(Click here for explanation in English.)
（この節の記述は，2026.9 に作成した最新版です．なお，参考のために6月に書いた古い暫定版は 3.3′[Old …..] に残してあります．）
2年前の2024年に，本WEBサイトの第3.1節と第3.2節において，円管内乱流速度分布のDSMC計算を行って，その可能性を論じたが，そこには解決すべき幾つかの問題が存在していた．すなわち，(a) 乱流速度分布を得るために，壁面近傍の複数の離れた位置に流速の局所ピークを発生させ，それによるレイリー・テイラー不安定を元にして管内全体を乱流状態に導いたが，実際の流れでは，壁面近傍の複数の離れた場所にそのような流速ピークが点在するわけではない．(b) 壁面粗さを，鈍い鋸刃形状で近似（DSMC法における宇佐美の疑似表面粗さ）したが，鋸刃の間隔と角度をレイノルズ数Reに応じて変化させるにあたり，物理的に意味を持たない近似式を作って対応させた．しかし，物理的な意味を説明できる数値，あるいは，レイノルズ数Reに依らない一定値を用いることが望まれる．(c) 専らUsys法を用いた計算に限定して論じていたが，Bird法に関する説明を加える必要がある，ということが挙げられる．
今回は，まずこの第3.3節において，改良した計算法について説明を行って，その計算結果を示し，次いで第3.4節では，そのソフトウェアを公開して，それの使用方法を説明する．なお，Re数の範囲は，前回と同様に 1000～7000とし，やはり，平均流速を変更することによりそのRe数を作り出している．また，対象気体や計算領域の大きさ，さらに，管路の上流と下流で周期境界を用いる等も，前回と同様な考えを使っている．すでに 第3.1節や 第3.2節に述べたことで変更のない事項に関しては，ここでは記述を省略しているので，それらのことをあらかじめ熟知した上で，今回の解説を読んでもらいたい．
(a) 壁面近傍の離れた複数の位置に流速ピークを作らない方法
これは比較的簡単に実現できた．すなわち，鋸刃形状を一点に固定するのではなく，流れ方向にゆっくりと移動させる．ただし，速く動かしすぎると乱流構造が崩れてしまうので，今回は，一つの鋸刃間隔を45ステップで動かすことにした．ここで，１ステップとは，本DSMC計算の分子移動と分子間衝突を分離している最小時間単位 Δtm（DTM=1.1μs）の200倍（220μs）を意味する．なお平均流速に比例（すなわちRe数に比例）して鋸刃間隔は変わるので，１ステップの鋸刃形状の移動量はRe数によって異なる．またこの場合，必ず45ステップで途中の出力結果の平均を取らないとデータが歪むことになる．途中結果は，15ステップずつ平均化され，さらにそれを３個まとめて１回の結果とする．すなわち，最初は，1～45ステップの平均となるが，2回目は，16～60ステップの平均結果というように計算が繰り返される．なお，計算時に作成される途中データは子フォルダ DATx に保存され，それは大きなファイルであるが，必要なくなれば自動的に削除される．さらに，分子が壁面で反射される際，移動中の鋸刃形状は，その時点の位置を中心にして代表長さ（管路直径）の1/35 だけの揺らぎ（摂動）が与えられる．
(b) 壁面粗さを鈍い鋸刃形状（鋸刃の間隔と刃の角度）で近似する方法について
鋸刃の間隔に関しては，平均流速に比例する値を用いた（Re数にも比例する）．すなわち，代表長さ（管路直径）と鋸刃間隔の比を，分子の最大確率速さと流れの平均流速の比にある係数を掛けたものと等しくなるようにする．これは，物理的に無理のない仮定と考える．具体的なプログラムでは，最大確率速さと平均流速の比を求めた後，係数1.30を掛け，それで代表長さを割算して鋸刃間隔を求める．管路長は，代表長さを鋸刃間隔で割算して得た値の整数部分に鋸刃間隔を掛けたものである．なお，この係数1.30や，揺らぎ（摂動）の幅1/35，また，次に示す鋸刃角度35度を変えることにより，Re数に応じた流速分布を今回の結果とは異なるものに調整できる．ところで，Re数に応じて実際の壁面形状が変化することはないはずであるが，これについては次のように考える．今回提案する，「乱流を引き起こすための外乱」としての「宇佐美の擬似表面粗さ」は，現実の表面粗さが持つ極めて複雑な性質を，単純な鋸刃の形状を用いて模擬するものである．その場合，現実の壁面形状においては鋸刃の間隔が多岐にわたる一方，あるレイノルズ数（すなわち，ある平均流速）に対して特定の鋸刃間隔だけがそのRe数に共振して流体に影響を及ぼすと考えれば，広範なレイノルズ数域でこの簡略化された近似を適用するために，鋸刃の幅を平均流速に比例させて伸び縮みさせるのは不自然な仮定ではない．刃の角度はRe数に依らず一定で，今回は 35度とする．今回の改良した計算法では刃角一定が可能になった．ただしこれも，実際の壁面粗さでは多くの異なる角度が存在していて，その流速に共鳴する角度だけが表面化（顕在化）すると考えれば，角度を固定する必要はないのかも知れない．これは今後の課題である．管路の上流と下流では周期境界を用いるので，鋸刃の刃数は整数に限る必要がある．そのため，Re数に応じて多少管路長が変わり，当初想定したように管路直径（代表長さ）に合わせることは必ずしも出来ない．今回は，小数点以下を切り捨てて整数化したので，流れ場の長さ（管路長）はRe数に応じて代表長さより少し短くなっている．
(c) Bird法を用いた乱流計算
以前の3.1節の説明では，Bird法を用いると，層流の放物線分布を厳密に計算できないという理由で，それ以上は説明を省略していた．今回は，従来のBird法に，Pareschiらの運動量とエネルギの厳密保存の方法を追加したものを「修正Bird法」と呼ぶことにして，円管内流速の解析に適用して調査した．
以下の計算では，Re数 1000～7000 の範囲で，それなりの結果の得られるように，鋸刃間隔のための係数を 1.30，鋸刃角度 35度，鋸刃形状のゆらぎ幅を，管路直径（代表長さ）の 1/35と設定して計算を実行した．なお，これらを「3つのパラメータ」と呼ぶ．3つのパラメータを変更すると，Re数変化に応じた流速分布を今回の計算結果とは異なるものにすることができるが，広いRe数範囲に亘って良好な結果を得るためには，3つのパラメータを注意深く選定する必要があることを忘れないでほしい．
3.3.1 Re=5500 での計算
まずは，Re=5500 で計算を行った．Figure 1 と 2は，管内流速分布の時間経過（動画）であり，初期流速分布は層流の放物線分布からランダムに分子速度を定めて計算を開始したもので，そこから乱流速度分布が完成するまでの変化が計算されている．なお，今回の計算におけるセル数は約103万個，分子数は約1300万個（Re数により計算領域の大きさが変化するので分子数もReにより多少変化する）で，乱数を用いて最初の分子速度を決定するので，それによって平均流速すなわちRe数を算出すると必ずしも希望するRe数にぴったり一致する値にはならないことに注意されたい．1回目の結果出力は，45ステップ後（分子移動と分子間衝突を分離する最小時間単位 Δtm＝1.1μsの200倍のさらに45倍，すなわち 9.9 ms となるが，その間の平均値なので，約 5 ms 経過時点での結果となる），2回目以降は，そこから15ステップ毎に得られたものを３つ平均した出力結果で，いずれも平均流速で無次元化している．
Figure 3 と 4 は，それぞれ，計算開始から乱流分布が完成するまでに得られた密度分布（数密度分布）と温度分布である．いずれも初期密度と初期温度で無次元化している．最終的な密度および温度の変化の程度は，壁面近傍を除く流れ場全体では，どちらも1%以下にとどまっているが，壁面近傍では，独自の壁面反射処理の影響で，密度で6%程度，温度で2%程度の変化が生じている．
Figure 5(a) は，表面粗さのない拡散反射面（乱反射面）を用いて計算した流速分布で，結果は，ほぼ完全な放物線分布（層流速度分布）となっているが，Figure 5(b) のオレンジ色部分のように，中心軸近傍で平均流速の2倍（理論上の最大値）を超える値が観測される場合がある．中心軸付近の密度は平均より少し低いので，本来は，流速ではなく質量速度（質量フラックス）で表現すべきである．そこで，これで調べたところ多少は改善されたが，その傾向が大きく変わることはなかった （Figure 5(c)）．なお，この場合の計算の初期分布には，平均流速で一様な分布を用いた．
Figure 6 は，修正Bird法により，乱流速度分布を求めたものである．修正Bird法の計算時間は，今回のDSMC計算では Usys法の約７割しか必要ないので，高速で計算が可能であり，結果も，Usys法と殆ど変わっておらず，その意味では非常に優れている（正確を期すために修正Bird法を用いたが，単なるBird法でも結果に変わりは見られなかった）．しかし，修正Bird法を使って，拡散反射壁面（壁面粗さなし）により層流速度分布を求めてみると，Figure 7 のような結果となり完全な放物線分布は得られない．また，壁面粗さを有する場合でも，Re数が小さくなるほどUsys法との差異が僅かに見られるようになる．Usys法が修正Bird法と異なっているのは，「空間を分割しているセルが見かけ上小さくなる」ということなので，乱流速度分布の計算にはセルをそれほど細かくする必要がないということが言えるかも知れない．しかし，層流速度分布においては十分な結果が得られないことを考えると，修正Bird法には何かしら不足したものがあると思われる（従来のBird法でも同様）．すなわち，その理由を明確にしない限り，修正Bird法を信頼して使用できないのではないかという疑念を払拭できない．
3.3.2 Re=1000～7000における流速分布
Figure 8～20 は，Re=1000～7000（平均流速は，おおよそ 10～70 m/s）で得られた管内流速の定常状態到達後の分布を，500 おきに描いたものである．また，Figure 21 は，それらをコマ送りで描いたものである．図はすべて平均流速で無次元化している．なお，Re=3000以上は45ステップでの平均であるが，平均流速の小さい状態では，特に菅の中心付近で変動が大きくなるので，Re=1000と1500および2000では180ステップでの平均を，Re=2500 では135ステップでの平均によって平滑化している（Re=3000以上であっても滑らかな結果を希望するなら90ステップや135ステップの平均を使うとよい）．低速で変動が大きくなる理由は，DSMC法が分子速度（最大確率速さが約350m/s）を基盤に解析されるためである．Re＝1000～2000 では層流を示す放物線分布が現れ，Re=2500～3500 では，層流から乱流への遷移を示し，Re=4000～7000 では，ほぼ乱流分布に達している．ただし，Re=6500 と 7000 の速度分布は，他の乱流分布から離れ始めている．この理由が，私の提案した「表面粗さ近似」を中核とする今回の計算方法の欠陥なのか，あるいは流速が大きくなり過ぎて圧縮性の影響が出てきたためなのか等のことは，今後検討すべき課題である．なお先に述べたように，修正Bird法（Bird法も同じ）では，Re数が小さくなると，結果が Usys法からずれてくる．一例として，Figure 22 に，Re=2500 における比較を示しておく．
先に述べたように，使用した３つのパラメータについては，今回の値を使用せずに，これを変えてもっと望ましい乱流・層流分布を得ることも可能であろうし，さらに，新たなパラメータを追加してより現実により近い管内流速分布を導くこともできるかも知れない．しかしそのためには，大がかりで面倒な繰返し計算による調査が必要となるであろう．
管内流速における乱流と層流の存在とその遷移について，著者は，元来，単純な数学で導くことのできる層流が基本にあって，何らかの原因により，そこから乱流が生み出されると考えていた．しかし，今回のDSMC計算を行って感じたことであるが，「レイリー・テイラー不安定を作り出す能力はReの大小に依らず存在するはずなので，その意味で流れは本来乱流と考えられるが，Re数が小さくなると慣性力が弱まって渦を持続させることができなくなり，乱流分布が崩れて層流になる」 と考えるのも悪くないように思われる．また，層流の方が，DSMC計算で用いる「セル」は小さいことが要求される．これらのことをどう解釈するかは，今後の読者諸君の考察に任せたい．さらに，私は，手元に実験結果を持っておらずシミュレーションの良否を明確に判定できないので，このことは是非，読者諸君に調査を依頼したい．
流速分布を求めるこの解析には，比較的長い計算時間が必要である．具体的には，一つのRe数の計算完了までに要する時間は，CPUにRyzen9_9950X3D（9950Xでも9950X3D2でも計算時間はほぼ同じ）を使用すると約44時間（32 thread使用），Ryzen7_6800H を用いると倍の約88時間（16 thread使用）となる（ただし，共に，M.2のSSDを使用した場合）．なお，DSMC計算に9950X3Dを使用する場合，高性能グラフィックボードは不要であり，CPUの冷却は空冷が可能である（著者の経験では，安価な簡易水冷はかえって性能が悪いし，温度制御の応答が遅く故障しやすい）．Ryzen7_6800Hは少し以前に発売されたCPUであるが，価格の割に性能が良く，ノートPC用なので消費電力も低く，コストパフォーマンスに優れている（2025年末に，ほぼ同性能のCPUが Ryzen7_170 として発売された．性能は若干落ちるが消費電力は低いらしい）．しかし，長時間計算というハードな使い方をするので頑丈なノートPCを手に入れる必要がある．著者の場合，Lenovo製16インチLaptopを使ってディスプレイをOFFにしたら，計算実行中は80W程度の消費電力であった（空冷用の空気取り入れ口は，電気掃除機で月に1回は清掃することを心掛けて）．