董馨,劉小民

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

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

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

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

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

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

1.1 牛頓流體流動的描述

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

(1)

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

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

(2)

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

1.2 非牛頓流體流動的描述

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

(3)

式中:μ(γ1)為動力黏度;μ∞和μ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表達式如下

(4)

1.3 Navier-Stokes流動拓撲優(yōu)化問題

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

f=-α(γ)u

(5)

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

(6)

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

針對牛頓流體和非牛頓流體的Navier-Stokes流動,建立如下的拓撲優(yōu)化問題

(7)

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

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

J1(Ω)=

(8)

(9)

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

2 數(shù)值算例

2.1 最小化壓力降

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

(10)

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

(11)

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

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

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

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

(a)短管 (b)方管

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

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

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

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

(a)短管 (b)方管

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

(a)短管 (b)方管

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

(a)短管 (b)方管

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

(a)短管 (b)方管

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

(a)短管 (b)方管

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

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

(a)短管 (b)方管

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

(a)短管 (b)方管

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

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

2.2 最小化能量耗散

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

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

(a)短管 (b)方管

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

(a)短管 (b)方管

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

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

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

(a)短管 (b)方管

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

3 結(jié) 論

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

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

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

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

參考文獻:

[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è)計與制備 [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.