基于细观有限元分析的砂样等效渗透性预测.pdf

2023年6 月30 日铁路地质与路基2023年第1期总第10 5期基于细观有限元分析的砂样等效渗透性预测王梦辉（中铁第四勘察设计院集团有限公司，湖北武汉430 0 6 3）【摘要砂土的宏观渗透性和细观孔隙水平的真实流速对设计反滤层具有重要意义。采用有限元方法离散不可压缩流体流动的纳维一斯托克斯方程，采用不可压缩单元来模拟不可压缩流体，引入周期性边界条件（periodicboundaryconditions,PBC）来降低模型尺寸大小的影响。采用该方法计算了X射线微层析成像（X-ray microtomography，XM T）技术重建的砂样和随机生成的砂样模型的渗透性，并与实测值和经验值进行了对比。结果表明：在砂样渗流过程中，孔隙流体绕过固体颗粒流动，细观孔隙流体流速分布并不均匀，而是存在优势水力路径。【关键字有限元分析表征体元周期性边界条件孔隙流体流动在岩土工程中，砂土的宏观渗透性和细观孔隙水平的真实流速是设计反滤层的重要因素，反滤层对防止土体发生流土和管涌破坏具有重要意义。宏观分析通常用于研究土体中的平均孔隙压力和流体速度。然而，仅仅通过宏观研究无法获得渗流过程中孔隙细观层面上某个具体点的流体流速和压强，关于渗流机理的研究也需要对细观层面上的流体流动进行详细地分析。细观层面上，砂样是天然的异质材料，其主要由较密实的固体颗粒和连通孔隙组成。在梯度压力作用下，流体主要从连通的孔隙通过，对孔隙中流体流动规律的研究是了解砂样宏观渗透性的基础。砂样是一种典型的多孔介质。多孔介质在宏观尺度上常表现为连续介质，特征尺寸在mm量级及以上，细观尺度表现为孔隙和固体介质两种材料组成的异质结构，特征尺寸在mmm量级。多孔介质内部孔隙流体的流动特性涉及到多孔介质的细观结构以及细观尺度流体的运动规律。基于此，众多学者在这方面做了相关研究，研究结果表明在细观孔隙层次上多孔介质孔隙内部流体的力等效渗透性流动特性仍符合经典流体力学理论2 ，只要求解不可压缩粘性流体流动的纳维斯托克斯方程便可得到细观孔隙流体流动的规律。大量学者采用数值方法求解不可压缩粘性流体流动的纳维斯托克斯方程，研究复杂形态孔隙中细观流体的流动。Karim等3 采用XMT技术重建了砂样的三维真实微结构，并采用有限元法计算其渗透率。Pazdniakou等4采用LBM对平行板间Poiseuille流动进行了验证，得到了与解析解一致的结果；计算了重构和真实多孔介质（砂岩）在较大频率范围内的动态渗透率，验证了普遍的标度行为，实测值与理论预测有很好的一致性。咎鹏等5 采用LBM 的有限体积形式计算了完全多孔Poiseuille流和Couette流，数值计算结果与理论解吻合较好。Bignonnet等提出了一种计算复杂微结构渗透率的快速傅里叶变换（fastFouriertransform，FFT）方法，该方法直接与多孔介质样品的3D成像技术相耦合，不需要网格划分和定义等效孔隙网络，大大缩短了计算时间。田志国等7 采用离散单元法（D EM）模拟颗粒堆积结构冷压过程，再23用LBM进行渗流模拟，研究了球形颗粒在不同堆积形式下的渗透率。郑晓培8 采用COMSOLMultiphysics有限元软件求解纳维-斯托克斯方程，研究了多孔介质的孔径、孔隙率、固体颗粒骨架结构对渗透性的影响。在求解不可压缩粘性流体流动的纳维一斯托克斯方程的数值方法中，有限单元法具有求解稳定和精度高的特点，可以离散含复杂孔隙网络的砂样模型，并求解相应的细观孔隙流体流动，是预测砂样等效渗透性的可靠手段。在进行细观孔隙流体流动数值分析时，由于计算规模的限制，往往只能选取一定大小的表征体元（representative volume element,RVE）进行建模求解。然而，实际的孔隙结构是无限大的，在表征体元边界上速度自由，计算得到的流速往往与表征体元的大小有关。在表征体元的边界上施加周期性边界条件，计算结果则与表征体元的大小无关，在用有限大小的表征体元中的流体流动描述无限大的孔隙结构中的流体流动时更加精确。本文采用有限元方法离散不可压缩流体流动的纳维斯托克斯方程，采用ABAQUS软件自定义不可压缩单元来模拟不可压缩流体，引人了周期性边界条件来分析砂样的渗流，计算了XMT技术重建的砂样和随机生成的砂样模型的渗透性。1孔隙流体流动的有限元分析方法1.1控制方程渗流过程中，孔隙中的流体往往处于低速流动状态，忽略流体流动的惯性项，仅考虑粘性力作用，则不可压缩粘性流体稳态流动的纳维斯托克斯方程为1p+UVu.=0pax1+UVu,=0pay1Qp+UVu,=0Poz24式中，J，分别表示流体沿x，y，z方向的单位质量上所受的力，p为流体密度，p(x,y,z)为点(x,y,z)处的压强，=v(x,y,z)为流体或者固体介质的运动粘度，(x,y,2)=L(t,y,2)p(s,y,)A(,J,2)为流体或固体介质的动力粘度，在流体区域(x,J,z)=，在固体区域(x,y,z)=，u x，u y，u,分别表示流体沿x，J，z 方向的流速，t代表某一瞬时时刻。为拉普拉斯（Laplace）算子符号假定孔隙结构为周期性分布且已知（通过试验获得或者通过相关的算法生成），从已知的孔隙结构中提取周期性分布的表征体元。假定表征体元边界上的流速周期性分布，在宏观上施加压强梯度，孔隙流体流动引起的附加压强也呈周期性分布，不可压缩粘性流体稳态流动的纳维斯托克斯方程可化简为Vu-VP-Vp=0式中，VP为施加的压强梯度，Vp是周期性分布的附加压强梯度，u为周期性分布的流体速度。真实压强为p=P+p不可压缩流体运动的连续性微分方程为Vu=0(3)式中，u代表流速u的散度。1.2有限元离散由纳维斯托克斯方程（2）和不可压缩流体运动的连续性微分方程（3）组成的基本方程组可以求解不可压缩流体运动的速度场ux，u y，u,和压强场p。采用4节点（二维问题）或8 节点（三维问题）等参数单元对速度和周期性分布的附加压强进行插值，可得到单元内任一点的速度和附加压强的表达式，将离散的速度和附加压强代人到（2）和（3）方程的弱形式(1)中得KL-GT（2)F(4)式中，K为刚度矩阵，G为耦合矩阵、F为外力等效节点力向量，u为节点流速所组成的向量，为节点附加压强所组成的向量。其中刚度矩阵，耦合矩阵、外力等效节点力向量的具体表达式为K,=JVN,VN,dQaN-N,dQOxapF-Ndax.式中，K,为第j个节点对i个节点的刚度系数；N,为第i个节点速度和周期性压强的等阶插值函数；G,为第j个节点对i个节点的aP耦合系数；为第i个节点r方向的外力等Ox,aP效节点力；是作用在r方向的压强梯度；2为求解域。公式（4）的数值稳定形式为K-G-G-L-QP=0-Q-M（元）式中，L为周期性压强对应的刚度矩阵，Q为周期性压强和稳定变量的耦合刚度矩阵，M为稳定变量对应的刚度矩阵，元是为保证数值求解稳定而引人的一个伪变量。消去元，式（8）可以简化为K-GL-G(-L-QMQ)小 LP)L,=IN,VN,dQaN2NdoxM,=N,ON,dQ式中，L，为周期性压强对应第j个节点对i个节点的刚度系数；，为周期性压强和r方向稳定变量对应第j个节点对i个节点的刚度系数；M,为稳定变量对应第j个节点对i个节点的刚度系数；为稳定因子，大小取为单元的面积或体积除以流体的动力粘度。基于公式（9）（12），本文编写了(5)ABAQUS子单元UEL程序实现一维、二维和三维孔隙中流体流速的有限元求解。(6)1.3周期性边界条件的施加平行等宽裂隙介质和由规则排列的圆形(7)或方形颗粒组成的多孔介质是多孔介质细观结构最简单的假设。从多孔介质中选取周期性分布的表征体元，如图1和图2 所示。图1平行等宽裂隙介质0uF0(8)F(9)0(10)(11)(12)h-ddX图2 规则排列的圆形颗粒组成的多孔介质在表征体元的边界上，流速和附加压强也呈周期性分布。因此，在计算过程中施加满足式（13）（16）的周期性边界条件，即可用表征体元内的流体流速分布来表征整个多孔介质内的流体流速分布。在采用有限元软件ABAQUS的计算过程中，本文采用编写python脚本对模型相应节点施加多点约束方程的方法来施加周期性边界条件。u(x+h,y)=u(x,y)p(x+h,y)=p(x,y)u(x,y+h)=u(x,y)(13)(14)(15)25p(x,y+h)=p(x,y)1.4均质化理论宏观层面上往往将多孔介质看作是均质材料，即细观上由孔隙和固体介质组成的异质材料的物理量经过均质化（统计平均）处理后可得到宏观的均质物理量。达西定律本质上是一个宏观层面上的力学定律。求得表征体元内任一点的流速后，通过统计平均得到整个表征体元的平均流速，该过程即为均质化过程。再由平均流速和所施加的宏观压强梯度的关系便可确定多孔介质的宏观等效渗透性。方程（9）（12）的有限元系统求解完毕后，即可得到周期性表征体元内任意点的流速。将所有单元形心处的流速采用面积或体积进行加权平均得到表征体元各个方向的平均流速，然后即可通过达西定律预测i方向的绝对渗透率K,(m或 pixel)K,=-LhAP式中，i是流体运动的方向，h,是流体沿i方向流经路径的总长度（表征体元i方向的长度），,是表征体元i方向的流体平均流速，P是i方向的压强差。沿i方向的水力传导系数，即渗透系数k,(m/s)，可通过下式计算k,=K,之式中，为孔隙流体的重度。2有限元模型的验证2.1平行裂隙介质渗透率的有限元分析平行裂隙介质和规则排列的圆形颗粒所组成的多孔介质的绝对渗透率分别有相应的解析解。对于平行裂隙介质，如果裂隙的宽度是h-a，裂隙的间距为，如图1所示，Al-Omari等9给出的平行裂隙介质的绝对渗透率为26(16)(17)(18)d(h-a)2K=12式中，为平行裂隙介质的孔隙率h-ad=h取h为1cm，有限元法计算的平行裂隙介质的绝对渗透率结果与解析解的对比如图3所示，有限元法计算的绝对渗透率与解析解一致。10-1,有限元法10-210-3101010-610-710-80.0图3平行裂隙介质的绝对渗透率2.2圆形阵列介质渗透率的有限元分析对于规则排列的圆形颗粒组成的多孔介质，如果相邻圆形颗粒圆心的间距为h，颗粒的直径为d，如图2 所示，Masad等1I给出的规则排列的圆形颗粒组成的多孔介质的绝对渗透率为2/2(h-d)s2K=9元h/2hln g.-0.738+.-0.887g+2.038g,4元(21)式中，为规则排列的圆形颗粒组成的多孔介质的固体体积分数元d?4h2取h为1cm，采用有限元法计算的规则排列的圆形颗粒组成的多孔介质绝对渗透率结果与解析解的对比如图4所示，有限元计算的绝对渗透率与解析解一致。本文有限元(19)(20)解析解0.20.4孔隙率/0.60.8(h-d h)(22)1.0方法可用于预测多孔介质的宏观等效渗透性。10-1一解析解x有限元法10-21010-10610-7102图4 规则排列的圆形介质的绝对渗透率3砂样的渗透性计算3.1基于XMT技术重建的砂样渗透性预测XMT技术被越来越多的用于确定多孔介质内部的细观结构图像，如水泥基材料、砂岩、土等。Karim等3.1采用XMT技术获得的砂样细观固体介质和孔隙分布图像如图5所示，图像大小为4 0 0 4 0 0 pixels（像素），约为14 0 0 14 0 0 m。采用MATLAB软件读取该图像的灰度值进行建模，砂样模型如图6 所示。取流体粘度=1.010-3Pas，沿x方向施加压强梯度-10 0 0 0 Pa/m，采用本文有限元方法获得的砂样孔隙流体细观流速分布如图7 所示，为了便于观察孔隙流体流速，此处仅显示模型的孔隙流体部分。在砂样中，孔隙层面的流体流速分布并不均匀，而是存在优势水力路径。在压力差的驱动下，孔隙流体绕过固体颗粒主要沿优势水力路径发生流动。孔隙中流体最大流速位置位于优势水力路径中，优势水力路径的位置则与孔隙的细观结构有关。图5采用XMT技术重建的砂样3 0.40.6孔隙率/0.81.0图6 砂样有限元模型流速/（mms)1.3971.2801.1641.0480.9310.8150.6980.5820.4660.3490.2330.1160.000图7 砂样孔隙内流体细观流速分布求得砂样模型内任一点的流体流速后，采用均质化理论计算砂样的绝对渗透率，计算结果为4.2 9 10-pixel。K a r i m等3 所给出的计算结果和实测值分别为3.5510-2pixel?和7.0 0 10 pixel，本文计算结果介于二者之间，并且本文所预测的结果更接近实测值，数值预测的结果与实测值在同一数量级上。3.2随机生成的砂样模型渗透性预测二维模型只考虑了流体沿平面上两个维度的流动情况，然而在实际的渗流问题中流体是在空间三个维度上发生流动，三维模型更加符合砂样的真实情况。基于此本文建立了三维砂样模型计算其渗透性，以便更加真实的模拟通过砂样的渗流。真实的砂土固体颗粒和孔隙的分布往往错综复杂，需要通过离散元或者异形颗粒离散元模拟沉积作用生成。本文忽略该过程，对真实的砂样进行了一定程度的简化。假设27砂样中的固体颗粒为球形颗粒，允许球形固体颗粒可以相互重叠用来表示形状更加复杂的固体颗粒，从而得到砂样的细观固体颗粒骨架和孔隙结构。为了解决ABAQUS软件在建立三维砂样有限元模型中面临的困难，对三维砂样渗透性的计算采用MATLAB程序生成ABAQUS计算输人文件，然后在ABAQUS软件内提交计算。建立边长为1mm的立方体表征体元，把立方体表征体元划分成10 0100 100的网格，共计10 0 0 0 0 0 个单元，每个单元的尺寸为0.0 1mm0.01mm0.01mm。采用MATLAB软件随机生成满足周期性边界条件的表征体元中球体颗粒的球心坐标，这些球体颗粒的直径都为0.2mm。通过计算单元形心到任一球形固体颗粒球心的距离来判断该单元是固体单元还是流体单元。对于单元形心在任一球形固体颗粒内部的单元认为该单元是固体颗粒单元，对该单元指派固体颗粒参数；对于单元形心不在任一球形固体颗粒内部的单元认为该单元是流体单元，对该单元指派孔隙流体参数。计算流体单元个数占表征体元中单元总个数的比率，该数值即为表征体元的孔隙率。通过对表征体元边界六个面上的相对应节点施加约束方程使表征体元边界上的流体流速满足周期性边界条件。生成的由随机分布的球形固体颗粒形成的三维砂样表征体元如图8 所示，该表征体元的孔隙率为0.3 7 5。隙流体重度y=9.8kN/m，计算该表征体元的渗透系数。该表征体元的渗透系数矩阵计算结果如下7.70-0.15-0.48k=-0.136.19-0.9010m/sL-0.47-0.935.16本模型的固体颗粒直径为0.2 mm，该粒径符合余钟波等12 所给出的细砂粒径范围。在该孔隙率下，用本模型近似表示细砂，余钟波等2 3 给出的细砂渗透系数经验值为2 10-7 2 10-4 m/s，由本模型计算所得渗透系数矩阵可见本模型渗透系数符合细砂渗透系数经验值范围。从表征体元的渗透系数矩阵可以看出，该砂样表征体元的渗透性具有轻微各向异性。沿x方向施加压强梯度-10 0 0 0 Pa/m，该表征体元的孔隙流体流速分布如图9 所示，为了清晰显示表征体元内部流体流速分布，将表征体元沿x,y，z 三个方向从中部切开。由图9 可见，孔隙流体流速的最小值位于球形固体颗粒的表面，孔隙流体流速的较大值位于表征体元中远离固体颗粒的位置，孔隙流体流速最大值的位置则与随机生成的砂样表征体元中球形固体颗粒的位置分布情况有关。流速/(mms)0.0000.001图9 三维砂样表征体元内流体流速分布4结论图8 随机生成的三维砂样表征体元采用有限单元法离散不可压缩粘性流取孔隙流体粘度u=1.010-3Pas，孔28体稳态流动的纳维一斯托克斯方程，使用ABAQUS自定义不可压缩单元来模拟不可压缩流体，引人周期性边界条件提高计算结果的精度，求得砂样表征体元内的渗流速度，采用均质化理论计算砂样的渗透性，得到以下结论：（1）砂样表征体元的渗透性具有轻微各向异性。（2）在砂样的渗流中，细观孔隙流体流速分布并不均匀，而是存在优势水力路径。孔隙中流体最大流速位于优势水力路径中，优势水力路径的位置与孔隙的细观结构有关。（3）在砂样表征体元中，固体颗粒表面孔隙流体流速最小，远离固体颗粒的位置孔隙流体流速较大。参考文献：1陈玉丽，马勇，潘飞，等.多尺度复合材料力学研究进展J.固体力学学报.2 0 18,3 9(1);1-68.2 李战华，周兴贝，朱善农.非极性小分子有机液体在微管道中的流量特性J力学学报。2002,34(3):432-438.3 Karim M R,Krabbenhoft K,Lyamin AV.Permeability determination of porous media usinglarge-scale finite elements and iterative solverJJ.International Journal for Numerical and AnalyticalMethods in Geomechanics.2014,38:991-1012.4Pazdniakou A,Adler P M.Dynamic permeability ofporous media by the lattice Boltzmann methodJ.Advances in Water Resources.2013,62:292-302.5 咎鹏，谢新生，陈燕萍.基于改进FV-LBM的多孔介质渗流模拟J.人民长江.2 0 2 0,51(2)：172-178.6Bignonnet F,Dormieux L.FFT-based bounds onthe permeability of complex microstructuresJ.International Journal for Numerical and AnalyticalMethods in Geomechanics.2014,38:1707-1723.7田志国，王沫然：冷压过程对球堆积结构渗透率的影响.空间控制技术与应用.2 0 2 0,46(4):72-75.8郑晓培.水泥土的渗透性及其计算模型D.南京：东南大学，2 0 18.9Al-Omari A,Masad E.Three dimensionalsimulation of fluid flow in X-ray CT imagesof porous mediaJ.International Journalfor Numerical and Analytical Methods inGeomechanics.2004,28:1327-1360.10 Masad E,Muhunthan B,Crowe C.Numericalmodelling of fluid flow in microscopic imagesof granular materialsJ.International Journalfor Numerical and Analytical Methods inGeomechanics.2002,26:53-74.1l KarimMR,Krabbenhoft K.Extraction ofEffective Cement Paste Diffusivities from X-rayMicrotomography ScansJJ.Transport in PorousMedia.2010,84:371-388.12余钟波，黄勇，SchwartzF.W.地下水水文学原理M.北京：科学出版社，2 0 0 8.29