LBM流体解析ツールXLBの活用【10:XLB独自例題の作成と分析】
この文書では、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計算を比較的容易に実装できます。
このXLBのgithubリポジトリは、以下にあります。
Autodesk/XLB
ここでは、最も基本となるCPUで動作する解析環境「JAX-CPU」を、以下の手順で準備した状態を用いて、独自例題の作成と分析を行います。
qiita.com
対象とする独自例題の概要
ここでは、「２次元問題として風洞中にある矩形周りの流れ」 を分析することを、目標とします。例えば、以下の論文では、格子ボルツマン法と有限体積法を用いた解析が報告されており、これを参考にします。
「Accurate computations of the laminar flow past a square cylinder based on two different methods: lattice-Boltzmann and finite-volume」
この論文中の図２では、以下の結果が示されています。
!image.png
Fig. 2. Streamlines around the square cylinder for different Reynolds numbers (a) Re = 1; (b) Re = 30; (c) Re = 60; (d) Re = 200
解析モデルと解析条件と境界条件などの設定
これより、レイノルズ数Reの違いによって、正方形周りの流れの変化を、比較検討できます。解析モデルと境界条件は、以下の図のような定義です。
!image.png
独自例題の、格子ボルツマン法LBMをAutodesk XLBで実施するために関する設定は、以下の通りです。
| 項目 | 設定 | 項目 | 設定 |
|--|--|--|--|
| 計算領域 | 2次元長方形| 流体モデル | D2Q9 |
| 衝突モデル | BGK| 流入 | 一定速度 |
| 流出 | ExtrapolationOutflowBC| 上下壁 | HalfwayBounceBackBC |
| 正方形物体 | HalfwayBounceBackBC| 計算バックエンド | JAX |
| GPU(Jax,Warp,Neon) | 使用しない| CPU | JAX-CPU |
独自例題の具体的な設計
以下では、独自例題のスクリプトを作るうえで、設定する項目をまとめておきます。
スクリプト内の変数名と上記設定値との対応
独自例題での設定値とスクリプト内の 「変数名との対応」 を、下表に示します。問題を簡単にするため、格子や物体などの寸法は、無次元化した数値を用いています。
| 項目 | 説明図の内容 | スクリプトの変数 | 設定値 | 単位・意味 |
|---|---|---|---:|---|
| 計算領域・X方向 | 領域の長さ | Nx | 300 | 格子点 |
| 計算領域・Y方向 | 領域の高さ | Ny | 120 | 格子点 |
| 流入速度 | X方向流入速度 | u_in | 0.05 | 格子単位/step |
| レイノルズ数 | 流れの無次元数 | Re | 100.0 | 無次元 |
| 正方形一辺 | 物体の大きさ | square_size | 40 | 格子点 |
| 正方形中心 X | 物体中心のX座標 | square_cx | 75 | 格子座標 |
| 正方形中心 Y | 物体中心のY座標 | square_cy | 60 | 格子座標 |
| 正方形回転角 | 物体の傾き | square_angle | 30.0 | 度 |
| 動粘性係数 | LBMの粘性 | visc | 0.020 | 格子単位 |
| 緩和パラメータ | BGKの緩和係数 | omega | 約1.7857 | 無次元 |
| 計算ステップ数 | 時間発展回数 | num_steps | 5000 | step |
| 出力間隔 | 結果を保存する間隔 | post_process_interval | 500 | step |
さらに重要なのは、「計算で決まる数値」 です。例えば、動粘性係数：visc と BGK緩和パラメータ：omega は直接入力しているのではなく、上の数値から計算されています。
visc = u_in * square_size / Re
omega = 1.0 / (3.0 * visc + 0.5)
したがって、$\nu=\frac{0.05\times40}{100}=0.020$、そして、$\omega=\frac{1}{3\times0.020+0.5}
\approx1.7857$です。
つまり、上の方のユーザーが変更する基本パラメータ*（設定値）と、そこから自動的に計算されるパラメータ（計算値）を分けると理解しやすくなります。この表にすると、「どの数値を変更すると、解析モデルの何が変わるのか」がかなり分かりやすくなります。
| 種類 | 変数 | 値 | 役割 |
|---|---|---:|---|
| 計算値 | visc | 0.020 | 動粘性係数 |
| 計算値 | omega | 1.7857 | BGK緩和パラメータ |
以上の設定により、下図の８の解析結果のイメージを目標とします。
!image.png
独自例題を理解する上でのポイント
この独自例題は、以下の【04:例題flow_past_sphere_3d.pyの分析】に対して、２次元に簡略化した内容です。ただし、流体内部の四角形の「位置・規模・角度」が任意で設定できるので、様々な条件を確認することができます。これを参考にして、「左から流入 → 中央の正方形物体 → 右へ流出」 という、外部流れの最小例にできます。
qiita.com
今回作るモデル
まずは非常に単純化します。設定する流体空間中に、正方形を任意の位置・規模・角度に、境界条件：HalfwayBounceBackとして置くことができます。２次元XY座標において、流体はX方向の＋向きに流れ、左辺が流入（Velocity BC）で右辺が流出（Outflow BC）です。
このプログラムの重要なところ
今回の目的は、XLBの構造をできるだけ単純に理解することです。特に重要なのは次の部分を解説します。
【D2Q9】
2次元LBMの最も基本的な速度モデルです。スクリプトでは以下の定義です。中央の1方向＋周囲8方向なので D2Q9 です。
velocity_set = xlb.velocity_set.D2Q9(
precision_policy=precision_policy,
compute_backend=compute_backend,
)
【JAX-CPU】
今回は計算規模も小さく、どの計算環境でも動作させるため、CPU動作 としています。従って、「Python → XLB → JAX → CPU」のような計算の仕組みになります。
compute_backend = ComputeBackend.JAX
os.environ["JAX_PLATFORMS"] = "cpu"
【正方形】
今回の解析対象の扱いの中心部分です。回転させた正方形内部を座標返還して、正方形内部の格子を判定しています。詳しくはスクリプトを見ながら、説明します。
【正方形の境界条件】
ここがLBMでは重要です。スクリプトでは以下の定義です。以下の設定で、正方形に到達した流体粒子を反射させる というモデルです。XLBの公式2次元キャビティ例でも、固体壁には HalfwayBounceBackBC が使われています。
bc_square = HalfwayBounceBackBC(
indices=square
)
【流入口】
ここでは、流入として「X方向にu_in=0.05・Y方向に0.0」としています。つまり一様流を左辺から入れます。
bc_inlet = RegularizedBC(
"velocity",
prescribed_value=(u_in, 0.0),
indices=inlet,
)
【流出口】
ここでは、流入として、以下のような定義とします。右側では流体を外へ出します。
bc_outlet = ExtrapolationOutflowBC(
indices=outlet
)
【斜めの境界条件】
この例題では、正方形の中心位置はそのままにして、指定した角度だけ回転させるように変更できます。ここでは、正方形そのものを回転させるのではなく、格子上の各セルについて「そのセルが回転した正方形の内部か？」を判定しています。
したがってXLBでは、「数学的な正方形 → 回転座標系で判定 → square_mask → 固体セル → HalfwayBounceBackBC」という流れになります。
なお、LBMの格子は整数セルなので、30°や45°にすると境界が階段状（staircase） になります。これは今回のような基本的なXLB学習例としては重要なポイントです。LBMにおける斜めの境界条件は、階段状のギザギザ（ボクセル表現）による誤差を減らし、滑らかな傾斜や曲面を表現するために非常に重要です。
【Reynolds数】
今回は、「Re = 100.0」としています。代表長さを正方形の一辺「square_size = 20」として、$Re=\frac{UL}{\nu}$から、動粘性係数 visc = u_in * square_size / Re としています。
さらにBGKの緩和係数 $\omega=\frac{1}{3\nu+1/2}$ を、「omega = 1.0 / (3.0 * visc + 0.5)」で計算しています。
なお、正方形の寸法Lを変化させた場合は、本来ならば、レイノルズ数Reが変化するところですが、このスクリプトでは、動粘性係数viscを計算しなおすので、レイノルズ数は変化しない設定となっています。
【公式例題との関係】
特に、現在のXLB公式例では3次元の球まわり流れについて、HalfwayBounceBackBC、RegularizedBC、ExtrapolationOutflowBC、IncompressibleNavierStokesStepper を組み合わせる構成が使われています。今回のスクリプトは、それを2次元・正方形に最小化したものと考えると分かりやすいです。GitHub
なお、現在のXLB main版の公式2D例では save_fields_vtk() と save_image() による結果保存も使われていますので、今回も同じ方式にしています。GitHub
独自例題のスクリプトの解説
以下では、スクリプトの記述に従って、ブロック毎に区分して、それぞれの機能をコメントの形で記入したので、これを用いて説明します。
以下の説明では、ブロックごとに分けて解説しますが、全体をつなげると「２次元空間での正方形周囲の流れ」を、XLB用いてLBM計算するスクリプト「flow_square_2d-1.py」としておきます。
（最後の「-1」は、最後に説明するスクリプトの改良を期待してのバージョン１と考えています）
スクリプト「flow_square_2d-1.py」の構成
このスクリプトは、以下の 「32のブロック」 に、機能が分かれています。各ブロックについて、スクリプトを提示して、コメント部分に説明を追加します。
① JAXをCPUに設定 → ② 計算条件を設定
→ ③ 正方形のサイズ・位置・角度を設定 → ④ XLBのバックエンドをJAXに設定
→ ⑤ D2Q9を設定 → ⑥ 粘性係数・BGK係数を計算
→ ⑦ XLBを初期化 → ⑧ 計算格子を作成
→ ⑨ 境界を取得 → ⑩ 左側 → 流入口
→ ⑪ 右側 → 流出口 → ⑫ 上下 → 固体壁
→ ⑬ 格子座標を作成 → ⑭ 正方形の角度を設定
→ ⑯ 正方形セルを判定 → ⑰ 固体セルのインデックスを作成
→ ⑱ 流入口 BC → ⑲ 上下壁 BC
→ ⑳ 正方形 BC → ㉑ 流出口 BC
→ ㉒ 境界条件を統合 → ㉓ Navier-Stokes Stepper
→ ㉔ LBM分布関数を初期化 → ㉕ 密度・速度計算
→ ㉖ 計算条件を表示 → ㉗ 結果保存関数
→ ㉘ 計算ステップ設定 → ㉙ タイマー開始
→ ㉚ LBM時間発展 → ㉛ 計算終了
→ ㉜ 計算時間表示
#ブロック【① JAXをCPUに設定】の解説：CPUによるXLBによるLBM計算設定
# ============================================================
# XLB : 2D Flow Around a Rotated Square : flow_square_2d-1.py
#
# 2次元の長方形領域の中に、正方形の物体を置き、
# 左から右へ流体を流す外部流れ問題を計算する。
# この正方形は、任意の位置と寸法と角度で置ける。
# 計算方法： LBM (Lattice Boltzmann Method)
# 速度モデル： D2Q9
# 衝突モデル： BGK
# バックエンド： JAX
# 計算デバイス： CPU
# 重要： このプログラムでは GPU(Jax,Warp,Neon) を使用しない。
# ============================================================
# ============================================================
# ① JAXをCPUで使用する設定
# ============================================================
# JAXが利用する計算デバイスをCPUに限定する。
# "cpu"を指定すると、
# Python → XLB → JAX → CPU という計算経路になる。
# NVIDIA GPUやCUDAは、この例では使用しない。
# この環境変数はJAXをimportする前に設定する必要がある。
import os
# 念の為、環境変数JAX_PLATFORMSをcpuに設定する。
os.environ["JAX_PLATFORMS"] = "cpu"
# ============================================================
# Pythonで使用する基本ライブラリ
# ============================================================
# 時間計測に使用する。
# 解析全体に何秒かかったかを最後に表示するために使用する。
import time
# NumPyは、
# ・格子座標の作成 ・正方形の領域判定 ・角度計算
# などに使用する。
import numpy as np
# JAX版NumPy。
# 計算結果の速度や密度など、
# XLB内部でJAX配列として扱われるデータを処理するときに使用する。
import jax.numpy as jnp
# XLB本体を読み込む。
import xlb
# XLBが使用する計算バックエンドを指定するためのクラス。
# 今回は ComputeBackend.JAX を使用する。
from xlb.compute_backend import ComputeBackend
# 数値精度を指定するためのクラス。
# 今回はFP32を使用する。
from xlb.precision_policy import PrecisionPolicy
# 計算格子を作成するための関数。
from xlb.grid import grid_factory
# CFDに対応するNavier-Stokes方程式をLBMとして時間発展させる
# XLBの主要なStepper。これでステップ毎の処理を進める。
from xlb.operator.stepper import (
IncompressibleNavierStokesStepper
)
# 境界条件を指定するクラスを読み込む。
# HalfwayBounceBackBC：
# 固体壁で流体を反射させる。
# RegularizedBC：
# 流入口などで速度を指定する。
# ExtrapolationOutflowBC：
# 流出口から流体を外へ流す。
from xlb.operator.boundary_condition import (
HalfwayBounceBackBC,
RegularizedBC,
ExtrapolationOutflowBC,
)
# 密度や速度などの巨視的物理量を
# 分布関数から計算するクラス。
from xlb.operator.macroscopic import Macroscopic
# 計算結果をVTKや画像として保存する関数。
# VTK： ParaViewなどで解析結果を確認できる。
# PNG： 速度分布などを画像として確認できる。
from xlb.utils import (
save_fields_vtk,
save_image,
)
#ブロック【② 計算条件を設定】の解説：解析空間や流入速度やレイノルズ数の設定
# ============================================================
# ② 計算条件
# ============================================================
# x方向の格子数。今回は300セル。
# 流れは基本的にx方向、つまり「左 → 右」に進む。
# 流入から物体、特に物体から流出においては、十分な空間を確保する。
Nx = 300
# y方向の格子数。今回は120セル。
# 対象の物体の上下に、十分な空間を確保する。
Ny = 120
# XLBに渡す計算領域のサイズ。
# (Nx, Ny) つまり x方向 = 300・y方向 = 120 の2次元格子になる。
# 対象とする空間内の物体の大きさに対応して、計算領域を定める
grid_shape = (Nx, Ny)
# ------------------------------------------------------------
# 流入速度
# ------------------------------------------------------------
# 左側から流入する流体の速度。
# 今回は無次元速度として0.05を使用する。
# x方向： 0.05 y方向： 0
# したがって流体は水平に流入する。
u_in = 0.05
# ------------------------------------------------------------
# Reynolds数
# ------------------------------------------------------------
# 流れの状態を表す無次元数。
# Re = U L / ν
# U： 代表速度 L： 代表長さ ν： 動粘性係数
# 今回は正方形の一辺を代表長さとする。例題ではReを前提条件とする。
Re = 100.0
#ブロック【③ 正方形のサイズ・位置・角度を設定】の解説：中心位置と回転角度の設定
# ============================================================
# ③ 正方形物体の条件
# ============================================================
# 正方形の一辺の長さ。
# 40 × 40セルの正方形を作る。流れに直交方向の格子数Nyと比較する
# 流体空間の格子数 Nx,Ny に適切に入る大きさが限界となる。
# 閉塞率=物体の断面積／空間の断面積が、50%以下になるように注意する。
square_size = 40
# ------------------------------------------------------------
# 正方形の中心位置
# ------------------------------------------------------------
# x方向の中心。
# Nx = 300なので、300 / 4 = 75 となる。
# したがって正方形の中心はx=75。流体空間の流入側の 1/4 に位置する。
# なお // は整数除算なので、割切れないときは、小数点切り捨ての整数化
square_cx = Nx // 4
# y方向の中心。
# Ny = 120なので、120 / 2 = 60 となる。
# したがって正方形の中心はy=60。上下方向の中央に位置する。
square_cy = Ny // 2
# ------------------------------------------------------------
# 正方形の回転角度
# ------------------------------------------------------------
# 正方形を反時計回りに回転させる角度。
# 0度： 通常の正方形
# 30度： 30度回転
# 45度： ひし形に近い形
# 90度： 見た目は0度と同じ
# この値を変更することで、
# 同じ計算プログラムで様々な角度を試すことができる。
# この角度は、細かな数値の任意の数値で設定可能です。
square_angle = 30.0
#ブロック【④ XLBのバックエンドをJAXに設定】の解説：併せて数値制度をFP32に設定
# ============================================================
# ④ XLBの計算バックエンド
# ============================================================
# XLBの計算バックエンドとしてJAXを選択する。
# JAXにはCPU/GPUなどの実行環境があるが、
# 先ほど、os.environ["JAX_PLATFORMS"] = "cpu"
# と指定しているため、今回はCPUで計算される。
compute_backend = ComputeBackend.JAX
# ------------------------------------------------------------
# 数値精度
# ------------------------------------------------------------
# FP32を使用する。
# 最初のFP32： 主な計算の精度
# 2番目のFP32： 内部計算などで使用する精度
# CFDやLBMの学習用として扱いやすい設定。
precision_policy = PrecisionPolicy.FP32FP32
#ブロック【⑤ D2Q9を設定】の解説：２次元の標準的なLBMの速度モデル
# ============================================================
# ⑤ LBMの速度モデル D2Q9
# ============================================================
# D2Q9速度モデルを作成する。
# D2： 2次元 Q9： 9個の離散速度
# 中心： 0 周囲： 8方向 の合計9方向となる。
# なお、D2Q5 として流速が非常に遅い問題用もあるかXLBでは未実装
velocity_set = xlb.velocity_set.D2Q9(
precision_policy=precision_policy,
compute_backend=compute_backend,
)
#ブロック【⑥ 粘性係数・BGK係数を計算】の解説：寸法とレイノルズ数からこれらを計算
# ============================================================
# ⑥ 粘性係数とBGK緩和係数
# ============================================================
# Reynolds数の定義
# Re = U L / ν から、
# ν = U L / Re と変形して、動粘性係数viscを求める。
# ここで、代表速度：U = u_in, 代表長さ：L = square_size とする。
visc = u_in * square_size / Re
# ------------------------------------------------------------
# BGKの緩和係数
# ------------------------------------------------------------
# LBMでは衝突過程にBGKモデルを使用する。
# BGKモデルでは緩和係数omegaが必要になる。
# XLBで使用する格子単位系では、
# omega = 1 / (3ν + 0.5) として計算する。
# 緩和係数omegaと動粘性係数viscの関係式から、D2Q9の場合の関係式となる。
omega = 1.0 / (3.0 * visc + 0.5)
#ブロック【⑦ XLBを初期化】の解説：LBM計算をXLBで行う計算手法の基本設定
# ============================================================
# ⑦ XLBの初期化
# ============================================================
# XLBに、
# ・D2Q9: velocity_set
# ・JAX: compute_backend
# ・FP32: precision_policy
# を標準設定として登録する。
# 以降のXLBオブジェクトが、この設定を利用できるようになる。
xlb.init(
velocity_set=velocity_set,
default_backend=compute_backend,
default_precision_policy=precision_policy,
)
#ブロック【⑧ 計算格子を作成】の解説：設定した解析空間の計算格子を作成する
# ============================================================
# ⑧ 計算格子を作成
# ============================================================
# XLBのgrid_factoryを使用して、
# 300 × 120の2次元計算格子を作成する。
# イメージ：
# y
# ↑ 300 × 120
# ┌─────────────┐
# │ |
# │ |
# └─────────────┘ → x
grid = grid_factory(
grid_shape,
compute_backend=compute_backend,
)
#ブロック【⑨ 境界を取得】の解説：２次元計算領域の４辺の外周格子を取得する
# ============================================================
# ⑨ 計算領域の境界を取得
# ============================================================
# 計算領域の外周セルを取得する。
# bounding_box_indices()は、
# 格子の ・left ・right ・top ・bottom
# などの境界位置を取得するために使用する。
box = grid.bounding_box_indices()
# ------------------------------------------------------------
# 外周端を除いた境界
# ------------------------------------------------------------
# remove_edges=Trueにすると、
# 境界の角などの重複セルを除いた境界を取得する。
# 流入口・流出口の境界条件を設定するときに利用する。
box_no_edge = grid.bounding_box_indices(
remove_edges=True
)
#ブロック【⑩ 左側 → 流入口】の解説：解析空間の左側を流入部として定義する
# ============================================================
# ⑩ 流入口
# ============================================================
# 左側境界を流入口として使用する。
# 流体は → → → → →
# のように左から右へ流れる。
inlet = box_no_edge["left"]
#ブロック【⑪ 右側 → 流出口】の解説：解析空間の右側を流出部として定義する
# ============================================================
# ⑪ 流出口
# ============================================================
# 右側境界を流出口として使用する。
# 流体は右側から計算領域の外へ出ていく。
outlet = box_no_edge["right"]
#ブロック【⑫ 上下 → 固体壁】の解説：解析空間の上下において固定壁として処理する
# ============================================================
# ⑫ 上下の壁
# ============================================================
# 上側と下側の境界を固体壁として設定する。
# velocity_set.dは空間の次元数。
# D2Q9なのでd=2。
# x方向・y方向の境界情報を組み合わせる。
walls = [
box["bottom"][i] + box["top"][i]
for i in range(velocity_set.d)
]
# 同じセルが重複している場合に整理する。
walls = np.unique(
np.array(walls),
axis=-1
).tolist()
#ブロック【⑬ 格子座標を作成】の解説：設定した大きさの解析空間に格子座標を設定する
# ============================================================
# ⑬ 格子座標を作る
# ============================================================
# x方向の格子番号を作る。
# 0, 1, 2, ..., 299
x = np.arange(Nx)
# y方向の格子番号を作る。
# 0, 1, 2, ..., 119
y = np.arange(Ny)
# x-yの2次元座標を作る。
# X[i,j]： x方向の位置 Y[i,j]： y方向の位置
# これによって、すべての格子セルについて
# 「正方形の内部かどうか」を判定できる。
X, Y = np.meshgrid(
x,
y,
indexing="ij"
)
#ブロック【⑭ 正方形の角度を設定】の解説：角度で設定された回転をラジアン単位に変換する
# ============================================================
# ⑭ 正方形を回転させる
# ============================================================
# 角度を度からラジアンに変換する。
# NumPyのsin/cosはラジアンを使用するため、
# 入力したデグリーの角度を変換する必要がある。
theta = np.deg2rad(square_angle)
# ------------------------------------------------------------
# 正方形中心を原点とした座標に変換
# ------------------------------------------------------------
# すべての格子点について、
# 正方形の中心からの相対座標を計算する。
# Xc： x方向の中心からの距離 Yc： y方向の中心からの距離
Xc = X - square_cx
Yc = Y - square_cy
#ブロック【⑮ 回転座標に変換】の解説：回転した正方形を座標変換して準備する
# ============================================================
# ⑮ 回転座標系へ変換
# ============================================================
# 正方形をangleだけ回転した状態を、
# 数学的に判定するために座標変換を行う。
# 座領変換のために、回転行列を利用している。
# Xr = [ cosθ sinθ ][ Xc ]
# [ -sinθ cosθ ][ Yc ]
# この変換によって、回転した正方形を
# 「回転していない正方形」として判定できる。
Xr = (
np.cos(theta) * Xc
+ np.sin(theta) * Yc
)
# y方向についても同じ座標変換を行う。
Yr = (
-np.sin(theta) * Xc
+ np.cos(theta) * Yc
)
#ブロック【⑯ 正方形セルを判定】の解説：格子上での正方形に対応する位置を判定する
# ============================================================
# ⑯ 正方形内部のセルを判定
# ============================================================
# 正方形の半分の長さを計算する。
# square_size = 40 なので、
# half = 20 となる。
half = square_size / 2.0
# ------------------------------------------------------------
# 正方形内部の判定
# ------------------------------------------------------------
# 回転後の座標Xr,Yrについて、
# -20 <= Xr <= 20
# -20 <= Yr <= 20
# を満たすセルを正方形内部とする。
# つまり、回転座標系で見ると
# 普通の40×40の正方形になっている。
square_mask = (
(np.abs(Xr) <= half)
& (np.abs(Yr) <= half)
)
#ブロック【⑰ 固体セルのインデックスを作成】の解説：格子のなかで対象物体の部分の指標を設定する
# ============================================================
# ⑰ 正方形セルのインデックスを取得
# ============================================================
# square_maskがTrueになったセルの
# x,yインデックスを取得する。
# これがXLBに、「ここは固体です」
# と教えるための情報になる。
indices = np.where(square_mask)
# XLBが利用できる境界インデックス形式に変換する。
# 正方形のセルについて、
# x方向、y方向の位置をリストとして保存する。
square = [
indices[i].astype(np.int32).tolist()
for i in range(velocity_set.d)
]
#ブロック【⑱ 流入口 BC】の解説：流入部の一様流の速度境界条件を設定する
# ============================================================
# ⑱ 流入口の境界条件
# ============================================================
# 左側から速度を指定して流体を入れる。
# "velocity"： 速度境界条件
# prescribed_value=(u_in, 0.0) は、
# u_x = 0.05
# u_y = 0 を意味する。
bc_inlet = RegularizedBC(
"velocity",
prescribed_value=(u_in, 0.0),
indices=inlet,
)
#ブロック【⑲ 上下壁 BC】の解説：２次元解析空間の上下の境界条件を設定する
# ============================================================
# ⑲ 上下壁の境界条件
# ============================================================
# 上下の壁を固体壁として扱う。
# HalfwayBounceBackは、壁に向かって進んできた流体粒子を
# 壁で反射させる境界条件。
# イメージ
# 流体 →
# │壁
# │
# ← 反射 │
bc_walls = HalfwayBounceBackBC(
indices=walls
)
#ブロック【⑳ 正方形 BC】の解説：空間内の物体の対する境界条件を設定する
# ============================================================
# ⑳ 正方形物体の境界条件
# ============================================================
# 正方形も固体として扱う。
# 流体が正方形にぶつかった場合、HalfwayBounceBackによって反射する。
# その結果、以下のような流れが発生する。
# → → → →
# ■■■
# → → → ■■■
# ■■■
# → → →
# LBMで固体壁を表現するための境界条件で、
# 壁の「半分の位置」で粒子を跳ね返す方法です。
bc_square = HalfwayBounceBackBC(
indices=square
)
#ブロック【㉑ 流出部の境界条件】の解説：流出部において出口境界条件を設定する
# ============================================================
# ㉑ 流出口の境界条件
# ============================================================
# 右側から流体が計算領域の外へ出ていく出口境界条件を設定する。
# ExtrapolationOutflowBCは、出口付近の流れを外側へ外挿して
# 流体を計算領域から出すために使用する。
bc_outlet = ExtrapolationOutflowBC(
indices=outlet
)
#ブロック【㉒ 境界条件をまとめる】の解説：解析空間の４種類の境界条件をまとめる
# ============================================================
# ㉒ 境界条件をまとめる
# ============================================================
# ここまでで作った境界条件を、1つのリストにまとめる。
# これをStepperに渡すことで、XLBが各境界で適切な処理を行う。
boundary_conditions = [
bc_walls,
bc_inlet,
bc_outlet,
bc_square,
]
#ブロック【㉓ Navier-Stokes Stepper】の解説：CFDのナビエストークス流れのLBM対応
# ============================================================
# ㉓ Navier-Stokes Stepper
# ============================================================
# XLBの時間発展計算を担当するStepperを作成する。
# IncompressibleNavierStokesStepper：
# 非圧縮Navier-Stokes流れを
# LBMによって時間発展させる。
# collision_type="BGK"：
# 衝突モデルとしてBGKを使用する。
stepper = IncompressibleNavierStokesStepper(
grid=grid,
boundary_conditions=boundary_conditions,
collision_type="BGK",
)
#ブロック【㉔ LBM分布関数を初期化】の解説：粒子の挙動を定める分布関数の初期化
# ============================================================
# ㉔ 初期分布関数を作成
# ============================================================
# LBMでは、f_i # という「各方向へ移動する粒子の分布関数」を計算する。
# f_0： 現在の分布関数
# f_1： 次の計算で使用する分布関数
# bc_mask： 境界条件の情報
# missing_mask： 不足している隣接セルなどの情報
# を取得する。
f_0, f_1, bc_mask, missing_mask = (
stepper.prepare_fields()
)
#ブロック【㉕ 密度・速度計算】の解説：LBMで密度や速度などのマクロな物理量の計算
# ============================================================
# ㉕ 巨視的物理量を計算するオブジェクト
# ============================================================
# LBMでは、分布関数 f から、
# 密度 rho 速度 u
# などの通常の流体力学的な物理量を計算する。
# その計算を担当するのがMacroscopic。
macro = Macroscopic(
compute_backend=ComputeBackend.JAX,
precision_policy=precision_policy,
velocity_set=xlb.velocity_set.D2Q9(
precision_policy=precision_policy,
compute_backend=ComputeBackend.JAX,
),
)
#ブロック【㉖ 計算条件を表示】の解説：解析の全体設定条件をまとめて出力する
# ============================================================
# ㉖ 計算条件を画面に表示
# ============================================================
print()
print("======================================")
print(" XLB : Flow Around Rotated Square")
print("======================================")
# 計算格子を表示する。
print(f"Grid : {Nx} x {Ny}")
# 今回の計算バックエンドを表示する。
print(f"Backend : JAX-CPU")
# 流入速度を表示する。
print(f"Velocity : {u_in}")
# Reynolds数を表示する。
print(f"Reynolds : {Re}")
# 正方形のサイズを表示する。
print(f"Square size : {square_size} x {square_size}")
# 正方形の中心位置を表示する。
print(
f"Square center : "
f"({square_cx}, {square_cy})"
)
# 正方形の回転角度を表示する。
print(f"Square angle : {square_angle} deg")
# 動粘性係数を表示する。
print(f"Viscosity : {visc}")
# BGK緩和係数を表示する。
print(f"Omega : {omega}")
print("======================================")
print()
#ブロック【㉗ 結果保存関数】の解説：必要量を計算してVTK,PNG形式で保存する
# ============================================================
# ㉗ 結果を保存する関数
# ============================================================
# 一定時間ごとに呼び出して、
# その時点の流れ場を保存する関数。
# step： 現在の時間ステップ
# f： 現在のLBM分布関数
def post_process(step, f):
# --------------------------------------------------------
# 密度と速度を計算
# --------------------------------------------------------
# LBMの分布関数fから、
# rho： 密度 u： 速度 を計算する。
rho, u = macro(f)
# --------------------------------------------------------
# 外周1セルを除去
# --------------------------------------------------------
# 結果表示用に外周セルを除く。
# rhoは、[密度] のデータ。
rho = rho[:, 1:-1, 1:-1]
# uは2次元速度なので、
# u[0]： x方向速度 u[1]： y方向速度 となる。
u = u[:, 1:-1, 1:-1]
# --------------------------------------------------------
# 速度の大きさ
# --------------------------------------------------------
# 速度ベクトル u = (u_x, u_y) の大きさ
# |u| = sqrt(u_x^2 + u_y^2)
# を計算する。
u_magnitude = jnp.sqrt(
u[0] ** 2
+ u[1] ** 2
)
# ========================================================
# 保存する物理量
# ========================================================
# ParaViewなどで確認できるように、物理量を辞書にまとめる。
fields = {
# 密度
"rho": rho[0],
# x方向速度
"u_x": u[0],
# y方向速度
"u_y": u[1],
# 速度の大きさ
"u_magnitude": u_magnitude,
}
# ========================================================
# VTK形式で保存
# ========================================================
# VTKファイルとして保存する
# ParaViewを使用すれば、速度分布などをカラー表示できる。
save_fields_vtk(
fields,
timestep=step,
prefix="square_flow"
)
# ========================================================
# PNG画像として保存
# ========================================================
# 速度の大きさを画像として保存する。
# 解析結果を簡単に確認するときに便利。
save_image(
fields["u_magnitude"],
timestep=step,
prefix="square_flow"
)
# ========================================================
# 計算状況を表示
# ========================================================
# 現在のステップと最大速度を表示する。
print(
f"step = {step:6d} : "
f"max velocity = "
f"{float(jnp.max(u_magnitude)):.6f}"
)
#ブロック【㉘ 計算ステップ設定】の解説：全体の解析ステップ数と結果保存間隔
# ============================================================
# ㉘ 時間発展計算
# ============================================================
# 総計算ステップ数。
# LBMでは、この時間ステップを1回ずつ進めながら
# 流れを計算する。
num_steps = 5000
# ------------------------------------------------------------
# 結果を保存する間隔
# ------------------------------------------------------------
# 500ステップごとに結果を保存する。
# 5000ステップ計算するので、
# 0 → 500 → 1000 → ... → 5000 の結果が保存される。
post_process_interval = 500
#ブロック【㉙ タイマー開始】の解説：解析時間測定のため開始時間を記録する
# ============================================================
# ㉙ 計算時間測定開始
# ============================================================
# 計算開始時刻を記録する。
start_time = time.time()
#ブロック【㉚ LBM時間発展】の解説：指定ステップまで解析を進めて途中結果を保存する
# ============================================================
# ㉚ LBM時間発展ループ
# ============================================================
# 0から4999まで5000回、
# LBMの時間発展を行う。
for step in range(num_steps):
# --------------------------------------------------------
# LBMを1ステップ進める
# --------------------------------------------------------
# XLBのStepperに
# f_0： 現在の分布関数
# f_1： 作業用分布関数
# bc_mask： 境界条件情報
# missing_mask： 不足セル情報
# omega： BGK緩和係数
# step： 現在の時間ステップ
# を渡して計算する。
f_0, f_1 = stepper(
f_0,
f_1,
bc_mask,
missing_mask,
omega,
step,
)
# --------------------------------------------------------
# 分布関数のバッファを交換
# --------------------------------------------------------
# 今回計算したf_1を、次のステップのf_0として使用する。
# これにより、現在 → 次の時刻 の計算を効率的に進める。
f_0, f_1 = f_1, f_0
# ========================================================
# 結果保存
# ========================================================
# 500ステップごと、または最終ステップで結果を保存する。
if (
step % post_process_interval == 0
or step == num_steps - 1
):
# 現在の分布関数から、密度・速度・速度の大きさを計算して
# ファイルに保存する。
post_process(
step,
f_0
)
#ブロック【㉛ 計算終了】の解説：time関数の現時時点から開始時点を引いて経過時間
# ============================================================
# ㉛ 計算終了
# ============================================================
# 計算終了時刻を取得する。
elapsed = time.time() - start_time
#ブロック【㉜ 計算時間表示】の解説：解析の正常完了表示と解析時間の秒単位表示
# ============================================================
# ㉜ 計算時間を表示
# ============================================================
print()
print("======================================")
print(" Simulation finished")
print("======================================")
# 計算にかかった秒数。
print(f"Elapsed time : {elapsed:.3f} s")
print("======================================")
このスクリプト全体の機能と設定の１３項目について、対応するスクリプトと解析条件や解析結果と組み合わせて、下図にまとめます。上記の32ブロックの説明をあわせて、スクリプトの改良に活用してください。
!image.png
スクリプト「flow_square_2d-1.py」の実行
まず、以下の【01:JAXcpu導入方法と動作確認】で準備した解析環境を、以下の手順で準備します。XLBの作業用フォルダ「~/XLBjaxcpuWork」において、独自例題フォルダ「FLOW2D」を作っておきます。
qiita.com
PowerShell 7.6.6
PS C:\Windows\System32> wsl -l -v
NAME STATE VERSION
* Ubuntu-24.04 Stopped 2
Ubuntu-24.04-XLBjaxCPU Stopped 2
PS C:\Windows\System32> wsl -d Ubuntu-24.04-XLBjaxCPU
dalab@dalab2312nb01:/mnt/c/Windows/System32$ cd
dalab@dalab2312nb01:~$ ls
XLBjaxCPU XLBjaxcpuWork
dalab@dalab2312nb01:~$ source ~/XLBjaxCPU/bin/activate
(XLBjaxCPU) dalab@dalab2312nb01:~$ cd ~/XLBjaxcpuWork/
(XLBjaxCPU) dalab@dalab2312nb01:~/XLBjaxcpuWork$ ls
LDC2D XLB-git
(XLBjaxCPU) dalab@dalab2312nb01:~/XLBjaxcpuWork$ mkdir FLOW2D
(XLBjaxCPU) dalab@dalab2312nb01:~/XLBjaxcpuWork$ ls
FLOW2D LDC2D XLB-git
このなかで、「２次元問題として風洞中にある矩形周りの流れ」を、XLBで解析するために、以上で説明したスクリプト「flow_square_2d-1.py」を、これらの32ブロックを繋げて１つのファイルとして作ります。
(XLBjaxCPU) dalab@dalab2312nb01:~/XLBjaxcpuWork/TEST2D$ cd FLOW2D/
(XLBjaxCPU) dalab@dalab2312nb01:~/XLBjaxcpuWork/FLOW2D$ vi flow_square_2d-1.py
(XLBjaxCPU) dalab@dalab2312nb01:~/XLBjaxcpuWork/FLOW2D$ ls
flow_square_2d-1.py
これを、pythonで実行します。最初にGPUは未設定なので、CPUで計算を進めると表示され、４か所の境界条件の定義が登録されます。
そのあとに、独自例題の解析条件設定が表示されます。そのあとで、解析全体5000ステップで、500ステップ毎に途中結果の保存のログが表示されます。ここでは、500ステップでの解析時間と最大速度も表示されます。
ここで、流入口速度は$u_{\mathrm{in}}=0.05\quad [\mathrm{lu/ts}]$ と考えます。ここで lu は lattice unit（格子数）、ts は time step（時間ステップ）です。つまり、１時間ステップあたり、進む格子数で定義します。この単位で最大速度を見ると、最終5000ステップで0.134なので約2.7倍になっており、正方形の上下の狭い空間を、速く流れていることがわります。
最後に、解析時間が約10秒と表示されます。格子数300x120=36000で、全5000ステップを、Intel Corei9で実行した結果になります。
なお、Windwsのタスクマネージャーで確認すると、CPUで並列処理していることがわかります。
(XLBjaxCPU) dalab@dalab2312nb01:~/XLBjaxcpuWork/FLOW2D$ python flow_square_2d-1.py
No GPU support is available; CPU fallback will be used.
registered bc RegularizedBC_8332170304619 with id 1
registered bc HalfwayBounceBackBC_8332175202663 with id 2
WARNING! Assuming no-slip condition for BC type = HalfwayBounceBackBC!
registered bc HalfwayBounceBackBC_8332163779297 with id 3
WARNING! Assuming no-slip condition for BC type = HalfwayBounceBackBC!
registered bc ExtrapolationOutflowBC_8332163779354 with id 4
Warp 1.17.0 initialized:
CUDA Toolkit 12.9, Driver 13.2
Devices:
"cpu" : "x86_64"
"cuda:0" : "NVIDIA GeForce RTX 3070 Ti Laptop GPU" (8 GiB, sm_86, mempool enabled)
Kernel cache:
/home/dalab/.cache/warp/1.17.0
======================================
XLB : Flow Around Rotated Square
======================================
Grid : 300 x 120
Backend : JAX-CPU
Velocity : 0.05
Reynolds : 100.0
Square size : 40 x 40
Square center : (75, 60)
Square angle : 30.0 deg
Viscosity : 0.02
Omega : 1.7857142857142856
======================================
Saved ./square_flow_0000000.vtk in 0.020215 seconds.
step = 0 : max velocity = 0.000000
Saved ./square_flow_0000500.vtk in 0.002383 seconds.
step = 500 : max velocity = 0.156543
Saved ./square_flow_0001000.vtk in 0.002592 seconds.
step = 1000 : max velocity = 0.129826
Saved ./square_flow_0001500.vtk in 0.002317 seconds.
step = 1500 : max velocity = 0.130439
Saved ./square_flow_0002000.vtk in 0.002421 seconds.
step = 2000 : max velocity = 0.133740
Saved ./square_flow_0002500.vtk in 0.002411 seconds.
step = 2500 : max velocity = 0.135603
Saved ./square_flow_0003000.vtk in 0.002438 seconds.
step = 3000 : max velocity = 0.132879
Saved ./square_flow_0003500.vtk in 0.002409 seconds.
step = 3500 : max velocity = 0.132424
Saved ./square_flow_0004000.vtk in 0.007832 seconds.
step = 4000 : max velocity = 0.138548
Saved ./square_flow_0004500.vtk in 0.007650 seconds.
step = 4500 : max velocity = 0.137304
Saved ./square_flow_0004999.vtk in 0.002434 seconds.
step = 4999 : max velocity = 0.133921
======================================
Simulation finished
======================================
Elapsed time : 10.336 s
======================================
解析結果として、500ステップ毎の「可視化用vtkファイル」と「速度分布pngファイル」が、保存されています。
(XLBjaxCPU) dalab@dalab2312nb01:~/XLBjaxcpuWork/FLOW2D$ ls
flow_square_2d-1.py square_flow_0002000.vtk square_flow_0004999.vtk square_flow_3000.png
square_flow_0000.png square_flow_0002500.vtk square_flow_0500.png square_flow_3500.png
square_flow_0000000.vtk square_flow_0003000.vtk square_flow_1000.png square_flow_4000.png
square_flow_0000500.vtk square_flow_0003500.vtk square_flow_1500.png square_flow_4500.png
square_flow_0001000.vtk square_flow_0004000.vtk square_flow_2000.png square_flow_4999.png
square_flow_0001500.vtk square_flow_0004500.vtk square_flow_2500.png
速度分布のPNG画像ファイルを、ステップ0500-4999の10枚分を確認します。前半で流れが成長し、後半でカルマン渦が確認できます。これは、正方形が30度回転しており、上下の流れが非対称になっており、カルマン渦が誘導されたと思います。
!image.png
次に、可視化用のVTKデータファイル「square_flow_...vtk」を、ParaViewを用いて確認します。
!image.png
この可視化用データでは、密度rhoを表示することもできます。
!image.png
さらに、MP4形式の動画ファイルを作成することもできます。ParaViewの手順（save animation）を確認して試してください。
以上で説明したスクリプト「flow_square_2d-1.py」において、変更可能なパラメータを、下図にまとめます。まず最初は、スクリプトの設定通りで解析を実行して、動作確認をします。そのあとは、この図に示すいろいろなパラメータを変更して、比較検討してみてください。
!image.png
例えば、正方形の回転角度０度とすると、傾かない状態の速度分布は、以下のようになります。流れは成長しますが、上下対称なので、カルマン渦は発生していません。
!image.png
コンパクト版スクリプト「flow_square_2d-0.py」
スクリプト全体の把握し、機能拡張などを目的として、コメントなどを削除して、ログの画面表示も圧縮して、コンパクトにしたスクリプト「flow_square_2d-0.py」を以下に示します。こちらを用いて一括コピーしても、同じ解析結果が得られます。なお、対応するブロック①～㉜について、該当する行末に記してあります。
import os
os.environ["JAX_PLATFORMS"] = "cpu"
import time
import numpy as np
import jax.numpy as jnp
import xlb
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, RegularizedBC, ExtrapolationOutflowBC
from xlb.operator.macroscopic import Macroscopic
from xlb.utils import save_fields_vtk, save_image #①
Nx, Ny = 300, 120
u_in, Re = 0.05, 100.0 #②
square_size, square_cx, square_cy, square_angle = 40, Nx // 4, Ny // 2, 30.0 #③
compute_backend = ComputeBackend.JAX
precision_policy = PrecisionPolicy.FP32FP32 #④
velocity_set = xlb.velocity_set.D2Q9(precision_policy=precision_policy, compute_backend=compute_backend) #⑤
visc = u_in * square_size / Re
omega = 1.0 / (3.0 * visc + 0.5) #⑥
xlb.init(velocity_set=velocity_set, default_backend=compute_backend, default_precision_policy=precision_policy) #⑦
grid = grid_factory((Nx, Ny), compute_backend=compute_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"][i] + box["top"][i] for i in range(velocity_set.d)]
walls = np.unique(np.array(walls), axis=-1).tolist() #⑫
X, Y = np.meshgrid(np.arange(Nx), np.arange(Ny), indexing="ij") #⑬
theta = np.deg2rad(square_angle)
Xc, Yc = X - square_cx, Y - square_cy
Xr = np.cos(theta) * Xc + np.sin(theta) * Yc
Yr = -np.sin(theta) * Xc + np.cos(theta) * Yc #⑭
half = square_size / 2.0 #⑮
square_mask = (np.abs(Xr) <= half) & (np.abs(Yr) <= half) #⑯
indices = np.where(square_mask)
square = [indices[i].astype(np.int32).tolist() for i in range(velocity_set.d)] #⑰
bc_inlet = RegularizedBC("velocity", prescribed_value=(u_in, 0.0), indices=inlet) #⑱
bc_walls = HalfwayBounceBackBC(indices=walls) #⑲
bc_square = HalfwayBounceBackBC(indices=square) #⑳
bc_outlet = ExtrapolationOutflowBC(indices=outlet) #㉑
boundary_conditions = [bc_walls, bc_inlet, bc_outlet, bc_square] #㉒
stepper = IncompressibleNavierStokesStepper(grid=grid, boundary_conditions=boundary_conditions, collision_type="BGK") #㉓
f_0, f_1, bc_mask, missing_mask = stepper.prepare_fields() #㉔
macro = Macroscopic(compute_backend=ComputeBackend.JAX, precision_policy=precision_policy,
velocity_set=xlb.velocity_set.D2Q9(precision_policy=precision_policy, compute_backend=ComputeBackend.JAX)) #㉕
print(f"Grid={Nx}x{Ny}, Backend=JAX-CPU, U={u_in}, Re={Re}, Square={square_size}x{square_size}, Center=({square_cx},{square_cy}), Angle={square_angle}°, ν={visc}, ω={omega}") #㉖
def post_process(step, f): #㉗
rho, u = macro(f)
rho, u = rho[:, 1:-1, 1:-1], u[:, 1:-1, 1:-1]
u_magnitude = jnp.sqrt(u[0] ** 2 + u[1] ** 2)
fields = {"rho": rho[0], "u_x": u[0], "u_y": u[1], "u_magnitude": u_magnitude}
save_fields_vtk(fields, timestep=step, prefix="square_flow")
save_image(fields["u_magnitude"], timestep=step, prefix="square_flow")
print(f"step={step:6d}, max velocity={float(jnp.max(u_magnitude)):.6f}")
num_steps, post_process_interval = 5000, 500 #㉘
start_time = time.time() #㉙
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 % post_process_interval == 0 or step == num_steps - 1:
post_process(step, f_0)
elapsed = time.time() - start_time #㉛
print(f"Simulation finished: {elapsed:.3f} s") #㉜
LBMの考え方を改めて確認する
以上で、「２次元問題として風洞中にある矩形周りの流れ」 を、XLBのスクリプトを用いて、LBM解析で求めることが出来ました。ここで、LBMを初めて学ぶ場合、まず 「流体そのものを直接動かすのではなく、各格子点にある“粒子の分布”を動かして、結果として流体の速度や密度を得る」 として考えると理解しやすいです。
今回のXLBの例題を理解するために、D2Q9 → 分布関数 → Collision → Streaming → 速度・密度の順に、できるだけ具体的に説明します。
1. LBMの基本的な考え方
通常のCFDでは、Navier-Stokes方程式を直接解いて、$u_x,\quad u_y,\quad p$ などを求めます。
一方、LBMでは少し違います。以下の手順、つまり、「分布関数を計算することで、流体の速度・密度を間接的に求める」 という方法です。
LBM(XLB) → 分布関数 f_i を計算 → Collision（衝突）
→ Streaming（移動） → 分布関数 f_i の更新 → 密度 ρ、速度 u を計算
2. 「分布関数」とは何か？
ここがLBMで最も重要な部分です。ある格子点を考えると、「D2Q9」 では、1つの格子点に9方向の分布関数があります。
!image.png
一般に、$f_0,f_1,f_2,\cdots,f_8$ と表します。例えば、以下のように、それぞれの方向に対応しています。
f0 → その場にとどまる成分
f1 → 右 f2 → 上 f3 → 左 f4 → 下
f5 → 右上 f6 → 左上 f7 → 左下 f8 → 右下
3. 分布関数は「粒子の量」と考える
例えば、ある格子点で、以下の状態だったとします。
右方向：f1 = 0.20・左方向：f3 = 0.10・上方向：f2 = 0.12・下方向：f4 = 0.12
右向きの分布が左向きより大きいので、「右 →→」という流れが生じます。つまり、各方向にどれだけ“流体の分布”が存在しているか を表しているのが f_i です。ここで重要なのは、f_i は普通の意味での「実際の粒子1個1個」ではないことです。LBMでは、流体の状態を9方向の分布として表現していると考える方が適切です。
4. 分布関数から密度を求める
9個の分布関数を全部足すと、密度になります。$\rho=\sum_i f_i$ つまり、「f0 + f1 +...+ f8：全部足すと、密度ρ」になります。
5. 分布関数から速度を求める
速度は、各方向の分布に「方向」を掛けて足し合わせます。$\rho\mathbf{u}=\sum_i f_i\mathbf{c}_i$ したがって、$\mathbf{u}=\frac{1}{\rho}\sum_i f_i\mathbf{c}_i$ となります。ここで、$\mathbf{c}_i$ は各方向の速度です。例えば、「右方向：+x 左方向：-x 上方向：+y 下方向：-y」です。
6. では「Collision」とは何か？
ここからLBMの核心です。英語のCollisionは、日本語では「衝突」です。ただし、流体解析のLBMでは、 分布関数を局所的に平衡状態へ近づける計算 と考える方が分かりやすいです。
7. Collisionのイメージ
例えば、ある格子点で、方向別の値が、以下のような状態だったとします。流体としては少し「乱れた」分布になっています。
現在の分布
右 ████████ 左 ██ 上 ████ 下 ███ 斜め ███...
そこで、「現在の分布 → Collision → 平衡状態に近づける」という操作をします。
8. BGKモデル
今回のXLBのスクリプトでは、「collision_type="BGK"」 を使用しています。BGKでは、非常に簡単に書くと、$f_i^{*}=f_i+\omega(f_i^{eq}-f_i)$ です。
ここで、$f_i$：現在の分布関数、$f_i^{eq}$：平衡分布関数、$f_i^*$：Collision後の分布関数、$\omega$：緩和係数、です。
9. この式を日本語にすると
式 $f_i^{*}=f_i+\omega(f_i^{eq}-f_i)$ を見ると、$f_i^{eq}-f_i$ は、「現在の分布と平衡状態との差」です。したがって、以下のような流れになります。
現在の分布 → 平衡分布との差 → 現在の分布を修正 → Collision後の分布
[差を計算] [ωを掛ける]
10.緩和係数 ωは何を意味する？
今回のスクリプトでは、$\omega=\frac{1}{3\nu+0.5}$ としていました。例えば、Re=100 の場合、ν=0.02 なので、$\omega=\frac{1}{3(0.02)+0.5}\approx1.4286$ となります。つまり、「Re → 粘性 ν
→ 緩和係数 ω → Collisionの強さ」という関係になっています。
11. Collisionだけでは流体は流れない
ここが初心者には特に重要です。Collisionをすると、「格子点A ● → 分布が変化します。」しかし、それだけでは隣の格子点には移動していません。そこで次に、Streaming を行います。
12. Streamingとは？
Streamingは、分布関数を、その分布が向いている隣の格子点へ移動させる 処理です。例えば、以下のように、左図の３つの格子点があるとします。A点に存在する右向き分布 $f_1$ は、Streamingによって、右図のように移動します。
A B C [Streaming] A f1 B
● → ● → ● ● → ●
つまり、$f_i(\mathbf{x}+\mathbf{c}_i,t+1)=f_i^*(\mathbf{x},t)$ です。
13. これが「流れる」の正体
例えば右向きの分布だけを見ると、それぞれの格子点にある右向き分布が、以下のように隣へ移動します。これを大量の格子点で繰り返します。すると、クロな流れが現れます。
時刻 t 時刻 t+1
● →1 ● →2 ● →3 ● ● →0 ● →1 ● →2 ●
14. CollisionとStreamingの違い
ここを明確に区別するとLBMがかなり分かりやすくなります。
| 処理 | 何をする？ | イメージ |
|---|---|---|
| Collision | 分布を平衡状態に近づける | その場で分布を混ぜる |
| Streaming | 分布を隣の格子へ移動する | 分布を運ぶ |
15. LBMの1ステップ
LBMの基本的な1ステップを非常に単純化すると、以下のような流れになります。
!image.png
16. 今回の正方形まわり流れで考える
今回のモデルは、２次元の長方形の流体空間の内部に、正方形の物体があります。左から流体が入ってきて、LBMでは、各格子点で以下の処理を繰り返します。
分布関数 f_i → Collision → Streaming → 境界条件 → 次の時刻の f_i
17. 正方形にぶつかるとどうなる？
Streamingをすると、分布が正方形の内部へ入ろうとする場合があります。しかし、正方形は固体なので、そのまま通過させません。そこで、Bounce Back などの境界条件を使います。
イメージとしては、流体の格子からの分布を反射させます。今回のスクリプトでは、「HalfwayBounceBackBC」を使っています。
18. そして後流ができる
正方形の後ろでは、分布関数の移動とCollisionを何度も繰り返すことで、複雑な流れが形成されます。これをマクロに見ると、後流・再循環・渦 として現れます。
19. 最後に「速度」を取り出す
LBMが直接計算している中心的な量は分布関数 $f_i$ です。そこから、$\rho=\sum_i f_i$と、$\mathbf{u}=\frac{1}{\rho}\sum_i f_i\mathbf{c}_i$ を計算します。今回のXLBスクリプトの、「rho, u = macro(f)」は、まさにこの処理です。
!image.png
20. XLBのスクリプトとの対応
今回のスクリプトをLBMの考え方に対応させると、次のような、説明図になります。
!image.png
21. f_0 と f_1 は何？
今回のプログラムには、「f_0, f_1 = stepper.prepare_fields()」がありました。ここは少し注意が必要です。f_0とf_1はD2Q9の9方向そのものではありません。 D2Q9の方向は、$f_0,\ldots,f_8$ ですが、XLBのスクリプトに出てくる「f_0 f_1」は、時間発展の前後の分布関数を格納するための2つの配列（バッファ） として使われています。
そして、以下の処理によって、計算した結果を次のステップへ渡しています。ここはLBMを初めて学ぶと非常に混乱しやすいところです。
f_0, f_1 = stepper(
f_0, f_1,
bc_mask,
missing_mask,
omega,
step
)
f_0, f_1 = f_1, f_0
22. 一番簡単な覚え方
LBMはまず、次の3語を覚えるとよいです。これらを何千回も繰り返すと、最終的に 「流体の流れ」 が現れます。
【分布関数：「各方向にどれだけ流体が分布しているか」】
↓【Collision：「その場所で分布を平衡状態に近づける」】
↓【Streaming：「その分布を方向に従って隣の格子へ運ぶ」】
以上で説明してきた、以下の流れで理解すると、XLBのプログラムが何をしているのかがかなり見通しやすくなります。
分布関数 → Collision → Streaming → 境界条件 → 次の時刻
→ .... → 正方形まわりの流れ → ux, uy, |u| → 解析結果
これを説明図としてまとめると、以下のようになります。
!image.png
2 likes
Mentions
LBM流体解析ツールXLBの活用【10:XLB独自例題の作成と分析】 この文書では、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計算を比較的容易に実装できます。 このXLBのgithubリポジトリは、以下にあります。 Autodesk/XLB ここでは、最も基本となるCPUで動作する解析環境「JAX-CPU」を、以下の手順で準備した状態を用いて、独自例題の作成と分析を行います。 qiita.com 対象とする独自例題の概要 ここでは、「２次元問題として風洞中にある矩形周りの流れ」 を分析することを、目標とします。例えば、以下の論文では、格子ボルツマン法と有限体積法を用いた解析が報告されており、これを参考にします。 「Accurate computations of the laminar flow past a square cylinder based on two different methods: lattice-Boltzmann and finite-volume」 この論文中の図２では、以下の結果が示されています。 !image.png Fig. 2. Streamlines around the square cylinder for different Reynolds numbers (a) Re = 1; (b) Re = 30; (c) Re = 60; (d) Re = 200 解析モデルと解析条件と境界条件などの設定 これより、レイノルズ数Reの違いによって、正方形周りの流れの変化を、比較検討できます。解析モデルと境界条件は、以下の図のような定義です。 !image.png 独自例題の、格子ボルツマン法LBMをAutodesk XLBで実施するために関する設定は、以下の通りです。 | 項目 | 設定 | 項目 | 設定 | |--|--|--|--| | 計算領域 | 2次元長方形| 流体モデル | D2Q9 | | 衝突モデル | BGK| 流入 | 一定速度 | | 流出 | ExtrapolationOutflowBC| 上下壁 | HalfwayBounceBackBC | | 正方形物体 | HalfwayBounceBackBC| 計算バックエンド | JAX | | GPU(Jax,Warp,Neon) | 使用しない| CPU | JAX-CPU | 独自例題の具体的な設計 以下では、独自例題のスクリプトを作るうえで、設定する項目をまとめておきます。 スクリプト内の変数名と上記設定値との対応 独自例題での設定値とスクリプト内の 「変数名との対応」 を、下表に示します。問題を簡単にするため、格子や物体などの寸法は、無次元化した数値を用いています。 | 項目 | 説明図の内容 | スクリプトの変数 | 設定値 | 単位・意味 | |---|---|---|---:|---| | 計算領域・X方向 | 領域の長さ | Nx | 300 | 格子点 | | 計算領域・Y方向 | 領域の高さ | Ny | 120 | 格子点 | | 流入速度 | X方向流入速度 | u_in | 0.05 | 格子単位/step | | レイノルズ数 | 流れの無次元数 | Re | 100.0 | 無次元 | | 正方形一辺 | 物体の大きさ | square_size | 40 | 格子点 | | 正方形中心 X | 物体中心のX座標 | square_cx | 75 | 格子座標 | | 正方形中心 Y | 物体中心のY座標 | square_cy | 60 | 格子座標 | | 正方形回転角 | 物体の傾き | square_angle | 30.0 | 度 | | 動粘性係数 | LBMの粘性 | visc | 0.020 | 格子単位 | | 緩和パラメータ | BGKの緩和係数 | omega | 約1.7857 | 無次元 | | 計算ステップ数 | 時間発展回数 | num_steps | 5000 | step | | 出力間隔 | 結果を保存する間隔 | post_process_interval | 500 | step | さらに重要なのは、「計算で決まる数値」 です。例えば、動粘性係数：visc と BGK緩和パラメータ：omega は直接入力しているのではなく、上の数値から計算されています。 visc = u_in * square_size / Re omega = 1.0 / (3.0 * visc + 0.5) したがって、$\nu=\frac{0.05\times40}{100}=0.020$、そして、$\omega=\frac{1}{3\times0.020+0.5} \approx1.7857$です。 つまり、上の方のユーザーが変更する基本パラメータ*（設定値）と、そこから自動的に計算されるパラメータ（計算値）を分けると理解しやすくなります。この表にすると、「どの数値を変更すると、解析モデルの何が変わるのか」がかなり分かりやすくなります。 | 種類 | 変数 | 値 | 役割 | |---|---|---:|---| | 計算値 | visc | 0.020 | 動粘性係数 | | 計算値 | omega | 1.7857 | BGK緩和パラメータ | 以上の設定により、下図の８の解析結果のイメージを目標とします。 !image.png 独自例題を理解する上でのポイント この独自例題は、以下の【04:例題flow_past_sphere_3d.pyの分析】に対して、２次元に簡略化した内容です。ただし、流体内部の四角形の「位置・規模・角度」が任意で設定できるので、様々な条件を確認することができます。これを参考にして、「左から流入 → 中央の正方形物体 → 右へ流出」 という、外部流れの最小例にできます。 qiita.com 今回作るモデル まずは非常に単純化します。設定する流体空間中に、正方形を任意の位置・規模・角度に、境界条件：HalfwayBounceBackとして置くことができます。２次元XY座標において、流体はX方向の＋向きに流れ、左辺が流入（Velocity BC）で右辺が流出（Outflow BC）です。 このプログラムの重要なところ 今回の目的は、XLBの構造をできるだけ単純に理解することです。特に重要なのは次の部分を解説します。 【D2Q9】 2次元LBMの最も基本的な速度モデルです。スクリプトでは以下の定義です。中央の1方向＋周囲8方向なので D2Q9 です。 velocity_set = xlb.velocity_set.D2Q9( precision_policy=precision_policy, compute_backend=compute_backend, ) 【JAX-CPU】 今回は計算規模も小さく、どの計算環境でも動作させるため、CPU動作 としています。従って、「Python → XLB → JAX → CPU」のような計算の仕組みになります。 compute_backend = ComputeBackend.JAX os.environ["JAX_PLATFORMS"] = "cpu" 【正方形】 今回の解析対象の扱いの中心部分です。回転させた正方形内部を座標返還して、正方形内部の格子を判定しています。詳しくはスクリプトを見ながら、説明します。 【正方形の境界条件】 ここがLBMでは重要です。スクリプトでは以下の定義です。以下の設定で、正方形に到達した流体粒子を反射させる というモデルです。XLBの公式2次元キャビティ例でも、固体壁には HalfwayBounceBackBC が使われています。 bc_square = HalfwayBounceBackBC( indices=square ) 【流入口】 ここでは、流入として「X方向にu_in=0.05・Y方向に0.0」としています。つまり一様流を左辺から入れます。 bc_inlet = RegularizedBC( "velocity", prescribed_value=(u_in, 0.0), indices=inlet, ) 【流出口】 ここでは、流入として、以下のような定義とします。右側では流体を外へ出します。 bc_outlet = ExtrapolationOutflowBC( indices=outlet ) 【斜めの境界条件】 この例題では、正方形の中心位置はそのままにして、指定した角度だけ回転させるように変更できます。ここでは、正方形そのものを回転させるのではなく、格子上の各セルについて「そのセルが回転した正方形の内部か？」を判定しています。 したがってXLBでは、「数学的な正方形 → 回転座標系で判定 → square_mask → 固体セル → HalfwayBounceBackBC」という流れになります。 なお、LBMの格子は整数セルなので、30°や45°にすると境界が階段状（staircase） になります。これは今回のような基本的なXLB学習例としては重要なポイントです。LBMにおける斜めの境界条件は、階段状のギザギザ（ボクセル表現）による誤差を減らし、滑らかな傾斜や曲面を表現するために非常に重要です。 【Reynolds数】 今回は、「Re = 100.0」としています。代表長さを正方形の一辺「square_size = 20」として、$Re=\frac{UL}{\nu}$から、動粘性係数 visc = u_in * square_size / Re としています。 さらにBGKの緩和係数 $\omega=\frac{1}{3\nu+1/2}$ を、「omega = 1.0 / (3.0 * visc + 0.5)」で計算しています。 なお、正方形の寸法Lを変化させた場合は、本来ならば、レイノルズ数Reが変化するところですが、このスクリプトでは、動粘性係数viscを計算しなおすので、レイノルズ数は変化しない設定となっています。 【公式例題との関係】 特に、現在のXLB公式例では3次元の球まわり流れについて、HalfwayBounceBackBC、RegularizedBC、ExtrapolationOutflowBC、IncompressibleNavierStokesStepper を組み合わせる構成が使われています。今回のスクリプトは、それを2次元・正方形に最小化したものと考えると分かりやすいです。GitHub なお、現在のXLB main版の公式2D例では save_fields_vtk() と save_image() による結果保存も使われていますので、今回も同じ方式にしています。GitHub 独自例題のスクリプトの解説 以下では、スクリプトの記述に従って、ブロック毎に区分して、それぞれの機能をコメントの形で記入したので、これを用いて説明します。 以下の説明では、ブロックごとに分けて解説しますが、全体をつなげると「２次元空間での正方形周囲の流れ」を、XLB用いてLBM計算するスクリプト「flow_square_2d-1.py」としておきます。 （最後の「-1」は、最後に説明するスクリプトの改良を期待してのバージョン１と考えています） スクリプト「flow_square_2d-1.py」の構成 このスクリプトは、以下の 「32のブロック」 に、機能が分かれています。各ブロックについて、スクリプトを提示して、コメント部分に説明を追加します。 ① JAXをCPUに設定 → ② 計算条件を設定 → ③ 正方形のサイズ・位置・角度を設定 → ④ XLBのバックエンドをJAXに設定 → ⑤ D2Q9を設定 → ⑥ 粘性係数・BGK係数を計算 → ⑦ XLBを初期化 → ⑧ 計算格子を作成 → ⑨ 境界を取得 → ⑩ 左側 → 流入口 → ⑪ 右側 → 流出口 → ⑫ 上下 → 固体壁 → ⑬ 格子座標を作成 → ⑭ 正方形の角度を設定 → ⑯ 正方形セルを判定 → ⑰ 固体セルのインデックスを作成 → ⑱ 流入口 BC → ⑲ 上下壁 BC → ⑳ 正方形 BC → ㉑ 流出口 BC → ㉒ 境界条件を統合 → ㉓ Navier-Stokes Stepper → ㉔ LBM分布関数を初期化 → ㉕ 密度・速度計算 → ㉖ 計算条件を表示 → ㉗ 結果保存関数 → ㉘ 計算ステップ設定 → ㉙ タイマー開始 → ㉚ LBM時間発展 → ㉛ 計算終了 → ㉜ 計算時間表示 #ブロック【① JAXをCPUに設定】の解説：CPUによるXLBによるLBM計算設定 # ============================================================ # XLB : 2D Flow Around a Rotated Square : flow_square_2d-1.py # # 2次元の長方形領域の中に、正方形の物体を置き、 # 左から右へ流体を流す外部流れ問題を計算する。 # この正方形は、任意の位置と寸法と角度で置ける。 # 計算方法： LBM (Lattice Boltzmann Method) # 速度モデル： D2Q9 # 衝突モデル： BGK # バックエンド： JAX # 計算デバイス： CPU # 重要： このプログラムでは GPU(Jax,Warp,Neon) を使用しない。 # ============================================================ # ============================================================ # ① JAXをCPUで使用する設定 # ============================================================ # JAXが利用する計算デバイスをCPUに限定する。 # "cpu"を指定すると、 # Python → XLB → JAX → CPU という計算経路になる。 # NVIDIA GPUやCUDAは、この例では使用しない。 # この環境変数はJAXをimportする前に設定する必要がある。 import os # 念の為、環境変数JAX_PLATFORMSをcpuに設定する。 os.environ["JAX_PLATFORMS"] = "cpu" # ============================================================ # Pythonで使用する基本ライブラリ # ============================================================ # 時間計測に使用する。 # 解析全体に何秒かかったかを最後に表示するために使用する。 import time # NumPyは、 # ・格子座標の作成 ・正方形の領域判定 ・角度計算 # などに使用する。 import numpy as np # JAX版NumPy。 # 計算結果の速度や密度など、 # XLB内部でJAX配列として扱われるデータを処理するときに使用する。 import jax.numpy as jnp # XLB本体を読み込む。 import xlb # XLBが使用する計算バックエンドを指定するためのクラス。 # 今回は ComputeBackend.JAX を使用する。 from xlb.compute_backend import ComputeBackend # 数値精度を指定するためのクラス。 # 今回はFP32を使用する。 from xlb.precision_policy import PrecisionPolicy # 計算格子を作成するための関数。 from xlb.grid import grid_factory # CFDに対応するNavier-Stokes方程式をLBMとして時間発展させる # XLBの主要なStepper。これでステップ毎の処理を進める。 from xlb.operator.stepper import ( IncompressibleNavierStokesStepper ) # 境界条件を指定するクラスを読み込む。 # HalfwayBounceBackBC： # 固体壁で流体を反射させる。 # RegularizedBC： # 流入口などで速度を指定する。 # ExtrapolationOutflowBC： # 流出口から流体を外へ流す。 from xlb.operator.boundary_condition import ( HalfwayBounceBackBC, RegularizedBC, ExtrapolationOutflowBC, ) # 密度や速度などの巨視的物理量を # 分布関数から計算するクラス。 from xlb.operator.macroscopic import Macroscopic # 計算結果をVTKや画像として保存する関数。 # VTK： ParaViewなどで解析結果を確認できる。 # PNG： 速度分布などを画像として確認できる。 from xlb.utils import ( save_fields_vtk, save_image, ) #ブロック【② 計算条件を設定】の解説：解析空間や流入速度やレイノルズ数の設定 # ============================================================ # ② 計算条件 # ============================================================ # x方向の格子数。今回は300セル。 # 流れは基本的にx方向、つまり「左 → 右」に進む。 # 流入から物体、特に物体から流出においては、十分な空間を確保する。 Nx = 300 # y方向の格子数。今回は120セル。 # 対象の物体の上下に、十分な空間を確保する。 Ny = 120 # XLBに渡す計算領域のサイズ。 # (Nx, Ny) つまり x方向 = 300・y方向 = 120 の2次元格子になる。 # 対象とする空間内の物体の大きさに対応して、計算領域を定める grid_shape = (Nx, Ny) # ------------------------------------------------------------ # 流入速度 # ------------------------------------------------------------ # 左側から流入する流体の速度。 # 今回は無次元速度として0.05を使用する。 # x方向： 0.05 y方向： 0 # したがって流体は水平に流入する。 u_in = 0.05 # ------------------------------------------------------------ # Reynolds数 # ------------------------------------------------------------ # 流れの状態を表す無次元数。 # Re = U L / ν # U： 代表速度 L： 代表長さ ν： 動粘性係数 # 今回は正方形の一辺を代表長さとする。例題ではReを前提条件とする。 Re = 100.0 #ブロック【③ 正方形のサイズ・位置・角度を設定】の解説：中心位置と回転角度の設定 # ============================================================ # ③ 正方形物体の条件 # ============================================================ # 正方形の一辺の長さ。 # 40 × 40セルの正方形を作る。流れに直交方向の格子数Nyと比較する # 流体空間の格子数 Nx,Ny に適切に入る大きさが限界となる。 # 閉塞率=物体の断面積／空間の断面積が、50%以下になるように注意する。 square_size = 40 # ------------------------------------------------------------ # 正方形の中心位置 # ------------------------------------------------------------ # x方向の中心。 # Nx = 300なので、300 / 4 = 75 となる。 # したがって正方形の中心はx=75。流体空間の流入側の 1/4 に位置する。 # なお // は整数除算なので、割切れないときは、小数点切り捨ての整数化 square_cx = Nx // 4 # y方向の中心。 # Ny = 120なので、120 / 2 = 60 となる。 # したがって正方形の中心はy=60。上下方向の中央に位置する。 square_cy = Ny // 2 # ------------------------------------------------------------ # 正方形の回転角度 # ------------------------------------------------------------ # 正方形を反時計回りに回転させる角度。 # 0度： 通常の正方形 # 30度： 30度回転 # 45度： ひし形に近い形 # 90度： 見た目は0度と同じ # この値を変更することで、 # 同じ計算プログラムで様々な角度を試すことができる。 # この角度は、細かな数値の任意の数値で設定可能です。 square_angle = 30.0 #ブロック【④ XLBのバックエンドをJAXに設定】の解説：併せて数値制度をFP32に設定 # ============================================================ # ④ XLBの計算バックエンド # ============================================================ # XLBの計算バックエンドとしてJAXを選択する。 # JAXにはCPU/GPUなどの実行環境があるが、 # 先ほど、os.environ["JAX_PLATFORMS"] = "cpu" # と指定しているため、今回はCPUで計算される。 compute_backend = ComputeBackend.JAX # ------------------------------------------------------------ # 数値精度 # ------------------------------------------------------------ # FP32を使用する。 # 最初のFP32： 主な計算の精度 # 2番目のFP32： 内部計算などで使用する精度 # CFDやLBMの学習用として扱いやすい設定。 precision_policy = PrecisionPolicy.FP32FP32 #ブロック【⑤ D2Q9を設定】の解説：２次元の標準的なLBMの速度モデル # ============================================================ # ⑤ LBMの速度モデル D2Q9 # ============================================================ # D2Q9速度モデルを作成する。 # D2： 2次元 Q9： 9個の離散速度 # 中心： 0 周囲： 8方向 の合計9方向となる。 # なお、D2Q5 として流速が非常に遅い問題用もあるかXLBでは未実装 velocity_set = xlb.velocity_set.D2Q9( precision_policy=precision_policy, compute_backend=compute_backend, ) #ブロック【⑥ 粘性係数・BGK係数を計算】の解説：寸法とレイノルズ数からこれらを計算 # ============================================================ # ⑥ 粘性係数とBGK緩和係数 # ============================================================ # Reynolds数の定義 # Re = U L / ν から、 # ν = U L / Re と変形して、動粘性係数viscを求める。 # ここで、代表速度：U = u_in, 代表長さ：L = square_size とする。 visc = u_in * square_size / Re # ------------------------------------------------------------ # BGKの緩和係数 # ------------------------------------------------------------ # LBMでは衝突過程にBGKモデルを使用する。 # BGKモデルでは緩和係数omegaが必要になる。 # XLBで使用する格子単位系では、 # omega = 1 / (3ν + 0.5) として計算する。 # 緩和係数omegaと動粘性係数viscの関係式から、D2Q9の場合の関係式となる。 omega = 1.0 / (3.0 * visc + 0.5) #ブロック【⑦ XLBを初期化】の解説：LBM計算をXLBで行う計算手法の基本設定 # ============================================================ # ⑦ XLBの初期化 # ============================================================ # XLBに、 # ・D2Q9: velocity_set # ・JAX: compute_backend # ・FP32: precision_policy # を標準設定として登録する。 # 以降のXLBオブジェクトが、この設定を利用できるようになる。 xlb.init( velocity_set=velocity_set, default_backend=compute_backend, default_precision_policy=precision_policy, ) #ブロック【⑧ 計算格子を作成】の解説：設定した解析空間の計算格子を作成する # ============================================================ # ⑧ 計算格子を作成 # ============================================================ # XLBのgrid_factoryを使用して、 # 300 × 120の2次元計算格子を作成する。 # イメージ： # y # ↑ 300 × 120 # ┌─────────────┐ # │ | # │ | # └─────────────┘ → x grid = grid_factory( grid_shape, compute_backend=compute_backend, ) #ブロック【⑨ 境界を取得】の解説：２次元計算領域の４辺の外周格子を取得する # ============================================================ # ⑨ 計算領域の境界を取得 # ============================================================ # 計算領域の外周セルを取得する。 # bounding_box_indices()は、 # 格子の ・left ・right ・top ・bottom # などの境界位置を取得するために使用する。 box = grid.bounding_box_indices() # ------------------------------------------------------------ # 外周端を除いた境界 # ------------------------------------------------------------ # remove_edges=Trueにすると、 # 境界の角などの重複セルを除いた境界を取得する。 # 流入口・流出口の境界条件を設定するときに利用する。 box_no_edge = grid.bounding_box_indices( remove_edges=True ) #ブロック【⑩ 左側 → 流入口】の解説：解析空間の左側を流入部として定義する # ============================================================ # ⑩ 流入口 # ============================================================ # 左側境界を流入口として使用する。 # 流体は → → → → → # のように左から右へ流れる。 inlet = box_no_edge["left"] #ブロック【⑪ 右側 → 流出口】の解説：解析空間の右側を流出部として定義する # ============================================================ # ⑪ 流出口 # ============================================================ # 右側境界を流出口として使用する。 # 流体は右側から計算領域の外へ出ていく。 outlet = box_no_edge["right"] #ブロック【⑫ 上下 → 固体壁】の解説：解析空間の上下において固定壁として処理する # ============================================================ # ⑫ 上下の壁 # ============================================================ # 上側と下側の境界を固体壁として設定する。 # velocity_set.dは空間の次元数。 # D2Q9なのでd=2。 # x方向・y方向の境界情報を組み合わせる。 walls = [ box["bottom"][i] + box["top"][i] for i in range(velocity_set.d) ] # 同じセルが重複している場合に整理する。 walls = np.unique( np.array(walls), axis=-1 ).tolist() #ブロック【⑬ 格子座標を作成】の解説：設定した大きさの解析空間に格子座標を設定する # ============================================================ # ⑬ 格子座標を作る # ============================================================ # x方向の格子番号を作る。 # 0, 1, 2, ..., 299 x = np.arange(Nx) # y方向の格子番号を作る。 # 0, 1, 2, ..., 119 y = np.arange(Ny) # x-yの2次元座標を作る。 # X[i,j]： x方向の位置 Y[i,j]： y方向の位置 # これによって、すべての格子セルについて # 「正方形の内部かどうか」を判定できる。 X, Y = np.meshgrid( x, y, indexing="ij" ) #ブロック【⑭ 正方形の角度を設定】の解説：角度で設定された回転をラジアン単位に変換する # ============================================================ # ⑭ 正方形を回転させる # ============================================================ # 角度を度からラジアンに変換する。 # NumPyのsin/cosはラジアンを使用するため、 # 入力したデグリーの角度を変換する必要がある。 theta = np.deg2rad(square_angle) # ------------------------------------------------------------ # 正方形中心を原点とした座標に変換 # ------------------------------------------------------------ # すべての格子点について、 # 正方形の中心からの相対座標を計算する。 # Xc： x方向の中心からの距離 Yc： y方向の中心からの距離 Xc = X - square_cx Yc = Y - square_cy #ブロック【⑮ 回転座標に変換】の解説：回転した正方形を座標変換して準備する # ============================================================ # ⑮ 回転座標系へ変換 # ============================================================ # 正方形をangleだけ回転した状態を、 # 数学的に判定するために座標変換を行う。 # 座領変換のために、回転行列を利用している。 # Xr = [ cosθ sinθ ][ Xc ] # [ -sinθ cosθ ][ Yc ] # この変換によって、回転した正方形を # 「回転していない正方形」として判定できる。 Xr = ( np.cos(theta) * Xc + np.sin(theta) * Yc ) # y方向についても同じ座標変換を行う。 Yr = ( -np.sin(theta) * Xc + np.cos(theta) * Yc ) #ブロック【⑯ 正方形セルを判定】の解説：格子上での正方形に対応する位置を判定する # ============================================================ # ⑯ 正方形内部のセルを判定 # ============================================================ # 正方形の半分の長さを計算する。 # square_size = 40 なので、 # half = 20 となる。 half = square_size / 2.0 # ------------------------------------------------------------ # 正方形内部の判定 # ------------------------------------------------------------ # 回転後の座標Xr,Yrについて、 # -20 <= Xr <= 20 # -20 <= Yr <= 20 # を満たすセルを正方形内部とする。 # つまり、回転座標系で見ると # 普通の40×40の正方形になっている。 square_mask = ( (np.abs(Xr) <= half) & (np.abs(Yr) <= half) ) #ブロック【⑰ 固体セルのインデックスを作成】の解説：格子のなかで対象物体の部分の指標を設定する # ============================================================ # ⑰ 正方形セルのインデックスを取得 # ============================================================ # square_maskがTrueになったセルの # x,yインデックスを取得する。 # これがXLBに、「ここは固体です」 # と教えるための情報になる。 indices = np.where(square_mask) # XLBが利用できる境界インデックス形式に変換する。 # 正方形のセルについて、 # x方向、y方向の位置をリストとして保存する。 square = [ indices[i].astype(np.int32).tolist() for i in range(velocity_set.d) ] #ブロック【⑱ 流入口 BC】の解説：流入部の一様流の速度境界条件を設定する # ============================================================ # ⑱ 流入口の境界条件 # ============================================================ # 左側から速度を指定して流体を入れる。 # "velocity"： 速度境界条件 # prescribed_value=(u_in, 0.0) は、 # u_x = 0.05 # u_y = 0 を意味する。 bc_inlet = RegularizedBC( "velocity", prescribed_value=(u_in, 0.0), indices=inlet, ) #ブロック【⑲ 上下壁 BC】の解説：２次元解析空間の上下の境界条件を設定する # ============================================================ # ⑲ 上下壁の境界条件 # ============================================================ # 上下の壁を固体壁として扱う。 # HalfwayBounceBackは、壁に向かって進んできた流体粒子を # 壁で反射させる境界条件。 # イメージ # 流体 → # │壁 # │ # ← 反射 │ bc_walls = HalfwayBounceBackBC( indices=walls ) #ブロック【⑳ 正方形 BC】の解説：空間内の物体の対する境界条件を設定する # ============================================================ # ⑳ 正方形物体の境界条件 # ============================================================ # 正方形も固体として扱う。 # 流体が正方形にぶつかった場合、HalfwayBounceBackによって反射する。 # その結果、以下のような流れが発生する。 # → → → → # ■■■ # → → → ■■■ # ■■■ # → → → # LBMで固体壁を表現するための境界条件で、 # 壁の「半分の位置」で粒子を跳ね返す方法です。 bc_square = HalfwayBounceBackBC( indices=square ) #ブロック【㉑ 流出部の境界条件】の解説：流出部において出口境界条件を設定する # ============================================================ # ㉑ 流出口の境界条件 # ============================================================ # 右側から流体が計算領域の外へ出ていく出口境界条件を設定する。 # ExtrapolationOutflowBCは、出口付近の流れを外側へ外挿して # 流体を計算領域から出すために使用する。 bc_outlet = ExtrapolationOutflowBC( indices=outlet ) #ブロック【㉒ 境界条件をまとめる】の解説：解析空間の４種類の境界条件をまとめる # ============================================================ # ㉒ 境界条件をまとめる # ============================================================ # ここまでで作った境界条件を、1つのリストにまとめる。 # これをStepperに渡すことで、XLBが各境界で適切な処理を行う。 boundary_conditions = [ bc_walls, bc_inlet, bc_outlet, bc_square, ] #ブロック【㉓ Navier-Stokes Stepper】の解説：CFDのナビエストークス流れのLBM対応 # ============================================================ # ㉓ Navier-Stokes Stepper # ============================================================ # XLBの時間発展計算を担当するStepperを作成する。 # IncompressibleNavierStokesStepper： # 非圧縮Navier-Stokes流れを # LBMによって時間発展させる。 # collision_type="BGK"： # 衝突モデルとしてBGKを使用する。 stepper = IncompressibleNavierStokesStepper( grid=grid, boundary_conditions=boundary_conditions, collision_type="BGK", ) #ブロック【㉔ LBM分布関数を初期化】の解説：粒子の挙動を定める分布関数の初期化 # ============================================================ # ㉔ 初期分布関数を作成 # ============================================================ # LBMでは、f_i # という「各方向へ移動する粒子の分布関数」を計算する。 # f_0： 現在の分布関数 # f_1： 次の計算で使用する分布関数 # bc_mask： 境界条件の情報 # missing_mask： 不足している隣接セルなどの情報 # を取得する。 f_0, f_1, bc_mask, missing_mask = ( stepper.prepare_fields() ) #ブロック【㉕ 密度・速度計算】の解説：LBMで密度や速度などのマクロな物理量の計算 # ============================================================ # ㉕ 巨視的物理量を計算するオブジェクト # ============================================================ # LBMでは、分布関数 f から、 # 密度 rho 速度 u # などの通常の流体力学的な物理量を計算する。 # その計算を担当するのがMacroscopic。 macro = Macroscopic( compute_backend=ComputeBackend.JAX, precision_policy=precision_policy, velocity_set=xlb.velocity_set.D2Q9( precision_policy=precision_policy, compute_backend=ComputeBackend.JAX, ), ) #ブロック【㉖ 計算条件を表示】の解説：解析の全体設定条件をまとめて出力する # ============================================================ # ㉖ 計算条件を画面に表示 # ============================================================ print() print("======================================") print(" XLB : Flow Around Rotated Square") print("======================================") # 計算格子を表示する。 print(f"Grid : {Nx} x {Ny}") # 今回の計算バックエンドを表示する。 print(f"Backend : JAX-CPU") # 流入速度を表示する。 print(f"Velocity : {u_in}") # Reynolds数を表示する。 print(f"Reynolds : {Re}") # 正方形のサイズを表示する。 print(f"Square size : {square_size} x {square_size}") # 正方形の中心位置を表示する。 print( f"Square center : " f"({square_cx}, {square_cy})" ) # 正方形の回転角度を表示する。 print(f"Square angle : {square_angle} deg") # 動粘性係数を表示する。 print(f"Viscosity : {visc}") # BGK緩和係数を表示する。 print(f"Omega : {omega}") print("======================================") print() #ブロック【㉗ 結果保存関数】の解説：必要量を計算してVTK,PNG形式で保存する # ============================================================ # ㉗ 結果を保存する関数 # ============================================================ # 一定時間ごとに呼び出して、 # その時点の流れ場を保存する関数。 # step： 現在の時間ステップ # f： 現在のLBM分布関数 def post_process(step, f): # -------------------------------------------------------- # 密度と速度を計算 # -------------------------------------------------------- # LBMの分布関数fから、 # rho： 密度 u： 速度 を計算する。 rho, u = macro(f) # -------------------------------------------------------- # 外周1セルを除去 # -------------------------------------------------------- # 結果表示用に外周セルを除く。 # rhoは、[密度] のデータ。 rho = rho[:, 1:-1, 1:-1] # uは2次元速度なので、 # u[0]： x方向速度 u[1]： y方向速度 となる。 u = u[:, 1:-1, 1:-1] # -------------------------------------------------------- # 速度の大きさ # -------------------------------------------------------- # 速度ベクトル u = (u_x, u_y) の大きさ # |u| = sqrt(u_x^2 + u_y^2) # を計算する。 u_magnitude = jnp.sqrt( u[0] ** 2 + u[1] ** 2 ) # ======================================================== # 保存する物理量 # ======================================================== # ParaViewなどで確認できるように、物理量を辞書にまとめる。 fields = { # 密度 "rho": rho[0], # x方向速度 "u_x": u[0], # y方向速度 "u_y": u[1], # 速度の大きさ "u_magnitude": u_magnitude, } # ======================================================== # VTK形式で保存 # ======================================================== # VTKファイルとして保存する # ParaViewを使用すれば、速度分布などをカラー表示できる。 save_fields_vtk( fields, timestep=step, prefix="square_flow" ) # ======================================================== # PNG画像として保存 # ======================================================== # 速度の大きさを画像として保存する。 # 解析結果を簡単に確認するときに便利。 save_image( fields["u_magnitude"], timestep=step, prefix="square_flow" ) # ======================================================== # 計算状況を表示 # ======================================================== # 現在のステップと最大速度を表示する。 print( f"step = {step:6d} : " f"max velocity = " f"{float(jnp.max(u_magnitude)):.6f}" ) #ブロック【㉘ 計算ステップ設定】の解説：全体の解析ステップ数と結果保存間隔 # ============================================================ # ㉘ 時間発展計算 # ============================================================ # 総計算ステップ数。 # LBMでは、この時間ステップを1回ずつ進めながら # 流れを計算する。 num_steps = 5000 # ------------------------------------------------------------ # 結果を保存する間隔 # ------------------------------------------------------------ # 500ステップごとに結果を保存する。 # 5000ステップ計算するので、 # 0 → 500 → 1000 → ... → 5000 の結果が保存される。 post_process_interval = 500 #ブロック【㉙ タイマー開始】の解説：解析時間測定のため開始時間を記録する # ============================================================ # ㉙ 計算時間測定開始 # ============================================================ # 計算開始時刻を記録する。 start_time = time.time() #ブロック【㉚ LBM時間発展】の解説：指定ステップまで解析を進めて途中結果を保存する # ============================================================ # ㉚ LBM時間発展ループ # ============================================================ # 0から4999まで5000回、 # LBMの時間発展を行う。 for step in range(num_steps): # -------------------------------------------------------- # LBMを1ステップ進める # -------------------------------------------------------- # XLBのStepperに # f_0： 現在の分布関数 # f_1： 作業用分布関数 # bc_mask： 境界条件情報 # missing_mask： 不足セル情報 # omega： BGK緩和係数 # step： 現在の時間ステップ # を渡して計算する。 f_0, f_1 = stepper( f_0, f_1, bc_mask, missing_mask, omega, step, ) # -------------------------------------------------------- # 分布関数のバッファを交換 # -------------------------------------------------------- # 今回計算したf_1を、次のステップのf_0として使用する。 # これにより、現在 → 次の時刻 の計算を効率的に進める。 f_0, f_1 = f_1, f_0 # ======================================================== # 結果保存 # ======================================================== # 500ステップごと、または最終ステップで結果を保存する。 if ( step % post_process_interval == 0 or step == num_steps - 1 ): # 現在の分布関数から、密度・速度・速度の大きさを計算して # ファイルに保存する。 post_process( step, f_0 ) #ブロック【㉛ 計算終了】の解説：time関数の現時時点から開始時点を引いて経過時間 # ============================================================ # ㉛ 計算終了 # ============================================================ # 計算終了時刻を取得する。 elapsed = time.time() - start_time #ブロック【㉜ 計算時間表示】の解説：解析の正常完了表示と解析時間の秒単位表示 # ============================================================ # ㉜ 計算時間を表示 # ============================================================ print() print("======================================") print(" Simulation finished") print("======================================") # 計算にかかった秒数。 print(f"Elapsed time : {elapsed:.3f} s") print("======================================") このスクリプト全体の機能と設定の１３項目について、対応するスクリプトと解析条件や解析結果と組み合わせて、下図にまとめます。上記の32ブロックの説明をあわせて、スクリプトの改良に活用してください。 !image.png スクリプト「flow_square_2d-1.py」の実行 まず、以下の【01:JAXcpu導入方法と動作確認】で準備した解析環境を、以下の手順で準備します。XLBの作業用フォルダ「~/XLBjaxcpuWork」において、独自例題フォルダ「FLOW2D」を作っておきます。 qiita.com PowerShell 7.6.6 PS C:\Windows\System32> wsl -l -v NAME STATE VERSION * Ubuntu-24.04 Stopped 2 Ubuntu-24.04-XLBjaxCPU Stopped 2 PS C:\Windows\System32> wsl -d Ubuntu-24.04-XLBjaxCPU dalab@dalab2312nb01:/mnt/c/Windows/System32$ cd dalab@dalab2312nb01:~$ ls XLBjaxCPU XLBjaxcpuWork dalab@dalab2312nb01:~$ source ~/XLBjaxCPU/bin/activate (XLBjaxCPU) dalab@dalab2312nb01:~$ cd ~/XLBjaxcpuWork/ (XLBjaxCPU) dalab@dalab2312nb01:~/XLBjaxcpuWork$ ls LDC2D XLB-git (XLBjaxCPU) dalab@dalab2312nb01:~/XLBjaxcpuWork$ mkdir FLOW2D (XLBjaxCPU) dalab@dalab2312nb01:~/XLBjaxcpuWork$ ls FLOW2D LDC2D XLB-git このなかで、「２次元問題として風洞中にある矩形周りの流れ」を、XLBで解析するために、以上で説明したスクリプト「flow_square_2d-1.py」を、これらの32ブロックを繋げて１つのファイルとして作ります。 (XLBjaxCPU) dalab@dalab2312nb01:~/XLBjaxcpuWork/TEST2D$ cd FLOW2D/ (XLBjaxCPU) dalab@dalab2312nb01:~/XLBjaxcpuWork/FLOW2D$ vi flow_square_2d-1.py (XLBjaxCPU) dalab@dalab2312nb01:~/XLBjaxcpuWork/FLOW2D$ ls flow_square_2d-1.py これを、pythonで実行します。最初にGPUは未設定なので、CPUで計算を進めると表示され、４か所の境界条件の定義が登録されます。 そのあとに、独自例題の解析条件設定が表示されます。そのあとで、解析全体5000ステップで、500ステップ毎に途中結果の保存のログが表示されます。ここでは、500ステップでの解析時間と最大速度も表示されます。 ここで、流入口速度は$u_{\mathrm{in}}=0.05\quad [\mathrm{lu/ts}]$ と考えます。ここで lu は lattice unit（格子数）、ts は time step（時間ステップ）です。つまり、１時間ステップあたり、進む格子数で定義します。この単位で最大速度を見ると、最終5000ステップで0.134なので約2.7倍になっており、正方形の上下の狭い空間を、速く流れていることがわります。 最後に、解析時間が約10秒と表示されます。格子数300x120=36000で、全5000ステップを、Intel Corei9で実行した結果になります。 なお、Windwsのタスクマネージャーで確認すると、CPUで並列処理していることがわかります。 (XLBjaxCPU) dalab@dalab2312nb01:~/XLBjaxcpuWork/FLOW2D$ python flow_square_2d-1.py No GPU support is available; CPU fallback will be used. registered bc RegularizedBC_8332170304619 with id 1 registered bc HalfwayBounceBackBC_8332175202663 with id 2 WARNING! Assuming no-slip condition for BC type = HalfwayBounceBackBC! registered bc HalfwayBounceBackBC_8332163779297 with id 3 WARNING! Assuming no-slip condition for BC type = HalfwayBounceBackBC! registered bc ExtrapolationOutflowBC_8332163779354 with id 4 Warp 1.17.0 initialized: CUDA Toolkit 12.9, Driver 13.2 Devices: "cpu" : "x86_64" "cuda:0" : "NVIDIA GeForce RTX 3070 Ti Laptop GPU" (8 GiB, sm_86, mempool enabled) Kernel cache: /home/dalab/.cache/warp/1.17.0 ====================================== XLB : Flow Around Rotated Square ====================================== Grid : 300 x 120 Backend : JAX-CPU Velocity : 0.05 Reynolds : 100.0 Square size : 40 x 40 Square center : (75, 60) Square angle : 30.0 deg Viscosity : 0.02 Omega : 1.7857142857142856 ====================================== Saved ./square_flow_0000000.vtk in 0.020215 seconds. step = 0 : max velocity = 0.000000 Saved ./square_flow_0000500.vtk in 0.002383 seconds. step = 500 : max velocity = 0.156543 Saved ./square_flow_0001000.vtk in 0.002592 seconds. step = 1000 : max velocity = 0.129826 Saved ./square_flow_0001500.vtk in 0.002317 seconds. step = 1500 : max velocity = 0.130439 Saved ./square_flow_0002000.vtk in 0.002421 seconds. step = 2000 : max velocity = 0.133740 Saved ./square_flow_0002500.vtk in 0.002411 seconds. step = 2500 : max velocity = 0.135603 Saved ./square_flow_0003000.vtk in 0.002438 seconds. step = 3000 : max velocity = 0.132879 Saved ./square_flow_0003500.vtk in 0.002409 seconds. step = 3500 : max velocity = 0.132424 Saved ./square_flow_0004000.vtk in 0.007832 seconds. step = 4000 : max velocity = 0.138548 Saved ./square_flow_0004500.vtk in 0.007650 seconds. step = 4500 : max velocity = 0.137304 Saved ./square_flow_0004999.vtk in 0.002434 seconds. step = 4999 : max velocity = 0.133921 ====================================== Simulation finished ====================================== Elapsed time : 10.336 s ====================================== 解析結果として、500ステップ毎の「可視化用vtkファイル」と「速度分布pngファイル」が、保存されています。 (XLBjaxCPU) dalab@dalab2312nb01:~/XLBjaxcpuWork/FLOW2D$ ls flow_square_2d-1.py square_flow_0002000.vtk square_flow_0004999.vtk square_flow_3000.png square_flow_0000.png square_flow_0002500.vtk square_flow_0500.png square_flow_3500.png square_flow_0000000.vtk square_flow_0003000.vtk square_flow_1000.png square_flow_4000.png square_flow_0000500.vtk square_flow_0003500.vtk square_flow_1500.png square_flow_4500.png square_flow_0001000.vtk square_flow_0004000.vtk square_flow_2000.png square_flow_4999.png square_flow_0001500.vtk square_flow_0004500.vtk square_flow_2500.png 速度分布のPNG画像ファイルを、ステップ0500-4999の10枚分を確認します。前半で流れが成長し、後半でカルマン渦が確認できます。これは、正方形が30度回転しており、上下の流れが非対称になっており、カルマン渦が誘導されたと思います。 !image.png 次に、可視化用のVTKデータファイル「square_flow_...vtk」を、ParaViewを用いて確認します。 !image.png この可視化用データでは、密度rhoを表示することもできます。 !image.png さらに、MP4形式の動画ファイルを作成することもできます。ParaViewの手順（save animation）を確認して試してください。 以上で説明したスクリプト「flow_square_2d-1.py」において、変更可能なパラメータを、下図にまとめます。まず最初は、スクリプトの設定通りで解析を実行して、動作確認をします。そのあとは、この図に示すいろいろなパラメータを変更して、比較検討してみてください。 !image.png 例えば、正方形の回転角度０度とすると、傾かない状態の速度分布は、以下のようになります。流れは成長しますが、上下対称なので、カルマン渦は発生していません。 !image.png コンパクト版スクリプト「flow_square_2d-0.py」 スクリプト全体の把握し、機能拡張などを目的として、コメントなどを削除して、ログの画面表示も圧縮して、コンパクトにしたスクリプト「flow_square_2d-0.py」を以下に示します。こちらを用いて一括コピーしても、同じ解析結果が得られます。なお、対応するブロック①～㉜について、該当する行末に記してあります。 import os os.environ["JAX_PLATFORMS"] = "cpu" import time import numpy as np import jax.numpy as jnp import xlb 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, RegularizedBC, ExtrapolationOutflowBC from xlb.operator.macroscopic import Macroscopic from xlb.utils import save_fields_vtk, save_image #① Nx, Ny = 300, 120 u_in, Re = 0.05, 100.0 #② square_size, square_cx, square_cy, square_angle = 40, Nx // 4, Ny // 2, 30.0 #③ compute_backend = ComputeBackend.JAX precision_policy = PrecisionPolicy.FP32FP32 #④ velocity_set = xlb.velocity_set.D2Q9(precision_policy=precision_policy, compute_backend=compute_backend) #⑤ visc = u_in * square_size / Re omega = 1.0 / (3.0 * visc + 0.5) #⑥ xlb.init(velocity_set=velocity_set, default_backend=compute_backend, default_precision_policy=precision_policy) #⑦ grid = grid_factory((Nx, Ny), compute_backend=compute_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"][i] + box["top"][i] for i in range(velocity_set.d)] walls = np.unique(np.array(walls), axis=-1).tolist() #⑫ X, Y = np.meshgrid(np.arange(Nx), np.arange(Ny), indexing="ij") #⑬ theta = np.deg2rad(square_angle) Xc, Yc = X - square_cx, Y - square_cy Xr = np.cos(theta) * Xc + np.sin(theta) * Yc Yr = -np.sin(theta) * Xc + np.cos(theta) * Yc #⑭ half = square_size / 2.0 #⑮ square_mask = (np.abs(Xr) <= half) & (np.abs(Yr) <= half) #⑯ indices = np.where(square_mask) square = [indices[i].astype(np.int32).tolist() for i in range(velocity_set.d)] #⑰ bc_inlet = RegularizedBC("velocity", prescribed_value=(u_in, 0.0), indices=inlet) #⑱ bc_walls = HalfwayBounceBackBC(indices=walls) #⑲ bc_square = HalfwayBounceBackBC(indices=square) #⑳ bc_outlet = ExtrapolationOutflowBC(indices=outlet) #㉑ boundary_conditions = [bc_walls, bc_inlet, bc_outlet, bc_square] #㉒ stepper = IncompressibleNavierStokesStepper(grid=grid, boundary_conditions=boundary_conditions, collision_type="BGK") #㉓ f_0, f_1, bc_mask, missing_mask = stepper.prepare_fields() #㉔ macro = Macroscopic(compute_backend=ComputeBackend.JAX, precision_policy=precision_policy, velocity_set=xlb.velocity_set.D2Q9(precision_policy=precision_policy, compute_backend=ComputeBackend.JAX)) #㉕ print(f"Grid={Nx}x{Ny}, Backend=JAX-CPU, U={u_in}, Re={Re}, Square={square_size}x{square_size}, Center=({square_cx},{square_cy}), Angle={square_angle}°, ν={visc}, ω={omega}") #㉖ def post_process(step, f): #㉗ rho, u = macro(f) rho, u = rho[:, 1:-1, 1:-1], u[:, 1:-1, 1:-1] u_magnitude = jnp.sqrt(u[0] ** 2 + u[1] ** 2) fields = {"rho": rho[0], "u_x": u[0], "u_y": u[1], "u_magnitude": u_magnitude} save_fields_vtk(fields, timestep=step, prefix="square_flow") save_image(fields["u_magnitude"], timestep=step, prefix="square_flow") print(f"step={step:6d}, max velocity={float(jnp.max(u_magnitude)):.6f}") num_steps, post_process_interval = 5000, 500 #㉘ start_time = time.time() #㉙ 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 % post_process_interval == 0 or step == num_steps - 1: post_process(step, f_0) elapsed = time.time() - start_time #㉛ print(f"Simulation finished: {elapsed:.3f} s") #㉜ LBMの考え方を改めて確認する 以上で、「２次元問題として風洞中にある矩形周りの流れ」 を、XLBのスクリプトを用いて、LBM解析で求めることが出来ました。ここで、LBMを初めて学ぶ場合、まず 「流体そのものを直接動かすのではなく、各格子点にある“粒子の分布”を動かして、結果として流体の速度や密度を得る」 として考えると理解しやすいです。 今回のXLBの例題を理解するために、D2Q9 → 分布関数 → Collision → Streaming → 速度・密度の順に、できるだけ具体的に説明します。 1. LBMの基本的な考え方 通常のCFDでは、Navier-Stokes方程式を直接解いて、$u_x,\quad u_y,\quad p$ などを求めます。 一方、LBMでは少し違います。以下の手順、つまり、「分布関数を計算することで、流体の速度・密度を間接的に求める」 という方法です。 LBM(XLB) → 分布関数 f_i を計算 → Collision（衝突） → Streaming（移動） → 分布関数 f_i の更新 → 密度 ρ、速度 u を計算 2. 「分布関数」とは何か？ ここがLBMで最も重要な部分です。ある格子点を考えると、「D2Q9」 では、1つの格子点に9方向の分布関数があります。 !image.png 一般に、$f_0,f_1,f_2,\cdots,f_8$ と表します。例えば、以下のように、それぞれの方向に対応しています。 f0 → その場にとどまる成分 f1 → 右 f2 → 上 f3 → 左 f4 → 下 f5 → 右上 f6 → 左上 f7 → 左下 f8 → 右下 3. 分布関数は「粒子の量」と考える 例えば、ある格子点で、以下の状態だったとします。 右方向：f1 = 0.20・左方向：f3 = 0.10・上方向：f2 = 0.12・下方向：f4 = 0.12 右向きの分布が左向きより大きいので、「右 →→」という流れが生じます。つまり、各方向にどれだけ“流体の分布”が存在しているか を表しているのが f_i です。ここで重要なのは、f_i は普通の意味での「実際の粒子1個1個」ではないことです。LBMでは、流体の状態を9方向の分布として表現していると考える方が適切です。 4. 分布関数から密度を求める 9個の分布関数を全部足すと、密度になります。$\rho=\sum_i f_i$ つまり、「f0 + f1 +...+ f8：全部足すと、密度ρ」になります。 5. 分布関数から速度を求める 速度は、各方向の分布に「方向」を掛けて足し合わせます。$\rho\mathbf{u}=\sum_i f_i\mathbf{c}_i$ したがって、$\mathbf{u}=\frac{1}{\rho}\sum_i f_i\mathbf{c}_i$ となります。ここで、$\mathbf{c}_i$ は各方向の速度です。例えば、「右方向：+x 左方向：-x 上方向：+y 下方向：-y」です。 6. では「Collision」とは何か？ ここからLBMの核心です。英語のCollisionは、日本語では「衝突」です。ただし、流体解析のLBMでは、 分布関数を局所的に平衡状態へ近づける計算 と考える方が分かりやすいです。 7. Collisionのイメージ 例えば、ある格子点で、方向別の値が、以下のような状態だったとします。流体としては少し「乱れた」分布になっています。 現在の分布 右 ████████ 左 ██ 上 ████ 下 ███ 斜め ███... そこで、「現在の分布 → Collision → 平衡状態に近づける」という操作をします。 8. BGKモデル 今回のXLBのスクリプトでは、「collision_type="BGK"」 を使用しています。BGKでは、非常に簡単に書くと、$f_i^{*}=f_i+\omega(f_i^{eq}-f_i)$ です。 ここで、$f_i$：現在の分布関数、$f_i^{eq}$：平衡分布関数、$f_i^*$：Collision後の分布関数、$\omega$：緩和係数、です。 9. この式を日本語にすると 式 $f_i^{*}=f_i+\omega(f_i^{eq}-f_i)$ を見ると、$f_i^{eq}-f_i$ は、「現在の分布と平衡状態との差」です。したがって、以下のような流れになります。 現在の分布 → 平衡分布との差 → 現在の分布を修正 → Collision後の分布 [差を計算] [ωを掛ける] 10.緩和係数 ωは何を意味する？ 今回のスクリプトでは、$\omega=\frac{1}{3\nu+0.5}$ としていました。例えば、Re=100 の場合、ν=0.02 なので、$\omega=\frac{1}{3(0.02)+0.5}\approx1.4286$ となります。つまり、「Re → 粘性 ν → 緩和係数 ω → Collisionの強さ」という関係になっています。 11. Collisionだけでは流体は流れない ここが初心者には特に重要です。Collisionをすると、「格子点A ● → 分布が変化します。」しかし、それだけでは隣の格子点には移動していません。そこで次に、Streaming を行います。 12. Streamingとは？ Streamingは、分布関数を、その分布が向いている隣の格子点へ移動させる 処理です。例えば、以下のように、左図の３つの格子点があるとします。A点に存在する右向き分布 $f_1$ は、Streamingによって、右図のように移動します。 A B C [Streaming] A f1 B ● → ● → ● ● → ● つまり、$f_i(\mathbf{x}+\mathbf{c}_i,t+1)=f_i^*(\mathbf{x},t)$ です。 13. これが「流れる」の正体 例えば右向きの分布だけを見ると、それぞれの格子点にある右向き分布が、以下のように隣へ移動します。これを大量の格子点で繰り返します。すると、クロな流れが現れます。 時刻 t 時刻 t+1 ● →1 ● →2 ● →3 ● ● →0 ● →1 ● →2 ● 14. CollisionとStreamingの違い ここを明確に区別するとLBMがかなり分かりやすくなります。 | 処理 | 何をする？ | イメージ | |---|---|---| | Collision | 分布を平衡状態に近づける | その場で分布を混ぜる | | Streaming | 分布を隣の格子へ移動する | 分布を運ぶ | 15. LBMの1ステップ LBMの基本的な1ステップを非常に単純化すると、以下のような流れになります。 !image.png 16. 今回の正方形まわり流れで考える 今回のモデルは、２次元の長方形の流体空間の内部に、正方形の物体があります。左から流体が入ってきて、LBMでは、各格子点で以下の処理を繰り返します。 分布関数 f_i → Collision → Streaming → 境界条件 → 次の時刻の f_i 17. 正方形にぶつかるとどうなる？ Streamingをすると、分布が正方形の内部へ入ろうとする場合があります。しかし、正方形は固体なので、そのまま通過させません。そこで、Bounce Back などの境界条件を使います。 イメージとしては、流体の格子からの分布を反射させます。今回のスクリプトでは、「HalfwayBounceBackBC」を使っています。 18. そして後流ができる 正方形の後ろでは、分布関数の移動とCollisionを何度も繰り返すことで、複雑な流れが形成されます。これをマクロに見ると、後流・再循環・渦 として現れます。 19. 最後に「速度」を取り出す LBMが直接計算している中心的な量は分布関数 $f_i$ です。そこから、$\rho=\sum_i f_i$と、$\mathbf{u}=\frac{1}{\rho}\sum_i f_i\mathbf{c}_i$ を計算します。今回のXLBスクリプトの、「rho, u = macro(f)」は、まさにこの処理です。 !image.png 20. XLBのスクリプトとの対応 今回のスクリプトをLBMの考え方に対応させると、次のような、説明図になります。 !image.png 21. f_0 と f_1 は何？ 今回のプログラムには、「f_0, f_1 = stepper.prepare_fields()」がありました。ここは少し注意が必要です。f_0とf_1はD2Q9の9方向そのものではありません。 D2Q9の方向は、$f_0,\ldots,f_8$ ですが、XLBのスクリプトに出てくる「f_0 f_1」は、時間発展の前後の分布関数を格納するための2つの配列（バッファ） として使われています。 そして、以下の処理によって、計算した結果を次のステップへ渡しています。ここはLBMを初めて学ぶと非常に混乱しやすいところです。 f_0, f_1 = stepper( f_0, f_1, bc_mask, missing_mask, omega, step ) f_0, f_1 = f_1, f_0 22. 一番簡単な覚え方 LBMはまず、次の3語を覚えるとよいです。これらを何千回も繰り返すと、最終的に 「流体の流れ」 が現れます。 【分布関数：「各方向にどれだけ流体が分布しているか」】 ↓【Collision：「その場所で分布を平衡状態に近づける」】 ↓【Streaming：「その分布を方向に従って隣の格子へ運ぶ」】 以上で説明してきた、以下の流れで理解すると、XLBのプログラムが何をしているのかがかなり見通しやすくなります。 分布関数 → Collision → Streaming → 境界条件 → 次の時刻 → .... → 正方形まわりの流れ → ux, uy, |u| → 解析結果 これを説明図としてまとめると、以下のようになります。 !image.png