第 26 回数値流体力学シンポジウム
講演番号 B11-4

格子ボルツマン法を用いた高効率な乱流の数値計算に関する研究
A study on highly efficient numerical simulation for turbulent flow using lattice Boltzmann model○ 安田 孝宏, 滋賀県立大, 滋賀県彦根市八坂町 2500, E-mail:tkyasuda@mech.usp.ac.jp
橋本 知久, 近畿大, 東大阪市小若江 3-4-1, E-mail:hasimoto@mech.kindai.ac.jp
丹野 格, 筑波技術大, つくば市天久保 4-3-15, E-mail: itaru_t@a.tsukuba-tech.ac.jp
田中 嘉宏, 東洋ゴム工業㈱, 兵庫県伊丹市藤ノ木 2-2-13, E-mail: y-tanaka@toyo-rubber.co.jp
南川 久人, 滋賀県立大, 滋賀県彦根市八坂町 2500, E-mail:minagawa@mech.usp.ac.jp
森西 晃嗣, 京工繊大, 京都市左京区松ヶ崎御所海道町, E-mail:morinisi@kit.ac.jp
里深 信行, 京工繊大/滋賀県立大名誉教授, E-mail: satofuka.n@office.usp.ac.jp
TakahiroYasuda, The University of Shiga Prefecture, 2500 Hassaka-cho, Hikone-shi, Shiga 522-8533
Tomohisa Hashimoto, Kinki University, 3-4-1 Kowakae, Higashi-Osaka, Osaka, 577-8502
Itaru Tanno, Tsukuba University of Technology, 4-3-15 Amakubo, Tsukuba, Ibaraki, 305-8520
YoshihiroTanaka, TOYO TIRE&RUBBER CO.,LTD., 2-2-13 Fujinoki, Itami, Hyogo, 664-0847
Hisato Minagawa, The University of Shiga Prefecture, 2500 Hassaka-cho, Hikone-shi, Shiga 522-8533
Koji Morinishi, Kyoto Institute of Technology, Matsugasaki, Sakyo-ku, Kyoto, 606-8585
Nobuyuki Satofuka, Nobuyuki Satofuka, Professor emeritus, Kyoto Institute of Technology/The
University of Shiga Prefecture
Quasi-equilibrium Lattice Boltzmann model (QELBM) has been developed in order to get stable solution at highReynolds number flow simulation in recent years. Since QELBM uses only two relaxation times to stable thecalculation, the algorithm is simpler than Multi-Relaxation Time model (MRT), which is the model often used to stable
the solution, and 13 velocity model for 3-dimensional calculation, which cannot introduce in Entropic Lattice
Boltzmann model, can be used. However, the Application of QELBM for turbulence flow has not been done yet, thus
QELBM can capture the turbulence structure by using coarse mesh. In this study, we applied QELBM to isotropic
homogeneous turbulence at Reynolds number Re = 62831 in various grid resolutions, and investigated the stability,accuracy and computational effort by comparing with LBGK model and ELBM.
１．緒論
格子ボルツマン法 (LBM)は流体を仮想粒子の集合体とみなし，速度分布関数の並進と衝突を解くことで流れ場を計算する手法であり，アルゴリズムが簡単，ポアソン方程式が不要，並列化に向いているという特長から，計算負荷の大きいポアソン方程式を解く従来の差分法等の手法と比べて，非定常・非圧縮流れの高効率な計算に有利である．その一方で，格子ボルツマン法の衝突演算子として広く用いられている単一の緩和時間を用いた LBGK モデルでは高レイノルズ数領域の流れ場で十分な格子解像度がない場合に計算が不安定になることが知られており，工学的にも重要な高レイノルズ数流れへの適用の妨げとなっている．この欠点を克服するために，多緩和時間モデル(MRT) (1)，エントロピック格子ボルツマンモデル(ELBM) (2)，準平衡モデル(QELBM)(3)(4)等が開発されている．特に，QELBM は 2 つの緩和時間を用いて計算を安定化させる手法で， MRT よりもアルゴリズムがシンプルである．また，QELBM は ELBM のように緩和時間を決定するために非線形方程式を計算する必要が無く，さらに，単一の緩和時間を用いる LBGK モデルや ELBM では使用不可能な速度の数が少な
い 3 次元 13 速度（D3Q13）モデルを用いることが可能であるため(4)，高レイノルズ数流れの高効率な数値計算手法として有望である．しかし，QELBM を格子解像度が十分でない乱流場に適用した場合に，格子解像度以下の乱流がどの程度正確に計算可能かどうかについての知見は調査されていない．そこで，本研究では
LBGK モデル，ELBM および QELBM で 2 次元一様等方性乱流の計算を行い，計算の安定性や渦度場および計算速度について比較した．

２．数値計算手法
２－１ 準平衡格子ボルツマンモデル
本研究で用いる QELBM(3)は，広く用いられている衝突項に単一の緩和時間を用いる LBGK モデルとは異なり，式(1)のようにボルツマン方程式の衝突項を 2 段階の緩和メカニズムを用いて表す．

 t f  v  f  1
 f  fC   1  fC  f M fs(1)ここで， f は速度分布関数であり，v はマクロの速度を表している．右辺第 1 項は緩和時間f を用いた制約平衡分布関数 fC への緩和過程であり，右辺第 2 項は緩和時間s を用いた fC から平衡分布関数 fM への緩和過程である．ここで，平衡分布関数は質量と運動量の保存を満たす条件下で H 関数を最小にするように決定され，制約平衡分布関数は質量と運動量に加えて, 応力テンソルの対角成分の保存を満たす条件下でH関数を最小にするように決定される。ここで，式(1)は式(2)で定義される準平衡分布関数 fQE を用ると，式(3)のように LBGK モデルと同様な形式で表される。

f QE f f f M  1   f Cs s (2)1 f  f QE f(3) t f  v  f  

ここで，緩和時間f およびs と動粘性係数 および体積粘性係数
は音速 cs を介してそれぞれ式(4)のような関係があるため，緩和時間f とs の比は式(5)のように粘性比に等しくなる。

   f c s2 ,    s c s21(4)Copyright © 2012 by JSFMfs



第 26 回数値流体力学シンポジウム
講演番号 B11-4
で格子点数 32×32，64×64，128×128 における LBGK モデル，ELBM および QELBM の渦度分布を比較する．図には参照解として，1024×1024 の高解像度格子を用いた LBGK モデルの結果も併せて示す．
表 1 で示したように，32×32 の低い格子解像度では，ELBM でのみ解が得られるが，図 1(a)の渦度分布から分かるように，格子解像度が小さいため，渦構造がほとんど捉えられていない．格子解像度が増加すると図 1(b)，図 1(c)のように徐々に細かい構造が捉えられ，参照解に近づくのが確認できる．しかし，せん断の強い部分の渦構造についてはこの格子解像度ではまだ捉えられていない．定量的にどの程度の渦構造が捉えられているかはエネルギースペクトル等を用いた調査が必要である．

(5)

ここで，は H 関数が単調に減少する条件を満たすために，よりも大きく設定しなければならないが，その条件を満たす範囲で安定な解が得られるように調節できる．粘性比が1 の場合はQELBM は LBGK モデルと同じモデルとなる。
２－２ 計算条件
本計算では，2 次元一様等方性乱流を取り扱う．速度モデルには D2Q9 モデルを用い，境界条件には周期境界条件を与えた．初~k , k  が以下の条件を満たすように期条件は波数空間の渦度 1 2与えた．
1~k1 , k2 222(6)E k   1 k 2 , k 2  k1  k22k k ３－３ 計算速度
計算速度を調査するために,LBG モデル，ELBM および QELBMの各モデルの CPU 時間を表 2 で比較する．ここで，計算は無次元の物理時間 t*=tU/L が t*= 10 に達するまで行った．計算のプラッ
トフォームはクロック周波数 2.8GHz の CPU を持つ PC であり，並列計算は行っていない．
表 2 から分かるように，同じ格子点数の場合，ELBM や QELBMの計算時間は LBGK モデルの約 2.5 倍，約 1.7 倍となる．これはELBMでは緩和時間を各格子点で調整するのに非線形方程式を解くためであり，QLBM では，計算に 2 段階の緩和過程を用いるためである．しかし，QELBM や ELBM では LBGK モデルで解が得られない 64×64 の格子解像度で計算を行うことが可能であり，その場合には表 2 に示すように LBGK モデルの 16％程度の時間で計算可能である．ただし，図 1 で示したように，格子解像度が低くなると細かい乱流構造を捕獲できないのでその兼ね合いで解法を選ぶ必要がある．
また，緒論で述べたようにQELBMはD3Q13モデルを適用できるため，さらなる計算時間の短縮が期待できる．

2

ここで，E(k)は波数 k におけるエネルギースペクトル，k1 およびk2 は x1 方向および x2 方向の波数をそれぞれ表している．なお，初期の E(k)は以下のように与えた．
2 2 E k   k exp   k 3
 3 (7)格子は格子間隔 1 の等間隔直交格子を用い，格子点数は格子解像度の影響をみるために 32×32，64×64，128×128，1024×1024の 4 種類の格子で計算を実行した．また，流れのマッハ数 Ma は
Ma = U/Cs = 0.087 とし，代表速度 U および計算領域長さに基づくレイノルズ数 Re は Re =62831 である．ここで，Cs は音速で Cs =1/(30.5)である．
３．計算結果および考察
３－１ 計算の安定性計算の安定性を調べるために，表１に LBGK モデル，ELBM お
よび QELBM の安定な解が得られる最大格子点数を示す．表 1 に示すように，ELBM, QELBM の計算の安定性が LBGK よりも向上していることが分かる．特に，ELBM は無条件安定な解法である
ため，本計算で用いた最も少ない格子点数で Re=1.0×105 の計算が行える．ただし，ELBM は各格子点毎に緩和時間を調整し，局所的に数値粘性を調整するため，本来の高レイノルズ数流れが計算できているかどうかについては検証が必要となる．QELBM はELBM に比べて安定性の向上率が少ないが，これは QELBM が無条件安定な解法では無いことに加えて，QELBM が体積粘性を調整し，音波等をダンピングすることにより安定性を向上させる手法であることが考えられる，すなわち，キャビティ流れのような特異点があり，そこからの強い音波が生じやすい問題では安定性の向上率が高いが，本研究のような周期境界を有する流れ場ではそのような圧力振動が少なく，安定性の向上が少ないと考えられる．

４．結論
本研究では LBGK モデル，ELBM および QELBM を 2 次元一様等方性乱流の計算に用い，計算の安定性や精度および計算速度について検証した．その結果，ELBM や QLBM で計算の安定性が向上することが分かった．また，低い格子解像度での計算結果は大まかな渦構造は捉えられるが，強いせん断層の部分では細かい格子の場合と差がみられることが分かった．計算時間は LBGK
モデルで解が得られない 64×64 の格子解像度で計算を行えばQELBM が最も少ないが，必要な精度との兼ね合いで解法を選ぶ必要があることが分かった．今後は，MRT も含めた各手法の比較やエネルギースペクトル等をで，低い格子解像度で得られた結果がどの程度渦構造を捉えられているか定量評価を行う予定である．

Table 1 Minimum number of grids required for obtaining stablesolution
ReLBGKELBM

QELBM

4

1.0×1064×64

32×32

32×32

4

5.0×10128×128

32×3264×64

5128×128

32×32

128×128

1.0×10LBGK(1024×1024)
ELBM(a) 32×32
Fig. 1 Vorticity distribution obtained by LBGK, ELBM and QELBM invarious grid resolution (Continued).

３－２ 渦度分布
計算結果の格子解像度による影響について調べるために，図 12Copyright © 2012 by JSFM

第 26 回数値流体力学シンポジウム
講演番号 B11-4
参考文献(1) D'Humieres,D., Bouzidi, M., and Lallemand, P., “Thirteen-velocity
three-dimensional lattice Boltzmann model”, Physical Review E,Vol. 63, (2001), pp. 066702-1-066702-7.
(2) Karlin , I.V., Gorban, A., Succi, S., Boffi,V.,, “Maximum entropy
principle for lattice kinetic equations”, Physical Review Letters, Vol.
81, No. 1 (1998), pp. 6-9.
(3) Asinari, P, Karlin, I. V., “Quasiequilibrium Lattice Boltzmann
Models with Tunable Bulk Viscosity for Enhancing Stability”,
Physical Review E, Vol. 81, (2010), pp. 016702-1-016702-15.
(4) Karlin, I.and Asinari, P, “Factorization Symmetry in Lattice
Boltzmann Method”, Physica A: Statistical Mechanics and its
Applications, Vol. 389, No. 8 (2010), pp. 1530-1548.

LBGK(1024×1024)

ELBM

QELBM
(b) 64×64

LBGK (1024×1024)

LBGK

ELBM

QELBM
(c) 128×128
Fig. 1 Vorticity distribution obtained by LBGK, ELBM and QELBM invarious grid resolutions.
Table 2 Comparison of CPU time
MethodCPU TimeRatio

ELBM 32×32

4.00.03

ELBM 64×64

29.70.24

QELBM 64×64

19.80.16

LBGK 128×128

121.01.0

ELBM 128×128

302.3

2.50

QELBM 128×128

202.4

1.673Copyright © 2012 by JSFM