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

格子ボルツマン法を用いた曲がり管内流れにおける物体の挙動シミュレーション
Lattice Boltzmann Simulation of Motion of a Body in Curved Square Pipe Flows

地村啓，信大院, 〒 380-8553 長野市若里 4-17-1, E-mail：10ta138a@shinshu-u.ac.jp

吉野正人, 信大工, 〒 380-8553 長野市若里 4-17-1, E-mail：masato@shinshu-u.ac.jp

鈴木崇司, 信大院, 〒 380-8553 長野市若里 4-17-1

Hiroshi CHIMURA, Shinshu University, 4-17-1 Wakasato, Nagano 380-8553, Japan
Masato YOSHINO, Shinshu University, 4-17-1 Wakasato, Nagano 380-8553, Japan
Takashi SUZUKI, Shinshu University, 4-17-1 Wakasato, Nagano 380-8553, JapanThe lattice Boltzmann method (LBM) with elastic model is applied to simulations of two-phase ﬂows containing
a deformable body with a viscoelastic membrane. By using this method, we investigate the motion of the body incurved square pipe ﬂows. From the results, it is seen that the body laterally migrates to an equilibrium position
between the centerline of the square pipe and the inner or outer wall. Also, it is found that bodies released fromthe same initial position migrate to diﬀerent equilibrium positions according to the body diameters and Dean
numbers.

1.

はじめに曲がり管内流れは，理学・工学における基本的な流れとして様々な分野で見ることができる．特に近年では，曲がり管内流れの特徴である二次流れに着目し，流体中に含まれる固体粒子の分離や分級が試みられている (1, 2, 3) ．
Ookawara ら (1) は，円弧状に曲がったマイクロチャネル
内における粒子を含む流れの数値計算を行い，粒子の挙動が直径により異なることを報告している．しかしながら，曲がり管内を流れる物体について，その詳細な挙動や平衡位置に関する報告例は少ない．
また，医学の分野でも曲がり管内流れは曲率をもった血管の流れとして見られる．血液中に含まれている赤血球は，粘弾性膜が流体を包む袋状の物体であることから，血流による流体力を受け変形を伴いながら流れる．そのため，血液流れは非常に複雑な流れである．このように，変形する物体を含む流れは，医学のみならず多くの分野で応用されており，変形する物体を含む流れ場の詳細を調べることは非常に重要である．和田ら (4) や Tsubota ら(5)は，ばねモデルによって赤血球を再現し，それらを粒子法を用いて解析することで，赤血球の変形と流れについて報告している．Ramanujan–Pozrikidis ら (6) は，境界要素法を用いて膜の変形と流れの連成解析を行い，せん断流れ場において物体の初期形状が挙動に与える影響について調べている．しかしながら，実際は曲率をもった血管が多いにも関わらず，そのような流路内での変形する物体を含む流れの解析はほとんど報告されていない．
曲がり管内を流れる物体の挙動は，流れ場のレイノルズ数，物体の直径および初期位置，流路の幾何形状など多くのパラメータが影響し，血管のようなスケールの小さな流路においては，実験・観測による解析は困難である．
このことから，物体の詳細な挙動を解析するには数値計算によるアプローチが有効であると考えられる．
二 相 系 格 子 ボ ル ツ マ ン 法 (7) (Lattice Boltzmann
Method，以下 LBM) は，デカルト座標系を用いること
によりアルゴリズムが簡単で，界面を陽に追跡する必要がなく，質量・運動量の保存性に優れているという長所をもつ数値計算法の一つである．このことから，変形す

1

c

c

9

8

c
c
c
c
c c
c
c
c
c
c
c
6

3

1

4

c

11

10

2

5

13

12
7

14

15

Fig. 1: 3D15V model.

る物体を含む流れ場の LBM による解析も近年行われている．Dupin ら (8) は，ラグランジュ的な粒子で構築された膜モデルによる物体を複数個含んだ平行平板間流れの解析を行っている．また，Wu ら (9) は，ばねモデルによって赤血球を再現し，LBM に様々な手法を組み込むことで計算負荷を抑えた手法を提案している．著者らのグループ (10, 11) も質点ばねモデルに基づく弾性力を既存の液液二相系 LBM(12) に組み込んだ固液二相流のための数値計算法を提案した．しかしながら，いずれも平行平板間やダクトなど曲率をもたない流路形状における解析であり，物体の挙動についても定性的な議論にとどまっている．
そこで本研究では，鈴木ら (13) が提案した曲率をもつ壁面のすべりなし境界条件を組み込んだ液液二相系 LBM(12)を用いて，変形を伴いながら流れる物体の曲がり管内における挙動を調べた．

2.

数値計算法

2.1 二相系格子ボルツマン法使用される物理量はすべて，代表長さ L，粒子の代表
速さ c，時間スケール t0 = L/U (U：流れの代表速さ)，基準 order parameter φ0 ，および基準密度 ρ0 を用いて無次元化したものである (7) ．本研究では，血液流れにおけるc 2011 by JSFM
Copyright 第 25 回数値流体力学シンポジウム
E03-5

赤血球と血漿の関係に見られるように，固体と流体の密度差が非常に小さい問題を対象とする．ここでは，既存の液液二相系 LBM(12) を基にした固体を含む二相流解析のための数値計算法について説明する．格子気体モデル
として，本研究では Fig. 1 に示す 3 次元 15 速度モデルを用いた．この速度モデルの粒子速度 ci (i = 1, 2, 3, . . . , 15)は以下のように与えられる．



Fluid A(FA)

Fluid B(FB)

[c1 , c2 , c3 , c4 , c5 , c6 , c7 , c8 , c9 , c10 , c11 , c12 , c13 , c14 , c15 ] =0 1 0 −1 0 0 0 1 −1 −1 1 1 −1 −1 1
0 0 1 0 −1 0 0 1 1 −1 −1 1 1 −1 −1
0 0 0 0 0 1 −1 1 1 1 1 −1 −1 −1 −1



(M)

(1)

2.2 物体のモデル
Fig. 2 に示すように，流体が粘弾性膜によって包み込まれた環状の物体のモデルを考える（以下では，周囲の外部流体を流体相 A，物体内部の流体を流体相 B とする）．
ここで，膜はばね接合した質点によって構成され，流体相 A と流体相 B を隔てる界面に膜を配置する．本モデル
では，質点間のばねの伸縮から質点に作用する弾性力を求め，膜を構成する質点を含む立方体格子を粘弾性をもつ膜とみなし，その立方体格子の頂点である格子点に質点に作用する弾性力を与えた．本研究において物体のモ
デルは，Fig. 3 に示すような Np 個の質点 Pj を頂点とし，ばねによって接合された三角形要素によって構成される三次元構造をもつ．三角形要素による膜モデルの構築法
は Ramanujan–Pozrikidis(6) および Sui ら (14) の方法を参考に次のように行った．まず，(i) 重心を原点とする正八面体を生成し，各頂点に質点を配置する．次に，(ii) 各
質点と質点に結ばれた線分の中心に新たに質点を配置す
る．さらに，(iii) 各面の隣接する質点同士を直線で結ぶことによって，各面が 4 個の三角形要素に分割される．この (ii) および (iii) の工程を n 回繰り返すことにより，各面が 4n 個の三角形要素によって分割された正八面体を構成した．また，このときの要素数 Ne ，質点数 Np および各質点を結ぶ線要素 (ばね) の総数 Ns は以下となる.
Ne
Np
Ns

=

4n × 8,

=

Ne
+ 2,
2

(3)

3
Ne .
2

(4)

=

(2)

R
xj ,
|r j |

yjs =

R
yj ,
|rj |

zjs =

Fig. 2: Model of a body with viscoelastic membrane.

PjkFig. 3: Diagram of triangular mesh (upper) and springmodel (lower).
小になるように決定する．本モデルでは膜の伸縮に対する弾性抵抗を考え，この抵抗力を質点間のばね k の伸縮より求めた．ばね k の弾性係数を Klk ，自然長を lk0 ，伸縮後の長さを lk とすると，膜の伸縮によって蓄えられる
弾性エネルギ El は次式で与えられる．
N

s
1
Klk
2



k=1

R
zj .
|r j |

P j'

lk

El =

lk − lk0
lk0

2
，

(6)

ここで，k 本目のばねの両端の質点を Pj および Pj  とし，その位置ベクトルをそれぞれ r j および r j  とすると，
lk = |rj − rj  | となる.
弾性エネルギ El が最小になるように仮想仕事の原理を用いて質点 Pj の位置ベクトル r j に関して解くと，質点 Pj に作用する弾性力 F j は次式となる．

質点 Pj の位置ベクトルを rj = (xj , yj , zj ) とすると，各質点を位置 (xj , yj , zj ) より球面上の質点位置 (xsj , yjs , zjs )へ以下の式によって投影し，半径 R の球形状の膜モデルを構築した．

xsj ={Membrane.

Fj = −

(5)

∂El
.
∂rj

(7)

2.3 質点の運動膜を構成する質点は上記の弾性力とともに流れ場による流体力を受け，LBM における計算格子内をラグラン
ジュ的に移動する．時間刻み Δt の間に移動する質点 Pjの位置 r j (t + Δt) は次の微分方程式

膜にはたらく弾性力は，ばね接合した質点に作用する弾性力を求めることで表した．本研究では，和田ら (4) やTsubota ら (5) のアイデアを参考に，物体の変形による
弾性エネルギを求め，そのエネルギが最小となるように物体の形状が決定する問題を扱う．各質点に作用する弾性力は仮想仕事の原理に基づき物体の弾性エネルギが最

drj (t)dt
2=

uj (t),(8)

c 2011 by JSFM
Copyright 第 25 回数値流体力学シンポジウム
E03-5

G

H

度分布関数 fi (x, t) および gi (x, t) の時間発展は，LBMを拡張させた Lattice Kinetic スキーム (7) を用いてそれぞれ次式により計算される．

VA
F

E
Δx

Pj
C

D
Δx

A

Δx

φ=

fic

15


=

8

15

fi ,

p=

1
gi ,
3 i=1

u=

15


gi ci .

(14)

i=1

κf
|∇φ|2
6
(15)
+3Ei φciα uα + Ei κf Gαβ (φ)ciα ciβ ,Hi φ + Fi p0 (φ) − κf φ∇2 φ −

9

10

15

H1 = 1, H2 = H3 = H4 = · · · = H15 = 0, ⎪⎪⎪⎪F = −7/3, F = 3E (i = 2, 3, 4, · · · , 15), ⎭
1

i

(16)

i

および

1
− aφ2 ,
1 − bφp0 (φ)=φTGαβ (φ)=

9 ∂φ ∂φ
3 ∂φ ∂φ−δαβ ,2 ∂xα ∂xβ
2 ∂xγ ∂xγ

(17)

(18)

である．式 (15)〜(18) において，α, β, γ = x, y, z（添字 α, β, γ はデカルト座標を表し，総和規約に従う），κf
は界面の厚さを決めるパラメータであり，a, b, T は orderparameter φ の最大値および最小値を決定するパラメー
タ，δαβ はクロネッカーのデルタである．
一方，式 (13) 中の関数 gic は，弾性力を表す外力項を加え次式で与えられる．
39cgi = Ei 3p + 3ciα uα − uα uα + ciα ciβ uα uβ
2
2



∂uβ
3
∂uα
+ Δx
+
ciα ciβ
4
∂xα
∂xβ
 

∂uβ
∂
∂uα
+3Ei ciα
+
μ
Δx
∂xβ
∂xα
∂xβ
(19)
+3Ei δΠ M Fα ciα ,ω(r j , xΨ )F j (t)δ(r j , xΨ ), (10)

ここで，δ(r j , xΨ ) は格子点 Ψ に寄与する質点 Pj を識別する関数であり，本モデルでは，格子点 Ψ を頂点にもつ
隣接する 8 個の立方体格子からなる参照領域内に含まれる質点を考える．つまり，関数 δ(r j , xΨ ) は以下のように定義した(8) ．
⎧⎪⎪⎨ 1 if the particle Pj is in

the reference domain,(13)⎫
E1 = 2/9, E2 = E3 = E4 = · · · = E7 = 1/9,⎪⎪⎪⎪⎬
E = E = E = · · · = E = 1/72,

j=1

⎪
⎪
⎩ 0

gic (x, t),

ここで，for Ψ = A, B, C, . . . , H,(9)
ここで，ω(r j , xΨ ) は Ψ ω(rj , xΨ ) = 1 となる重み係数であり，xΨ は格子点 Ψ の位置ベクトルである．ω(rj , xΨ )は図 4 に示すように，質点 Pj を含む立方体格子 (格子間隔を Δx とする) の格子点 Ψ に対し，対角反対側に位置する格子点と質点 Pj を頂点とする直方体の体積によって決まる（例えば，Ψ = A に対する重み係数 ω(rj , xA )は，図中の灰色部の体積 VA を立方体格子の体積 (Δx)3で割った値に等しい）．なお，質点の位置は一般に格子点と一致しないため，速度の場合と同様に，質点 Pj に作用する弾性力 F j (t) に重み係数をかけて質点 Pj が含まれる立方体格子の 8 個の格子点 Ψ に分配する．したがって，格子点 Ψ における弾性力 F (xΨ , t) は，周囲の質点に作用する弾性力からの寄与の総和として以下の式より求める．

=

gi (x + ci Δx, t + Δt) =式 (12) 中の関数 fic は，次式で与えられる．

Ψδ(r j , xΨ )

(12)

i=1

を 4 次の Runge–Kutta 法を用いて解くことで求めた．ここで，uj (t) は質点 Pj の速度であり，質点が含まれる立方体格子の頂点である格子点 Ψ (= A, B, C, . . . , H) の流速を用いて線形補間で以下の式より求める．
uj (t) =ω(r j , xΨ )u(xΨ , t)

F (xΨ , t) =fic (x, t),

上式において，fic および gic は後述の巨視的変数のみで定義される関数である．二相を表す order parameter φ，
流体の圧力 p および流速 u は，速度分布関数を用いてそれぞれ次のように定義する．

B

Fig. 4: Cubic lattice with sides Δx where particle Pj is
included. The weighting coeﬃcient ω(rj , xA ) is equal
to the shaded volume VA divided by volume (Δx)3 of
the cubic lattice.

Np

fi (x + ci Δx, t + Δt) =(11)otherwise.

ここで，μ は粘度を表し，右辺最終項は弾性力を表す外力項であり，添字の Π は FA，FB，M のいずれかを表す．また，物体の界面張力 κg に関して，弾性力の影響を調べるため κg = 0 にすることで，二相系 LBM に考慮されている界面張力を 0 とした．式 (15)，(18) および (19)2.4 物体を含む二相系 LBM
以下の物理量における添字 FA, FB, M は，それぞれ流体相 A, 流体相 B，膜を表すものとする．本計算では，格子点xおよび時刻 t における速度ci をもつ仮想粒子の速3c 2011 by JSFM
Copyright 第 25 回数値流体力学シンポジウム
E03-5

D

Inflow y'

s
z'
y
x
z

Outflow /W = −0.25)(a) Case 1 (zc0

H

y' s
z'

Flow

W

R

murayamaθ /W = 0)(b) Case 2 (zc0Fig. 5: Computational domain of curved square pipe :overall view (upper) and cross section (lower).
における微分には，テーラー級数展開を用いて導かれる
以下の中心差分近似を使用した．

∂ψ
∂xα

∇2 ψ

3.

15

≈

≈

1 
ciα ψ(x + ci Δx),
10Δx i=1

15(Δx)2

 15


(20)
 /W = 0.25)(c) Case 3 (zc0

ψ(x + ci Δx) − 14ψ(x) . (21)Fig. 6: Motion of a body for various initial positions at
Dn = 5.30; top view (left), right-side view (right). Solidline indicates the trajectory of the body.

i=2

計算結果および考察

3.1 計算条件
Fig. 5 に示すように，矩形断面 (H × W ) をもち，流入
口から流出口までの管軸長さ L の曲がり管内に直径 D の粘弾性皮膜をもった変形する物体（以下では単に物体と呼ぶ）を一つ配置し，流入口と流出口の圧力差 Δp で駆動される流れ場における物体の挙動計算を行った．計算領
域は，H × W = 40Δx × 40Δx, L ∼ 354Δx，曲がり部の曲率半径は R = 40Δx とした．以下では，流入口断面の中心を原点とし管軸に沿う軸を s 軸と定義する．また，y 
軸および z  軸を Fig. 5 のように定義し，θ を曲がり部開始地点から流れ方向に 0◦ ≤ θ ≤ 180◦ の範囲で定義する．
境界条件は，物体がじゅうぶん平衡状態に達した後，時刻 t = 0 で流入口 (s = 0) および流出口 (s = L) に圧力差
Δp = 1.00 × 10−4 を与えた．その他の壁面にはすべりなし境界条件を用いた．物体の直径は D = 12Δx とし，Dとダクト幅 W との比は D/W = 0.30 である．また，物体の密度および粘性係数をそれぞれ ρB = 1, μB = 0.01Δxとし，界面厚さに関するパラメータを κf = 0.01(Δx)2 で一定とした．その他の計算パラメータは，a = 9/49, b =2/21, T = 0.55, φmax = 4.937, φmin = 2.251 とした．曲がり管内流れの無次元数として，レイノルズ数 Re および
曲率比 C をそれぞれ Re = ρF uin Dh /μF , C = Dh /R と定義した．ここで，周囲流体の密度は ρF = 1 とし，uin は各4

時刻における入口断面平均流速を，物体が初期位置より平衡位置に到達するまで時間平均したものとする．また，Dhは水力直径（矩形断面においては Dh = 2HW/(H + W )となる）であり，さらに，Re と曲率比 C の積である Dean
数は Dn = Re C と定義した．

3.2

z  軸上の初期位置より運動する物体の挙動

物体の挙動に関する初期位置の影響を調べるため，物



, zc0
) を sc0 /L = 0.17, yc0/H = 0体の初期位置 (sc0 , yc0
で固定し，zc0 /W =−0.25, 0, 0.25 の 3 ケースに対して計算を行った．それぞれのケースを Case 1,Case 2,Case 3 とする．本計算においては曲率比 C = 1.00 より，Re = Dnとなる. Dn = 5.30 における各初期位置より運動を開始した物体の挙動を Fig. 6 に示す．なお，Fig. 6 の左図は計算領域を上側 (y = +∞) から見た図，右図は計算領域を右側 (x = +∞) から見た図である．図より，Case 1 では，物体は管軸 (y  = z  = 0) と外側管壁 (z  /W = −0.50)
間を流れている様子が見られる．また，Case 2 およびCase 3 では，物体は曲がり部において外側管壁方向へ移動していき，Case 2 は管軸を越えて管軸と外側管壁間を，Case 3 では管軸を越えることなく管軸と内側管壁(z  /W = 0.50) 間をそれぞれ流れている様子が見られる．

c 2011 by JSFM
Copyright 第 25 回数値流体力学シンポジウム
E03-5

θ = 0 θ = 90 θ =180

0.50

zc' W

0.25
0

-0.25
Case 1

-0.50
0

0.2

Case 2

0.4

s L

0.6

Case 30.8

1.0(a) D/Dh = 0.20

(a)

yc' H

0.50

0

(b) D/Dh = 0.30

-0.50
-0.30

0

zc' W

0.30

(b)
Fig. 7: Body trajectories for various initial positions
at Dn = 5.30; ○ and × indicate initial and equilibrium
positions, respectively. (a) Trajectories projected on thez  –s plane. (b) Trajectories projected on the y  –z  plane.

Outer side

Inner side(c) D/Dh = 0.50

Fig. 9: Motion of a body released from the same initial
position on the z  -axis for diﬀerent body diameters at
Dn = 5.30; top view (left), right-side view (right). Solidline indicates the trajectory of the body.

y'sでは外側管壁から内側管壁に向かう二次流れが見られる．
その結果，上下で対称な一対の渦（Dean 渦）が形成されることがわかる．Gossett ら (3) は，曲がり管内を流れる固体粒子（固体粒子の直径 DP とダクト幅 W との比
DP /W = 0.10）について，(i) ダクト中央の内側管壁近傍の固体粒子は，内側管壁から外側管壁に向かう二次流
れの影響を受けて外側管壁方向へと運動し，(ii) 外側管壁
近傍の固体粒子は，外側管壁から内側管壁に向かう二次
流れの影響を受けて，上下のいずれかを経由して内側管壁方向へと運動していくと報告している．本計算結果と比較すると，物体のダクト断面内における移動量は小さ
いが，Case 2 および Case 3 における曲がり部での物体の外側管壁方向への運動は上記 (i) と対応していると考えられる．また，Case 1 での内側管壁方向への運動は上記 (ii) と対応していると考えられる．以上の結果から，曲がり管内の単一物体は固体粒子の場合と同様に Dean 渦の影響を受け，物体は変形しながら流れることがわかった．初期位置において管軸と外側管壁の間に物体を配置したケースでは，内側管壁方向へ流れることがわかった．
また，管軸および管軸と内側管壁の間に物体を配置したケースでは，物体は変形を伴いながら外側管壁方向へ流れる様子が見られた．

z'Fig. 8: Velocity vectors on cross section for θ = 90◦ at
Dn = 5.30.

また，全てのケースで物体は変形しながら流れていることがわかる．さらに，このときの物体の重心の軌跡を Fig. 7に示す．図 (a) は，物体の重心の軌跡を z  –s 平面に投影した図，図 (b) はダクト断面（y  –z  平面）に投影した図である．図 (a) から，各初期位置（図中の○）より運動を開始した物体は，曲がり部 θ = 180◦ を超えた時点で最も外側管壁側へ至り，それ以降の直線部では s 軸と平行に流れることがわかる．また，図 (b) から Case 2 および Case 3 は大きく外側管壁方向へ流れていることがわかる．なお，各ケースにおける物体の平衡位置はそれぞれ(yc /H, zc /W )=(0, −0.24), (0, −0.09), (0, 0.13) である．
本計算結果における物体の挙動は，曲がり部における二次流れの影響によるものと考えられる．また，Dn = 5.30
における物体を含まない場合の θ = 90◦ の断面の流速ベクトル図を Fig. 8 に示す．図よりダクト中央では内側管壁から外側管壁に向かう二次流れが見られ，上下壁近傍

5c 2011 by JSFM
Copyright 第 25 回数値流体力学シンポジウム
E03-5

θ = 0 θ = 90 θ =180

0.50

zc' W

0.25
0
D Dh =0.2-0.25

D Dh =0.3

-0.50
0

0.2

D Dh =0.5(a) Dn = 0.24

0.4

s L

0.6

0.8

1.0

(a)

yc' H

0.50

(b) Dn = 1.46

0

-0.50
-0.30

0

zc' W

0.30

(b)
(c) Dn = 5.30

Fig. 10: Body trajectories for the same initial positionon the z  -axis for diﬀerent body diameters at Dn = 5.30;○ and × indicate initial and equilibrium positions, respectively. (a) Trajectories projected on the z  –s plane.
(b) Trajectories projected on the y  –z  plane.
3.3 物体の挙動におよぼす管径比の影響
物体の挙動におよぼすダクト断面の管径比の影響を調べるため，Dean 数が Dn ∼ 5.30 において，管径比を0.20 ≤ D/Dh ≤ 0.50 の範囲で変化させて計算を行った．
はじめに，計算結果の一例として， Dn = 5.30 におけ/H, zc0/W )=(0.17, 0, 0) より運動る初期位置 (sc0 /L, yc0を開始した. 管径比が D/Dh = 0.20, 0.30, 0.50 の物体の挙動および重心の軌跡をそれぞれ Fig. 9 および Fig.10 に示す．図より，すべてのケースにおいて物体は曲がり部において管軸と外側管壁方向へ流れ，曲がり部 θ = 180◦を超えた時点で最も外側管壁側へ至っている．z  方向の移動に関して，管径比 D/Dh = 0.20, 0.30 のケースはほぼ同じ移動量を示したのに対し，管径比 D/Dh = 0.50 のケースは他のケースに比べて移動量が小さい結果が得られた．このことから，管径比 D/Dh が大きくなるにつれて物体は外側管壁方向へ流れにくくなることがわかった．
Bhagat ら (15) は，曲がり管内を流れる固体粒子に作用する力として，ストークス抵抗を仮定した二次流れによる力 FD ，および主流せん断流れによる力 FＬ を考えた．二つの力 FD および FL は，二次流れによる力 FD が固体粒子の直径 DP に比例する (FD ∝ DP ) のに対して，主流のせん断流れによる力 FL は固体粒子の直径 DP の 4 乗に比例する (FL ∝ DP 4 )．この固体粒子に作用する力 FD は二次流れの方向であり，FL は壁面に向かう方向である．したがって，Dn = 5.30，D/Dh = 0.50 に対する物体の挙動の計算結果は，物体に作用する主流のせん断流れによる力が二次流れによる力より大きい (FD < FL ) ため，二次流
れの方向へ移動しづらく本計算結果のように管軸上を流れ3.4 物体の挙動におよぼす Dean 数の影響
次に，物体の挙動におよぼす Dean 数 Dn の影響を調べ
/H, zc0/W )=(0.17, 0, 0)
た．物体を初期位置 (sc0 /L, yc0に配置し，Dean 数を 0.24 ≤ Dn ≤ 1786 の範囲で変化させて計算を行った．計算結果の一例として，Dn =0.24, 1.46, 5.30, 32.8 における物体の挙動および重心の軌跡をそれぞれ Fig. 11 および Fig. 12 に示す．図より，
Dean 数の違いにより物体が管軸を越えないケース (Dn =0.24, 1.48) と管軸を越えるケース (Dn = 5.30, 32.8) に分類されることがわかる．Dean 数の低い Dn = 0.24, 1.46
では，物体は曲がり部を通過する際に管軸を越えること6

c 2011 by JSFM
Copyright (d) Dn = 32.8

Fig. 11: Motion of a particle released from the same initial position on the z  -axis at various Dean numbers; topview (left), right-side view (right). Solid line indicates
the trajectory of the body.
たと考えられる．すなわち，物体の流れ場に対する追従性
が，その直径により変化することがわかる．なお，管径比D/Dh = 0.20, 0.30, 0.50 における物体の平衡位置はそれ
ぞれ (yc /H, zc /W ) = (0, −0.11), (0, −0.09), (0, −0.016)
である．以上の結果から，いずれの管径比においても物
体は曲がり部において外側管壁方向へ流れる様子が見られ，管径比が大きくなるにつれて物体の平衡位置が +z 
方向に移動することがわかった．

第 25 回数値流体力学シンポジウム
E03-5
θ = 0 θ = 90 θ =180

0.50

zc' W

0.25

Dn =0.24

Dn =5.30

Dn =1.46

Dn =32.8が物体の挙動におよぼす影響について調べるため，物体の直径を変化させて計算を行い，物体の平衡位置について調べた．さらに，Dean 数が 0.24 ≤ Dn ≤ 1786 の範囲において計算を行い，物体の挙動におよぼす Dean 数の影響を調べた．
今後は，物体の形状や個数，変形能の違いが物体の挙動に与える影響について調べる予定である．

0

-0.25
-0.50
0

0.2

0.4

s L

0.6

0.8

謝辞
本研究の一部は，日本学術振興会 学術研究助成基金助
成金（基盤研究 (C) No.23560192）により実施したものである．ここに記して謝意を表します．

1.0

(a)

(1) Ookawara, S., Street, D. and Ogawa, K., “Numerical study on development of particle concentration
proﬁles in a curved microchannel,” Chem. Eng.
Sci., 61 (2006), pp. 3714–3724.

yc' H

0.50

(2) Ookawara, S., Oozeki, N., Ogawa, K., Löb, P. andHessel, V., “Process intensiﬁcation of particle separation by lift force in arc microchannel with bifurcation,” Chem. Eng. Process., 49 (2010), pp.
697–703.

0

-0.50
-0.50

0

zc' W

0.50

(3) Gossett, D.R. and Di Carlo, D., “Particle focusing mechanisms in curving conﬁned ﬂows,” Anal.
Chem., 35 (2009), pp. 8459–8465.

(b)(4) 和田，小林，“膨潤赤血球の体積減少に伴う形状変化の数値シミュレーション,” 日本機械学会論文集 (A編), 69 (2003), pp. 14–21.

Fig. 12: Particle trajectories at various Dean numbers;○ and × indicate initial and equilibrium positions, respectively. (a) Trajectories projected on the z  –s plane.
(b) Trajectories projected on the y  –z  plane.

(5) Tsubota, K. and Wada, S., “Particle method forcomputer simulation of red blood cell motion inblood ﬂow,” Comput. Methods Programs Biomed.,
83 (2006), pp. 139–146.

なく管軸と内側管壁間を流れている．一方，Dean 数の高
い Dn = 5.30, 32.8 では，物体は曲がり部において遠心力の効果が大きく，外側管壁方向へ移動していき，管軸を越えて管軸と外側管壁間を流れている．これは，曲がり部においてダクト中央における外側管壁方向の二次流れの影響と，スリップ速度 (主流速度と物体の速度の差)と流速分布形状による内側管壁方向へのせん断流れの影
響の大小関係が，Dean 数によって変化するためであると考えられる．物体が管軸を越えない Dean 数では，曲がり部においてダクト中央における外側管壁方向への二次流
れの影響よりも，スリップ速度と流速分布形状による内側管壁方向へのせん断流れの影響が大きいため，管軸を越えない．これに対し，物体が管軸を越える Dean 数では，反対に二次流れの影響のほうが大きくなるため，物体は管軸を越えて流れるという結果が得られたと考えられる．なお，Dn = 0.24, 1.46, 5.30, 32.8 における物体の平衡位置はそれぞれ (yc /H, zc /W ) = (0, 0.05), (0, 0.001),(0, −0.09), (0, −0.29) であった．以上の結果から，曲がり管内を流れる物体に関して，Dean 数により平衡位置が異なることがわかった．

4.

(6) Ramanujan, S. and Pozrikidis, C., “Deformationof liquid capsules enclosed by elastic membranes in
simple shear ﬂow: Large deformation and the eﬀectof ﬂuid viscosities,” J. Fluid Mech., 361 (1998), pp.
117–143.
(7) Inamuro, T., “Lattice Boltzmann methods for viscous ﬂuid ﬂows and for two-phase ﬂuid ﬂows,”
Fluid Dyn. Res., 38 (2006), pp. 641–659.
(8) Dupin, M.M., Halliday, I., Care, C.M. and Munn,
L.L., “Lattice Boltzmann modelling of blood cell
dynamics,” Int. J. Comput. Fluid Dyn., 22 (2008),
pp. 481–492.
(9) Wu, J. and Aidun, C.K., “Simulating 3D deformable particle suspensions using lattice Boltzmann method with discrete external boundary
force,” Int. J. Numer. Meth. Fluids, 62 (2009),
pp. 765–783.

おわりに曲率をもつ壁面に対する境界条件を組み込んだ液液二相系 LBM を用いて，正方形断面の曲がり管内を流れる物
体の挙動計算を行った．Dean 数が Dn = 5.30 において，物体の初期位置をダクト中央の内側管壁と外側管壁の間で変化させて計算を行った．また，ダクト断面の管径比7(10) Yoshino, M. and Murayama, T., “A lattice Boltzmann method for a two-phase ﬂow containing solid
bodies with viscoelastic membranes,” Eur. Phys.
J. Special Topics, 171 (2009), pp. 151–157.
(11) Murayama, T., Yoshino, M. and Hirata, T.,“Three-dimensional lattice Boltzmann simulationof two-phase ﬂow containing a deformable bodyc 2011 by JSFM
Copyright 第 25 回数値流体力学シンポジウム
E03-5with a viscoelastic membrane,” Commun. Comput. Phys., 9 (2011), pp. 1397–1413.

(14) Sui, Y., Chew, Y.T., Roy, P. and Low, H.T., “Ahybrid method to study ﬂow-induced deformationof three-dimensional capsules,” J. Comput. Phys.,
227 (2008), pp. 6351–6371.

(12) Inamuro, T., Tomita, R. and Ogino, F., “Lattice
Boltzmann simulations of drop deformation andbreakup in shear ﬂows,” Int. J. Mod. Phys. B,
17 (2003), pp. 21–26.

(15) Bhagat, A.A.S., Kuntaegowdanahalli, S.S. and Papautsky, I., “Continuous particle separation in spiral microchannels using dean ﬂows and diﬀerential
migration,” Lab on a Chip, 8 (2008), pp. 1906–
1914.

(13) 鈴木, 村山, 中西, 吉野, “格子ボルツマン法を用いた曲がり管内における二相流シミュレーション,” 数値
流体力学シンポジウム講演論文集, 24 (2010), E6-5.

8c 2011 by JSFM
Copyright