第 25 回数値流体力学シンポジウム
D12-2

準平衡格子ボルツマンモデルを用いた高レイノルズ数流れの数値計算の高効率化
Study on Highly Efficient Numerical Simulation of Flow in High Reynolds Number Region using
Quasie-Equilibrium Lattice Boltzmann Model○ 安田孝宏, 滋賀県立大, 滋賀県彦根市八坂町 2500, E-mail:tkyasuda@mech.usp.ac.jp
南川久人, 滋賀県立大, 滋賀県彦根市八坂町 2500, E-mail:minagawa@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
森西 晃嗣, 京工繊大, 京都市左京区松ヶ崎御所海道町, 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
Hisato Minagawa, 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
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 3-dimensonal flow has not been done yet,thus the detail of the advantage of QELBM for other methods is not known. In this study, we applied QELBM to flow
in a duct with square cross section at Reynolds number Re = 100, and investigated the stability, accuracy andcomputational effort by comparing with theoretical solution and 13 velocity MRT. As a result, it was found thatQELBM can efficiently calculate the flow without losing the computational stability and accuracy.
度を損なうこと無く,高効率な計算が可能であることを確認した(5)．
本研究では，QELBM が入口・出口境界を伴う場合にも有効であるかどうかを調査するため， 3 次元正方ダクト流れの計算を行い，計算結果の安定性・妥当性，計算速度について検証した．

１．緒論
格子ボルツマン法 (LBM)は流体を仮想粒子の集合体とみなし，速度分布関数の並進と衝突を解くことで流れ場を計算する手法であり，アルゴリズムが簡単，ポアソン方程式が不要，並列化に向いているという特長から，計算負荷の大きいポアソン方程式を解く従来の差分法等の手法と比べて，非定常・非圧縮流れの高効率な計算に有利である．その一方で，格子ボルツマン法の衝突演算子として広く用いられている単一の緩和時間を用いた LBGK モ
デルでは高レイノルズ数領域の流れ場で十分な格子解像度がない場合に計算が不安定になることが知られており，工学的にも重要な高レイノルズ数流れへの適用の妨げとなっている．この欠点を克服するために，多緩和時間モデル(MRT) (1)，エントロピック格子ボルツマンモデル(ELBM) (2)，準平衡モデル(QELBM)(3)(4)と呼ばれるモデルが開発されている．特に，QELBM は 2 つの緩和時間を用いて計算を安定化させる手法で， MRT よりもアルゴリズムがシンプルである．また，QELBM は ELBM のように緩和時間を決定するために非線形方程式を計算する必要が無く，さらに，
ELBM では用いることが出来ない速度の数が最小の 3 次元 13 速度（D3Q13）モデルを用いることが可能であるため，高レイノルズ数流れの高効率な数値計算手法として有望である．しかし，現在までに QELBM を用いて 3 次元計算を行った例は無く，計算の安定性や精度および計算速度がどの程度向上するのかについては不明である．そこで，最近，著者らは D3Q13 準平衡モデルを用いて立方キャビティ流れの数値解析を行い，計算の安定性，計算の精度，計算速度について MRT や LBGK モデルと比較し，QELBMは MRT や LBGK モデルと比べて，計算の安定性や計算結果の精２．数値計算手法
２－１ 準平衡格子ボルツマンモデル
以下に QELBM の計算手法の概略を示す．詳細については文献(4)を参照されたい．
QELBM では式(1)に示すように衝突項を 2 段階の緩和メカニズ
ムを用いて表す．右辺の第 1 項は緩和時間1 での準平衡分布関数f *への緩和であり，第 2 項は緩和時間2 での準平衡状態から平衡分布関数 f M への緩和である．この 2 つの緩和時間はそれぞれ，非平衡の圧力テンソルの対角成分と非対角成分を調節する作用があり，LBGK モデルを D3Q13 モデルに適用する際の弊害となる低マッハ数極限でのエントロピー不一致の問題を回避できる．
11(1) t f  v  f    f  f *    f *  f M 1
2
ここで， f *は，圧力テンソルの対角成分，流体の密度およびマクロな速度の方向成分 uを用いて以下のように表わされる．
1
f *0 ,0 ,0     1   xx   yy   zz 2
1*f  , ,0     xx u x 1  u y    yy u y 1  u x 8
  zz 1  u x  u y 1

Copyright © 2011 by JSFM

第 25 回数値流体力学シンポジウム
D12-2
1f * ,0 ,     xx u x 1  u z    zz u z 1  u x 8  yy 1  u x  u z (2)３．計算結果３－１ 計算の安定性計算の安定性を調べるために，表１に U = 0.05 の場合の安定な解が得られる最大 Re 数を示す．表には比較のために，同じ格子点数の 13 速度の MRT の結果も併せて示す．表 1 に示すように，計算の安定性は QELBM の方が優れている．この結果は著者らが最近行った 3 次元キャビティ流れの結果(5)と異なる傾向であり，この原因として，出口境界の影響が考えられる．しかし，QELBMや MRT では自由度を持つパラメータが含まれていることから，厳密な比較を行うためには，それらの影響を考慮する必要がある．



1f *0 , ,     yy u y 1  u z    zz u z 1  u y 8
  xx 1  u y  u z 

尚，D3Q13 の場合，とには 1 あるいは-1 が用いられる．また，平衡分布関数 f M は式(2)のに以下の Mを代入することで求められる．
1M(3)  u23
本研究では，2 つの緩和時間1 および2 に以下の値を用いる．
 1   , 2  2Table 1 Maximum Reynolds number required obtaining stable solution(4)

ここで，緩和時間 と動粘性係数の間には以下の関係がある．
(5)2式(1)を離散化して式を整理すると，新しい時間ステップでの分布関数は以下の式で与えられる．

g  x  vt ,t  t   1   g  x ,t M*

300

4 H 2  dp  3  dx (7) coshiz H   cos iy H i3  1  1   coshi 2 i 1 2

ここで， u，p およびはそれぞれ x 方向の速度成分，圧力および粘性係数を示し，右辺の無限級数は反復回数を 200 として数値解により近似した．
図 3 に示すように，QELBM で得られた速度分布は下流側に移動するにつれて若干の増加が見られる．この傾向は特に U = 0.1の場合で顕著である．これは，ダクト壁面による摩擦の影響で圧力損失が生じ，密度が減少したため，質量保存の法則に従って速度が増加したものと考えられる．故に，図 4 に示すように縦軸を質量流量で整理すると，一様流入速度に関わらず流れ方向のどの位置においても値がほぼ理論値と一致する．また，圧力勾配を理論値と比較すると，U = 0.05 の場合，QELBM で 1.38×10-4，理論値の圧力勾配 1.28×10-4，U = 0.1 の場合，QELBM で 7.07×10-4，理論値の圧力勾配 5.12×10-4 となり，一様流入速度が小さくなる程，両者の値は近づくことが分かった．この原因としては一様流入速度が大きい場合にはマッハ数が大きくなるため，計算誤差が大きくなることが考えられる．また，U = 0.05 の場合の速度分布や圧力勾配の理論値からの若干のずれは，格子解像度が影響していることが考えられる．以上の結果から非圧縮性流れのようなマッハ数が低い条件では，QELBM は出口・入口を有する流れ場においても妥当な解を与えることが分かった．

15Hx500i 1,3 ,515H=15U=0.1U

255×15×15

２－２ 計算条件
本計算では，図 1 に示すように，一辺の長さが H の 3 次元正方ダクトの流れ場を取り扱う．D3Q13 モデルでは図 2 に示すように静止粒子および格子点を中心とする各平面の対角方向の速度で構成される．境界条件には x = 0.に位置する入口平面（y-z 平面）で速度 U = 0.05, 0.1 の一様流入速度を与え，ダクトの壁面では速度 0の条件を与えた．また，流出境界条件として圧力を 0.9 で固定し，速度を上流側のマクロな速度を外挿入する条件を用いた．格子は格子間隔 1 の等間隔直交格子を用い，格子点数 255×15×15 である．計算を実行したレイノルズ数 Re は特に指定しない限り，Re=100 で行った．なお，本研究で用いた Re は U, L および に基づいている．

yMRT (D3Q13)u y , z  2t
1Mg ,   4    g     2  t4  t

QELBM (D3Q13)

３－２ 速度分布及び圧力勾配計算で得られた結果の信頼性を調べるために，図 3 に U = 0.05
および U = 0.1 の場合の，z/H = 0 における異なる流れ方向位置での x-y 平面の速度分布の結果を示す．また，図には比較のために次式で表される理論値の結果も示す．

(6) f  ,u  f  ,u, 2

Number of Gridz

L=255127Fig. 1 Definition of the computational geometry.

３－３ 計算速度
QELBM の計算速度を調査するために,QELBM および MRT の各モデルの CPU 時間を表 2 で比較する．ここで，計算ステップ数は 1.0×105，U = 0.05 で，計算のプラットフォームはクロック周波数2.26GHz のCPU を持つPC であり，並列計算は行っていない．
表 2 より QELBM の計算時間は MRT の約 60％であり，著者らが最近行った3次元キャビティ流れの結果(5) とほぼ同じ結果が得られた．3 次元キャビティ流れから本研究のダクト流れのプログラ
ム上の変更点は格子点数及び境界条件のみであるため，計算時間Fig. 2 Definition of the velocity of D3Q13 model.
2Copyright © 2011 by JSFM

第 25 回数値流体力学シンポジウム
D12-2にも大きな相違がないものと考えられる．尚，今回用いた計算プ
ログラムは高速化のために特別なチューニングを行っていないので，その点については今後，調査する必要がある．

3４．結論
本研究ではD3Q13QELBM を用いて3 次元正方ダクト流れの数値解析を行い，計算の安定性，計算の精度，計算速度について理論値，MRT と比較した．その結果，QELBM モデルは 3 次元立方キャビティ流れの場合と同様に，入口，出口境界を有する場合でも計算の安定性や計算結果の妥当性を損なうこと無く, 高効率な計算が可能であることが分かった．今回の計算で対象としたレイノルズ数は Re =100 と低いため，今後は QELBM を一様等方性乱流に適用して， QELBM の解，特に低い格子解像度を用いた場合の解が，どの程度高レイノルズ数の乱流渦構造を捉えているかについて調査する予定である．

x/H=5.0x/H=10.0

U=0.05

2.52u/uave 1.5
10.50-1-0.5

00.5

10.5

1y/Hx/H=15.0
Theoryx/H=15.0
Theoryx/H=5.0x/H=10.0

3

U=0.1

32.5

U=0.052.5x/H=15.0
Theoryx/H=5.0x/H=10.0

22u/uave 1.5u/U 1.5

110.5

0.5

0‐1‐0.50y/H0-1-0.5

00.5

1Fig. 4 Mass flow rate profile at different position in flow direction(z/H = 0).

y/Hx/H=5.0
x/H=10.0

Table 2 Comparison of CPU time (Number of time step = 1.0×105)

x/H=15.0
TheoryQELBM (D3Q13)

3U=0.1
2.53MRT (D3Q13)

CPUtime [sec]5.87×10

1.02×104

Ratio

1.00

1.75

2

参考文献
(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.
(5) 安田，里深，南川, "格子ボルツマン法を用いた高レイノルズ数流れの数値計算手法に関する研究," 日本機械学会 2011 年
度年次大会講演論文集, (2011), pp. G050011.

u/U 1.5
10.5
0‐1‐0.50y/H

0.5

1Fig. 3 Velocity profile at different position in flow direction ((z/H = 0).

3Copyright © 2011 by JSFM