高 笑，李宗利，李 洋

(西北农林科技大学 水利与建筑工程学院，陕西 杨凌 712100)

裂缝对渗流场影响的数值模拟方法研究

高 笑，李宗利，李 洋

(西北农林科技大学 水利与建筑工程学院，陕西 杨凌 712100)

【目的】 研究裂缝对渗流场影响的数值模拟方法，为岩石、混凝土等材料裂缝渗流场的变化研究提供支持。【方法】 在等效连续介质模型框架下提出渗透系数张量修正法，并基于流量等效原理推导出修正单元的渗透系数张量矩阵。应用有限元分析软件，分别建立裂缝位于构件边界和构件内部时的分析模型，应用渗透系数张量修正法得到的渗透系数张量表达式定义裂缝单元的渗透特性，将模拟所得渗流场与裂缝真实存在的数值分析结果进行比较，并通过改变修正单元的尺寸分析渗透系数张量修正法对单元尺寸的依赖性。【结果】 修正单元为10倍裂缝宽度时，与精确解相比，各分析模型所得渗流量及相应断面孔隙水压的误差均小于0.5%。就单元尺寸依赖性而言，为满足工程精度要求，从渗流量角度分析，修正单元的尺寸应不大于600倍裂缝宽度；从孔隙水压角度分析，修正单元的尺寸应不大于200倍裂缝宽度。【结论】 所建立的等效处理方法合理可行，当修正单元的尺寸不大于200倍裂缝宽度时，无论是渗流量还是孔隙水压，该处理方法均能满足工程精度要求。

裂缝；渗流场；数值模拟；渗透系数张量修正法；单元尺寸

高坝、地下洞室等结构的渗流场数值分析是目前研究的重要课题之一。岩体和混凝土自身存在着裂缝，这些裂缝控制着其周围的渗流场。目前处理裂缝效应的模型[1-2]主要有3种，分别是等效连续介质模型、离散介质模型和双重介质模型，这3种模型都只是考虑已有裂缝的渗流问题，却忽略了岩体或混凝土在力的作用过程中新生裂缝对其渗流的影响。尤其是在渗流-应力或损伤分析中，裂缝是动态发展的，势必引起渗流场的局部变化，严格地讲需要重新划分网格，但实际上很难做到。目前国内外就动态裂缝扩展对应力场影响的研究已有成熟方法，最具代表性的有Rots[3]基于旋转裂缝模型的劲度矩阵修正法及Bazant等[4]基于裂纹带模型的损伤折减法等，但是考虑动态扩展裂缝对渗流场影响的研究成果却很少。杨天鸿等[5-6]从细观力学角度，针对岩体裂缝引入突变系数修正裂缝单元的渗透性能，但并不能体现裂缝所产生渗流的各向异性，突变系数的取值也缺乏理论或试验支撑。

等效连续介质模型[7-9]将裂缝中的水流等效到研究体中，在分析过程中可不用重新划分网格。而渗透系数张量是将研究体作为等效连续渗透介质的重要参数，在很大程度上决定了渗流场分析的正确与否。对裂隙岩体渗透系数张量的研究现已较为成熟，Long等[10]应用连续介质理论得到了裂缝性岩体的等效渗透率张量,但却忽略了基岩的渗透性。李亚军等[11]依据流量等效原理，得到含有多条裂缝的等效渗透率张量的数学模型。刘建军等[12]基于流量等效原理，推导出含裂缝岩块的渗透系数张量。荣冠等[13]应用能量叠加原理得到节理岩体的渗透系数张量，虽强调了节理的连通性，但却忽略了基岩以及非连通节理对节理岩体渗透性的影响。姚军等[14]基于流量等效原理，应用边界元算法得到了裂缝性多孔介质网格块的等效渗透率张量。以上理论都是将裂缝的渗透性等效到整个研究体或表征体元中，并未从数值模拟的单元层次考虑裂缝对渗透性的影响。可见需结合考虑裂缝与基质体的渗透特性，研究裂缝单元的渗透系数张量，通过修正裂缝单元的渗透性，则可体现出裂缝对混凝土或岩体渗流场的影响。为此，本研究基于等效连续介质模型，利用流量等效原理，综合考虑基质与裂缝的渗透特性，修正裂缝单元的渗透系数张量，而数值分析模型并不重新划分网格，以此体现裂缝的存在及其对渗流场的影响，以期为考虑裂缝存在及其发展对岩石、混凝土等材料渗流场影响的数值分析研究提供参考。

1 裂缝单元的渗透系数张量

含有裂缝的单元，即修正单元，在产生裂缝之前其渗透系数张量为各向同性，之后变为各向异性，且以裂缝的扩展方向为渗流主方向(下文中的主方向均为裂缝扩展方向)。由于垂直裂缝方向的研究体结构变化很小，认为垂直渗流主方向的渗透系数保持不变，因此只需修正裂缝单元主方向的渗透系数。

基于流量等效原理，通过修正单元的流量，等于通过裂缝的流量与通过完整单元的流量除去裂缝流量所剩部分(下称残缺单元)流量之和。裂缝与单元相比很小，且完整单元的渗透系数也远小于裂缝的渗透系数，则可以用通过完整单元的流量代替通过残缺单元的流量，因此通过修正单元的流量等于通过裂缝与完整单元的流量之和。单元的流量以流出流量计算，任意单元的流量等于主方向(裂缝方向)和垂直主方向流量之和。又因垂直主方向的渗透系数保持不变，所以修正单元主方向流量等于裂缝流量与完整单元主方向流量之和，等效原理示意图如图1所示，图中1、2、3、4为单元节点编号，5、6为裂缝两端位置点编号。

裂缝修正单元的流量等效原理为：

(1)

式中：kζ′为修正单元等效之后主方向的渗透系数，m/s；l14、l34分别为节点1与4、3与4的距离，m；h14、h34分别为单元边界l14、l34位置处的平均水头，m；nζ为渗流主方向；-∂h14/∂nζ、-∂h34/∂nζ分别为单元边界l14、l34在主方向的水力梯度；θ1、θ2分别为主方向与l14、l34的夹角(≤90°)。

②如图1(b)所示，完整单元主方向ζ流出的流量qζ″为：

(2)

式中：kζ″为完整单元主方向的渗透系数，m/s；其余同上。

③如图1(c)所示，裂缝的流量qf为：

(3)

(4)

式中：kf为裂缝渗透系数，由立方定律[15-16]得到，m/s；h5、h6分别为由节点水头插值所得裂缝扩展末端、扩展始端的水头，m；l56为裂缝长度，m；b为裂缝宽度，m；γ为水的重度，N/m3；μ为水的黏度，(N·s)/m2。

④由qζ′=qζ″+qf，得：

kζ′=kζ″+kf[b(h6-h5)/l56]/

(l14sinθ1·∂h14/∂nζ+l34sinθ2·∂h34/∂nζ)。

(5)

由于修正单元在等效之后，可视为均匀材料，且单元相对较小，修正单元各位置在裂缝扩展方向的水力梯度以及裂缝在其扩展方向的水力梯度相差很小，可取相等，则(5)式可表示为：

kζ′=kζ″+kfb/L。

(6)

式中：L=l14sinθ1+l34sinθ2，为位于裂缝两侧且距裂缝最远两节点到裂缝的距离之和。

⑤由(6)式求出修正单元渗流主方向渗透系数kζ′，因垂直裂缝扩展方向的材料结构特性未发生变化，则垂直渗流主方向渗透系数kζ′不变，取混凝土渗透系数，即kη′=kζ″，则可得修正单元在局部坐标系(ζ，η)下的渗透系数张量如(7)式所示。经坐标转化，则可得修正单元在整体坐标系(x，y)下的渗透系数张量如(8)式所示。

(7)

(8)

式中：kxy=kyx。λ可按下式计算：

(9)

式中：α为整体坐标系x轴与局部坐标系ζ轴的夹角。

2 算例验证

为验证以上理论的正确性与精确度，用有限元大型商业软件ABAQUS对图2两算例进行数值模拟分析。算例1裂缝位于混凝土构件边界；算例2裂缝位于混凝土构件内部。裂缝宽为0.08 mm，长为0.5 m，与水平方向夹角α=36.87°，裂缝位置(A、B)、荷载(p1、p2)及分析断面Ⅰ-Ⅰ(x=0.2 m)、Ⅱ-Ⅱ (x=0.4 m)、Ⅲ-Ⅲ(x=0.7 m)、Ⅳ-Ⅳ(x=0.3 m)、Ⅴ-Ⅴ(x=0.7 m)的位置等详见图2。

算例1，2模型参数：构件尺寸1 m×1 m；混凝土密度ρ=2 400 kg/m3，弹性模量E=20 GPa，泊松比υ=0.167，除修正单元之外的混凝土渗透系数k=1×10-10m/s，各向同性；裂缝渗透系数各向异性，主方向渗透系数kf=4.589×10-3m/s(由(4)式求得)，取15°水的黏度μ=0.114×10-4(N·s)/m2、重度γ=9 810 N/m3，垂直主方向的渗透系数取混凝土的渗透系数(考虑了裂缝与混凝土接触面的渗流，取混凝土渗透系数更为合理)；构件左侧、右侧分别作用1 000和100 Pa的均匀水压力，并以孔隙水压力形式施加于模型相应边界；左、右侧为透水边界，上、下侧为不透水边界；固定构件的水平位移、竖向位移及角位移；修正单元、较大规则单元均采用孔压/位移耦合的CPE8RP四边形八节点平面应变单元，过渡单元采用孔压/位移耦合的CPE6MP三角形六节点平面应变单元；用Soils(固结)分析步进行稳态渗流分析。本研究只分析裂缝对渗流场的影响，不考虑重力对渗流场的影响。

首先，分别建立算例1、2裂缝真实存在的模型，将裂缝视为一种嵌入在非裂缝单元边界的均质材料单元，其主方向渗透系数按式(4)计算，单元的尺寸为裂缝宽度，这样处理实际上是精确计算，可作为精确解，为下文中算例1、2含有修正单元模型的模拟结果的检验标准。算例1有42 020个单元、86 669个节点，算例2有37 902个单元、78 423个节点。然后，分别建立算例1、2修正单元为10倍裂缝宽度的模型，即N=L/b=10，修正单元主方向渗透系数为4.589×10-6m/s，算例1有1 626个单元、3 625个节点，算例2有1 910个单元、4 189个节点。算例1、2精确解模型与相应含有修正单元模型的网格划分、单元分布以及模拟所得渗流场见图3、4。

算例1含有修正单元的模型模拟所得单宽流量为1.182×10-11m3/s，相应的精确解为1.179×10-11m3/s。算例2含有修正单元的模型模拟所得单宽流量为1.049×10-11m3/s，相应的精确解为1.047×10-11m3/s。经计算可得算例1、2修正单元为10倍裂缝宽度的模型模拟所得单宽渗流量q与其相应精确解相比，误差分别为0.29%和0.21%，都远小于工程允许误差(5%)。从图3(a)、(b)和图4(a)、(b)可以看出，算例1、2中修正单元为10倍裂缝宽度模型模拟所得渗流场与其相应精确解整体上差别很小。算例1中断面Ⅰ-Ⅰ、Ⅱ-Ⅱ和Ⅲ-Ⅲ孔隙水压的对比见图5，算例2中断面Ⅳ-Ⅳ和 Ⅴ-Ⅴ 孔隙水压的对比见图6。

由图5可以看出，在断面Ⅰ-Ⅰ、Ⅱ-Ⅱ和Ⅲ-Ⅲ，算例1修正单元为10倍裂缝宽度模型模拟所得孔隙水压与精确解相应位置的孔隙水压相比，最大误差依次约为0.12%，0.47%,0.32%，均远小于工程允许误差。由图6可以看出，在断面Ⅳ-Ⅳ和Ⅴ-Ⅴ，算例2修正单元为10倍裂缝宽度模型模拟所得孔隙水压与精确解相应位置的孔隙水压相比，最大误差分别约为0.29%和0.38%，远小于工程允许误差。可见，无论裂缝位于构件边界还是内部，应用式(5)～(9)对裂缝单元的渗透系数张量进行等效处理，用以分析裂缝对渗流场的影响，该方法合理可行。

3 单元大小局限性分析

以上两算例虽验证了应用式(5)～(9)对修正单元的渗透系数张量进行等效处理的可行性，但其修正单元为裂缝宽度10倍这一数量级，修正单元相对较小。为此，将算例1、2中的修正单元依次放大到裂缝宽度的100，200，300，400，600，700和800倍，分别从渗流量和孔隙水压角度，研究单元大小对该等效处理方法的影响。

算例1、2中相应修正单元主方向渗透系数见表1，各模型建立方法(包括网格划分、单元分布)与修正单元为10倍裂缝宽度类似，且渗流场分布与相应精确解相似，这里不再重复。依据分析的需要，图7仅给出算例1部分模型断面，即断面Ⅰ-Ⅰ、Ⅱ-Ⅱ和Ⅲ-Ⅲ的孔隙水压对比结果，图8仅给出算例2部分模型断面，即断面Ⅳ-Ⅳ和Ⅴ-Ⅴ的孔隙水压对比结果。算例1、2各模型模拟所得的渗流量如表1所示。

由表1可得，算例1修正单元为100，200，300，400，600，700，800倍裂缝宽度时，各模型的渗流量与其精确解相比，误差分别为1.29%，2.05%，2.67%，2.86%，4.37%，4.93%和5.44%；算例2修正单元为100，200，300，400，600，700倍裂缝宽度时，各模型的渗流量与其精确解相比，误差分别为1.28%，2.05%，2.46%，3.10%，4.55%和5.10%。可见仅从渗流量角度考虑，要满足工程精度(误差<5%)要求，当裂缝位于构件边界时，修正单元应不大于700倍裂缝宽度；当裂缝位于构件内部时，修正单元应不大于600倍裂缝宽度；当构件边界和内部都有裂缝时，修正单元应不大于600倍裂缝宽度。

由图7可以看出，算例1修正单元为300，400倍裂缝宽度的模型与其相应位置的精确解相比，在断面Ⅰ-Ⅰ,最大误差依次约为1.19%，1.32%；在断面Ⅱ-Ⅱ，最大误差依次约为4.86%，5.86%，且均在裂尖附近位置但并不在裂尖；在断面Ⅲ-Ⅲ，最大误差依次约为2.88%，3.15%。由图8可知，算例2修正单元为200，300倍裂缝宽度的模型与其相应位置的精确解相比，在断面Ⅳ-Ⅳ，最大误差依次约为3.52%，4.11%；在断面Ⅴ-Ⅴ，最大误差依次约为4.78%，5.61%。可见仅从孔隙水压角度考虑，要满足工程精度(误差<5%)要求，当裂缝位于构件边界时，修正单元应不大于300倍裂缝宽度；当裂缝位于构件内部时，修正单元应不大于200倍裂缝宽度；当构件边界和内部都有裂缝时，修正单元应不大于200倍裂缝宽度。

综上可知，同时考虑渗流量和孔隙水压，应用式(5)～(9)对不大于200倍裂缝宽度的修正单元进行等效处理，可以满足工程精度要求。

4 结 论

1)在荷载作用下，构件边界或内部产生裂缝，在等效连续介质模型不变网格框架下，依据流量等效原理修正裂缝单元的渗透系数张量以体现裂缝对构件渗流场的影响，是合理可行的。

2)本研究提出的渗透系数张量修正法受到修正单元大小的影响。仅从渗流量角度，修正单元尺寸应不大于600倍裂缝宽度；仅从孔隙水压角度，修正单元尺寸应不大于200倍裂缝宽度；二者同时考虑，则要求修正单元不大于200倍裂缝宽度。

[1] 宋晓晨,徐卫亚.裂隙岩体渗流概念模型研究 [J].岩土力学,2004,2(2):226-231.

Song X C,Xu W Y.A study on conceptual models of fluid flow in fractured rock [J].Rock and Soil Mechanics,2004,2(2):226-231.(in Chinese)

[2] 祝云华,刘新荣,梁宁慧,等.裂隙岩体渗流模型研究现状与展望 [J].工程地质学报,2008,16(2):178-183.

Zhu Y H,Liu X R,Liang N H,et al.Current research and prospects in modeling seepage field in fractured rock mass [J].Chinese Journal of Engineering Geology,2008,16(2):178-183.(in Chinese)

[3] Rots J G.Computational modeling of concrete fracture [M].Delft:Technische Hogeschool Delft,1988.

[4] Bazant Z P,Oh B H.Crack band theory for fracture of concrete [J].Materials and Structure,1983,16:155-177.

[5] 杨天鸿,唐春安,朱万成,等.岩石破裂过程渗流与应力耦合分析 [J].岩土工程学报,2001,23(4):489-493.

Yang T H,Tang C A,Zhu W C,et al.Coupling analysis of seepage and stress in rock failure process [J].Chinese Journal of Geotechnical Engineering,2001,23(4):489-493.(in Chinese)

[6] 杨天鸿,屠晓利,於 斌,等.岩石破裂与渗流耦合过程细观力学模型 [J].固体力学学报,2005,26(3):333-337.

Yang T H,Tu X L,Yu B,et al.A micromechanical model for simulating the coupling of fracture and flow of rock [J].Acta Mechanica Solida Sinica,2005,26(3):333-337.(in Chinese)

[7] 张有天,张武功.裂隙岩石渗透特性渗流数学模型及系数量测 [J].岩石力学,1982(8):41-52.

Zhang Y T,Zhang W G.Seepage mathematical model and coefficient measurement of fractured rock [J].Chinese Journal of Rock Mechanics,1982(8):41-52.(in Chinese)

[8] Oda M.An equivalent continuum model for coupled stress and fluid flow analysis in jointed rock masses [J].Water Resources Research,1986,22(13):1845-1856.

[9] Snow D T.Anisotropic permeability of fractured media [J].Water Resources Research,1969,5(6):1273-1289.

[10] Long J C S,Remer J S,Wilson C R,et al.Porous media equivalents for networks of discontinuous fractures [J].Water Resources Research,1982,18(3):645-658.

[11] 李亚军,姚 军,黄朝琴,等.裂缝性油藏等效渗透率张量计算及表征单元体积研究 [J].水动力学研究与进展,2010,25(1):1-7.

Li Y J,Yao J,Huang Z Q,et al.Calculation of equivalent permeability tensor and study on representative element volume for modeling fractured reservoirs [J].Chinese Journal of Hydrodynamics,2010,25(1):1-7.(in Chinese)

[12] 刘建军,刘先贵,胡雅礽,等.裂缝性砂岩油藏渗流的等效连续介质模型 [J].重庆大学学报：自然科学版,2000,23(S1):158-160.

Liu J J,Liu X G,Hu Y R,et al.The equivalent continuum media model of fracture sandstone reservoir [J].Journal of Chongqing University：Natural Science Edition，2000,23(S1):158-160.(in Chinese)

[13] 荣 冠,周创兵,王恩志.裂隙岩体渗透系数张量计算及其表征单元体积初步研究 [J].岩石力学与工程学报,2007,26(4):740-746.

Rong G,Zhou C B,Wang E Z.Preliminary study of permeability tensor calculation of fractured rock mass and its representative elementary volume [J].Chinese Journal of Rock Mechanics and Engineering,2007,26(4):740-746.(in Chinese)

[14] 姚 军,李亚军, 黄朝琴,等.裂缝性油藏等效渗透率张量的边界元求解方法 [J].油气地质与采收率,2009,16(6):80-83.

Yao J,Li Y J,Huang Z Q,et al.Calculation of equivalent permeability tensors of fractured reservoirs using boundary element method [J].Petroleum Geology and Recovery Efficiency,2009,16(6):80-83.(in Chinese)

[15] Bear J.Dynamics of fluids in porous media [M].New York:Elsevier,1972.

[16] 田开铭,万 力.各向异性裂隙介质渗透性的研究与评价 [M].北京:学苑出版社,1989.

Tian K M,Wan L.Research and evaluation of the permeability of anisotropic fractured media [M].Beijing:Academy Press,1989.(in Chinese)

Numerical simulation method for effect of cracks on seepage field

GAO Xiao,LI Zong-li,LI Yang

(CollegeofWaterResourcesandArchitecturalEngineering,NorthwestA&FUniversity,Yangling,Shaanxi712100,China)

【Objective】 This study focused on numerical simulation method for the impact of cracks on seepage field to support the research of variations in seepage field of rock or concrete caused by cracks.【Method】 Based on the equivalent continuum model and equivalent flow principle, modified permeability tensor method was proposed and the permeability tensor matrix of the modified element was derived.Then,using finite element software,numerical model with a boundary crack or an inner crack was built.The seepage fields of established model were compared to models with real cracks and element size was modified to analyze the effect of element size on this method.【Result】 Compared with the exact solution,the error of the seepage flux and the corresponding pore water pressure of the models with modified elements of 10 times of the crack width was less than 0.5%.To meet the precision requirement in engineering,from the view of seepage flux,fixed element size should be no more than 600 times of the crack width.From the view of pore water pressure,fixed element size should be no more than 200 times of the crack width.【Conclusion】 The proposed method was reasonable and feasible.Considering both seepage flux and pore water pressure,the proposed method can meet the engineering precision requirements when modified element size was no more than 200 times of the crack width.

cracks;seepage field;numerical simulation;modified permeability tensor method;element size

2013-12-13

国家自然科学基金项目(51379178，50779057)；中央高校基本科研业务费专项(ZD2012015)

高 笑(1987-),男，陕西佳县人，在读硕士，主要从事水利水电工程设计研究。E-mail:gaoxiao1119@sina.cn

李宗利(1967-),男，陕西凤翔人，教授，博士生导师，主要从事水工结构与岩土工程稳定分析理论、水工金属结构设计理论与方法研究。E-mail:zongli02@163.com

时间:2015-04-13 12:59

10.13207/j.cnki.jnwafu.2015.05.017

TV223.4

A

1671-9387(2015)05-0222-07

网络出版地址:http://www.cnki.net/kcms/detail/61.1390.S.20150413.1259.017.html