計算数理工学論文集 Vol. 17 (2017 年 12 月), 論文 No. 05-171215

JASCOME

熱流動解析に対する境界条件強制型

Smoothed Profile-Lattice Boltzmann Method の提案
BOUNDARY CONDITION-ENFORCED SMOOTHED PROFILE-LATTICE BOLTZMANNMETHOD FOR INCOMPRESSIBLE THERMAL FLOWS瀬田剛 1) ，内山 知実 2) ，高野

登 3)

Takeshi SETA, Tomomi UCHIYAMA and Noboru TAKANO
1) 富山大学大学院理工学研究部 (工学)(〒 930-8555

富山市五福 3190,

2) 名古屋大学(〒 464-0814

名古屋市千種区不老町,

E-mail: uchiyama@is.nagoya-u.ac.jp)(〒 930-8555

富山市五福 3190,

E-mail: takano@eng.u-toyama.ac.jp)

未来材料・システム研究所

3) 富山大学大学院理工学研究部 (工学)E-mail: seta@eng.u-toyama.ac.jp)

In this study, the smoothed proﬁle-lattice Boltzmann method is proposed to enforce theDirichlet type boundary conditions for thermal-hyperbolic systems. A smoothed-proﬁlewith a smoothly spreading solid-ﬂuid interface layer identiﬁes the solid domain for non-slip
and non-penetration conditions. To satisfy the Dirichlet boundary condition, the present
scheme computes the interaction force and the heat transfer based on the high-order accurate deﬁnitions of the ﬂuid velocity and the temperature. The analytical and numericalsolutions for the symmetric shear ﬂows demonstrate that the proposed method eﬀectively
reduces the boundary-value deviation. In several two-dimensional numerical investigations including those determining the Couette ﬂows, ﬂow around a circular cylinder, andnatural convection, the numerical results of streamlines, isotherms show good agreementwith those of the previous studies or the existing theoretical work.
Key Words : Smoothed Proﬁle Method, Lattice Boltzmann Method, Dirichlet Boundary
Condition, Thermo-Hydraulics
1. はじめに題を解くことができる (4) ．本論文では，Guo による流速の格子ボルツマン法 (Lattice Boltzmann Method, LBM) は，

定義式 (7) に基づき定式化された Smoothed Proﬁle-Lattice

独立した分布関数の運動により流体挙動が再現され，圧力にBoltzmann Method (SP-LBM) (6) を温度方程式を解析出来

関するポアソン方程式を解く必要がないため，並列効率に優るように改良する．熱源を含む温度の定義式 (8) に基づき，れ，高速計算に適している(1)．直交デカルト格子系で，非圧縮性 Navier-Stokes (NS) 方程式が計算される LBM を用い，ソース項を計算した場合の SP-LBM の境界条件の計算精度への影響を検証する．

複雑な幾何形状を有し，大規模計算が必要な，蓮の葉表面の超撥水作用や，電子機器の冷却現象を解明するためには，埋め込み境界法 (Immersed Boundary Method, IBM) (2, 3) や，Smoothed Proﬁle Method (SPM) (4, 5, 6) の適用が効果的である．Guo によって提案された流速の定義式 (7) に基づき，Immersed Boundary-Lattice Boltzmann Method (IB-LBM)の外力を導出することで，Direct forcing method (2) で観察された固体壁面への流体の漏れの除去に成功している (3) ．

SPM では，任意の界面厚さを有するプロファイルを定義する
ことにより，デカルト格子のみで，流体と構造体との連成問

2. 計算手法
2.1. Lattice Boltzmann Method
構造物と流体との熱輸送を含む相互作用がある非圧縮性
流体に対する連続の式，NS 方程式，温度方程式，∇ · u = 0,∂u1
+ (u · ∇)u = − ∇p + ν∇2 u + F ,∂tρ∂T+ (u · ∇)T = χ∇2 T + Q,
∂t

(1)
(2)
(3)

を，LBM により解析する．ここで ρ は流体の密度，u は流2017 年 9 月 21 日受付，2017 年 11 月 8 日受理

速，T は流体の温度，p は圧力，ν は動粘性係数，χ は熱拡

1.5

Periodic Boundary

1.5
theory○ τ= 1+ τ= 5
τ = 10● τ = 20

T0 > T1

up

Tp = T0

0.5

Ti > To
To

Ti

u/u p

Tp = T1

Ro
Ri

0.5

0

-0.5 1.00-0.5 1.00

up

-1.5

0.98

-1.0

0.96

θ

48

0

49

50

51

50

52

100

150

y

-1.5

200

1.5

Schematic diagrams.

-0.5 1.00

が時間発展することで，流体運動が求められる．式 (1)-(3)-1.0に対応する LBM の発展方程式は,

-1.5

τf

+ δt Fα (x, t),(4)eqgα (x, t) − gα(x, t)+ δt Qα (x, t),
−
τg

51

52

100

y

150

200

theory○ τ= 1+ τ= 5
τ = 10● τ = 200-0.5 1.00

0.98

0.98

-1.0

0.96
48

0

49

50

50

51

52

100

150

y

-1.5

200

0.96
48

0

49

50

51

50

52

100

y

150

200

(b-1) velocity(b-2) temperature(b) present SP-LBM
Fig. 2gα (x + cα δt , t + δt ) = gα (x, t)

50

50

0.5

0

LBM では，離散速度 cα に対応した粒子速度分布関数 fα ，gα

−

49

1.0

T /Tp

u/u p

0.5

fα (x + cα δt , t + δt ) = fα (x, t)

48

1.5
theory○ τ= 1+ τ= 5
τ = 10● τ = 201.0

散率，F は構造物と流体との相互作用力，Q は熱源である．

0.96

0

(a-1) velocity(a-2) temperature(a) conventional SP-LBM (4)(a) symmetric shear flows (b) cylindrical Couette flowfα (x, t) − fαeq (x, t)

0

0.98

-1.0

Periodic Boundary

Fig. 1

theory○ τ= 1+ τ= 5
τ = 10● τ = 20

1.0

T /Tp

Periodic Boundary

−u p

h

Periodic Boundary

1.0

Velocity and temperature proﬁles obtained throughcalculation of the symmetric shear ﬂows.
(5)
粒子境界近傍に格子幅に相当する境界幅 ξ を仮定する (4) ．

で与えられる．緩和時間 τf ，τg は，動粘性係数と熱拡散率に対し ν = c2 δt (τf − 0.5)/3，χ = c2 δt (τg − 0.5)/3 の関係があeq
る．D2Q9 モデルに対する平衡分布関数 fαeq ，gαは，[3cα · u9(cα · u)23u2 ]fαeq = ωα ρ 1 +
+
−
,
c2
2c4
2c2
[
2
2]
3cα · u9(cα · u)3ueqgα= ωα T 1 +
+
− 2 ,
c2
2c4
2c

(6)
(7)

1
で
で定義される．重み係数は ω0 = 49 ，ω1−4 = 19 ，ω5−8 = 36ある．式 (4)，(5) に，Chapman-Enskog 展開を適用すると，

φi (x, t) = s(ai − |x − X i (t)|),
x < −ξ/2,

 0,(
)
1
πx
s(x) =sin ξ + 1
|x| < ξ/2,
2



1.
x > ξ/2.

ここで，添え字 i は粒子番号，ai は粒子半径，X i (t) は時間する計算領域全体に対する粒子のプロファイル φ は，

φ(x, t) =方程式 (3) が導出される．圧力は p = c2 ρ/3 で与えられる．
ため，Guo は流速の定義式と外力項を，

1∑
F
fα cα ,
u = u ∗ + δt ,
ρ=
fα ,
u =
ρ
2
α
α
(
)[
]
1
3(cα − u)9(cα · u)
Fα = ωα ρ 1 −+cα · F ,
2τf
c2
c4
∑

∗

(8)

(12)

t における粒子の中心位置を示す．総数 Np 個の粒子が存在

式 (4) から連続の式 (1) と NS 方程式 (2) が，式 (5) から温度
運動量フラックスに対する外力の非物理的な影響を除去する

(11)

Np
∑

φi (x, t),

(13)

i=1

で与えられる．粒子の速度場 up ，温度場 Tp は，

φ(x, t)up (x, t) =Np∑φi (x, t)V i (t)

i=1

(9)

+

Np
∑

φi (x, t)Ωi (t)(x − X i (t)),

(14)

i=1

のように提案した (7) ．同様に，温度の時間に関する微分項
による温度方程式への非物理的な影響を除去するために温度とソース項は，

(
)
1
Q
T =
Q, (10)gα , T = T + δt , Qα = ωα 1 −
2
2τg
α
∗

∑

のように定義される∗(8)．なお，ρ，u，T 等変数の無次元化

は，文献 (9) に従う．

2.2. Smoothed Proﬁle Method
次式で示されるように，SPM では流体領域において φi = 0，
固体領域において φi = 1 となるプロファイル φi を定義し，

φ(x, t)Tp (x, t) =Np∑φi (x, t)Ti (t),

(15)

i=1

で表される．V i ，Ωi はそれぞれ，時間 t における i 番目の粒子の並進速度と角速度を表す．ここで，式 (8) に，u = upを代入すると，SP-LBM における外力は，

F (x, t) = φ(x, t)2(up (x, t) − u∗ (x, t))
,
δt

(16)

となる．式 (16) を式 (8) に代入すると，u(x, t) = φ(x, t)up (x, t) + (1 − φ(x, t))u∗ (x, t),

(17)

1.01

1.0

1.0

1.00
Analytical solutions (25), (29)Temperature (Conventional)○ Temperature (Present)+ Velocity (Conventional)● Velocity (Present)

0.98

0.8
0.9

0.6

u/up

Tj /Tp, uj /up
0
0

0.99

58

59

60

61

0.4

Present (ξi = 2δx)
Present (ξi = 0)
+ Hu et al. (ξi = 2δx)
● Hu et al. (ξi = 0)
Jafari et al. (ξi = 2δx)
○ Jafari et al. (ξi = 0)

0.97

0.2

0.96
0.95
0.5

0
1

5

10

20

30

10

τ

20

30

(a) boundary-value deviation

40

x

50

60

70

(a) velocity

0.8

0.6

Analytical solutions (25), (26), (29)Temperature (Conventional)○ Temperature (Present)+ Velocity (Conventional)● Velocity (Present)

1.0

0.9

0.6

0.08

0.4

1.0

0.8

T /Tp

(Tj -Tj -2 )/Tp , (uj -uj -2 )/u p
0
0
0
0

1.0

58

59

60

61

0.4

Present (ξi = 2δx)
Present (ξi = 0)
+ Hu et al. (ξi = 2δx)
● Hu et al. (ξi = 0)
Jafari et al. (ξi = 2δx)
○ Jafari et al. (ξi = 0)

0.04

0.2

0.5

1

0.2

2

0

0
0.5

1

5

10

20

30

τ

10

(b) boundary slip
Fig. 3

20

30

40

x

50

60

70

(b) temperature

Comparison of analytical and numerical solutions de-

Fig. 4

Tangential velocity and temperature proﬁles as cal-pending on the relaxation time.

culated by the SP-LBMs in cylindrical Couette ﬂows.

が得られる．流体領域において φ = 0 であるから，式 (16) よで与えられ，流体の密度，流速，温度は次式で求められる．

り，外力は F = 0 となる．式 (8) に F = 0 を代入すると，u∗ρ=は流体の速度 uf と一致する．よって，式 (17) は，Jafari に

∑

fα ,

α

よって示された SPM における流速の定義式 (4) と一致する．

u=1∑cα fα ,
ρ α

T =

∑

gα .

(24)

α

Hu によって提案された SP-LBM では，本モデルと同様に，
u = φup + (1 − φ)uf .

(18)

流速，外力項，温度，ソース項に対し，式 (8)，(9)，(10) が用いられるが，外力と熱源は，Jafari によって提案された式式 (18) の両辺の発散を取り，∇ · uf = 0 とすれば，∇ · u = ∇φ · (up − uf ),(21)，(22) が用いられる (5) ．式 (12) の φi に対し，様々なプ
ロファイルが提案されており，文献 (8) では一般的な tanh 関(19)数を用いたが，本研究では Jafari のモデル (4) と同様に sinが得られる．境界において up と uf が一致し，流体の固体関数を用いる．

壁への漏れがなければ，u に対する連続の式 (1) が満足される．同様にして，熱源は，3. 計算結果∗

Q(x, t) = φ(x, t)2(Tp (x, t) − T (x, t)),δt

3.1. Symmetric shear ﬂows

(20)

Fig. 1(a) に示す symmetric shear ﬂows に対し，200 × 200

によって求められる．本手法では，式 (16)，(20) から得られの格子点を用い，y = 50δx ，y = 150δx の位置に，SPM によ

た F ，Q を式 (8)，(9)，(10) に代入することで，熱輸送を伴り，up = 0.01，up = −0.01 の流速を，Tp = 1，Tp = −1 の温う流体運動が計算される．

度を，それぞれ，設定する．ここで，δx = cδt であり，全ての

なお，Jafari は，粒子を剛体球と考える Direct forcing methodに基づき，相互作用力 F を，

F (x, t) = φ(x, t)up (x, t) − u(x, t)
,
δt

(21)

Tp (x, t) − T (x, t).
δt

3cα · F,c2

y = 150δx の位置でのみ，φ = 1 となり，それ以外の格子点では φ = 0 となる．Fig. 2(a) に Jafari によって提案された従来型の SP-LBM による計算結果を，Fig. 2(b) に，本手法による計算結果を示す．SP-LBM においても，IB-LBM と同(22)

様に，緩和時間 τ が増加すると，境界近傍で流速と温度に歪が生じている．本手法では，式 (16)，(20) により強制的に境

この場合，式 (4)，(5) の外力項とソース項は，

Fα = ωα ρ

計算で，δx = δt = c = 1 とする．後述する解析解の導出手順の単純化のため，界面幅を ξ = 0 とする．つまり，y = 50δx ，で与えた．同様に，ソース項 Q は次式で与えることが出来る．

Q(x, t) = φ(x, t)(2)

界値が設定されるため，Dirichlet 条件が常に満足されること

Qα = ωα Q,(23)が，Fig. 2(b) から分かる．

0.2

Error

0.1

Slope = -1(a-1) Re = 20(a-2) Re = 40(a) IB-LBM (Direct forcing method) (2)

0.01

0.001
0.0003Present (ξi = 2δx)
Present (ξi = 0)
+ Hu et al. (ξi = 2δx)
● Hu et al. (ξi = 0)
Jafari et al. (ξi = 2δx)
○ Jafari et al. (ξi = 0)
50

100

200

400

(b-1) Re = 20(b-2) Re = 40(b) IB-LBM (Implicit correction method) (3)Number of grid points(a) error in velocity
0.2

Error

0.1

Slope = -1(c-1) Re = 20(c-2) Re = 40(c) SP-LBM (Jafari et al.) (4)

0.01

0.001
0.0003Present (ξi = 2δx)
Present (ξi = 0)
+ Hu et al. (ξi = 2δx)
● Hu et al. (ξi = 0)
Jafari et al. (ξi = 2δx)
○ Jafari et al. (ξi = 0)
50

100

200

(d-1) Re = 20(d-2) Re = 40(d) SP-LBM (Hu et al.) (5)400

Number of grid points(b) error in temperature
Fig. 5

Relative error versus number of grid points in thecalculation of the cylindrical Couette ﬂows.
文献 (8) の手法を用い，Fig. 2 で観察された境界近傍の歪に対する解析解を導出する．境界点 j0 における値 ϕj0 と，境界点から 2 格子分離れた点 j0 − 2 における値 ϕj0 −2 に関する

Jarari の手法の解析解は，以下のように与えられる．
ϕj0(2τ − 1)B(AC + 1) + 4τ 2 − 5τ + 3=,
ϕp4τ 2 − 3τ + 2
ϕj0 −2= BC.
ϕp

20 and 40. The thick line indicates the zero streamline. The

内円の半径は Ri = 40δx ，外円の半径は Ro = 80δx ，内円の(26)

温度は Ti = 1，外円の温度は To = 0 とする．外円は静止しており，内円は一定の速度 uθp = 0.01c で回転している．緩和
時間は τf = τg = 0.75 とする．厳密解 ûθ ，T̂ は，

4τ − 3τ + 2
3
,
B= 2
,
4τ 2 − 3τ + 5
4τ − 3τ + 5
48
∑
1
1
C=
,
(k + 1 − Ak) (k − A(k − 1))

A=

(27)

ûθ (R) = uθpT̂ (R) =である．本手法に対する解析解は，

ϕj0 −2
= BC,
ϕp

(29)

となる．C は式 (28) で与えられ，A，B は以下のようになる．

B=

68τ 2 − 8τ + 11

R/Ro − Ro /R,
Ri /Ro − Ro /Ri

(31)

(28)

k=1

8τ 2 − 8τ + 5,
8τ 2 − 8τ + 11

Streamlines for ﬂow past a circular cylinder at Re =

(25)

2

A=

Fig. 6gray line indicates the surface of the circular cylinder.

ϕ は流速または温度を表す．ここで，

ϕj0
= 1,
ϕp

(e-1) Re = 20(e-2) Re = 40(e) present SP-LBM

To log(R/Ri ) − Ti log(R/Ro ),log(Ro /Ri )(32)で与えられる．Jafari の手法 (4) ，Hu の手法 (5) ，本手法に対する y = 100δx の断面における流速分布と温度分布を Fig. 4に示す．仮想的な固液界面幅 ξ = 2δx を設定した場合，温度分布，速度分布共に，ξ = 0 の場合よりも厳密解からのずれ

.

(30)

が大きくなる．Fig. 5 に計算領域の一辺の格子点数の数と誤Fig. 3 に，ϕj0 と ϕj0 − ϕj0 −2 に関する数値解と解析解と差との関係を示す．ξ = 2δx の場合，固液界面における φ のを示す．解析解と数値解とは良く一致しており，緩和時間がプロファイルの影響により，本手法の計算精度は向上しない．

大きくなると，境界近傍の歪が大きくなることが分かる．本

固液界面幅が ξ = 0 であれば，Dirichlet 境界条件が厳密に満手法にでは，緩和時間に関わらず，境界値が設定値と常に一足される本手法の誤差が，Fig. 5 で最も小さくなっている．

致することが解析解と数値解から証明された．

3.3. Flow past a circular cylinder

3.2. Cylindrical Couette ﬂows

レイノルズ数 Re = 20 と 40 における円柱周り流れの計

Fig. 1(b) に示す cylindrical Couette ﬂows の計算を行い，算 を 行 う．緩 和 時 間 は τf = 0.65 に 設 定 し ，円 柱 の 直 径 は

SP-LBM の解析精度を検証する．計算領域は 200δx × 200δx ，

D = 100δx ，代表速さは u0 = Reν/D で与える．格子点数

Table 1

1.5○ Dennis and Chang (Re = 40)● Dennis and Chang (Re = 20)
Present (Re = 40)
Present (Re = 20)

1.0

Comparison of drag coeﬃcient Cd , wake length

2L/D, separation angle θs , pressure coeﬃcient at the front
stagnation point Cp (π), and pressure coeﬃcient at the rear

0.5

Cp

stagnation point Cp (0) as calculated by the proposed SP0

LBM and in previous studies.

-0.5
-1.0
-1.5

0

30

60

120

90

150

180

θ

Fig. 7

Pressure distribution on the cylinder.

40D × 40D の計算領域の (16D, 20D) の位置に円柱を設置
し，密度は ρ = 1.0 とする．数値的安定性のため，ξ = 2δxとした SPM により円柱の境界条件を設定する．外周の 4 辺における fα の境界条件は，密度が ρ = 1.0，流速が u0 に対応した平衡分布関数 fαeq により与える．Fig. 6 に，流れ関数
分布を示す．Fig. 6(a)，(c)，(d) に示されるように，Direct

forcing method に基づき，外力を求める手法では，流速の漏れが発生している．Fig. 6(b) に示される Implicit correction

Re References
Cd 2L/D θs Cp (π) −Cp (0)(10)
20 FDM
2.045 1.88 43.7 1.269 0.589
SAM (11)2.053 1.786 43.37 1.274 0.582
LBM (12)2.152 1.842 42.96 1.233 0.563
IB-LBM (3) 2.103 1.900 40.89 1.220 0.580
SP-LBM (4) 2.092 1.900 41.62 1.250 0.587
SP-LBM (5) 2.091 1.900 41.55 1.250 0.588
Present
2.095 1.920 41.64 1.250 0.587
40 FDM (10)1.522 4.69 53.8 1.144 0.509(11)SAM
1.550 4.357 53.34 1.117 0.554
LBM (12)1.499 4.490 52.84 1.133 0.487(3)
IB-LBM
1.566 4.600 50.70 1.105 0.505
SP-LBM (4) 1.562 4.640 51.23 1.133 0.512
SP-LBM (5) 1.562 4.640 51.18 1.133 0.512
Present
1.565 4.660 51.29 1.135 0.512

method の結果と同様に，定義式 (8) と，式 (9) の外力項を用いる本手法では，流速の漏れが除去されることが Fig. 6(e)から分かる．
本手法に対する抗力係数 Cd と圧力係数 Cp とを求める．

ps − p0
Cp =.
(1/2)ρu20

Fd
Cd =
,
(1/2)ρu20 D(33)

抗力 Fd は，運動量保存則に基づき，F の x 成分 Fx から，∫
Fd = −(34)
Fx (x, t)dx,∀pのように求められる．円柱表面上の圧力 ps は IBM で用いられる補間関数を用い，境界から 3δx 離れた点で評価され
2

る．一様流の圧力は p0 = c ρ/3 で与えられる．Fig. 7 において，本手法で得られた圧力係数 Cp の分布は参照解 (10) と間を τf = 3ν + 0.5 から導出する．浮力は，ブジネスク近似

G = βG(T − Tm )j により与える．β は体積膨張率，G は重力加速度の大きさ，j は重力加速度の向きと逆向きの単位ベ
クトル，Tm は基準温度を表し，Tm = 0.5 とする．体積膨張率は，β = Raνχ/h3 で与える．本計算により得られた流れ関数および温度分布は，参照解 (13, 14) と良い一致を示しており，本 SP-LBM は自然対流解析に有効であることが明らかになった．本手法から求められた円柱上の平均ヌセルト数

N̄u と流れ関数の最大値 ψmax とが，参照解 (13, 14) と良い一致を示すことが Table 2 から分かる．

4. おわりに
SP-LBM において外力と熱源の導出方法を修正することに一致しており，圧力分布も正確に計算出来ることが分かる．

より，Dirichlet 境界条件の計算精度を向上させた．Symmetric

Table 1 にまとめられた抗力係数 Cd ， 後流長さ L，剥離点

shear ﬂows の計算において，本 SP-LBM においても，緩和θs ，圧力係数 Cp が参照解 (10, 11, 12) と良い一致を示した．

時間が大きいと境界近傍に流速や温度の歪が発生するが，境

なお，Table 1 の IB-LBM (3) ，SP-LBM (4) ，SP-LBM (5) の界値は設定値と常に一致することが解析解と数値解から証

値は，文献値ではなく，著者らにより実施された Fig. 6 の計明された．円柱周り流れと自然対流の計算において，固体壁

算結果から導出された値である．

への流体の漏れが除去され，流れ関数分布，温度分布とも適3.4. Natural convection

切に計算されたことから，本手法の熱流動解析への有効性が最後に，円柱周りの自然対流解析を行い，ヌセルト数，最確認された．本論文では単純な定常問題のみが検証されて大流れ関数，流速分布，温度分布について検証する．温度おり，今後，流れや温度の発達過程を含む非定常問題への本

To = 0 の正方形の領域の中心に温度 Ti = 1 の円柱を置く．格

手法の有効性を検証する必要がある．本研究は JSPS 科研費

子点数は 221 × 221，ξ = 0 とする．プラントル数を Pr = 0.71

16K06070 の助成を受けたものである．ここに記して謝意をとし，レイリー数 Ra と，計算領域の幅 L と半径 R との比表す．

R/L を変化させた場合の温度分布と流速分布を Fig.8 に示す．Fig.8 の破線で示すように，正方形領域の端から幅 10δx

参考文献

までの固体領域を SPM により設定する．これにより，代表
長さは h = 200δx となる．緩和時間を τg = 0.6 とし，熱拡散(1) Chen, S.，Doolen, G. D.： Lattice Boltzmann method

率 χ を求めた後，ν = χPr から動粘性係数を求め，緩和時

for ﬂuid ﬂows， Annu. Rev. Fluid Mech.， 30(1998)，

Table 2

Surface-averaged Nusselt number and maximum

absolute value of stream function.

Ra R/L0.1
104(a-1) R/L = 0.1(a-2) R/L = 0.2(a) Ra = 104

0.2(a-3) R/L = 0.3
0.3
0.1
105

0.2
0.3

(b-1) R/L = 0.1(b-2) R/L = 0.2(b) Ra = 105(b-3) R/L = 0.3

0.1
106

0.2
0.3

N̄u
ψmax
N̄u
ψmax
N̄u
ψmax
N̄u
ψmax
N̄u
ψmax
N̄u
ψmax
N̄u
ψmax
N̄u
ψmax
N̄u
ψmax

Present SIMPLE (13) DQM (14)
2.019
2.071
2.08
1.792
1.73
1.71
3.234
3.331
3.24
1.016
1.02
0.97
5.601
5.826
5.40
0.510
0.50
0.49
3.672
3.825
3.79
10.16
10.15
9.93
4.918
5.08
4.86
8.391
8.38
8.10
6.490
6.212
6.21
5.180
5.10
5.10
5.887
6.107
6.11
20.81
25.35
20.98
8.860
9.374
8.90
23.95
24.07
24.13
12.54
11.62
12.00
20.47
21.30
20.46

(7) Guo, Z., Zheng, C., Shi, B.： Discrete lattice eﬀects on(c-1) R/L = 0.1(c-2) R/L = 0.2(c) Ra = 106(c-3) R/L = 0.3

the forcing term in the lattice Boltzmann method, Phys.
Rev. E, 65(2002), pp. 046308．

Streamlines and isotherm contours for diﬀerent val-

(8) Seta, T., Hayashi, K., Tomiyama, A.：Analytical andues of R/L at Ra = 104 , Ra = 105 , and Ra = 106 . Stream-

numerical studies of the boundary slip in the immersedlines and isotherm contours are shown on the left- and right-

boundary-thermal lattice Boltzmann method， Int. J.

hand sides, respectively.

Numer. Methods Fluids，(2017)，DOI: 10.1002/ﬂd.4462

Fig. 8(published online)．
pp. 329–364．

(9) Inamuro, T., Yoshino, M., Inoue, H., Mizuno, R.,(2) Feng, Z.-G., Michaelides, E. E.： Proteus: a direct forc-

Ogino, F. ： A lattice Boltzmann method for a bi-

ing method in the simulations of particulate ﬂows, J.

nary miscible ﬂuid mixture and its application to a

Comput. Phys., 202(2005), pp. 20–51．

heat-transfer problem， J. Comput. Phys., 179(2002)，(3) Wu, J.，Shu, C.：Implicit velocity correction-based im-

pp. 201–215．

mersed boundary-lattice Boltzmann method and its ap-(10) Dennis S. C. R, Chang G. Z.： Numerical solutions forplications, J. Comput. Phys., 228(2009), pp. 1963–

steady ﬂow past a circular cylinder at Reynolds numbers

1979．

up to 100， J. Fluid Mech.， 42(1970)， pp. 471–489．

(4) Jafari, S., Yamamoto, R., Rahnama, M.： Lattice-(11) Nieuwstadt F., Keller, H. B.：Viscous ﬂow past circular

Boltzmann method combined with smoothed-proﬁlecylinders， Comput. Fluids， 1(1973)， pp. 59–71．

method for particulate suspensions， Phys. Rev. E，(12) He X, Doolen G. D.： Lattice Boltzmann method on

83(2011)， pp. 026702．

curvilinear coordinates system: ﬂow around a circular(5) Hu, Y, Li, D., Shu, S., Niu, X,：An eﬃcient smoothedcylinder， J. Comput. Phys., 134(1997)， pp. 306–15．

proﬁle-lattice Boltzmann method for the simulation of(13) Moukalled, F.，Acharya, S.： Natural convection inforced and natural convection ﬂows in complex geome-the annulus between concentric horizontal circulartries， Int. Commun. Heat Mass， 68(2015)， pp. 188–and square cylinders， J. Thermophysics Heat Tr.，

199．

10(1996)， pp. 524–531．

(6) Seta, T., Uchiyama, T., Takano, N.： Smoothed proﬁle-(14) Shu, C.，Zhu, Y. D.： Eﬃcient computation of natural

lattice Boltzmann method for non-penetration and wet-

convection in a concentric annulus between an outerting boundary conditions in two and three dimensions，square cylinder and an inner circular cylinder， Int. J.

Comput. Fluids， 159(2017)， pp. 64–80．

Numer. Methods Fluids， 38(2002)， pp. 429–445．