微流控芯片通道結構的拓撲優化研究

董馨,劉小民

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

微流控芯片的微流動通道加工技術和結構布置一直為學術界與工業界所重視,它除了連接流體器件的出口和入口、起到連通作用之外,還具有流量分配、試劑混合等功能,同時也希望流體流經微流通道時能夠將某些目標值降到最小,如能量耗散、阻力、壓降等,以滿足微流控芯片能耗低的優點。

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

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

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

1 數學描述

1.1 牛頓流體流動的描述

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

(1)

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

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

(2)

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

1.2 非牛頓流體流動的描述

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

(3)

式中:μ(γ1)為動力黏度;μ∞和μ0分別代表著無限剪切黏性和零剪切黏性。參數取值[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流動拓撲優化問題

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

f=-α(γ)u

(5)

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

(6)

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

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

(7)

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

這里,分別選擇壓降和能量耗散為優化目標函數,相應的函數表達式如下

J1(Ω)=

(8)

(9)

圖1 雙進-單出微通道設計區域示意圖

2 數值算例

2.1 最小化壓力降

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

(10)

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

(11)

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

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

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

表1 壓力降函數的目標函數值

(a)短管 (b)方管

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

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

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

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

(a)短管 (b)方管

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

(a)短管 (b)方管

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

(a)短管 (b)方管

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

(a)短管 (b)方管

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

(a)短管 (b)方管

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

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

(a)短管 (b)方管

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

(a)短管 (b)方管

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

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

2.2 最小化能量耗散

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

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

(a)短管 (b)方管

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

(a)短管 (b)方管

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

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

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

(a)短管 (b)方管

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

3 結 論

本文針對兩種雷諾數條件下的雙進-單出型微流道結構進行拓撲優化研究,獲得的主要結論如下。

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

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

(3)在非牛頓流體的Stokes流動中,當目標函數為壓力降函數時,雙進-單出微通道的匯合點不隨通道長度的改變產生變化;當目標函數為能量耗散函數時,雙進-單出微通道的匯合點隨著微通道長度的增加而后延。這表明,微流通道壓力降主要體現在雙管流動時的沿程損失,匯合變為單管能夠有效減小損失,從而降低壓力降;而能量耗散與雙管結構的角度、長度、匯合點都有關,主要是降低微通道內流體流動時的黏性耗散與匯合時流體間的摩擦耗散。

參考文獻:

[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] 王志鑫. 用于混合及微液柱制備的微流控芯片設計與制備 [D]. 北京: 北京化工大學,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.