多孔介质渗流现象广泛存在于自然界及工程应用中，如土壤中的地下水渗流、土石坝工程的稳定性、石油天然气的开采等^[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 < {R}{e} < 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（a）所示的实验装置开展微流控芯片模型渗流实验. 该实验装置主要由流量控制子系统和微观模型观测子系统2大部分组成.

1）流量控制子系统由高精度流量注射泵（PHD ULTRA 4400， Harvard Apparatus，US）和压力传感器(uPS01-T116，LabSmith，US)组成，注射泵的量程在 0.01~10.00 $ \mathrm{m}\mathrm{L}/\mathrm{m}\mathrm{i}\mathrm{n} $，压力传感器的测量为−100~150 kPa. 如图1(a)所示，注射泵以恒定流量向微观模型入口注入液体，压力传感器实时监测芯片两端压力变化.

2） 微观模型观测子系统由微观模型安装平台、倒置的光学显微镜 (Ti2，Nikon Corp，Japan)和内置高速相机（Zyla 4.0，Andor Technology，Belfast，UK）组成. 微观模型长 $ L=20\;\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{m}\mathrm{m} $，宽 $ W=10\;\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{m}\mathrm{m} $，深度 $ h=50 $ µm，内部设计有不同直径D和排列方式的微柱. 本次实验所用芯片尺寸如表1所示. 芯片模型前后出口均为环形树枝分叉状，分叉通道的宽为800 μm，其余尺寸如图1(c)所示. 由于二维聚二甲基硅氧烷（PDMS）弹性模量较小，在孔隙压力较大时将发生较大弹性变形，对实验压力的测量分析带来较大干扰^[23], 因此本实验采用玻璃制作微流控芯片模型.

实验过程包括无气水的制备、芯片模型排气和渗流平衡3个阶段. 每种类型芯片试验重复3次，取平均值作为对应流量下的稳定压力值，具体步骤如下.

1）无气水的制备. 在实验开始前，将去离子水放入真空抽气机中抽气30 min，排除水中的气体，避免气泡吸附在芯片内部，从而对试验结果造成影响.

2）芯片模型排气. 将实验装置按照图1所示搭建完毕后，先将压力传感器清零，然后以较高流量(1 000 µL/min)通入无气水，利用高流速水流驱除芯片内部残余的气泡. 当在显微镜下观测芯片内部无明显气泡时，可认为此时芯片内部已经完全饱和.

3）渗流平衡. 在芯片达到饱和后，调节压力传感器，将流量调节至100 µL/min，然后持续通入无气去离子水，当压力传感器显示的压力稳定后进行记录，然后逐级加大流量至2 000 µL/min，重复上述步骤，整个试验过程在室温下(约20 ℃)进行.

基于实验芯片模型，数值模拟的物理模型尺寸设置为如图2(a)所示，即 $ L=20\;\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{m}\mathrm{m} $， $ W= 10\;\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{m}\mathrm{m} $， $ h=50$ µm. 圆形为微柱，2个圆之间的区域为内部孔隙通道，上下边界固定，液体流向水平向左. 液相介质设为水，在20 ℃条件下，密度为998.2 $（\mathrm{k}\mathrm{g}/{{\rm{m}}}^{3}）$，黏度为0.001 $ \mathrm{k}\mathrm{g}/\mathrm{m}\cdot \mathrm{s} $，流固界面设置为无滑移边界条件. 基于芯片模型尺寸，选取3D建模. 应用结构化网格方式对图2(a)中的芯片模型进行网格划分，图2(b)为选取的单元示意图.

多孔介质中流体流动过程满足连续性假设，可采用连续性方程和 Navier-Stokes方程描述.即连续性方程为

(1) $ \frac{{\partial u}}{{\partial x}}+\frac{{\partial v}}{{\partial y}} = 0. $

动量守恒方程为

(2) $ u\frac{{\partial u}}{{\partial x}}+v\frac{{\partial u}}{{\partial y}} = - \frac{1}{\rho }\frac{{\partial p}}{{\partial x}}+\mu \left(\frac{{{\partial ^2}u}}{{\partial {x^2}}}+\frac{{{\partial ^2}u}}{{\partial {y^2}}}\right), $

(3) $ u\frac{{\partial v}}{{\partial x}}+v\frac{{\partial v}}{{\partial y}} = - \frac{1}{\rho }\frac{{\partial p}}{{\partial y}}+\mu \left(\frac{{{\partial ^2}v}}{{\partial {x^2}}}+\frac{{{\partial ^2}v}}{{\partial {y^2}}}\right). $

式中： $ u $为 $ x $方向流体速度， $ v $为 $ y $方向流体速度, p为流体压力, $ \;\rho $为流体密度， $ \;\mu $为流体黏度.

在实际工程应用中，通常将多孔介质内部流动在雷诺数较小时视为层流. 为了简化数值计算，针对上述方程作如下假定：1)内部流态为稳定层流；2)多孔介质的固体表面及通道内壁面均为光滑无滑移界面；3)流体的物性参数为常数；4)忽略质量力.

入口边界条件设置为

(4) $ \left. \begin{gathered} x = 0, \\ u = {u_{{\text{in}}}}, \\ v = 0. \\ \end{gathered} \right\} $

出口边界条件设置为

(5) $ \left. \begin{gathered} x = L, \\ \frac{{\partial u}}{{\partial x}} = \frac{{\partial v}}{{\partial x}} = 0. \\ \end{gathered} \right\} $

由于设置的流固界面为光滑无滑移边界，因此界面处的流速还满足如下关系：

(6) $ u = v = 0. $

Ansys Fluent是一款得到国内外广泛应用的计算流体力学的商业软件，尤其在单相流渗流和多相流驱替方面.为了提高网格质量，网格类型主要采取以六面体为主、四面体为辅的网格单元. 通过设置自动捕捉模式，可以很好地捕捉复杂结构几何特性的几何小面. 为了确保最后计算结果的可靠性，还对网格进行独立性检测. 在流速为 $ 3.3\times {10}^{-4}\;\mathrm{m}/\mathrm{s} $ 时，图3(a)展示芯片模型Sq 0.60-1000在不同网格尺寸下的入口压力( $\Delta {P}$)-网格节点数( $ {N} $)图. 随着网格最大尺寸的减小，物理模型网格节点数增大，模型入口压力值增大. 当物理模型网格节点数大于 $ 1.1\times {10}^{6} $时，模型的入口压力值基本保持不变，此时网格的最大尺寸设置为0.02 mm. 为了加快模拟计算过程，物理模型网格的最大尺寸设置为0.02 mm.

为了得到三维多孔介质内部流体的流动特性，采用有限体积法对控制方程进行求解，同时调用SIMPLEC算法对模型进行压力-速度耦合运算. 模型计算过程采取双精度，动量方程的离散采用二阶迎风格式，松弛因子设为 0.5，2个迭代步之间的残差设置为 $ {10}^{-6} $. 当残差曲线达到 $ {10}^{-6} $且出口处的平均流速保持不变时，可认为此次迭代计算达到收敛状态.

图3(b)展示了芯片模型Sq 0.60-1 000的入口压力实验值与数值模拟值对比结果，2条曲线趋势保持一致. 尽管模拟值略小于实验值，但最大偏差为17.6%，此误差在允许范围之内，因此该数值模拟方法可用于模拟准二维微流控芯片模型的渗流问题.

雷诺数的定义为

(7) $ {Re} = \frac{{\rho ud}}{\mu }. $

式中：u为流速，d为管径， $ \mu $为流体黏度系数.

根据研究对象的不同，主要有基于孔隙尺寸、颗粒直径和渗透率的3种雷诺数表达式^[24]. 在准二维微流控芯片模型中，参考文献[22]， $ d $ 取微柱直径D. 本次实验条件下的雷诺数为0.1~210.0.

由于压力传感器测量的压力值为芯片模型入口和出口处的压力差值，并非芯片模型内部两端，需要对所测压力值进行修正^[12]：

(8) $ \Delta P = \Delta {P_{\rm{m}}} - (\Delta {P_{\rm{t}}}+\Delta {P_{\rm{c}}}+\Delta {P_{\rm{{in}}}}+\Delta {P_{\rm{d}}}+\Delta {P_{\rm{{ex}}}}+2\Delta {P_{\rm{b}}}+\Delta {P_{\rm{{ev}}}}). $

式中： $ \Delta P $为芯片模型内部压降， $ \Delta {P}_{\mathrm{m}} $为实测压降， $ \Delta {P}_{\mathrm{t}} $为连接入口和传感器的圆形导管产生的压降, $ \Delta {P}_{\mathrm{c}} $为芯片模型内部出口和入口处树形分叉通道产生的压降, $ \Delta {P}_{\mathrm{i}\mathrm{n}} $和 $ \Delta {P}_{\mathrm{e}\mathrm{x}} $为芯片模型入口和出口界面变化产生的压降， $\Delta {P}_{\mathrm{d}}$为不充分发展流部分产生的压降, $ \Delta {P}_{\mathrm{b}} $为液体直角转弯产生的压降, $ \Delta {P}_{\mathrm{e}\mathrm{v}} $为电黏滞效应产生的压降.

对于 $ \Delta {P}_{\mathrm{c}} $，Bahrami等^[25]基于矩形断面的Poiseuille流，结合达西公式提出相应的解析解：

(9) $ \Delta {P_{\rm{c}}} = \frac{{2{f_f}\rho {u^2}L}}{{{D_{\rm{h}}}}}, $

(10) $ {f_f}{{\rm{Re}}} = \frac{{12}}{{\left( {1 - \dfrac{{192}}{{{{\text{π}} ^5}}}\xi \tanh \left( {{\text{π}} /2\xi } \right)} \right)\left( {1+\xi } \right){\xi ^{1/2}}}}, $

(11) $ {{{Re}}}A^ {{1/2}} = \frac{{\rho u{A^{1/2}}}}{\mu }. $

式中： $ {f}_{f} $为扩散摩擦因子；L为通道长度； $ {D}_{\mathrm{h}} $为矩形截面的特征长度； $ \xi $为通道矩形截面的纵横比，并且 $\xi < 1.0$； ${{{Re}}}A^ {{1/2}}$为雷诺数. Yovanovich^[26]通过研究发现，对于非圆形截面通道， $ {A^{1/2}} $作为特征长度用于计算雷诺数更为合适，A为矩形通道截面面积，在本实验中， $D_{\rm{h}}= {A^{1/2}} $.

由图1(c)中芯片入口和出口处的树形分叉尺寸，并联立式(9)~(11)可求得 $ \Delta {P}_{\mathrm{c}} $. 对于 $ \Delta {P}_{\mathrm{t}} $，可直接根据圆形管道中的阻力计算：

(12) $ \Delta {P_t} = \frac{{32\mu {L_1}u}}{{d_1^2}}. $

式中： $ {L}_{1} $为圆形导管的长度； $ {d}_{1} $为导管内直径， $d_1 $=0.8 mm.

对于直角转弯处和管道中截面变化产生的压降 $ \Delta {P}_{\mathrm{b}} $、 $ \Delta {P}_{\mathrm{i}\mathrm{n}} $和 $ \Delta {P}_{\mathrm{e}\mathrm{x}} $，文献[12]中给相应的计算公式：

(13) $ \Delta {P_{\rm{{in}}}}+\Delta {P_{\rm{{ex}}}}+2\Delta {P_{\rm{b}}} = \frac{{\rho {Q^2}}}{{2A}}\left[ {{K_{\rm{c}}}+{K_{\rm{e}}}+2{K_{\rm{b}}} \times {{\left( {A/{A_{\rm{t}}}} \right)}^2}} \right]. $

式中： $ \mathrm{Q} $为流量； $ {A}_{\mathrm{t}} $为圆形管道截面积； $ {K}_{c}\mathrm{、}{K}_{e} \mathrm{和}{K}_{b} $均为损失系数，均取1.2. Akbari 等^[12]经过计算，发现在低雷诺数下，直角转弯处和管道中截面变化产生的压降之和占模型测量压降值的百分比小于1%；在较高雷诺数下（ ${R}{e}\ge 100$）,直角转弯处和管道中截面变化产生的压降之和占比小于5%. 实验的研究范围主要集中在低雷诺数流速区间，为了简化计算，忽略该值的影响.

对于 $ \Delta {P}_{\mathrm{e}\mathrm{v}} $和 $\Delta {P}_{\mathrm{d}}$，由于芯片为玻璃制作，变更且内部流态为层流，在不同流速和不同孔隙率模型下，与测量压降值的占比小于0.3%^[12]，因此可以忽略不计. 结合芯片模型尺寸及式(9)~(13)求解式（8），得到芯片模型内部主体的压降修正值 $ \Delta P $. 图4给出在孔隙率均为0.60时，微柱直径分别为500 、1 000 和2 000 µm的3种芯片模型在微柱方形排列方式下的实验单位压降随流速的变化规律. 在整个实验流速阶段（ $0.001 \; \mathrm{ }67\mathrm{ }\sim 0.067\;\mathrm{ }00\;\mathrm{ }\mathrm{m}/\mathrm{s}$）, 3种芯片模型的单位压降与流速的线性拟合度分别为0.999 4、0.998 2和0.996 9，均具有良好的线性关系. 由Darcy渗透公式可知，在该流速范围内，芯片模型中的流态为层流. 在相同排列方式下，微柱直径小的芯片模型在相同流量下表现出更大的压降.

流动阻力S又称摩擦阻力，简称流阻，是黏性流体在运动过程中产生的阻碍流动的反作用力，定义为

(14) $ S = \frac{{\Delta P}}{Q}. $

图5所示在孔隙率 $ \varepsilon =0.60 $，不同微柱直径和排列方式下流阻随流速变化规律. 当流速大于 $ 0.005\;\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{m}/\mathrm{s} $时，流阻S随着流速u的增加基本保持不变；在相同流速和微柱排列方式条件下，微柱直径越小，芯片模型的流阻更大；在相同流量和微柱直径条件下，错开排列芯片模型比方形排列芯片模型具有更大的流阻. 然而在流速较小时（ $ u < 0.005\;\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{m}/\mathrm{s} $），4种芯片模型的结果均表现为先减小后增大，并且初始最大值大于后期稳定值. Kundu等^[24]对颗粒直径分别为 $2.5\;\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{}\mathrm{}\mathrm{、} 3.5\;\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{}\mathrm{}\mathrm{和} 5\;\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{m}\mathrm{m}$的土柱进行渗流实验，发现在流速处于 $2.88\times {10}^{-5} \sim 1.90\times {10}^{-4}\;\mathrm{m}/\mathrm{s}$时，流阻随着流速的增大而减小. 当流速持续增大时，流阻随着流速的增加保持不变，即由前达西流流态进入达西流流态. 在本实验中，当流速 $ u < 0.003\;\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }5\;\mathrm{ }\mathrm{m}/\mathrm{s} $时，前达西流流态比较明显. 在低流速时，孔隙通道中由电渗力和毛细力产生的流动势能在与流体方向相反的壁面上产生小的逆流^[27]. 当流速处于 $ 0.003\;\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }5\sim 0.005\;\mathrm{ }0\mathrm{ }\mathrm{m}/\mathrm{s} $时，注射泵动力系统的脉冲效应比较明显，这对实验结果造成一定的影响.

在低雷诺数下，孔隙中的流体主要受黏滞力控制，水力梯度与流速呈线性关系. 可由Darcy渗透公式表示：

(15) $ - \frac{{{\rm{d}}p}}{{{\rm{d}}L}} = \frac{\mu }{K}u. $

式中：K为渗透率.

对于多孔介质材料渗透特性的预测最早由Kozeny开展，通过一系列等长度和直径的平行毛细管模型计算流体流动时多孔介质两端的压降，进而分析多孔介质的渗透率. Carman对该方法进行修正，得到较合理的Kozeny-Carman（KC）公式. 尽管KC公式是一个半经验公式，但是由于含有的参数简单并且容易获取，因此被广泛应用于地下水渗流和油气开采等工程领域，表达式为

(16) $ K = \frac{1}{{c{{S_0}^2}}}\frac{{{\varepsilon ^3}}}{{{{\left( {1 - \varepsilon } \right)}^2}}}. $

式中： $ c $为Kozeny常数， ${{S}}_0$为比表面积. 随后引入迂曲度 $ \tau $进行修正：

(17) $ K = \frac{1}{{c{\tau ^2}{{S_0}^2}}}\frac{{{\varepsilon ^3}}}{{{{\left( {1 - \varepsilon } \right)}^2}}} = \frac{1}{{k{{S_0}^2}}}\frac{{{\varepsilon ^3}}}{{{{\left( {1 - \varepsilon } \right)}^2}}}. $

式中： $ k=c{\tau }^{2} $为KC常数，影响大小的因素主要有孔隙率、孔隙形状和大小、迂曲度等^[28-29].

在大量实验基础上，Kaviany^[30]提出当迂曲度 $ \tau $为 $\sqrt 2$时，球形颗粒组成的多孔介质中的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公式存在较大误差.

芯片模型孔隙率越大，内部流体流态越接近平板流. 平板流的渗透率表达式为^[33]

(18) $ K = \frac{{{a^2}}}{{12}}. $

式中：a为两平行板的开度.

图7展示孔隙率从0.327增大至0.900，芯片模型渗透率实验值与数值模拟值、KC预测值及平板流渗流理论值的对比. 由图中曲线可知，KC公式仅在孔隙率为0.50~0.60时，大致适用于准二维微流控芯片模型的渗透率预测；当孔隙率小于0.50时，KC公式预测值小于芯片模型模拟值；当孔隙率大于0.60时，KC公式预测值大于芯片模型模拟值，并且差值趋势增大. 当孔隙率逐渐增大（接近1）时，芯片模型模拟值收敛于平板流渗流理论值. 基于此现象，提出适用于预测准二维微流控芯片模型渗透率的修正平板流渗流公式：

(19) $ K = \frac{{{h^2}}}{{12}}{\varepsilon ^2}. $

式中： $ h $为芯片模型的微柱高度. 由图7中数据可知，修正平板流渗流公式可以准确预测准二维微流控芯片模型渗流值.

迂曲度是影响多孔介质渗透率的重要参数，也是定量描述颗粒排列方式对多孔介质渗透性影响的最重要指标. 迂曲度越大，表明流体在孔隙中的流动路径越曲折，进而造成更大的流阻. Carman^[6]将迂曲度定义为有效流线长度与模型长度的比值. 基于上述理论，采取数值模拟方法求解流线长度，取 $ h=0.025\;{\rm{mm}} $处的剖面作为特征平面，得到每个芯片模型的平均迂曲度值，对多孔介质的渗透性影响进行量化：

(20) $ {l_1} = \sum\nolimits_{i} {} {\sqrt {{{\left( {{x_i} - {x_{i - 1}}} \right)}^2}+{{\left( {{y_i} - {y_{i - 1}}} \right)}^2}} } , $

(21) $ \overline \tau =\sum_1^i \frac{{l_i}}{{iL}} .$

式中： $ {l}_{1} $为一条流线长度，如图8所示.

根据上述迂曲度的定义，在流速为0.016 7 m/s时，利用数值模拟方法求得的芯片模型平均迂曲度值如图9所示. 不同流速下流线形态发生变化，平均迂曲度值将发生细微变化，但波动幅值在5%以内. 由图中曲线可知，当孔隙率固定时，平均迂曲度的值随微柱直径增大而减小，且下降趋势减弱. 例如，当孔隙率为0.60时，方形排列的芯片模型在微柱直径分别是 $ 500$、 $ 1\;000 $和 $ 2\;000$ µm时，平均迂曲度值分别为1.1670、1.156 5和1.150 1. 在微柱排列方式相同下，方形排列的芯片模型的平均迂曲度值为1.05左右，明显小于错开排列的芯片模型平均迂曲度值1.13. 在相同排列方式下，孔隙率越小，平均迂曲度越大. 相比于排列方式，孔隙率对平均迂曲度的影响更大. 这从微观视角解释了图6(a)的实验结果. 上述模拟结果表明，迂曲度并非常数2^1/2，而是一个和孔隙率、颗粒排列方式、颗粒直径等因素相关的参数. 上述结果也显示，在相同孔隙率条件下，颗粒直径的影响小于2%. 因此在实际工程应用中，颗粒直径对迂曲度的影响可以忽略，影响迂曲度的因素主要为孔隙率、颗粒排列、孔隙大小和分布^[34].

在相同孔隙率下，芯片模型微柱排列方式的各向异性将对芯片模型渗透率产生影响. 如图10所示，在孔隙率为0.60时，旋转角 $ \theta $分别是20°、30°、40°、50°和60°这5种微柱错开排列模式下的芯片模型数值模拟流速场云图. 随着旋转角的增大，微柱孔喉间的最大流速逐渐较小；在流动方向上，微柱颗粒的前后两端出现三角形状的“低速区”；在垂直流动方向上出现“高速区”、“低速区”面积远小于“高速区”面积； “高速区”随着旋转角的增大而减小，形状也由三角形变成“香蕉形”. 当孔隙中发生细颗粒迁移时，“低速区”极易形成颗粒堆积^[35].

图11(a)展示利用CFD建模计算，在微柱错开排列模式下，各向异性（旋转角 $ \theta $）对芯片模型渗透率的影响. 由图中曲线可以得到，在相同孔隙率下，旋转角 $ \theta $对芯片模型渗透率的有一定影响. 当 $ \theta =20{\text{°}} $时，微柱排列对流体的流阻最大，此时芯片模型的渗透率最小；反之当 $ \theta =70{\text{°}} $时，芯片模型的渗透性最大；在中间区域，旋转角 $ \theta $对芯片模型的渗透率几乎没有影响. 以孔隙率为0.60的模型为例，当 $ \theta =20{\text{°}} $时，芯片模型的渗透率为 $ 5.69\times {10}^{-11}\;{\rm{m}}^{2} $；当 $ \theta =70{\text{°}} $时，芯片模型的渗透率为 $ 9.27\times {10}^{-11}\;{\rm{m}}^{2} $；当 $ \theta \mathrm{在}30{\text{°}} \sim 60{\text{°}} $时，芯片模型的渗透率为 $ 8.0\times {10}^{-11}\;{\rm{m}}^{2} $左右，该值一般作为特征值为采用. 3种孔隙率下旋转角 $ \theta $产生的渗透率相对于特征值的幅值分别为 $28.87{\text{%}} \left({\varepsilon }=0.60\right)\mathrm{、} 15.39{\text{%}}\left({\varepsilon }=0.70\right) \mathrm{和}7.2{\text{%}}\left({\varepsilon }=0.80\right)$. 孔隙率越小，幅值越大，即旋转角在较大和较小时对芯片模型的渗透率影响越大，这也为在实际工程中，颗粒填充提供新的视角.

相对于颗粒排列，孔喉直径也常被用于研究对多孔介质渗透率的影响. 在旋转角 $ \theta $变化的过程中，芯片模型微柱间的孔喉直径也相应发生变化，随着 $ \theta $的增大而增大. 图11(b)展示芯片模型渗透率随孔喉直径的变化规律，拟合曲线为对数关系. 3种孔隙率下的曲线数据均表明，芯片模型渗透率随着孔喉直径的增大而增大；拟合曲线的系数为正，表明渗透率随着孔喉增大而增大的趋势减小. 在不考虑优势流通道的前提下，相对于大孔喉、小孔喉对多孔介质的渗透性更具有决定性作用.

为了进一步探究流体在芯片模型中的流动特性，图12(a)展示芯片模型平均迂曲度随旋转角 $ \theta $的变化规律. 曲线表明迂曲度在 $ \theta =20{\text{°}} $时最大，在 $ \theta =70{\text{°}} $时最小. 在旋转角度相同情况下，孔隙率越小，迂曲度越大. 图12(b)所示为芯片模型平均迂曲度值与相关文献中的数据对比分析. Comiti等^[36]提出的模型与旋转角较小时芯片模型的平均迂曲度值吻合较好.当旋转角较大时，芯片模型平均迂曲度值落在最下方，表明3种模型都不能单独预测该排列情况下的平均迂曲度值；当旋转角在 $ 30{\text{°}}\sim 60{\text{°}} $时，芯片模型平均迂曲度值均落在文献[37-38]中，具有较强的一致性. 每个模型的假设条件不同，参考文献中提出的迂曲度模型的适用对象并不具有普适性. 后期的研究方向可以聚焦于无序排列和三维颗粒堆积方面的实验和模拟.

针对均匀颗粒周期性排列的准二维微流控芯片模型的渗透性开展一些列实验和数值模拟研究，得到以下结论.

1）在相同孔隙率下，相比于微柱方形排列的芯片模型，微柱错开排列的芯片模型的渗流路径更长，微观上表现在迂曲度值更大，宏观上表现为更大的流阻和压降，最终体现为渗透率更小.

2）KC公式不适用于预测准二维微流控芯片模型的渗透率，而提出的修正平板流渗流公式能准确预测准二维微流控芯片模型渗透率.

3）实验数据和数值模拟结果表明，微柱直径和排列方式（各向异性）也是影响芯片模型渗透率大小的因素.

4）数值模拟结果表明，在各向异性下，旋转角在 $ 20\text{°}\mathrm{和}70\text{°} $时，对均匀颗粒等间距排列的芯片模型的渗透率影响最大，而在中间角度（ $ 30\text{°}\sim 60\text{°} $）下芯片模型的渗透率基本一致. 数值模拟结算得到的芯片模型平均迂曲度值也与文献中的参考数据吻合较好.