SENAC Vol. 55, No. 2（2022. 4）― 25 ―［共同研究成果］

直交格子 法による 移動境界 問題の解 法に関す る研究高橋俊：東海大学工学部動力機械工学科
佐々木大輔：金沢工業大学工学部航空システム工学科
山田剛治 ：東海大学工学部機械工学科

Euler 方程式と断熱壁面，滑り壁面境界条件を用いた埋め込み境界法により，広範囲の速度域を安定に解析できる解析手法を構築した．極超音速流中における円柱と球の解析を実施して，先行研究の数値解析結果や，実験結果との定性的な一致を確認した．また複雑形状の例として２球から生じる衝撃波が干渉する流れ場の解析を行って，本手法の適用性について確認した．今後は移動境界問題に適用すべく本手法の拡張を行い，最終的には多数の球と衝撃波が干渉した際の運動量，エネルギーの変換や，流れ中において様々な相対移動を行う物体周りの解析と現象解明に取り組む．

1.

はじめに流体中を伝搬する衝撃波が物体と干渉する現象は以前から幅広く研究されており，爆発や衝撃による物体の破壊から，衝撃波が誘起する相変化[1]や，体内の結石の破砕[2, 3]など様々な問題設定の研究が行われてきた．衝撃波が物体表面を通過する時に物体に作用する瞬間的な非定常抗力は衝撃波と物体の干渉現象において特徴的な現象の一つとして挙げられる．Tanno[4]や Sun[5]らは単一の球に衝撃波を負荷した時に生じるこの非定常抗力について実験と数値計算によって計測，調査し，非定常抗力係数は定常時の抗力係数の数倍に達することを明らかにした．またこれらの研究では非定常抗力が生じる際の衝撃波の可視化から，衝撃波が球体を通過して後流側よどみ点で収束した際に抗力が大きく低下する現象や，Re 数によっては瞬間的に抵抗が負になる現象を示した．また数値解析では Re 数依存性の議論から Re 数が大きい場合には非粘性解析でも十分に非定常抗力が議論できると結論付けられた．彼らの研究はこれまで多くの非定常抗力の研究者から引用され，現在まで世界中で詳細なメカニズムが研究されている．例としては衝撃波を微粒子に当てた際の軌道解析[6]，衝撃波管の境界層の影響なども考慮して計測を高精度化した実験的な研究[7]や，複数の粒子の流体力学的相互作用に関する CFD 研究[8, 9]などが挙げられる．これらの研究の応用の一つが粒子幕（particle curtain，particle cloud）のモデル化と現象解明である．質点モデルを CFD に組み入れた研究[10-13]では，質点モデルの高度化とそれが運動量交換，エネルギー交換に及ぼす影響について検討されている．実験での粒子幕の研究には[14, 15]などが存在し，CFDによる粒子幕の研究では埋め込み境界法（Immersed Boundary Method: IBM）を用いて，乱流エネルギー生成と散逸，運動量交換，エネルギー交換などが詳細に議論されてきた[16-18]．実際に非定常抗力や定常抗力を受けた微粒子は流れの中で加速度運動して流れ場も同時に変化する．IBMを用いた非粘性 CFD 解析で粒子の移動を調査した研究例も存在する[8, 19]が，CFD で粒子の運動の影響を議論した研究例は少ない．これまでは単一球周りの高 Ma 数，低 Re 数の流れ場の詳細な数値解析によるデータベース構築の研究[20-23]や，実験的研究[24]，移動物体周りの CFD による研究[25]，CFD を用いた抵抗モデルの高精度化[26]と，継続的に傾向を調査してきた．本研究ではさらなる粒子幕に関する研究の深化に向けて，1 つもしくは 2 つの粒子周りの流れ場を IBM により解析する計算手法の構築と検証計算を行った．

― 26 ―

2.

SENAC Vol. 55, No. （2 2022. 4）

計算方法

2.1 支配方程式
支配方程式には無次元化された三次元圧縮性 Euler 方程式を用いる．(1)，(2)，(3)はそれぞれ連続の式，運動量保存式，エネルギー保存式である．
𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕(𝜕𝜕𝜕𝜕𝑢𝑢𝑢𝑢𝑖𝑖𝑖𝑖 )+=0𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕
𝜕𝜕𝜕𝜕𝑥𝑥𝑥𝑥𝑖𝑖𝑖𝑖
𝜕𝜕𝜕𝜕(𝜕𝜕𝜕𝜕𝑢𝑢𝑢𝑢𝑖𝑖𝑖𝑖 ) 𝜕𝜕𝜕𝜕�𝜕𝜕𝜕𝜕𝑢𝑢𝑢𝑢𝑖𝑖𝑖𝑖 𝑢𝑢𝑢𝑢𝑗𝑗𝑗𝑗 �𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕
+=−𝜕𝜕𝜕𝜕𝑥𝑥𝑥𝑥𝑗𝑗𝑗𝑗𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝑥𝑥𝑥𝑥𝑖𝑖𝑖𝑖
𝜕𝜕𝜕𝜕(𝜕𝜕𝜕𝜕𝜌𝜌𝜌𝜌)𝜕𝜕𝜕𝜕+�(𝜕𝜕𝜕𝜕𝜌𝜌𝜌𝜌 + 𝜕𝜕𝜕𝜕)𝑢𝑢𝑢𝑢𝑗𝑗𝑗𝑗 � = 0𝜕𝜕𝜕𝜕𝜕𝜕𝜕𝜕
𝜕𝜕𝜕𝜕𝑥𝑥𝑥𝑥𝑗𝑗𝑗𝑗

(1)

(2)

(3)

ここで位置 xi は球形粒子の直径 d，流速 ui は基準音速 c0，時間 t は無次元時間 d/c0，総エネル
ギーE は基準音速の二乗 c02，密度 ρ は基準密度 ρ0，圧力 p は ρ0c02 でそれぞれ無次元化されている．また Euler 方程式内で用いる諸量を求めるための無次元化された状態方程式は以下の式で表される．
1
𝜕𝜕𝜕𝜕 = (𝛾𝛾𝛾𝛾 − 1) �𝜕𝜕𝜕𝜕𝜌𝜌𝜌𝜌 − 𝜕𝜕𝜕𝜕𝑢𝑢𝑢𝑢𝑖𝑖𝑖𝑖 𝑢𝑢𝑢𝑢𝑖𝑖𝑖𝑖 �2(4)

計算格子には等間隔直交格子を用いて，非粘性流束計算には SLAU スキーム[31]を用い，時間積分は 3 次精度の Total Variation Diminishing (TVD) Runge-Kutta 法を用いて行う．物体境界はレベルセット関数と埋め込み境界法を用いて表現する．これらの手法の計算精度は著者らの先行研究で既に実施済みである[27, 28]．また本研究では粒子の運動方程式との連成解析を行うべく進めている．その際，流体力の算出には先行研究により提案されたポリゴン形状を不要とする簡易な方法を用いる[29]．Navier-Stokes 方程式を用いた際の，この手法を応用した運動方程式との連成解析の精度については既に検証済みである[25, 30]．

2.2 埋め込み境界法
これまでの研究では滑り無し壁面境界条件を用いた検証結果を報告していたが，本研究では高Re 数の流れ場を模擬した解析を高効率に実施するために滑り壁面境界条件の埋め込み境界法を適用する．図 1 が滑り壁面用の埋め込み境界法の模式図で，黒，黒縁灰色．灰色，白抜き緑，黄色の丸印がそれぞれ計算セル，ゴーストセル，物体セルの中心位置，壁面上の点，イメージポイントを示す．黄色，緑，紫，赤，青の線と矢印がそれぞれ壁面からの距離，物体境界線，速度ベクトル，速度ベクトルの壁面法線方向成分，接線方向成分である．φg，φi はそれぞれ着目しているゴーストセルと壁面上の点までの最小距離と，イメージポイントと壁面上の点までの最小距離（プローブ長さ）で，今回はイメージポイントにおける再帰参照を防止するため φi=1.75Δx（Δx は格子幅）の固定値とした．イメージポイントには trilinear 内挿で周囲の 8 点から物理量を内挿する．

直交格子法による移動境界問題の解法に関する研究

図1𝒗𝒗𝒗𝒗𝑔𝑔𝑔𝑔,𝑛𝑛𝑛𝑛 = −― 27 ―

滑り壁のゴーストセルの速度の定義方法

𝜙𝜙𝜙𝜙𝑔𝑔𝑔𝑔
�𝒗𝒗𝒗𝒗 ∙ 𝒏𝒏𝒏𝒏�𝒏𝒏𝒏𝒏
𝜙𝜙𝜙𝜙𝑖𝑖𝑖𝑖 𝑖𝑖𝑖𝑖,𝑛𝑛𝑛𝑛(5)

𝒗𝒗𝒗𝒗𝑔𝑔𝑔𝑔,𝑡𝑡𝑡𝑡 = 𝒗𝒗𝒗𝒗𝑖𝑖𝑖𝑖,𝑛𝑛𝑛𝑛 − �𝒗𝒗𝒗𝒗𝑖𝑖𝑖𝑖,𝑛𝑛𝑛𝑛 ∙ 𝒏𝒏𝒏𝒏�𝒏𝒏𝒏𝒏(6)

𝜕𝜕𝜕𝜕𝑔𝑔𝑔𝑔 = 𝜕𝜕𝜕𝜕𝑖𝑖𝑖𝑖(8)𝜌𝜌𝜌𝜌𝑔𝑔𝑔𝑔 =𝒗𝒗𝒗𝒗𝑔𝑔𝑔𝑔 = 𝒗𝒗𝒗𝒗𝑔𝑔𝑔𝑔,𝑛𝑛𝑛𝑛 + 𝒗𝒗𝒗𝒗𝑔𝑔𝑔𝑔,𝑡𝑡𝑡𝑡(7)

𝜕𝜕𝜕𝜕𝑔𝑔𝑔𝑔 = 𝜕𝜕𝜕𝜕𝑖𝑖𝑖𝑖

(9)

1 𝜕𝜕𝜕𝜕𝑔𝑔𝑔𝑔 1+ �𝒗𝒗𝒗𝒗 ∙ 𝒗𝒗𝒗𝒗 �𝛾𝛾𝛾𝛾 𝛾 1 𝜕𝜕𝜕𝜕𝑔𝑔𝑔𝑔 2 𝑔𝑔𝑔𝑔 𝑔𝑔𝑔𝑔

(10)

下添え字 i，g，n，t はそれぞれイメージポイント，ゴーストセル，壁面法線方向成分，壁面接線方向成分の諸量を指す．イメージポイントに内挿した速度 vi，密度 ρi，圧力 pi を用いて，式(5)(10)を順に計算して断熱条件と滑り壁条件を満足するようにゴーストセルの諸量である速度 vg，密度 ρg，圧力 pg，エネルギーEg を求める．本研究では物体境界をレベルセット関数によって表現しているため，壁面との単位法線ベクトル n はレベルセット関数の勾配から容易に求められる．
これまでの研究では滑り無し壁面境界条件を適用していたが，本研究では高 Re 数の流れ場を模擬した解析を高効率に実施するために，滑り壁面境界条件の埋め込み境界法を適用した．また今回は物体境界からプローブ長さの 2 倍程度まで物体の内部に位置しているセルをゴーストセルとした．

― 28 ―

3.

SENAC Vol. 55, No. （2 2022. 4）

結果と考察

3.1 円柱周りの解析（M∞=8.1）
構築した手法の検証として一様流マッハ数 M∞=8.1 の中に置かれた円柱周りの解析を行った．
等間隔直交格子を用いて主流方向と主流垂直方向にそれぞれ 100，200 セル，合計 20000 セルを用いて，円柱の直径を 80 セルで分割した．図 2 は用いた格子と解析により得られた等圧力線図を示しており，この結果が構造格子によって得られた先行研究[31]と良好な一致を示した．このことか
ら本解析手法の妥当性と有効性を確認した．

図2

M∞=8.1 における格子と円柱周りの等圧力線図

3.2 球周りの解析（M∞=6.7）
一様流マッハ数 M∞=6.7 の流れ中の球周りの流れ場を，直径を 40 セルと 80 セルで分割して解析を行った．図 3，図 4 はそれぞれの解析から得られた中央断面における等圧力線図と密度勾配の絶対値の分布図である．両者とも離脱衝撃波が安定して捉えられ，80 セルの方が等圧力線が滑らかで，また密度勾配の絶対値の濃淡が明瞭であった．衝撃風洞で本解析と同様の条件で可視化実験を実施した結果，衝撃波離脱距離が両者とも可視化結果とほぼ一致することを確認しており，これらの解析結果が妥当であることを確認した．

直交格子法による移動境界問題の解法に関する研究

図3

図4― 29 ―M∞=6.7 における球周りの等圧力線図（直径 40 分割，80 分割）

M∞=6.7 における球周りの密度勾配の絶対値分布図（直径 40 分割，80 分割）

3.3 ２球周りの解析（M∞=6.7）
一様流マッハ数 M∞=6.7 の流れ中の２球周りの流れ場を，直径を 40 セルで分割して解析を行った．図 5，図 6 はそれぞれ中央断面における等圧力線図と，密度勾配の絶対値の分布図である．前方の球周りで生じた離脱衝撃波が後方の球に当たり流れ場中で最大の圧力を生じた．図 5 の等圧力分布図はこの最大圧力で無次元化して示してある．球の後流では最大圧力の 1%未満の領域も確認される広範囲の速度と圧力が生じる運動量変化が大きな流れ場であったが，解析は安定に実施できたことから本手法の堅牢性が確認できた．図 6 に示した密度勾配の絶対値の分布には２球から生じる衝撃波の干渉と衝撃波の三次元的な分布が見られた．本解析で採用した埋め込み境界法
はこのような三次元任意形状の複数物体周りの計算格子作成が容易に行えることから，これらの三次元的な流れ場構造の把握に対する適用性が確認された．

― 30 ―

SENAC Vol. 55, No. （2 2022. 4）

図5

図64.

M∞=6.7 における２球周りの無次元圧力分布と無次元等圧力線図

M∞=6.7 における２球周りの密度勾配の絶対値分布（中央断面，前方から見た等値面）

結論

Euler 方程式と断熱壁面，滑り壁面境界条件を用いた埋め込み境界法により，広範囲の速度域を安定に解析できる解析手法を構築した．極超音速流中における円柱と球の解析を実施して，先行研究の数値解析結果や，実験結果との定性的な一致を確認した．また複雑形状の例として２球から生じる衝撃波が干渉する流れ場の解析を行って，本手法の適用性について確認した．今後は移動境界問題に適用すべく本手法の拡張を行い，最終的には多数の球と衝撃波が干渉した際の運動量，エネルギーの変換や，流れ中において様々な相対移動を行う物体周りの解析と現象解明に取り組む．

直交格子法による移動境界問題の解法に関する研究― 31 ―

謝辞
本研究は東北大学サイバーサイエンスセンターのスーパーコンピュータを利用することで実現することができた．また研究にあたって同センター関係各位に有益なご指導とご協力をいただいた．この場をお借りして謝意を表する．
また本解析手法の開発において，東京大学大学院工学系研究科航空宇宙工学専攻航空宇宙システム学講座 今村太郎 准教授には貴重な意見をいただいた．ここに謝意を表する．

参考文献[1] Duvall, G. E., et al. “Phase transitions under shock-wave loading.” Reviews of Modern Physics 49.3
(1977): 523.
[2] Chaussy, C., et al. “First clinical experience with extracorporeally induced destruction of kidney stones
by shock waves.” The Journal of urology 127.3 (1982): 417-420.
[3] Takayama, K., et al. “Applications of shock wave research to medicine.” WIT Transactions on Modellingand Simulation 41 (2005).
[4] Tanno, H., et al. “Interaction of a shock with a sphere suspended in a vertical shock tube.” Shock Waves
13.3 (2003): 191-200.
[5] Sun, M., et al. “Unsteady drag on a sphere by shock wave loading.” Shock waves 14.1 (2005): 3-9.
[6] Jourdan, G., et al. “Drag coefficient of a sphere in a non-stationary flow: new results.” Proceedings of theRoyal Society A: Mathematical, Physical and Engineering Sciences 463.2088 (2007): 3323-3345.
[7] Wagner, Justin L., et al. “Shock tube investigation of quasi-steady drag in shock-particle interactions.”
Physics of Fluids 24.12 (2012): 123301.
[8] Sridharan, P., et al. “Shock interaction with one-dimensional array of particles in air.” Journal of Applied
Physics 117.7 (2015): 075902.
[9] Mehta, Y., et al. “Numerical investigation of shock interaction with one-dimensional transverse array ofparticles in air.” Journal of Applied Physics 119.10 (2016): 104901.
[10] Ling, Y., et al. “Interaction of a planar shock wave with a dense particle curtain: Modeling andexperiments.” Physics of Fluids 24.11 (2012): 113301.
[11] Regele, J. D., et al. “Unsteady effects in dense, high speed, particle laden flows.” International Journalof Multiphase Flow 61 (2014): 1-13.
[12] Jiang, L.-J., et al. “Particle-resolved simulations of shock-induced inviscid flow through particle-curtain
at initial stage.” Computers & Fluids (2021): 105196.
[13] Shallcross, G. S., et al. “A volume-filtered description of compressible particle-laden flows.”
International Journal of Multiphase Flow 122 (2020): 103138.
[14] Bordoloi, A. D., et al. “Relaxation drag history of shock accelerated microparticles.” Journal of Fluid
Mechanics 823 (2017).
[15] DeMauro, E. P., et al. “Unsteady drag following shock wave impingement on a dense particle curtain
measured using pulse-burst PIV.” Physical review fluids 2.6 (2017): 064301.
[16] Das, P., et al. “A sharp interface Cartesian grid method for viscous simulation of shocked particle-ladenflows.” International Journal of Computational Fluid Dynamics 31.6-8 (2017): 269-291.
[17] Hosseinzadeh-Nik, Z., et al. “Investigation and quantification of flow unsteadiness in shock-particle
cloud interaction.” International Journal of Multiphase Flow 101 (2018): 186-201.
[18] Osnes, A. N., et al. “Computational analysis of shock-induced flow through stationary particle clouds.”
International Journal of Multiphase Flow 114 (2019): 268-286.
[19] Luo, K., et al. “Studies on shock interactions with moving cylinders using immersed boundary method.”
Physical Review Fluids 2.6 (2017): 064302.

― 32 ―SENAC Vol. 55, No. （2 2022. 4）[20] Nagata, T., et al. “Investigation on subsonic to supersonic flow around a sphere at low Reynolds numberof between 50 and 300 by direct numerical simulation.” Physics of Fluids 28.5 (2016): 056101.
[21] Nagata, T., et al. "Direct numerical simulation of flow past a transversely rotating sphere up to aReynolds number of 300 in compressible flow." Journal of Fluid Mechanics 857 (2018): 878-906.
[22] Nagata, T., et al. “Direct numerical simulation of subsonic, transonic and supersonic flow over an
isolated sphere up to a Reynolds number of 1000.” Journal of Fluid Mechanics 904 (2020).
[23] Nagata, T., et al. “Direct numerical simulation of flow around a heated/cooled isolated sphere up to aReynolds number of 300 under subsonic to supersonic conditions.” International Journal of Heat andMass Transfer 120 (2018): 284-299.
[24] Nagata, T., et al. “Experimental investigation of transonic and supersonic flow over a sphere forReynolds numbers of 103–105 by free-flight tests with schlieren visualization.” Shock Waves 30 (2020):
139-151.
[25] Mizuno, Y., et al. “Direct numerical simulation of gas–particle flows with particle–wall collisions using
the immersed boundary method.” Applied Sciences 8.12 (2018): 2387.
[26] Loth, E., et al. “Supersonic and Hypersonic Drag Coefficients for a Sphere.” AIAA Journal 59.8 (2021):
3261-3274.
[27] Takahashi, S., et al. “A numerical scheme based on an immersed boundary method for compressible
turbulent flows with shocks: application to two-dimensional flows around cylinders.” Journal of Applied
Mathematics 2014 (2014): 252478.
[28] Mizuno, Y., et al. “A simple immersed boundary method for compressible flow simulation around a
stationary and moving sphere.” Mathematical Problems in Engineering 2015 (2015): 438086.
[29] Nonomura, T., et al. “A comparative study on evaluation methods of fluid forces on Cartesian grids.”
Mathematical Problems in Engineering 2017 (2017): 8314615.
[30] Nagata, T., et al. “A simple collision algorithm for arbitrarily shaped objects in particle‐resolved flow
simulation using an immersed boundary method.” International Journal for Numerical Methods in Fluids
92.10 (2020): 1256-1273.
[31] Kitamura, K., et al. “Towards shock-stable and accurate hypersonic heating computations: A newpressure flux for AUSM-family schemes.” Journal of Computational Physics 245 (2013): 62-83.