Flow simulation in micro channel based on stabilized bubble function element

安定化気泡関数要素によるマイクロチャネル内の流れ解析

Flow simulation in micro channel based on stabilized bubble function element

○ 田井 茂俊（長岡技科大院） 正 倉橋 貴彦（長岡技科大院）

Shigeotoshi Tai, Graduate school of Nagaoka University of Technology, 1603-1, Kamitomioka, Nagaoka, Niigata, Takahiko Kurahashi, Nagaoka University of Technology, 1603-1 Kamitomioka, Nagaoka, Niigata,

Key Words: Microchannel, FEM, Bubble function, SUPG method, droplet generation

1. 序論

近年マイクロ化学チップを利用した液滴生成技術が注目 されている．例えば，医療分野におけるドラッグデリバリー システム(DDS)では病巣への薬剤の標的指向などで薬物を 包んでいるカプセル形成微粒子の単分散化が重要な点とし て挙げられている．

マイクロスケール流路内において生成される液滴サイズ 及び液滴間隔は，流入する二液体と流量条件の組み合わせに より変化する．そこで，事前に予測できるツールの開発を目 指して有限要素解析による検討を行った．流れ解析を行う 際,SUPG 法における安定化パラメータにおいて時間増分量

∆𝑡を考慮した安定化パラメータ^[1][2]と等しくなる様に気泡関 数要素における安定化パラメータを設定し,数値実験による 考察を行う.

本論文では,SUPG 法の安定化パラメータを使用した安定 化気泡関数有限要素法に関する考察をRotating Cone問題に 対して行う．また，その結果を踏まえて，下記に示す支配方 程式に対してRotating Cone問題において使用した安定化パ ラメータを適用し，マイクロチャネル内の流れ解析を行う．

流れ解析では，本研究室で実施されたマイクロ化学チップに よる液滴生成実験（図1）の再現を行う．数値実験ではマイ クロチャネルの流路が円環であると仮定し，ポアズイユ流れ を考えて流入させると共に，マイクロ化学チップに二液体を 流入させる際に使用するシリンジポンプの送液特性を考慮 する．

2. 安定化気泡関数有限要素法による離散化 本検討では流れ場を表現するための支配方程式として式 (1)に示す非圧縮性粘性流体に対するナビエ・ストークスの 運動方程式と式(2)に示す連続式,また,界面位置を表す式(3) の輸送方程式を用いる.式(1)～(3)は総和規約による表示を しており，本研究では平面二次元を対象とする．

𝑉_𝑖̇ + 𝑉_𝑗𝑉_𝑖,𝑗+1

𝜌𝑃_,𝑖− 𝜐(𝑉_𝑖,𝑗+ 𝑉_𝑗,𝑖)

,𝑗= 𝑓_𝑖^𝑉 𝑉𝑖,𝑗 = 0

𝜙̇ + 𝑉_𝑖𝜙_,𝑖= 0

ここで𝑉_𝑖は流速, 𝑃は圧力である.また,νは動粘性係数,ρは

密度,fiVは外力（チップを水平に置いているため,体積力とし ての重力は無視し,表面(界面)張力のみ考慮する.）𝜙は界面 形状を表す指標関数（以後、指標関数）を示す.またマイク ロスケール流路内におけるに流体の表面(界面)張力を表す ため，CSF モデル^[3]を適用し，体積力fiVは式(4)の様に表さ れる.離散化手法としては，SUPG 法における安定化気泡関 数要素を用いた有限要素法を適用し，時間方向の離散化に対 しては θ 法を適用する．式(4)における𝜅(𝑥)は表面(界面)曲 率を表しており,式(5)により与えられる.

𝑓_𝑖^𝑉 = σκ 𝜌_,𝑖

𝜌

̂

𝜌 𝜌

̃

κ(𝑥) = 1

|𝒏|[𝑛_𝑖

|𝒏||𝒏|_,𝑖− 𝑛_𝑖,𝑖]

ここで σは表面(界面)張力係数であり〈𝜌〉,[ρ]及び界面上に おける法線ベクトル niは, 〈𝜌〉 = 0.5 × (𝜌₁+ 𝜌₂), [𝜌] = 𝜌₂− 𝜌₁(𝜌₁< 𝜌₂), ni=𝜙_,𝑖と与えられる.ここにθは0≦𝜙≦1で設定 される値であり,指標関数𝜙の値が0.5≦𝜙あるいは𝜙＜0.5に より流体の種別を分ける．

流速 Vi,圧力 P,指標関数𝜙に対する補間関数を式(6),(7),(8) に示す．気泡関数要素を流速および指標関数の補間に適用し，

圧力の補間には三角形一次要素を用いる．

𝑉_𝑖= 𝑁₁𝑉_𝑖1+ 𝑁₂𝑉_𝑖2+ 𝑁₃𝑉_𝑖3+ 𝑁₄𝑉̃_𝑖4= {𝑁_𝐵}^𝑇{𝑉_𝑖} 𝑃 = 𝑁1𝑃1+ 𝑁2𝑃2+ 𝑁3𝑃3= {𝑁}^𝑇{𝑃}

𝜙 = 𝑁₁𝜙₁+ 𝑁₂𝜙₂+ 𝑁₃𝜙₃+ 𝑁₄𝜙̃₄= {𝑁_𝐵}^𝑇{𝜙}

ここで Nは面積座標(𝜂₁,𝜂2,𝜂3)により表される面積座標を 示しており,𝑁₁=𝜂₁,𝑁₂=𝜂₂,𝑁₃=𝜂₃,𝑁₄=27𝜂₁𝜂₂𝜂₃と表される．気 泡関数要素を用いた流体解析においては，数値粘性が適切に 与えられないことが知られている．運動方程式(式(1))に対し て気泡関数要素を適用した場合，安定化有限要素法における 安定化パラメータは式(9)のように表される．

𝜏_{𝐵𝑢𝑏𝑏𝑙𝑒}^{𝑀𝑜𝑚.} ＝ (∫ 𝑁_{Ω 4}

𝑒 𝑑Ω)²

(𝜐 + 𝜐′) ∫ 𝑁Ω𝑒 4,𝑖𝑁_4,𝑖𝑑Ω𝐴_𝑒

ここにν′は人工粘性を表す係数である．ここで，この安定 化パラメータ𝜏_{𝐵𝑢𝑏𝑏𝑙𝑒}^{𝑀𝑜𝑚.} はSUPG法による安定化パラメータ

𝜏_{𝑆𝑈𝑃𝐺}(式(10))と等価となるように設定する．

𝜏_{𝑆𝑈𝑃𝐺}= ((2 Δ𝑡)

2

+ (2|𝑈|

ℎ )

2

+ (4𝜐 ℎ²)

2

)

−1 2

ここで|𝑈|は要素内の流速の大きさ，hは要素の代表寸法を表 しており, ℎ = √2𝐴_𝑒により計算される．最終的に運動方程式 を解く際，要素重心点において，人工粘性を表す係数𝜐′を 含めて式(11)に示す粘性が与えられることになる．

Distilled water : 10[μl/min]

Ethyl acetate 0.5[μl/min]

Fig1. Experimental result

(1) (2) (3)

(4) (5)

(6) (7) (8)

(9)

(10)

(𝜐 + 𝜐^′) ∫ 𝑁_4,𝑖𝑁_4,𝑖𝑑Ω

Ω_𝑒

= 1

𝐴_𝑒(∫ 𝑁₄𝑑Ω

Ω𝑒

)

2

((2 Δ𝑡)

2

+ (2|𝑈|

ℎ )

2

+ (4𝜐 ℎ²)

2

)

1 2

上記内容は，輸送方程式にも適用でき，輸送方程式に対して も要素重心点において，式(12)に示す数値粘性を与える．

𝜐^′∫ 𝑁4,𝑖𝑁4,𝑖𝑑Ω

Ω_𝑒

= 1

𝐴𝑒(∫ 𝑁4𝑑Ω

Ω_𝑒

)

2

((2 Δ𝑡)

2

+ (2|𝑈|

ℎ )

2

)

1 2

3. 安定化気泡関数有限要素法による流れ場の解析 3.1 Rotating Cone 問題に対する考察

非定常問題の数値解析例として，移流方程式の解の精度を 検証するためにRotating Cone問題を取り扱い，SUPG法に おける安定化パラメータを用いた気泡関数要素の精度を調 べる．空間領域(-1≦x≦1，-1≦y≦1)とし，次式に示される

「コーン形状」(図3)をした変数𝜙に対する空間分布を与え る(式(13),(14))．

𝜙 = {0.5(𝑐𝑜𝑠(2𝜋𝑟) + 1) 𝑟 ≤ 0.5 0 𝑟 > 0.5

𝑟 = √𝑥²+ (𝑦 + 0.5)²

これを初期条件として原点を中心に回転する定常な移流速 度場𝑢 = −𝑦,𝑣 = 𝑥を与えてコーン状の𝜙の分布を移流させ る問題を考える．

図2のようなメッシュを用意し，case1:気泡関数要素を使 用し，安定化パラメータを付与しない場合，case2:気泡関数 要素を使用し，𝜏_{𝑆𝑈𝑃𝐺}= ((²

∆𝑡)²+ (^{2√𝑢2+𝑣2}

ℎ )

2

)

−¹

2(SUPG法にお ける安定化パラメータ)と等価となる安定化パラメータを要 素重心点に入れた場合について比較を行う．

計算条件は，節点数：1,681，要素数：3,200，時間増分量：

∆𝑡=0.001，Time step=6,280，時間方向の離散化にはクランク・

ニコルソン法(𝜙 = 0.5)，空間方向の離散化には気泡関数要 素(𝑁₄の形状関数𝑁₄= 27𝜂1𝜂2𝜂3)を使用したガラーキン法を 適用する．

表1にコーンの各周ごとの中心位置高さを示す．表1 よ り，case 1の条件での解析結果では「コーン形状」が時間経 過とともに減衰していくことが分かる．しかし，case 2の条 件の場合，コーンの先端形状の減衰は小さく，初期の「コー ン形状」に近くなっていることが分かり，高精度に計算が行 えていることを確認できる．

3.2 マイクロチャネル内の流れ解析

本節では，SUPG法による安定化パラメータを用いた気泡 関数要素を使用し，マイクロチャネル内の十字型接合部（オ リフィス部）に着目し，二流体の有限要素解析を行う．

3.2.1 解析モデルと解析条件

有限要素法解析メッシュは図4に示す．二流体は図4に示 す通りに流入し，流量設定は二流体 Q(1)（酢酸エチル）を 0.5μl/min(1.124mm/s)，Q(2)(蒸留水)を10μl/min(22.48mm/s)と する．流量は断面積により割ることで平均流量を算出し，チ ャネル内流路が円環と仮定して，ポアズイユ流れにて流入さ せている．θ法におけるθは0.5を適用し，各流体に対する 物性値及び，その他の計算条件は表2に示す．また，流入流 量の多い蒸留水のシリンジポンプによる送液特性を考え，蒸 留水側の流量にsin波で流量の変化を与える(図4参照)．

3.2.2 解析結果

液滴生成の数値実験結果を図6から図9に示す．この結果 から，液滴が作成される際に実験で見られる二流体の挙動を 再現することができ，流入した酢酸エチルの液滴が生成され，

流れ方向に運ばれていく様子が確認できる．

Total number of nodes 5,105

Total number of elements 9,720

Time increment Δ𝑡 (sec.) 10^-6

Time steps 280,000

Density ρ(1) (kg/mm³) 8.984×10^-7 Kinematic eddy viscosity ν(1) (mm²/sec) 0.557 Density ρ(2) (kg/mm³) 9.982×10^-7 Kinematic eddy viscosity ν(2) (mm²/sec) 1.004 (12)

(13)

(14)

Fig.2 Computational domain Fig.3 Initial value distribution

0 1.00 1.00

1 0.96 1.00

2 0.91 1.00

3 0.83 0.99

4 0.75 0.99

φ(x=0.0,y=-0.5)

Round Coordinates Cone height

in case 1 in case 2 Table 1 Comparison of cone height in cases 1 and 2

Table 2 Computational conditions

Fig. 4 Finite element mesh and definition of boundary condition Fig.4 Finite element mesh and definition of boundary condition

Fig.6 Distribution of index function(T = 0.18sec) (11)

4. 結論

本研究では，CSFモデルを考慮した非圧縮粘性流体に対す るナビエ・ストークス方程式に対してSUPG法における安定 化パラメータを用いた気泡関数要素による有限要素法の有 用性の確認と，マイクロ化学チップ内流路における二相流解 析において液滴が球状になる前の分離を再現した．以下に得 られた知見を示す．

・Rotating Cone問題にてSUPG法における安定化パラメー タを用いた気泡関数要素により安定に，また，精度良く計 算できることが分かった．

・マイクロチャネル内の流れ解析において，流入した酢酸エ チルの液滴が生成される二流体の挙動が再現できた．

現状では分離した酢酸エチルは細い棒状になり進んでいく が，実際の実験では，図1のように球状になることが分かっ ている．今後はこの原因を究明し，現象の再現性を向上させ る．

参考文献

[1]

有限要素法による流れのシミュレーション，丸善出版，

2012

pp58-61

[2] Tayfun E. Tezduyar and Yasuo Osawa, Finite element stabilization parameters computed from element matrices and vectors, Comput. Methods Appl. Mech. Engrg, 190, pp411- 430, 2000

[3] J.U. Brackbill, D.B.Kothe and C.Zemach, A Continuum method for mokiling surface tension , Journal of

Computational physics, 100, pp.335-354, 1992

Fig.7 Distribution of index function(T = 0.19sec)

Fig.8 Distribution of index function (T = 0.20sec)

Fig.9 Distribution of index function (T = 0.21sec)