第 27 回数値流体力学シンポジウム
E01-2

DSMC 計算でセルを大きくできる分子間衝突法 U-system のその後
（テーパ円柱後流 3 次元解析と円管内乱流速度分布再現に向けて）
Subsequent development of the intermolecular collision scheme, U-system, for DSMC method○ 宇佐美 勝, 名城大, 〒468-8502 名古屋市天白区塩釜口 1-501, E-mail: usamim@meijo-u.ac.jp
野村 将志，村木 卓馬， 名城大院，〒468-8502 名古屋市天白区塩釜口 1-501
Masaru Usami, Meijo University, 1-501 Shiogamaguchi, Tempaku-ku, Nagoya 468-8502, Japan
Masashi Nomura and Takuma Muraki, Meijo University
The direct simulation Monte Carlo (DSMC) method(1)(2) has been developed for a numerical method in a rarefied gas
flow, and now has been expanding its feasibility to continuum flow region. Although Meiburg(3) has pointed out that itis difficult for the DSMC method to deal with flows including vortices since a collision pair of molecules is selected ina cell in disregard of their locations, the DSMC method using the new collision scheme (U-system)(4) has been verified
to be effective for flows with periodic vortex shedding at the RGD24(5) and the RGD26(6). One problem to resolve
associated with U-system was such that the momentum and the energy accompanied with a collision pair were notconserved simultaneously during U-system collision process. Sun et al. (7) have found that the correction scheme
evolved by Pareschi and Trazzi(8) for momentum and energy is useful for the conservation process in U-system. The
purpose of this study is to verify Pareschi's correction scheme and to calculate the three dimensional flow field behinda circular taper cylinder. U-system is also applied to a velocity profile in a circular tube under the condition of aturbulent flow even though it is still at a testing stage.
２．新分子間衝突計算法 (新衝突法U-system)(4)(5)(6)(11)(12)
従来の衝突計算法では，セルが大きくなるとそれに比例して二つの衝突分子の位置は離れるため，シミュレーションが忠実に行われなくなる．本来，衝突する２つの分子は同じ位置で選ばれるのが最も好ましいが，有限個の分子しか用いないシミュレーショ
ンでは，たとえセルを小さくしたとしても全く同じ位置で２つの分子を見つけるのは困難である．従来のDSMC法では，この問題を少しでも解決するために，セルの大きさを平均自由行程以下に保つことを必要な条件とした．さて，衝突する２つの分子の位置座標を一致させるため，片方はシミュレーション上の分子から実際に選び，もう片方の分子は，その速度を，1 つ目の分子の位置の速度分布関数から仮想的に作り出すという方法が考えられる．
ただし，この方法は，任意の場所の速度分布関数を正確に表現する必要があるし，仮想的に作り出した分子は衝突後に削除される
ので，運動量やエネルギの保存が困難になり，得られる結果に不規則で大きな変動が含まれてしまう．その不利益を極力小さくするために，本新衝突法は，なるべく近くに存在する分子の速度を修正して仮想分子を作り出すという考えに立ち，これを発展させた(宇佐美(11)の補遺参照)．
最終的に考案された新衝突法では，衝突計算直前に一方の分子Pを，あたかも，もう一方の分子Qの位置に一致させたかのように分子速度の補正を行う．いま，流れ場中の分子は，全て，ある温度と流れ速度を持つ局所平衡状態の速度分布関数（Maxwell分布）に従うものと仮定する．なお，以下では，２つの分子位置における速度分布関数をマクスウェル分布で仮定して補正式を導いているが，得られた補正式自体は，２つの位置の分布関数が必ずしも
マクスウェル分布でなくとも有効である．唯一必要な条件は，「分
布中心が，流速に応じて左右に平行移動し，分子の熱運動速度が，温度の平方根に比例するような分布形状」というものである．大
きいセルであるといっても同じセル内にある２つの分子であるから，マクスウェル分布でないとしても，上記の条件を満たすのはそれほど困難なことではないと考える．さて，分子Pの位置の分布関数をf1，分子Qの位置の分布関数をf2とし，分子Pの速度がf1上の一点Aにあるものとしよう(Fig. 1)．衝突計算に際して，このA
点を分布関数f2上のどこかに移せばよいことになるが，最も自然な方法は，f1上のA点位置とf2上の位置が相対的に同じになるよう移動させることである(B点)．

１．緒言
直接シミュレーションモンテカルロ法(DSMC法)(1)(2) は希薄気体流の解析法として開発されたが，計算機技術の発展に伴って，その適用範囲は，連続流状態まで広がる傾向にある．
DSMC法は，計算の発散といった不安定要素がなく，また，境界条件も比較的与え易い．一方，DSMC法で最も問題となるのは，流れ場のセル
分割を，平均自由行程と同程度の細かさで行う必要があるという点である．従来の衝突計算法では，同一セルに含まれる任意の分子２個は，その相対位置に関係なく衝突ペアとして選択されるた
め，セルを大きくすることはできなかった．
物体の後流に生じる渦とその剥離流れは，流体力学において古くから存在する問題であるが，現在でもまだ解析が続けられている(9)．Meiburg(3) は，斜め平板後流に生じる渦の剥離流れに関連させて，DSMC法の分子間衝突に関する欠点を指摘している．Bird(10)は，その指摘に対して強く反論しているが，渦が関与する解析では，セル分割を平均自由行程より十分小さく（例えばその 1/3 以下に）する必要性を特に強調せざるを得なかった．
著者らは，衝突する分子の位置に応じて分子速度の補正を行った上で分子間衝突を行うという算法 (U-system) により，DSMC計算においてセル長が大きくとも良好な結果を得る方法を考案し，様々な問題に適用してその手法の改善を行っている(4)(5)(6)(11)(12)．今
回は，Pareschi(8)らが提案した「分子の運動量と並進エネルギをセルごとに厳密保存する手法」をU-systemに適用することを試み，それにより，直径が変化するテーパ円柱の後方に生じる渦を解析したものである．直径一定の円柱の後方に生じる渦の解析はすでに多くなされており，豊富な知見が得られているが，円柱の直径が変化するテーパ円柱については，その後方にどのような渦が発生し，また直径一定の円柱と場合と，結果にどのような差異が生じるのか等について，まだ調査すべき点が多く残されている．
DSMC法を用いることにより，この種の渦を含んだ高速流れや低密度流れの解析が，今後，比較的容易に行えるようになると考える．また，純粋なBirdの衝突計算法では再現できないとされている円管内乱流速度分布について，U-systemを適用した場合の結果について示し，分子シミュレーションによる乱流再現の可能性について検討した．
1

Copyright © 2013 by JSFM

第 27 回数値流体力学シンポジウム
E01-2
時間は，従来の計算法より増加するが，セルを大きくできる点を考慮すれば，計算全体での時間比較では新衝突法が有利となる．
さて，分子速度の補正を行うと，各衝突分子ペアの持つ運動量とエネルギが衝突前後で厳密には保存されない．もちろん流れ場
全体としてあるいは長時間にわたって平均をとれば，統計的ばらつきの範囲内で保存されるのであるが，多くの計算例を検討した
結果，一回の衝突において，衝突前後のエネルギに大きな差が生
じると，特に衝撃波が生じるような変化の激しい流れ場では，ひずみを成長させる原因となることが判明した．そこで，衝突直前の２つの分子が持つエネルギ和と衝突直後のエネルギ和を求め，その比の平方根を用いて衝突後の分子速度成分を修正し，衝突前
後でエネルギが保存されるような改良を考えた．例えば，同じ質量の２つの分子P, Qの衝突前の速度成分をcPi , cQi ( i = 1~3) とし，衝突後の値をcPi*, cQi* ( i = 1~3) としたとき ( cP1 = uA, cP1* = uD )，式(7) で得られるHを，衝突後の両分子の各速度成分cPi*, cQi* に掛算して，最終的に修正された衝突後の分子速度とした．

Fig. 1 Modification of velocity in the new collision scheme.
このB点位置を求める計算は簡単である．いま，x方向速度成分uの分布関数について考えると，f1とf2は次のように書ける．

ここで，Rは気体定数，Tx1, Tx2およびU1, U2は，x方向温度および流速のx成分である．f1およびf2をそれぞれ最大値で無次元化し，それらを等しいとおくと，この方法はU-systemの計算安定化に大きく貢献したが，運動量に関する保存が保証されてないとの指摘があることと，また，計算条件によっては温度上昇が予想外に生じるといった事例も発生して，更なる改善が望まれていた．Sunら(7) は，分子間衝突を無次元速度空間で行うことを考え，速度を無次元化するためのスケール変数に流速と温度を用いるという概念を導入してDSMC解析を行っている（LCDSMC法）．その計算法は，基本的には，著者らの新衝突法U-systemと同じ算法ではあるが，特筆すべきは，分子の運動量と並進エネルギを各セルごとに厳密保存するために，
Pareschi ら(8) が提案した補正方法を採用したことである．
流速と温度が与えられた平衡状態の速度分布関数から，乱数を用いて，複数の分子の速度を再現した場合，再現された分子の速度から，逆に流速と温度を求めても，厳密には元に戻らない．これを厳密に元に戻すために，Pareschi らは，運動量とエネルギの
厳密保存の方法を考案した．いま，この方法をU-systemに適用することを考える．
衝突前と衝突後の分子速度成分をu , v , w およびu*, v*, w* とし，これを，Pareschi らの補正法を使って運動量とエ
ネルギを厳密保存して u+, v+, w+ となるものとしよう．この場合，補正は次式のように，4 つのパラメータによって行うことになる．

したがって，

この修正された速度uB（なおuAは修正前の速度）を用いて，これ以降の計算は一般的な衝突計算法（No-time-counter法(1) あるいは
Null-collision法(13)）で行う．結果として得られた衝突後の速度uC (C点) は，分子Pを元の位置に戻すのに伴い，式(5)を用いて再び修正される(D点)．

さて，式(4)(5)を用いるためには，分子PとQそれぞれの位置における流速Uと温度Tのデータが必要である．したがって，最初は，従来の計算法でセル中心の流速と温度を仮に求めておく．隣接するセル中心の値から各セルの四隅（三次元の場合は八隅）の値を導いた後，セル内の任意位置における値は，セル四隅の値を用いて双線形補間（三次元の場合なら三線形補間)で求める．例えば二次元の場合，長方形セルの四隅の座標を (x, y), (x+A, y), (x, y+B),(x+A, y+B) とし，
そこにおける温度をそれぞれ T1, T2, T3, T4 とした
とき，セル内の任意位置 (x+a, y+b) の温度Tは，式(6)により求めることができる(双線形補間法)．

(8)

ここで，
を求めるために，運動量保存式 (9) とエネルギ保存式 (10) を用いる．

(9)
DSMC 法の長い繰返し計算の中で，従来の方法による計算は最初の１回だけでよく，以降は，新衝突計算法によって流速と温度のデータが逐次更新される．なお，新衝突法の計算１回に要する
2

Copyright © 2013 by JSFM

第 27 回数値流体力学シンポジウム
E01-2
３．完全平衡法(14)
流れ場の密度が大きくなるにつれて分子間衝突の回数は増加し，例えば，大気圧状態に近づくとその計算量は膨大なものとなって
DSMC 法の適用が困難になる．一方，分子間衝突が増加すると，分子の速度分布関数は平衡状態に達し，それ以上，分子間衝突を行っても分布関数に変化は生じなくなる．もし非平衡状態を無視してもよい流れ場であるなら（(注) 一般の CFD 計算では平衡状態しか取り扱わない），
必ずしも忠実に分子間衝突の計算を行わずとも，絶えず分布関数を平衡状態のものに保ちつつ DSMC 法に則った分子シミュレーションを行って，その結果，計算負荷を減らすことができると考えられる．
DSMC 法では，同一セルの中にある分子を２つ任意に選んだ後，その相対速度に応じて衝突確率を求め，衝突可と判断されたものについて分子間衝突の計算を行う．この方法では，衝突になかなか関与できない分子がいつまでも残る可能性があり，速度分布が平衡状態に達するまでに，相当の計算量が必要となる．いま，分子の相対速度に関係なく整然と分子間衝突を行わせ，全ての分子が必ず衝突に関与するような算法を用いれば，非平衡状態の計算はできなくとも，大気圧に近い密度の流れ場で効率的な分子シミュレーションができる．
Fig. 2 は，１つのセルに存在する N 個の分子から任意に２個を選んで衝突させた後，片方の分子を除いて，残った分子群から再び分子を２個選び衝突させるという算法で，セル内の全ての分子を最低１回は衝突させ（計 N－1 回の衝突計算），この一連の操作
を 5 回繰り返して速度分布の変化の様子を調べたものである（横軸は最大確率速度で無次元化した速度）．
そこにおける最初の速度分布は，分子速度の大きさに無関係に一様な確率分布を持つものとした．計算開始時の分布関数が，一様分布といった極端なものであるにも関わらず，3 回の繰り返しで非常にきれいな平衡状態の速度分布関数が得られている．
なお，この完全平衡法の適用を試みるにあたっては，新衝突法の利用が大前提となっている．
すなわち，通常の DSMC 衝突法（従来法）であるとすれば，完全な連続流状態ではセルを極端に小さくする必要があり，
この完全平衡法がほとんど有効に機能しない．
新衝突法と組み合わせて使うことにより完全平衡法が効果的に利用できるのである．

(10)
結果として，

は，下記の一連の式で求められる．

(11)

ここで，

(12)

なお，τ を求める際に，平方根の前の符号には必ず正を用いる必要がある．負であっても保存則は満足するが，負であると，補正
前後の速度方向を一致させるために，修正用パラメータ λ が非常に大きな値になってしまう（補正値 λ と流速の絶対値がほぼ同じ
大きさになる）．一方，正の場合，λ は流速の 1/100 以下である．
結果として，平方根の前の符号に負を用いると，保存則は成立す
るが U-system の効果は殆どなくなってしまう．
DSMC法では，分子移動の計算と分子間衝突の計算が時間ステップdt で分離され交互に計算が繰り返される．本新衝突法ではセ
ルを大きくするのに伴い，平均的速度の分子がセルを 5 × dt で横切るように，すなわちdt を大きく設定して計算速度の増加を図っている．なお，そのために，分子情報のサンプリングは，分子間衝突の前後で１回ずつ行うのが有効である．なお，Sunら(7) のLCDSMC法では，時間ステップは，相変わらず平均自由時間より十分小さく取って計算しているので，解析に要する時間はかなり多く必要になると思われる．

４．複数スレッドCPUによる並列計算(14)近年の計算機技術の急速な発展は，
パソコン CPU についても同様で，筆者の経験では，計算速度と記憶容量のどちらについても2 年でほぼ 2 倍，すなわち，20 年で約 1000 倍の性能向上が実現されている．最近は計算速度に多少の足踏みが見られるが，その代わり，１つの CPU の中に複数の演算処理機能を持たせた，複数コ
ア・複数スレッドを持つ CPU が出現し，それに対応できるコンパ
イラもようやく機能するようになってきた．
DSMC 計算を最も効率よく並列化するためには，流れ場をスレ
ッド数と同じ数の領域に分割して各領域を一つのスレッドに受け持たせることである．これにより各領域の計算は，他の領域とは完全に独立する（メモリアクセス時の待機時間は別）．ただし，分子が領域の境界を横切ったなら，その分子の情報は，移動先の領域の計算に受け渡されなければならないので，その手続きを並列計算に伴って追加する必要がある．一方，別の並列計算の方法として，流れ場を分割せずに，DSMC 法の各計算ルーチンごとに並列化を実行することが考えられる．この方法では，計算を完全に独立させて並列化することは難しく，並列化の効果を十分期待できない面もあるが，その計算法を一度確立しておけば，プログラ
ムの簡単な書き換えで，どんな DSMC 計算の高速化にも適用でき
る．以下では，後者の立場に立ち，8 個あるいは 12 個のスレッド

Fig. 2 Transition to the equilibrium distribution.
3

Copyright © 2013 by JSFM

第 27 回数値流体力学シンポジウム
E01-2
ものである．ここで，U-system-5 とあるのは，今回のPareschi ら(8)の補正法を組み込む一つ前のU-system 版を用いて得た結果で，
U-system-8 とあるのが，Pareschi らの補正を加えた計算によるものである．純粋なBird法による結果が非常に鈍化した分布となっているのに比較して，U-system-5 もU-system-8 も鋭い変化を再現した結果となっている．なお，Bird法でセル長を平均自由行程の150 倍と細かくした結果も表示しているが，600 倍のU-systemの結果にはまだ及ばない．なお，Fig. 4 に示すようにU-system で 150倍にすると，さらに鋭い分布が得られる．なお，先に述べた補正
パラメータτに負を用いた場合の結果もFig. 5 に示しているが，
U-systemの効果はほとんど消えている．Fig. 6 は，同様に，温度分
布T/T∞ を比較したものであり，密度分布と同様に，U-system の効
果は明白である．なお，Pareschi らの補正の影響（U-system-5 とU-system-8 の比較）については，この噴流構造解析にはそれほど明確に現れてはいないが，他の解析の比較によって，安定性の十分な向上が認められており，今後は，U-system-8 を採用することにする．

を持つCPU（Intel Core i7）とIntel Fortran コンパイラ（Composer XE）
を用いて，DSMC 計算プログラムの各ルーチンを並列化する方法を紹介する．
4.1 並列計算プログラミングにおける一般的事項 Intel Fortran では，並列計算を実行するために，DO ループ直前に次のような記述をする．なお，PRIVATE の後のカッコ内で，並列計算で独立して用いる変数を宣言し，REDUCTION の後のカッコ内で，並列計算終了時にスレッド間のデータ集計を行いたい変数を宣言する．
!$OMP PARALLEL DO PRIVATE(…, …, …) REDUCTION(+: …,…, …)
4.2 一様乱数の作成 DSMC計算の並列プログラミングにおいては，8 個（12 個）のスレッドで並列計算する際に，独立した 8 個（12 個）の一様乱数を使う必要がある．また，プログラムの中には並列化できない部分もあり，そこでも乱数を使用するので，それも考慮して 9 個（13 個）の乱数列を準備した．具体的には，周期 246の独立した 2 つの乱数列において，一周期を複数等分した位
置に初期値を定め，それらから出発する 9 個（13 個）の一様乱数列を作成して用いた．なお，一つの一様乱数を用いて並列計算することも可能ではあるが，高速化の効果は薄れる．
4.3 分子移動（Molecular movement）の並列化 この計算ルーチン
では，基本的に個々の分子は独立して空間移動するので，全分子数を 8 個（12 個）のグループに分割して並列化すれば良い．しかし，計算領域外に流出した分子を捨て去るとき，分子情報を記憶している配列を整理する計算処理が必要で，そこは単純には並列化できないので特殊な処理が必要である．
4.4 分子所属セル調査（Cell identification）の並列化 この計算ルー
チンは分子ごとに独立であり，簡単に並列化して高速化が実現できる．
4.5 分子並べ替え処理（Molecular sorting）の並列化 このルーチン
は，計算が相互依存しているため，簡単には並列化できない．以前，ベクトル計算機で高速化を図った際には，強制ベクトル化という確率的な手法を用いたが，今回のスレッド数は 8 個（12 個）に過ぎないのでそれは不可能である．しかし，DO ループ内で若干複雑な計算処理を行うことにより，高速化が大きく期待できな
いまでも，並列化できないわけではない．
4.6 分子間衝突（Intermolecular collision）の並列化 このルーチンは，セルごとに計算が独立しているので，並列化が可能で，しかも密度が高くなるぼど全計算時間に占める分子間衝突の割合が増加するので，高速化に最も寄与できる部分である．
4.7 流れ場サンプリング（Sampling of flow properties）の並列化 このルーチンは，計算規模が大きくなると全体の計算時間に占める割合が無視できなくなるので高速化したい部分であるが，分子情報とセル情報が交差して現れるので並列化は難しい．参照用配列（cross referencing array）を用いると並列化できなくはないが，高速化はあまり期待できない．

Bird 600 times
n/n∞U-system-5 600 times
n/n∞U-system-8 600 times
n/n∞Bird 150 times

５．超音速自由噴流の構造比較による U-system の効果
テーパ円柱後流の 3 次元渦構造についてのDSMC解析に入る前に，セルが粗いDSMC計算において新衝突法U-systemがどの程度の効果を発揮するかについて，超音速自由噴流の構造解析を用いて明らかにしておく．この噴流解析については，最初にU-systemを提案した際(4) に取り扱ったもので，そのときは不安定現象が現れたため速度補正の程度を制限せざるを得なかったという非常に敏感な計算対象である．上流よどみ点圧力 760 torrのアルゴン単一
気体が，直径 2.5 mmの円形オリフィスから圧力比 1/50 の背圧の下流側へ流出する際に生じる噴流構造を求める．Fig. 3 は，上流
よどみ点の密度で無次元化して表した密度分布n/n∞ で，上流よどみ点の平均自由行程の 600 倍を基本セルの一辺長として解析した

n/n∞Fig. 3 Comparison of density profiles for supersonic free jet
between various collision schemes.
4

Copyright © 2013 by JSFM

第 27 回数値流体力学シンポジウム
E01-2
６．テーパ円柱後流に生じる渦はく離流れの 3 次元 DSMC 解析
二次元流れ場の渦はく離流れ解析については，すでに，完全平衡法も含めた比較が行われており，粗いセルにおけるU-systemの効果が明確に示されている(15)．さらに，三次元流れ場のテーパ円柱後流の渦はく離流れについてもU-system-5 を用いた解析により
3 次元性を持つ渦発生の様子が再現されている(15)．今回は，
U-system-8 を用いて，直径を変化させたテーパ円柱の後方に生じる渦をDSMC法により解析したものである．
Fig. 7 に示すように，流れ場としての計算領域は，主流方向（x軸）14 mm，縦方向（y軸）6 mm，奥行方向（z軸）7.5 mm であり，
テーパ円柱の中心軸は，主流の流入境界面より 3 mm後方に置かれている．テーパ円柱の長さは，流れ場の奥行きzと同じで 7.5 mmとし，円柱の最大直径dを 2 mm，最小直径d2を最大径の 0.75 倍(d2＝1.5 mm）として解析した．直径変化の割合は一定である．なお，
y-z平面の中心でx＝11.2 mmの位置に監視点を設け，その位置での流速の時間変化を記録した．

U-system-8 150 times
n/n∞Fig. 4 Density profile for supersonic free jet by U-system-8
with cell length of 150 times mean free path.

n/n∞Fig. 5 Density profile for supersonic free jet by U-system-8
with a negative parameter λ.
Fig. 7 Computational domain of the DSMC method.
T/T∞

分子モデルはアルゴン Ar に相当する VHS 分子を使用し，総分
子数は約4 億個である．
圧力は最初1340 Pa で始め，
途中で2680 Paに変化させる．
セルの一辺の長さは最大直径 d を 60 分割したもの
で，平均自由行程 λ の 8.5 倍（途中からは 17 倍）である．セル分
割数は，x 方向に 420 個，y 方向に 180 個，z 方向に 225 個とし，計 1701 万個のセルを用いる．Ar 単一気体の初期条件は，温度 288
K，圧力 1340 Pa で流速ゼロある．レイノルズ数は Re＝190（途中から Re＝380），クヌッセン数は Kn＝0.00119（途中から Kn＝0.00060）とする．したがって，代表長さに平均直径を用いると，主流速度は U = 103 m/s (M= 0.326) となる．
計算領域を取り囲む境
界条件については，z 軸に垂直な境界では，分子は鏡面反射する
ものとし（対称境界），それ以外の境界は，主流速度 U と同じ流速の局所平衡状態（温度288 K，圧力は最初1340 Pa，途中から 2680
Pa）を仮定する．また，円柱壁面において分子は拡散反射する．
Fig. 8 は，横軸に時間経過（データサンプル回数で表示．1 回が0.85 μs に相当）
，縦軸に，監視点における流速の主流方向成分（青色）と垂直方向成分（桃色）の時間変化を示したものである．特に，主流に垂直な速度成分(桃色)に注目すると，最初，主流が 1340
Pa で流れ込んだ後，しばらくすると振動を生じ始め，時間とともに垂直振動は増加してほぼ一定の振幅に達する．その後（2866 回
目のサンプル以降．
計算開始より2.4 ms 以降），
主流圧力を2680 Paに増加すると，垂直振動の振幅は少し増加するが，必ずしも振幅一定で安定しているわけではなく，わずかに変動を繰り返して推移する．そして，時刻 6.3 ms（7392 回目のサンプル）において，衝突計算を Bird の従来法に変更すると，たちまち振動は減衰してカルマン渦は消えてしまうが，さらに，再度，U-system に戻すと振動は復活する．

Bird 600 times

T/T∞

U-system-8 600 times

T/T∞

U-system-8 150 times
Fig. 6 Comparison of temperature profiles for supersonic free
jet between various collision schemes.
5

Copyright © 2013 by JSFM

第 27 回数値流体力学シンポジウム
E01-2
Flow velocity components

m/s

Time

Fig. 8 Oscillation of velocity components of flow by the new collision scheme.

× 0.85 μs

Fig. 10 Vorticity contours around a taper circular cylinder onthree different planes.

Fig. 9 Velocity vectors (direction only) around a taper cylinder on threedifferent planes.
6

Copyright © 2013 by JSFM

第 27 回数値流体力学シンポジウム
E01-2
Fig. 9 は，Re = 380 の場合のテーパ円柱後流の異なる 3 断面（最大直径断面 z/d = 0，中心部断面 z/d = 1.9，最小直径断面 z/d = 3.75）で観測した「流れ方向を示すベクトル（流線ベクトル）」で，各軸の目盛は，円柱最大径 d の倍数で示されている．z = 0 の位置でテーパ円柱の直径は最大で，z の値が大きくなるほど直径は小さくなる．円柱背後に渦のはく離流れが得られているが，各断面の渦中心は主流方向にずれ，テーパ円柱の径が小さい断面ほど，渦はより速く下流へ流れており，径が小さい場所ほど流速の大きいことがわかる．同様な状況は，渦度等値線を描いた Fig. 10 にも表現されている（暖色が正の渦度，寒色が負の渦度）．
流れが，主流速度一定で計算領域に流入しているにもかかわらず，円柱背後において円柱軸方向に異なる流速となるためには，
低速部分から高速部分へ主流を横切るような流れが存在しなければならない（連続の法則）．Fig. 11 は，最大直径断面のベクトル図の中に，後流渦の内部の一点を通る流線を赤色で描いたものであり，Fig. 12 は，赤色の流線だけを抜き出して示したものである．
また，Fig. 13 は，同じ流線を異なる角度から眺めたものである．
流線は流跡線と異なるので，この流線が必ずしも同じ流体部分の移動した軌跡とは言えないが，渦の中を円柱が細くなる方向に回転しながら進む流線が頻繁に出現しており，渦の中を通って，流れが主流を横切る方向に移動するのであれば，連続の法則をうまく説明できる．
Fig. 14 は，同時刻において別の流線を描いたもの，
Fig. 15 は，異なる時刻を描いたものであるが，やはり，円柱が細くなる方向に渦巻いた流線が得られている．
Re = 190 の場合，以前の報告(15) では，円柱上側にできる渦と下側の渦は，最大直径の背後では離れたまま流れ去るとしていたが，今回の解析では，
流れるにつれて接近する様子が見られた．
また，
Re = 380 の場合と同様，渦内部の一点を通る流線は，Fig. 16 示す
ごとく，主流を横切って円柱が細くなる方へ進むことが頻繁に見受けられた．

Fig. 13 Stream line behind a taper cylinder (View from different direction).

Fig. 14 Different stream line behind a taper cylinder at the same time.

Fig. 15 Different stream lines behind a taper cylinder at another time.

Fig. 11 Velocity vectors and a stream line (red line).

Fig. 12 Stream line behind a taper cylinder for Re = 380.
7

Fig. 16 Stream lines behind a taper cylinder for Re = 190.
Copyright © 2013 by JSFM

第 27 回数値流体力学シンポジウム
E01-2
７．DSMC 法による円管内乱流速度分布の再現に向けてレイノルズの管内乱流の実験以来，乱流の研究は 130 年にわたって流体力学の分野をにぎわしているが，乱流モデルのような人為的な操作を加えないで管内乱流を計算機内に再現することは至難の業であり，分子シミュレーションの分野でもDSMC法では再現不可能というのが一般的な見方である．著者らも，分子運動に人為的な乱れを加えることによる流速分布の乱流化は可能である
ものの，そのような操作なしで，レイノルズ数Reがある限界値を超えて増加したとき，流速分布が放物線層流分布から乱流分布に変化するという現象は，これまでDSMC法で再現することはできなかった．例えば，以前U-system の旧版で求めた流速分布(5)は，温度と密度の極端な変化を伴っていたため，信頼できるものではなかった．今回，運動量と並進エネルギの厳密保存の実現で，新衝突法U-systemの安定性が格段に増進したのを機に，円管内流速分布（Fig. 17）のReによる変化について再び挑戦してみた．
円管断面のセル分割は，まず，一辺 4 mm の正方形を 400 × 400の正方セルで等分割した上で，これに直径 4 mm の円を内接させ（Fig. 18 参照，実際の格子数は 400 × 400）
，円の内部にある格子が実質の管内セルとなる．
円に削られていないセルは 1 辺 0.01mmの立方体であるが，円で削られたセルは，分子間衝突の計算に必要となるのでセル体積が正確に計算されている．
計算領域としては，本来，長い円筒状の計算空間を考えるべきであるが，今回は計算を単純化するために，管軸方向にはセルを一つしか配置しないという「二次元解析」で行っている．すなわ
ち，流れ方向（管軸方向）の計算空間は 0.01 mm の長さしかない（Fig. 19）
．ただし，今回の解析プログラムでは 3 次元計算にも対応できるように，軸方向の分子位置座標も取り扱っている．円管の流入出断面は周期境界とし，一方の管断面から流出した分子は反対側の断面から流入することになる．円管壁面で分子は拡散反射するので，そのままでは，時間経過とともに軸方向の運動量が失われて流れが停止していく．そこで，それを防ぐため，分子移動の計算終了時（分子間衝突の計算の開始前）に，失った総運動
量を全分子に均等に配分する．この操作は，流れが停止するのを防ぐために管断面に圧力差を設けていることに相当する．なお，
重力の影響は考えてない．
有効なセルの総数は約 12 万 5 千個，総分子数 125 万個，アルゴン単一気体を VHS 分子モデルで近似して取り扱う．
初期条件は，
平均流速が分子の最大確率速度の 0.1 倍（35 m/s）となるような層流分布，もしくは一様分布で，初期条件の違う両者によりほぼ同
じ流速分布が安定して得られるまで計算を続けた．圧力は Re =5000 の場合 344 Pa，温度は 288 K としており，平均自由行程は，344 Pa のとき 0.015 mm である．Bird 法と U-system 法とで結果の比較を行った．調査した Re は，Re＝2, 10, 50, 200, 1000, 5000（このときのクヌッセン数はKn＝9.4, 1.9, 0.38, 0.094, 0.019, 0.0038）
で，Re を変更する場合，
例えば Re = 1000 のとき 68.8 Pa といったように，圧力設定を Re に比例させる．乱流とは無縁の小さな Re を計算した理由は，クヌッセン数 Kn が 1 以上になると希薄化の効果が大きくなり壁面で滑り現象を生じるので，それと乱流との違いを明確にしたかったのと，その条件では，Bird 法と U-system で差
異のないことを確認したかったためである．

Fig. 18 Dividing method into a network of cells for across-section of a circular tube.
0.01

4

Fig. 19 Simulated region for velocity profile in a circular tube.
v/vmaxr/r0
Fig. 20 Velocity profile of flow in a circular tube calculatedwith the intermolecular collision scheme of Bird's original.
v/vmaxr/r0
Fig. 21 Velocity profile of flow in a circular tube calculatedwith the intermolecular collision scheme of U-system.

Fig. 17 Velocity profile of flow in a circular tube.
8

Copyright © 2013 by JSFM

第 27 回数値流体力学シンポジウム
E01-2
運動量とエネルギの厳密保存が可能な，従来より安定性が増した U-system（U-system-8）を用いて，Re 変化に伴う円管内流速分
布の変化を DSMC 法により再現した．今回は，まだテスト計算の域を超えてはいないが，Bird 法では得られなかった層流分布からのずれが U-system を用いた計算により得られたことは，今後，分
子シミュレーションによる乱流速度分布再現の可能性に道を開くものと考える．

Fig. 20 に，Birdオリジナルの分子間衝突計算法によって計算した流速分布を示す．横軸は，管中心からの位置rを管半径r0で無次元化した値を示し，縦軸は，流速vを，初期条件で与えた流速から得られるであろう層流放物線分布の最大値vmで無次元化した値で示している．Re=2 およびRe=10 では，Knが 1 より大きくなる
ので，壁面上での滑り現象に代表される「希薄化の影響」によって，流速分布は，層流分布とは異なるものが得られている．しかし，Re=50 以上では，いずれの結果も，放物線分布となって，層流速度分布が出現している（Re=5000 の最大値が 1.0 を少し超えている理由は，最初に与えた分子流速の偏りが原因である可能性が高い）．一方，Fig. 21 は，U-systemによる計算結果で，まずRe=2とRe=10 では，希薄化効果によりBird法と同じ結果となっている．
そして，Re=50 とRe=200 は層流分布となっているが，Re=1000 では放物線分布から少しずれ，Re=5000 では，そのずれはさらに大きくなって，乱流分布に近づく傾向にある．ただし，今回は，層流分布からのずれは確認できたものの，1/7 乗法則のような乱流分布にはほど遠い結果しか得られていない．
さて，Bird 法では，
Re=50 以上のいずれでも層流分布が得られ，
一方，U-system では，Re=1000 と 5000 で放物線分布からのずれが見られたのは，どういう理由からであろうか．U-system がセルを極端に小さくしたものに近い計算だからであろうか．そうであるなら，Bird 法であっても，セルを極端に小さくすれば，層流分
布からのずれが現れるであろうか．セルを細かくすることにより計算時間が多くかかることは言うまでもないが，もし流れ場中の小さな不規則な乱れが乱流を発達させる原因とするなら，セルを小さくすればするほど，いったん生じた微細な乱れも，殆どの場
合，
途中でつぶされてしまい，
それが大きな渦に発達するのには，途方もなく時間がかかる可能性がある．
数値計算そのものには乱れ要素が含まれていないのに，流れの中に乱れが発生して層流速度分布からのずれが生じたとしたら，
DSMC 計算で用いている乱数のゆらぎ，あるいは，壁面境界でのセル配置の不規則性等が原因である可能性が考えられる．いずれにしても，今回の流速分布は乱流分布にはほど遠いものであり，今後は，計算に，壁面粗さの影響を加える等の外的な乱れを加えることにより，流速分布の乱流化が可能になるかどうかを調べる必要がある．二次元解析から三次元解析に拡張することも必要であろう．
現在，Re=10000 は計算中でまだ結果が出ていないが，途中経過を見る限りでは Re=5000 の結果と差異はあまりないようである．
なお，完全平衡法も使用してみたが，この方法では Re を特定できないし，また，期待するほどの結果も得られなかった．今回，他に，断面格子数 1600×1600 についても，また，平均流速 69 m/sと 17.3 m/s についても計算を行っているが，まだ報告できるだけの十分な結果としてはまとめられていない．

参考文献(1) Bird,G.A., Molecular Gas Dynamics and the Direct Simulation of Gas
Flows, (1994), Clarendon.
(2) 日本機械学会編，計算力学ハンドブック第Ⅲ巻（4.1 DSMC 法の基礎），(2009)，pp. 130-149，日本機械学会．
(3) Meiburg, E., Physics of Fluids, 29-10 (1986), pp. 3107-3113.
(4) Usami, M. and Nakayama, T., Rarefied Gas Dynamics, AIP
Conference Proceedings, 663 (2003), pp. 374-381.
(5) Usami, M. and Mizuguchi, K., Rarefied Gas Dynamics, AIP
Conference Proceedings, 762 (2005), pp. 686-691.
(6) Usami, M., Kamiya, T., Maeda, H. and Furukawa, Y., Rarefied Gas
Dynamics, AIP Conference Proceedings, 1084 (2009), pp. 365-370.
(7) Sun, Q., Liu H. and Fan J., Rarefied Gas Dynamics, AIP Conference
Proceedings, 1084 (2009), pp. 311-316.
(8) Pareschi, L. and Trazzi, S., Int. J. Numerical Methods in Fluids, 48-9
(2005), pp. 947-983.
(9) 日本鋼構造協会編，構造物の耐風工学，（1997），東京電気大学
出版会．
(10) Bird, G. A., Physics of Fluids, 30-2 (1987), pp. 346-366.
(11) 宇佐美，ながれ，26-4 (2007)，pp. 273-282．
(12) 宇佐美・井村，機論，B 編，75-749 (2009)，pp. 1-10.
(13) Koura, K., Progress in Astronautics and Aeronautics (AIAA), 117(1989), pp. 25-39.
(14) 宇佐美ほか，
日本流体力学会年会2010講演要旨集およびCD，(2010)，p. 112.
(15) 宇佐美ほか，物体後流の高速渦はく離流れに関するDSMC分子シミュレーション，日本流体力学会年会2012講演論文集USB，(2012)．

８．まとめ
テーパ円柱の後流に生じる渦はく離流れを，分子シミュレーシ
ョン法である DSMC 法を用いて解析した．
DSMC 法においては，
流れ場空間をセルと呼ばれる微小体積要素に分割する際，通常，セルの大きさを平均自由行程以下にする必要がある．しかし本解析では，分子間衝突の前後に，衝突に関与する分子の速度をその位置座標に応じて修正するという手法 U-system により，セルを大きくすることを可能にし，これにより，三次元場における渦はく離流れの解析を実現した．なお，U-system において，分子の運動量と並進エネルギを厳密保存するためにPareschi らの補正法を使うことは有効である．テーパ円柱後流の流れは，著しい三次元性を持ち，後流渦の内部を通る流線は，主流を横切って円柱直径の細くなる方向に進む傾向が見られる．
9

Copyright © 2013 by JSFM