確定性橫向遷移中液固兩相流的三維格子玻爾茲曼模擬*

彭浩宇，陳英懷，郭鐘寧，黃志剛

（廣東工業大學機電工程學院，廣東廣州 510006）

0 引言

利用微流控芯片對微觀物質進行快速分離是一種前景廣闊的新興技術，可用于生物醫療、環境檢測和精細化工等科技領域[1]。近十幾年發展起來的確定性橫向遷移技術（Deterministic Lateral Displacement，DLD）采用微柱陣列對微觀顆粒進行分選，由于其清潔、簡便和可連續分離等優點，在微流控分選技術中占據重要位置[2]。目前，國內外對DLD已經開展了比較廣泛的研究工作，然而關于流體和顆粒在芯片微通道中的具體運動情況，仍需深入研究。

目前關于DLD的研究主要集中于微柱形狀、排列、粒子類型和邊界的影響等方面。早期的DLD裝置多采用圓形微柱陣列，近期有學者研究了多邊形、流線形[3]甚至“工”字形[4]等微柱陣列的粒子分選特點。傳統DLD中的微柱多采用規則的斜向陣列，也有研究人員采用不規則排布[5]，達到特別的粒子分選目標。在實際使用中，DLD微柱陣列封裝于微流控芯片之中，流道邊界的影響需要加以考慮并適合處理[6]，以提升分離效果。關于粒子類型，除傳統的剛性粒子外，近來有不少研究關注細胞和病毒等生物粒子的分離[7]。由于生物粒子柔性大、可變形和易受損傷等特點，使其分離過程中存在一些獨特的現象，有學者專門對此進行了設計[8]，保證分離的效率和生物使用性。

確定性橫向遷移裝置的單元流道尺寸微小，流體速度相對較低，早期有學者根據層流思想，提出了其中的粒子分離理論[9]。然而該理論并未考慮粒子與流體之間的相互作用（多相流），也忽略了邊界的影響，可能使理論預測結果出現較大偏差。針對這一問題，隨后研究人員根據實驗結果對分離理論進行了修正，提出了多種經驗模型[10，11]。另一方面，則有學者利用數值仿真技術對DLD的分離機制進行了深入研究。普遍的做法是采用有限單元法或有限體積法，對DLD微柱陣列之間的流體動力學進行計算[7，12]。為探討粒子與流體的相互作用，Quek等人采用浸入邊界法研究變形粒子在DLD中的分離情況[13]。Krueger等則采用浸入邊界-格子玻爾茲曼法，研究了變形紅細胞在DLD中的運動[14]；國內韋建輝等采用類似方法，對變形粒子在圓形、多邊形及“工”字形微柱中的運動進行了分析[15]。喻文廣等采用直接數值模擬方法對DLD中剛性顆粒的運動進行了計算[16]。數值方法對DLD中流體行為的研究取得了很好的研究成果，不過目前的計算模型多為二維模型，很少考慮芯片頂底邊界的影響。關于粒子與流體之間的相互作用，也仍需要繼續深入探討。

格子玻爾茲曼法（Lattice Boltzmann Method，LBM），從介觀層面上描述流體狀態，在計算上基于粒子描述，特別適合處理多相微流體和復雜邊界等問題[17]。本文將采用三維格子玻爾茲曼方法對DLD裝置中的微流體和粒子運動進行分析計算，分析不同形狀和排列的微柱間的流體運動，研究邊界對流場的影響，探討粒子分離和阻塞等的詳細機理。

1 計算模型

1.1 確定性橫向遷移基本原理

確定性橫向遷移的基本原理如圖1所示，分選裝置由規則排列的微柱組成，每列微柱與前一列稍微錯開，形成斜向排布的陣列。假定陣列的周期為N，微柱間距為λ，則陣列的傾斜度為ε=1/N，相鄰列的橫向偏移量為ε?λ。流體在微柱間流動時，自然形成傾斜的主流道和橫向的分流道。考察縱向的流動，由于微柱的“柵欄”效果，流體被分成N條具有相同流量的流束，其寬度為RC。當粒子尺寸顯著大于RC時，微柱會連續將其推向相鄰的流束，從而使粒子沿傾斜的主流道運動，稱為“L”形軌跡。而當粒子尺寸小于RC時，將始終保持在初始流束內，呈現繞微柱而過的“Z”字形軌跡，其總體沿縱向運動。由此，粒子的運動軌跡按臨界半徑RC區分，可實現不同尺寸粒子的分離[9]。在實際使用中，可以將具有不同臨界半徑的微柱陣列進行組合，實現更多類別的粒子分離。

圖1 確定性橫向遷移的粒子分離原理（以周期N＝3示例）

1.2 格子玻爾茲曼模型

格子玻爾茲曼方法（LBM）根據流體微觀分布函數的變化，求解流體動力學問題。模型中將將流體離散為一系列的微團，微團以特定方式在規則的格子上碰撞和遷移，如圖2所示，對所有格點上的微團運動進行統計，便可得出流體的宏觀運動狀態。微團的碰撞運動可根據微觀分布函數的變化進行描述[18]，如下式：

式中： fi表示微團i的微觀分布函數；x表示格子上的格點；c表示流體微團的離散速度集合；Δt表示離散時間步長；Ω表示碰撞矩陣； feq表示平衡態分布函數。

碰撞方程可在動量空間中，采用碰撞和遷移兩個基本步驟進行數值求解。

圖2 格子玻爾茲曼法計算原理

由基本原理可知，LBM是的一種基于拉格朗日粒子描述的計算方法，特別適合多相流動和復雜界面流動的模擬。在LBM的計算模型中，格子的安排有多種方式，本文進行DLD中多相微流動的三維計算，因此選擇了三維的D3Q27格子，每個微團可以向鄰近的27個格子進行遷移，具有很高的自由度。

關于DLD中固相粒子運動的計算，本文采用的是離散相（Discrete Phase Model，DPM）方法，粒子表示為離散的小球，以拉格朗日方法追蹤軌跡，可以對粒子的密度、半徑、質量、與邊界的碰撞方式等參數進行定義。粒子運動方程如下：

式中：up、 ρp，u和 ρ分別為粒子的運動速度、密度和流體的速度與密度；g為重力；FD為曳力，與粒子直徑、流體粘度和雷諾數等參數有關；FA為外力。

1.3 仿真模型及參數

確定性橫向遷移的仿真計算模型如圖3所示，高度方向（Z軸）與橫向（Y軸）兩端設為固定邊界，橫向（X軸）入口設置流速邊界條件，速度0.1 m/s?～0.5m/s變化。出口設為壓力出口，表壓為0。對兩種形狀的微柱陣列進行了對比計算：“圓柱形”和“鵝蛋形”微柱，其幾何外形見圖4。圓形微柱的主要參數為：直徑50μm，主流道寬度a1=25μm，副流道寬度b1=25μm，陣列周期：N=19；“鵝蛋形”微柱的主要參數：短軸×長軸17×24μm，主流道寬度 a2=32μm，副流道寬度b2=15μm，陣列周期：N=59。微柱高度均為50μm，其他參數詳見文獻[7][19]，選擇這組參數方便與文獻中的實驗數據進行對比，計算流域包含一個完整的微柱陣列周期。

圖3 DLD的仿真計算模型

圖4 微柱的幾何形狀及陣列結構

采用均勻格子玻爾茲曼網格，單元的尺寸為2μm，時間步長采用自動調節模式，courant數設置為0.001，采用多松馳方法進行碰撞計算。初步計算表明，采用這些參數，可使穩定性參數維持在1以下，計算結果穩定可靠。臨界偏轉半徑理論值近似為RC～ε?a，固相粒子的半徑以此為基準，在0.25ε?a～4ε?a的范圍內變化，以精確界定真實臨界偏轉半徑。基礎流體設定為液體水，密度為1 000 kg/m3，粘度為0.001Pa·s，每次在微柱陣列入口處釋放不同半徑的固相粒子，其中固相粒子密度為1 050 kg/m3。

2 仿真結果及討論

2.1 微柱間的流場分布

對DLD裝置中的流場分布進行了分析，圓形微柱間的局部典型流場如圖5所示，鵝蛋形微柱的結果類似。從圖5（a）水平方向的流場可以看出流體在主流道中的流速最快，粒子將主要沿主流道運動。而副流道中的流速相對較低，而且總體沿橫向方向。粒子進行橫向遷移，走“Z”形軌跡，必須通過副流道。由于邊界的存在，邊界附近的微柱后方出現了渦流，有可能使流入其中的粒子發生阻塞。從圖5（b）高度方向的流場看到，沿流道的頂邊和底邊也存在低流速區。粒子主要在中部高流速區運動，少數進入邊界附近的粒子有可能發生阻塞。從三維仿真的結果得出結論，不僅是兩側邊界，還有頂底邊界，都可能對DLD的粒子收集率產生影響。為提高收集率，在實際芯片設計中，可以在兩側邊界處開設輔助流道或者增加柱列的總寬度；而為降低頂底邊界的影響，可以適當提高微柱的高度，增加有效工作區。

圖5 圓形微柱DLD裝置中的典型流場（入流速度u=0.5m/s）

相鄰兩個微柱之間的流場直接決定了粒子的偏轉情況，本文對圓形和鵝蛋形相鄰微柱間的流速進行了對比，結果如圖6所示。速度與坐標都經過了規整化處理，u=u/umax，y=y/w，w為微柱陣列總寬。可以看到，當采用圓形柱列時，兩側邊界對流場均勻性的影響更為顯著，邊界附近微柱間的流速相比中心流速降低率達到15.2%。兩鵝蛋形柱列，由于其定向導流作用，流場更加均勻，邊界流速降低率低于1%。此外，鵝蛋形微柱間的流場并非對稱分布，而是向大端方向傾斜，從而使粒子能更順利地通過微柱間的流道和副流道，起到減少堵塞的作用。本文的仿真結果證實了非圓微柱對改善DLD阻塞情況的作用，與文獻[20]的結論相符。更重要的是，首次揭示鵝蛋形（或三角形，通常鵝蛋形微柱可認為是帶有過渡圓弧的三角形微柱，在實際加工中由于銳邊加工的困難，實際制備的三角形微柱最終都類似于蛋形）對改善流場的均勻性具有顯著效果，可以提高DLD的有效工作區，對DLD裝置的優化設計具有指導作用。

圖6 圓形和鵝蛋形微柱間流速的對比

2.2 臨界偏轉半徑

對粒子在圓形和鵝蛋形微柱陣列中的粒子運動軌跡進行追蹤，典型結果如圖7、8所示。由理論公式估算，圓形和鵝蛋形微柱陣列中的粒子偏轉半徑分別約為3μm和2μm，因此計算了半徑在1～6μm （間隔0.5μm）之間的粒子的運動情況，分析了不同尺寸粒子的運動軌跡。在圓形DLD中，發現當粒子半徑小于3.5μm時，運動軌跡為Z模式，而粒徑大于4μm的粒子運動模式為L模式。為更精確界定臨界偏轉半徑RC，進一步細分粒子直徑，發現RC約為3.5μm。

而在鵝蛋形微柱陣列中，也計算了半徑在1～6μm（間隔0.5μm）之間的粒子的運動情況，發現當粒子半徑小于2μm時，運動軌跡為Z模式，而粒徑大于2.5μm時，粒子運動模式為L模式。進一步細分發現RC約為2.4μm。本文數值計算的臨界偏轉半徑與文獻[7][19]中的研究結果相符。

圖7 圓柱形微柱陣列粒子運動軌跡

圖8 鵝蛋形微柱陣列粒子運動軌跡

為對比圓形與鵝蛋形微柱陣列的粒子通過性，計算了各自的主流道寬度－臨界偏轉半徑比率 γ?=?a?/?Rc。圓形陣列 γ?=?7.1，而鵝蛋形陣列 γ?=?13.3，接近為圓形陣列的2倍。以上結果表明：為達到相同的偏轉效果，鵝蛋形陣列允許采用更寬的主流道，而寬流道意味著更好的粒子通過性，可顯著改善粒子阻塞現象。

本文還研究了速度對臨界偏轉半徑的影響，發現不同速度下的臨界偏轉半徑幾乎沒有變化。其原因是由于微流道的特征尺寸很小（10-6m量級），在通常工作流速下（10-1m/s量級），由雷諾數計算公式可知其雷諾數為10-1量級，遠小于2 300，故其流動呈層流狀態，流場特性沒有隨著速度變化發生明顯變化，因此臨界偏轉半徑不變。

3 結論

確定性橫向遷移（DLD）是一種應用廣泛的粒子分選技術，本文采用三維格子玻爾茲曼法，對DLD中的流場特性和粒子臨界偏轉半徑進行了研究。對比了圓柱形和鵝蛋形DLD中的流場特性，仿真發現非圓微柱能改善DLD的阻塞，與之前文獻中的結論相符。相較于圓形陣列，鵝蛋形陣列在改善流場的均勻性方面具有顯著效果，可以提高DLD的有效工作區。仿真還計算了粒子的臨界偏轉半徑，計算結果與前人的實驗結果一致。這些結論對DLD裝置的優化設計具有指導作用。