固気二相流の離散粒子シミュレーション
田中敏嗣（大阪大学大学院工学研究科 機械システム工学専攻）
１．はじめに
周知のように、流体単一相のみの流れ（単相流）の問題は数値シミュレーションが高度に発達している分野である。これに対して、固相、気相、液相のうち複数の相が混ざった流れ（混相流）は、幅広いスケールの内部構造をもち、数値シミュレーションの基礎となる数理モデルや計算の手法などに関しても解決されなければならない問題が多く残されている。筆者らのグループでは、気流中に固体粒子を含む固気二相流の問題に対して、個々の粒子運動に着目した粒子運動モデルである離散粒子モデルを用いて数値シミュレーションを行う手法やモデルなどに関する研究を行ってきた。また、このような手法により流動層内で形成される気泡や、粒子分散流中で形成される粒子クラスターなどのメゾスケール構造形成の予測と制御に関する研究を行っている[1]。ここでは、これまでサイバーメディアセンターの計算機を利用して行った研究のうちから、粒子クラスターの3次元構造形成の数値シミュレーションに対して、希薄気体運動の計算手法であるDSMC (Direct Simulation Monte-Carlo)法[2]を応用することにより計算時間の短縮を行ったことについて述べることにする。
２．粒子間相互作用モデル
離散粒子モデルは、基本的には、固体、液体、気体などの問題を分子運動レベルで扱う分子運動モデルと同様のものである。離散粒子モデルにも分子運動モデルと同様に
Hard sphere モデルと Soft sphere モデルが存在し、Hard sphere モデルにより衝突支配運動を、Soft
sphereモデルにより多体接触運動が記述できる（図1）。離散粒子モデルが分子運動モデルと本質的に異なる点は、粒子間の近接相互作用である衝突や接触がエネルギー散逸を伴う点にある。衝突や接触などの相互作用を伴う粒子群の運動を計算機により数値解析する場合に、粒子間衝突や接触を検出する手続きは非常に大きな計算負荷となる。ある粒子と別の粒子との間で衝突や接触が起こるかどうかは、これらの粒子の相対運動に対して幾何学的な判定を行えばよい。N個の粒子の運動を扱うとすると、単純にすべての粒子の組み合わせに対して判定するとすれば、その数はN(N-1)/2、つまりN2に比例するものとなり、Nが大きくなると計算負荷が非常に増大する。
この問題を回避するのは以下に示すように容易である。図2に示すように計算領域をセルに分割し、何番の粒子がどのセルに含まれているかを帳簿に記帳するのである。粒子間衝突を扱う場合、セルの寸法を1回の時間進行での粒子間の最大相対移動量以下に選んでおけば、あるセル(I,J)に含まれる粒子が衝突可能な粒子は、セル(I,J)とそれに隣接セル内の粒子だけである。図2の2次元の問題の例では、中央のセル(I,J)に含まれる粒子に対する衝突判定領域は、ハッチングを施した9個のセルである。このような工夫により、必要な衝突判定の回数はNに比例したものとすることができる。３．ベクトル機における計算の問題点
筆者らのグループでは、離散粒子モデルを用いて、粒子運動と流体運動が強い相互作用を行う流れの計算を行っている。流体運動解析に対するベクトル機の利点を活用するためSX5を用いているのであるが、その際、問題となるのが粒子間衝突の計算に要する計算時間である。前節で述べた粒子間衝突の計算手法では、衝突判定領域内に存在するすべての粒子との間で決定論的に衝突判定を行い、衝突判定という条件判断に基づいて衝突する場合には衝突の手続きを行う。そのため、計算負荷の大きな粒子間衝突の処理に関してベクトル化率を上げることができず、粒子間衝突に要する時間が計算時間の大きな部分を占めるようになっていた。ある例では、計算時間の91.7%を粒子間衝突計算が占めるものとなっていた。計算領域の大規模化、粒子濃度の高濃度条件への展開を行うためには、粒子間衝突計算の高速化は避けては通れない問題となった。
４．DSMC法の導入と計算の高速化
分子運動モデルによる希薄気体運動の計算手法の一つであるDSMC法では、分子間衝突確率に基づく確率論的な手法により、分子運動を計算する。つまり、ある分子が1回の時間進行中に衝突するかどうか、衝突するとすればどの分子とどのような衝突をするかが、乱数のサイコロを振るモンテカルロ法により取り扱われる。それに加えて、実在の膨大な数の分子から抽出したと標本と見なされる少数の分子運動のみを計算することにより計算負荷が低減される。筆者らの研究では、その粒子濃度分布の計算結果が図3に示されるような、粒子濃度が大きく変化する流れを取り扱い、粒子濃度が非常に希薄な領域も存在する。このため、計算条件に対して計算粒子数の低減は行わず、衝突判定の手続きだけにDSMC法を応用することにした。△t時間の時間進行あたりの粒子iの衝突確率Piは次式で与えられる。
図3は粒子間衝突計算にDSMC法を応用して計算を行った粒子クラスターの3次元計算の結果である[5]。計算はすべての境界を周期境界とする直方体状の計算領域内で、重力の作用下での粒子運動と気流の運動を求めた。約149万個の粒子運動が取り扱われており、粒子の体積分率は1×10-3である。図は一つの鉛直面内における粒子濃度の分布を示している。この図に示されるように、粒子濃度は一様にはならず、空間的に大きな非一様性を生じ、時間的および空間的に大きな乱れを生み出す。
図3の計算結果は粒子間衝突に決定論的手法を用いた結果と比較され、よく一致することが確認されている。また、検証実験が現在も引き続き行われており、空間パターンの特徴がよく一致することが確認されている。
粒子間衝突に関する計算時間については、DSMC法の採用により決定論的手法を用いるものに比べて13.6%まで短縮できた。DSMC法の採用以外にも、ベクトル化率が必ずしも高くなかった粒子運動計算のベクトル化率を上げる努力をした結果、ベクトル化率を99.49%まで上げることができた。
５．おわりに
DSMC法の採用により、これまで全計算時間の約92%を占めていた粒子間衝突の計算時間は大幅に短縮し、全計算時間は11.5%とすることができた。本研究では、大阪大学工学研究科機械物理工学専攻修士課程の茶園史也君の多大なる協力を得た。記して謝意を表する。
文 献
[1]Bird, G. A., Molecular Gas Dynamics, Oxford Univ. Press, (1976), 118.[2]田中敏嗣,ながれ, Vol.21, No.3 (2002), 250.
[3]Nanbu, K.,J. Phys. Soc. Japan, Vol.49, No.5 (1980), 2042.
[4]Illner, R., and Neunzert, H., Transport Theory and Statistical Phys., Vol.16 (1987), 141.
[5]茶園史也,田中敏嗣,辻 裕, 日本機械学会関西支部第78期定時総会講演会講演論文集, (2003), 14-27.