DLD裝置中剛性顆粒運動的直接數(shù)值模擬研究*

喻文廣，余釗圣，葉尚軍，邵雪明，王博康

(浙江大學力學系，浙江 杭州310027)

0 引 言

顆粒分離技術(shù)是生物、醫(yī)學檢測等領(lǐng)域中的一項關(guān)鍵技術(shù)。確定性側(cè)向位移分離技術(shù)(DLD)是一種基于粒子尺寸實現(xiàn)高效、連續(xù)分離顆粒的技術(shù)。

DLD分離技術(shù)自2004年被 Huang等人［1］提出后，由于其具有高效、高分辨率、連續(xù)分離等優(yōu)點，已受到廣泛的關(guān)注，目前文獻中已有不少的實驗和應用研究報道［2-9］。關(guān)于DLD分離技術(shù)的數(shù)值研究工作還相對缺乏。Loutherback等人［10］采用COMSOL軟件求解了無顆粒存在時的單相流場，然后根據(jù)速度剖面和前述近似理論計算了臨界半徑，并采用該方法研究了三角形微柱的形狀和尺寸對分離臨界半徑的影響。AIFandi等人［11］采用COMSOL軟件求解了菱形和翼形微柱陣列的流場，然后研究了球形顆粒在其中的分離特性(不考慮顆粒對流場的影響)。Quek等人［12］采用浸沒邊界法數(shù)值模擬了彈性顆粒在圓柱陣列中的運動。

至今為止，關(guān)于DLD裝置中剛性顆粒的運動和臨界半徑的預測，都基于前述的近似理論。Inglis等人［13］通過實驗研究了圓柱陣列中剛性顆粒的臨界尺寸，發(fā)現(xiàn)前述近似理論的臨界半徑比實驗值要小。本研究首次采用高精度直接數(shù)值模擬方法對剛性顆粒的運動和臨界半徑進行研究。另外，由于實際DLD裝置中側(cè)向圓柱數(shù)目較多(上百或者更多)，數(shù)值模擬如果采用與實際相同數(shù)量的柱列數(shù)目則計算量太大，實際的計算只能取有限的側(cè)向區(qū)域大小，側(cè)向的邊界為人工邊界，對于人工邊界如何給邊界條件是非常重要的問題。本研究將首先研究該邊界條件問題。

本研究首先通過計算采用側(cè)向采用固壁邊界條件時不同數(shù)量圓柱的流場和周期性邊界下流場，并對流場傾斜角度進行分析，以論證側(cè)向邊界條件;然后，對剛性顆粒在DLD中的運動情況進行模擬分析其臨界分離直徑，并與實驗結(jié)果進行對比，驗證了該方法的精確性。

1 數(shù)值模型

1.1 虛擬區(qū)域方法

針對流固兩相流動的典型數(shù)值模擬方法［14-15］有雙流體模型、拉格朗日點粒子模型和直接數(shù)值模擬。在雙流體模型中，固體顆粒被當成連續(xù)相處理，而后兩個模型中顆粒是離散的，以拉格朗日方法來追蹤軌跡。在點粒子模型中，顆粒和流體的相互作用力是按經(jīng)驗公式給的，對于顆粒很小并且濃度很稀的情況，這樣的模型是合適的。在針對流固兩相流的直接數(shù)值模擬方法中，考慮了顆粒的體積以及顆粒和流體界面上的運動學和動力學邊界條件，是一種沒有引入近似模型的方法，所以這種方法被稱為“直接數(shù)值模擬”。本研究采用的直接力/虛擬區(qū)域(Direct Forcing/Fictitious Domain，簡寫為DF/FD)方法［16］屬于直接數(shù)值模擬方法。DF/FD方法是通過改進最早由Glowinski等［17］提出的分布式拉格朗日乘子/虛擬區(qū)域(DLM/FD)方法得到的。該方法的中心思想是假設固體顆粒的內(nèi)部也充滿流體，然后通過虛擬體積力(即拉格朗日乘子)來使得這部分流體的運動符合剛體顆粒的運動。這樣，流場區(qū)域就變?yōu)楹唵螀^(qū)域，可以采用固定的結(jié)構(gòu)化網(wǎng)格來求解。

為了描述上的簡要，以下公式中只考慮一個顆粒情況。假設顆粒的密度、體積、轉(zhuǎn)動慣量、速度和角速度分別為 ρs，Vp，J，U 和 ωp。用 P(t)和 Ω 分別表示顆粒所占的區(qū)域和整個計算區(qū)域。采用下述特征量來無量綱化方程:特征長度Lc，特征速度Uc，特征時間Lc/Uc，特征壓力，特征擬體力。無量綱化后的虛擬區(qū)域法控制方程為:

式中:u—流體的速度，p—流體的壓力，λ—定義在粒子區(qū)域的擬體力(拉格朗日乘子)，r—顆粒中心為原點的位置矢徑，ρr—顆粒密度 ρs與流體密度 ρf的比值。

ρr定義為:ρr= ρs/ρf，雷諾數(shù)定義為:Re= ρfUcLc/μ，弗勞德數(shù)定義為，無量綱的體積定義為:，無量綱的轉(zhuǎn)動慣量定義為:J*=J/。

數(shù)值離散格式的主要特點是:采用算子分裂格式將耦合方程組式(1～5)分解為兩個子問題:流體子問題和拉格朗日乘子子問題。流體子問題采用基于投影格式的二階有限差分方法求解;拉格朗日乘子問題采用直接力格式求解。

1.2 控制參數(shù)

DLD裝置的示意圖如圖1所示。該裝置由周期性布置的微柱列構(gòu)成，假定圓柱間距為S，沿著流動方向每排豎直柱列與前一排柱列側(cè)向錯開一定位移量ΔS，則柱列的周期為N=S/ΔS，傾斜柱列的斜率為ε=1/N。DLD分離技術(shù)的原理(近似理論)如下:當流體流過周期為N的柱列時，豎直方向兩相鄰圓柱間的流體可被分成N條具有相同流量的窄流體束如圖1所示。流體束的相對位置(或序號)在流經(jīng)每列圓柱時會發(fā)生改變，但在經(jīng)過N列圓柱之后又回其原來位置。假定與圓柱相鄰的流束的寬度為Rc，當顆粒的半徑小于Rc時，可以近似認為顆粒(中心)跟隨流體在流體束內(nèi)部運動，呈現(xiàn)“之”字形(zigzag)運動軌跡;當顆粒半徑大于Rc時，如顆粒中心位于第2條流束，則可以近似認為顆粒在流經(jīng)每列圓柱時都會處于兩圓柱間的第2條流束位置，因此顆粒會沿著傾斜通道一直運動，導致側(cè)向位移。根據(jù)不同尺寸的粒子運動軌跡不同的特性，就可以實現(xiàn)分離。在實際應用過程中，可以將不同臨界半徑的柱列組合起來，實現(xiàn)更大尺寸范圍的粒子的分離。

本研究只考慮固體和流體密度相等的情況，即ρr=1。也不考慮重力的作用，即Fr=0。不考慮顆粒的布朗運動，在控制方程式(1～5)的推導中已經(jīng)做了忽略。

如圖1所示，(初始)顆粒直徑為d(半徑為a)，固定圓柱的直徑為D，兩個方向的圓柱中心距均為S，間隙為Lg，柱列的周期為N，斜率為ε=1/N。本研究只考慮D=Lg=S/2的情形。

圖1 DLD裝置示意圖(周期N=3)

取均勻來流速度U0為特征速度，圓柱直徑D為特征長度。因而:

對于實際DLD裝置中的流動，雷諾數(shù)比較低，流體慣性可以忽略。本研究取Re=0.1，驗證算例表明在Re＜1時，雷諾數(shù)取值不影響顆粒的運動。

2 邊界條件和數(shù)值參數(shù)

關(guān)于邊界條件，流向入口為均勻來流，出口為無反射出流條件，側(cè)向取周期性邊界條件，計算域以及邊界條件如圖2所示。計算區(qū)域大小為56D×28D，圓柱數(shù)量為24×14。流場網(wǎng)格數(shù)為2 048×1 024，一個圓柱直徑覆蓋約37個網(wǎng)格。時間步長為0.00 025D/U0。對于d/D=0.6的較大顆粒，驗證算例表明:把上面給定的網(wǎng)格尺寸和時間步長加大一倍，對結(jié)果的影響很小。限于計算條件，本研究所研究的顆粒尺寸d＞0.2D，圓柱陣列周期N≤7。

圖2 計算域以及邊界條件

對于側(cè)向邊界條件，實際的DLD裝置在上、下兩側(cè)為固壁，且側(cè)向的圓柱數(shù)量較多(如數(shù)百個)，如果采用實際的邊界條件和圓柱數(shù)量，計算量很大。周期性邊界條件相當于側(cè)向的圓柱數(shù)量無窮多，如果采用周期性邊界計算量將大大減小。但是周期性邊界與實際裝置的固壁邊界條件又會有所差異，下面就該問題進行探討。

本研究分別模擬了上下邊界采用固壁邊界條件，圓柱數(shù)目分別為:24×14、24×28、24×57、24×115的柱列無顆粒流動，和上、下邊界采用周期性邊界條件時柱列數(shù)量為(24×28)柱列無顆粒流動。

通過分析發(fā)現(xiàn)，DLD中流場并非完全像上文提到的近似理論那樣:假定DLD裝置中的流線是周期性的，其周期等于圓柱陣列周期，因而平均意義上的流線是水平的，沒有側(cè)向偏移。事實上，由于傾斜通道的存在，流線有總體的傾斜，且傾斜角度隨上、下壁面間距增大而增加(傾斜角度α定義為:與流線相切的直線與水平方向的夾角，傾斜流線圖如圖3所示)。靠近上壁面處的流體傾斜角度較小，離壁面越遠傾斜角度增大。不同側(cè)向位置的流線傾斜角度與位置關(guān)系如圖4所示，其中橫坐標為無量綱化側(cè)向位置，縱坐標為不同位置流線的傾斜角度。從圖4可以看出，在靠近壁面處由于壁面的存在使得流線傾斜較小，但隨著離遠離壁面增大，流線的傾斜角度迅速增大，當超過壁面距離一段距離后增長平緩，在中間區(qū)域內(nèi)緩慢變化并達到最大值。因為在中間區(qū)域，壁面影響作用最小，傾斜度最大，本研究將選區(qū)中間區(qū)域的流場傾斜角度來分析邊界的影響。

圖3 周期N=3、上下邊界為固壁邊界條件、柱列數(shù)目為24×115的無顆粒流場中的傾斜流線

圖4 不同側(cè)向位置處流線的傾斜角度

下面將分析不同側(cè)向柱列數(shù)量的分離裝置中間區(qū)域流場的傾斜程度，并以此來表征側(cè)向柱列數(shù)目以及邊界條件對流場的影響。

模擬結(jié)果如圖5所示，從圖5可以看出，當側(cè)向柱列數(shù)目為14時流線傾斜角度接近于0°，當隨著側(cè)向圓柱數(shù)量的增多，流線傾斜角度越來越大。側(cè)向圓柱柱列數(shù)目為115時，其偏移角度α≈3.4°。而當采用周期性邊界條件下(N=3)的流線傾斜角度α≈3.7°，兩者非常接近。當圓柱數(shù)目繼續(xù)增加時，傾斜角度會越來越接近周期性邊界條件的計算結(jié)果。因此，當DLD裝置中側(cè)向圓柱數(shù)目大于100時，采用周期性邊界條件是合理的。

圖5 固壁邊界下不同側(cè)向柱列數(shù)目對傾斜角度的影響

3 剛性顆粒在DLD中的運動

Inglis等人測量了剛性顆粒在DLD分離裝置中臨界直徑，其結(jié)果表明近似理論給出的臨界直徑顯著低于實驗值，而至今尚未有采用直接模擬方法計算剛性顆粒臨界直徑的報道。本研究將計算剛性顆粒的臨界直徑，以驗證虛擬區(qū)域方法的可靠性。

在模擬中，較小的顆粒會沿著流線運動，當顆粒直徑接近臨界直徑時，顆粒相鄰兩次下拐的間距會隨著顆粒直徑增加而增加，直至沒有下拐。本研究定義顆粒在計算區(qū)域內(nèi)沒有出現(xiàn)下拐為側(cè)向遷移模態(tài)，否則皆為“之”字形模態(tài)。兩種運動模態(tài)如圖6所示，由圖6可知，小于臨界尺寸的顆粒沿“之”字形路線運動，大于臨界尺寸的顆粒沿傾斜通道做側(cè)向遷移運動。

圖6 不同尺寸顆粒的運動模態(tài)

圖7 DLD裝置中剛性顆粒臨界直徑數(shù)值結(jié)果和實驗結(jié)果的對比

本研究把計算得到的臨界半徑和Inglis等人的實驗結(jié)果做了對比，其結(jié)果如圖7所示。圖7中橫坐標為柱列的斜率ε，即柱列周期的倒數(shù)(ε=1/N);縱坐標為顆粒直徑與圓柱間距的比值。從圖7可看出，數(shù)值結(jié)果和實驗結(jié)果吻合較好。需要說明的是:圖中給出的實驗擬合曲線主要根據(jù)ε＜0.2(N＞5)的數(shù)據(jù)，對于ε＞0.2(N＜5)，側(cè)向遷移模態(tài)的數(shù)據(jù)很少，因而無法確定較精確的臨界半徑。對于周期3(ε=1/3)陣列，模擬中沒有觀察到側(cè)向遷移模態(tài)(d/D≤0.98)，而ε≈0.32的實驗數(shù)據(jù)也沒有顯示存在側(cè)向遷移模態(tài)。

4 結(jié)束語

本研究用虛擬區(qū)域法對剛性顆粒在DLD裝置中的運動進行模擬，并分析了臨界半徑，得到了與Inglis的實驗結(jié)果非常貼近的臨界直徑分布，驗證了該方法在計算該問題時的準確性。此外，還對側(cè)向邊界條件問題進行了研究，其結(jié)果顯示:當上、下邊界采用周期性邊界條件時的模擬結(jié)果，與上、下邊界采用固壁邊界的側(cè)向數(shù)目較多(N＞100)柱列的計算結(jié)果非常接近。所以對于模擬數(shù)量較多的實際DLD裝置，可以采用周期性邊界條件簡化計算量。

在驗證了該方法的準確性之后，筆者將在以后的工作中模擬彈性顆粒在DLD中的運動，研究分析顆粒的變形特性、初始形狀以及初始取向等因素對臨界半徑的影響。

［1］HUANG L R，COX E C，AUSTIN R H，et al.Continuous particle separation through deterministic lateral displacement［J］.Science，2004，304(5673):987-990.

［2］LOUTHERBACK K，PUCHALLA J，AUSTIN R H，et al.Deterministic microfluidic ratchet［J］.Physical Review Letter，2009，102(4):45301-45305.

［3］BEECH J P，HOLM S H，ADOLFSSON K，et al.Sorting cells by size，shape and deformability［J］.Lab on a Chip，2012，12(6):1048-1051.

［4］DAVIS J A，INGLIS D W，MORTON K J，et al.Deterministic hydrodynamics:taking blood apart［J］.Proceedings of the National Academy of Sciences of the United States of America，2006，103(40):14779-14784.

［5］NAGRATH S，SEQUIST L V，MAHESWARAN S，et al，I-solation of rare circulating tumour cells in cancer patients by microchip technology［J］.Nature，2007，450(7173):1235-1239.

［6］INGLIS D W，MORTON K J，DAVIS J A，et al.Microfluidic device for label-free measurement of platelet activation［J］.Lab on a Chip，2008，8(6):925-931.

［7］HOLM S H，BEECH J P，BARRETT M P，et al.Separation of parasites from human blood using deterministic lateral displacement［J］.Lab on a Chip，2011，11(7):1326-1332.

［8］FRECHETTE J，DRAZER G.Directional locking and deterministic separation in periodic arrays［J］.Journal of Fluid Mechanics，2009(627):379-401.

［9］BALVIN M，SOHN E，IRACK T，et al.Directional locking and the role of irreversible interactions in deterministic hydrodynamics separations in microfluidic devices［J］.Physical Review Letters，2009，103(7):78301-78304.

［10］LOUTHERBACK K，CHOU K，NEWMAN J，et al.Improved performance of deterministic lateral displacement arrays with triangular posts［J］.Microfluidics and Nanofluidics，2010，9(6):1143-1149.

［11］Al-FANDI M，Al-ROUSAN M，JARADAT M，et al.New design for the separation of microorganisms using microfluidic deterministic lateral displacement［J］.Robotics and Computer-Integrated Manufacturing，2011，27(2):237-244.

［12］QUEK R，LE D V，CHIAM K H.Separation of deformable particles in deterministic lateral displacement devices［J］.Physical Review E，2011，83(5):56301-56307.

［13］INGLIS D W，DAVIS J A，AUSTIN R H，et al.Critical particle size for fractionation by deterministic lateral displacement［J］.Lab on a Chip，2006，6(5):655-658.

［14］解婧陶.板料成形數(shù)值模擬中回彈問題的處理［J］.機電技術(shù)，2011(6):91-93.

［15］楊 陽，趙建平.EPS聚合反應釜內(nèi)混合過程的數(shù)值模擬［J］.輕工機械，2012，30(6):1-4.

［16］YU Z，SHAO X.A direct-forcing fictitious domain method for particulate flows［J］.Journal of Computational Physics，2007，227(1):292-314.

［17］GLOWINSKI R，PAN T W，HESLA T I，et al.A distributed Lagrange multiplier/fictitious domain method for particulate flows［J］.International Journal of Multiphase Flow，1999，25(5):755-794.