Autodesk

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