鄭馨瑤 楊 剛 張會(huì)臣

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

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

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

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

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

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

1.1 模擬構(gòu)型

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

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

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

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

1.2 模擬方法

首先采用格子Boltzmann方法求解流體速度場(chǎng)，在引入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

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

(2)

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

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

(3)

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

(4)

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

(5)

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

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

(6)

(7)

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

(8)

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

(9)

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

2 結(jié)果及分析

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

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

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

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

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

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

圖6 不同入口速度下分叉點(diǎn)截面處的 壓力降(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ù)一定時(shí)，分叉角度為120°時(shí)，壓力降的最大值相差2.32%，分叉角度為105°時(shí)，壓力降的最大值相差5.92%。從圖7(a)、(c)中可以看出，當(dāng)分叉角度一定時(shí)，溶液質(zhì)量分?jǐn)?shù)為0.1%的情況下，壓力降的最大值相差2.32%，溶液質(zhì)量分?jǐn)?shù)為0.3%時(shí)，壓力降的最大值相差18.42%?？梢?jiàn)，隨著溶液質(zhì)量分?jǐn)?shù)的增大，使入口速度對(duì)出口1處的壓力降影響更大。所以相比于分叉角度，溶液質(zhì)量分?jǐn)?shù)使入口速度對(duì)出口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)分叉角度一定時(shí)，溶液質(zhì)量分?jǐn)?shù)為0.1%時(shí)，不同入口速度下的微通道出口2截面處的壓力降趨勢(shì)一致，而溶液質(zhì)量分?jǐn)?shù)為0.3%時(shí)，不同入口速度下壓力降的最大值相差2.69%，且壓力降隨著入口速度的增大而減小?？梢?jiàn)隨著溶液質(zhì)量分?jǐn)?shù)的增大使入口速度對(duì)出口1處的壓力降影響更大，所以相比于分叉角度，溶液質(zhì)量分?jǐn)?shù)使入口速度對(duì)出口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ù)對(duì)微通道內(nèi)壓力降的影響

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

圖9 不同溶液質(zhì)量分?jǐn)?shù)下分叉點(diǎn)截面 處的壓力降(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%時(shí)，壓力降最大，隨著溶液質(zhì)量分?jǐn)?shù)的降低，壓力降減小。由于水是非牛頓流體，沒(méi)有流變特性，導(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處壓力值隨時(shí)間的變化趨勢(shì)如圖11所示。根據(jù)圖11將模擬過(guò)程分為3個(gè)階段，在模擬初期(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)，壓力值迅速降低，但是在水的模擬過(guò)程中，模擬開(kāi)始不久就達(dá)到穩(wěn)定，沒(méi)有壓力值降低的階段。所以剪切稀化流體在流動(dòng)過(guò)程中產(chǎn)生了壓力值大幅度降低的現(xiàn)象，分析其原因是由于剪切稀化流體的流變特性導(dǎo)致的，因?yàn)槠漯ざ入S著剪切速率的增加而減小，導(dǎo)致壓力也隨之減小。第三階段即模擬穩(wěn)定后(3 000～5 000)，壓力值的降低幅度隨著溶液質(zhì)量分?jǐn)?shù)的增大而增大，與出口截面處壓力降結(jié)果一致。

圖11 不同溶液質(zhì)量分?jǐn)?shù)下出口1處壓力隨時(shí)間的 變化趨勢(shì)(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%時(shí)，壓力降最大；隨著溶液質(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é)論

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

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

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

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