この文書では、WSL2のUbuntu24.04上でLBM（格子ボルツマン法）流体解析ツールXLBを活用する方法を紹介します。
「XLB（Accelerated Lattice Boltzmann）」は、Autodeskが公開している、GPUによる高速計算と機械学習との連携を重視したPython製LBM（格子ボルツマン法）ライブラリです。2次元・3次元の流体解析に対応し、標準的なBGKモデルや高レイノルズ数流れに適したKBCモデルなどを利用できます。
最大の特徴は、NVIDIA Warp、Google JAX、Autodesk Neonなどの計算バックエンドを利用できる点です。いずれも、PythonコードからGPUを活用した高速なLBM計算を比較的容易に実装できます。
以下の【07:Warp導入方法と動作確認】では、バックエンドにWarpを GPUで利用する環境を構築し、例題rotating_sphere_3d.pyの動作確認まで進めました。
この【08】では、【07】の手順で構築したWarpのGPU活用の環境をもとに、例題 rotating_sphere_3d.pyを分析してゆきます。
https://github.com/Autodesk/XLB/blob/main/examples/cfd/rotating_sphere_3d.py
XLBのGitHub公開例題の中で、最初に理解する例題としては lid_driven_cavity_2d.py（2次元 Lid-Driven Cavity）があり【02】で分析しました。次に3次元の例題として、【04】では flow_past_sphere_3d.py を取り上げ、XLBによる3次元流体解析の基本とスクリプトの設定について分析しました。さらに【06】では multires_windtunnel_3d.py を取り上げ、マルチ解像度格子によるXLBによる3次元流体解析のスクリプトの設定について分析しました。
本記事【08】では、バックエンドとしてWarpを用いた例題として rotating_sphere_3d.py を分析します。この例題では、回転する球周りの3次元流れを対象として、XLBのWarpバックエンドを用いたLBM解析を行います。
さらに、D3Q27速度モデル、STL形状からの計算格子生成、球表面の回転速度を考慮した境界条件、レイノルズ数200による流れ条件の設定、抗力・揚力などの流体力の計算、計算結果の可視化など、XLBを3次元の流体解析へ応用するための基本的な機能が含まれています。
そのため、本例題は、flow_past_sphere_3d.py などで学んだXLBの基本からさらに進み、任意の3次元形状をSTLデータとして扱い、回転する固体表面と流体との相互作用をLBMで解析する方法を学ぶための例題として位置づけることができます。特に、XLBのWarpバックエンドを利用してGPU上で3次元LBM計算を実行する流れを理解する上でも重要な例題です。
例題「3次元 回転する球周りの流れ」の概要
この例題は、3次元の流体中に回転する球を配置し、その周囲を流れる流体の挙動をXLBで解析する例題です。例題スクリプトでは、3次元格子、D3Q27速度モデル、Warpバックエンド、FP32計算精度を用いてLBMによる流体解析を行います。球の形状にはSTLファイルを使用し、STLから球の表面形状を読み込んで計算領域内に配置します。
また、入口には流入速度を与える境界条件、出口には流出境界条件を設定し、球表面には回転する固体表面の速度を考慮した境界条件を設定しています。球はy軸を中心として回転し、流体が球の周囲を通過することで生じる流れ場を解析します。流入速度やレイノルズ数などの条件を設定し、球に作用する流体力を計算することで、回転する物体周りの流れの特徴を評価します。
さらに、計算ではLBMの格子単位（無次元単位） を用いて流速や時間発展を表現しており、prescribed_velocity = 0.04 のような値によって流入速度を設定します。このようなLBM単位と物理量との関係を理解することも、本例題を学ぶ重要なポイントです。
この例題では、単純な2次元流れの解析からさらに進み、3次元流れ、D3Q27モデル、STL形状の読み込み、回転する固体表面の境界条件、レイノルズ数による流れ条件の設定、流体力の評価、GPUを利用したWarpバックエンドによる計算という、XLBを3次元の工学的な流体解析へ応用するための基本的な機能を学ぶことができます。
例題の概念図を、以下に示します。
はい。先ほどのAhmed Body版の構成を維持しながら、現在の rotating_sphere_3d.py の実際の処理内容に合わせて全面的に置き換えると、以下の構成が適切です。
特に、rotating_sphere_3d.py はAhmed Body版とは異なり、マルチ解像度ではなく単一解像度、NeonではなくWarp、実験値との比較ではなく、抗力係数・揚力係数の時間変化を評価する例題です。現在のXLB main のスクリプト内容とも照合しています。:chatgpt-content-reference{index="0"}
例題「3次元 回転する球周りの流れ」の解析
以下では、スクリプトのコードそのものは示さず、対象とする例題 rotating_sphere_3d.py が「何を計算しているのか」を、「物理モデル → 3次元LBM → 計算格子 → 物理条件 → 球形状 → 回転境界条件 → 抗力・揚力 → 出力」 の順に説明します。
この例題は、XLBのWarpバックエンドを利用して、回転する球の周囲を流れる3次元流れをLBMで解析するものです。
以下の説明では、【06:例題multires_windtunnel_3d.pyの分析】の内容を用いています。詳しく知りたい方は、以下の【06】も参考にしてください。
https://qiita.com/ShibataRyoichi/items/5c0bf7d533a201e30bf3
1. この例題で何を解析するのか
この例題では、3次元の流体中に球を配置し、その球を回転させながら、球の周囲を流れる流体の挙動を解析します。球に対して流体が一定速度で流入すると、球の前方では流れがよどみ、球の側面では流れが加速します。また、球の後方では流れが剥離し、複雑な後流が形成されます。
さらに、この例題では球を単純に固定するのではなく、y軸を中心として回転させます。球の回転によって球表面の速度が変化するため、球周辺の流れも変化します。その結果、球には流れ方向の抗力だけでなく、回転による揚力も発生します。
したがって、この例題では、「球周辺の3次元速度場・球後方の後流・球に作用する抗力・球に作用する揚力・ 抗力係数 $C_D$・揚力係数 $C_L$」を、時間発展に沿って計算します。
2. 計算領域は3次元
rotating_sphere_3d.py では、3次元のLBM格子を使用します。計算格子は、「$320\times224\times224$」の3次元格子です。つまり、「x方向：320格子・y方向：224格子・z方向：224格子」となっています。
概念的には、風洞型の計算領域を考えます。x方向が主流方向であり、流体は左から右へ流れます。球は計算領域の内部に配置され、y軸を中心として回転します。
3. LBMでは何を計算しているのか
この例題では、通常のNavier-Stokes方程式を直接離散化するのではなく、格子ボルツマン法（LBM：Lattice Boltzmann Method） によって流れを計算します。
LBMでは、各格子点において、$f_i(\mathbf{x},t)$という粒子分布関数を計算します。
ここで、$i$：粒子の移動方向、$\mathbf{x}$：格子上の位置、$t$：時間、$f_i$：方向 $i$ の分布関数を表します。この分布関数を時間発展させることで、最終的に流体の密度と速度を求めます。
4. D3Q27とは
この例題では、3次元LBMの速度モデルとして、D3Q27 を使用します。D3Q27を分解すると、「D3：3次元・Q27：27個の離散速度」という意味です。
5. 衝突モデルによる時間発展
LBMでは、「衝突 → ストリーミング → 境界条件」 を繰り返して流れを時間発展させます。
この例題では、衝突モデルとして、KBC（Karlin-Bösch-Chikatamarla） を使用しています。
したがって、flow_past_sphere_3d.py などで使用される衝突モデルとは異なり、本例題ではKBCによって分布関数を時間発展させます。
6. 分布関数から密度と速度を求める
LBMで直接計算しているのは分布関数 $f_i$ ですが、流体解析で最終的に必要になるのは密度と速度です。
密度は、$\rho=\sum_i f_i$として求められます。速度は概念的に、$\rho\mathbf{u}=\sum_i f_i\mathbf{e}_i$から求められます。ここで、$\mathbf{e}_i$はD3Q27の各離散速度です。
この例題では、これらのマクロ量から、$u_x,\quad u_y,\quad u_z$および速度の大きさ$|\mathbf{u}|=\sqrt{u_x^2+u_y^2+u_z^2}$を求めます。
7. 流入速度はLBMの格子単位
この例題の流入速度は、$U=0.04$です。これはSI単位の$0.04\ {\rm m/s}$ではありません。LBMで使用する**格子単位（lattice units）**です。つまり、$U_{\mathrm{LBM}}=0.04$という無次元化された速度を与えています。
したがって、この値だけから実際のm/sに変換することはできません。物理単位へ対応させる場合には、格子間隔と時間刻みを対応させる必要があります。この例題では、LBM上の速度 $0.04$ と球の直径、レイノルズ数などの関係から流れの条件を設定します。
8. レイノルズ数は200
この例題では、$Re=200$を設定しています。レイノルズ数は、$Re=\frac{UD}{\nu}$で定義されます。ここで、$U$：代表速度、$D$：球の直径、$nu$：動粘性係数、です。
したがって、本例題では、$Re=200$となるようにLBM上の粘性係数を設定しています。
コード内部では、$\nu=\frac{UD}{Re}$という関係から粘性を決定し、その粘性からLBMの緩和パラメータを求めています。
9. 球の形状をSTLから読み込む
この例題では、球をPythonの数式で生成するのではなく、STLファイルから読み込みます。
概念的には、「sphere.stl → 三角形表面 → 球の3次元形状 → LBM格子上に配置 → 固体境界」という流れになります。
この方法を使用することで、単純な球だけでなく、将来的にはCADなどから作成した複雑な3次元形状も同じ考え方で扱うことができます。
10. STL形状を計算領域に配置する
読み込んだSTLの頂点座標は、そのまま使用するのではありません。まず球の座標を整理し、計算格子の単位に合わせてスケーリングします。概念的には、「$\text{STL座標}\rightarrow\text{LBM格子単位}\rightarrow\text{計算領域内へ移動}$」という変換を行います。そして、球が計算領域の中央付近に配置されるように座標を移動します。
11. 球はy軸を中心として回転する
この例題の重要な特徴が、球の回転です。球はy軸を中心として回転します。角速度ベクトルを、$\boldsymbol{\omega}=(0,\omega,0)$とすると、球表面のある点 $\mathbf{r}$ における表面速度は、$\mathbf{u}_s=\boldsymbol{\omega}\times\mathbf{r}$で表されます。
つまり、球の表面上の場所によって速度の方向と大きさが変わります。この回転速度を境界条件として流体に与えることで、回転球周りの流れを計算します。
12. 回転速度の無次元パラメータ
球の回転は、単純に角速度を指定するのではなく、流入速度と球径に対して無次元化された回転速度を用いて設定されています。本例題では、$\Omega_{\mathrm{nondim}}=-0.2$が使用されています。
また、実際のLBM計算で使用する角速度は、流入速度 $U$と球径 $D$ を利用して、$\omega=\frac{2U\Omega_{\mathrm{nondim}}}{D}$という関係から設定されます。
したがって、この例題では、「流入速度 + 球径 + レイノルズ数 + 回転速度」 によって回転球周りの流れ条件を決めています。
13. 入口には速度境界条件を与える
計算領域の左側が入口です。入口では、x方向に一定速度の流体が流入するようにします。速度は、$\mathbf{u}_{in}=(U,0,0)$です。本例題では、$U=0.04$を使用します。
入口の境界条件には、Regularized Boundary Conditionを使用します。
14. 出口ではDo Nothing境界条件を使用する
計算領域の右側は出口です。出口では流速を固定するのではなく、Do Nothing Boundary Conditionを使用します。これによって、球の後方で形成された流れが計算領域の外へ流出できるようにします。
15. 上下面・前後面は周期境界として扱われる
本例題では、上面・下面・前面・背面について明示的な境界条件を割り当てていません。そのため、XLBではこれらの面が**周期境界（Periodic Boundary）**として扱われます。
つまり、解析空間の上下や、z方向の前後が周期的につながるような計算領域になります。これは、無限に広がる流体中の球を有限の計算領域で近似するための設定として理解できます。
16. 球表面の境界条件
球表面には、流体と固体の間の相互作用を表現する境界条件を設定します。現在のスクリプトでは、Hybrid Boundary Conditionを使用しています。その内部では、Nonequilibrium Regularizedを境界処理として使用しています。
さらに、STL形状との距離情報を利用し、「STL表面 → Voxelization → LBM格子上の固体境界 → Hybrid Boundary Condition」という処理を行います。STLのボクセル化には、Ray方式が使用されています。
以上の、13,14,15,16の境界条件の内容を図示してまとめます。
17. 球表面に回転速度を与える
球表面では、単純な「速度0」の固定壁条件にはしません。球の回転によって発生する表面速度、$\mathbf{u}_s=\boldsymbol{\omega}\times\mathbf{r}$を境界条件として使用します。
したがって、以下の違いがあります。この点が、単純な flow_past_sphere_3d.py との大きな違いです。
固定球：表面速度 = 0 → 通常の球周り流れ
回転球：表面速度 = Ω × r → 回転によって流れが変化 → 揚力が発生
18. 球後方には後流が形成される
流体が球に衝突すると、球の前方では流れが分岐し、球の側面を通過します。その後、球の後方では流れが剥離し、渦を伴う後流が形成されます。さらに球が回転しているため、球表面の速度分布が上下方向の流れを変化させます。
その結果、「速度の低下・流れの剥離・渦構造・非対称な速度分布」などが発生します。これが回転球周りの流れの重要な特徴です。
19. 分布関数を時間発展させる
LBM計算では、初期状態から時間ステップを繰り返します。概念的には、以下の処理を繰り返します。
初期状態 → Collision → Streaming → Boundary Condition
→ 次の時間ステップ → Collision → Streaming → Boundary Condition → ……
本例題では最大、$100000$ステップまで計算します。また、一定間隔で計算状況を表示し、一定間隔で流れ場と流体力を評価します。
20. 流体の速度場を計算する
計算された分布関数から、$u_x,\quad u_y,\quad u_z$を求めます。さらに、$|\mathbf{u}|=\sqrt{u_x^2+u_y^2+u_z^2}$として速度の大きさを計算します。
この速度の大きさを利用すると、「球 ● → 球周辺の加速 → 球後方の速度低下 → 後流・渦構造」を可視化できます。本例題では、計算領域中央の断面における速度分布を画像として保存します。
21. 球に作用する力を計算する
この例題では、流れ場だけではなく、球に作用する力も計算します。そのために、Momentum Transfer（運動量移動法） を使用します。流体と球表面との間で交換される運動量を利用して、球に作用する流体力を求めます。
球に作用する力を、$\mathbf{F}=(F_x,F_y,F_z)$とすると、$F_x$：流れ方向の力、$F_y$：y方向の力、$f_z$：z方向の力、となります。
22. 抗力と揚力
この例題では、$F_D=F_x$を抗力として扱います。一方、$F_L=F_z$を揚力として扱います。回転していない球では上下方向の力は対称になりやすいですが、球が回転すると流れが非対称になるため、揚力が発生します。これは回転球に特徴的なMagnus効果を理解するためにも重要です。
23. 抗力係数を計算する
抗力そのものは、速度や球の大きさによって変化します。そこで、無次元化した抗力係数、$C_D$を使用します。本例題では、$C_D=\frac{2F_D}{U^2A}$として計算します。ここで、$F_D$：抗力、$U$：LBM上の流入速度、$A$：球の投影面積、です。球の投影面積は、$A=\frac{\pi D^2}{4}$です。
したがって、「球表面 → Momentum Transfer → 抗力 FD → 投影面積 A と速度 U で無次元化 → 抗力係数 CD」という流れになります。
24. 揚力係数を計算する
揚力についても同様に無次元化します。$C_L=\frac{2F_L}{U^2A}$です。したがって、本例題から、$C_D,\quad C_L$という2つの重要な流体力学的指標を得ることができます。特に回転球では、回転によって揚力がどのように発生するかを (C_L) によって評価できます。
25. 抗力・揚力を時間変化として記録する
この例題では、最終的な1個の (C_D) や (C_L) を求めて終了するのではありません。計算途中でも一定間隔で、$C_D(t)$、および$C_L(t)$を計算します。そして、「時間ステップ → CD, CLを計算 → 時系列データとして保存 → 移動平均を計算 → CD・CLの時間変化を可視化」をします。さらに、複数の移動平均幅を用いることで、瞬間的な変動と長時間平均の両方を見ることができます。
26. 計算結果を画像として出力する
本例題では、計算途中および計算終了時に流れ場を後処理します。主な処理は、以下の通りです。
- 分布関数から密度・速度を計算
- 速度の大きさを計算
- 球に作用する力を計算
- $C_D$ を計算
- $C_L$ を計算
- 速度分布画像を保存
- $C_D、C_L$ の時系列グラフを保存
つまり、単にLBM計算を実行するだけではなく、流れ場 + 流体力 + グラフまで自動的に評価する例題になっています。
27. VTKによる3次元結果出力
この例題にはVTK形式で流れ場を保存するための処理も用意されています。ただし、現在のスクリプトでは、**VTK保存処理はコメントアウトされています。
**したがって、標準状態では主に速度分布の画像と $_D、C_L$ のグラフを出力します。必要に応じてVTK出力を有効化すれば、ParaViewなどを利用して3次元流れ場を可視化することもできます。
28. WarpバックエンドによるGPU計算
この例題では、ComputeBackend.WARPを使用しています。したがって、「rotating_sphere_3d.py → XLB → Warp Backend → CUDA → NVIDIA GPU」という計算経路になります。今回の環境では、RTX A4500が認識されており、WarpによってGPU計算を実行できます。このため、この例題は、「XLB → Warp → CUDA → GPU」 というXLBのGPU計算環境を確認するための代表的な3次元例題としても利用できます。
29. この例題で最終的に得られるもの
このプログラムの結果を整理すると、大きく3種類に分けられます。
① 流れ場：$u_x,\quad u_y,\quad u_z,\quad |\mathbf{u}|$ です。
② 球に作用する流体力：$F_D,\quad F_L$ です。
③ 無次元流体力係数：$C_D,\quad C_L$ です。
全体を整理すると、以下のような構成となります。
rotating_sphere_3d.py
┌──────────────┼──────────────┐
流れ場 流体力 係数
velocity Drag / Lift CD / CL
速度画像 Momentum Transfer 時系列グラフ
30. この例題の本質
この rotating_sphere_3d.py は、XLBを理解するうえで、基本的な3次元物体周り流れからさらに一歩進んだ例題です。単純な球周りの流れだけではなく、「球を回転させることで、固体表面の速度を境界条件として流体に与え、その結果として発生する抗力・揚力を計算する」 ことが本例題の中心です。
全体を整理すると、以下の流れになります。
$\text{STL}\rightarrow\text{3D LBM格子} \rightarrow f_i\rightarrow\text{KBC衝突}\rightarrow\text{Streaming}\rightarrow\text{回転球境界}\rightarrow(\rho,\mathbf{u})\rightarrow\mathbf{F}\rightarrow(C_D,C_L)$
つまり、この例題を理解すると、「XLBでは、STLで定義された3次元形状を計算領域に配置し、WarpによるGPU上のLBM計算によってその周囲の流れを求め、さらに物体表面との運動量交換から抗力・揚力を評価できる」 ということが理解できます。
特に重要なのは、以下の連の流れです。
① D3Q27による3次元LBM
② KBCによる時間発展
③ STLによる球形状の取り込み
④ y軸周りの球の回転
⑤ 回転表面速度を考慮した境界条件
⑥ Momentum Transferによる流体力計算
⑦ $C_D, C_L$ の時間変化
⑧ Warp + CUDA + GPUによる計算
31. flow_past_sphere_3d.py から何が発展したのか
この例題を、通常の固定球周り流れと比較すると理解しやすくなります。
flow_past_sphere_3d.py：流入流れ → ● 球 → 後流を計算
rotating_sphere_3d.py：流入流れ → ↻ ● 回転する球 → 表面速度を考慮
→ 3次元流れを計算 → 運動量交換を計算 → Drag / Lift → CD / CL
へ発展しています。したがって、この例題は、「XLBで3次元流れを計算する」だけでなく、「動く・回転する物体と流体との相互作用を計算し、物体に作用する力まで評価する」 ための例題として位置づけることができます。
なお、現在のXLB main の実装では、冒頭のdocstringに「halfway bounce-back」とありますが、実際の球境界には HybridBC（nonequilibrium_regularized、mesh distance、RAY voxelization、回転速度profile）が設定されています。記事では実際のコードで使われている境界条件を優先して「HybridBC」と記述するのが正確です。
例題スクリプトの具体的な解説
元のスクリプトは変更せず、学習用の詳細コメントだけを追加し、①〜⑫の12ブロックに整理しています。元スクリプトが D3Q27、Warp、STL、回転速度プロファイル、HybridBC、MomentumTransfer、時間発展を含む構成であることもそのまま維持しています。
12ブロックの構成
- ライブラリの読み込み
- 計算条件・物理条件の設定
- XLBの初期化と計算格子の生成
- STL形状の読み込み
- STL形状の座標変換と配置
- 回転する球の表面速度プロファイル
- 境界条件の設定
- LBM時間発展器・初期状態・フィールドの準備
- 流体力係数を可視化する補助関数
- 計算結果の後処理
- 流体力計算とマクロ量計算の準備
- LBM時間発展ループと結果出力
特に⑦では、実際のスクリプトにある RegularizedBC、DoNothingBC、HybridBC、RAY方式、回転速度profile をそのまま説明しています。
スクリプトファイル：
# =============================================================================
# ① ライブラリの読み込み
# -----------------------------------------------------------------------------
# このブロックでは、XLBによるLBM計算、WarpによるGPU計算、STL形状の読み込み、
# NumPy/JAXによる配列処理、VTK・画像出力、グラフ作成などに必要なライブラリを
# 読み込みます。
#
# 特に重要なのは、XLBを中心に、
# XLB → ComputeBackend.WARP → NVIDIA Warp → CUDA → GPU
# という計算経路を構成する点です。
#
# なお、このブロックでは計算そのものはまだ開始していません。
# =============================================================================
"""
Rotating sphere 3-D example (single-resolution).
Simulates flow past a sphere rotating about the y-axis using the
halfway bounce-back BC with a prescribed rotational-velocity profile.
Computes drag and lift coefficients over time and saves VTK snapshots.
"""
import xlb
import trimesh
import time
import warp as wp
import numpy as np
import jax.numpy as jnp
from typing import Any
from xlb.compute_backend import ComputeBackend
from xlb.precision_policy import PrecisionPolicy
from xlb.grid import grid_factory
from xlb.operator.stepper import IncompressibleNavierStokesStepper
from xlb.operator.boundary_condition import (
HalfwayBounceBackBC,
FullwayBounceBackBC,
RegularizedBC,
DoNothingBC,
HybridBC,
)
from xlb.operator.force.momentum_transfer import MomentumTransfer
from xlb.operator.macroscopic import Macroscopic
from xlb.utils import save_fields_vtk, save_image, warp_array_to_jax
import matplotlib.pyplot as plt
from xlb.operator.equilibrium import QuadraticEquilibrium
from xlb.operator import Operator
from xlb.velocity_set.velocity_set import VelocitySet
from xlb.operator.boundary_masker import MeshVoxelizationMethod
# =============================================================================
# ② 計算条件・物理条件の設定
# -----------------------------------------------------------------------------
# ここでは、3次元計算格子の大きさ、計算バックエンド、精度、D3Q27速度モデル、
# 流入速度、計算ステップ数、Reynolds数など、解析全体の基本条件を設定します。
#
# diam は球の代表長さをLBM格子単位で表した値です。
# grid_shape は球の大きさを基準として、x方向10倍、y/z方向7倍の計算領域を
# 構成します。
#
# Re = U D / ν から粘性係数 visc を求め、さらにLBMの緩和パラメータ omega を
# 計算します。
#
# また、rot_rate_nondim は球の回転を無次元化して表した値であり、これを
# 流入速度と球径から実際のLBM計算で使う回転速度 rot_rate に変換します。
# =============================================================================
# -------------------------- Simulation Setup --------------------------
# Grid parameters
wp.clear_kernel_cache()
diam = 32
grid_size_x, grid_size_y, grid_size_z = 10 * diam, 7 * diam, 7 * diam
grid_shape = (grid_size_x, grid_size_y, grid_size_z)
# Simulation Configuration
compute_backend = ComputeBackend.WARP
precision_policy = PrecisionPolicy.FP32FP32
velocity_set = xlb.velocity_set.D3Q27(precision_policy=precision_policy, compute_backend=compute_backend)
wind_speed = 0.04
num_steps = 100000
print_interval = 1000
post_process_interval = 1000
# Physical Parameters
Re = 200.0
visc = wind_speed * diam / Re
omega = 1.0 / (3.0 * visc + 0.5)
# Rotational speed parameters (see [1] which discusses the problem in terms of 2 non-dimensional parameters: Re and Omega)
# [1] J. Fluid Mech. (2016), vol. 807, pp. 62–86. c© Cambridge University Press 2016 doi:10.1017/jfm.2016.596
# \Omega = \omega * D / (2 U_\infty) where Omega is non-dimensional and omega is dimensional.
rot_rate_nondim = -0.2
rot_rate = 2.0 * wind_speed * rot_rate_nondim / diam
# Print simulation info
print("\n" + "=" * 50 + "\n")
print("Simulation Configuration:")
print(f"Grid size: {grid_size_x} x {grid_size_y} x {grid_size_z}")
print(f"Backend: {compute_backend}")
print(f"Velocity set: {velocity_set}")
print(f"Precision policy: {precision_policy}")
print(f"Prescribed velocity: {wind_speed}")
print(f"Reynolds number: {Re}")
print(f"Max iterations: {num_steps}")
print("\n" + "=" * 50 + "\n")
# Initialize XLB
# =============================================================================
# ③ XLBの初期化と計算格子の生成
# -----------------------------------------------------------------------------
# ここでXLBを指定したD3Q27・Warp・FP32設定で初期化します。
#
# その後、grid_factory() によって3次元LBM格子を生成します。
# さらに、計算領域の各面をインデックスとして取得し、
# ・入口（left）
# ・出口（right）
# ・上下面、前後面（walls）
# を後の境界条件設定で利用できるようにします。
#
# walls に明示的な境界条件を設定していない点は、後のブロックで重要になります。
# =============================================================================
xlb.init(
velocity_set=velocity_set,
default_backend=compute_backend,
default_precision_policy=precision_policy,
)
# Create Grid
grid = grid_factory(grid_shape, compute_backend=compute_backend)
# Bounding box indices
box = grid.bounding_box_indices()
box_no_edge = grid.bounding_box_indices(remove_edges=True)
inlet = box_no_edge["left"]
outlet = box["right"]
walls = [box["bottom"][i] + box["top"][i] + box["front"][i] + box["back"][i] for i in range(velocity_set.d)]
walls = np.unique(np.array(walls), axis=-1).tolist()
# Load the mesh (replace with your own mesh)
# =============================================================================
# ④ STL形状の読み込み
# -----------------------------------------------------------------------------
# 解析対象となる球をSTLファイルから読み込みます。
# STLは三角形の集合として球表面を表しており、mesh.vertices に各頂点座標が
# 格納されます。
#
# この例題では、単純な数式で球を生成するのではなく、STL形状をXLBのLBM計算
# 領域へ取り込むことで、実際の3次元形状を境界として扱う方法を示しています。
# =============================================================================
stl_filename = "../stl-files/sphere.stl"
mesh = trimesh.load_mesh(stl_filename, process=False)
mesh_vertices = mesh.vertices
# Transform the mesh points to be located in the right position in the wind tunnel
# =============================================================================
# ⑤ STL形状の座標変換と計算領域への配置
# -----------------------------------------------------------------------------
# 読み込んだSTLの座標をLBM格子で扱えるように変換します。
#
# まず最小座標を引いて形状の基準位置をそろえ、STLの物理的な大きさと
# LBM格子上で必要な大きさから dx を求めます。
#
# その後、球を計算領域の所定位置へ移動し、最終的な球表面座標 sphere を
# 作ります。
#
# 最後に、LBM上の球径 diam と球の投影面積 sphere_cross_section を求めます。
# 後の抗力係数・揚力係数の計算で、この投影面積が基準面積として使われます。
# =============================================================================
mesh_vertices -= mesh_vertices.min(axis=0)
mesh_extents = mesh_vertices.max(axis=0)
length_phys_unit = mesh_extents.max()
length_lbm_unit = grid_shape[1] / 7
dx = length_phys_unit / length_lbm_unit
mesh_vertices = mesh_vertices / dx
shift = np.array([grid_shape[0] / 3, (grid_shape[1] - mesh_extents[1] / dx) / 2, (grid_shape[2] - mesh_extents[2] / dx) / 2])
sphere = mesh_vertices + shift
diam = np.max(sphere.max(axis=0) - sphere.min(axis=0))
sphere_cross_section = np.pi * diam**2 / 4.0
# Define rotating boundary profile
# =============================================================================
# ⑥ 回転する球の表面速度プロファイル
# -----------------------------------------------------------------------------
# 球がy軸を中心に回転するため、球表面の各格子点に与える速度を定義します。
#
# 角速度ベクトルは、
# ω = (0, rot_rate, 0)
# とし、球中心から表面点までの位置ベクトルを r とすると、
#
# u_s = ω × r
#
# によって、その位置における球表面速度を求めます。
#
# この関数はWarpのデバイス側関数として定義され、球表面の境界条件を計算する
# 際に利用されます。
# =============================================================================
def bc_profile():
"""Build a Warp function returning the rotational wall velocity at a voxel."""
dtype = precision_policy.compute_precision.wp_dtype
_u_vec = wp.vec(velocity_set.d, dtype=dtype)
angular_velocity = _u_vec(0.0, rot_rate, 0.0)
origin_np = shift + diam / 2
origin_wp = _u_vec(origin_np[0], origin_np[1], origin_np[2])
@wp.func
def bc_profile_warp(index: wp.vec3i):
x = dtype(index[0])
y = dtype(index[1])
z = dtype(index[2])
surface_coord = _u_vec(x, y, z) - origin_wp
return wp.cross(angular_velocity, surface_coord)
return bc_profile_warp
# Define boundary conditions
# =============================================================================
# ⑦ 境界条件の設定
# -----------------------------------------------------------------------------
# ここでは入口・出口・球表面の境界条件を定義します。
#
# 入口：
# RegularizedBC により、x方向へ wind_speed = 0.04 の一定速度を与えます。
#
# 出口：
# DoNothingBC により、流体が計算領域外へ流出できるようにします。
#
# 球：
# HybridBC を使用し、nonequilibrium_regularized を境界処理として指定します。
# STL形状をRAY方式でボクセル化し、mesh distance を利用して固体境界を
# 定義します。さらに、⑥で定義した回転速度プロファイルを適用します。
#
# 上下面・前後面には境界条件を割り当てていません。このXLBの設定により、
# それらは周期境界として扱われます。
# =============================================================================
bc_left = RegularizedBC("velocity", prescribed_value=(wind_speed, 0.0, 0.0), indices=inlet)
bc_do_nothing = DoNothingBC(indices=outlet)
# bc_sphere = HalfwayBounceBackBC(mesh_vertices=sphere, voxelization_method="ray", profile=bc_profile())
bc_sphere = HybridBC(
bc_method="nonequilibrium_regularized",
mesh_vertices=sphere,
use_mesh_distance=True,
voxelization_method=MeshVoxelizationMethod("RAY"),
profile=bc_profile(),
)
# Not assining BC for walls makes them periodic.
boundary_conditions = [bc_left, bc_do_nothing, bc_sphere]
# Setup Stepper
# =============================================================================
# ⑧ LBM時間発展器・初期状態・フィールドの準備
# -----------------------------------------------------------------------------
# IncompressibleNavierStokesStepper を生成し、計算格子・境界条件・衝突モデル
# をまとめてLBMの時間発展処理へ渡します。
#
# 衝突モデルには KBC を指定しています。
#
# その後、CustomInitializer により初期速度を一様流として設定し、
# stepper.prepare_fields() によって分布関数 f_0, f_1 と、境界条件・欠損情報を
# 表すマスクを準備します。
#
# ここまでで、時間発展を開始するための計算状態が整います。
# =============================================================================
stepper = IncompressibleNavierStokesStepper(
grid=grid,
boundary_conditions=boundary_conditions,
collision_type="KBC",
)
# Make initializer operator
from xlb.helper.initializers import CustomInitializer
initializer = CustomInitializer(
bc_id=bc_do_nothing.id,
constant_velocity_vector=(wind_speed, 0.0, 0.0),
velocity_set=velocity_set,
precision_policy=precision_policy,
compute_backend=compute_backend,
)
# Prepare Fields
f_0, f_1, bc_mask, missing_mask = stepper.prepare_fields(initializer=initializer)
# -------------------------- Helper Functions --------------------------
# =============================================================================
# ⑨ 流体力係数を可視化する補助関数
# -----------------------------------------------------------------------------
# この関数は、計算中に蓄積した抗力係数または揚力係数の時系列をグラフ化します。
#
# 生の係数だけでなく、10、100、1,000、10,000、100,000ステップを窓とする
# 移動平均を計算することで、瞬間的な変動と長時間的な傾向を確認できます。
#
# prefix が drag の場合は抗力係数、lift の場合は揚力係数の画像として保存します。
# =============================================================================
def plot_coefficient(time_steps, coefficients, prefix="drag"):
"""
Plot the drag coefficient with various moving averages.
Args:
time_steps (list): List of time steps.
coefficients (list): List of force coefficients.
"""
# Convert lists to numpy arrays for processing
time_steps_np = np.array(time_steps)
coefficients_np = np.array(coefficients)
# Define moving average windows
windows = [10, 100, 1000, 10000, 100000]
labels = ["MA 10", "MA 100", "MA 1,000", "MA 10,000", "MA 100,000"]
plt.figure(figsize=(12, 8))
plt.plot(time_steps_np, coefficients_np, label="Raw", alpha=0.5)
for window, label in zip(windows, labels):
if len(coefficients_np) >= window:
ma = np.convolve(coefficients_np, np.ones(window) / window, mode="valid")
plt.plot(time_steps_np[window - 1 :], ma, label=label)
plt.ylim(-1.0, 1.0)
plt.legend()
plt.xlabel("Time step")
plt.ylabel("Drag coefficient")
plt.title("Drag Coefficient Over Time with Moving Averages")
plt.savefig(prefix + "_ma.png")
plt.close()
def post_process(
# =============================================================================
# ⑩ 計算結果の後処理
# -----------------------------------------------------------------------------
# post_process() は、一定ステップごとに流れ場と流体力を評価する中心的な関数です。
#
# まず分布関数をJAX配列へ変換可能な形にし、Macroscopic演算子から密度 rho と
# 速度 u を計算します。
#
# 次に境界セルを除外して速度の大きさ |u| を求め、中央y断面の速度分布画像を
# 保存します。
#
# さらに MomentumTransfer により球表面で交換された運動量から抗力・揚力を求め、
#
# CD = 2 FD / (U^2 A)
# CL = 2 FL / (U^2 A)
#
# として無次元化します。
#
# 求めた CD、CL、時刻をリストへ保存し、移動平均グラフも更新します。
# =============================================================================
step,
f_0,
f_1,
grid_shape,
macro,
momentum_transfer,
missing_mask,
bc_mask,
wind_speed,
car_cross_section,
drag_coefficients,
lift_coefficients,
time_steps,
):
"""Compute macroscopic fields, force coefficients, and save VTK output."""
"""
Post-process simulation data: save fields, compute forces, and plot drag coefficient.
Args:
step (int): Current time step.
f_current: Current distribution function.
grid_shape (tuple): Shape of the grid.
macro: Macroscopic operator object.
momentum_transfer: MomentumTransfer operator object.
missing_mask: Missing mask from stepper.
bc_mask: Boundary condition mask from stepper.
wind_speed (float): Prescribed wind speed.
car_cross_section (float): Cross-sectional area of the car.
drag_coefficients (list): List to store drag coefficients.
lift_coefficients (list): List to store lift coefficients.
time_steps (list): List to store time steps.
"""
wp.synchronize()
# Convert to JAX array if necessary
if not isinstance(f_0, jnp.ndarray):
f_0_jax = warp_array_to_jax(f_0)
else:
f_0_jax = f_0
# Compute macroscopic quantities
rho, u = macro(f_0_jax)
# Remove boundary cells
u = u[:, 1:-1, 1:-1, 1:-1]
u_magnitude = jnp.sqrt(u[0] ** 2 + u[1] ** 2 + u[2] ** 2)
fields = {"ux": u[0], "uy": u[1], "uz": u[2], "u_magnitude": u_magnitude}
# Save fields in VTK format
# save_fields_vtk(fields, timestep=step)
# Save the u_magnitude slice at the mid y-plane
mid_y = grid_shape[1] // 2
save_image(fields["u_magnitude"][:, mid_y, :], timestep=step)
# Compute lift and drag
boundary_force = momentum_transfer(f_0, f_1, bc_mask, missing_mask)
drag = boundary_force[0] # x-direction
lift = boundary_force[2]
cd = 2.0 * drag / (wind_speed**2 * car_cross_section)
cl = 2.0 * lift / (wind_speed**2 * car_cross_section)
print(f"CD={cd}, CL={cl}")
drag_coefficients.append(cd)
lift_coefficients.append(cl)
time_steps.append(step)
# Plot drag coefficient
plot_coefficient(time_steps, drag_coefficients, prefix="drag")
plot_coefficient(time_steps, lift_coefficients, prefix="lift")
# Setup Momentum Transfer for Force Calculation
# =============================================================================
# ⑪ 流体力計算とマクロ量計算の準備
# -----------------------------------------------------------------------------
# 球の境界条件を対象として MomentumTransfer を生成します。
# これにより、球表面で流体と固体の間に交換された運動量から流体力を評価できます。
#
# また、Macroscopic 演算子をJAXバックエンドで準備します。
# ここでは、LBMの時間発展そのものはWarpで行いながら、密度・速度などの
# マクロ量の計算にはJAXを利用する構成になっています。
#
# 最後に、抗力係数・揚力係数・時間ステップを保存するためのリストを初期化します。
# =============================================================================
bc_car = boundary_conditions[-1]
momentum_transfer = MomentumTransfer(bc_car, compute_backend=compute_backend)
# Define Macroscopic Calculation
macro = Macroscopic(
compute_backend=ComputeBackend.JAX,
precision_policy=precision_policy,
velocity_set=xlb.velocity_set.D3Q27(precision_policy=precision_policy, compute_backend=ComputeBackend.JAX),
)
# Initialize Lists to Store Coefficients and Time Steps
time_steps = []
drag_coefficients = []
lift_coefficients = []
# -------------------------- Simulation Loop --------------------------
# =============================================================================
# ⑫ LBM時間発展ループと結果出力
# -----------------------------------------------------------------------------
# ここから実際のLBM計算を開始します。
#
# 各ステップで、
# ① stepper による Collision / Streaming / Boundary の時間発展
# ② 分布関数バッファの入れ替え
# ③ 一定間隔で計算進捗を表示
# ④ 一定間隔で post_process() を実行
# を繰り返します。
#
# post_process() では、速度場、抗力、揚力、CD、CLを評価し、画像やグラフを
# 更新します。
#
# num_steps = 100000 なので、最大10万ステップまで計算します。
# 最後に "Simulation completed successfully." を表示して計算を終了します。
# =============================================================================
start_time = time.time()
for step in range(num_steps):
# Perform simulation step
f_0, f_1 = stepper(f_0, f_1, bc_mask, missing_mask, omega, step)
f_0, f_1 = f_1, f_0 # Swap the buffers
# Print progress at intervals
if step % print_interval == 0:
elapsed_time = time.time() - start_time
print(f"Iteration: {step}/{num_steps} | Time elapsed: {elapsed_time:.2f}s")
start_time = time.time()
# Post-process at intervals and final step
if (step % post_process_interval == 0) or (step == num_steps - 1):
post_process(
step,
f_0,
f_1,
grid_shape,
macro,
momentum_transfer,
missing_mask,
bc_mask,
wind_speed,
sphere_cross_section,
drag_coefficients,
lift_coefficients,
time_steps,
)
print("Simulation completed successfully.")
例題スクリプトの結果の分析
以下の記事【07:Warp導入方法と動作確認】で解説した手順を用いて、環境を準備して例題を実行すると、多数の解析結果が作られます。画面に表示されるログについては、【07】で概要を確認しました。
ここでは、ディレクトリ「ROTS3D」に保存されたファイルについて分析します。ファイルが多数ありますが、現在の状態では、CD（drag_ma.png）とCL（lift_ma.png）のグラフのpng画像ファイルと、速度分布の結果画像pngファイルが、1000ステップごとに、99999ステップまでの101段階で保存されます。
CD（抗力係数）とCL（揚力係数）のグラフ
いずれのグラフも、「Raw」が解析データそのもの、「MA 10」が直近10ステップの移動平均、「MA 100」が直近100ステップの移動平均、です。
移動平均により、細かな時間変化ではなく、滑らかな変動を見ることができます。
まず、抗力係数CD（drag_ma.png）のグラフは以下の通りです。
グラフのステップ０から暫く値が上下に動いていますが、これは、球後方の後流に渦が形成され、流体力が時間的に変化するためです。
次に、揚力係数CL（lift_ma.png）のグラフは以下の通りです。
速度分布の結果画像の分析
静止状態直後1000から、最初の3000ステップまでは、以下の状態です。流れが進んでいる様子が分かります。
途中の41000から43000までは、大きな変化はないようです。
最後の98000から99999までも、大きな変化はないようです。
このように、速度分布が「すぐに変化しなくなったように見える」のは、計算開始時から一様流が設定されており、球の周囲だけが短時間で流れ場を調整するためです。ただし、実際に流れ全体が完全に定常になったとは限りません。
【球の周囲の変化は比較的速い】
この例では、速度は$U=0.04$、球径は$D=32$です。したがって、流体が球１個分の距離を進む時間は、$t_D=\frac{D}{U}=\frac{32}{0.04}=800$ ステップ程度です。
つまり、約800ステップで流体が球径分だけ移動することになります。スクリプトでは1000ステップごとに速度場を出力しています。
そのため、画像を見るとすぐに、球周辺の基本的な流れ構造がかなり形成されているように見えます。
【計算領域全体が落ち着くにはもっと時間が必要】
計算領域のx方向は、$L_x=320$です。流速が0.04なので、入口から出口まで流体が通過する時間は、$t_{\mathrm{domain}}=\frac{320}{0.04}=8000$ ステップです。
さらに、球は入口から約$\frac{320}{3}\approx107$付近に置かれています。したがって、球から出口までの距離は約$320-107=213$格子程度です。
これを流体が通過するには、$\frac{213}{0.04}\approx5330$ステップ程度かかります。
つまり、1000～2000ステップ程度では、球周辺の局所的な流れは形成されても、計算領域全体の流れが完全に落ち着いたとは言えません。
以上により、**「速度分布がすぐ変化しなくなる」＝「計算がすぐ定常状態になった」**とは考えない方がよいです。
この例題では、
① 約800ステップで球径1個分の流れが進む
② 球周辺の基本的な流れ構造は比較的早く形成される
③ その後も後流や流体力には時間変動が残る可能性がある
④ 速度画像ではその変化が見えにくい
⑤ $C_D$、$C_L$ の時系列を見ることで変化を確認できる**
という理解が適切です。