溶鋼内の介在物移動及び気泡運動に対する埋め込み格子法の発展と応用
代表研究者 産業技術短期大学 機械工学科 講師 浅尾 慎一
共同研究者 産業技術短期大学 機械工学科 教授 樋口 善彦
1．研究背景鉄鋼製品の安全性の基準の 1 つとして，非金属介在物や気泡の含有率が挙げられる．非金属介在物とは溶鋼接触面の耐火物やスラグ，溶鋼中に含まれる金属成分が空気中の酸素と反応して生じた生成物である．非金属介在物や気泡は鋼よりも密度が低いため，通常は溶鋼表面に向けて浮いてくるが，あまりにも小さいため浮上速度が非常に遅く時間がかかる．そのため，浮上しきらない非金属介在物や気泡（以下，介在物と呼ぶ）はタンディッ
シュから鋳型に取り込まれて鋼中に取り込まれていく．
浮遊速度に関しては，理論的にストークスの式によって導き出すことができる．
しかし，これは静止している溶鋼中に対する介在物や気泡の浮遊速度である．現実の溶鋼は流動しており，その流動によって，介在物の移動が行われるため，浮遊速度はストークスの式と大きく異なる．
今，溶鋼内を自由に動く介在物のような，大きい計算領域を物体が自由に動く様子を計算する場合，手法としては大きく２つに分かれる．一つは計算格子を用いて物体境界を表現する方法であり，物体適合座標を用いる方法や重畳格子法[1]を用いる方法である．一方，計算格子を用いずに物体境界を表現する方法として仮想境界法[2]や，カットセル法[3]などがある．これらは物体表面に関係なく計算格子を生成し，物体付近の格子点に仮想的な外力を与えることで物体表面を表現する方法である．前者の方法は物体表面を鋭く捉えるため，物体表面に対する計算の信頼性が高い．しかしながら，移動する物体の格子生成は難しく，格子生成に時間がかかってしまう．また，計算格子を物体にあわせて張るため，格子の形状がいびつになる部分が発生し，その部分では計算精度が落ちる可能性がある．重畳格子法の場合は計算領域全体に張るメイン格子と物体周りに生成するサブ格子が独立しているため，上記のようないびつな格子が生成する可能性は低い．しかしながら，互いの格子の境界での物理量を内挿するため，計算精度の低下，流れの保存則を満足しない．一方，後者の方法では格子生成をする必要がないため，物体表面の生成は容易である．しかしながら，物体表面を計算格子で表さないため，物体を鋭く捉えることができず計算の信頼性が低い．さらに物体を動かす際，流れの保存則に加えて幾何保存則[4]も考慮する必要がある．ここで幾何保存則とは格子の移動・変形が流れに影響を及ぼさないことをいう．
以上の背景より，研究者は重畳格子法のように格子生成に時間がかからず，かつ，流れの保存則と幾何保存則を同時に満足する手法である物体透過格子法[5]を開発してきた．この手法は移動格子有限体積法[6]という流れの保存則と幾何保存則を同時に満足する手法と，計算格子の追加・削除を行うことでスムーズに物体を移動させるものである．本研究－ 129 －

ではこの物体透過格子法を用いることで介在物の浮遊の様子を正確に捉えることができると考える．
以上より，本研究では上記の物体透過格子法を用いて流動解析を行うが，特にタンディッシュ内の堰の有無と介在物の動きについて調べる．ただし，タンディッシュ内の介在物の流れの特徴を端的に捉えるため，2 次元ベースで計算することにする．また介在物の個数も数個程度で行うことにする．
2．基礎方程式
2.1 流体に関する基礎方程式
基礎方程式は，連続の式と 2 次元非圧縮性ナビエ・ストークス方程式である．これらの式を無次元化したのち，保存形で記述すると，連続の式は式(1)，2 次元非圧縮性ナビエ・ストークス方程式は式(2)，(3)として記述される．
𝜕𝑢 𝜕𝑣
+= 0,
𝜕𝑥 𝜕𝑦(1)𝜕𝒒 𝜕𝑬 𝜕𝑭++= 𝟎,𝜕𝑡 𝜕𝑥 𝜕𝑦(2)𝑬 = 𝑬 + 𝑬 − 𝑬 ,𝑭 = 𝑭 + 𝑭 − 𝑭 ,
𝒒=

𝑝0𝑢𝑢𝑣
,𝑬 =,𝑭 =,, 𝑬𝒂 = 𝑢 , 𝑭𝒂 =𝑣
𝑝𝑣
0𝑢𝑣
𝜕𝑢𝜕𝑢
1 𝜕𝑥
1 ⎛𝜕𝑦⎞, 𝑭𝒗 =𝑬𝒗 =.
𝑅𝑒 𝜕𝑣
𝑅𝑒 ⎜𝜕𝑣 ⎟𝜕𝑥⎝𝜕𝑦⎠(3)

ここで，式中の 𝒒 は保存量ベクトル，𝑬𝒂 , 𝑭𝒂 は 𝑥, 𝑦 方向の移流項ベクトル，𝑬 , 𝑭は圧力項ベクトル，𝑬𝒗 , 𝑭𝒗 は 𝑥, 𝑦方向の粘性項ベクトルを示す．また，𝑢, 𝑣 は 𝑥, 𝑦 方向の速度，𝑝 は圧力，𝑅𝑒 はレイノルズ数を示している．また，これらの式は式(4)で表される関係式を用いて無次元化を行った結果である．
𝑥=𝑥̅
𝑦𝑢𝑣̅
𝑡̅
𝑝̅,𝑦 = ,𝑢 =,𝑣 =,𝑡 =,𝑝 =,𝐿𝑈
𝐿𝑈𝐿 ⁄𝑈𝜌̅ 𝑈𝑈𝐿𝜌̅ 𝑈 𝐿𝑅𝑒 ==.
𝜇̅𝜈̅(4)

ここで， は有次元量を表し，𝐿 , 𝑈 , 𝜌̅ ，𝜇̅ ，𝜈̅ はそれぞれ代表長さ，代表速度，密度，粘性係数，動粘性係数を示している．

－ 130 －

2.2 物体運動の支配方程式
介在物の運動として 2 次元の並進と回転を考慮した 3 自由度の運動方程式を用いる．流
体の支配方程式と同様の，空間に固定された空間固定座標系 𝑥, 𝑦 を用いた並進運動方程式は式(5)で記述される．
⎧𝑀 𝑑𝑢 = 𝐹⎪ 𝑑𝑡.
𝑑𝑣
⎨𝑀=𝐹⎪ 𝑑𝑡⎩
ここで，𝑀 は物体の重量，𝑢 , 𝑣(5)はそれぞれ 𝑥, 𝑦 方向の物体重心の速度，𝑥 , 𝑦 はそれ

ぞれ 𝑥, 𝑦 座標での物体重心，𝐹 , 𝐹 はそれぞれ 𝑥, 𝑦 方向の物体が受ける力を表す．
また，物体の回転運動方程式は式(6)で記述される．
𝐼𝑑𝜔=𝑇.
𝑑𝑡(6)

ここで，𝐼 は物体の慣性モーメント，𝜔 は物体の角速度，𝑇 は物体が受けるトルクを表す．ただし，物体は剛体とする．
また，これらの式は式(7)で表される関係式を用いて無次元化を行った結果である．
𝑀=𝐹𝑀𝑢𝑣̅
𝐹,𝑢 =,𝑣 =,𝐹 =,𝐹 =,𝑈𝑈𝜌̅ 𝐿 𝐿𝜌̅ 𝑈 𝐿 𝐿𝜌̅ 𝑈 𝐿 𝐿𝐼̅𝜔
𝑇𝐼 =,𝜔 =,𝑇 =.
𝑈 ⁄𝐿𝜌̅ 𝐿 𝐿𝜌̅ 𝑈 𝐿 𝐿

ここで， は有次元量を表す．また 𝐿(7)は 𝑧 方向における単位長さである．

3．数値計算法
3.1 移動格子有限体積法
移動格子有限体積法は，Fig. 1 のような，時間・空間の 3 次元の検査体積に有限体積法を適用することで，格子が移動・変形をしても幾何保存則を厳密に満足しながら計算を行うことが可能な手法である．𝑛 時間段階での計算セル(青色で示す)が ∆𝑡 後の 𝑛 + 1 段階で別の場所(赤色で示す)に移動したとする．格子点の時間を含んだ位置ベクトルを 𝑹 =(𝑥, 𝑦, 𝑡) とし，右下の添え字 𝑖, 𝑗 が格子点番号，右上の添え字 𝑛, 𝑛 + 1 が時間段階を表す．以下， は時間と空間の 3 次元領域を示す．移動格子有限体積法の定式化の基礎となる検査体積 𝛺 は，3 次元空間 (𝑥, 𝑦, 𝑡) において 𝑛 時間段階での計算セル面と 𝑛 + 1 段階の計算セル面を結んで形成される六面体である．この検査体積の側面は計算セルの 4 つの表面がそれぞれの時刻 𝑡での面と時刻 𝑡での面とを結んで形成される．この時間

𝑛 段階と 𝑛 + 1 段階の計算格子により形成される 3 次元の検査体積を用いて基礎方程式の離散化を行う．また，物理量の定義はセル中心型を採用する．

－ 131 －Fig. 1 検査体積 𝛺式(2)を検査体積 𝛺 に対して積分する．
𝜕𝑬 𝜕𝑭 𝜕𝒒++𝑑𝛺 = 0.
𝜕𝑥 𝜕𝑦 𝜕𝑡(8)𝜕 𝜕 𝜕(𝑬, 𝑭, 𝒒) d𝛺 = 0.
, ,𝜕𝑥 𝜕𝑦 𝜕𝑡(9)式(8)を発散形で表すと，3 次元空間に対してガウスの発散定理を適用し，離散化する．
𝑬𝑛 + 𝑭𝑛 + 𝒒𝑛= 0.

(10)

ここで，𝑛 , 𝑛 , 𝑛 は Fig. 2 に示すような検査体積における外向き単位法線ベクトルの𝑥, 𝑦, 𝑡 方向成分を表す．

Fig. 2 検査体積 𝛺 の法線ベクトル－ 132 －検査表面 𝑙 = 5, 6 は時間方向にのみベクトルを持つため，𝑛 = 𝑛 = 0 であり，𝑛 は計算セルの面積に相当する．𝑙 = 1~6の各検査表面の流束ベクトル 𝑬, 𝑭, 𝒒 の評価は，各面の時間段階に合わせて評価を行う．𝑙 = 1~4 では 𝑛 + 1⁄2 時間段階，𝑙 = 5では 𝑛 時間段階，
𝑙 = 6 では 𝑛 + 1 時間段階である．これを考慮すると，式(10)は次式となる．
𝒒(𝑛 ) + 𝒒 (𝑛 ) +𝑬

⁄

𝑛 +𝑭⁄

𝑛 +𝒒⁄

𝑛= 0.

(11)

これが移動格子有限体積法の基礎式となる．
3.2 物体透過格子法
Fig. 3 に物体透過格子法による物体の移動方法を記載する．

Fig. 3 物体透過格子法による物体の移動方法
まず，Step1 のように物体周りに格子が存在する場合を考える．物体が右に移動すると，Step2 のように物体と物体の右側の格子間隔が狭くなり，左側では格子間隔が広くなる．
そこで，物体の右側にある格子線を削除し，物体の左側に格子線を追加する．同時に物体の上下部分で格子連結線の付け替えを行うと Step3 のようになる．以上より本手法を用いることで，計算領域内を自由に物体が動くことが可能になる．物体透過格子法の定式化に関しては参考文献[5]に譲る．
3.3 流束評価法と物体・流体の連成手法
移動格子有限体積法で基礎方程式を離散化すると，幾何保存則[4]を満足するのに必要な移動格子項 𝒒⁄𝑛 が現れる．移動格子項と移流項の評価には，QUICK 法を，粘性項と圧力項の評価には中心差分型離散近似を用いた．流れに対する基礎方程式（連続の式とNavier-Stokes 方程式）の解法には SMAC 法[7]を用いた．SMAC 法における第 1 段階の解法には LU-SGS 法[8]，第 2 段階の圧力修正方程式の解法には Bi-CGstab 法[9]を用い，－ 133 －

前処理として不完全 LU 分解法[10]を用いた．ところで，SMAC 法をコロケーション格子に対して行うと，通常，圧力修正方程式の解は振動するが，本論文では，Rhie and Chowが提案した方法[11]により振動を回避した．
物体・流体の連成手法は SMAC 法の第 1 段階のときに行うことにする．計算手順としては Fig. 4 のようなフローチャートのように行った．

Fig. 4 物体・流体の連成手法に対するフローチャート

4．タンディッシュ内の流れのシミュレーション
4.1 計算モデルと計算条件
Fig. 5 に介在物が 2 つ含まれるタンディッシュの計算モデルを示す．

Fig. 5 タンディッシュの計算モデル－ 134 －

ここで，タンディッシュの大きさは 𝐿 = 3000 mm，𝐻 = 1000 mm あり，参考文献[12]のものに近い値を設定した．流入口の大きさは 𝐿 = 100 mm，流出口の大きさは 𝐿 =100 mm とし，流出口と壁からの距離は 𝐿 = 150 mm とした．また介在物の直径は100 mm とした．2 次元計算のため奥行きは単位長さである．流入速度は 𝑈 = 3.04 ×10 m/s としている．これは重量流量に換算すると 0.1t/min に相当する．参考文献[12]によれば，重量流量は 4~18t/min である．しかしながら，本計算では 2 次元計算であることを考慮すると大きなレイノルズ数で計算しても乱流は表現できず，むしろ細かい流れによって大局的な流れが捉えられないことを想定したため，重量流量を本来の 1/10 とした．物体の密度，流体の密度はともに鉄の融点での密度 6980kg/m とした．流体の動粘性係数は 0.7 × 10m /s とした．レイノルズ数 Re は，代表長さを介在物の直径，代表速度を流入速度とし計算すると Re = 4340 である．
次に，介在物の位置について，介在物①を (𝑥, 𝑦) = (1050 mm, 600 mm) に，介在物②を (𝑥, 𝑦) = (1950 mm, 600 mm) の位置に設定した．堰は幅が 100 mm，高さが 200 mmであり，タンディッシュの下中央に設置した．計算は堰があるケースと堰がないケースの2 ケースを計算した．
次に計算条件について格子点数は 451 × 151 とした．Fig. 6 に初期格子を示す．

Fig. 6 計算モデルの計算格子計算領域内に等間隔の構造格子を形成した．なお介在物周りの格子は楕円形格子生成法を用いて格子を生成した．
初期条件については，圧力はゲージ圧として 0Pa，速度は 𝑢 = 𝑣 = 0 m/s とした．一定
時間物体を固定した状態で，流入口から流速を与えた後に物体の拘束をはずした．境界条
件について，流入面は圧力ノイマン型境界，速度固定 (𝑢 = 0, 𝑣 = 𝑈)，流出面は圧力固定(𝑝 = 0.0)，速度は線形外挿，壁面は，圧力ノイマン型境界，速度固定 (𝑢 = 𝑣 = 0)，移動物
体は圧力：ノイマン型境界，速度：物体の速度とする．計算は 𝑡 = 33.4 s まで行った．時間刻み幅は代表時間の 1/100 とした．

－ 135 －

4.2 計算結果
Fig. 7 にそれぞれの介在物の軌跡を示す．

Fig. 7 介在物の軌跡 (左：堰がある場合，右：堰がない場合)
ここで，Body1 は介在物①，Body2 は介在物②を示す．また図中の点は 3.04s ごとに打っている．点と点との間隔が広いほど，介在物の速度が速いことになる．
まず， ①の介在物について，はじめは堰がある場合，ない場合ともに同じような軌跡を辿ることがわかる．しかし後半では堰がある方が，ない方よりも早い段階で介在物が高い位置に来ている．これより，介在物の浮遊という点では堰がある方が早く浮遊すると言える．また，堰のない場合における①の介在物は， 𝑥座標の1.5 m の位置付近を通過したあと，一度戻り，再度進む挙動をみせている．これは①の介在物付近の流れ場が複雑になっていることを示唆している．一方，②の介在物については，①の介在物と同様に，はじめは堰がある場合，ない場合ともに同じような軌跡を辿る．しかしながら，堰がある方が，介在物の移動速度が速くなっていることがわかる．
次に Fig. 8 に介在物の速さの時間履歴を示す．

Fig. 8 介在物の速さの時間履歴 (左：堰がある場合，右：堰がない場合)
ここで，Body1 は介在物①，Body2 は介在物②を示す．また介在物の速さは物体の速度の絶対値である．

－ 136 －

速度の時間履歴の傾向は堰がある場合，ない場合ともに同じである．𝑡 = 5 s から 𝑡 =10 s まで一度加速し，𝑡 = 10 s から 𝑡 = 15 s まで減速，その後加速する様子が伺える．
しかし，2 回目の加速以降，傾向は似ているものの，速度の絶対値については堰のない場合の方が大きい．このことから，速度の絶対値の小さい方が，介在物の挙動の制御が容易であると言えるのであれば，堰のある方が介在物の挙動制御は堰のない場合に比べ容易であると言える．また，速さが小さい方が介在物がタンディッシュに衝突した場合，ダメージが小さいと推測できるので，堰がある方がタンディッシュの寿命という点では有利であると言える．
次に Fig. 9 に 𝑡 = 6.08 s，𝑡 = 12.16 s，𝑡 = 21.28 s，𝑡 = 27.36 s における渦度分布を示す．

Fig. 9 𝑡 = 6.08 s，𝑡 = 12.16 s，𝑡 = 21.28 s，𝑡 = 27.36 s における渦度分布(左：堰がある場合，右：堰がない場合)
赤色が反時計回りの渦，青色が時計回りの渦を示す．まず 𝑡 = 6.08 s において，物体周りの流れに差異はなかった．𝑡 = 12.16 s については下中央付近で上向きの流れが確認できる．ただし，上向きの流れの発生原因は堰がある場合とない場合で異なる．堰がある場合，－ 137 －

堰に流体がぶつかることによって上に向きの流れが発生している．
一方，堰がない場合 𝑡 =6.08 s の左下にある渦が壁に沿って移流するが，その渦が剥離して上に移動する．その移動によって上向きの流れが発生している．𝑡 = 21.28 s の時には，堰がある場合，ない場合の両方において，上向きの流れが壁に当たっている．その衝突によって上中央付近の壁では左右に流れが発生していることが分かる．𝑡 = 27.36 s においては堰がある場合は，上向きの流れに沿って介在物①が上に移動している．一方，堰がない場合，上向きの流れには乗らず，下中央付近の青色の渦，すなわち，時計回りの渦に介在物①が巻き込まれることがわかる．これより，堰がない場合の介在物①は一度 𝑥 の負の方向に移動する．この移動が介在物①の浮遊するタイミングが遅くなっている理由となる．
以上より，タンディッシュ内の流れを可視化することによって，堰の有無と介在物の動きの関係を確認した．
5．結論
本研究では上記の物体透過格子法を用いて流動解析を行うが，特にタンディッシュ内の堰の有無と介在物の動きについて調べた．その結果，介在物の軌跡，介在物の速度履歴，及び，各時間におけるタンディッシュ内の流れを確認することができ，介在物の浮遊するタイミングのしくみを知ることができた．
謝辞
本研究は，公益財団法人 JFE21 世紀財団の 2018 年度技術研究助成による支援を受けて実施した．ここに深く感謝の意を表する．
参考文献[1] Steger, J. L., Dougherty, F. C. and Benek, J. A.: A Chimera Grid Scheme; Advancesin Grid Generation, American Society of Mechanical Engineers Fluids Engineering
Division, Vol. 5, pp. 55–70 (1983).
[2] Saiki, E. M. and Bringen, S.: Numerical Simulation of a Cylinder in Uniform Flow:
Application of a Virtual Boundary Method, Journal of Computational Physics, Vol.
123, No. 2, pp. 450–465 (1996).
[3] Ye, T., Mittal, R., Udaykumar, H. S. and Shyy, W.: An Accurate Cartesian Grid
Method for Viscous Incompressible Flows with Complex Immersed Boundaries,
Journal of Computational Physics, Vol. 159, pp. 209–240 (1999).
[4] Obayashi, S.: Freestream capturing for moving coordinates in three dimensions,
AIAA Journal, Vol. 30, No. 4, pp. 1125–1128 (1992).
[5] Asao, S., Matsuno, K. and Yamakawa, M.: Trans-Mesh Method and Its Application－ 138 －

to Simulations of Incompressible Fluid-Rigid Bodies Interaction, Journal ofComputational Science and Technology, Vol. 5, No. 3, pp. 163– 174 (2011).
[6] 三原清孝, 松野謙一, 里深信行：移動格子有限体積法 (第 1 報, 基礎的定式化と検証),
日本機械学会論文集 B 編 , Vol. 65, No. 637, pp. 1279–1285 (1999).
[7] Amsden, A. and Harlow, F.: A Simplified MAC Technique for Incompressible Fluid
Flow Calculations, Journal of Computational Physics, Vol. 6, pp. 322–325 (1970).
[8] Yoon, S. and Jameson., A.: Lower-Upper Symmetric-Gauss-Seidel Method for theEuler and Navier-Stokes Equations, AIAA Journal, Vol. 26, pp. 1025–1026 (1988).
[9] van der Vorst, H. A.: Bi-CGSTAB: A Fast and Smoothly Converging Variant of BiCG for the Solution of Nonsymmetric Linear Systems, SIAM Journal on ScientificComputing, Vol. 13, No. 2, pp. 631–644 (1992)[10] 村田健郎, 名取亮, 唐木幸比古：大型数値シミュレーション, pp. 51–153, 岩波書店
(1990).
[11] Rhie, C. H. and Chow, W. L.: Numerical Study of the Turbulent Flow Past an Airfoilwith Trailing Edge Separation, AIAA Journal, Vol. 21, No. 11, pp. 1525–1532 (1983).
[12] 古米 孝平, 荒牧 則親, 三木 祐司, 村井 剛: 高清浄度鋼錶片の製造方法及びタンディッシュ 特願 2013-190799.

－ 139 －