董馨,劉小民

(西安交通大學(xué)能源與動(dòng)力工程學(xué)院,710049,西安)

微流控芯片的微流動(dòng)通道加工技術(shù)和結(jié)構(gòu)布置一直為學(xué)術(shù)界與工業(yè)界所重視,它除了連接流體器件的出口和入口、起到連通作用之外,還具有流量分配、試劑混合等功能,同時(shí)也希望流體流經(jīng)微流通道時(shí)能夠?qū)⒛承┠繕?biāo)值降到最小,如能量耗散、阻力、壓降等,以滿足微流控芯片能耗低的優(yōu)點(diǎn)。

文獻(xiàn)[1]采用電學(xué)比擬法研究了Stokes流動(dòng)的流道流量分配設(shè)計(jì),但電學(xué)比擬法只能作為一個(gè)近似求解方法,不能準(zhǔn)確地設(shè)計(jì)高精度的流體網(wǎng)絡(luò)。文獻(xiàn)[2]設(shè)計(jì)了一種用于試劑混合的微通道,由T型微通道和Y型微通道組合而成,其中T型微通道用于制備微液柱,Y型微通道和T型微通道組合用于兩種試劑混合。現(xiàn)有的這些微流設(shè)備的設(shè)計(jì)大部分都只是基于設(shè)計(jì)者的經(jīng)驗(yàn)和直覺,或者依賴于研究者的嘗試改進(jìn)方法進(jìn)行的。這些設(shè)計(jì)方法帶有較大的隨機(jī)性,很難獲得具有最佳性能的設(shè)計(jì)構(gòu)型。

部分研究者基于數(shù)值分析技術(shù)的尺寸優(yōu)化、形狀優(yōu)化和拓?fù)鋬?yōu)化方法對(duì)微流器件進(jìn)行設(shè)計(jì)和改進(jìn)。文獻(xiàn)[3]首次將拓?fù)鋬?yōu)化設(shè)計(jì)引入流體力學(xué)領(lǐng)域。文獻(xiàn)[4-6]對(duì)中低雷諾數(shù)下的二維定常不可壓縮黏性層流流動(dòng)進(jìn)行拓?fù)鋬?yōu)化設(shè)計(jì),并將Stokes流動(dòng)擴(kuò)展到Navier-Stokes流動(dòng),以此對(duì)文獻(xiàn)[3]的工作進(jìn)行延伸。拓?fù)鋬?yōu)化作為更高層次的優(yōu)化方法,得到的設(shè)計(jì)結(jié)果更加合理,且不受原有微流通道拓?fù)浣Y(jié)構(gòu)的限制。文獻(xiàn)[7]首次采用密度類拓?fù)鋬?yōu)化方法,通過對(duì)不同雷諾數(shù)下牛頓流體和非牛頓流體進(jìn)行拓?fù)鋬?yōu)化,研究對(duì)象是Carreau-Yasuda非牛頓本構(gòu)模型。Jensen等采用變密度方法研究了記憶型非牛頓流體的拓?fù)鋬?yōu)化[8-9]。文獻(xiàn)[10]運(yùn)用變密度方法,研究了剪切漸薄效應(yīng)的非牛頓流體拓?fù)鋬?yōu)化,以達(dá)到壁面應(yīng)力最小化的目標(biāo)。文獻(xiàn)[11]研究了冪率型非牛頓流體拓?fù)鋬?yōu)化,并比較了剪切漸薄和剪切漸厚非牛頓流體拓?fù)鋬?yōu)化的設(shè)計(jì)結(jié)果。

基于以上分析,本文采用基于真實(shí)血液改進(jìn)的Cross非牛頓流體模型,對(duì)雙進(jìn)口-單出口的微流道進(jìn)行結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計(jì),分別以能量耗散函數(shù)和壓降函數(shù)為目標(biāo),比較分析了非牛頓效應(yīng)和流道長度對(duì)最優(yōu)拓?fù)浣Y(jié)構(gòu)的影響。

1 數(shù)學(xué)描述

1.1 牛頓流體流動(dòng)的描述

在微流控芯片試樣中,大部分流體屬于牛頓流體且滿足不可壓縮條件,特征是應(yīng)力和應(yīng)變滿足本構(gòu)方程

(1)

式中:σ為應(yīng)力張量;ε為應(yīng)變張量;u為流體的流動(dòng)速度;p為流體壓強(qiáng);μ為流體的動(dòng)力黏度。

設(shè)Ω是一個(gè)具有Lipschitz連續(xù)邊界Γ=?Ω的有界開區(qū)域,代表的是流體區(qū)域,則二維不可壓縮定常流動(dòng)的Navier-Stokes方程可表示為

(2)

式中:ρ為流體密度;f為流體微元所受體積力。

1.2 非牛頓流體流動(dòng)的描述

微流控芯片在生物醫(yī)學(xué)行業(yè)具有廣泛應(yīng)用,這些應(yīng)用包括遺傳和單細(xì)胞分析、蛋白質(zhì)研究、細(xì)胞遷移、藥物篩選、干細(xì)胞和神經(jīng)細(xì)胞培養(yǎng)等[12-14],其中常涉及血液這一典型的非牛頓流體。非牛頓流體的應(yīng)力和應(yīng)變率之間不遵循牛頓內(nèi)摩擦定律,考慮血液的剪切漸薄效應(yīng),本文將動(dòng)力黏度選擇為改進(jìn)的Cross模型[15]

(3)

式中:μ(γ1)為動(dòng)力黏度;μ∞和μ0分別代表著無限剪切黏性和零剪切黏性。參數(shù)取值[16-17]:μ∞=3.5 mPa·s,μ0=0.16 Pa·s,λ=0.82 s,a=1.23,b=0.64,血液密度ρ=1.058 g/cm3。其中,剪切率γ1表達(dá)式如下

(4)

1.3 Navier-Stokes流動(dòng)拓?fù)鋬?yōu)化問題

文獻(xiàn)[3]提出了定常Stokes流動(dòng)的密度法拓?fù)鋬?yōu)化模型,隨后文獻(xiàn)[5]提出定常Navier-Stokes流動(dòng)的密度法拓?fù)鋬?yōu)化模型。該模型是通過在Navier-Stokes方程中加入人工Darcy摩擦力進(jìn)行的,將式(2)的體力項(xiàng)設(shè)為

f=-α(γ)u

(5)

式中:α可以用來表征多孔介質(zhì)內(nèi)流動(dòng)的流體滲透率,它是設(shè)計(jì)變量(即優(yōu)化變量)γ的一個(gè)插值函數(shù)

(6)

式中:αmin和αmax分別為α(γ)的最大值和最小值,αmin一般取為0,αmax取值越大則黏滯力越大,固體滲透率越小(原則上滲透越小越接近真實(shí)的情況)。但是,由于拓?fù)鋬?yōu)化中密度法依賴于流體的滲透,在αmax取值太大時(shí)容易出現(xiàn)數(shù)值不穩(wěn)定的情況,一般取足夠大但有限的值以保證優(yōu)化問題的數(shù)值穩(wěn)定性且接近真實(shí)情況。設(shè)計(jì)變量γ在0與1之間取值,γ取0對(duì)應(yīng)于固相材質(zhì),γ取1對(duì)應(yīng)于液相材質(zhì)。q為正實(shí)數(shù),用來調(diào)節(jié)插值曲線的凹凸性α(γ),q的取值會(huì)影響管道收斂的灰度,一般取值為1。

針對(duì)牛頓流體和非牛頓流體的Navier-Stokes流動(dòng),建立如下的拓?fù)鋬?yōu)化問題

(7)

計(jì)算區(qū)域示意圖如圖1所示,設(shè)計(jì)域的長、寬分別為L和H,長度單位為mm。在本文中,長度分別選取L=0.5H(短管),L=H(方管),L=2H(長管),以這3種長度為代表進(jìn)行優(yōu)化設(shè)計(jì),分析微流通道長度對(duì)通道結(jié)構(gòu)最優(yōu)形狀的影響。

這里,分別選擇壓降和能量耗散為優(yōu)化目標(biāo)函數(shù),相應(yīng)的函數(shù)表達(dá)式如下

J1(Ω)=

(8)

(9)

圖1 雙進(jìn)-單出微通道設(shè)計(jì)區(qū)域示意圖

2 數(shù)值算例

2.1 最小化壓力降

選擇最小化進(jìn)出口壓力降為優(yōu)化目標(biāo),即以降低微流通道的沿程壓力損失、增大流體的流動(dòng)速度為目標(biāo),改善流體在微流體器件中的流動(dòng)特性。微通道系統(tǒng)的壓降接近于系統(tǒng)總壓降,可近似為系統(tǒng)總的機(jī)械能損失,即流動(dòng)阻力損失。微通道系統(tǒng)的壓降表達(dá)如下

(10)

式中:γ2為流體的比重,γ2=ρg,g為重力加速度。高度變化引起的相對(duì)壓降和速度能可以忽略不計(jì),式(10)簡化為

(11)

考慮到微流通道進(jìn)出口流量的不均勻性,因此目標(biāo)函數(shù)為進(jìn)出口積分值之差,即以目標(biāo)函數(shù)為壓降的表達(dá)式見式(8)。

目標(biāo)函數(shù)為壓力降函數(shù)的非牛頓流體在Re=0.1(Stokes流動(dòng))和Re=100(Navier-Stokes流動(dòng))時(shí),對(duì)應(yīng)的3種模型(短管、方管、長管)的拓?fù)鋬?yōu)化設(shè)計(jì)結(jié)果見圖2和圖3。如圖2所示,當(dāng)Re=0.1時(shí),在3個(gè)模型中兩條通道雙管均在設(shè)計(jì)域的1.5 mm處匯合為單管,隨后從出口流出,表明在小雷諾數(shù)下(近似為Stokes流動(dòng)),微管道長度對(duì)于非牛頓流體雙管的匯合點(diǎn)的位置沒有影響。

由表1給出的目標(biāo)函數(shù)值可知,由于壓降相同,即系統(tǒng)的流動(dòng)損失相同,近似為沿程阻力損失相同,表明隨著微通道長度增加,雙管流動(dòng)部分的管徑必須增加,以保證沿程阻力損失不變。

表1 壓力降函數(shù)的目標(biāo)函數(shù)值

(a)短管 (b)方管

(c)長管圖2 Re=0.1時(shí)最小化壓力降的非牛頓流體微通道最優(yōu)拓?fù)浣Y(jié)構(gòu)及速度矢量分布

如圖3所示,當(dāng)Re=100時(shí),設(shè)計(jì)域中進(jìn)口處流動(dòng)的角度隨著流道長度增加而減小,在長管中進(jìn)口角為0°,3個(gè)模型中兩條通道雙管分別在設(shè)計(jì)域的1.9、3、5.1 mm處匯合為單管。計(jì)算結(jié)果表明,在較大雷諾數(shù)下,雙管的匯合點(diǎn)會(huì)隨著管道的長度增加向后延遲,同時(shí)雙管流動(dòng)部分的管徑也隨之增加。原因是:Stokes流體流動(dòng)速度較小,黏性力占優(yōu),在雙管流動(dòng)中黏性摩擦較大,減少雙管流動(dòng)部分可有效減小壓力損耗;而Navier-Stokes流體流速較大,慣性力占優(yōu),黏性摩擦的作用較小,慣性力促使雙管流動(dòng)部分增加。

將圖2與圖3相同模型的最優(yōu)結(jié)構(gòu)進(jìn)行對(duì)比分析,結(jié)果表明雷諾數(shù)的增大導(dǎo)致匯合點(diǎn)后移。這是由于隨著雷諾數(shù)的增大,黏性力與慣性力相比逐漸可以忽略不計(jì),流動(dòng)速度增大會(huì)使流體在雙管中流動(dòng)時(shí)間更長,從而延后了匯合點(diǎn)。

圖4和圖5分別給出了圖2和圖3對(duì)應(yīng)的動(dòng)力黏度云圖,表明通道中心處流體黏性最大,流動(dòng)時(shí)由流體本身產(chǎn)生的摩擦阻力最大。Re=0.1比Re=100時(shí)的流體的黏性梯度變化大,這會(huì)導(dǎo)致非牛頓效應(yīng)在Stokes流動(dòng)中更加突出;在Navier-Stokes流動(dòng)的最優(yōu)拓?fù)浣Y(jié)構(gòu)中,由于非牛頓流體的黏性梯度變化不明顯,使得它與牛頓流體的最優(yōu)結(jié)構(gòu)相似,此時(shí)可以忽略非牛頓流體黏性效應(yīng)的作用。

(a)短管 (b)方管

(c)長管圖3 Re=100時(shí)最小化壓力降的非牛頓流體微通道最優(yōu)拓?fù)浣Y(jié)構(gòu)及速度矢量分布

(a)短管 (b)方管

(c)長管圖4 Re=0.1時(shí)最小化壓力降的非牛頓流體黏性分布

(a)短管 (b)方管

(c)長管圖5 Re=100時(shí)最小化壓力降的非牛頓流體黏性分布

(a)短管 (b)方管

(c)長管圖6 Re=0.1時(shí)最小化壓力降的牛頓流體微通道最優(yōu)拓?fù)浣Y(jié)構(gòu)及速度矢量分布

(a)短管 (b)方管

(c)長管圖7 Re=100時(shí)最小化壓力降的牛頓流體微通道最優(yōu)拓?fù)浣Y(jié)構(gòu)及速度矢量分布

目標(biāo)函數(shù)為壓力降的牛頓流體在Re=0.1(Stokes流動(dòng))和Re=100(Navier-Stokes流動(dòng))時(shí),對(duì)應(yīng)的3種模型(短管、方管、長管)的拓?fù)鋬?yōu)化設(shè)計(jì)結(jié)果見圖6和圖7。對(duì)比分析可以看到:當(dāng)Re=0.1和Re=100時(shí),隨著模型長度的增加,兩條流道匯合點(diǎn)的水平距離都在增加;同一種模型下,不同雷諾數(shù)下的最優(yōu)流道結(jié)構(gòu)也基本相同。由此可得出結(jié)論:當(dāng)流動(dòng)介質(zhì)為牛頓流體時(shí),在雙-單管通道的結(jié)構(gòu)中,層流小雷諾數(shù)下的最優(yōu)拓?fù)浣Y(jié)構(gòu)不隨雷諾數(shù)產(chǎn)生較明顯的改變。

(a)短管 (b)方管

(c)長管圖8 Re=0.1時(shí)非牛頓流體與牛頓流體微通道最優(yōu)結(jié)構(gòu)的交叉對(duì)比

(a)短管 (b)方管

(c)長管圖9 Re=100時(shí)非牛頓流體與牛頓流體微通道最優(yōu)結(jié)構(gòu)的交叉對(duì)比

圖8和圖9分別交叉對(duì)比了非牛頓流體與牛頓流體在Re=0.1和Re=100時(shí)同一工況下的最優(yōu)結(jié)構(gòu)。當(dāng)Re=0.1時(shí),非牛頓流體微通道的匯合點(diǎn)都比牛頓流體匯合點(diǎn)靠前,且隨著微通道長度的增加,2個(gè)匯合點(diǎn)的距離也隨之增加;當(dāng)Re=100時(shí),兩條線基本重合,匯合點(diǎn)都隨著微通道長度的增加而產(chǎn)生后移。這說明:在以壓力降為目標(biāo)函數(shù)的雙進(jìn)-單出管微流通道最優(yōu)拓?fù)浣Y(jié)構(gòu)中,非牛頓效應(yīng)主要體現(xiàn)在黏性力占優(yōu)的Stokes流動(dòng)中;對(duì)于慣性力不可忽略的Navier-Stokes流動(dòng)來說,當(dāng)介質(zhì)為非牛頓流體時(shí),可將其近似作為牛頓流體處理,以減少計(jì)算和設(shè)計(jì)的工作量。

2.2 最小化能量耗散

微流控芯片的特點(diǎn)之一為低能耗,這就要求流體在流經(jīng)通道時(shí)的能量耗散盡可能小,能量耗散的目標(biāo)函數(shù)表達(dá)見式(9)。

由上節(jié)的結(jié)論可知,非牛頓流體在Stokes流動(dòng)中的非牛頓效應(yīng)比Navier-Stokes流動(dòng)中明顯,而且在微流控芯片中,流動(dòng)速度非常低,相較于較大的黏性力可忽略掉慣性項(xiàng),近似為Stokes流動(dòng)。因此,在下面的算例中,只選取Re=0.1為例進(jìn)行優(yōu)化設(shè)計(jì)。

(a)短管 (b)方管

(c)長管圖10 最小化能量耗散的非牛頓流體微通道最優(yōu)拓?fù)浣Y(jié)構(gòu)及速度矢量分布

(a)短管 (b)方管

目標(biāo)函數(shù)為能量耗散函數(shù)的非牛頓流體在Re=0.1(Stokes流動(dòng))時(shí),對(duì)應(yīng)的3種模型拓?fù)鋬?yōu)化設(shè)計(jì)結(jié)果見圖10,圖11為黏性分布云圖。由此可知,隨著微流通道變長,匯合點(diǎn)不斷延后。當(dāng)優(yōu)化目標(biāo)為最小化壓力降函數(shù)時(shí),要求流體降低沿程損失,以增加流體流動(dòng)速度,因此兩條流道提前匯合可以減小沿程損失,達(dá)到壓力降最小的目標(biāo)。當(dāng)目標(biāo)函數(shù)為能量耗散函數(shù),考慮到兩條流道匯合時(shí)所產(chǎn)生的黏性耗散和摩擦耗散,模型長度增加使匯合點(diǎn)產(chǎn)生一定的延后,有助于使能量耗散更小。黏性分布圖表明,通道中心處的黏性較大,匯合后的黏性比雙管流動(dòng)的黏性有所減小。短管、方管和長管的目標(biāo)函數(shù)值分別為7.4×10-7、8.2×10-7和1.2×10-6W/m,表明通過最優(yōu)拓?fù)湓O(shè)計(jì)只能使能量耗散增幅減小,但隨著長度增加必然會(huì)引起能量耗散值的增大。

(c)長管圖11 最小化能量耗散的非牛頓流體微通道黏性分布云圖

當(dāng)目標(biāo)函數(shù)為能量耗散函數(shù)的牛頓流體在Re=0.1(Stokes流動(dòng))時(shí),對(duì)應(yīng)的3種模型拓?fù)鋬?yōu)化設(shè)計(jì)結(jié)果見圖12。短管和方管最優(yōu)結(jié)構(gòu)的雙管在設(shè)計(jì)域中沒有匯合,而長管的最優(yōu)結(jié)構(gòu)是在約1/3處匯合為單管,而后在約2/3處又分流為上下對(duì)稱的兩條流道,最后在設(shè)計(jì)域的出口處經(jīng)兩條流道流出。長管匯合為單通道后又分流為雙通道的原因是為了降低流動(dòng)速度,從而減小流動(dòng)中的黏性摩擦力,且雙管中總的耗散比繼續(xù)保持單管流動(dòng)狀態(tài)的耗散小,因此最優(yōu)拓?fù)浣Y(jié)構(gòu)如圖12c所示。目標(biāo)函數(shù)值分別為5.1×10-9、6×10-9和1.4×10-8W/m。

(a)短管 (b)方管

(c)長管圖12 最小化能量耗散函數(shù)的牛頓流體微通道最優(yōu)拓?fù)浣Y(jié)構(gòu)及其速度矢量分布

3 結(jié) 論

本文針對(duì)兩種雷諾數(shù)條件下的雙進(jìn)-單出型微流道結(jié)構(gòu)進(jìn)行拓?fù)鋬?yōu)化研究,獲得的主要結(jié)論如下。

(1)非牛頓效應(yīng)對(duì)Stokes流體的敏感性較高,在雙進(jìn)-單出通道的最小化壓力降最優(yōu)拓?fù)浣Y(jié)構(gòu)中,雙管的匯合點(diǎn)不隨微通道長度改變發(fā)生變化;介質(zhì)為牛頓流體時(shí),隨著流道長度增加,匯合點(diǎn)也隨之后延。上述差別在Navier-Stokes流動(dòng)中并沒有體現(xiàn),因此在基于Stokes流動(dòng)的微流控芯片通道結(jié)構(gòu)設(shè)計(jì)中,要充分考慮非牛頓效應(yīng)。

(2)在最小化壓力降為目標(biāo)的牛頓流體微流道拓?fù)鋬?yōu)化設(shè)計(jì)中,層流小雷諾數(shù)(包括Re為0.1和100)下,相同雙進(jìn)-單出微通道模型的最優(yōu)拓?fù)浣Y(jié)構(gòu)基本相同,不隨雷諾數(shù)產(chǎn)生改變。

(3)在非牛頓流體的Stokes流動(dòng)中,當(dāng)目標(biāo)函數(shù)為壓力降函數(shù)時(shí),雙進(jìn)-單出微通道的匯合點(diǎn)不隨通道長度的改變產(chǎn)生變化;當(dāng)目標(biāo)函數(shù)為能量耗散函數(shù)時(shí),雙進(jìn)-單出微通道的匯合點(diǎn)隨著微通道長度的增加而后延。這表明,微流通道壓力降主要體現(xiàn)在雙管流動(dòng)時(shí)的沿程損失,匯合變?yōu)閱喂苣軌蛴行p小損失,從而降低壓力降;而能量耗散與雙管結(jié)構(gòu)的角度、長度、匯合點(diǎn)都有關(guān),主要是降低微通道內(nèi)流體流動(dòng)時(shí)的黏性耗散與匯合時(shí)流體間的摩擦耗散。

參考文獻(xiàn):

[1] OH K W,LEE K,AHN B,et al. Design of pressure-driven microfluidic networks using electric circuit analogy [J]. Lab on A Chip,2012,12(3): 515-545.

[2] 王志鑫. 用于混合及微液柱制備的微流控芯片設(shè)計(jì)與制備 [D]. 北京: 北京化工大學(xué),2014: 1-9.

[3] BORRVALL T,PETERSSON J. Topology optimization of fluids in stokes flow [J]. International Journal for Numerical Methods in Fluids,2010,41(1): 77-107.

[4] EVGRAFOV A. Topology optimization of slightly compressible fluids & dagger [J]. Journal of Applied Mathematics and Mechanics,2010,86(1): 46-62.

[5] GERSBORG-HANSEN A,SIGMUND O,HABER R B. Topology optimization of channel flow problems [J]. Structural & Multidisciplinary Optimization,2005,30(3): 181-192.

[6] LIU Z,GAO Q,ZHANG P,et al. Topology optimization of fluid channels with flow rate equality constraints [J]. Structural & Multidisciplinary Optimization,2011,44(1): 31-37.

[7] PINGEN G,MAUTE K. Optimal design for non-Newtonian flows using a topology optimization approach [J]. Computers & Mathematics with Applications,2010,59(7): 2340-2350.

[8] EJLEBJERG J K,SZABO P,OKKELS F. Topology optimization of viscoelastic rectifiers [J]. Applied Physics Letters,2012,100(23): 234102.

[9] JENSEN K E,SZABO P,OKKELS F. Optimization of bistable viscoelastic systems [J]. Structural & Multidisciplinary Optimization,2014,49(5): 733-742.

[10] HYUN J,WANG S,YANG S. Topology optimization of the shear thinning non-Newtonian fluidic systems for minimizing wall shear stress [J]. Computers & Mathematics with Applications,2014,67(5): 1154-1170.

[11] ZHANG B,LIU X M,SUN J J. Topology optimization design of non-Newtonian roller-type viscous micropumps [J]. Structural & Multidisciplinary Optimization,2016,53(3): 409-424.

[12] BLAZEJ R G,KUMARESAN P,MATHIES R A. Microfabricated bioprocessor for integrated nanoliter-scale Sanger DNA sequencing [J]. Proceedings of the National Academy of Sciences of the United States of America,2006,103(19): 7240-7245.

[13] GAO J,YIN X F,FANG Z L. Integration of single cell injection,cell lysis,separation and detection of intracellular constituents on a microfluidic chip [J]. Lab on A Chip,2004,4(1): 47-52.

[14] LAGALLY E T,SCHERER J R,BLAZEJ R G,et al. Integrated portable genetic analysis microsystem for pathogen/infectious disease detection [J]. Analytical Chemistry,2004,76(11): 3162-3170.

[15] ZHANG B,LIU X. Topology optimization study of arterial bypass configurations using the level set method [J]. Structural & Multidisciplinary Optimization,2014,51(3): 1-26.

[16] ABRAHAM F,HEINKENSCHLOSS M B. Shape optimization in unsteady blood flow: a numerical study of non-Newtonian effects [J]. Computer Methods in Biomechanics & Biomedical Engineering,2005,8(2): 127-137.

[17] CHEN H,SU J,LI K,et al. A characteristic projection method for incompressible thermal flow [J]. Numerical Heat Transfer: Part B Fundamentals,2014,65(6): 554-590.