Tips

どんなアプリ？

• 2次元の流体計算をブラウザにさせるウェブアプリです。この画面が表示された後は、サーバーとの間で通信をしません。
• 計算を開始すると、ブラウザがお使いのデバイスの CPU コアを 1 個占領すると思います。メモリは数百 MB 程度使うかもしれません。
• 計算結果の更新スピードが「計算スピード」欄に表示されます。新しめのパソコンならば、5 - 6 fps 程度の値が出ると予想されます。最大値は 8.0 fps です。
• お絵かきモードで自由に描いた物体を、一様流中に置くことができます。
• いろいろなボタンやつまみをイジって楽しんでください。
• 過激な設定の計算をすると、途中で結果が発散してしまうかもしれません。その過程も含めてお楽しみください。発散しても、初期化ボタン ( ) を押せば計算をやり直せます。

計算結果について

• 流速、渦度、圧力の計算結果はカラーマップとして図示されます。無限遠での値をグレー (ゼロ) として、それより大きければ赤色 (正の値)、小さければ青色 (負の値) になります。ただし、スケールや原点は変更できます。
• カラーマップのスケールとは、図示された色が頭打ちになるときの値と原点の差です。小さく設定すれば色が濃くなります。流速と圧力について、値 1.0 は無限遠での絶対値を指します。スケールは一時停止中のみ変更できます。
• カラーマップの原点は、グレーで図示されるときの値です。原点も一時停止中のみ変更できます。
• 渦度とは、流れの回転の度合いを表す物理量です。その付近の流れに左回りの傾向があれば赤、右回りの傾向があれば青になります。
• ベルヌーイ式とは、ベルヌーイの定理において保存される式 (動圧 + 静圧 / 全水頭 / 全エンタルピーの変化量) のことを指します。
• レイノルズ数とは、流体の粘性が現象に与える影響力の指標です。レイノルズ数が小さいほどドロドロの流体 (あるいは小さな領域)、大きいほどサラサラの流体 (あるいは大きな領域) を計算していることになります。
• 抗力、揚力とは、流体が物体に与える合力のことです。一様流方向 (今回は右向き) の成分を抗力、それに垂直な方向 (今回は上向き) の成分を揚力と言います。抗力係数、揚力係数はそれらの大きさの指標です。

結果から分かることの例

• 流線との関係性が可視化されるので、渦度の意味を理解することができます。
• 渦の中心は低圧になることから、渦における遠心力と圧力傾度力のつりあいがイメージできます。
• レイノルズ数による流れ方の違いを大雑把に理解できます。
• ベルヌーイの定理の成立条件を理解できます。ベルヌーイの定理は、定常流かつ粘性が無視できる場合に、ある瞬間の同じ流線上でベルヌーイ式の値が一定値を取るというものです。例えば、レイノルズ数を高く設定した場合、物体の上流側は定常流に近くなりますから、ベルヌーイ式が一定値を取ります。逆に、下流側は定常とはみなせませんから、ベルヌーイ式は変動します。また、レイノルズ数を低く設定した場合は、非粘性の仮定が崩れますからベルヌーイの定理は成り立ちません。
• その他、上手く物体壁を描き、レイノルズ数を低く設定すれば、圧力降下などの現象の理解に役立つかもしれません。
• 面白い現象を見つけてみてください。ただし、その現象が現実にも起きうるのか、数値計算の設定による人工物なのかを判別するのは難しい場合が多いです。

結果を見るときの注意点

• ブラウザでリアルタイムに計算できるように、格子は粗く設定してあります。
• それに関連して、物体表面付近での流れの計算には大きな誤差が伴います。このため、抗力や揚力の大きさについて、特にレイノルズ数が大きい場合には、上手く計算できません。
• レイノルズ数もいくらか過大評価していると思われます。
• 図示しているのは計算領域のうちの一部であり、図の外における流れも裏で計算しています。しかし、図の境界付近に物体を描くと、その先にある数値境界の人工的な影響が表れる可能性があります。
• 流出側 (計算領域の右辺) の数値境界条件には壊れやすいものを用いているため、高レイノルズ数での計算の破綻はこの境界が原因である可能性が高いです。
• 物体表面では流速がゼロになるような境界条件が課されています。
• 計算開始直後の状態硬直を避けるために、初期条件には一様流からの微小な擾乱が加えてあります。
• 計算直後には、物体を突然置いたことによる激しい音波や定在波が発生します。落ち着くまでにはしばらく計算を進める必要があります。定在波の影響は、特に圧力に強く表れ、しぶとく残ります。
• 用いている計算スキームには、乱流モデルは組み込まれていません。また、計算格子も粗いので、乱流に関する現象は再現できません。
• 実際の車や飛行機、船などの場合、レイノルズ数は 10 の 6 乗から 9乗のような、ずっと大きい値になります。このアプリの結果をそのままこれらの系に当てはめられるわけではありません。実用的なシミュレーションには相応の計算機または計算時間が必要です。

計算の詳細

物理設定

• 流体は非圧縮性を想定しています。つまり、空気というよりは水です。しかし、今回の問題設定・計算精度では、両者に共通する大雑把な性質を再現するのが精一杯なので、気にする必要は無いかもしれません。
• 図に書かれている目盛りを単位にして長さを表すことにします。
• 一様流の速度はマッハ 0.1 です。これは、水にしてはかなり大きいです。陽解法スキームを用いているため、計算速度をなるべく上げるための処方です。一般に、レイノルズ数を一定に保ちながら流速の設定を変えると、圧力変化の大きさは変わりますが、流れの様子は定性的には変わらないことが多いです。
• 一様流の速度が 1 になるように物理時間の単位を取ることにします。このとき、一度の図の更新で物理時間が 0.05 だけ進みます。
• レイノルズ数 Re は次のように定義されます。図の 1 目盛りが表す長さを L、一様流の速度を U、流体の動粘性率を ν としたとき、Re = LU/ν です。例えば、常温の水を想定し、L = 10 cm、U = 1 cm/s としたときに、Re = 1000 となります。

計算領域

• 図には x 方向に 8、y 方向に 4 の領域が示されていますが、計算領域はその外側にも広がります。図の左下を原点 (0, 0) としたとき、およそ −4 ≤ x ≤ 16, −8 ≤ y ≤ 12 の正方形領域が計算されます。
• 図に示された領域には 179 × 89 個の計算格子が存在します。図の外側の計算領域には、図の領域に比べて 2 倍の粗さ (面密度にして 1/4) の計算格子が設定されています。計算格子の総数は 61344 個です。
• 図を一度更新するために、図の外側の領域は 9 回、図の領域は 18 回、後述する時間積分スキームを実行します。

計算スキーム

• 格子ボルツマン法 (LBM, e.g. [1]) を用いています。この方法では、セルオートマトン風に離散化されたボルツマン方程式が計算されます。
• LBM の利点は、比較的計算が速いことと、境界条件の設定が容易なことです。お絵かきモードの実装のために選びました。
• LBM は擬似圧縮性解法の一種です。解く方程式は、近似的に非圧縮性の質量保存則&ナビエ-ストークス方程式に帰着するものですが、解く際に若干の圧縮性を認めます。
• 空間 2 次精度の解法を用いています。科学研究に使われるような有限差分法や有限体積法と比べると、精度は落ちます。また、物体形状に対応して計算格子の形を変化させることができないのも欠点です。因みに、3 次元計算に関しては、空間 4 次精度の LBM スキームも存在し、乱流遷移の再現を行った研究もあります [2]。
• 衝突項の計算には、セントラルモーメントを用いています。具体的には、文献 [3] の方法です。これにより、高レイノルズ数でも安定に計算を進められます。しかし、今回は格子点数が粗いので、精度には期待できません。
• 物体表面での固定境界条件には、ハーフウェイバウンスバックスキーム (e.g. [1]) を用いています。
• 計算領域の上辺、下辺、左辺の数値境界では、局所平衡分布を仮定して未知量を補完しています。
• 計算領域の右辺の数値流出境界では、文献 [4] の "convective condition" を用いています。領域内に与える影響が少ないのが利点ですが、安定性は他の境界より悪いです。
• 細かい計算格子の領域 (図示された領域) と粗い計算格子の領域は、空間的に 3 次スプライン補間、時間的に 2 次ラグランジュ補間を用いて滑らかに接続しています。

可視化手法

• 流線の可視化には LIC (Line Integral Convolution) 法 (e.g. [5]) を用いています。用いる入力ノイズとカーネルには文献 [6] のものを採用しました。

アプリの構成

• UI とお絵かきモードは JavaScript、流体計算と可視化は WebAssembly によって実装されています。
• 並列化なしの CPU による計算です。

更新情報

更新履歴

• 2023 年 2 月 16 日：公開
• 2023 年 2 月 16 日：大きめのメモリリークが起きて、30 分くらい計算を続けるとクラッシュしてしまう問題を解決。
• 2023 年 2 月 24 日：ベルヌーイ式も図示できるようにした。
• 2023 年 2 月 27 日：線を描きすぎると図が真っ暗になってしまう問題を解決。図の端に線を描くとフリーズする問題を解決。

現在分かっている問題

やり残したこと

• GPU を使うと計算速度がどう変わるのか試してみたい。

参考文献

• [1] 瀬田剛, (2021),『格子ボルツマン法』(森北出版).
• [2] M. Geier, A. Pasquali and M. Schönherr, (2017), Parametrization of the cumulant lattice Boltzmann method for fourth order accurate diffusion part II: Application to flow around a sphere at drag crisis, Journal of Computational Physics, 348, 889-898, doi: 10.1016/j.jcp.2017.07.004.
• [3] A. De Rosis, (2016), Non-orthogonal central moments relaxing to a discrete equilibrium: A D2Q9 lattice Boltzmann model, Europhysics Letters, 116, 44003, doi: 10.1209/0295-5075/116/44003.
• [4] Z. Yang, (2013), Lattice Boltzmann outflow treatments: Convective conditions and others, Computers and Mathematics with Applications, 65, 160-171, doi: 10.1016/j.camwa.2012.11.012.
• [5] B. Cabral and L. C. Leedom, (1993), Imaging vector fields using line integral convolution, SIGGRAPH '93: Proceedings of the 20th annual conference on Computer graphics and interactive techniques, 263-270, doi: 10.1145/166117.166151.
• [6] R. Wegenkittl, E. Groller and W. Purgathofer, (1997), Animating flow fields: rendering of oriented line integral convolution, Proceedings. Computer Animation '97, 15-21, doi: 10.1109/CA.1997.601035.

注意事項

• アプリの著作権は作者に帰属します。他所で再配布を行ったり、他のコンテンツの一部として用いることはできません。
• 何かございましたら以下の連絡先にお願いします。
• 峰田竜二: contact@solarphys.com

太陽の科学館トップページへ