基于微流控芯片模型的滲流實驗與數值模擬

聶紹凱，劉鵬飛，巴特，陳云敏

(浙江大學 巖土工程研究所，浙江 杭州，310058)

多孔介質滲流現象廣泛存在于自然界及工程應用中，如土壤中的地下水滲流、土石壩工程的穩定性、石油天然氣的開采等[1-2].長期以來，許多學者對多孔介質的滲流特性開展大量的研究，傳統的多孔介質滲流實驗往往基于宏觀尺度的土柱實驗和箱型模型試驗等[3-4].最早描述滲流規律的公式是由法國工程師Darcy于1856年通過實驗總結得到的Darcy模型，后期不斷有學者對Darcy模型進行修正，得到其他一些半經驗模型，如Darcy-Forchheimer修正模型[5]等.在理論方面，對于多孔介質滲流特性的研究將多孔介質的微細觀結構參數與滲流行為關聯起來，得到相應的解析解，目前在實際工程中應用較為廣泛的是Kozeny-Carman方程[6].現有的研究表明，影響多孔介質滲流特性的主要因素有輸送流體的物性、孔隙率、孔隙連通性、顆粒形狀、比表面積等[7-8].由于多孔介質結構的復雜性及微米尺度孔隙間的動態滲流過程難以捕捉，精確描述多孔介質內部滲流機理仍是一個富有挑戰的課題.

基于激光刻蝕技術制備而成的微觀孔隙模型，可以在準二維尺度上復刻多孔介質的微細觀孔隙結構，且采用顯微鏡和照相機直接觀測微米尺度下的滲流狀態，能夠采集到比三維透明土土柱實驗更高分辨率的流體運動圖像[9].因此，利用簡化的二維模型研究多孔介質的滲透特性也是很多學者的研究熱點，尤其是應用于地下油氣驅替開采[10]、二氧化碳封存[11]等方向.Akarbi等[12]使用微模型研究了流體流進矩形通道后的壓力變化，并與任意截面通道壓降理論解相比較，發現兩者能較好地符合.Joseph等[13]對原位的儲層巖心進行取樣掃描，提取巖心的孔隙結構構建二維裂隙通道模型，發現巖心的滲透率隨多孔介質的孔隙和喉道數量增加而明顯上升.Tamayol 等[14]構建1D、2D和3D特征單元的纖維多孔結構，并利用Navier-Stokes方程，求解較大雷諾數（0.01＜Re＜200）下模型的滲透率，結果表明模型滲透率是孔隙率和纖維排列的函數關系.現有的文獻結果表明，微流控模型的實驗壓降數據與Darcy模型和Hagen-Poiseuille方程預測結果具有良好的一致性.

數值模擬在微米尺度下的多孔介質滲流機理研究方面發揮重要作用.常用的數值模擬方法有計算流體力學方法（CFD）[15-16]、格子玻爾茲曼方法（lattice Boltzmann method, LBM）[17]、孔隙網絡模型方法（pore network model, PNM）[18]等.Tamayol等[15]利用計算流體方法對矩形微通道的滲流特性進行模擬，并將模擬結果與實驗數據進行對比，結果吻合良好.Yazdchi等[19]通過聯合數學解析和有限單元模擬，探究孔隙微結構對規則顆粒排列的微模型滲透性的影響.黃永平等[16]利用計算流體力學模型建立層流條件下多孔介質內部滲流的理論模型，并進行數值模擬，分析流體在多孔介質內部的流動特性.研究結果表明，多孔介質內部的孔隙結構對滲流特性有著重要的影響.相比于其他的數值模擬方法，計算流體力學方法的優點在于模型參數較少，無須標定，但是對計算資源的要求較高.

由于受到微流控模型尺寸、形狀與邊壁摩擦等因素影響，已有的研究表明傳統的宏觀滲流模型不能準確預測微流控模型的滲流特征.比如，Kosar等[20]采用低長寬比的微模型探究層流條件下的壓降和摩擦因數,發現宏觀條件下試驗模型不能預測微細觀尺度下的試驗結果.Yeom等[21]考慮多孔介質的內部結構對滲透率的影響，但是忽略了邊界效應.Gunda等[22]測量低雷諾數下帶有集成微柱的微通道兩端的壓降和流阻，通過Darcy方程計算出模型的滲透率，但是實驗結果和理論值存在一定偏差.

本研究基于微流控芯片加工技術，通過制作準二維微流控芯片模型來模擬多孔介質內部的骨架及孔隙結構，通過測量芯片兩端的壓降并進行修正，計算芯片模型的滲透率.采用直接求解Navier-Stokes方程模擬單相流滲流過程，將微柱顆粒排列方式對芯片模型的滲透性影響量化，并與實驗結果進行對比研究.結合芯片孔隙參數，系統性研究微流控模型滲透率和各向異性（微柱排列方式）、微柱直徑、孔隙率的關系，并探究Kozeny-Carman公式的適用性，最后提出適用于預測準二維微流控芯片模型滲透率的修正平板流滲流方程.

1 實驗裝置與實驗流程

1.1 顯微鏡-微觀模型實驗裝置

采用如圖1（a）所示的實驗裝置開展微流控芯片模型滲流實驗.該實驗裝置主要由流量控制子系統和微觀模型觀測子系統2大部分組成.

圖1 顯微鏡-微觀模型實驗裝置Fig.1 Microscopy-micromodel experimental system

1）流量控制子系統由高精度流量注射泵（PHD ULTRA 4400， Harvard Apparatus，US）和壓力傳感器(uPS01-T116，LabSmith，US)組成，注射泵的量程在 0.01～10.00 m L/min，壓力傳感器的測量為-100～150 kPa.如圖1(a)所示，注射泵以恒定流量向微觀模型入口注入液體，壓力傳感器實時監測芯片兩端壓力變化.

2） 微觀模型觀測子系統由微觀模型安裝平臺、倒置的光學顯微鏡 (Ti2，Nikon Corp，Japan)和內置高速相機（Zyla 4.0，Andor Technology，Belfast，UK）組 成.微 觀 模 型 長L=20mm ，寬W=10mm，深度h=50 μm，內部設計有不同直徑D和排列方式的微柱.本次實驗所用芯片尺寸如表1所示.芯片模型前后出口均為環形樹枝分叉狀，分叉通道的寬為800 μm，其余尺寸如圖1(c)所示.由于二維聚二甲基硅氧烷（PDMS）彈性模量較小，在孔隙壓力較大時將發生較大彈性變形，對實驗壓力的測量分析帶來較大干擾[23], 因此本實驗采用玻璃制作微流控芯片模型.

表1 實驗芯片模型類型及內部結構參數1)Tab.1 Geometrical dimensions of different microfluidic chip models considered in present work

1.2 實驗流程

實驗過程包括無氣水的制備、芯片模型排氣和滲流平衡3個階段.每種類型芯片試驗重復3次，取平均值作為對應流量下的穩定壓力值，具體步驟如下.

1）無氣水的制備.在實驗開始前，將去離子水放入真空抽氣機中抽氣30 min，排除水中的氣體，避免氣泡吸附在芯片內部，從而對試驗結果造成影響.

2）芯片模型排氣.將實驗裝置按照圖1所示搭建完畢后，先將壓力傳感器清零，然后以較高流量(1000 μL/min)通入無氣水，利用高流速水流驅除芯片內部殘余的氣泡.當在顯微鏡下觀測芯片內部無明顯氣泡時，可認為此時芯片內部已經完全飽和.

3）滲流平衡.在芯片達到飽和后，調節壓力傳感器，將流量調節至100 μL/min，然后持續通入無氣去離子水，當壓力傳感器顯示的壓力穩定后進行記錄，然后逐級加大流量至2000 μL/min，重復上述步驟，整個試驗過程在室溫下(約20 ℃)進行.

2 數值模擬

2.1 物理模型

基于實驗芯片模型，數值模擬的物理模型尺寸設置為如圖2(a)所示，即L=20mm ，W=10mm，h=50 μm.圓形為微柱，2個圓之間的區域為內部孔隙通道，上下邊界固定，液體流向水平向左.液相介質設為水，在20 ℃條件下，密度為998.2（ kg/m3）， 黏 度為0.001 k g/m·s，流 固 界 面設置為無滑移邊界條件.基于芯片模型尺寸，選取3D建模.應用結構化網格方式對圖2(a)中的芯片模型進行網格劃分，圖2(b)為選取的單元示意圖.

圖2 數值模擬物理模型示意圖Fig.2 Schematic diagram of numerical simulation physical model

2.2 控制方程

多孔介質中流體流動過程滿足連續性假設，可采用連續性方程和 Navier-Stokes方程描述.即連續性方程為

動量守恒方程為

式中：u為x方向流體速度，v為y方 向 流 體 速 度,p為流體壓力, ρ為流體密度，μ 為流體黏度.

在實際工程應用中，通常將多孔介質內部流動在雷諾數較小時視為層流.為了簡化數值計算，針對上述方程作如下假定：1)內部流態為穩定層流；2)多孔介質的固體表面及通道內壁面均為光滑無滑移界面；3)流體的物性參數為常數；4)忽略質量力.

2.3 邊界條件

入口邊界條件設置為

出口邊界條件設置為

由于設置的流固界面為光滑無滑移邊界，因此界面處的流速還滿足如下關系：

2.4 模擬方法

Ansys Fluent是一款得到國內外廣泛應用的計算流體力學的商業軟件，尤其在單相流滲流和多相流驅替方面.為了提高網格質量，網格類型主要采取以六面體為主、四面體為輔的網格單元.通過設置自動捕捉模式，可以很好地捕捉復雜結構幾何特性的幾何小面.為了確保最后計算結果的可靠性，還對網格進行獨立性檢測.在流速為3.3×10-4m/s 時，圖3(a)展示芯片模型Sq 0.60-1000在不同網格尺寸下的入口壓力(ΔP)-網格節點數(N)圖.隨著網格最大尺寸的減小，物理模型網格節點數增大，模型入口壓力值增大.當物理模型網格節點數大于1 .1×106時，模型的入口壓力值基本保持不變，此時網格的最大尺寸設置為0.02 mm.為了加快模擬計算過程，物理模型網格的最大尺寸設置為0.02 mm.

圖3 數值模擬物理模型的網格獨立性檢驗和正確性驗證Fig.3 Grid independence test and validation of numerical microfluidic chip model

為了得到三維多孔介質內部流體的流動特性，采用有限體積法對控制方程進行求解，同時調用SIMPLEC算法對模型進行壓力-速度耦合運算.模型計算過程采取雙精度，動量方程的離散采用二階迎風格式，松弛因子設為 0.5，2個迭代步之間的殘差設置為10-6.當殘差曲線達到10-6且出口處的平均流速保持不變時，可認為此次迭代計算達到收斂狀態.

圖3(b)展示了芯片模型Sq 0.60-1000的入口壓力實驗值與數值模擬值對比結果，2條曲線趨勢保持一致.盡管模擬值略小于實驗值，但最大偏差為17.6%，此誤差在允許范圍之內，因此該數值模擬方法可用于模擬準二維微流控芯片模型的滲流問題.

3 實驗數據分析

雷諾數的定義為

式中：u為流速，d為管徑，μ為流體黏度系數.

根據研究對象的不同，主要有基于孔隙尺寸、顆粒直徑和滲透率的3種雷諾數表達式[24].在準二維微流控芯片模型中，參考文獻[22]，d取微柱直徑D.本次實驗條件下的雷諾數為0.1～210.0.

3.1 壓降分析

由于壓力傳感器測量的壓力值為芯片模型入口和出口處的壓力差值，并非芯片模型內部兩端，需要對所測壓力值進行修正[12]：

式中：ΔP為芯片模型內部壓降，ΔPm為實測壓降，ΔPt為連接入口和傳感器的圓形導管產生的壓降,ΔPc為芯片模型內部出口和入口處樹形分叉通道產生的壓降, ΔPin和 ΔPex為芯片模型入口和出口界面變化產生的壓降，ΔPd為不充分發展流部分產生的壓降, ΔPb為液體直角轉彎產生的壓降,ΔPev為電黏滯效應產生的壓降.

對于 ΔPc，Bahrami等[25]基于矩形斷面的Poiseuille流，結合達西公式提出相應的解析解：

式中：ff為擴散摩擦因子；L為通道長度；Dh為矩形截面的特征長度；ξ為通道矩形截面的縱橫比，并且 ξ ＜1.0 ；ReA1/2為雷諾數.Yovanovich[26]通過研究發現，對于非圓形截面通道，A1/2作為特征長度用于計算雷諾數更為合適，A為矩形通道截面面積，在本實驗中，Dh=A1/2.

由圖1(c)中芯片入口和出口處的樹形分叉尺寸，并 聯立 式(9)～(11)可求得ΔPc.對于ΔPt，可直接根據圓形管道中的阻力計算：

式中：L1為 圓形導管的長度；d1為導管內直徑，d1=0.8 mm.

對于直角轉彎處和管道中截面變化產生的壓降 ΔPb、 ΔPin和 ΔPex， 文獻[12]中給相應的計算公式：

式中：Q 為流量；At為圓形管道截面積；Kc、Ke和Kb均為損失系數，均取1.2.Akbari 等[12]經過計算，發現在低雷諾數下，直角轉彎處和管道中截面變化產生的壓降之和占模型測量壓降值的百分比小于1%；在較高雷諾數下（Re≥100）,直角轉彎處和管道中截面變化產生的壓降之和占比小于5%.實驗的研究范圍主要集中在低雷諾數流速區間，為了簡化計算，忽略該值的影響.

對 于ΔPev和ΔPd，由于芯片為玻璃制作，變 更且內部流態為層流，在不同流速和不同孔隙率模型下，與測量壓降值的占比小于0.3%[12]，因此可以忽略不計.結合芯片模型尺寸及式(9)～(13)求解式（8），得到芯片模型內部主體的壓降修正值 ΔP.圖4給出在孔隙率均為0.60時，微柱直徑分別為500、1000和2000μm的3種芯片模型在微柱方形排列方式下的實驗單位壓降隨流速的變化規律.在整個實驗流速階段（0.00167～0.06700m/s）,3種芯片模型的單位壓降與流速的線性擬合度分別為0.9994、0.9982和0.9969，均具有良好的線性關系.由Darcy滲透公式可知，在該流速范圍內，芯片模型中的流態為層流.在相同排列方式下，微柱直徑小的芯片模型在相同流量下表現出更大的壓降.

圖4 當孔隙率為 0 .60時不同微柱直徑芯片模型的實驗單位壓降隨流速的變化規律Fig.4 Variation of unit experimental pressure drop in microfluidic chip models with different velocities under different micropillar diameters when porosity was 0.60

3.2 流動阻力分析

流動阻力S又稱摩擦阻力，簡稱流阻，是黏性流體在運動過程中產生的阻礙流動的反作用力，定義為

圖5所示在孔隙率 ε =0.60，不同微柱直徑和排列方式下流阻隨流速變化規律.當流速大于0.005m/s時，流阻S隨著流速u的增加基本保持不變；在相同流速和微柱排列方式條件下，微柱直徑越小，芯片模型的流阻更大；在相同流量和微柱直徑條件下，錯開排列芯片模型比方形排列芯片模型具有更大的流阻.然而在流速較小時（u＜0.005m/s），4種芯片模型的結果均表現為先減小后增大，并且初始最大值大于后期穩定值.

圖5 當孔隙率 為 0.60時不同微柱直徑和排列方式下實驗流阻隨流量變化規律Fig.5 Variation of experimental flow resistance in microfluidic chips with different flow rate under different micro pillar diameter and arrangements when porosity was 0.60

Kundu等[24]對顆粒直徑分別為 2 .5、3.5和5mm的土柱進行滲流實驗，發現在流速處于 2.88×10-5～1.90×10-4m/s時，流阻隨著流速的增大而減小.當流速持續增大時，流阻隨著流速的增加保持不變，即由前達西流流態進入達西流流態.在本實驗中，當流速u＜0.0035m/s時，前達西流流態比較明顯.在低流速時，孔隙通道中由電滲力和毛細力產生的流動勢能在與流體方向相反的壁面上產生小的逆流[27].當流速處于 0 .0035～0.0050m/s時，注射泵動力系統的脈沖效應比較明顯，這對實驗結果造成一定的影響.

4 滲透性分析

在低雷諾數下，孔隙中的流體主要受黏滯力控制，水力梯度與流速呈線性關系.可由Darcy滲透公式表示：

式中：K為滲透率.

4.1 Kozeny-Carman公式預測芯片模型滲透率

對于多孔介質材料滲透特性的預測最早由Kozeny開展，通過一系列等長度和直徑的平行毛細管模型計算流體流動時多孔介質兩端的壓降，進而分析多孔介質的滲透率.Carman對該方法進行修正，得到較合理的Kozeny-Carman（KC）公式.盡管KC公式是一個半經驗公式，但是由于含有的參數簡單并且容易獲取，因此被廣泛應用于地下水滲流和油氣開采等工程領域，表達式為

式中：c為Kozeny常數，S0為比表面積.隨后引入迂曲度τ進行修正：

式中：k=cτ2為KC常數，影響大小的因素主要有孔隙率、孔隙形狀和大小、迂曲度等[28-29].

在大量實驗基礎上，Kaviany[30]提出當迂曲度τ 為時，球形顆粒組成的多孔介質中的Kozeny常數c值大致為2.5，因此KC常數為5.0.對于泥炭層介質，Mathava等[31]發現KC常數為3.4.

根據實驗所測壓降值，結合達西滲透式(15)可求得芯片模型內部通道的滲透率.圖6(a)給出不同孔隙率和排列方式下的芯片模型實驗結果對比圖.從圖中可以看出，芯片模型滲透率隨著微柱直徑的增大而增大，但是增大最大差值在1倍之內.在相同孔隙率和微柱直徑條件下，錯開排列芯片模型具有更小的滲透率，這與前面錯開排列模型中具有較大流阻的結果一致；在相同排列方式和微柱直徑條件下，孔隙率越小，芯片模型滲透率更小，并且滲透率差值隨著微柱直徑的增大而增大.

為了探究KC公式對微流控芯片模型滲透率預測的適用性，圖6(b)給出在孔隙率為0.60時，不同芯片模型滲透率的實驗值、模擬值和KC預測值的對比.對于KC預測值，由于微流控芯片模型在尺寸比例上更接近煤炭層介質結構，故取KC常數為3.4.圖中數據表明，模擬值和KC預測值均大于實驗值.模擬值比實驗值的誤差范圍為9.78%～28.43%，KC預測值比實驗值的誤差范圍為9.67%～73.97%，并且誤差隨著微柱直徑的增大而減小.造成模擬值與實驗值誤差的原因主要有：1）在芯片模型主體壓降值的計算過程中，采取部分簡化處理；2）在數值模擬運算過程中，做部分假設；3）加工工藝無法保證芯片模型內壁面完全光滑；4）流量泵的脈沖影響.KC公式是在三維實驗的基礎上提出的半經驗公式，KC常數也是一個受孔隙率和孔隙結構影響的系數，并且隨孔隙率的增大而增大[32]，并非固定不變，因此對于準二維微流控芯片模型的滲透率預測，KC公式存在較大誤差.

圖6 不同類型芯片模型的滲透率對比Fig.6 Comparison of permeability obtained by different microfluidic chip models

芯片模型孔隙率越大，內部流體流態越接近平板流.平板流的滲透率表達式為[33]

式中：a為兩平行板的開度.

圖7展示孔隙率從0.327增大至0.900，芯片模型滲透率實驗值與數值模擬值、KC預測值及平板流滲流理論值的對比.由圖中曲線可知，KC公式僅在孔隙率為0.50～0.60時，大致適用于準二維微流控芯片模型的滲透率預測；當孔隙率小于0.50時，KC公式預測值小于芯片模型模擬值；當孔隙率大于0.60時，KC公式預測值大于芯片模型模擬值，并且差值趨勢增大.當孔隙率逐漸增大（接近1）時，芯片模型模擬值收斂于平板流滲流理論值.基于此現象，提出適用于預測準二維微流控芯片模型滲透率的修正平板流滲流公式：

式中：h為芯片模型的微柱高度.由圖7中數據可知，修正平板流滲流公式可以準確預測準二維微流控芯片模型滲流值.

圖7 芯片模型預測滲透率與實驗值、模擬值的對比Fig.7 Comparison of predicted permeability of microfluidic chip models with experimental data and numerical data

4.2 迂曲度

迂曲度是影響多孔介質滲透率的重要參數，也是定量描述顆粒排列方式對多孔介質滲透性影響的最重要指標.迂曲度越大，表明流體在孔隙中的流動路徑越曲折，進而造成更大的流阻.Carman[6]將迂曲度定義為有效流線長度與模型長度的比值.基于上述理論，采取數值模擬方法求解流線長度，取h=0.025mm處的剖面作為特征平面，得到每個芯片模型的平均迂曲度值，對多孔介質的滲透性影響進行量化：

式中：l1為一條流線長度，如圖8所示.

圖8 芯片模型流線示意圖Fig.8 Schematic diagram of streamline for chip models

根據上述迂曲度的定義，在流速為0.0167 m/s時，利用數值模擬方法求得的芯片模型平均迂曲度值如圖9所示.不同流速下流線形態發生變化，平均迂曲度值將發生細微變化，但波動幅值在5%以內.由圖中曲線可知，當孔隙率固定時，平均迂曲度的值隨微柱直徑增大而減小，且下降趨勢減弱.例如，當孔隙率為0.60時，方形排列的芯片模型在微柱直徑分別是 500、1000和 2000 μm時，平均迂曲度值分別為1.1670、1.1565和1.1501.在微柱排列方式相同下，方形排列的芯片模型的平均迂曲度值為1.05左右，明顯小于錯開排列的芯片模型平均迂曲度值1.13.在相同排列方式下，孔隙率越小，平均迂曲度越大.相比于排列方式，孔隙率對平均迂曲度的影響更大.這從微觀視角解釋了圖6(a)的實驗結果.上述模擬結果表明，迂曲度并非常數21/2，而是一個和孔隙率、顆粒排列方式、顆粒直徑等因素相關的參數.上述結果也顯示，在相同孔隙率條件下，顆粒直徑的影響小于2%.因此在實際工程應用中，顆粒直徑對迂曲度的影響可以忽略，影響迂曲度的因素主要為孔隙率、顆粒排列、孔隙大小和分布[34].

圖9 不同類型芯片模型平均迂曲度對比Fig.9 Comparison of average tortuosity in different microfluidic chip models

4.3 各向異性對滲透率的影響

在相同孔隙率下，芯片模型微柱排列方式的各向異性將對芯片模型滲透率產生影響.如圖10所示，在孔隙率為0.60時，旋轉角 θ分別是20°、30°、40°、50°和60°這5種微柱錯開排列模式下的芯片模型數值模擬流速場云圖.隨著旋轉角的增大，微柱孔喉間的最大流速逐漸較小；在流動方向上，微柱顆粒的前后兩端出現三角形狀的“低速區”；在垂直流動方向上出現“高速區”、“低速區”面積遠小于“高速區”面積； “高速區”隨著旋轉角的增大而減小，形狀也由三角形變成“香蕉形”.當孔隙中發生細顆粒遷移時，“低速區”極易形成顆粒堆積[35].

圖10 當流量為10 μL/min、孔隙率為0.60時，旋轉角 θ對芯片模型(D=1000 μm)速度分布的影響Fig.10 Effect of rotation angle θ on velocity distribution of microfluidic chip models obtained by numerical simulation when flow rate was 10 μL/min, micropillar diameter was 1000 μm and porosity was 0.60

圖11(a)展示利用CFD建模計算，在微柱錯開排列模式下，各向異性（旋轉角 θ）對芯片模型滲透率的影響.由圖中曲線可以得到，在相同孔隙率下，旋轉角 θ對芯片模型滲透率的有一定影響.當 θ =20°時，微柱排列對流體的流阻最大，此時芯片模型的滲透率最小；反之當 θ =70°時，芯片模型的滲透性最大；在中間區域，旋轉角 θ對芯片模型的滲透率幾乎沒有影響.以孔隙率為0.60的模型為 例，當θ =20°時，芯片模型的滲透率為5.69×10-11m2；當θ =70°時，芯片模型的滲透率為9.27×10-11m2；當θ 在30°～60°時，芯片模型的滲透率為 8 .0×10-11m2左右，該值一般作為特征值為采用.3種孔隙率下旋轉角 θ產生的滲透率相對于特征值的幅值分別為 28.87%(ε=0.60)、15.39%(ε=0.70)和7.2%(ε=0.80).孔隙率越小，幅值越大，即旋轉角在較大和較小時對芯片模型的滲透率影響越大，這也為在實際工程中，顆粒填充提供新的視角.

相對于顆粒排列，孔喉直徑也常被用于研究對多孔介質滲透率的影響.在旋轉角 θ變化的過程中，芯片模型微柱間的孔喉直徑也相應發生變化，隨著 θ的增大而增大.圖11(b)展示芯片模型滲透率隨孔喉直徑的變化規律，擬合曲線為對數關系.3種孔隙率下的曲線數據均表明，芯片模型滲透率隨著孔喉直徑的增大而增大；擬合曲線的系數為正，表明滲透率隨著孔喉增大而增大的趨勢減小.在不考慮優勢流通道的前提下，相對于大孔喉、小孔喉對多孔介質的滲透性更具有決定性作用.

圖11 旋轉角和孔喉直徑對微柱錯開排列芯片模型滲透率的影響Fig.11 Effect of rotation angle and pore throat diameter on permeability of microfluidic chip models obtained by numerical simulation under staggered arrangement

為了進一步探究流體在芯片模型中的流動特性，圖12(a)展示芯片模型平均迂曲度隨旋轉角 θ的變化規律.曲線表明迂曲度在 θ =20°時最大，在θ=70°時最小.在旋轉角度相同情況下，孔隙率越小，迂曲度越大.圖12(b)所示為芯片模型平均迂曲度值與相關文獻中的數據對比分析.Comiti等[36]提出的模型與旋轉角較小時芯片模型的平均迂曲度值吻合較好.當旋轉角較大時，芯片模型平均迂曲度值落在最下方，表明3種模型都不能單獨預測該排列情況下的平均迂曲度值；當旋轉角在30°～60°時，芯片模型平均迂曲度值均落在文獻[37-38]中，具有較強的一致性.每個模型的假設條件不同，參考文獻中提出的迂曲度模型的適用對象并不具有普適性.后期的研究方向可以聚焦于 無序排列和三維顆粒堆積方面的實驗和模擬.

圖12 不同旋轉角下芯片模型的平均迂曲度值Fig.12 Average tortuosity of microfluidic chip models under different rotation angles

5 結論

針對均勻顆粒周期性排列的準二維微流控芯片模型的滲透性開展一些列實驗和數值模擬研究，得到以下結論.

1）在相同孔隙率下，相比于微柱方形排列的芯片模型，微柱錯開排列的芯片模型的滲流路徑更長，微觀上表現在迂曲度值更大，宏觀上表現為更大的流阻和壓降，最終體現為滲透率更小.

2）KC公式不適用于預測準二維微流控芯片模型的滲透率，而提出的修正平板流滲流公式能準確預測準二維微流控芯片模型滲透率.

3）實驗數據和數值模擬結果表明，微柱直徑和排列方式（各向異性）也是影響芯片模型滲透率大小的因素.

4）數值模擬結果表明，在各向異性下，旋轉角在 20°和70°時，對均勻顆粒等間距排列的芯片模型的滲透率影響最大，而在中間角度（30°～60°）下芯片模型的滲透率基本一致.數值模擬結算得到的芯片模型平均迂曲度值也與文獻中的參考數據吻合較好.