(Click here for explanation in English.)
（以下に示す記述は，暫定版です．現在，プログラムを少し改良して計算を繰返していますが，その結果自体は，暫定版のものと殆ど変わりません．最新版の公開は，今年の秋ごろを予定しています．）
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/30 だけの揺らぎ（摂動）が与えられる．
(b) 壁面粗さを鈍い鋸刃形状（鋸刃の間隔と刃の角度）で近似する方法について
鋸刃の間隔に関しては，平均流速に比例する値を用いた（Re数にも比例する）．すなわち，管路長（代表長さとしての管路直径に近い値）と鋸刃間隔の比を，分子の最大確率速さと流れの平均流速の比に，ある係数を掛けたものと等しくなるようにする．これは，物理的に無理のない仮定と言える．具体的なプログラムでは，最大確率速さと平均流速の比を求めた後，係数1.39を掛け，それを整数化（切捨て）した値で管路長を割算して鋸刃間隔を求める．なお，この係数を変えれば，境界層厚さを調整できることが判明している．Re数に応じて実際の壁面形状が変化することはないはずであるが，これについては次のように考える．今回提案している「宇佐美の擬似表面粗さ」は，現実の表面粗さが持つ極めて複雑な性質を，単純な鋸刃の形状を用いて模擬するものである．その場合，実際の壁面形状においては鋸刃の間隔が多岐にわたる一方，あるレイノルズ数（すなわち，ある平均流速）に対して特定の鋸刃間隔だけがそのレイノルズ数に共振して流体に影響を及ぼすと考えれば，広範なレイノルズ数域でこの簡略化された近似を適用するために，鋸刃の幅を平均流速に比例させて伸び縮みさせるのは妥当なことである．刃の角度はRe数に依らず一定で，今回は 35度とする．今回の改良した計算法では刃角一定が可能になった．ただしこれも，実際の壁面粗さでは多くの異なる角度が存在していて，その流速に共鳴する角度だけが表面化（顕在化）すると考えれば，角度を固定する必要はないのかも知れない．管路の上流と下流では周期境界を用いるので，鋸刃の刃数は整数に限る必要がある．そのため，Re数に応じて多少管路長が変わり，本来の設定である菅直径（代表長さ）に等しくすることが出来なくなる．今回は，小数点以下を切り捨てて整数化したので，流れ場の長さ（管路長）はRe数に応じて代表長さより少し短くなってくる．
(c) Bird法を用いた乱流計算
以前の3.1節の説明では，Bird法では，層流の放物線分布を厳密に計算できないという理由で，それ以上は説明を省略していた．今回は，従来のBird法に，Pareschiらの運動量とエネルギの厳密保存の方法を追加したものを「修正Bird法」と呼ぶことにして，円管内流速の解析に適用して調査した．
以下の計算では，Re数 1000～7000 の範囲で，それなりの結果の得られるように，鋸刃間隔のための係数を 1.39，鋸刃角度 35度，鋸刃形状のゆらぎ幅を，管路長（計算領域長さ）の 1/30 と設定して計算を実行した．なお，これらを「３つのパラメータ」と呼ぶ．係数1.39 を変更すると境界層厚さを変化させることができるが，変更した値を用いて広いRe数の範囲で良好な結果を得るためには，それ以外のパラメータを注意深く選定する必要があることを忘れないでほしい．
3.3.1 Re=5500 での計算
まずは，Re=5500 で計算を行った．Figure 1 と 2は，管内流速分布の時間経過（動画）であり，初期流速分布は層流の放物線分布からランダムに分子速度を定めて計算を開始したもので，そこから乱流速度分布が完成するまでの変化が計算されている．なお，今回の計算におけるセル数は約103万個，分子数は約1530万個（Re数により計算領域の大きさが変化するので分子数もReにより多少変化する）で，乱数を用いて最初の分子速度を決定するので，それによって平均流速すなわちRe数を算出すると必ずしも希望するRe数にぴったり一致する値にはならないことに注意されたい．1回目の結果出力は，45ステップ後（分子移動と分子間衝突を分離する時間ステップΔtm＝1.1μsの200倍のさらに45倍，すなわち 9.9 ms となるが，その間の平均値なので，約 5 ms 経過時点での結果となる），2回目以降は，そこから15ステップ毎に得られたものを３つ平均した出力結果で，いずれも平均流速で無次元化している．
Figure 3 と 4 は，それぞれ，計算開始から乱流分布が完成するまでに得られた密度分布（数密度分布）と温度分布である．いずれも初期密度と初期温度で無次元化している．流れ場全体では，どちらも1%～2%の変化にとどまっているが，壁面近傍における密度には，独自の壁面反射処理の影響で5%程度の変化が生じている．
Figure 5 は，表面粗さのない拡散反射面（乱反射面）を用いて計算した結果で，結果は，完全な放物線分布（層流速度分布）となっている．なお，この場合の計算の初期分布には，平均流速で一様な分布を用いた．
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=2500～7000 は，45ステップでの平均であるが，平均流速の小さい状態では，特に菅の中心付近で変動が大きくなるので，Re=1500と2000では，90ステップでの平均を，Re=1000 では，180ステップでの平均によって平滑化している．低速で変動が大きくなる理由は，DSMC法が分子速度（最大確率速さが約 350m/s）を基盤に解析されるためである．Re＝1000～2000 は層流を示す放物線分布であり，Re=2500～3500 は，層流から乱流への遷移を示しており，Re=4000～7000 で，ほぼ乱流分布に達したと言うことができる．ただし，Re=6500 と 7000 の速度分布は，他の乱流分布から離れ始めている．この理由が，今回の計算方法における欠陥なのか，あるいは流速が大きくなり過ぎて圧縮性の影響が出てきたためなのか等のことは，今後検討すべき課題である．なお先に述べたように，修正Bird法（Bird法も同じ）では，Re数が小さくなると，結果が Usys法からずれてくる．一例として，Figure 22 に，Re=2500 における比較を示しておく．
先に述べたように，使用した３つのパラメータについては，今回の値を使用せずに，これを変えてもっと望ましい乱流・層流分布を得ることも，言い換えれば，別の境界層厚さについて挑戦することも可能であろうし，さらに，新たなパラメータを追加してより現実により近い管内流速分布を導くこともできるかも知れない．しかしそのためには，大変面倒で大がかりな繰返し計算の調査が必要となるであろう．
管内流速における乱流と層流の存在とその遷移について，著者は，元来，単純な数学で導くことのできる層流が基本にあって，何らかの原因により，そこから乱流が生み出されると考えていた．しかし，今回のDSMC計算を行って言えることであるが，「流れは本来乱流であり，それが何らかの理由で乱流成分が打ち消されて層流が作られる」と考えた方が自然のように思われる．さらに，層流の方が，DSMC計算で用いる「セル」は小さいことが要求される．これらのことをどう解釈するかは，今後の読者諸君の考察に任せたい．