郝鵬 張麗麗 丁明明2)?

1) (伊犁師范大學(xué)物理科學(xué)與技術(shù)學(xué)院,新疆凝聚態(tài)相變與微結(jié)構(gòu)實驗室,伊寧 835000)

2) (廣東工業(yè)大學(xué)輕工化工學(xué)院,廣州 510006)

采用基于流固耦合的有限元方法,對二維模型中高分子囊泡在微管流中慣性遷移現(xiàn)象進行了系統(tǒng)研究,分析了囊泡因受到流體作用力而形變并發(fā)生慣性遷移現(xiàn)象的機理.研究表明: 隨著雷諾數(shù)的增大,囊泡慣性遷移的平衡位置離其初始位置越來越遠;隨著阻塞比的增加,囊泡慣性遷移后的平衡位置越來越接近壁面.對于囊泡膜的模量和黏度以及膜厚,結(jié)果表明模量和黏度決定了囊泡的變形程度,模量對囊泡平衡位置影響較小,但增大黏度和膜厚會促進囊泡的平衡位置偏向管道中心.本研究有助于進一步明晰囊泡在慣性遷移過程中的形變和平衡位置,為囊泡在藥物輸運、化學(xué)反應(yīng)和生理過程的應(yīng)用提供可靠的計算依據(jù).

1 引言

囊泡是由兩親性分子有序組合且具有雙層膜結(jié)構(gòu)的一種組合體,其大小約為30 nm—100 μm,呈圓形或橢圓形.雙層膜內(nèi)外通常為簡單的黏性液體,在新陳代謝、物質(zhì)輸運、免疫識別、浮力調(diào)節(jié)、酶存儲與釋放等過程中發(fā)揮重要作用[1-3].由于囊泡膜與細(xì)胞膜的結(jié)構(gòu)非常相似,所以其作為簡單生物膜模型而得到廣泛的研究[4].

慣性遷移指隨機分布在微通道中的粒子與流體相互作用之后,達到指定的平衡位置,并平行于壁面移動.慣性微流控技術(shù)[5-6]是一種能夠?qū)崿F(xiàn)顆?；蚣?xì)胞高通量精確操控的新方法,廣泛應(yīng)用在生物、醫(yī)學(xué)、材料和環(huán)境等領(lǐng)域.慣性遷移原理是慣性微流控技術(shù)的核心,因此許多學(xué)者對微顆粒在流體中的動力學(xué)行為進行了廣泛的研究,例如剛性小球[7-9]、囊泡[10]、膠囊[11-12]或紅細(xì)胞[13-15]等.

最初,Segré和Silberberg[16]通過實驗發(fā)現(xiàn): 在圓管層流中,球形顆粒慣性遷移后的平衡位置約為管道半徑0.6 倍處的圓環(huán)區(qū)域,稱為“管縮效應(yīng)”.Segré和Silberberg 及之后的研究者認(rèn)為引發(fā)顆粒聚集的升力是由流場的慣性作用而產(chǎn)生的,因此被稱為顆粒慣性聚集現(xiàn)象.Carlo 等[17-18]采用模擬方法描述了較大中性懸浮顆粒的動力學(xué)行為,揭示了升力與通道和顆粒幾何形狀的復(fù)雜關(guān)系.本質(zhì)上,顆粒在微管流中慣性遷移是剪切梯度升力(指向通道壁面[19])和壁面誘導(dǎo)升力(指向通道中心[20])兩者相互平衡的結(jié)果.Carlo 等[21]研究了管道壁面、雷諾數(shù)及彎曲管道二次流等因素對細(xì)胞或顆粒慣性遷移的影響.Asmolov[7]研究了大雷諾數(shù)時剛性小球以及中性和非中性粒子的慣性遷移問題,采用匹配漸近展開的方法推測顆粒的平衡位置,研究表明曲率對無擾動速度分布的影響是顯著的,在大雷諾數(shù)下,升力的值與線性剪切流對應(yīng)的值不同.Matas 等[8]研究了雷諾數(shù)為67—1700 的泊肅葉流中懸浮粒子的遷移規(guī)律,并驗證了Segré-Silberberg[16]的理論,證明了粒子的管狀捏縮效應(yīng),并觀察到大雷諾數(shù)顆粒的遷移位于Segré-Silberberg 平衡位置[9].Morita 等[22]同樣驗證了Segré-Silberberg 提出的理論,并證明了隨著雷諾數(shù)增加,管道的入口長度對顆粒平衡位置的影響隨之增強.Yao 等[23]發(fā)現(xiàn)顆粒在高雷諾數(shù)圓管流體中會出現(xiàn)內(nèi)外兩個平衡位置.Nakayama 等[24]對圓管中懸浮顆粒進行數(shù)值模擬和實驗測量發(fā)現(xiàn)隨著雷諾數(shù)增加,不同顆粒慣性聚集平衡分布,存在三個區(qū)域與之相對應(yīng).

除了以上考慮的硬顆粒,國內(nèi)外學(xué)者對“軟顆?！?如囊泡等在微管流中的動力學(xué)行為也做了許多有意義的研究[25-28].Shin 等[29]采用二維數(shù)值模擬方法發(fā)現(xiàn): 當(dāng)囊泡膜的彈性確定時,其最終平衡位置隨雷諾數(shù)的變化而變化.然而,Alghalibi 等[30]發(fā)現(xiàn)軟顆粒直徑是管道直徑的0.2 倍時,顆粒的平衡位置與顆粒的柔軟程度或流速無關(guān),總是停留在管道中心.Kruger 等[31]發(fā)現(xiàn),當(dāng)雷諾數(shù)和毛細(xì)數(shù)足夠大時,軟粒子與軟粒子相互作用增強了其向通道中心的遷移.Hur 等[32]通過實驗發(fā)現(xiàn)不同彈性和黏度的顆粒和液滴具有獨立的橫向動態(tài)平衡位置,即在通道中粒子可以彼此分離.

盡管“軟”“硬”顆粒在微通道中的動力學(xué)行為受到了學(xué)者們的廣泛關(guān)注,但目前關(guān)于囊泡的研究缺少對膜黏性的考慮.例如,聚合物構(gòu)成的囊泡,依賴于聚合物網(wǎng)絡(luò)的交聯(lián)密度[33],可增加幾何性和機械性.此外,膜的黏彈性和收縮性取決于膜的成分和制備工藝[34].彈性囊泡忽略了膜黏度,其平衡位置由壁效應(yīng)和剪切梯度效應(yīng)之間的平衡決定[35].但是聚合物囊泡的黏性特征尤為顯著,勢必會影響囊泡在慣性遷移中的形變和平衡位置,因此,研究囊泡的黏彈性對其慣性遷移的影響具有重要意義.

本文采用基于流固耦合的有限元方法對黏彈性囊泡在微管流中的慣性遷移現(xiàn)象進行了系統(tǒng)研究,計算發(fā)現(xiàn),囊泡受到流體的擠壓會發(fā)生形變,從圓形變?yōu)闄E圓狀.當(dāng)雷諾數(shù)Re＜50 時,囊泡慣性遷移現(xiàn)象不明顯,但是隨著雷諾數(shù)的增加,囊泡慣性遷移的徑向位移越來越大,達到平衡位置后,保持一定的形狀平行于壁面移動.囊泡膜黏度的增大使得其慣性遷移后的平衡位置向管道中心偏移,但隨著囊泡阻塞比的增加,囊泡慣性遷移的平衡位置越來越接近管道邊緣.控制影響囊泡慣性遷移的不同因素可以有效地確定其最終平衡位置,有助于實現(xiàn)精準(zhǔn)分離.

2 仿真模型與有限元分析

如圖1 所示,在前期工作基礎(chǔ)上[36-42],通過二維模型計算囊泡在微管流中慣性遷移行為.當(dāng)流體流過囊泡時,會在囊泡表面施加一個垂直于流動方向的分力,稱為升力.升力分為壁面誘導(dǎo)升力和剪切梯度升力.剪切梯度升力指顆粒在管道流中由于流體速度分布引起的力,即當(dāng)流體繞過顆粒時,由于顆粒表面上下速度差所產(chǎn)生的力(通常指向通道壁面或者剪切率增加的方向).壁面誘導(dǎo)升力指顆?？拷诿鏁r,由于流體壓強會對顆粒產(chǎn)生遠離壁面方向的力,其隨著顆粒與壁面之間距離增加而減小[43].

圖1 囊泡受力分析圖(剪切梯度升力指向通道壁面,壁面誘導(dǎo)升力指向通道中心)Fig.1.Vesicle force analysis.(Shear gradient lift points to the channel wall,wall-induced lift points to the channel center).

圖2 是囊泡在微管流中慣性遷移的示意圖,管道內(nèi)流體假定為水,密度ρ為1000 kg/m3,黏度μ為0.001 Pa·s,且認(rèn)為是不可壓縮牛頓流體.水的流動采用湍流k-ε模型描述,矩形上下壁面為無滑動條件,左右兩邊為開放邊界,左邊為入口,初始流速為V,右邊為出口,由此形成微管內(nèi)水的流動.可變形囊泡假設(shè)具有黏彈性,置于管道內(nèi)部.圖2(b)—(d)展示了不同時刻囊泡周圍流速的分布以及囊泡自身的形變.

圖2 (a) 囊泡慣性遷移示意圖;(b)—(d) 不同時刻囊泡周圍流速圖 (管道寬為H=150 μm、長為D=1300 μm,囊泡半徑為a=20 μm,囊泡膜厚1 μm,囊泡內(nèi)外均為水.膜的楊氏模量為5000 Pa.管道入口速度為V,囊泡表面到管道壁面的距離為L)Fig.2.(a) Schematic representation of the vesicle inertial migration;(b)—(d) flow velocity around vesicles at different times (The channel width is H=150 μm and length is D=1300 μm.The vesicle radius is a=20 μm and the vesicle membrane is 1 μm thick.Water is both inside and outside the vesicles.The Young’s modulus of the membrane is 5000 Pa.The inlet speed is V,L is distance from the vesicle surface to the channel wall).

微管流中不可壓縮流的連續(xù)方程和運動方程為[36]

式中,F為外力;?表示(?/?x,?/?y);Ufluid為流速;I為單位矩陣;P為壓力.

假設(shè)囊泡膜具有黏彈性,其運動描述如下[42]:

其中ρsolid為膜密度,Usolid為位移,σ為應(yīng)力張量,Fv為單位體積的力.用Kelvin-Voigt 模型描述囊泡膜的黏彈性行為,其彈性應(yīng)變率和應(yīng)力張量的關(guān)系為[39]

式中η=0.022 Pa·s 為膜黏度;ε為膜應(yīng)變張量[39],

其中τ為松弛時間[39],

G為剪切模量[39],

E=5000 Pa 是楊氏模量,ν=0.33 是泊松比.

對于囊泡膜的變形與遷移,變形指數(shù)γ和水平速度vx、垂直速度vy分別為[39]

其中Lmax是囊泡質(zhì)心到囊泡表面的最大距離;S是囊泡膜包圍的面積;Usolidx和vx是囊泡的水平位移和速度,相對而言Usolidy和vy是囊泡的垂直位移和速度.

考慮到囊泡運動后形狀變化,其在不同時刻的傾向角θ變化可表示為[39]

式中,(x1,y1)表示囊泡質(zhì)心坐標(biāo);(x2,y2)表示囊泡膜上離質(zhì)心最遠點的坐標(biāo).

采用流固耦合計算流體和囊泡膜的相互作用,邊界處的相互作用可以表示為[36]

其中vw是囊泡膜的速度,n為流固界面的單位法向量,Γ為流體施加在囊泡邊界上的總力.因此,流體受到囊泡的反作用力為反方向.

體系雷諾數(shù)為Re=ρvd/u.根據(jù)以上參數(shù),本文雷諾數(shù)范圍為1＜Re＜400.在COMSOL Multiphysics 5.3a 軟件中,利用有限元方法來求解非結(jié)構(gòu)網(wǎng)格上的控制方程.采用任意拉格朗日-歐拉方法計算囊泡膜和流體網(wǎng)格節(jié)點的移動變化以及求解流固耦合有關(guān)方程[35].為了提高計算精度和收斂性,減少計算時間,采用自由三角網(wǎng)格來提高網(wǎng)格質(zhì)量并且對網(wǎng)格進行細(xì)化,且自動重新劃分網(wǎng)格以避免由流固相互作用導(dǎo)致的網(wǎng)格質(zhì)量問題[44,45].

3 結(jié)果及分析討論

3.1 阻塞比、雷諾數(shù)和初始位置對囊泡慣性遷移的影響

計算過程中,首先觀察到囊泡形狀的變化,從圓形變?yōu)榘朐滦?再由半月形變?yōu)闄E圓形,并保持橢圓狀以一定的傾向角平行于壁面移動,如圖2(b)—(d)所示,這與之前的研究結(jié)果一致[46-47].第一次形變由于流體擠壓囊泡促使其形變,且形變隨著雷諾數(shù)的增加而增大.第二次形變由囊泡第一次形變后,囊泡膜積累了彈性應(yīng)變能,促使囊泡擁有恢復(fù)原先形狀的能力.此外,囊泡運動過程中,由于管道流中囊泡表面上下速度差,使得囊泡向受力大的方向偏移,直到達到平衡(即垂直速度為零).

如圖3 所示,計算了囊泡初始質(zhì)心距壁距離40 μm 時的慣性遷移行為,其中囊泡阻塞比κ=a/h定義為囊泡半徑與管道寬度一半之比(ɑ為囊泡半徑,h為管道寬度一半).結(jié)果表明,不同阻塞比下,隨著雷諾數(shù)的增大,囊泡表現(xiàn)出完全不同的遷移行為.當(dāng)k=0.1 時,隨著雷諾數(shù)的增加,囊泡慣性遷移后的平衡位置逐漸趨于管道中心.然而,當(dāng)k=0.3 時,隨著雷諾數(shù)的增加,囊泡的平衡位置開始趨向管道中心,但隨著雷諾數(shù)的進一步增大,開始向管道邊緣靠近.此外,當(dāng)k=0.5 時也取得了相似的結(jié)果.因此,改變雷諾數(shù)或顆粒阻塞比都可以顯著改變囊泡的最終平衡位置.此外,Matas等[8]實驗報道中顆粒的阻塞比范圍約為0.0238—0.125 時,如圖3 虛線所示,隨著雷諾數(shù)的增加,顆粒偏向管道中心,這與本文軟球顆粒阻塞比為0.1時,囊泡慣性遷移平衡位置規(guī)律相似.

圖3 不同阻塞比下,雷諾數(shù)對慣性遷移平衡位置的影響(r 代表囊泡達到平衡位置后質(zhì)心的縱坐標(biāo).黑色虛線代表Matas 等[8]的實驗報道結(jié)果,其顆粒直徑為190 μm—1 mm,管道寬度為8 mm,即顆粒的阻塞比范圍為0.0238—0.125)Fig.3.Effect of Reynold numbers on the equilibrium position of inertial migration with different blocking ratios (r represents the ordinate of the centroid of the vesicle after reaching the equilibrium position.The black dashed line represents the experimentally reported results of Matas et al.[8]with particle diameters of 190 μm—1 mm and pipe widths of 8 mm,i.e.the particle blocking ratios ranged from 0.0238 to 0.125).

理論上,顆粒在微管流中慣性遷移是剪切梯度升力[19]和壁面誘導(dǎo)升力[20]兩者相互平衡的結(jié)果.為了能夠進一步闡釋不同阻塞比和雷諾數(shù)時囊泡平衡位置的差異,計算囊泡升力隨時間的變化圖.如圖4(a)所示: 當(dāng)雷諾數(shù)Re=100 時,隨著顆粒阻塞比的增加,囊泡受到的升力越大.因此,囊泡尺寸會改變囊泡受到的升力大小,進而影響平衡位置.如圖4(b)所示: 當(dāng)阻塞比k=0.3 時,囊泡受到的升力隨著雷諾數(shù)的增加而增加.因此,雷諾數(shù)也是改變囊泡平衡位置的重要因素.此外,無論初始升力存在多大差異,最終都趨于零,即囊泡達到平衡位置.其次,發(fā)現(xiàn)隨著雷諾數(shù)和顆粒阻塞比的增加,囊泡達到平衡的時間越久,這也為不同顆粒進行篩選和分離提供了時間窗口.

圖4 不同雷諾數(shù)和阻塞比下囊泡升力隨時間的變化圖(a) Re=100,k=0.1,0.3,0.5;(b) Re=50,100,250,k=0.3 (F 為升力,即是壁面誘導(dǎo)升力與剪切梯度升力的總和,方向為管道徑向)Fig.4.Plot of vesicle lift over time under different Reynolds numbers and blocking ratios: (a) Re=100,k=0.1,0.3,0.5;(b) Re=50,100,250,k=0.3 (Lift is the sum of wall-induced lift and shear gradient lift,in the tube radial).

進一步計算了囊泡初始位置對其慣性遷移規(guī)律的影響.計算結(jié)果發(fā)現(xiàn): 囊泡初始位置對其最終平衡位置影響顯著.如圖5 所示,阻塞比k=0.32、模量E=5000 Pa 時,不同初始位置導(dǎo)致囊泡慣性遷移后平衡位置發(fā)現(xiàn)明顯變化.當(dāng)L=51 μm時,即位于管道中心時,隨著雷諾數(shù)的變化,其平衡位置始終處于管道中心,沒有發(fā)生偏移.然而,當(dāng)初始位置L=31 μm 時,也就是處于管道一側(cè)時,結(jié)果表明囊泡隨著雷諾數(shù)的增加,呈現(xiàn)先趨于管道中心聚集而后開始向管道壁側(cè)遷移.隨著初始位置向管道壁側(cè)移動,這一現(xiàn)象更加明顯,例如L=1 μm 時,即初始時刻靠近管道壁,隨著雷諾數(shù)的增加,其最終平衡位置變化顯著.此外,發(fā)現(xiàn)當(dāng)L=11 μm 時,隨著雷諾數(shù)的增加,囊泡遷移后的平衡位置趨壁現(xiàn)象最明顯,甚至大于L=1 μm 的情況.因此,分析了不同初始位置時囊泡慣性遷移過程中所受力的變化.結(jié)果表明,囊泡初始位置距壁距離L=11 μm (L/H≈0.1,即約為管道寬度的1/10)時,產(chǎn)生的升力最大,從而導(dǎo)致囊泡慣性遷移的平衡位置更靠近管道壁面.

圖5 囊泡初始位置對慣性遷移的影響Fig.5.Effect of the initial vesicle position on the inertial migration.

為了驗證囊泡初始距壁距離為11 μm 時,升力最大這一觀點,描述了雷諾數(shù)為50 和200 時,囊泡的初始位置不同,其升力隨時間的變化圖,如圖6(a)和圖6(b)所示,結(jié)果表明囊泡初始距壁距離L=11 μm 時,升力最大.

圖6 不同初始位置時囊泡升力隨時間的變化 (a) Re=50;(b) Re=200Fig.6.Variations of vesicle lift with time at different initial positions: (a) Re=50;(b) Re=200.

以上結(jié)果表明,囊泡的初始位置、阻塞比和雷諾數(shù)是影響其平衡位置的關(guān)鍵因素,通過聯(lián)合調(diào)整相應(yīng)參數(shù),可實現(xiàn)對平衡位置的控制,進而實現(xiàn)精確篩選.

3.2 囊泡膜的模量、黏度和膜厚對慣性遷移的影響

囊泡膜的模量和黏度是表征膜特性的兩個最重要參數(shù),因此也勢必會影響其慣性遷移規(guī)律.模量決定著囊泡膜的變形難易程度,因此,選取囊泡膜初始位置距壁距離L=16 μm,阻塞比k=0.32,膜厚1 μm 進行分析計算,如圖7 所示,隨著模量的增加,囊泡平衡位置幾乎沒有任何變化.然而,也注意到當(dāng)Re=350 時,不同模量的囊泡的平衡位置有了差異.但由于模型收斂性的限制,沒能進一步計算更大雷諾數(shù)的情況.

圖7 囊泡膜模量對慣性遷移的影響Fig.7.Effect of the modulus of vesicle membrane on the inertial migration.

對于許多囊泡,尤其是高分子類的囊泡,膜的黏性是一個無法忽略的屬性.已有研究表明,利用高分子囊泡模擬細(xì)胞的動力學(xué)過程中,最重要的一點是細(xì)胞膜的黏度,增加膜的黏度可以減弱軟顆?；蛏锛?xì)胞的變形[48].分析了囊泡形變及其慣性遷移的軌跡,囊泡形變越大表示其受力越大,促使其平衡位置朝向管道壁面移動.圖8 描述了不同黏度η時,囊泡形變及其遷移的軌跡(補充文件圖S1和圖S2 (online)給出了不同膜量和模厚的計算結(jié)果).在圖8 中描述了黏度η=1 和100 Pa·s 時,囊泡慣性遷移過程的形變和平衡位置軌跡圖以及囊泡所受最大升力.從圖8 可以明顯地看出,增大膜黏度可以抵抗因受流體作用力而發(fā)生的形變,并使得囊泡升力減小導(dǎo)致其平衡位置降低.

圖8 不同膜黏度時囊泡的形變及其遷移位置軌跡圖Fig.8.Deformation of vesicles and their migration location trajectories at different membrane viscosities.

其次,計算了模量E=10000 Pa,阻塞比為0.32,囊泡膜初始距壁距離L=16 μm 時,不同膜黏度η對囊泡慣性遷移的影響.如圖9 所示,在所計算的黏度和雷諾數(shù)區(qū)間,囊泡慣性遷移后的平衡位置向管道中央偏移.這主要可歸因于黏度會改變囊泡在流體中的形變,進而影響受力情況,從而出現(xiàn)了規(guī)律性偏向管道中心的現(xiàn)象.

圖9 囊泡膜的黏度對慣性遷移的影響Fig.9.Effect of vesicle membrane viscosity on the inertial migration.

最后,還計算了模量5000 Pa,阻塞比為0.32,囊泡膜初始距壁距離L=16 μm 時,不同膜厚度I對囊泡慣性遷移的影響.如圖10 所示,在計算的雷諾數(shù)區(qū)間,觀測到隨著囊泡膜厚度的增加,囊泡慣性遷移后的平衡位置趨向管道中心.

圖10 囊泡膜厚對慣性遷移的影響Fig.10.Effects of vesicle membrane thickness on inertial migration.

4 結(jié)論

高分子囊泡在微管流中慣性遷移現(xiàn)象的規(guī)律是應(yīng)用慣性微流控技術(shù)實現(xiàn)其分離的重要基礎(chǔ).本文采用基于流固耦合的有限元方法對高分子囊泡在微管流中慣性遷移現(xiàn)象進行了系統(tǒng)研究,分析了囊泡因受到流體作用力而形變并發(fā)生慣性遷移現(xiàn)象的機理.研究表明: 隨著雷諾數(shù)的增大,囊泡慣性遷移的平衡位置離其初始位置越來越遠;隨著阻塞比的增加,囊泡慣性遷移后的平衡位置越來越接近壁面.此外,模擬還發(fā)現(xiàn)當(dāng)初始位置約為管道寬度的1/10 時,囊泡受到的升力最大,從而導(dǎo)致囊泡慣性遷移的平衡位置更靠近管道壁面.對于囊泡膜的模量和黏度,結(jié)果表明模量和黏度決定了囊泡的變形程度,模量對囊泡平衡位置影響較小,但增大黏度和膜厚會促進囊泡的平衡位置偏向管道中心.本研究有助于進一步明晰囊泡在慣性遷移過程中的形變和平衡位置,為囊泡在藥物輸運、化學(xué)反應(yīng)和生理過程的應(yīng)用提供可靠的計算依據(jù).

然而也注意到,由于本文采用的模型在大雷諾數(shù)時收斂性不佳,導(dǎo)致本文計算的囊泡膜的模量和黏度變化區(qū)間,尤其是雷諾數(shù)大小受到了限制.如何進一步提高模型的收斂性和計算效率,以便更好地考慮更大參數(shù)范圍內(nèi)囊泡的慣性遷移規(guī)律,是未來工作的重點.