第 2 章 渦法と粒子追跡法による混相流解析手法と数値モデル
2.5 相間相互作用の計算方法
2.5.2 粒子の並進運動による流れ場の変化の計算方法
(1) 粒子の並進運動に対するTwo-wayモデル
本研究で用いる渦法と粒子追跡法によるラグランジュ・ラグランジュ解析に限らず，
オイラー・ラグランジュ型の解析においても，粒子の並進運動に作用する流体力の反作 用として，粒子が流体へ与える影響を流体の運動方程式（運動量の保存式）に導入する 方法が提案されている．ラグランジュ・ラグランジュ解析に対するTwo-wayモデルでは，
個々の粒子が受ける流体力の反作用から，流れ場解析における渦要素の渦度変化を時々 刻々算出する計算方法が用いられている．
Joia ら[18]は，液滴流に対する Tang ら[85]の渦法解析に修正を加えて，粒子による渦要 素の渦度変化を考慮した二次元解法を提案し，固気二相噴流を解析した．ただし，渦度 変化の計算法にいくつか不明な点があり，実験結果等との比較検証がないため，解法の 妥当性が十分検討されていなかった．
これに対して，内山ら[15] - [17]は，渦法を用いた固気二相自由乱流のTwo-way解法を提 案している．流れ場を計算格子に分割し，各格子において粒子運動に起因する渦要素の 渦度変化を，面積や体積による配分方法によって求める Two-way モデルを用いている．
固気二相混合層の解析に適用し，気相組織渦と粒子運動との関係や，粒子分布・速度分 布等を求め，実験結果と比較し，高精度な解析ができることを示している．
本研究では，上記のような外部流れの固気混相流を対象に内山ら[15] - [17], [20], [21]が提案
したTwo-wayモデルを用いる．流体相および粒子相ともにラグランジュ解析を用いると
Two-way couplingの計算の際に，離散点数のべき乗で演算量が増大する問題がある．そ
のため，上記研究例で行われているように，オイラー的な計算格子を援用することで計 算負荷を低減させることを試みた．
(2) 計算方法（二次元）
本研究では，Two-wayモデルにオイラー的な計算格子を援用する．この手法では，先 ず初めに，ラグランジュ的に計算された個々の粒子が受ける流体力 fp（抗力（定常力お よび非定常力），揚力等）を，単位体積あたりの力Fpとしてそれぞれ計算格子に配分す る．次に，格子点間でFpが線形変化するものと仮定することで微分演算を行い，二相間 での運動量授受を算出し，流れ場解析における渦要素の渦度変化として時々刻々フィー ドバックする．以下に計算方法の詳細を記す．
粒子が流体へ与える外力を考慮した非圧縮粘性流れにおける流体相の質量および運動
量の保存式は以下のようになる．ただし，粒子－流体の相互作用が考慮されても，粒子 体積率が十分低い仮定のもとで，流体相の質量保存式には粒子の混入による質量変化を 考慮していない．
0
u (2.5.1)
p f f
t u u p u F
u 1 2 1
(2.5.2) ここで，Fpは単位体積あたりの流体が粒子から受ける力であり，粒子に働く流体力（並 進運動に作用する力）の反作用力として運動方程式に導入される．
上式から，粒子の並進運動による相互作用を考慮した渦度輸送方程式は以下のように なる．
p
t u ω ω u ω f F
ω 2 1
(2.5.3) 粒子との相互作用（運動量の授受）による渦度の変化を考慮するため，流れ場に固定 した任意の閉曲線まわりの循環 の時間変化を考える．
ω dS Dt D dt
d (2.5.4)
ここで，dSは面積ベクトルである．
二次元の渦度輸送方程式において粘性拡散項を無視した式を，上式に代入し，Stokes の定理を用いて変形すれば次式を得る．
r F S
F d d
dt d
p f p
f
1
1 (2.5.5)
drは線素ベクトルである．
ここで，オイラー的手法を援用するため，解析領域を任意の四角形格子に分割する．
一つの格子を図 2.5.2 に示す．格子点 1～4 における Fpの値Fp （ =1～4）に対し，x およびy 方向成分Fpx ，Fpy は，次のように格子点に配分する．図 2.5.2に示した一つ の格子に含まれる粒子数をNp，粒子jに作用する流体力をfp jとすれば，四つの格子点に おける値Fp は面積配分によって，次式で求められる．
j p N
j j p
p
y x A y
x f
F
1
1 （ =1, 2, 3, 4） (2.5.6)
ただし，A jは格子点 を見込む四角形の面積であり，A1j A2j A3j A4j x yであ る．なお，各格子点は複数の格子に共有され，解析領域内部では四つの格子に関係する．
このため，各格子点におけるFp の値は関係する格子に上式を適用して得られる値の総 和として与える．
次に，格子点間で Fpが線形変化するものと仮定すれば，この格子における の時間 変化率 tは，上式から導かれる次式で与えられる．
1 4 3 2
4 3 2 1
2 2 1
py py py py
px px px px f
F F F y F
F F F x F
t (2.5.7)
格子に Nv個の渦要素が存在する場合，時間間隔 tにおける渦要素一つ当たりの循環 の変化量を Nvとする．渦要素が存在しない場合には，循環 の渦要素一つを格子 中央から新しく発生させる．
(3) 計算方法（三次元）
二次元での計算方法と同様に，粒子の並進運動による相互作用（運動量の授受）を考 慮した渦度輸送方程式は，以下のようになる．
p
t u ω ω u ω f F
ω 2 1
(2.5.8)
ここで，Fpは単位体積あたりの流体が粒子から受ける力であり，粒子に働く流体力（並 進運動に作用する力）の反作用力である．
粒子と流体との相間相互作用による渦度の変化を考慮するため，流れ場の任意の体積 における渦強度αの時間変化率は，Reynoldsの輸送定理と流れ場の連続の式から得られ る式に，伸縮項と粘性拡散項を無視した渦度輸送方程式を代入すれば，以下のように求 められる．
Grid(i) : 1-2-3-4
fpi
A3 j
Fp4
Fp1
1
4 3
2 A2 j
Δx Δy
Fp2
Fp3
A4 j
A1 j
Solid particle
Vortex element Grid(i) : 1-2-3-4
fpi
A3 j
Fp4
Fp1
1
4 3
2 A2 j
Δx Δy
Fp2
Fp3
A4 j
A1 j
Solid particle
Vortex element
図2.5.2. 粒子の並進運動に対するTwo-wayモデル（二次元）
S F n V
α F
d dt d
d
p f
p f
1
1 (2.5.9)
ここで，nは外向き単位法線ベクトルである．
解析領域を図2.5.3 (a)に示すような直方体格子に分割する．格子の各面においてFpの 値が既知であれば，格子内における渦強度αの時間変化率 α tは上式から求められる．
例えば，x方向成分 x tは次式で表される．
z D py C py y A pz B pz f
x F F S F F S
t
1 (2.5.10)
ここで，FpzAは面AにおけるFpzの値， Syはy軸に垂直な格子面の面積 x zを表す．
また，FpyCは面CにおけるFpyの値， Szはz軸に垂直な格子面の面積 x yを表す．
格子に Nv個の渦要素が存在する場合，時間間隔 tにおける渦要素一つ当たりの渦強 度の変化量を α Nv とする．渦要素が存在しない場合には，強度 αの渦要素一つを格
図2.5.3. 粒子の並進運動に対するTwo-wayモデル（三次元）
(b) 格子点上での流体力の算出 (a) 格子の概略
Fp1
Fp2
Fp3
Fp4
Fp5
Fp7
Fp8
fp j 1
2
3 4
5
7 8
6
x y
z
A
B
C D
x y z
fp j
Fp6
子中央から発生させる．
一方，渦強度の時間変化率 α tの算出に用いた格子面における Fpの値（たとえば FpzA）は，面を構成する四つの格子点でのFpの平均値として求める．格子に含まれる粒 子数をNp，粒子jに作用する流体抗力をfp jとし，二次元の解法を拡張すれば，格子点 における値Fp （ =1～8）次式で求められる．
j p N
j
j p
p
z y x
V z
y
x f
F
1
1 （ =1, 2, …, 8） (2.5.11)
ここで，Vは格子体積 x y z，V jは格子点 を見込む直方体の体積である（図2.5.3 (b) 参照）．なお，各格子点は複数の格子に共有され，解析領域内部では八つの格子に関 係する．このため，Fp の値は関係する格子に上式を適用して得られる値の総和として 与える．
以上のように，粒子が流体に及ぼす影響は，粒子に作用する流体力の反作用から，渦 度の変化として求める．