基于格子Boltzmann方法的分叉微通道內(nèi)非牛頓流體流動特性研究*

鄭馨瑤 楊 剛 張會臣

(大連海事大學(xué)船舶與海洋工程學(xué)院 遼寧大連 116026)

微通道是指水力直徑在微米量級的流動通道，具有高的表面體積比、高的傳熱效率及較少的流體使用量等特點，廣泛應(yīng)用于微流體系統(tǒng)[1-3]。由于尺度的限制，微通道內(nèi)流體的流動特性很難通過實驗完全揭示清楚，并且有些實驗現(xiàn)象的分析也亟待深入。對微通道內(nèi)流體流動進(jìn)行數(shù)值模擬，是對實驗中難以觀測記錄現(xiàn)象和對實驗現(xiàn)象加以分析的有效手段[4]。這些數(shù)值模擬方法主要包括計算流體動力學(xué)、格子玻爾茲曼方法以及分子動力學(xué)。相比于宏觀尺度下的計算流體動力學(xué)方法與微觀尺度下的分子動力學(xué)方法，格子玻爾茲曼方法是一種介于流體連續(xù)性假設(shè)與分子動力學(xué)之間的方法，節(jié)省計算資源的同時，可揭示微通道內(nèi)流體流動的機制[5]。

早期一些學(xué)者借助格子Boltzmann方法模擬微通道內(nèi)氣體的流動。TANG等[6]采用格子玻爾茲曼方法模擬了微通道中的等溫氣體流動，將反彈邊界和鏡面反射格式組合模擬滑移邊界條件[7]，推導(dǎo)出了松弛時間與努森數(shù)的關(guān)系。LIU和GUO[8]采用格子玻爾茲曼方法研究了長矩形微通道內(nèi)壓力驅(qū)動的氣體流動，并與其他方法所得數(shù)據(jù)進(jìn)行了對比。此外，由于格子玻爾茲曼方法易于模擬復(fù)雜邊界，還被廣泛應(yīng)用于液滴動力學(xué)特性的模擬。梁宏等人[9]采用格子 Boltzmann 方法模擬了液滴在分叉微通道中的遷移過程,發(fā)現(xiàn)當(dāng)毛細(xì)數(shù)足夠大且壁面具有親水性時，液滴會發(fā)生二次分裂；而毛細(xì)數(shù)足夠小時，液滴在分叉處滯留并且不發(fā)生分裂，通過調(diào)節(jié)出口流量比可以使液滴實現(xiàn)非對稱分裂。LIM、KEIVAN等[10-11]采用格子玻爾茲曼方法模擬了微通道中液滴的形成，研究表明液滴的大小以及液滴之間的距離與微通道的幾何形狀關(guān)系密切。

應(yīng)用格子Boltzmann方法，通過局部計算剪切應(yīng)力可模擬非牛頓流體的流動。DONG等[12]將流變方程引入格子Boltzmann方法中，建立了非牛頓流體兩相流的格子Boltzmann模型，模擬了存在矩形粗糙單元的微通道內(nèi)冪律流體非牛頓流體的位移，發(fā)現(xiàn)冪律指數(shù)的降低導(dǎo)致位移速度的增加，相對粗糙度的減小和粗糙單元間距的增大可加大位移速度，從而減少位移時間。謝馳宇等[13]將非牛頓流體模型引入格子Boltzmann方法中，研究多孔復(fù)雜結(jié)構(gòu)內(nèi)多相非牛頓流體的驅(qū)替過程，發(fā)現(xiàn)了界面張力在驅(qū)替過程中起阻礙作用，并分析了壓差、潤濕性對驅(qū)替過程的影響。SHI和TANG[14]采用格子Boltzmann方法模擬黏性指進(jìn)現(xiàn)象，發(fā)現(xiàn)剪切增稠流體中的黏性指進(jìn)增強，剪切稀化流體中的黏性指進(jìn)被抑制，重力加速度導(dǎo)致“手指”形狀不對稱，高毛細(xì)數(shù)與高黏度比導(dǎo)致“手指”狹窄。該研究從介觀的角度解釋了黏性指進(jìn)現(xiàn)象。

盡管通過實驗研究與數(shù)值模擬，對微通道內(nèi)流體流動的過程有了一定程度的了解，但對于流動過程中的細(xì)節(jié)仍不夠清晰。本文作者采用格子Boltzmann方法，模擬不同角度的分叉微通道內(nèi)剪切稀化流體的流動，通過流動過程中密度隨時間的變化趨勢以及穩(wěn)態(tài)流動下的密度，得到微通道內(nèi)壓力的分布以及流動區(qū)間的壓力降，著重分析溶液質(zhì)量分?jǐn)?shù)、入口速度與分叉角度對非牛頓流體流動特性的影響，闡明流體特性和微通道幾何構(gòu)型對非牛頓流體流動行為的影響機制。

1 模型構(gòu)型與模擬方法

1.1 模擬構(gòu)型

文中模擬的是高度和寬度均為300 μm的矩形截面多角度分叉微通道，分叉角度分別為60°、75°、90°、105°和120°。圖1所示為分叉角度為120°的微通道結(jié)構(gòu)示意圖。微通道具有2個出口，非牛頓流體從左側(cè)入口以固定速度流入，流經(jīng)分叉點處后由2個出口流出。

圖1 微通道結(jié)構(gòu)示意Fig 1 Structure of microchannel

選取3種不同質(zhì)量分?jǐn)?shù)的聚丙烯酰胺(PAM)水溶液進(jìn)行模擬，探討流變特性對其流動特性的影響。聚丙烯酰胺水溶液為剪切稀化流體，其流變特性采用冪律模型描述，稠度系數(shù)K和冪律指數(shù)n2個參數(shù)根據(jù)流變儀測得的實驗數(shù)據(jù)擬合得到，K和n及其他特性參數(shù)如表1所示。

表1 流體的特性參數(shù)

1.2 模擬方法

首先采用格子Boltzmann方法求解流體速度場，在引入BGK碰撞算子[15]后，具有二次精度的格子Boltzmann方程離散為

(1)

平衡態(tài)分布函數(shù)的具體形式一般與離散速度模型的構(gòu)造有關(guān)，因此離散速度模型的構(gòu)造至關(guān)重要。文中采用DnQm模型[16]，D2Q9是最具代表性的格子BGK模型，二維D2Q9模型如圖2所示。

圖2 D2Q9模型Fig 2 D2Q9 model

對于二維D2Q9，其局部平衡分布函數(shù)的表達(dá)式為

(2)

式中：wk表示k方向的權(quán)重因子；cs為格子聲速；ck為沿著格子遷移方向的速度矢量。

流體的宏觀密度經(jīng)平衡分布函數(shù)求和得到

(3)

壓力由宏觀密度得到[17]

(4)

廣義牛頓流體中的冪律模型形式簡單，廣泛使用，且與剪切稀化流體的真實特性基本相符，所以選用冪律模型描述非牛頓流體的動力黏度

(5)

式中：K為稠度系數(shù)；n為冪律指數(shù)；γ為剪切速率。

采用根據(jù)剪切黏度調(diào)整松弛時間的直接方法，進(jìn)行格子玻爾茲曼方法中非牛頓流體特性的模擬[18-20]。首先計算格子的黏度，再根據(jù)黏度調(diào)整松弛時間。其中，通過式(6)迭代計算形變速率張量，根據(jù)式(7)對剪切速率進(jìn)行計算。

(6)

(7)

最后通過格子玻爾茲曼方法中松弛時間與黏度的關(guān)系調(diào)整松弛時間，實現(xiàn)剪切稀化特性的模擬。

(8)

在格子玻爾茲曼方法的模擬中，一般采用格子單位，與實際有關(guān)的問題模擬時需要進(jìn)行物理單位與格子單位的轉(zhuǎn)化。對于不可壓縮流體流動問題，需要匹配宏觀雷諾數(shù)和格子雷諾數(shù)。宏觀雷諾數(shù)經(jīng)驗公式為

(9)

通常，格子速度Ulattice的數(shù)量級保持在0.1或者0.2，取值與宏觀速度無關(guān)。文中Ulattice選擇0.1，微通道截面的實際尺寸為300 μm，其格子單位取值為14，匹配宏觀雷諾數(shù)和格子雷諾數(shù)后得到密度ρ的初始值。其中，根據(jù)微通道截面形狀，a的取值為0.212 1，b的取值為0.676 6。

2 結(jié)果及分析

2.1 分叉角度對微通道內(nèi)流體壓力降的影響

圖3所示是模擬趨于穩(wěn)定后入口截面處與分叉點處的壓力差值。數(shù)據(jù)整體分布在負(fù)值區(qū)域，表明從入口處流經(jīng)到分叉點處時壓力增大。從圖中可以看出，90°分叉角度壓力差值最?。划?dāng)分叉角度大于90°時，壓力差值隨著分叉角度的增大而減??；當(dāng)分叉角度小于90°時，壓力差值隨著分叉角度的增大而增大；壓力差值最大值出現(xiàn)在中線偏向分叉點壁面一側(cè)，越靠近壁面壓力差值越小。從圖中還可以發(fā)現(xiàn)，當(dāng)分叉角度互補時，其壓力差值結(jié)果相近，相差分別為10.72%和12.51%，且分叉角度小于90°的微通道壓力差值大于分叉角度大于90°的微通道。

圖3 不同分叉角度下分叉點截面處的 壓力降(w=0.1%,v=2 m/s)Fig 3 Pressure drop at the cross section of the bifurcation point in different bifurcation angles(w=0.1%,v=2 m/s)

出口1截面處的壓力降如圖4所示，整體分布趨勢呈拋物線形，壓力降最大值出現(xiàn)在微通道中間截面處，越靠近壁面，其壓力降值越小;90°分叉角度時壓力降最小，壓力降最大值出現(xiàn)在分叉角度為75°和105°處，最大值相差4.96%。從圖4(a)、(b)中可以看出，入口速度為2 m/s時，壓力降的最大值相差4.96%，入口速度為4 m/s時，壓力降的最大值相差8.94%。從圖4(a)、(c)中可以看出，溶液質(zhì)量分?jǐn)?shù)為0.1%時，壓力降的最大值相差4.96%，溶液質(zhì)量分?jǐn)?shù)為0.3%時，壓力降的最大值相差22.57%，說明隨著溶液質(zhì)量分?jǐn)?shù)的增大，分叉角度對出口1處的壓力降影響更大。因此，相比于入口速度，溶液質(zhì)量分?jǐn)?shù)使得分叉角度對出口1壓力降的影響更為顯著。

出口2截面處的壓力降如圖5所示，不論分叉角度如何，其壓力降趨勢一致。

2.2 入口速度對微通道內(nèi)壓力降的影響

圖6所示是模擬趨于穩(wěn)定后入口截面處與分叉點處的壓力差值，可以看出，入口速度為10 m/s時壓力差值最小，為2 m/s時壓力差值最大，即入口速度越大，壓力差值越小；而截面處壓力差值最大值出現(xiàn)在中線偏向分叉點壁面一側(cè)，越靠近壁面壓力差值越小。

圖6 不同入口速度下分叉點截面處的 壓力降(w=0.1%,θ=120°)Fig 6 Pressure drop at the cross section of the bifurcation point at different inlet velocities(w=0.1%,θ=120°)

出口1截面處的壓力降如圖7所示。從圖7(a)、(b)中可以看出，當(dāng)溶液質(zhì)量分?jǐn)?shù)一定時，分叉角度為120°時，壓力降的最大值相差2.32%，分叉角度為105°時，壓力降的最大值相差5.92%。從圖7(a)、(c)中可以看出，當(dāng)分叉角度一定時，溶液質(zhì)量分?jǐn)?shù)為0.1%的情況下，壓力降的最大值相差2.32%，溶液質(zhì)量分?jǐn)?shù)為0.3%時，壓力降的最大值相差18.42%。可見，隨著溶液質(zhì)量分?jǐn)?shù)的增大，使入口速度對出口1處的壓力降影響更大。所以相比于分叉角度，溶液質(zhì)量分?jǐn)?shù)使入口速度對出口1壓力降的影響更為顯著。

圖7 不同入口速度下出口1截面處的壓力降Fig 7 Pressure drop at the cross section of the outlet 1 at different inlet velocities (a) w=0.1%,θ=120°； (b)w=0.1%,θ=105°;(c)w=0.3%,θ=120°

出口2截面處的壓力降如圖8所示。從圖8(a)、(b)中可以看出，不同入口速度條件下出口2截面處壓力降基本相同。從圖8(a)、(c)中可以看出，當(dāng)分叉角度一定時，溶液質(zhì)量分?jǐn)?shù)為0.1%時，不同入口速度下的微通道出口2截面處的壓力降趨勢一致，而溶液質(zhì)量分?jǐn)?shù)為0.3%時，不同入口速度下壓力降的最大值相差2.69%，且壓力降隨著入口速度的增大而減小。可見隨著溶液質(zhì)量分?jǐn)?shù)的增大使入口速度對出口1處的壓力降影響更大，所以相比于分叉角度，溶液質(zhì)量分?jǐn)?shù)使入口速度對出口2壓力降的影響更為顯著。

圖8 不同入口速度下出口2截面處的壓力降Fig 8 Pressure drop at the cross section of the outlet 2 at different inlet velocities (a) w=0.1%,θ=120°； (b)w=0.1%,θ=105°;(c)w=0.3%,θ=120°

2.3 溶液質(zhì)量分?jǐn)?shù)對微通道內(nèi)壓力降的影響

圖9所示是模擬趨于穩(wěn)定后入口截面處與分叉點處壓力差值，可以看出，水的壓力差值比任一質(zhì)量分?jǐn)?shù)剪切稀化流體的壓力差值都小；溶液質(zhì)量分?jǐn)?shù)為0.5%時壓力差值最大，隨著溶液質(zhì)量分?jǐn)?shù)增大，壓力差值增大。溶液質(zhì)量分?jǐn)?shù)為0.3%和0.5%時，微通道分叉點截面處的壓力差值基本相同，分析其原因是在這2種質(zhì)量分?jǐn)?shù)下，描述其非牛頓特性的稠度系數(shù)K和冪律指數(shù)n的值相近，其流變特性的相似導(dǎo)致了流動過程中產(chǎn)生的壓力差值相近。

圖9 不同溶液質(zhì)量分?jǐn)?shù)下分叉點截面 處的壓力降(v=2 m/s,θ=120°)Fig 9 Pressure drop at the cross section of the bifurcation point at different solution concentrations(v=2 m/s,θ=120°)

出口1截面處的壓力降如圖10所示，壓力降的最大值相差3.01%;溶液質(zhì)量分?jǐn)?shù)為0.5%時，壓力降最大，隨著溶液質(zhì)量分?jǐn)?shù)的降低，壓力降減小。由于水是非牛頓流體，沒有流變特性，導(dǎo)致水的壓力降遠(yuǎn)遠(yuǎn)小于剪切稀化流體的壓力降。

圖10 不同溶液質(zhì)量分?jǐn)?shù)下出口1截面處的 壓力降(v=2 m/s,θ=120°)Fig 10 Pressure drop at the cross section of the outlet 1 at different solution concentrations(v=2 m/s,θ=120°)

出口1處壓力值隨時間的變化趨勢如圖11所示。根據(jù)圖11將模擬過程分為3個階段，在模擬初期(0～1 000)，壓力值隨著溶液質(zhì)量分?jǐn)?shù)的增大而增大，且溶液質(zhì)量分?jǐn)?shù)為0.3%與0.5%的微通道內(nèi)壓力峰值差值為27.5%，遠(yuǎn)遠(yuǎn)小于與溶液質(zhì)量分?jǐn)?shù)為0.1%的壓力峰值差值55.21%。在這一階段，水的壓力值已經(jīng)達(dá)到穩(wěn)定。隨著模擬的進(jìn)行(1 000～3 000)，溶液質(zhì)量分?jǐn)?shù)大的微通道內(nèi)，壓力值迅速降低，但是在水的模擬過程中，模擬開始不久就達(dá)到穩(wěn)定，沒有壓力值降低的階段。所以剪切稀化流體在流動過程中產(chǎn)生了壓力值大幅度降低的現(xiàn)象，分析其原因是由于剪切稀化流體的流變特性導(dǎo)致的，因為其黏度隨著剪切速率的增加而減小，導(dǎo)致壓力也隨之減小。第三階段即模擬穩(wěn)定后(3 000～5 000)，壓力值的降低幅度隨著溶液質(zhì)量分?jǐn)?shù)的增大而增大，與出口截面處壓力降結(jié)果一致。

圖11 不同溶液質(zhì)量分?jǐn)?shù)下出口1處壓力隨時間的 變化趨勢(v=2 m/s,θ=120°)Fig 11 The trend of pressure over time at the outlet 1 at different solution concentrations(v=2 m/s,θ=120°)

出口2截面處的壓力降如圖12所示。不同溶液質(zhì)量分?jǐn)?shù)條件下，壓力降的最大值相差3.55%；溶液質(zhì)量分?jǐn)?shù)為0.1%時，壓力降最大；隨著溶液質(zhì)量分?jǐn)?shù)的增大，壓力降下降。

圖12 不同溶液質(zhì)量分?jǐn)?shù)下出口2截面處 的壓力降(v=2 m/s,θ=120°)Fig 12 Pressure drop at the cross section of the outlet 2 at different solution concentrations(v=2 m/s,θ=120°)

3 結(jié)論

通過將非牛頓流體冪律模型引入牛頓流體格子Boltzmann模型的方法，數(shù)值模擬了矩形截面微通道內(nèi)剪切稀化流體的流動特性，探討了流體流變特性和微通道幾何構(gòu)型對非牛頓流體流動行為的影響，得到如下結(jié)論：

(1)分叉角度為90°的微通道壓力降最小，當(dāng)分叉角度大于90°時，壓力降隨著分叉角度的增大而減小，當(dāng)分叉角度小于90°時，壓力降隨著分叉角度的增大而增大，且當(dāng)微通道分叉角度在數(shù)學(xué)概念上互補時，其壓力降相近。

(2)微通道內(nèi)各截面處壓力降呈拋物線形，最大值出現(xiàn)在微通道中間截面處，越靠近壁面壓力降越小;在溶液質(zhì)量分?jǐn)?shù)與分叉角度一定時，隨著入口速度的增大，壓力降增大；相比于分叉角度，溶液質(zhì)量分?jǐn)?shù)使得入口速度對出口處壓力降的影響更為顯著。

(3)在入口速度與分叉角度一定時，隨著溶液質(zhì)量分?jǐn)?shù)的增大，壓力降增大，且當(dāng)流體的流變特性相近時，其流動特性相近;而相比于入口速度，溶液質(zhì)量分?jǐn)?shù)使得分叉角度對出口處壓力降的影響更為顯著。