複雑な地形や建物を含んだ大規模3次元爆風シミュレーション - J

複雑な地形や建物を含んだ大規模３次元爆風シミュレーション
青木 尊之（東工大・学術国際情報センター）
Takayuki AOKI, Tokyo Institute of Technology
FAX: 03-5734-3276
E-mail: [email protected]
１．はじめに
最近爆発事故のニュースが多いように感じるが、社会開
発において爆薬は欠くことができない存在である。強大な
破壊力を持つために貯蔵庫の安全性を検証することは極め
て重要である。しかし大規模な爆発実験には多くの危険が
伴うため、爆発をコンピュータ上で仮想的に実験する数値
シミュレーションが注目されている。
爆薬は化学反応の速度が非常に早く、放出されるエネル
ギーも膨大である。TNT などの強力な固体爆薬は起爆の瞬
間に爆轟が進み、一瞬のうちに化学反応が爆薬全体に広が
る。ほぼ固体状態のままで爆轟は終了するが、温度・圧力・
内部エネルギーは放出されたエネルギーのために極限とも
言える状態に上昇している。爆轟により生成された物質は
超音速で空気中を膨張して行き爆風となる。膨張とともに
次第に密度が低下してガス状になるが、依然として温度が
高いために炭素などが発光を続け燃えているように見える。
我々が目にする爆発の火炎面は爆轟生成ガスの表面である。
爆風による爆轟生成ガスの広がる範囲は限られているが、
爆轟生成ガスが空気中に超音速で膨張するために強い衝撃
波を発生し、遠方にまで伝播して被害をもたらす。爆轟生
成ガスと空気の界面は重力不安定（レーリー・テーラー不
安定性）であり、ずり流れの状態になるためにケルビン・
ヘルムホルツ不安定性も発生するなど流体力学的に非常に
不安定である。
２．IDO 法による高精度数値計算
爆風を詳細に解析するには、爆轟生成ガスの複雑な運動
に対して細かい計算格子を用いる必要があり、さらに遠方
にまで伝播する衝撃波の強さを精度よく評価することも重
要である。従って爆発（爆風）解析の計算領域も広くなる
ため、非常に大規模な計算が必要となる。
実験を行わずに爆発（爆風伝播）現象を計算機の中で忠
実に再現させるために、圧縮性流体方程式を時間・空間で
離散化し数値シミュレーションを行う必要がある。通常は
差分法などが用いられるが、爆轟直後の密度は空気の 1000
倍以上あり、ミリ秒オーダーで膨張して密度が空気の 100
分の 1 以下に低下する部分もある。このような急激な密度
変化に対して精度の高い計算を行うために局所補間微分オ
ペレータ(IDO)法 1,2) を用いる。通常の数値計算で用いる密
度や圧力などの物理量の他に、その空間勾配を従属変数と
して独立に計算し、それらを用いて局所空間に高次エルミ
ート補間を形成する。この補間関数から導かれる高階空間
微係数を有効に利用してコンパクトなステンシルで安定・
高精度な計算を実現している。また、密度や圧力が数値誤
差の原因で、少しでも負の値を取ると計算が続行不可能に
なるため、方程式中の対流項に対しては解の単調整を保障
する有理関数型のエルミート補間 3) が用いられている。一
方、時間発展に関しては、十分な空間精度と計算の安定性
を確保するために４段４次精度のルンゲクッタ法が用いら
れている。
３．極端な密度変化を伴う衝撃波計算の検証
解析解のある衝撃波管問題で、初期密度比 1000 倍と圧
力比 10000 倍の状態の隔壁を取り外して計算を開始させ
た結果が図 1 である。高密度側を爆轟生成ガスと同じ状態
にし、低密度側は空気を想定している。密度・圧力のオー
ダーが変わる現象で、空気中を伝播する衝撃波に対しても
解析解と良く一致する結果が得られている。
Fig.1 Shock tube problem with large density and
pressure jump
爆薬の地表爆発を想定した２次元円筒座標系の計算を行
った結果を示す。円筒状の TNT 32kg の爆轟直後を初期状
態とし、地面は剛体として扱っている。
Fig.2 Blast wave driven by the explosion on the ground.
図 2 は密度分布を示している。左図は球状衝撃波が先に広
がり、中心付近に不安定な爆轟生成ガスが漂うようすを示
している。右図は中心付近の詳細計算であり、高圧の爆轟
生成ガスが膨張し、中心付近が非常に低密度になる。この
領域の圧力は次第に低下し、周囲からガスが流れ込んでき
て中心付近で衝突することにより 2 次衝撃波が発生する。
４．模擬火薬庫からの爆風計算
実際は山の斜面などにトンネル状の火薬庫を作ること
が多い。実験でも図 3 のような筒状の火薬庫を想定するこ
とが多く、実験結果と比較するために同じ条件で 600×1100
格子を用い、円筒 2 次元計算を行った。
Fig.3
Schematic view of the explosive magazine.
火薬庫の出口か等距離で方向が 1.8 度、28 度、60 度、85
度の地点で測定された圧力波形（時間履歴）と計算結果の
比較を図 4 に示す。衝撃波の強さは方向で異なり、衝撃波
５．将来展望
ここまで爆風計算のみを紹介してきたが、安全性という
観点から最も重要なことは爆発により破壊された物質がど
のように飛散するかである。これは圧縮性流体運動と構造
物の変形・破壊・運動との連成問題であり、今後の大きな
課題と考えられる。
破壊された複雑な形状の物体と流体の相互作用を詳細に
計算するためには十分細かい空間解像度が必要となる。全
空間に細かい格子を用いると計算資源がいくらあっても足
りない。物体の近傍にのみ詳細な格子を配置する必要があ
るが、物体は移動や変形するので境界適合の構造格子では
対応できない。AMR(Adaptive Mesh Refinement)法は有望な
格子法であり、図6はAMR法を用いてレーリーテーラー不
安定性の成長を計算した結果である。重い流体と軽い流体
の界面に細かい格子を配置して複雑に変形する界面を詳細
に記述している。
-0.4
-0.2
Y
0
0.2
0.4
0
0.2
0.4
0.6
0.8
1
X
Fig.6
Fig.4 Comparison between the experimental pressure histories
and the numerical results.
の到達時間を除いて数値計算は実験を非常によく再現して
いることが分かる。爆轟生成ガスに対しては、JWL 状態方
程式、空気に対しては理想気体の状態方程式を適用し、圧
縮性流体の速度を元に体積分率を計算し、分圧で爆轟生成
ガスと空気の混合領域の圧力を決めている。
実際の複雑な地形や建物を考慮した爆風の影響を評価す
るためには 3 次元計算が必要であり、400×400×400 程度の
格子点で計算を行うために（独）産業技術総合研究所・爆
発安全研究センターの Myrinet2000 接続の PC クラスター
（128CPU:Linpack 492TFlops）を用い、領域分割による並
列化を行った。図 5 は爆轟生成ガスが流体力学的に不安定
な状態で噴出して行くようすをフォト・リアリスティック
に可視化したものである。5)
Fig.5
Growth of the Rayleigh-Taylor instability with AMR
method.
参考文献
1) T. Aoki: Interpolated Differential Operator (IDO) Scheme
for Solving Partial Differential Equations, Comput. Phys.
Commun., Vol.102, 132-146 (1997)
2) T. Aoki, S. Nishita, K. Sakurai: Interpolated Differential
Operator (IDO) Scheme and Application to Level Set
Method, Comput. Fluid Dynamics Journal , Vol.9, No.4,
406-417 (2001)
3) F. Xiao, T. Yabe and T. Ito: Constructing Oscillation
Preventing Scheme for the Advection Equation by a Rational
Function, Comp. Phys. Comm., Vol.93, 1-12 (1996)
4) 加藤香・青木尊之・久保田士郎・吉田正典 (2003a) : 爆
轟直後からの爆風伝播３次元数値解析, 平成 14 年度 衝
撃波シンポジウム講演論文集, Vol.14, 247-248
5) 青木尊之: 爆発シミュレーションとビジュアライゼー
ション, 画像ラボ, Vol.15, No.7, 7-9 (2004)
Eruption process of the detonation product gas by three-dimensional simulation