「気液自由界面を用いたマイクロミキサー」の性能評価に関するPhoenicsの応用事例

ナノ・マイクロ応用理工学研究室（小野研究室）所属
芝浦工業大学大学院博士（後期）課程 2年

山田 崇

緒言
• マイクロミキシングとは何か• マイクロ空間内において，複数の流体（液体，気体）を効率よく混合させる．
• ダウンサイジングの代表的活用例の一つ．（例；μ-TAS，有機化合，DNA検出）• マイクロ空間とは代表長さ1～1000マイクロメートルの範囲内にある空間．
レイノルズ数Re（空間内の流体の慣性力-粘性力の比）から

ul
Re 


(1)

但し，ρ；流体の密度，μ；流体の粘度，u；流体の速度l；空間の代表長さ
例；純水（流速1mm/s，空間の代表長さ１マイクロメートル）
流体の密度 ρ=0.998×103 (kg/m3) ,粘度 μ=1.002×10-3(Pa・s)

 Re  0.001 （層流，流れに乱れを起こすのが難しい）
この空間に生じる物理現象 → 慣性力よりも粘性力の支配が強い傾向を示す．

• ダウンサイジングとは何か• たとえば，空間内で二種類の流体が接触しているとき，その二流体の接触界面の表面積-界面を形成する流体の体積比は1/L（長さ）の次元で表される．
• 空間内の流体の代表長さは，空間の代表長さのオーダーを超えることはない．
→ 空間内の代表長さが１マイクロメートルのとき，二流体の表面積-体積比は1,000,000以上となる．（また，界面張力/慣性力，伝導伝熱/対流伝熱，拡散/対流などの比が桁違いに大きくなる．）→接触界面を通じた熱・物質移動，化学反応がマクロ空間内に比べ著しく向上する．
これにより，効果的な熱交換，物質の混合・分離，精密な化学反応の制御が可能となる．

流体１流体２流体２(a)(b)

Fig. 1 空間内の二液の分布例 ； (a) 連続相，(b) 分散相

流体１• マイクロミキサーとは何か
マイクロ空間内での流体運動・現象の一般的な特性• 慣性力よりも粘性力，対流よりも拡散の影響が支配的（流れに乱れが起きにくい）• 「表面積-体積比」の増加による二流体間の接触界面での熱・物質移動速度の高速化
 マイクロミキサーとは，以上のマイクロ空間内での流体の特性を考慮して，マイクロスケールの管径（代表長さ）の直線（曲線）流路内に二流体を，主に管内で連続相を形成するように流入させ，その接触界面上での逐次混合により，効率的な混合を達成させる装置．
として、一般に定義されている．
流体１
接触界面

流体２

Fig. 2 遂次混合のイメージ（左から右に混合が進展）• マイクロミキサーの問題点① １台当たりの生成量が少ない → 多数のミキサーを接続して生成量を確保すること（ナンバリング-アップ）が大前提．
② 代表長さ（管径）の極端な短縮 → 微細管内（マイクロミキサー内部）の圧力損失の大幅上昇• ミキサー内部の破損頻発 → メンテナンス困難• ナンバリングアップ時の膨大な圧力損失 → 装置大型化

ナンバリング-アップ+代表長さ（管径）の大幅な短縮 → 非現実的解決策• １台当たりの（ある程度の）生成量の確保 → （マイクロ空間としては）大きな代表長さを持つ空間内に微小な代表長さを持つ複数流体の層（=流体セグメント）を多数配置• マイクロミキサー内の流量増加 → 高Re数時に発生する流路断面内の二次流れの活用• ミキサー内部の圧力損失が低く，かつ高効率なマイクロミキサーの開発
気液自由界面（気泡表面）の活用

気液自由界面の活用• 気液自由界面とはある均一な液体が他の均一な気相と接している境界を指す．界面上には界面（表面）張力が働く．液体と接している気体（気泡）では，表面張力と気相-液相間の圧力差によって外形（曲率半径）が定まる．

2p  pliq  pgas R(2)Δp；ラプラス圧，pliq；液相側圧力，pgas；気相側圧力，γ；界面（表面）張力，R；界面の曲率半径
液相

気液自由界面気相0.1mm0.025mmFig. 3 当研究室既往研究における微細管内での気液自由界面形成例

気液自由界面活用の利点• （形状は限られるが）マイクロ空間内に微小構造物を構築することが比較的容易．
• 気液自由界面表面でのせん断応力がほぼゼロ．
• 界面上の流体の挙動は，周囲の流体の濃度差，温度差，密度，粘性の微小な差異により，大きな影響を受ける（界面不安定性）．
これらを積極的に用いることにより，ミキサー内部の圧力損失が低く，かつ高効率なマイクロ混合デバイスが開発できる可能性がある．
当研究室が長年扱ってきたテーマ
対流拡散の利用
当研究室では，気液自由界面活用に加えて，マイクロミキサーにおいても対流拡散による混合速度向上を期待できる混合デバイスを考案し，研究を継続してきた．
•三次元非定常移流拡散方程式

  2C  2C  2C C uC vC wC D 2  2  2     (3)
tx
yz
yz  x
移流項活用 → 混合速度高速化

どのようにして，対流拡散による混合速度向上を目指したか？

→ 局所的に速い対流を起こすことにより，移流項活用が可能となる．

当研究室での研究事例

界面振動型マイクロミキサー

マランゴニ対流を用いた
マイクロミキサー（1）マランゴニ対流を用いたマイクロミキサー（２）• これらのミキサー開発に関して以上の３例の混合デバイスは微細管内の非定常現象（界面振動やマランゴニ効果）を利用→ 混合性能やデバイスの機能の予測や混合原理の究明

CFDによる解析が有効

Phoenicsを用いたCFD解析を試みた．

熱流体解析ソフト「PHOENICS」• 当研究室での使用実績から，「PHOENICS」使用のメリットとして（１）計算モデルの作成が容易（２）流れ場の計算にSIMPLE法が用いられることにより，非定常計算においても時間間隔（時間ステップ）をあまり小さくせず，計算メッシュを余り細かく設定せずとも，比較的正確な計算結果（現象の傾向等）を得ることが出来る
また，計算時間を短縮できる．
（３）解析領域内の任意のOBJECTを解析計算内において，比較的簡単に移動境界解析を実現できる「MOFOR」機能を搭載している．
（４）グラフィック機能が充実しているため，解析結果を容易に可視化できる．
また，計算結果を動画として保存することが可能．
（５）CADなどで作成したDATAをSTLファイル化して解析領域内のOBJECTの形状として取り込むことが可能．
（６）Informにより，解析計算のカスタム化が容易 ．
以上の特色を活かした解析モデルを構築する．

• 「PHOENICS」による基本的な解析手法• 実験時においては…
気相内の圧力管理 → 界面形状決定（CFDで界面形状を解く蓋然性があまり無い）• 気液自由界面もひとつの「オブジェクト」と捉える（境界条件に注意）• 解析すべき領域となるのは「液相側」のみ• 解析モデルのベース
単相流による解析（tutorial 10 「濃度拡散の非定常解析」を参考）+• STLファイルによるCADデータ読み込み（事例１）• MOFORの活用（事例１）• INFORMを用いた特定領域へのソース項（速度境界）の追加（事例２，３）

それぞれの応用事例について，簡便な解析モデルを作成し，計算を実施した．

事例１；界面振動型マイクロミキサーモデル作成と計算の流れ
幾何的な計算モデル作成

境界条件（流入・流出速度，壁面条件）設定
従属変数（拡散係数）設定

計算ステップ数（非定常計算） ，繰り返し計算回数設定
MOFOR機能設定
初期条件作成，読み込み
計算・解析

事例１；界面振動型マイクロミキサー
計算モデル概形
解析流体の物性

解析領域
x 方向成分；7.241[mm](2液合流部から混合流路の長さ；7[mm])y 方向成分；0.3[mm]z 方向成分；0.1 [mm]IMAT67 純水 at 20℃
Bubble 20.141mm

m1m
0.

Mixture Channel

0.1mm

WALL 3

WALL 10.1mm

0.
1mm

WALL 20.1mm

INLET1（濃度1） 0.1mm

ＹINLET2（濃度0）Ｘ

Fig. 5 解析モデルの概形

Bubble 1

事例１；界面振動型マイクロミキサー（１）計算モデルの作成• 計算モデルは２次元モデルとして作成（理由；非定常計算時における計算負荷低減）• 計算メッシュ：x（水平）方向；150，y（鉛直）方向；24，z（奥行）方向；1• X方向では，流体の流入部と気泡運動部に関しては格子間隔を小さく，それ以外の場所では格子間隔を大きく取っている．
• 計算モデル中，Y字型合流部，気泡オブジェクト等の形状データはCADにて製作
Bubble2
INLET2（濃度1）
WALL2（CAD)

WALL4（BROCAGE)WALL3OUTLET（CAD)
WALL5（BROCAGE)

WALL１（CAD)
INLET1(濃度0)Bubble1 Bubble2

INLET2（濃度1）
WALL2（CAD)(a)
WALL4（BROCAGE)WALL3OUTLET（CAD)
WALL１（CAD)
INLET1(濃度0)

WALL5（BROCAGE)
Bubble1(b)Fig. 6 解析モデルの基本構成と3次元CADによる形状データ作成例，(a)気泡対向配置モデル(b)気泡交互配置モデル(c) 3次元CADデータ例(c)事例１；界面振動型マイクロミキサー（２）境界条件（流入・流出速度，壁面条件）• 各Inletには，実験時の体積流量である8.3×10-12[m3/s]を設定• 壁面に「198 Solid with smooth-wall friction.」を
気泡部には「199 Solid allowing fluid-slips at walls.」を用い，境界条件を構築する．
WALL Surface (198 Solid with smooth-wall friction. )Gas-Liquid Free Interface (199 Solid allowing fluid-slips at walls. )Bubble ( Solid allowing fluid-slips at walls.)
Flow Property

Fig. 7 気泡オブジェクト周辺部
での境界条件事例１；界面振動型マイクロミキサー（３）従属変数（及び交換係数）設定• 従属変数φ；「PHOENICS」中で「一般的な変数」として予約されている「C1」を用いる．
• 交換係数Γφ設定のため，拡散係数を実験時の値に相当する，3.0×10-10[m2/s]に設定．

Fig. 8 PHOENICS内で取り扱われる一般形微分方程式（PHOENICS ver. 2009 ユーザーマニュアル 第3章「PHOENICSの基礎」p.23より転載）

事例１；界面振動型マイクロミキサー（４）計算ステップ数（非定常計算） ，繰り返し計算回数，MOFOR機能設定• 計算時間；5.0 [s]，計算ステップ数；1600（0.03125秒間隔），繰り返し計算回数；3500回• 実験画像から，1周期内8分割時の気泡の振動変位を確認し，*mof ファイルに反映．

Fig. 9 Moforを用いた気泡オブジェクトによる移動境界流れ解析の様子と*mof ファイルの一例

事例１；界面振動型マイクロミキサー（５）初期条件作成，読み込み① Moferを適用していない解析モデルで定常計算での事前計算を行う．
② 計算結果のうち，Phiファイルと境界条件・壁面条件（pbcファイル）を非定常計算での初期条件として，Restart機能を用いてImportする．
③ 本計算を行う．

Fig. 10 解析モデルの事前計算結果（ = 非定常計算時の初期条件）

事例１；界面振動型マイクロミキサー（６）計算結果(a)(b)

Fig. 10 解析モデルの計算結果（ 気泡オブジェクト付近），(a)対向形，(b)交互型

事例１；界面振動型マイクロミキサー（7）解析結果• CFD解析結果は，実験値と良い相関を示した．
• 「交互型」は「対向型」に比べ，高い混合効果を示した．
• 気泡オブジェクト（移動境界面）の位置により，混合効果に変化が見られた．

Putting point ofFirst bubble
Putting point ofSecond bubble

Relative Concentration Difference(%)

1000(s)0.125(s)0.25(s)0.375(s)0.5(s)
0.625(s)0.75(s)0.875(s)8060

40

Most extended
point of bubble

2000Distance(mm)(a)

1234Distance(mm)56(b)

Fig. 11 CFDでの解析結果，(a) 実験結果との比較，(b)１周期内での混合効果の変化（交互型）

7

事例１；界面振動型マイクロミキサー（７）解析結果(a)(b)
0.125(s)0.25(s)0.375(s)0.5(s)
0.675(s)0.75(s)0.875(s)1.0(s)Y Velocity(μ m/s)

1000.00
0.00-1000.00-2000.00-3000.00

Position of Bubbles

4000.00

0.125(s)0.25(s)0.375(s)0.50(s)0.625(s)0.75(s)0.875(s)1.0(s)

Positions ofBubbles

3000.00
2000.00

Y Velocity(μ m/s)

2000.00

1000.00
0.00-1000.00-2000.00-4000.00-3000.00-5000.00
0.500.60

0.70

0.80

0.90

1.00

1.10

Distance(mm)

(c)

1.20

1.30

1.40

1.50

-4000.00
2.750

2.875

3.000

3.125

3.250

3.375

Distance(mm)

3.500

3.625

3.750

(d)

Fig. 12 気泡オブジェクト付近での垂直方向速度の解析結果，(a) 対向型，(b)交互型映像， (c) １周期内での速度変化（対向型），(d) １周期内での速度変化（交互型）

事例２；マランゴニ対流を用いたマイクロミキサー（１）モデル作成と計算の流れ
幾何的な計算モデル作成

境界条件（流入・流出速度，壁面条件）設定
従属変数（拡散係数）設定

計算ステップ数（非定常計算） ，繰り返し計算回数設定
Informを用いて気泡オブジェクト表面に速度境界を設置
初期条件作成，読み込み計算・解析

事例２；マランゴニ対流を用いたマイクロミキサー（１）マランゴニ効果とは• 気体と液体の接触面である気液界面において，液相内の濃度勾配が界面に対して平行に存在するとき，この界面で起こる表面張力（分子間力）の局所的な不均一によって，界面上にせん断力（マランゴニ応力）が生じる現象．

Liquid phase
Concentration； Highσ； small；Marangoni StressConcentration； Lowσ； large

ΔCx；Liquid Flow• Marangoni convection Gas phasexFree Surface• このせん断力（マランゴニ応力）により生じる対流をマランゴニ対流という．
• 対流発生条件に限界はなく，温度（濃度）勾配があれば対流が発生する．

Fig. 1 Marangoni effect (a) Definition of Marangoni effect, (b) Value of surface tension of acetic a事例２；マランゴニ対流を用いたマイクロミキサー（１）（1）計算モデル作成
解析領域x 方向成分；0.6[mm], y 方向成分；0.1[mm]• 解析流体の物性
z 方向成分；0.45 [mm]IMAT67 純水 at 20℃• 解析モデル；2次元モデル0.2mm• セル数 ；120×1×90 = 10800 cells

Inlet 10.45mm

0.2mm

Inlet 2• 流入条件；両Inlet共に2.5×10-11[m3/s]• 非定常計算：6秒間• 計算刻み幅Δt：0.05秒• 繰り返し計算回数：10000回

Fig. 13 解析モデルの概形

事例２；マランゴニ対流を用いたマイクロミキサー（１）（2） Informを用いて気泡オブジェクト表面に速度境界を設置• ナビエ-ストークス方程式
速度境界
Velocity
Boundarydydxdx
xu
 div( uu)  div(  grad u)  grad pliq  St

• 速度境界面に働く力（マランゴニ力）

f   dx dy  dy dxdy
x
x



生成項に代入（= f ）• 濃度計算値から算出される表面張力値dzzy

ST = -0.1571C13+0.2397C12-0.1378C1+0.0690y• 速度境界面に働く力 f はST を用いて以下のように表される．
xdx     ST [1]  ST
 x 
f   dx dy  dy dxdy
x
x


 ST [1]  ST(  )
xdx
ST [1]  STf dxdy  ST [1]  ST dydx

事例２；マランゴニ対流を用いたマイクロミキサー（１）（３）初期条件作成，読み込み① Informによる速度境界を適用していない解析モデルで定常計算での事前計算を行う．
② 計算結果のうち，Phiファイルと境界条件・壁面条件（pbcファイル）を非定常計算での初期条件として，Restart機能を用いてImportする．
③ 本計算を行う．

濃度勾配大大きなせん断力発生（本計算時）

Fig. 14 解析モデルの事前計算結果（ = 非定常計算時の初期条件）

事例２；マランゴニ対流を用いたマイクロミキサー（１）（４）計算・解析結果• 合流点付近に大きな渦流発生• 渦流は約1秒で完全に消滅（実験結果との整合性あり）• その理由として，界面オブジェクト上面の濃度分布が早い段階で飽和→せん断力低下

界面上の濃度分布を途切れさせない（維持・持続させる）ことが重要

0.2mm

Fig. 15 CFDでの解析結果，(a) 0.1秒毎の濃度分布の変化，(b) 0.1秒毎の速度分布の変化

事例３；マランゴニ対流を用いたマイクロミキサー（２）（１）解析モデルのコンセプト• CFD解析結果→実験流路作成（CFD主導）0.2mmOff set 0.1mm

Wall: Flow of pure water

Wall

INLET1

Bubble0.2mm

Wall0.1mm

INLET20.1mm

INLET3: Flow of acetic acidsolution: Flow of mixed fluid: Marangoni force

Wall

界面上に常に濃度勾配を供給

Fig. 16 新解析モデルのコンセプト

事例３；マランゴニ対流を用いたマイクロミキサー（２）（２）計算モデル概要• 解析流体の物性，IMAT67 純水 at 20℃• 解析モデル；2次元モデル，セル数 ；120×1×90 = 10800 cells]• 非定常計算：6秒間，計算刻み幅Δt：0.05秒，繰り返し計算回数：10000回→ 本計算までの流れは今までと同様
OutletOutlet0.15mm0.07mm0.10mm

0.30mm

Inlet 12.0μl/minConcentration
value; 0.000.30mm

OutletInterfacial object(Boundary condition; Full-slip on the object)

Fig. 17 解析モデルの概形

Inlet 32.0μl/minConcentration
value; 0.00

Inlet 21.0μl/minConcentration
value; 0.30

事例３；マランゴニ対流を用いたマイクロミキサー（２）（３）計算結果(a)(b)

Fig. 15 CFDでの計算結果，(a) 濃度分布の変化，(b) 速度分布の変化

事例３；マランゴニ対流を用いたマイクロミキサー（２）（４）計算・解析結果• 合流点付近に大きな渦流発生

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

(m)

(n)

(o)

(p)

• 渦流は解析時間中で維持（ある程度の実験結果との整合性あり）• 界面オブジェクト上面の濃度分布が界面右端で維持→せん断力保持• 実験時の流れの状態とは異なる挙動→２次元（CFDモデル）と３次元（実験時）との差異

Fig. 19 CFDでの解析結果，0.05秒毎の濃度分布の変化
（１秒間）
(q)

(r)

(s)

(t)

結言
芝浦工大工学部機械機能工学科ナノ・マイクロ応用理工学研究室（小野研究室）では，「気液自由界面を用いたマイクロ混合器」に対する一連の研究に対して，「PHOENICS」を用いたCFD解析を実施してきた．
その結果，以下のような成果を上げることが出来た．
• 気液自由界面の変動を利用したマイクロ混合デバイスの性能評価に対して，「MOFOR」機能を用いた移動境界面解析により，簡便なモデルでの性能評価，混合メカニズムの解明を行うことが可能となった．
• マランゴニ対流を用いたマイクロ混合デバイスの性能評価に対しては，「Inform」機能を用いて，解析モデルのカスタマイズを簡単に行うことが出来，その性能評価，混合メカニズム，問題点を把握することだけの計算結果が得られた．

• それを基に，CFDモデルにて改善を行い，それを実際の実験装置製作に活かすことで，混合デバイスの性能改善につなげることが可能となった．

御静聴ありがとうございました．