日航123便事故調査報告書付録「圧力隔壁からの与圧空気の流出の数値計算」を3次元流体解析で再現を試みる
こちらの日航123便事故調査報告書付録にはp53「圧力隔壁からの与圧空気の流出の数値計算による検討」の章があります。
この結果が事故原因が圧力隔壁開口と推定された大きな根拠となってます。
その計算が正しいのかどうかを現在の最新の3次元流体解析で解明します。
流体解析にはOpenFoam（バージョン13)というフリーのオープンソースライブラリを使用します。
このライブラリは産業界や研究でも広く利用されているもので、スパコンでの解析でも利用されてます。
ライブラリそのままでは使いにくいため、それを使いやすくした汎用的なオーサリングツール"CrazyEasyFluidAnalysisSystem"を独自開発しまして、それを利用して今回のシミュレーションを行っております。
OpenFoam自体は改変することなくオリジナルのまま使用してます。計算結果の一次データはOpenFoamから出力されたものとなってますが2次的に出力されたデータは独自の実装で算出したものもあります。
どの計算結果がOpenFoamによるものかというのはこちらのmakegraph.pyの説明の表に記載してます。出力されるCSVファイルの意味も説明があります。
今回の解析につきましては、だれでも無料で再現できるように全てダウンロードできるようになってます。
実行環境、各種解析ツール、実行スクリプトが揃ってますので実行の再現はできます。こちらのページからダウンロードできます。インストール手順はページ内で説明があります。
実行結果全体はサイズが大きいことからCSVファイルやグラフ画像のみを掲載します。
今回使用した3次元モデル定義のXMLファイルと設定の.configファイルと実行用のスクリプト.shファイルおよび解析用の.pyファイルらをzipにまとめました。こちらからjal123standard.zipはダウンロードできます。
解凍しますと４つのファイルがありますので、それらをC:/wcfdフォルダに配置したとしてdocker上のコマンドラインで実行するには、
（Dockerや"CrazyEasyFluidAnalysisSystem"および付属アプリは事前にインストールしておきます。OpenFoamは"CrazyEasyFluidAnalysisSystem"のDockerイメージアーカイブに含まれます。）
cd /home/user/wcfd ./jal123standard.shとするだけです。結果がjal123standardフォルダに作成されます。
なお、普通のWindowsノートPCで動作させてますので、特別高機能なハードウェアは必要ありません。
シミュレーション内容説明
 |  |  |
上図のような3次元モデルを作成しました。
後部圧力隔壁が開口した場合にAPU防火壁および垂直尾翼が破壊されるかどうかを検証します。
事故調査報告書付録p70では８つの区域に分けて計算を繰り返すことにより結果を得てます。
今回も8つの区域(room01からroom08および外気のroom00)に分けて、B747-100形状になるべく合わせて作成しました。
なおB747-100形状はこちらの62-2-JA8119-05.pdf事故調査報告書p140を参考にモデリングしたSTLファイルを利用してますが、
これは形状を分かりやすく重ねて表示するための素材であり、計算には影響してません。
また、各区域は空気移動ができる接続部としましてdoor（事故調査報告書のATに相当）を配置しました。
基本的には6面体で立体形状を作成してます。本来は曲面で実装すべきですが、複数の6面体を組み合わせることで曲面となるべく乖離が少ないようにしてます。
各room（事故調査報告書での第１から８室に相当）は単数または複数の６面体のboxで構成されてます。
roomおよびdoorの容量や耐圧等のデータは事故調査報告書付録に準じてます。垂直尾翼形状は事故調査報告書付録の標準的なケースで採用されているモード３を実装しました。
| p64 | p65 | p71 |
 |  |  |
容量や温度や開口部面積や耐圧等の一覧が参照できます。doorの名前になる数字4桁は前2桁と後ろ2桁が接続する部屋の番号を示してます。
door0102が圧力隔壁開口部に相当しまして、面積は報告書の基準ケースと同様の1.8m2となってます。
結論
流速の解析動画(3秒間を0.1倍速度で再生）はこちらです。3次元流体解析を行った結果、圧力隔壁開口(door0102)での空気噴流ではAPU防火壁(door0600)は破壊されたものの、垂直尾翼(door0300)は破壊されませんでした。
結果のCSVファイルはこちら(resultbox.csv)とこちら(resultdoor.csv)です。
各door（開口部を表現）での差圧が耐圧を超えたかどうかはresultdoor.csvのbrokenの列で出力されてます。brokenが1であれば破壊。0であれば破壊なしとなってます。
垂直尾翼(door0300)はresultdoor.csvで下記の強調部分で最大差圧とその時刻が記録されてます。
気圧に関してはOpenFoamのほうで計算された結果となってます。
これによりますと、開口から0.319秒後に最大差圧2.64psiとなっていることが確認できますが設定した耐圧4.0psiを超えてません。つまり破壊されておりません。
参考までに厳密には垂直尾翼耐圧は4.75psiなのですが差圧4.0psiで破壊可能と報告書で扱っているので今回のシミュレーションでは耐圧として4.0psiを設定してます。
差圧をグラフにしますと上記のようになります。青いグラフのdoor0600がAPU防火壁にかかる差圧で、オレンジのグラフのdoor0300が垂直尾翼にかかる差圧です。
APU防火壁は圧力隔壁開口直後に破壊されるため、そこから空気が抜けて一気に差圧が下がっていることが確認できます。
垂直尾翼のほうはAPU防火壁が破壊後にある程度は差圧が上昇するのですが、4.0psiの設定耐圧までは上昇しません。
この結果は実はある程度は予測できていたものです。
以前には今回の3次元解析ではなく事故調査報告書と同様の次元のない数値解析で試したことがありますが、今回と同様に垂直尾翼は破壊されませんでした。
以前の独自シミュレーションはこちら(ver1 20秒ほどで結果が表示されます)とこちら(ver2 1分ほどで結果が表示されます)でWebブラウザで実行できます。
以前の独自シミュレーションは報告書とは異なる計算式を独自に実装してますが報告書と同様に８つの場所の状態を計算してます。
今回の3次元シミュレーションはOpenFoamのほうで3次元メッシュ全体で計算を行ってます。
報告書での計算はあくまで気圧や気温の状態遷移を計算していて形状や乱流は考慮されてませんので、3次元シミュレーションと結果が異なる理由としましては、その違いもあるのだろうと思います。
一般的な解釈では現在の3次元流体シミュレーションのほうが実際の状況に近い結果が出るとされますから、当時の事故調査結果のほうが現実との乖離が大きいとみなされると思います。
事故原因は少なくとも圧力隔壁開口ではないということを今回の現代科学のシミュレーションは示してます。
分析
事故調査報告書付録p76での客室の気圧変化のグラフを示します。赤字を追加してます。縦軸単位をpsiにしてます。これで比較がしやすいです。
今回のシミュレーションでの客室後部での気圧変化をグラフで見てみます。3秒間の表示です
比較しますと、今回のシミュレーションのほうが減圧速度がやや速いことがわかります。
ところで、事故調査報告書付録p81の「付図6 耐圧限界値係数と隔壁開口面積の変化」のグラフは以下です。
グラフの斜線部分の範囲内であればAPU防火壁と垂直尾翼の両方が破壊されるという意味になってます。
横軸が隔壁開口面積ですから、開口が大きいければ大きいほど減圧速度が速くなり両方が破壊可能な確率は高くなるということです。
事故調査報告書の減圧速度よりも今回のシミュレーションのほうが減圧速度がやや速いということは、
それだけAPU防火壁と垂直尾翼の両方が破壊される確率は上がっているはずなのです。
しかし、結果としましては垂直尾翼は破壊に至りませんでしたので、仮に事故調査報告書と同様の減圧速度であればやはり破壊には至らないということになります。
次に、APU防火壁から後方へ抜ける空気噴流により機体の前方加速度を算出してみます。
resultdoor.csvにはvelocity_mpsおよびmassFlow_kgpsの列があります。APU防火壁を破壊後に通過する流速や質量流量はdoor0600の行にありますので、そのデータを抽出して、
脱落していない残りの機体重量を234264kgとしまして、空気噴流により生じた機体の前方加速度を求めます。
APU防火壁開口部（door0600）における流速(velocity_mps 単位m/sec)および質量流量(massFlow_kgps 単位kg/sec)のグラフを示します。
 |  |
jal123standard.zipには加速度算出のためのスクリプトacalcu.pyが含まれてます。
実行するには
python ./acalcu.py jal123standard/resultdoor.csv --name door0600 --mass 234264とします。resultdoor.csvファイルがacalcu.pyと同じフォルダに配置している場合はresultdoor.csvのファイル名のみの指定で結構です。
結果のグラフはこちらです。
下の事故調査報告書付録p92「付録5付図１ DFDR拡大図」のLNGG（前後方向加速度）のグラフに示されている通り前方加速度の突出は0.047gでしたが、その2倍以上の値が上記グラフではプロットされてます。
ただし、これはAPU防火壁から抜ける空気噴流がすべて前方加速度として使われたという前提ですので、理論上の最大値とみなすことが出来ます。
このグラフで分かることは、圧力隔壁開口の空気噴流が前方加速度を発生させたと仮定しますと、２秒間程度は有意な加速度突出が記録されるはずということです。
しかしながら、実際のDFDRでの前方加速度の突出は約0.25秒間(最大範囲で0.5秒間）の範囲で収まってます。
つまり、圧力隔壁開口の空気噴流では前方加速度の一時的な突出は説明が困難です。
もっと短い時間での一瞬の噴流のほうがDFDRの記録とは整合性がとれると考えられます。
予備試験
今回のシミュレーションの前に単純な形状モデルでの断熱膨張の実験を行いました。
こちらのページ内で説明しているoneroomというフォルダ内にtest1,test2,test3,test4の４つのテストの実験です。
室内気圧および外気圧は日航123便の状況と同様ですが、容積は異なります。
 |  |
気圧変化や温度変化のグラフからtest2の設定を参考に今回のシミュレーションを行いました。
なお、test2の温度グラフでは室内が-43度にまで下がっている様子となってます。
事故調査報告書付録p74の「付図4基準ケース」の客室温度低下とほぼ同じ（時間スケールは異なる）曲線となってます。
外気圧まで室内の気圧が下がった場合に、外気温度よりも低い-43度にまで下がる状況が計算上は確かであることを確認しました。
応用
 |  |
第一室(room01)は実際には天井で客室と天井裏で仕切られている構造です。また、客室の最後部には化粧室がありますし、座席もあります。
それらがどのように影響するかというのを簡易モデルを追加して確認します。
図の位置がbox0104という後部座席の計測位置となります。また、隔壁近くの天井穴としまして、赤い囲いの部分の面積を開けました。
他には、前方部分で天井と客室は通じている構造としました。実際には小さい通気口が至る所にあると考えられます。
実行するにはこちらのzipファイルにxml,config,shファイルがありますので、解凍してご利用ください。C:/wcfdに配置した場合dockerのLinuxのコマンドライン上では（他の実行環境はインストール済としまして）
./jal123detail.shと実行することでjal123detailフォルダ内に結果が出力されます。グラフや動画は付属ソフトで作成できます。結果確認はWindows用のOpenFoamの標準的なビューアであるParaViewで確認可能です。
以下がbox0104位置での3秒間のグラフです。気圧、温度、風速、湿度となってます。
 |  |
 |  |
湿度は（単純計算では）100%以上が霧となって見える状態の意味です。
風速のグラフを確認しますと、はじめの１秒間で、2m/secから9m/sec程度となってます。
このモデルでは通路無しとしてモデリングしてますが動画のほうで確認しますと狭い部分で乱流も見られますので、通路脇の位置ではもっと風速はあがると考えられます。
気温のグラフでは、やはり断熱膨張では外気温以下に下がることが確認できます。この３秒間以上のシミュレーションを行った場合はもっと下がるものと考えられます。
動画（３秒間を１０倍スロー）は以下のリンクです。
気圧の変化のアニメーション
温度の変化のアニメーション
流速の変化のアニメーション
こちらがCSVファイルです。垂直尾翼であるdoor0300は破壊されませんでした。door0300の最大差圧は0.319秒時に2.82psiでした。
resultbox.csv
resultdoor.csv