この文書では、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計算を比較的容易に実装できます。
以下の【05:Neon導入方法とマルチ解像度例題の分析】では、バックエンドにNeonを GPUで利用する環境を構築し、例題multires_windtunnel_3d.pyの動作確認まで進めました。
この【06】では、【05】の手順で構築したNeonのGPU活用の環境をもとに、例題 multires_windtunnel_3d.pyを分析してゆきます。
https://github.com/Autodesk/XLB/blob/main/examples/cfd/multires_windtunnel_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 を分析します。この例題では、Ahmed Body周りの3次元流れを対象として、マルチ解像度LBM（Multi-Resolution LBM） による計算を行います。さらに、Neonバックエンドの利用、STL形状からの計算格子生成、物理単位とLBM単位の変換、抗力・揚力係数の計算、実験結果との比較など、XLBを実際の工学的な流体解析へ応用するための重要な機能が含まれています。
そのため、本例題は、flow_past_sphere_3d.py で学んだXLBの基本からさらに進み、複雑な3次元形状に対するマルチ解像度LBM解析の考え方と、XLBによる実用的な流体解析スクリプトの構成を学ぶための例題として位置づけることができます。
例題「3次元 マルチ解像度での自動車モデル周囲の流れ」の概要
この例題は、3次元の流体中に自動車形状を模擬したAhmed Bodyを配置し、その周囲を流れる空気の流れをXLBで解析する例題です。例題スクリプトでは、3次元格子、D3Q27モデル、Neonバックエンド、FP32計算精度を用い、マルチ解像度LBM（Multi-Resolution LBM） によって計算を行います。Ahmed Bodyの形状にはSTLファイルを使用し、車体周辺を細かい格子、遠方を粗い格子とすることで、計算精度と計算量のバランスを図っています。
また、入口には速度境界条件、出口には流出境界条件を設定し、車体表面にはSTL形状に基づく境界条件を設定しています。さらに、車体に作用する力を運動量移動法（Momentum Transfer）によって求め、抗力係数（Cd）と揚力係数（Cl）を計算します。計算結果はHDF5形式などで出力するとともに、車体後方の速度分布を実験データと比較することで、解析結果を評価します。
この例題では、単純な3次元物体周り流れの解析からさらに進み、複雑な形状、マルチ解像度格子、Neonバックエンド、流体力の評価、実験結果との比較という、XLBを工学的な流体解析に応用するための重要な機能を学ぶことができます。
例題の概念図を、以下に示します。
以下の説明では、スクリプトのコードそのものは示さず、対象とする例題 multires_windtunnel_3d.py が「何を計算しているのか」を、物理モデル → マルチ解像度LBM → 計算格子 → 境界条件 → Ahmed Body → 抗力・揚力 → 出力・実験値との比較の順に、文章と数式で詳しく説明します。
1. この例題で何を解析するのか
この例題では、3次元の流路の中に自動車形状を模擬したAhmed Bodyを配置し、その周囲を流れる空気の流れをXLBで解析します。
Ahmed Bodyは、自動車周りの流れを研究するために広く利用されている簡略化された車体モデルです。本例題では、25度のスラント角を持つAhmed Bodyを対象としています。
流体が車体に衝突すると、車体前方では流れがよどみ、車体上面や側面では流れが加速します。また、車体後方では流れが剥離し、複雑な後流（wake） が形成されます。
この例題では、単に速度や密度の分布を求めるだけではありません。車体に作用する抗力（Drag）と揚力（Lift）を計算し、さらに実験結果と速度分布を比較するところまで行います。
したがって、本例題は単純なLBMの計算例から一歩進んだ、工学的な3次元流体解析の例題として位置づけることができます。
2. 計算領域は3次元
この例題 multires_windtunnel_3d.py では、2次元の lid_driven_cavity_2d.py とは異なり、3次元のLBM格子を使用します。
さらに、flow_past_sphere_3d.py のように計算領域全体を同じ細かさの格子で計算するのではなく、複数の解像度を持つ計算格子を使用します。
概念的には、以下のように、車体周辺を細かく、車体から離れた領域を粗くします。
これは、車体周辺では流れの変化が大きく、詳細な計算が必要である一方、車体から遠い領域では比較的粗い格子でも流れを表現できるためです。
つまり、マルチ解像度LBMの基本的な考え方は、必要な場所には細かい格子を使い、必要性の低い場所では粗い格子を使うということです。
3. LBMでは何を計算しているのか
この例題では、通常のNavier-Stokes方程式を直接離散化するのではなく、格子ボルツマン法（LBM：Lattice Boltzmann Method） によって流れを計算します。
LBMでは、各格子点において流体の密度や速度を直接時間発展させるのではなく、$f_i(\mathbf{x},t)$という、粒子分布関数（distribution function） を計算します。
ここで、$i$：粒子の移動方向、$\mathbf{x}$：格子上の位置、$t$：時間、$f_i$：方向 (i) の粒子分布、を表します。今回の例題では、3次元LBMの速度モデルとしてD3Q27を使用しています。
4. D3Q27とは
D3Q27は、3次元LBMで使用される代表的な速度モデルです。名前を分解すると、「D3 → ３次元・Q27 → 27個の離散速度」という意味です。
D3Q27では、1つの格子点に対して、$f_0,f_1,\ldots,f_{26}$、という27個の分布関数を扱います。
flow_past_sphere_3d.py のD3Q19に対して、今回の例題ではD3Q27を使用している点が重要です。
5. 衝突モデルによる時間発展
LBMでは、基本的に、「衝突（Collision）→ 移流・ストリーミング（Streaming）→ 境界条件」、を繰り返すことで流れを時間発展させます。
今回の multires_windtunnel_3d.py では、シミュレーション初期化時に衝突モデルとしてKBCが指定されています。
概念的には、$f_i\rightarrow\text{Collision} \rightarrow f_i^{*}\rightarrow\text{Streaming}\rightarrow
f_i^{n+1}$ という流れです。
ここで、$f_i$：現在の分布関数、$f_i^{*}$：衝突後の分布関数、$f_i^{n+1}$：次の時間ステップの分布関数です。
この例題では、KBC衝突モデルによって分布関数を更新しながら、Ahmed Body周辺の流れを計算します。なお、flow_past_sphere_3d.py で使用していたBGKモデルとは異なり、今回の例題ではKBCが使用されている点が重要です。
6. 分布関数から密度と速度を求める
LBMで直接計算しているのは $f_i$ ですが、流体解析で最終的に必要となるのは、密度 $\rho$ と速度 $\mathbf{u}$ です。
密度は、$\rho=\sum_i f_i$ として求めることができます。また、速度は概念的に、$\rho\mathbf{u}=\sum_i f_i\mathbf{e}_i$ から求められます。ここで $\mathbf{e}_i$ はD3Q27の離散速度です。
本例題では、これらのマクロ量を利用して速度分布や密度分布を出力します。
7. 物理単位とLBM単位を変換する
今回の例題では、実際の空気流れを想定して、入口速度を $U_{\mathrm{physical}}=38.0\ {\rm m/s}$ としています。一方、LBMの計算では、$U_{\mathrm{LBM}}=0.05$ という格子単位の速度を使用します。
したがって、「現実の物理世界 → 38.0 m/s → UnitConvertor → LBM計算空間 → 0.05」という変換を行います。
また、最も細かい格子の物理的な大きさとして、$\Delta x_{\mathrm{fine}}=0.005\ {\rm m}$ を設定しています。
このように、現実の物理量をLBMで計算可能な格子単位へ変換することが、この例題の重要な準備になります。
8. Ahmed Bodyの形状をSTLから読み込む
今回の例題では、球のようにプログラムの中で単純な数式によって物体形状を作るのではありません。STLファイルで定義されたAhmed Bodyの形状を読み込みます。 使用するSTLは、25度のスラント角を持つAhmed Bodyです。以下のような流れになります。
STLファイル → Ahmed Bodyの三角形表面 → XLBの計算領域 → LBM格子上の固体境界
これによって、CADなどで作成した複雑な3次元形状を、XLBによるLBM解析に利用できます。
9. マルチ解像度の計算格子を生成する
STL形状を読み込んだ後、計算領域をマルチ解像度化します。本例題では、４段階の解像度レベルを持つ格子を生成します。
今回のスクリプトでは、最も細かい格子のサイズを0.005 mとして、そこから解像度を粗くしていきます。つまり、$\Delta x_{\mathrm{coarse}}=\Delta x_{\mathrm{fine}}2^{L-1}$ という関係で、解像度レベルが変化します。これによって、Ahmed Body近傍では細かい格子を使い、遠方では粗い格子を使うことができます。
10. 入口には速度境界条件を与える
計算領域の左側は入口です。入口では、空気が一定の速度で計算領域へ流入するように設定します。入口速度は、$\mathbf{u}_{in}=(U,0,0)$ とし、$U=38.0\ {\rm m/s}$ という物理速度に対応するLBM速度を使用します。
今回の例題では、入口にRegularized Boundary Conditionを使用しています。
11. 上面・下面・前面・背面の境界条件
流路の上面・下面・前面・背面についても境界条件を設定します。今回のスクリプトでは、これらの境界に対してNonequilibrium Regularizedを利用したHybrid Boundary Conditionを使用しています。
これによって、Ahmed Bodyを囲む3次元の流路全体を定義しています。
12. Ahmed Bodyは固体として扱われる
この例題の中心となるのがAhmed Bodyです。STLから読み込んだ車体表面を、LBMの計算格子上で固体境界として扱います。球の例題では、数式によって球の内部を判定していました。
一方、今回の例題では、以下のような処理を行います。これにより、複雑な3次元形状をLBM格子上に配置することができます。つまり、車体表面には、Hybrid Boundary Conditionが使用され、STL形状との距離情報も利用しています。
STL → 三角形表面 → Voxelization → LBM格子上の固体領域
13. 出口では流体を外へ流す
計算領域の右側は出口です。入口では速度を指定しますが、出口では速度を直接指定するのではなく、Do Nothing Boundary Conditionを使用します。概念的には、流体が計算領域から外へ流出できるようにします。
14. Ahmed Body後方には後流が形成される
車体前方では流れが車体を避けるように分かれます。車体上面では流れが加速し、車体後方の25度スラント部分では流れが剥離します。このような複雑な流れが形成され、車体後方の後流構造を計算できることが、この例題の重要なポイントです。
15. 車体に作用する抗力と揚力を計算する
この例題では、単に流れ場を計算して終了するのではありません。車体表面で流体と固体の間に生じる運動量交換から、車体に作用する流体力を計算します。XLBでは、MultiresMomentumTransferを利用しています。
車体に作用する力を、$\mathbf{F}=(F_x,F_y,F_z)$ とします。今回の例題では、流れ方向である $x$ 方向の力を抗力として扱い、$F_D=F_x$ とします。また、鉛直方向である $z$ 方向の力を揚力として、$F_L=F_z$ として扱います。
16. 抗力係数
抗力そのものだけでは、流速や物体の大きさが変わった場合に比較しにくくなります。そこで、抗力を無次元化した抗力係数 (C_D) を使用します。
本例題では、$C_D=\frac{2F_D}{U_{\mathrm{LBM}}^2 A}$ という形で計算します。
ここで、$F_D$：抗力、$U_{\mathrm{LBM}}$：LBM上の代表速度、$A$：基準面積、です。
本例題では、車体の投影面積をLBM格子上から求め、それを基準面積として利用しています。
したがって、以下のような流れになります。
車体表面 → 運動量交換 → 抗力 FD → 基準面積・速度で無次元化 → 抗力係数 CD
17. 揚力係数
揚力についても同様に、$C_L=\frac{2F_L}{U_{\mathrm{LBM}}^2 A}$ として無次元化します。
したがって、本例題から、$C_D, C_L$ という、車体周りの流体力学的な性能を評価するための重要な指標を得ることができます。
18. 時間発展
LBMの計算は1回だけではありません。本例題では、計算領域の大きさとLBM速度から、流体が計算領域を通過する回数を基準として計算ステップ数を決定します。設定では、「flow_passes = 2」となっているため、概念的には流れが計算領域を2回程度通過する時間スケールまで計算します。
時間発展では、$f_i^n \rightarrow f_i^{n+1}$ という更新を繰り返します。
概念的には、以下のような処理を繰り返します。
初期状態 → Collision → Streaming → Boundary → 次の時刻
→ Collision → Streaming → Boundary → ...
19. 計算途中で流れ場と流体力を出力する
本例題では、計算終了後だけでなく、計算途中でも結果を出力します。
一定の割合で、以下の処理を行います。
- マクロ量を計算
- 車体に作用する力を計算
- $C_D$ を計算
- $C_L$ を計算
- 速度場・密度場を出力
- 速度分布を画像として保存
20. 速度の大きさを可視化する
速度は、$\mathbf{u}=(u_x,u_y,u_z)$ として求められます。その大きさは、$|\mathbf{u}|=\sqrt{u_x^2+u_y^2+u_z^2}$ です。
本例題では、この速度場を利用して、計算領域内の流れの状態を可視化します。
特にAhmed Body後方の後流を見ることで、「流れの剥離・速度の低下・後流領域・車体周辺の流れの変化」などを確認できます。
21. 実験結果との比較
今回の例題では、計算結果を可視化するだけでなく、実験データとの比較も行います。ahmed.json に保存された実験データを読み込み、車体後方の複数の位置において速度分布を比較します。
つまり、「実験で観測された速度分布を、XLBのマルチ解像度LBMがどの程度再現できるか」 を確認できます。
22. この例題で最終的に得られるもの
このプログラムの出力を整理すると、大きく4種類あります。
① 流れ場：$u_x,\quad u_y,\quad u_z$ および速度の大きさ$|\mathbf{u}|$ です。
② 密度場：$\rho$ です。
③ 車体に作用する流体力：$F_D,\quad F_L$ から、$C_D,\quad C_L$ を求めます。
④ 実験値との比較：車体後方の速度分布を実験データと比較します。
全体を整理すると、以下のようになります。
multires_windtunnel_3d.py
┌────────────────┼────────────────┐
流れ場 流体力 実験比較
velocity Drag / Lift 速度分布
density Cd / Cl XLB vs 実験
└────────────────┼────────────────┘
↓
HDF5 / XDMF
23. この例題の本質
この multires_windtunnel_3d.py は、XLBを理解するうえで、flow_past_sphere_3d.py からさらに一歩進んだ例題です。
flow_past_sphere_3d.py では、「LBMによって3次元の物体周り流れを計算する」 ことが中心でした。
一方、この例題では、「複雑な3次元形状をSTLから読み込み、マルチ解像度LBMによって周囲の流れを計算し、さらに車体に作用する流体力を求め、実験結果と比較する」 ところまで進んでいます。
全体を整理すると、以下のような流れになります。
$\text{STL}\rightarrow\text{マルチ解像度格子} \rightarrow f_i\rightarrow\text{LBM時間発展}\rightarrow(\rho,\mathbf{u})\rightarrow\mathbf{F}\rightarrow(C_D,C_L)$
]
つまり、この例題を理解すると、「XLBはLBMを計算するだけのライブラリではなく、STLで定義された実際の3次元形状を対象として、マルチ解像度による効率的な流体解析を行い、さらに物体に作用する流体力や実験結果との比較まで行えるCFD基盤である」 ということが理解できます。
特に本例題の重要なポイントは、「細かい格子を計算領域全体に配置するのではなく、Ahmed Body周辺に細かい格子を集中させることで、計算精度と計算コストのバランスを取る」 というマルチ解像度LBMの考え方です。
解析形状ファイルahmed_25deg_m.stlの確認
スクリプトでは解析形状ファイルとして「Ahmed_25_NoLegs.stl」が指定されています。この代わりに「ahmed_25deg_m.stl」を用いてもよいかを確認します。
まず、標準的なAhmed bodyの25°モデルは、以下のような幾何学的な条件です。
| 項目 | 標準的な値 |
|---|---|
| 全長 | 1.044 m |
| 幅 | 0.389 m |
| 高さ | 0.288 m |
| 後部スラントφ | 25° |
| 脚 | 0.05 m程度 |
他の記述として、COMSOLの公式Ahmed body例でも、25°モデルについて全長1.044 m、幅0.389 m、高さ0.288 mが示されています。
以下の「CFD-Online」の解析モデルの解説でも、解析形状の寸法が指定されています。
https://www.cfd-online.com/Wiki/Ahmed_body
ここで「ahmed_25deg_m.stl」について、以下のPyhonスクリプトで、大きさを確認してみます。
python - <<'PY'
import trimesh
mesh = trimesh.load("ahmed_25deg_m.stl")
print("vertices :", len(mesh.vertices))
print("faces :", len(mesh.faces))
print("bounds :")
print(mesh.bounds)
size = mesh.bounds[1] - mesh.bounds[0]
print("size [m] :", size)
PY
vertices : 14315
faces : 28626
bounds :
[[-0.1 -0.1945 0. ]
[ 0.944 0.1945 0.338 ]]
size [m] : [1.044 0.389 0.338]
最後の数値を見ると、STLモデルのBounding Boxの寸法としては、車体の高さ0.288mに脚部分の0.05mを加えた0.338mとなり、全長と幅も対応しています。よって車体の寸法としては、利用できると思います。
ただし、利用するファイル「ahmed_25deg_m.stl」は、下図のように脚があり、スクリプトに指定されたファイル「Ahmed_25_NoLegs.stl」は、「NoLegs：脚なし」なので、その部分は異なっていることに、注意します。
例題スクリプトの具体的な解説
multires_windtunnel_3d.py の処理内容そのものは変更せず、学習用の説明コメントだけを追加して、①〜⑫の12ブロックに整理しました。
今回の12ブロックは次の構成です。
| ブロック | 内容 |
|---|---|
| ① | 解析条件・XLBの基本設定 |
| ② | STL形状からマルチ解像度計算領域を生成 |
| ③ | 入口・出口・壁・Ahmed Bodyの境界条件 |
| ④ | マルチ解像度LBM・KBCの初期化 |
| ⑤ | 抗力・揚力・係数の計算 |
| ⑥ | 実験データとの速度分布比較 |
| ⑦ | XLB初期化と計算格子生成 |
| ⑧ | 物理単位とLBM単位の変換 |
| ⑨ | 出力・計算格子・初期流れの準備 |
| ⑩ | 運動量交換・計算統計の準備 |
| ⑪ | LBM時間発展・計算結果の出力 |
| ⑫ | 計算終了後の後処理・実験値との比較 |
特に⑥では、先ほど確認した ahmed.json の19か所の測定位置をすべて対象にすることを明記しています。スクリプトでは data の各キーを順番に読み込み、各位置で (z) 方向の速度分布を取得して実験値と比較しています。:chatgpt-content-reference{index="0"}
また、②では実際のスクリプトにある5段階のマルチ解像度格子、③では RegularizedBC、DoNothingBC、HybridBC、⑪では sim.step() による時間発展、というように、実際のコードの構造に対応したコメントとしています。
"""
Ahmed Body aerodynamics with multi-resolution LBM
XLBを用いて、Ahmed Body（25度スラント）周りの
3次元流れをマルチ解像度LBMで解析する例題。
このスクリプトでは、
① 解析条件の設定
② STLからマルチ解像度計算領域を作成
③ 境界条件の設定
④ マルチ解像度LBMの初期化
⑤ 抗力・揚力の計算方法を定義
⑥ 実験データとの速度分布比較方法を定義
⑦ XLBと計算格子の初期化
⑧ 物理単位とLBM単位の変換
⑨ 出力・初期条件・計算ステップの準備
⑩ 運動量交換と計算統計の準備
⑪ LBMの時間発展
⑫ 計算終了後の後処理
という流れで計算を行う。
"""
# ============================================================
# 必要なライブラリを読み込む
# ============================================================
import neon
import warp as wp
import numpy as np
import time
import os
import matplotlib.pyplot as plt
import trimesh
import shutil
import xlb
from xlb.compute_backend import ComputeBackend
from xlb.precision_policy import PrecisionPolicy
from xlb.grid import multires_grid_factory
from xlb.operator.boundary_condition import (
DoNothingBC,
HybridBC,
RegularizedBC,
)
from xlb.operator.boundary_masker import MeshVoxelizationMethod
from xlb.utils.mesher import (
prepare_sparsity_pattern,
make_cuboid_mesh,
MultiresIO,
)
from xlb.utils import UnitConvertor
from xlb.operator.force import MultiresMomentumTransfer
from xlb.helper.initializers import CustomMultiresInitializer
# Warpのカーネルキャッシュをクリアする。
#
# XLBのNeonバックエンドではWarpも利用されるため、
# 実行開始時に以前のカーネルキャッシュの影響を避ける。
wp.clear_kernel_cache()
# Warpからの不要なメッセージを抑制する。
wp.config.quiet = True
# ============================================================
# ① 解析条件・XLBの基本設定
# ============================================================
#
# このブロックでは、解析全体の基本条件を設定する。
# 主な設定は、
# ・入口速度
# ・空気の動粘度
# ・最も細かい格子幅
# ・解析する流れの長さ
# ・Ahmed BodyのSTLファイル
# ・計算バックエンド
# ・数値精度
# ・LBMの速度モデル
# である。
# ------------------------------------------------------------
# 物理速度とLBM速度
# ------------------------------------------------------------
#
# 実際の風洞実験では、入口速度は
# U = 38.0 m/s
# としている。
# 一方、LBMでは物理単位そのものを直接計算するのではなく、
# 無次元化された格子単位（lattice unit）で計算する。
# ここでは
# wind_speed_mps = 38.0
# wind_speed_lbm = 0.05
# として、物理速度とLBM速度を対応付ける。
# 後でUnitConvertorを使用して、
# 物理単位 ⇔ LBM単位
# の変換を行う。
# ------------------------------------------------------------
# 格子幅
# ------------------------------------------------------------
#
# voxel_size = 0.005 m
# は最も細かい格子の格子幅である。
# 今回は5段階のマルチ解像度格子を使用するため、
# Level 1 : 0.005 m = 5 mm
# Level 2 : 0.010 m = 10 mm
# Level 3 : 0.020 m = 20 mm
# Level 4 : 0.040 m = 40 mm
# Level 5 : 0.080 m = 80 mm
# という関係になる。
# 車体付近では細かい格子を使い、
# 車体から遠い領域では粗い格子を使うことで、
# 計算量を削減する。
# ============================================================
# --------------------------
# Physical and simulation parameters
# --------------------------
# LBM内部で使用する入口速度
wind_speed_lbm = 0.05
# 実際の物理的な入口速度 [m/s]
wind_speed_mps = 38.0
# 計算領域を流れが何回通過する程度まで計算するか
flow_passes = 2
# 空気の動粘度 [m^2/s]
kinematic_viscosity = 1.508e-5
# 最も細かい格子幅 [m]
voxel_size = 0.005
# --------------------------
# STL geometry
# --------------------------
# Ahmed BodyのSTLファイル
stl_filename = "./ahmed_25deg_m.stl"
# 出力ディレクトリなどに使用する解析名
script_name = "Ahmed"
# --------------------------
# I/O settings
# --------------------------
# 計算の進捗を何%ごとに表示するか
print_interval_percentage = 1
# ファイル出力の頻度を途中で変更する位置
#
# ※元スクリプトの値とコメントをそのまま維持している。
file_output_crossover_percentage = 10
# 前半のファイル出力数
num_file_outputs_pre_crossover = 20
# 後半のファイル出力数
num_file_outputs_post_crossover = 5
# --------------------------
# XLB calculation settings
# --------------------------
# 計算バックエンドとしてNeonを使用する。
#
# 今回の例題では、マルチ解像度計算をNeonバックエンドで実行する。
compute_backend = ComputeBackend.NEON
# 計算精度
precision_policy = PrecisionPolicy.FP32FP32
# 3次元LBMの速度モデルとしてD3Q27を使用する。
#
# D3Q27は3次元空間で、
# 6方向の面方向
# 12方向の辺方向
# 8方向の頂点方向
# に加えて静止方向を持つ27速度モデルである。
#
# 今回はKBC衝突モデルを使用するため、
# D3Q27を使用している。
velocity_set = xlb.velocity_set.D3Q27(
precision_policy=precision_policy,
compute_backend=compute_backend,
)
# ============================================================
# ② STL形状からマルチ解像度計算領域を作る
# ============================================================
#
# この関数では、Ahmed BodyのSTL形状を読み込み、
# その形状を囲むマルチ解像度の計算領域を作成する。
# 処理の流れは、
# STL
# ↓
# 車体の大きさを取得
# ↓
# 計算領域内に配置
# ↓
# 5段階の解像度を設定
# ↓
# XLB用のlevel_dataを生成
# となる。
# ------------------------------------------------------------
# domain_multiplier
# ------------------------------------------------------------
#
# 各行が1つの解像度レベルを表している。
# 5行あるため、今回の計算は5段階の解像度を持つ。
# 最細格子幅は0.005 mなので、
# Δx_L = 0.005 × 2^(L-1)
# と考えることができる。
# したがって、
# 5 mm → 10 mm → 20 mm → 40 mm → 80 mm
# と格子が粗くなる。
# ------------------------------------------------------------
# 重要な考え方
# ------------------------------------------------------------
#
# Ahmed Body表面付近では流れの変化が大きい。
# 特に、
# ・車体前方の流れの分岐
# ・車体表面近傍の速度変化
# ・25度スラント部での流れの剥離
# ・車体後方の後流
# を表現するには細かい格子が必要である。
# 一方、車体から十分離れた場所では、
# 同じ細かさの格子を使用する必要性が小さい。
# そこで、
# 重要な場所 → 細かい格子
# 遠い場所 → 粗い格子
# とする。
# これがマルチ解像度LBMの基本的な考え方である。
# ============================================================
def generate_cuboid_mesh(stl_filename, voxel_size):
"""
STL形状を読み込み、
マルチ解像度の計算領域を生成する。
"""
# --------------------------------------------------------
# 各解像度レベルの計算領域を指定
# 各行は
# [-x, +x, -y, +y, -z, +z]
# 方向の領域倍率を表している。
# 車体周辺を細かくし、外側を粗くするため、
# レベルごとに異なる領域を設定している。
# --------------------------------------------------------
domain_multiplier = [
[3.0, 4.0, 2.5, 2.5, 0.0, 4.0],
[1.2, 1.25, 1.75, 1.75, 0.0, 1.5],
[0.8, 1.0, 1.25, 1.25, 0.0, 1.2],
[0.5, 0.65, 0.6, 0.60, 0.0, 0.6],
[0.25, 0.25, 0.25, 0.25, 0.0, 0.25],
]
# --------------------------------------------------------
# STLを読み込む
# --------------------------------------------------------
mesh = trimesh.load_mesh(
stl_filename,
process=False,
)
# STLが空の場合はエラー
if mesh.is_empty:
raise ValueError(
"Loaded mesh is empty or invalid."
)
# --------------------------------------------------------
# STLの最小座標・最大座標を取得
# --------------------------------------------------------
#
# min_bound : x,y,zそれぞれの最小座標
# max_bound : x,y,zそれぞれの最大座標
# ここからAhmed Bodyの大きさを取得する。
# --------------------------------------------------------
min_bound = mesh.vertices.min(axis=0)
max_bound = mesh.vertices.max(axis=0)
# 車体のx,y,z方向の大きさ
partSize = max_bound - min_bound
# Ahmed Body後端のx座標
#
# 実験データでは、この車体後端を
# x = 0 の基準として扱う。
x0 = max_bound[0]
# --------------------------------------------------------
# STLを計算領域内へ移動するための移動量を計算
# --------------------------------------------------------
stl_shift = np.array(
[
domain_multiplier[0][0] * partSize[0]
- min_bound[0],
domain_multiplier[0][2] * partSize[1]
- min_bound[1],
domain_multiplier[0][4] * partSize[2]
- min_bound[2],
],
dtype=float,
)
# --------------------------------------------------------
# STLを移動
# --------------------------------------------------------
mesh.apply_translation(stl_shift)
# 頂点法線を計算させる
_ = mesh.vertex_normals
# STL頂点をNumPy配列として取得
mesh_vertices = np.asarray(mesh.vertices)
# 一時STLとして保存
mesh.export("temp.stl")
# --------------------------------------------------------
# XLB用のマルチ解像度メッシュを生成
# --------------------------------------------------------
#
# make_cuboid_mesh()によって、
# ・各レベルの格子
# ・各レベルの領域
# ・STLとの位置関係
# などがlevel_dataとして作られる。
# --------------------------------------------------------
level_data = make_cuboid_mesh(
voxel_size,
domain_multiplier,
"temp.stl",
)
# 実際に生成されたレベル数
num_levels = len(level_data)
# 最細格子基準での全体格子形状
grid_shape_finest = tuple(
[
int(i * 2 ** (num_levels - 1))
for i in level_data[-1][0].shape
]
)
print(
f"Full shape based on finest voxel size is "
f"{grid_shape_finest}"
)
# 一時STLは不要なので削除
os.remove("temp.stl")
# 後で使用する情報を返す
return (
level_data,
mesh_vertices,
tuple(
[int(a) for a in grid_shape_finest]
),
stl_shift,
x0,
)
# ============================================================
# ③ 入口・出口・壁・Ahmed Bodyの境界条件
# ============================================================
#
# LBMでは計算領域の境界で、
# 「粒子分布をどのように扱うか」
# を指定する必要がある。
# 今回は、
# 左側 → 入口
# 右側 → 出口
# 上下面 → 流路境界
# 前後面 → 流路境界
# 車体 → 固体境界
# としている。
# ------------------------------------------------------------
# 入口
# ------------------------------------------------------------
#
# RegularizedBCを使用し、
# u = (U, 0, 0)
# という速度を与える。
# ------------------------------------------------------------
# 出口
# ------------------------------------------------------------
#
# DoNothingBCを使用する。
# 入口のように速度を固定するのではなく、
# 流体が計算領域外へ流出できる境界として扱う。
# ------------------------------------------------------------
# 車体
# ------------------------------------------------------------
#
# Ahmed BodyはSTLから読み込まれているため、
# 数式で球や円柱の内部を判定するのではない。
# STLの三角形表面をLBM格子に対応させ、
# Voxelizationによって固体境界を作る。
# ============================================================
def setup_boundary_conditions(
grid,
level_data,
body_vertices,
wind_speed_mps,
):
"""
各境界の境界条件を設定する。
"""
# --------------------------------------------------------
# 物理速度[m/s]をLBM速度へ変換
# --------------------------------------------------------
wind_speed_lbm = (
unit_convertor.velocity_to_lbm(
wind_speed_mps
)
)
# --------------------------------------------------------
# 各面に対応する格子インデックスを取得
# --------------------------------------------------------
left_indices = (
grid.boundary_indices_across_levels(
level_data,
box_side="left",
remove_edges=True,
)
)
right_indices = (
grid.boundary_indices_across_levels(
level_data,
box_side="right",
remove_edges=True,
)
)
top_indices = (
grid.boundary_indices_across_levels(
level_data,
box_side="top",
remove_edges=False,
)
)
bottom_indices = (
grid.boundary_indices_across_levels(
level_data,
box_side="bottom",
remove_edges=False,
)
)
front_indices = (
grid.boundary_indices_across_levels(
level_data,
box_side="front",
remove_edges=False,
)
)
back_indices = (
grid.boundary_indices_across_levels(
level_data,
box_side="back",
remove_edges=False,
)
)
# --------------------------------------------------------
# 入口
# --------------------------------------------------------
#
# 流れ方向をx方向として、
# u = (wind_speed_lbm, 0, 0)
# を指定する。
# つまり、
#
# x方向 → 38 m/s
# y方向 → 0
# z方向 → 0
# という流入条件である。
# --------------------------------------------------------
bc_inlet = RegularizedBC(
"velocity",
prescribed_value=(
wind_speed_lbm,
0.0,
0.0,
),
indices=left_indices,
)
# --------------------------------------------------------
# 出口
# --------------------------------------------------------
#
# DoNothingBCを使用する。
# --------------------------------------------------------
bc_outlet = DoNothingBC(
indices=right_indices
)
# --------------------------------------------------------
# 上面・下面・前面・背面
# --------------------------------------------------------
#
# Nonequilibrium Regularizedを利用した
# Hybrid Boundary Conditionを使用する。
# --------------------------------------------------------
bc_top = HybridBC(
bc_method="nonequilibrium_regularized",
indices=top_indices,
)
bc_bottom = HybridBC(
bc_method="nonequilibrium_regularized",
indices=bottom_indices,
)
bc_front = HybridBC(
bc_method="nonequilibrium_regularized",
indices=front_indices,
)
bc_back = HybridBC(
bc_method="nonequilibrium_regularized",
indices=back_indices,
)
# --------------------------------------------------------
# Ahmed Body
# --------------------------------------------------------
#
# STLの頂点をLBM単位へ変換する。
# さらに、
# MeshVoxelizationMethod("AABB_CLOSE")
# を利用してSTL形状を格子上へ対応付ける。
# use_mesh_distance=True
# によって、メッシュまでの距離情報も利用する。
# --------------------------------------------------------
bc_body = HybridBC(
bc_method="nonequilibrium_regularized",
mesh_vertices=unit_convertor.length_to_lbm(
body_vertices
),
voxelization_method=MeshVoxelizationMethod(
"AABB_CLOSE",
close_voxels=4,
),
use_mesh_distance=True,
)
# 境界条件をリストとして返す
return [
bc_top,
bc_bottom,
bc_front,
bc_back,
bc_inlet,
bc_outlet,
bc_body,
]
# ============================================================
# ④ マルチ解像度LBMの計算モデルを初期化する
# ============================================================
#
# ここでは実際のLBM計算を管理する
# MultiresSimulationManagerを作成する。
# collision_type="KBC"
# により、KBC衝突モデルを使用する。
# LBMでは、
# ① 分布関数の衝突
# ② 分布関数の移流
# ③ 境界条件
# を繰り返して流れ場を時間発展させる。
# この部分では、その時間発展全体を管理する
# シミュレーションマネージャを構築する。
# ============================================================
def initialize_simulation(
grid,
boundary_conditions,
omega_finest,
initializer,
collision_type="KBC",
mres_perf_opt=xlb.MresPerfOptimizationType.FUSION_AT_FINEST,
):
"""
マルチ解像度LBM計算を初期化する。
"""
sim = xlb.helper.MultiresSimulationManager(
omega_finest=omega_finest,
grid=grid,
boundary_conditions=boundary_conditions,
collision_type=collision_type,
initializer=initializer,
mres_perf_opt=mres_perf_opt,
)
return sim
# ============================================================
# ⑤ 車体に作用する抗力・揚力を計算する
# ============================================================
#
# LBMでは、車体表面で流体と固体の間に運動量交換が発生する。
# MultiresMomentumTransferを利用すると、
# その運動量交換から車体に作用する力を計算できる。
# 今回は、
# F = (Fx, Fy, Fz)
# としたとき、
# Drag = Fx
# Lift = Fz
# としている。
# さらに、
# Cd = 2 FD / (U^2 A)
# Cl = 2 FL / (U^2 A)
# に対応する形で、LBM単位の抗力・揚力から
# 無次元の抗力係数Cd、揚力係数Clを計算する。
# ============================================================
def compute_force_coefficients(
sim,
step,
momentum_transfer,
wind_speed_lbm,
reference_area,
):
"""
抗力・揚力と、その係数を計算する。
"""
# 車体表面での運動量交換から力を計算
boundary_force = momentum_transfer(
sim.f_0,
sim.f_1,
sim.bc_mask,
sim.missing_mask,
)
# x方向の力 → 抗力
drag = boundary_force[0]
# z方向の力 → 揚力
lift = boundary_force[2]
# 抗力係数
cd = (
2.0 * drag
/ (
wind_speed_lbm**2
* reference_area
)
)
# 揚力係数
cl = (
2.0 * lift
/ (
wind_speed_lbm**2
* reference_area
)
)
# 数値計算が異常になった場合は停止
if np.isnan(cd) or np.isnan(cl):
print(
f"NaN detected in coefficients "
f"at step {step}"
)
raise ValueError(
f"NaN detected in coefficients "
f"at step {step}: "
f"Cd={cd}, Cl={cl}"
)
# 後で時間履歴を保存するために記録
drag_values.append([cd, cl])
return cd, cl, drag
# ------------------------------------------------------------
# Cd・Clの時間変化をグラフ化する関数
# ------------------------------------------------------------
def plot_force_coefficients(
drag_values,
output_dir,
print_interval,
script_name,
percentile_range=(15, 85),
use_log_scale=False,
):
"""
CdとClの時間変化をグラフとして保存する。
"""
drag_values_array = np.array(
drag_values
)
# 横軸：計算ステップ
steps = np.arange(
0,
len(drag_values) * print_interval,
print_interval,
)
# Cd
cd_values = drag_values_array[:, 0]
# Cl
cl_values = drag_values_array[:, 1]
# グラフの縦軸範囲を計算
y_min = min(
np.percentile(
cd_values,
percentile_range[0],
),
np.percentile(
cl_values,
percentile_range[0],
),
)
y_max = max(
np.percentile(
cd_values,
percentile_range[1],
),
np.percentile(
cl_values,
percentile_range[1],
),
)
padding = (
y_max - y_min
) * 0.1
y_min = y_min - padding
y_max = y_max + padding
if use_log_scale:
y_min = max(
y_min,
1e-6,
)
# グラフを作成
plt.figure(
figsize=(10, 6)
)
plt.plot(
steps,
cd_values,
label="Drag Coefficient (Cd)",
color="blue",
)
plt.plot(
steps,
cl_values,
label="Lift Coefficient (Cl)",
color="red",
)
plt.xlabel(
"Simulation Step"
)
plt.ylabel(
"Coefficient"
)
plt.title(
f"{script_name}: "
"Drag and Lift Coefficients Over Time"
)
plt.legend()
plt.grid(True)
plt.ylim(
y_min,
y_max,
)
if use_log_scale:
plt.yscale("log")
plt.savefig(
os.path.join(
output_dir,
"drag_lift_plot.png",
)
)
plt.close()
# ============================================================
# ⑥ 実験データとXLB計算結果を比較する
# ============================================================
#
# Ahmed Bodyでは、車体後方の流れが重要である。
# この関数では、
# ahmed.json
# に保存された実験データを読み込み、
# XLBで計算した速度分布と比較する。
# ------------------------------------------------------------
# 実験データ
# ------------------------------------------------------------
#
# ahmed.jsonには複数のx位置が登録されている。
# 今回のデータでは、
# -1.162
# -1.062
# ...
# -0.012
# 0.038
# 0.088
# 0.138
# 0.188
# 0.238
# 0.288
# などの位置が存在する。
# 車体後端をx=0とすると、
# x > 0
# が車体後方である。
# 特に、
# 0.038 ～ 0.288 m
# の6点が車体後方の速度分布測定位置である。
# ------------------------------------------------------------
# 比較する速度
# ------------------------------------------------------------
#
# y=0の対称面上で、
# z = 0 ～ 0.8 m
# の鉛直方向の速度分布を取得する。
# 取得するのは速度ベクトルのx成分、
# Ux
# である。
# そのため最終的なグラフは、
# 横軸：Ux [m/s]
# 縦軸：z [m]
# となる。
# 実験値とXLB計算値を同じグラフに描くことで、
# 後流の速度分布をどの程度再現できているか確認できる。
# ============================================================
def compute_voxel_statistics(
sim,
bc_mask_exporter,
sparsity_pattern,
boundary_conditions,
unit_convertor,
):
"""
各解像度レベルのセル数などを集計する。
"""
fields_data = (
bc_mask_exporter.get_fields_data(
{"bc_mask": sim.bc_mask}
)
)
bc_mask_data = (
fields_data["bc_mask_0"]
)
level_id_field = (
bc_mask_exporter.level_id_field
)
# --------------------------------------------------------
# 各レベルの固体セル数を計算
# --------------------------------------------------------
solid_voxels = []
for lvl in range(num_levels):
level_mask = (
level_id_field == lvl
)
solid_voxels.append(
np.sum(
bc_mask_data[level_mask]
== 255
)
)
# --------------------------------------------------------
# 有効な流体セル数を計算
# --------------------------------------------------------
active_voxels = [
np.count_nonzero(mask)
for mask in sparsity_pattern
]
active_voxels = [
max(
0,
active_voxels[lvl]
- solid_voxels[lvl],
)
for lvl in range(num_levels)
]
# 全レベルの流体セル数
total_voxels = sum(
active_voxels
)
# --------------------------------------------------------
# 1グローバルステップあたりの格子更新数
# --------------------------------------------------------
#
# 細かい格子ほど1セルあたりの計算負荷が高くなるため、
# レベルごとの更新量を考慮して合計する。
# --------------------------------------------------------
total_lattice_updates_per_step = sum(
active_voxels[lvl]
* (
2 ** (
num_levels - 1 - lvl
)
)
for lvl in range(num_levels)
)
# --------------------------------------------------------
# 抗力・揚力係数に使う基準面積を計算
# --------------------------------------------------------
finest_level = 0
mask_finest = (
level_id_field
== finest_level
)
bc_mask_finest = (
bc_mask_data[mask_finest]
)
active_indices_finest = np.argwhere(
sparsity_pattern[0]
)
# 最後に登録された境界条件が車体
bc_body_id = (
boundary_conditions[-1].id
)
solid_voxels_indices = (
active_indices_finest[
bc_mask_finest
== bc_body_id
]
)
# YZ平面へ投影した面積を取得
unique_jk = np.unique(
solid_voxels_indices[:, 1:3],
axis=0,
)
reference_area = (
unique_jk.shape[0]
)
# LBM単位の面積を物理単位へ変換
reference_area_physical = (
reference_area
* unit_convertor.reference_length**2
)
return {
"active_voxels": active_voxels,
"solid_voxels": solid_voxels,
"total_voxels": total_voxels,
"total_lattice_updates_per_step":
total_lattice_updates_per_step,
"reference_area":
reference_area,
"reference_area_physical":
reference_area_physical,
}
# ------------------------------------------------------------
# シミュレーション結果と実験結果の比較
# ------------------------------------------------------------
def plot_data(
x0,
output_dir,
delta_x_coarse,
sim,
IOexporter,
prefix="Ahmed",
):
"""
Ahmed Bodyの各x位置で、
実験値とXLB計算値の速度分布を比較する。
"""
# --------------------------------------------------------
# XLBから出力されたCSVを読み込む関数
# --------------------------------------------------------
def _load_sim_line(csv_path):
"""
IOexporter.to_line()で出力したCSVを読み込む。
戻り値：
z : 高さ
Ux : x方向速度
"""
data = np.genfromtxt(
csv_path,
delimiter=",",
names=True,
autostrip=True,
dtype=None,
encoding="utf-8",
)
if data.size == 0:
raise ValueError(
f"No data in {csv_path}"
)
z = np.asarray(
data["z"],
dtype=float,
)
ux = np.asarray(
data["value"],
dtype=float,
)
return z, ux
# --------------------------------------------------------
# 実験データを読み込む
# --------------------------------------------------------
import json
ref_data_path = "./ahmed.json"
with open(
ref_data_path,
"r",
) as file:
data = json.load(file)
# --------------------------------------------------------
# ahmed.jsonに登録された全x位置について繰り返す
# --------------------------------------------------------
for x_str in data["data"].keys():
# 実験値のx方向速度
refX = np.array(
data["data"][x_str]["x-velocity"]
)
# 実験値の高さ
refY = np.array(
data["data"][x_str]["height"]
)
# ----------------------------------------------------
# 車体後端x0を基準にして、
# 実験データの相対位置x_posを実座標へ変換
# ----------------------------------------------------
x_pos = float(x_str)
x1 = x0 + x_pos
print(
f" x1 is {x1}"
)
# ----------------------------------------------------
# LBMの分布関数から速度場を再計算
# ----------------------------------------------------
sim.macro(
sim.f_0,
sim.bc_mask,
sim.rho,
sim.u,
streamId=0,
)
filename = os.path.join(
output_dir,
f"{prefix}_{x_str}",
)
wp.synchronize()
# ----------------------------------------------------
# y=0の対称面上で、
# x=x1の鉛直線に沿って速度を取得
# ----------------------------------------------------
#
# start_point=(x1,0,0)
# end_point=(x1,0,0.8)
#
# なので、xとyを固定し、
# z方向に沿って速度をサンプリングする。
# ----------------------------------------------------
IOexporter.to_line(
filename,
{"velocity": sim.u},
start_point=(
x1,
0,
0,
),
end_point=(
x1,
0,
0.8,
),
resolution=250,
# velocityのx成分を取得
component=0,
# マルチ解像度格子でのモデル単位を考慮
radius=delta_x_coarse,
)
# ----------------------------------------------------
# CSVを読み込む
# ----------------------------------------------------
csv_path = (
filename
+ "_velocity_0.csv"
)
print(
f"CSV path is {csv_path}"
)
try:
sim_z, sim_ux = (
_load_sim_line(
csv_path
)
)
except Exception as e:
print(
f"Failed to read "
f"{csv_path}: {e}"
)
continue
# ----------------------------------------------------
# 実験値とXLB計算値を同じグラフに描画
# ----------------------------------------------------
plt.figure(
figsize=(4.5, 6)
)
# 実験値
plt.plot(
refX,
refY,
"o",
mfc="none",
label="Experimental)",
)
# XLB計算値
plt.plot(
sim_ux,
sim_z,
"-",
lw=2,
label="Simulation",
)
plt.xlim(
np.min(refX) * 0.9,
np.max(refX) * 1.1,
)
plt.ylim(
np.min(refY),
np.max(refY),
)
plt.xlabel(
"Ux [m/s]"
)
plt.ylabel(
"z [m]"
)
plt.title(
f"Velocity Plot at "
f"{x_pos:+.3f}"
)
plt.grid(
True,
alpha=0.3,
)
plt.legend()
plt.tight_layout()
# 各位置ごとにPNGを保存
plt.savefig(
filename + ".png",
dpi=150,
)
plt.close()
# ============================================================
# ⑦ XLBを初期化し、マルチ解像度格子を生成する
# ============================================================
#
# ここからmain処理に入る。
# まずxlb.init()で、
# ・D3Q27
# ・Neon
# ・FP32FP32
# という計算環境をXLBへ登録する。
# 次にgenerate_cuboid_mesh()を呼び出し、
# STLからマルチ解像度格子情報を生成する。
# その後prepare_sparsity_pattern()で、
# 実際に計算する必要のある格子セルを整理する。
# ============================================================
# ------------------------------------------------------------
# XLB initialization
# ------------------------------------------------------------
xlb.init(
velocity_set=velocity_set,
default_backend=compute_backend,
default_precision_policy=precision_policy,
)
# ------------------------------------------------------------
# STL → マルチ解像度格子
# ------------------------------------------------------------
(
level_data,
body_vertices,
grid_shape_zip,
stl_shift,
x0,
) = generate_cuboid_mesh(
stl_filename,
voxel_size,
)
# ------------------------------------------------------------
# 各レベルの有効セル情報を準備
# ------------------------------------------------------------
sparsity_pattern, level_origins = (
prepare_sparsity_pattern(
level_data
)
)
# ============================================================
# ⑧ 物理単位とLBM単位を対応付ける
# ============================================================
#
# LBM計算では、
# 長さ
# 速度
# 時間
# 粘性
# などをLBM単位で扱う。
# 一方、今回の入力条件は、
# 速度 = 38.0 m/s
# 格子幅 = 0.005 m
# 動粘度 = 1.508e-5 m^2/s
# という物理単位で与えている。
# UnitConvertorによって両者を対応付ける。
# ------------------------------------------------------------
# マルチ解像度
# ------------------------------------------------------------
#
# num_levels = 5 の場合、
# 最細 = 0.005 m
# 最粗 = 0.005 × 2^4
# = 0.080 m
# となる。
# ============================================================
unit_convertor = UnitConvertor(
velocity_lbm_unit=wind_speed_lbm,
velocity_physical_unit=wind_speed_mps,
voxel_size_physical_unit=voxel_size,
)
# レベル数
num_levels = len(
level_data
)
# 最粗格子の格子幅
delta_x_coarse = (
voxel_size
* 2 ** (num_levels - 1)
)
# 動粘度をLBM単位へ変換
nu_lattice = (
unit_convertor.viscosity_to_lbm(
kinematic_viscosity
)
)
# ------------------------------------------------------------
# 緩和係数omegaを計算
# ------------------------------------------------------------
#
# BGK型LBMなどでは緩和時間と粘性が関係する。
# 今回のMultiresSimulationManagerでは、
# 最細レベルのomegaを計算して渡している。
# omega = 1 / (3ν + 0.5)
# ------------------------------------------------------------
omega_finest = (
1.0
/ (
3.0 * nu_lattice
+ 0.5
)
)
# ============================================================
# ⑨ 出力管理・計算格子・初期条件を準備する
# ============================================================
#
# このブロックでは、実際の計算を開始する前の準備を行う。
# ・出力ディレクトリ
# ・HDF5/XDMF出力
# ・マルチ解像度格子
# ・計算ステップ数
# ・出力間隔
# ・初期速度
# ・LBMシミュレーションオブジェクト
# を順番に準備する。
# ============================================================
# ------------------------------------------------------------
# 出力ディレクトリを作成
# ------------------------------------------------------------
current_dir = os.path.join(
os.path.dirname(__file__)
)
output_dir = os.path.join(
current_dir,
script_name,
)
# 既存フォルダがあれば削除
if os.path.exists(
output_dir
):
shutil.rmtree(
output_dir
)
# 新しく作成
os.makedirs(
output_dir
)
# ------------------------------------------------------------
# 速度・密度の出力設定
# ------------------------------------------------------------
field_name_cardinality_dict = {
"velocity": 3,
"density": 1,
}
# HDF5/XDMF出力用
h5exporter = MultiresIO(
field_name_cardinality_dict,
level_data,
offset=-stl_shift,
unit_convertor=unit_convertor,
)
# 境界マスク出力用
bc_mask_exporter = MultiresIO(
{"bc_mask": 1},
level_data,
offset=-stl_shift,
unit_convertor=unit_convertor,
)
# ------------------------------------------------------------
# マルチ解像度計算格子を生成
# ------------------------------------------------------------
#
# ここで初めて、XLBの計算に使用するgridオブジェクトを作る。
# ============================================================
grid = multires_grid_factory(
grid_shape_zip,
velocity_set=velocity_set,
sparsity_pattern_list=sparsity_pattern,
sparsity_pattern_origins=[
neon.Index_3d(*box_origin)
for box_origin in level_origins
],
)
# ------------------------------------------------------------
# 最粗レベルの格子サイズを取得
# ------------------------------------------------------------
coarsest_level = (
grid.count_levels - 1
)
grid_shape_x_coarsest = (
grid.level_to_shape(
coarsest_level
)[0]
)
# ------------------------------------------------------------
# 計算ステップ数
# ------------------------------------------------------------
#
# flow_passesと入口速度から、
# 流体が計算領域を何回通過する程度まで
# 時間発展させるかを決める。
# ------------------------------------------------------------
num_steps = int(
flow_passes
* (
grid_shape_x_coarsest
/ wind_speed_lbm
)
)
# ------------------------------------------------------------
# 進捗表示の間隔
# ------------------------------------------------------------
print_interval = max(
1,
int(
num_steps
* (
print_interval_percentage
/ 100.0
)
),
)
# ------------------------------------------------------------
# ファイル出力頻度を設定
# ------------------------------------------------------------
crossover_step = int(
num_steps
* (
file_output_crossover_percentage
/ 100.0
)
)
file_output_interval_pre_crossover = (
max(
1,
int(
crossover_step
/ num_file_outputs_pre_crossover
),
)
if num_file_outputs_pre_crossover > 0
else num_steps + 1
)
file_output_interval_post_crossover = (
max(
1,
int(
(
num_steps
- crossover_step
)
/ num_file_outputs_post_crossover
),
)
if num_file_outputs_post_crossover > 0
else num_steps + 1
)
# ------------------------------------------------------------
# 境界条件を生成
# ------------------------------------------------------------
boundary_conditions = (
setup_boundary_conditions(
grid,
level_data,
body_vertices,
wind_speed_mps,
)
)
# ------------------------------------------------------------
# 初期速度場を設定
# ------------------------------------------------------------
#
# 最初の計算状態では、
# 計算領域全体に入口速度に対応する速度を与える。
# bc_idには出口の境界IDを渡している。
# ============================================================
wind_speed_lbm = (
unit_convertor.velocity_to_lbm(
wind_speed_mps
)
)
initializer = CustomMultiresInitializer(
bc_id=boundary_conditions[-2].id,
constant_velocity_vector=(
wind_speed_lbm,
0.0,
0.0,
),
velocity_set=velocity_set,
precision_policy=precision_policy,
compute_backend=compute_backend,
)
# ------------------------------------------------------------
# マルチ解像度LBMを初期化
# ------------------------------------------------------------
sim = initialize_simulation(
grid,
boundary_conditions,
omega_finest,
initializer,
)
# ------------------------------------------------------------
# 格子数・固体セル数・基準面積などを計算
# ------------------------------------------------------------
stats = compute_voxel_statistics(
sim,
bc_mask_exporter,
sparsity_pattern,
boundary_conditions,
unit_convertor,
)
active_voxels = (
stats["active_voxels"]
)
solid_voxels = (
stats["solid_voxels"]
)
total_voxels = (
stats["total_voxels"]
)
total_lattice_updates_per_step = (
stats[
"total_lattice_updates_per_step"
]
)
reference_area = (
stats["reference_area"]
)
reference_area_physical = (
stats[
"reference_area_physical"
]
)
# ------------------------------------------------------------
# 初期状態の境界マスクを保存
# ------------------------------------------------------------
filename = os.path.join(
output_dir,
f"{script_name}_initial_bc_mask",
)
try:
bc_mask_exporter.to_hdf5(
filename,
{
"bc_mask":
sim.bc_mask
},
compression="gzip",
compression_opts=0,
)
xmf_filename = (
f"{filename}.xmf"
)
hdf5_basename = (
f"{script_name}"
"_initial_bc_mask.h5"
)
except Exception as e:
print(
"Error during initial "
f"bc_mask output: {e}"
)
wp.synchronize()
# ============================================================
# ⑩ 車体表面の運動量交換と計算統計を準備する
# ============================================================
#
# MultiresMomentumTransferを作成する。
# これによって、計算中にAhmed Body表面へ作用する力を
# 運動量交換から計算できる。
# また、これまでに計算した格子数などを画面に表示する。
# ここで、
# ・最細格子幅
# ・最粗格子幅
# ・レベル数
# ・流体セル数
# ・固体セル数
# ・物理入口速度
# ・LBM入口速度
# ・基準面積
# を確認できる。
# ============================================================
# ------------------------------------------------------------
# 運動量交換法による力の計算器
# ------------------------------------------------------------
momentum_transfer = MultiresMomentumTransfer(
boundary_conditions[-1],
mres_perf_opt=
xlb.MresPerfOptimizationType.
FUSION_AT_FINEST,
compute_backend=compute_backend,
)
# ------------------------------------------------------------
# 計算条件を表示
# ------------------------------------------------------------
print(
"\n"
+ "=" * 50
+ "\n"
)
print(
f"Number of flow passes: "
f"{flow_passes}"
)
print(
f"Calculated iterations: "
f"{num_steps:,}"
)
print(
f"Finest voxel size: "
f"{voxel_size} meters"
)
print(
f"Coarsest voxel size: "
f"{delta_x_coarse} meters"
)
print(
f"Total voxels: "
f"{sum(np.count_nonzero(mask) for mask in sparsity_pattern):,}"
)
print(
f"Total active voxels: "
f"{total_voxels:,}"
)
print(
f"Active voxels per level: "
f"{[int(v) for v in active_voxels]}"
)
print(
f"Solid voxels per level: "
f"{[int(v) for v in solid_voxels]}"
)
print(
f"Total lattice updates per global step: "
f"{total_lattice_updates_per_step:,}"
)
print(
f"Number of refinement levels: "
f"{num_levels}"
)
print(
f"Physical inlet velocity: "
f"{wind_speed_mps:.4f} m/s"
)
print(
f"Lattice velocity (ulb): "
f"{wind_speed_lbm}"
)
print(
f"Computed reference area (bc_mask): "
f"{reference_area} lattice units"
)
print(
f"Physical reference area (bc_mask): "
f"{reference_area_physical:.6f} m^2"
)
print(
"\n"
+ "=" * 50
+ "\n"
)
# ============================================================
# ⑪ LBMの時間発展を実行する
# ============================================================
#
# ここがシミュレーションの中心部分である。
# forループの1回がLBMの1ステップに相当する。
# sim.step()
# によって、
# ・分布関数の衝突
# ・分布関数の移流
# ・境界条件
# ・マルチ解像度格子間の処理
# が進められる。
# これをnum_steps回繰り返すことで、
# 初期状態からAhmed Body周りの流れ場が発達していく。
# ------------------------------------------------------------
# 定期的に行う処理
# ------------------------------------------------------------
#
# 一定ステップごとに、
# ① 分布関数から速度・密度を計算
# ② Cd・Clを計算
# ③ HDF5/XDMFへ保存
# ④ 速度場のスライス画像を保存
# ⑤ 計算性能を表示
# する。
# ============================================================
wp.synchronize()
start_time = time.time()
compute_time = 0.0
steps_since_last_print = 0
# CdとClの時間履歴
drag_values = []
# ------------------------------------------------------------
# メインの時間発展ループ
# ------------------------------------------------------------
for step in range(num_steps):
step_start = time.time()
# --------------------------------------------------------
# LBMを1ステップ進める
# --------------------------------------------------------
sim.step()
wp.synchronize()
# 計算時間を記録
compute_time += (
time.time()
- step_start
)
steps_since_last_print += 1
# ========================================================
# 進捗表示のタイミング
# ========================================================
if (
step % print_interval == 0
or step == num_steps - 1
):
# ----------------------------------------------------
# 分布関数から密度・速度を計算
# ----------------------------------------------------
sim.macro(
sim.f_0,
sim.bc_mask,
sim.rho,
sim.u,
streamId=0,
)
wp.synchronize()
# ----------------------------------------------------
# 抗力・揚力を計算
# ----------------------------------------------------
cd, cl, drag = (
compute_force_coefficients(
sim,
step,
momentum_transfer,
wind_speed_lbm,
reference_area,
)
)
# ----------------------------------------------------
# 出力ファイル名
# ----------------------------------------------------
filename = os.path.join(
output_dir,
f"{script_name}_{step:04d}",
)
# ----------------------------------------------------
# 速度・密度をHDF5へ保存
# ----------------------------------------------------
h5exporter.to_hdf5(
filename,
{
"velocity": sim.u,
"density": sim.rho,
},
compression="gzip",
compression_opts=0,
)
# ----------------------------------------------------
# 速度場の断面画像を出力
# ----------------------------------------------------
#
# plane_normal=(0,1,0)
#
# なので、y方向に垂直な断面、
# すなわちx-z平面を取り出す。
#
# Ahmed Bodyの中心対称面付近の流れを
# 画像として確認することができる。
# ----------------------------------------------------
h5exporter.to_slice_image(
filename,
{
"velocity": sim.u
},
plane_point=(
1,
0,
0,
),
plane_normal=(
0,
1,
0,
),
grid_res=2000,
bounds=(
0.25,
0.75,
0,
0.5,
),
show_axes=False,
show_colorbar=False,
slice_thickness=
delta_x_coarse,
)
# ----------------------------------------------------
# 計算時間と性能を計算
# ----------------------------------------------------
end_time = time.time()
elapsed = (
end_time
- start_time
)
total_lattice_updates = (
total_lattice_updates_per_step
* steps_since_last_print
)
# MLUPS
#
# Million Lattice Updates Per Second
#
# LBMの計算性能を表す代表的な指標。
MLUPS = (
total_lattice_updates
/ compute_time
/ 1e6
if compute_time > 0
else 0.0
)
# 現在までに何回分の流れが
# 計算領域を通過したか
current_flow_passes = (
step
* wind_speed_lbm
/ grid_shape_x_coarsest
)
# 残りステップ
remaining_steps = (
num_steps
- step
- 1
)
# 残り時間の推定
time_remaining = (
0.0
if MLUPS == 0
else (
total_lattice_updates_per_step
* remaining_steps
)
/ (
MLUPS
* 1e6
)
)
hours, rem = divmod(
time_remaining,
3600,
)
minutes, seconds = divmod(
rem,
60,
)
time_remaining_str = (
f"{int(hours):02d}h "
f"{int(minutes):02d}m "
f"{int(seconds):02d}s"
)
# 進捗率
percent_complete = (
(step + 1)
/ num_steps
* 100
)
# ----------------------------------------------------
# 進捗を画面に表示
# ----------------------------------------------------
print(
f"Completed step "
f"{step}/{num_steps} "
f"({percent_complete:.2f}% complete)"
)
print(
f" Flow Passes: "
f"{current_flow_passes:.2f}"
)
print(
f" Time elapsed: "
f"{elapsed:.1f}s, "
f"Compute time: "
f"{compute_time:.1f}s, "
f"ETA: "
f"{time_remaining_str}"
)
print(
f" MLUPS: "
f"{MLUPS:.1f}"
)
print(
f" Cd={cd:.3f}, "
f"Cl={cl:.3f}, "
f"Drag Force "
f"(lattice units)="
f"{drag:.3f}"
)
# 次の表示区間のためにリセット
start_time = time.time()
compute_time = 0.0
steps_since_last_print = 0
# ========================================================
# ファイル出力
# ========================================================
# 計算前半と後半で出力頻度を変える
file_output_interval = (
file_output_interval_pre_crossover
if step < crossover_step
else file_output_interval_post_crossover
)
if (
step % file_output_interval == 0
or step == num_steps - 1
):
# 最新の速度・密度を取得
sim.macro(
sim.f_0,
sim.bc_mask,
sim.rho,
sim.u,
streamId=0,
)
filename = os.path.join(
output_dir,
f"{script_name}_{step:04d}",
)
try:
# HDF5として保存
h5exporter.to_hdf5(
filename,
{
"velocity": sim.u,
"density": sim.rho,
},
compression="gzip",
compression_opts=0,
)
# XDMFファイル名
xmf_filename = (
f"{filename}.xmf"
)
# 対応するHDF5ファイル名
hdf5_basename = (
f"{script_name}_{step:04d}.h5"
)
except Exception as e:
print(
f"Error during file "
f"output at step "
f"{step}: {e}"
)
wp.synchronize()
# ========================================================
# 最終ステップ
# ========================================================
if step == num_steps - 1:
# 最終流れ場について、
# 実験データとの速度分布比較を行う。
plot_data(
x0,
output_dir,
delta_x_coarse,
sim,
h5exporter,
prefix="Ahmed",
)
# ============================================================
# ⑫ 計算終了後の後処理・結果保存
# ============================================================
#
# 計算終了後は、
# ・Cdの時間履歴
# ・Clの時間履歴
# ・Cd/Clのグラフ
# ・平均Cd
# ・平均Cl
# ・実験Cdとの比較
# を行う。
# つまり、この例題は単に「流れを計算して画像を出す」
# だけではなく、
# 流れ場
# ↓
# 車体表面の力
# ↓
# Cd・Cl
# ↓
# 実験値との比較
# という工学的な解析まで行っている。
# ============================================================
# ------------------------------------------------------------
# Cd・ClをCSVへ保存
# ------------------------------------------------------------
if len(drag_values) > 0:
with open(
os.path.join(
output_dir,
"drag_lift.csv",
),
"w",
) as fd:
fd.write(
"Step,Cd,Cl\n"
)
for i, (cd, cl) in enumerate(
drag_values
):
fd.write(
f"{i * print_interval},"
f"{cd},{cl}\n"
)
# Cd・Clの時間変化をグラフ化
plot_force_coefficients(
drag_values,
output_dir,
print_interval,
script_name,
)
# ------------------------------------------------------------
# 計算後半50%の平均Cd・Clを計算
# ------------------------------------------------------------
drag_values_array = np.array(
drag_values
)
if len(drag_values) > 0:
# 後半50%の開始位置
start_index = (
len(drag_values)
// 2
)
# 後半50%だけを取り出す
last_half = (
drag_values_array[
start_index:,
:
]
)
# 平均抗力係数
avg_cd = np.mean(
last_half[:, 0]
)
# 平均揚力係数
avg_cl = np.mean(
last_half[:, 1]
)
print(
"Average Drag Coefficient "
"(Cd) for last 50%: "
f"{avg_cd:.6f}"
)
print(
"Average Lift Coefficient "
"(Cl) for last 50%: "
f"{avg_cl:.6f}"
)
# --------------------------------------------------------
# 実験値との比較
# --------------------------------------------------------
#
# この例題では実験Cdとして0.3088を使用している。
# これはahmed.jsonの速度分布そのものではなく、
# スクリプト末尾で比較対象として設定されている値である。
# --------------------------------------------------------
print(
"Experimental Drag "
"Coefficient (Cd): "
f"{0.3088}"
)
# 相対誤差
print(
"Error Drag Coefficient "
"(Cd): "
f"{((avg_cd - 0.3088) / 0.3088) * 100:.2f}%"
)
else:
print(
"No drag or lift data collected."
)
例題スクリプトの結果の分析
以下の記事【05:Neon導入方法とマルチ解像度例題の分析】で解説した手順を用いて、環境を準備して例題を実行すると、多数の解析結果が作られます。画面に表示されるログについては、【05】で概要を確認しました。
なお、正常に解析が進む場合は、ログファイル「Neon.log」には、何も書かれないようです。画面表示のログも保存されていないようです。
ここでは、ディレクトリ「Ahmed」に保存されたファイルについて分析します。ファイルが多数あるので、種類ごとにディレクトリを作り、ここに移動してから、ファイルを調べます。
(XLBneon) ryos@dalab2403sv01:~/XLBneonWork/MWIND3D/Ahmed$ mkdir D-h5+xmf
(XLBneon) ryos@dalab2403sv01:~/XLBneonWork/MWIND3D/Ahmed$ mv *.h5 D-h5+xmf/
(XLBneon) ryos@dalab2403sv01:~/XLBneonWork/MWIND3D/Ahmed$ mv *.xmf D-h5+xmf/
(XLBneon) ryos@dalab2403sv01:~/XLBneonWork/MWIND3D/Ahmed$ mkdir D-volmagpng
(XLBneon) ryos@dalab2403sv01:~/XLBneonWork/MWIND3D/Ahmed$ mv *magnitude.png D-volmagpng/
(XLBneon) ryos@dalab2403sv01:~/XLBneonWork/MWIND3D/Ahmed$ mkdir D-volplot
(XLBneon) ryos@dalab2403sv01:~/XLBneonWork/MWIND3D/Ahmed$ mv *.csv D-volplot/
(XLBneon) ryos@dalab2403sv01:~/XLBneonWork/MWIND3D/Ahmed$ mv *.png D-volplot/
(XLBneon) ryos@dalab2403sv01:~/XLBneonWork/MWIND3D/Ahmed$ ls
D-h5+xmf D-volmagpng D-volplot
拡張子 h5, xmfファイルの分析
XLBの Ahmedボディ多重解像度計算 で、D=h5+xmf ディレクトリに保存される拡張子 .h5 と .xmf は、解析結果をParaViewなどで可視化するための一組のファイルと考えると分かりやすいです。
XLB公式ドキュメントでも、Neonによるマルチレゾリューション計算では HDF5/XDMF形式でマルチ解像度データを出力できると説明されています。
XLBによるLBM計算を行い、解析結果データとして「速度・密度等」を計算します。
これを、「xxxx.h5：データ本体」「xxxx.xmf：データの説明書」として保存します。
これを、ParaViewの可視化ツールで読み込み、「速度・圧力・渦など」を３次元可視化します。
簡単に言うと、.h5 が実データ、.xmf がその実データをParaViewなどに読ませるための案内ファイルです。
【.h5 ファイルとは】
.h5 は HDF5（Hierarchical Data Format 5） 形式です。HDF5は、大量の数値データを効率よく保存するための形式です。Ahmedボディのような3次元LBMでは、格子点が非常に多くなります。そこで、効率的にデータを保存するための、階層構造を持つデータベースになっています。
【.xmf ファイルとは】
.xmf は XDMF（eXtensible Data Model and Format） 形式のファイルです。こちらは非常に重要な機能があります。つまり、XMF自身に巨大な解析データを全部入れるのではなく、以下のような形で、HDF5データへの「案内役」 をします。
.xmfファイル
├─ 格子の大きさ・格子の間隔・座標・どんな物理量があるか・データの次元
└─ そのデータはどのH5に入っているか ⇒ .h5ファイル
このように、２つのファイルに分けるのは、XML形式の.xmfファイルは、可読文字のアスキーデータですが、この形式で膨大なデータを直接記述すると、ファイルが非常に大きくなります。
そこで、バイナリデータとして効率よく記述できる.h5ファイルにデータそのものは記述し、その定義のみを.xmfファイルに記述しています。
【ParaViewではどう使うのか】
通常は **.xmfをParaViewで開きます。**つまり、WindowsのParaViewを起動して、XLBneonのUbuntuファイルシステムに接続して、「XLBneonWork/MWIND3D/Ahmed/D-h5+xmf」にある、ファイル「Ahmed_...xmf」を選択して、開けば可視化できるはずです。
しかし、自分の環境では、ParaViewが異常終了したり、表示ができなかったです。このファイルを使ったParaViewの可視化は、追ってご報告します。
実験との比較検討用csvデータとグラフ画像の分析
ディレクトリ「D-volplot」には、Ahmedボディ周りの速度場に対して、実験結果とこのXLBによるLBM解析の結果を比較します。
実験結果は、ディレクトリ「MWIND3D」にあるファイル「ahmed.json」で、スクリプトの中で比較用データとして、設定しています。
解析結果は、ディレクトリ「D-volplot」にある拡張子.csvのファイルで、指定する位置に対応して、複数のファイルがあります。
この２つの結果を、指定する位置に対応したグラフ画像を作り、拡張子.pngファイルとして保存してあります。
以下の11個のグラフが、ボディ後部の座標原点から、前方向（X座標ー方向）に進んだ評価点の、高さ方向の速度値の、「実験が水色丸、解析が橙色線」として、表示しています。ボディ部分からその前部は、実験と解析の相違が大きいです。
今回は解析実現を主目的にとしており、解析精度を実現することは今後の課題にしたいと思います。
以下の６個のグラフが、ボディ後部の座標原点から、後方向（X座標＋方向）に進んだ評価点の、高さ方向の速度値の、「実験が水色丸、解析が橙色線」として、表示しています。ボディ部分からその後部は、実験と解析の相違が小さく、傾向は類似しています。
どちらの結果もLBMの解析において、位置による速度値のブレが大きいので、解析条件などの工夫が必要と思われます。
実験との比較検討用csvデータとグラフ画像の分析
ディレクトリ「D-volmagpng」には、Ahmedボディ周りの速度場に対して、縦方向の中央断面の画像が保存してあります。
解析は０から4160ステップまで実行し、間隔40ステップで画像が作られています。ただし、最初の0000,0041ステップは、流体が挙動していないので、画像は無視してください。
初めのほうの、0082-0697ステップの状態は、以下の通りです。
終わりのほうの、3256-4141ステップの状態は、以下の通りです。
以上で、マルチ解像度格子を用いたNeonバックエンドのXLBによるLBM解析の、スクリプトの解説と解析結果の分析 を終わります。