AutodeskのXLBには、三角形メッシュを流体領域内の境界として扱う機能があります。この機能を使うと、CADから書き出したSTLを読み込み、モデル周りの3次元流れを格子ボルツマン法（LBM）で計算できます。
本稿では、次の一連の流れを説明します。
CADモデル
↓ STLへ書き出し
Trimeshで読み込み・検査
↓ 三角形ごとの頂点列へ変換
LBM格子単位へ拡大縮小・配置
↓ XLBのHalfwayBounceBackBCへ渡す
3次元風洞シミュレーション
↓
境界マスク・速度場をVTKへ出力
↓
Windows版ParaViewで確認
最初から大きな自動車モデルを使うと、STLの問題、配置の問題、GPUメモリ不足、計算条件の問題を切り分けにくくなります。そこで、まず小さな球STLで処理経路を確認し、その後に任意のCADモデルへ差し替えます。
この記事でできること
- CADから出力したSTLをTrimeshで読み込む
- STLが閉じた形状か、法線の向きが整っているかを確認する
- STLをXLBが扱える三角形頂点列へ変換する
- CAD座標をLBM格子座標へ拡大縮小し、風洞内へ配置する
- STL表面へ滑りなし境界条件を設定する
- 流入口、流出口、風洞壁を設定して3次元計算を実行する
- STLが境界として認識されたかをVTKで確認する
- 密度、速度成分、速度の大きさをVTKへ保存する
- ParaViewでSTL、境界マスク、流れ場を重ねて確認する
検証基準と重要な前提
本稿は2026年9月1日時点のXLB mainブランチ、次のコミットを基準にしています。
9470e54a8d7ccd68d8e5563ca7a573040841ea8c
この時点のXLBでは、STLなどのメッシュを境界マスクへ変換するMeshBoundaryMaskerは、次の条件で実装されています。
- 3次元速度格子のみ
-
ComputeBackend.WARPまたはComputeBackend.NEON -
ComputeBackend.JAXではメッシュ境界生成は未実装
本稿では構成を単純にするため、NVIDIA Warpを使う単一GPU構成にします。JAX/CPUだけで実行したLid-driven cavityの環境では、そのままSTL境界を使用できない点に注意してください。
また、XLBの公式windtunnel_3d.pyは、既定値が512 × 128 × 128格子、100,000ステップです。これは動作確認としては大きいため、本稿では128 × 64 × 64格子、500ステップへ縮小します。
本稿の条件は、CAD→STL→XLB→VTKという処理経路を確認するためのものです。実機相当の空力値を得るには、物理単位変換、格子収束性、計算領域サイズ、境界条件、乱流モデル、時間収束、実験値との比較などを別途検討する必要があります。
STLを使う前に理解しておきたいこと
STLはCADのソリッドそのものではない
STLに保存されるのは、CADの履歴、拘束、材質、曲面式ではなく、表面を近似する三角形です。曲面の滑らかさは、STLへ書き出すときの三角形分割精度で決まります。
ただし、STLだけを極端に細かくしても、XLB側の格子が粗ければシミュレーション上の形状精度は上がりません。STLの三角形数と、モデルを何個のLBMセルで表現するかは別の解像度です。
STLには単位情報がない
STLの座標は任意単位で、ファイル自体にはmm、cm、mといったスケール情報が入りません。本稿のコードはモデル最大寸法を指定セル数へ正規化するため、最初の形状確認では元の単位に依存しません。
物理単位と対応させる段階では、CAD書き出し時の単位を把握し、ボクセル寸法や代表速度を明示してください。XLBには物理単位とLBM単位を変換するUnitConvertorもあります。
できるだけ閉じたメッシュを使う
風洞内の固体として扱うSTLは、次の条件を満たすものが扱いやすくなります。
- 穴のない閉じた表面である
- 非多様体エッジがない
- 重複面や面積ゼロの三角形がない
- 三角形の向きが一貫している
- 外向き法線になっている
- 意図しない離れた部品が含まれていない
Trimeshではis_watertight、is_winding_consistent、is_volumeを確認できます。
なぜmesh.verticesをそのまま渡さないのか
ここは任意のCAD STLへ差し替えるうえで重要です。
現在のXLBのメッシュ境界実装は、受け取ったmesh_verticesに対して内部で次のようなインデックスを作ります。
mesh_indices = np.arange(mesh_vertices.shape[0])
つまり、入力配列の0、1、2番目が最初の三角形、3、4、5番目が次の三角形、という並びを前提としています。
一方、一般的なTrimesh.verticesは、複数の面から共有される頂点テーブルです。頂点0、1、2が必ず同じ面を構成するとは限りません。そのため本稿では、面情報を展開した次の配列を渡します。
model_vertices = np.asarray(mesh.triangles, dtype=np.float32).reshape(-1, 3)
mesh.trianglesは(面数, 3頂点, XYZ)なので、reshape(-1, 3)すると「三角形ごとに3頂点が連続する配列」になります。
process=Falseとwatertight判定
STLは各三角形が独立した3頂点を持つ形で保存されることがあります。process=Falseで生データを保持すると、座標が同じ頂点も統合されません。そのままis_watertightを調べると、本当は閉じた形状でもFalseになることがあります。
そこで、境界へ渡す元データは変更せず、検査用コピーだけ頂点を統合します。
mesh = trimesh.load_mesh(STL_FILE, process=False)
check_mesh = mesh.copy()
check_mesh.merge_vertices()
print(check_mesh.is_watertight)
この方法なら、STL内の三角形記録を保ったまま、形状の接続性を確認できます。
1. XLBのWarp環境を準備する
既にXLBをWSL2へインストール済みの場合も、記事と同じAPIへ固定したいときは別の仮想環境を作ると安全です。
cd ~/src
git clone https://github.com/Autodesk/XLB.git
cd XLB
git checkout 9470e54a8d7ccd68d8e5563ca7a573040841ea8c
python3 -m venv .venv
source .venv/bin/activate
python -m pip install --upgrade pip setuptools wheel
python -m pip install -e ".[warp]"
XLB、Warp、JAX、GPUデバイスを確認します。
python - <<'PY'
import jax
import warp as wp
import xlb
wp.init()
print("XLB :", xlb.__version__)
print("Warp:", wp.__version__)
print("JAX :", jax.__version__)
print("Warp device:", wp.get_device())
PY
NVIDIA GPUを使用できる場合、Warp deviceは環境に応じて次のように表示されます。
cuda:0
cpuになっている場合でもWarp自体は動く可能性がありますが、3次元風洞計算は非常に遅くなります。WSL2からnvidia-smiが実行できること、Windows側のNVIDIAドライバー、WSL2のGPU連携を先に確認してください。
2. 作業ディレクトリを作る
XLBをeditable installしていれば、計算スクリプトはXLBリポジトリ外へ置けます。
mkdir -p ~/src/xlb-cad-stl/models
cd ~/src/xlb-cad-stl
CADから出力したSTLを配置する場合は、たとえばWindowsのダウンロードフォルダからコピーします。
cp /mnt/c/Users/<Windowsユーザー名>/Downloads/my_model.stl models/
計算時に頻繁に読むファイルは、/mnt/cではなくWSL2側へコピーしておく方が扱いやすくなります。
3. 最初の確認用STLを生成する
任意のCADモデルを使う前に、Trimeshで閉じた球STLを生成します。
python - <<'PY'
from pathlib import Path
import trimesh
Path("models").mkdir(exist_ok=True)
mesh = trimesh.creation.icosphere(subdivisions=2, radius=1.0)
mesh.export("models/sample_sphere.stl")
print("written: models/sample_sphere.stl")
print("triangles:", len(mesh.faces))
print("watertight:", mesh.is_watertight)
PY
次のように表示されれば準備完了です。
written: models/sample_sphere.stl
triangles: 320
watertight: True
4. STL風洞計算スクリプトを作る
cad_windtunnel_minimal.pyを作成し、次のコードを保存します。
from pathlib import Path
import time
import numpy as np
import trimesh
import warp as wp
import xlb
from xlb.compute_backend import ComputeBackend
from xlb.grid import grid_factory
from xlb.operator.boundary_condition import (
ExtrapolationOutflowBC,
FullwayBounceBackBC,
HalfwayBounceBackBC,
RegularizedBC,
)
from xlb.operator.boundary_masker import MeshVoxelizationMethod
from xlb.operator.macroscopic import Macroscopic
from xlb.operator.stepper import IncompressibleNavierStokesStepper
from xlb.precision_policy import PrecisionPolicy
from xlb.utils import save_fields_vtk, save_image, warp_array_to_jax
# -----------------------------------------------------------------------------
# User settings
# -----------------------------------------------------------------------------
STL_FILE = Path("models/sample_sphere.stl")
OUTPUT_DIR = Path("output/cad_windtunnel")
GRID_SHAPE = (128, 64, 64)
MODEL_LENGTH_LBM = 20.0 # Model maximum extent in lattice cells
WIND_SPEED = 0.03
REYNOLDS_NUMBER = 50.0
NUM_STEPS = 500
OUTPUT_INTERVAL = 250
OUTPUT_DIR.mkdir(parents=True, exist_ok=True)
# -----------------------------------------------------------------------------
# 1. Load and inspect the STL
# -----------------------------------------------------------------------------
mesh = trimesh.load_mesh(STL_FILE, process=False)
if not isinstance(mesh, trimesh.Trimesh) or mesh.is_empty:
raise ValueError(f"A valid triangle mesh could not be loaded: {STL_FILE}")
# Raw STL often has duplicated vertices per triangle. Merge only on a copy for
# topology checks; keep the original triangle records for XLB.
check_mesh = mesh.copy()
check_mesh.merge_vertices()
print(f"triangles : {len(mesh.faces):,}")
print(f"watertight : {check_mesh.is_watertight}")
print(f"winding consistent : {check_mesh.is_winding_consistent}")
if not check_mesh.is_watertight:
raise ValueError("The STL is not watertight. Repair holes/non-manifold edges first.")
# Optional orientation adjustment, for example 90 degrees about X:
# mesh.apply_transform(
# trimesh.transformations.euler_matrix(np.deg2rad(90), 0, 0, axes="sxyz")
# )
# -----------------------------------------------------------------------------
# 2. Convert to triangle soup, scale, and place in lattice coordinates
# -----------------------------------------------------------------------------
# XLB currently constructs triangle indices as 0, 1, 2, ... internally.
# Therefore pass consecutive triples of vertices, not an indexed vertex table.
model_vertices = np.asarray(mesh.triangles, dtype=np.float32).reshape(-1, 3).copy()
model_vertices -= model_vertices.min(axis=0)
physical_extents = np.ptp(model_vertices, axis=0)
scale = MODEL_LENGTH_LBM / physical_extents.max()
model_vertices *= np.float32(scale)
lbm_extents = np.ptp(model_vertices, axis=0)
# Flow is +X. Leave an upstream region and center the model in Y/Z.
translation = np.array(
[
GRID_SHAPE[0] * 0.25,
(GRID_SHAPE[1] - lbm_extents[1]) / 2,
(GRID_SHAPE[2] - lbm_extents[2]) / 2,
],
dtype=np.float32,
)
model_vertices += translation
model_vertices = np.ascontiguousarray(model_vertices)
mesh_min = model_vertices.min(axis=0)
mesh_max = model_vertices.max(axis=0)
if np.any(mesh_min < 3) or np.any(mesh_max >= np.asarray(GRID_SHAPE) - 3):
raise ValueError(
f"The placed model does not fit in the domain: min={mesh_min}, max={mesh_max}"
)
# Export the transformed STL so it can be checked together with VTK in ParaView.
faces = np.arange(len(model_vertices), dtype=np.int64).reshape(-1, 3)
trimesh.Trimesh(model_vertices, faces, process=False).export(
OUTPUT_DIR / "placed_model_lbm.stl"
)
# -----------------------------------------------------------------------------
# 3. Initialize XLB and define the wind-tunnel boundary conditions
# -----------------------------------------------------------------------------
backend = ComputeBackend.WARP
precision = PrecisionPolicy.FP32FP32
velocity_set = xlb.velocity_set.D3Q27(
precision_policy=precision,
compute_backend=backend,
)
xlb.init(
velocity_set=velocity_set,
default_backend=backend,
default_precision_policy=precision,
)
grid = grid_factory(GRID_SHAPE, compute_backend=backend)
box = grid.bounding_box_indices()
box_no_edge = grid.bounding_box_indices(remove_edges=True)
inlet = box_no_edge["left"]
outlet = box_no_edge["right"]
walls = [
box["bottom"][axis]
+ box["top"][axis]
+ box["front"][axis]
+ box["back"][axis]
for axis in range(velocity_set.d)
]
walls = np.unique(np.asarray(walls), axis=-1).tolist()
bc_inlet = RegularizedBC(
"velocity",
prescribed_value=(WIND_SPEED, 0.0, 0.0),
indices=inlet,
)
bc_outlet = ExtrapolationOutflowBC(indices=outlet)
bc_walls = FullwayBounceBackBC(indices=walls)
bc_model = HalfwayBounceBackBC(
mesh_vertices=model_vertices,
voxelization_method=MeshVoxelizationMethod("RAY"),
)
stepper = IncompressibleNavierStokesStepper(
grid=grid,
boundary_conditions=[bc_walls, bc_inlet, bc_outlet, bc_model],
collision_type="KBC",
)
f_0, f_1, bc_mask, missing_mask = stepper.prepare_fields()
wp.synchronize()
# Confirm that the STL produced mesh-boundary cells. A zero count means that the
# geometry, placement, voxelization, or XLB/Warp combination must be checked.
bc_ids = np.asarray(bc_mask.numpy())[0]
model_bc_cells = int(np.count_nonzero(bc_ids == bc_model.id))
print(f"model BC id : {bc_model.id}")
print(f"model BC cells : {model_bc_cells:,}")
save_fields_vtk(
{"bc_id": bc_ids},
timestep=0,
output_dir=str(OUTPUT_DIR),
prefix="boundary_mask",
)
if model_bc_cells == 0:
raise RuntimeError("No STL-derived boundary cells were generated.")
# Re = U L / nu, and nu = (1 / omega - 0.5) / 3 in lattice units.
nu = WIND_SPEED * MODEL_LENGTH_LBM / REYNOLDS_NUMBER
omega = 1.0 / (3.0 * nu + 0.5)
print(f"nu={nu:.8g}, omega={omega:.8g}")
macro = Macroscopic(
compute_backend=ComputeBackend.JAX,
precision_policy=precision,
velocity_set=xlb.velocity_set.D3Q27(
precision_policy=precision,
compute_backend=ComputeBackend.JAX,
),
)
# -----------------------------------------------------------------------------
# 4. Run and write VTK/PNG output
# -----------------------------------------------------------------------------
def write_output(step: int) -> None:
wp.synchronize()
populations = warp_array_to_jax(f_0)
rho, velocity = macro(populations)
rho_np = np.asarray(rho[0])
velocity_np = np.asarray(velocity)
speed_np = np.sqrt(np.sum(velocity_np**2, axis=0))
save_fields_vtk(
{
"rho": rho_np,
"u_x": velocity_np[0],
"u_y": velocity_np[1],
"u_z": velocity_np[2],
"u_magnitude": speed_np,
},
timestep=step,
output_dir=str(OUTPUT_DIR),
prefix="cad_windtunnel",
)
mid_y = speed_np.shape[1] // 2
save_image(
speed_np[:, mid_y, :],
timestep=step,
prefix=str(OUTPUT_DIR / "cad_windtunnel_mid_y"),
)
started = time.perf_counter()
for step in range(NUM_STEPS):
f_0, f_1 = stepper(f_0, f_1, bc_mask, missing_mask, omega, step)
f_0, f_1 = f_1, f_0
if step % 100 == 0:
wp.synchronize()
print(f"step={step:5d}, elapsed={time.perf_counter() - started:.2f}s")
started = time.perf_counter()
if step % OUTPUT_INTERVAL == 0 or step == NUM_STEPS - 1:
write_output(step)
print("Simulation completed successfully.")
5. 球STLで動作確認する
仮想環境が有効であることを確認して実行します。
cd ~/src/xlb-cad-stl
source ~/src/XLB/.venv/bin/activate
python cad_windtunnel_minimal.py
最初にSTLの検査結果が表示されます。
triangles : 320
watertight : True
winding consistent : True
続いて、XLBがSTLから作成した境界IDとセル数が表示されます。
model BC id : <数値>
model BC cells : <0より大きい数値>
model BC cellsが0ではないことが重要です。
計算が進むと、100ステップごとに経過時間が表示され、0、250、499ステップでVTKとPNGが保存されます。
Warpの初回実行ではカーネルのロードやコンパイルが入るため、2回目以降と比べて準備に時間がかかることがあります。
6. 出力ファイルを確認する
正常終了すると、次のようなファイルが生成されます。
output/cad_windtunnel/
├── placed_model_lbm.stl
├── boundary_mask_0000000.vtk
├── cad_windtunnel_0000000.vtk
├── cad_windtunnel_0000250.vtk
├── cad_windtunnel_0000499.vtk
├── cad_windtunnel_mid_y_0000.png
├── cad_windtunnel_mid_y_0250.png
└── cad_windtunnel_mid_y_0499.png
それぞれの役割は次のとおりです。
| ファイル | 内容 |
|---|---|
placed_model_lbm.stl |
拡大縮小・移動後のSTL。座標はLBM格子単位 |
boundary_mask_0000000.vtk |
XLBが生成した境界IDの3次元分布 |
cad_windtunnel_*.vtk |
密度、XYZ速度成分、速度の大きさ |
cad_windtunnel_mid_y_*.png |
Y中央断面の速度の大きさ |
WSL2の出力フォルダをWindowsのエクスプローラーで開きます。
explorer.exe output/cad_windtunnel
7. ParaViewでSTL境界を確認する
流れを見る前に、STLがXLBの格子境界として認識されているかを確認します。
配置後STLを開く
Windows版ParaViewで次のファイルを開きます。
placed_model_lbm.stl
Applyを押し、モデルの向きと位置を確認します。流れ方向は+Xです。
境界マスクを開く
次に、次のファイルを開きます。
boundary_mask_0000000.vtk
Applyを押し、Coloringをbc_idへ変更します。
外壁、流入口、流出口、STL境界は異なるIDを持ちます。ターミナルへ表示されたmodel BC idを使い、ParaViewのThresholdフィルターでその値だけを抽出すると、STL表面に対応する境界セルを確認できます。
Filters → Threshold
Scalars: bc_id
Lower Threshold: model BC id
Upper Threshold: model BC id
境界セルが形状に沿って現れれば、STLの読み込みとボクセル境界生成は成功です。
8. ParaViewで速度場を確認する
次の連番VTKをファイルシリーズとして開きます。
cad_windtunnel_...vtk
Apply後、Coloringをu_magnitudeへ変更します。
3次元データ全体をSurface表示すると外側しか見えないため、Sliceフィルターを追加します。
Filters → Slice
たとえばY中央断面を見る場合は、法線をY方向へ設定し、原点を領域中央へ置きます。断面をu_magnitudeで着色すると、モデル前方の減速域や後方の後流を確認できます。
VTKには次の配列を保存しています。
| 配列 | 意味 |
|---|---|
rho |
密度 |
u_x |
X方向速度 |
u_y |
Y方向速度 |
u_z |
Z方向速度 |
u_magnitude |
速度の大きさ |
9. 自分のCADモデルへ差し替える
球モデルで一連の処理を確認できたら、スクリプト冒頭のSTLパスを変更します。
STL_FILE = Path("models/my_model.stl")
向きを変更する
CADの前後方向とXLBの流れ方向が一致しない場合は、mesh.trianglesへ変換する前に回転を適用します。
たとえばX軸周りへ90度回転する場合は、スクリプト内のコメントを外します。
mesh.apply_transform(
trimesh.transformations.euler_matrix(
np.deg2rad(90),
0,
0,
axes="sxyz",
)
)
複数軸を回転する場合は、X、Y、Zの角度を指定します。
mesh.apply_transform(
trimesh.transformations.euler_matrix(
np.deg2rad(90),
np.deg2rad(0),
np.deg2rad(-90),
axes="sxyz",
)
)
回転後は必ずplaced_model_lbm.stlをParaViewで開き、前後、上下、左右が意図どおりか確認してください。
モデルを大きくする
MODEL_LENGTH_LBM = 40.0
この値は、モデルの最大寸法を何個のLBMセルへ割り当てるかを表します。
- 20セル程度：処理経路を確認するための粗い設定
- 40～100セル程度：形状をもう少し表現したい場合の出発点
- それ以上：GPUメモリと計算時間を確認しながら調整
これは精度を保証する基準ではありません。結果に必要な解像度は形状、Reynolds数、評価量によって異なります。最終的には格子数を段階的に増やし、結果の変化が十分小さくなることを確認します。
計算領域も拡大する
モデルだけを大きくすると、壁との距離や下流長さが不足します。
GRID_SHAPE = (256, 128, 128)
MODEL_LENGTH_LBM = 48.0
D3Q27では、分布関数の主要2配列だけでも、おおよそ次のメモリを使います。
2 × 27 × NX × NY × NZ × 4 byte
実際には境界マスク、中間配列、Warp/JAX間の変換、VTK出力用配列も必要です。GPUメモリぎりぎりではなく、余裕を残して設定してください。
地面へ近づける
本稿の最小版は、モデルをY方向とZ方向の中央へ置きます。自動車のように地面へ近づけたい場合は、Z方向の移動量を変更します。
translation = np.array(
[
GRID_SHAPE[0] * 0.25,
(GRID_SHAPE[1] - lbm_extents[1]) / 2,
3.0,
],
dtype=np.float32,
)
ただし、モデルと底面境界が交差すると境界条件が重なります。最初は数セルの隙間を設け、境界マスクを確認してから調整してください。公式windtunnel_3d.pyでも、ボクセル化方式に応じて底面からのシフト量を変えています。
10. 計算条件の意味
流入速度
WIND_SPEED = 0.03
これはm/sではなくLBM単位の速度です。最初の確認では小さな値から始めます。値を大きくしすぎると、非圧縮近似や数値安定性に影響します。
Reynolds数と緩和係数
本稿ではモデル最大寸法を代表長さLとし、次の関係を使います。
Re = U L / ν
ν = U L / Re
ω = 1 / (3ν + 0.5)
コードでは次の部分です。
nu = WIND_SPEED * MODEL_LENGTH_LBM / REYNOLDS_NUMBER
omega = 1.0 / (3.0 * nu + 0.5)
omegaが2へ極端に近づく条件は不安定になりやすいため、発散する場合は次を試します。
- Reynolds数を下げる
- モデルを表すセル数を増やす
- 流入速度を下げる
-
BGKではなくKBCを使う - STLや境界マスクに穴がないか確認する
500ステップは完成した流れ場ではない
本稿の500ステップは、計算が開始し、VTKが生成されることを確認するための設定です。十分に発達した後流や、収束した抗力係数を求めるには、計算領域と流速に応じてステップ数を増やす必要があります。
11. XLBのボクセル化方式を切り替える
現在のXLBには次の方式が登録されています。
| 指定 | 概要 | 主な注意点 |
|---|---|---|
RAY |
格子方向へレイを飛ばし、メッシュとの交差を調べる | 公式風洞exampleの既定方式 |
AABB |
三角形と単位ボクセルの交差をAABB検索で調べる | おおよそ1セル厚の表面を作る |
AABB_CLOSE |
AABB結果へ膨張・収縮のclose処理を加える | 小さな穴を埋めやすいが形状も変わり得る |
WINDING |
winding numberを使って形状を判定する | 法線・面の向きが重要。Neonでは未実装 |
本稿はRAYを使います。
voxelization_method=MeshVoxelizationMethod("RAY")
AABBへ変更する場合は次のようにします。
voxelization_method=MeshVoxelizationMethod("AABB")
小さな隙間をclose処理で埋める場合は、処理幅も指定します。
voxelization_method=MeshVoxelizationMethod(
"AABB_CLOSE",
close_voxels=2,
)
方式を変更した後は、必ずboundary_mask_0000000.vtkを確認してください。穴を埋める処理は、意図した開口部まで閉じる可能性があります。
12. よくあるエラーと確認方法
MeshBoundaryMasker is only implemented for WARP and NEON
STL境界をJAXバックエンドで作ろうとしています。
backend = ComputeBackend.WARP
になっていることを確認します。JAXは本稿では後処理の巨視的量計算に使っていますが、STL境界生成と時間発展はWarpで実行します。
MeshBoundaryMasker is only implemented for 3D velocity sets
D2Q9などの2次元速度格子を使用しています。3次元STLにはD3Q19またはD3Q27が必要です。本稿はKBCと組み合わせるためD3Q27を使用します。
The mesh must be fully contained within the domain
配置後STLが格子外へ出ています。
-
MODEL_LENGTH_LBMを小さくする -
GRID_SHAPEを大きくする - X方向の配置位置を変更する
- 回転後の寸法を確認する
placed_model_lbm.stlをParaViewで開くと、配置を視覚的に確認できます。
watertight: False
CADまたはメッシュ修復ツールで次を確認します。
- 穴
- T字接続
- 非多様体エッジ
- 重複面
- 内部面
- 離れた小部品
- 面積ゼロの三角形
process=Falseによる未統合頂点だけが原因の場合もあるため、本稿のように検査用コピーへmerge_vertices()を適用してから判断します。
model BC cells : 0
STLから境界が作成されていません。次の順で確認します。
-
placed_model_lbm.stlが計算領域内にあるか -
mesh.triangles.reshape(-1, 3)を渡しているか - 三角形数が0ではないか
-
MODEL_LENGTH_LBMが小さすぎないか -
RAYからAABBへ切り替えるとどうなるか - XLBと
warp-langのバージョンを記録する -
boundary_mask_0000000.vtkにモデル境界IDがあるか
境界セル数を最初に検査することで、モデルがない状態のまま長時間計算してしまうことを防げます。
途中でNaNやInfが出る
- Reynolds数を下げる
- 流入速度を下げる
- 格子解像度を上げる
- STLの穴や境界の重なりを確認する
- モデルと流入口・流出口・外壁の距離を増やす
-
KBCを使用する
密度rhoの最小値と最大値も監視すると、発散の兆候を見つけやすくなります。
GPUメモリ不足になる
最初は次の設定へ戻します。
GRID_SHAPE = (128, 64, 64)
MODEL_LENGTH_LBM = 20.0
NUM_STEPS = 500
VTK出力時にはWarpの分布関数をJAX/CPU側へ変換するため、一時的なメモリも必要です。出力間隔を広げ、不要な配列をVTKへ保存しないことも有効です。
XLB 0.3.1でメッシュが消える
XLBのIssue #148には、XLB 0.3.1と新しいWarpの組み合わせでメッシュが消えるという報告があります。2026年5月14日のコメントでは、XLB 0.3.1に対してwarp-lang==1.7.2へ戻す切り分けが提案され、Warp 1.11.1では症状が再現したとされています。
ただし、本稿が基準にしているmainコミットは、0.3.1からAPIと依存条件が変わっています。本稿の環境へ0.3.1向け回避策をそのまま混ぜないでください。
旧リリースを検証するときは別の仮想環境を作り、公式0.3.1 exampleと組み合わせて切り分けます。本稿のスクリプトは、現在のMeshVoxelizationMethod APIを使うため、0.3.1向けではありません。
13. コマンドライン版スクリプトへ発展させる
複数のSTLや条件を試す場合は、設定値をソースコードへ直接書くより、次のようなコマンドライン引数へすると便利です。
python cad_windtunnel.py \
--stl models/my_model.stl \
--grid 192 96 96 \
--model-length 32 \
--rotate-x 90 \
--wind-speed 0.03 \
--reynolds 100 \
--steps 2000 \
--output-interval 500 \
--voxelization RAY \
--device cuda:0
計算前の形状検査だけを行うモードも有用です。
python cad_windtunnel.py \
--stl models/my_model.stl \
--grid 192 96 96 \
--model-length 32 \
--inspect-only
この段階で配置後STLをParaViewへ読み込み、向きとサイズを確認してからGPU計算へ進めると、試行錯誤を減らせます。
14. 抗力・揚力を計算するには
XLBの公式examples/cfd/windtunnel_3d.pyでは、MomentumTransferを使ってメッシュ境界へ作用する力を求めています。
概念的には次の流れです。
from xlb.operator.force.momentum_transfer import MomentumTransfer
momentum_transfer = MomentumTransfer(
bc_model,
compute_backend=ComputeBackend.WARP,
)
boundary_force = momentum_transfer(
f_0,
f_1,
bc_mask,
missing_mask,
)
ただし、抗力係数・揚力係数へ変換するには、代表面積、代表速度、密度、座標軸、時間平均区間を適切に定義する必要があります。まず本稿の手順で、形状と境界マスク、速度場が妥当であることを確認してから追加するのが安全です。
まとめ
CADのSTLをXLBへ渡す処理は、単にtrimesh.load_mesh()を呼ぶだけではありません。特に重要なのは次の点です。
- STLを検査し、閉じた三角形表面を用意する
-
mesh.triangles.reshape(-1, 3)で三角形ごとの頂点列へ展開する - CAD座標をLBM格子単位へ変換し、領域内へ配置する
- Warpバックエンドと3次元速度格子を使う
-
HalfwayBounceBackBCへSTL頂点とボクセル化方式を渡す - 計算前にモデル境界セル数と
boundary_maskを確認する - VTKと配置後STLをParaViewで重ねて確認する
最初は球などの単純形状・小規模格子・低いReynolds数で処理経路を確立し、その後にCADモデル、格子解像度、風洞寸法、物理条件を段階的に引き上げると、問題を切り分けやすくなります。