ABAQUS渗流应力耦合分析中渗透荷载施加问题探讨

，，

(1.三峡大学 水利与环境学院，湖北 宜昌 443002；2.西北农林科技大学 水利与建筑工程学院，陕西 杨凌 712100)

1 研究背景

目前，关于岩体或混凝土渗流应力耦合问题的研究较多[1-4]，鉴于理论分析的繁琐和试验模拟的局限，数值分析成为重要的研究手段，然而不同的软件对于水力边界处理方法不同。ABAQUS作为一款计算功能强大的有限元分析软件，提供了渗流应力耦合分析功能，能够求解多孔介质的饱和渗流、非饱和渗流，以及二者的混合问题[5-7]。该软件的渗流应力耦合分析理论是在Biot固结理论的基础上建立起来的，可以直接考虑骨架颗粒变形与孔隙水之间的相互作用，即“直接耦合”；用户还可以通过编写子程序的方式定义渗透系数与应力、应变等力学参数之间的关系，进而实现直接耦合过程中引起的水力参数和力学参数的变化，即“间接耦合”；因此在利用该软件进行完全渗流应力耦合分析时包含了“直接耦合”和“间接耦合”2种情况。ABAQUS中专门设有孔压单元，较常规应力位移单元多一个孔隙水压自由度，能够方便进行渗流应力耦合分析。然而在模型上下游边界指定随高程线性变化的孔压来满足水头条件的同时，是否还需要再次对该边界定义相应的静水压力荷载，一直是学者不断争论的问题之一。

一些学者认为只要定义了孔压边界，软件将自行根据水头计算出模型内渗透力，然后按照体力考虑水荷载的作用[8]；而有的学者则认为由于材料存在孔隙，故水荷载不可能百分之百作用在边界面上，因此从安全角度出发，认为应该再次施加静水压力荷载[9-10]。另外，利用该软件对混凝土重力坝进行渗流应力耦合分析时，如果考虑了坝基的渗透作用，是否需要再单独定义坝体扬压力荷载也成为学者争论的问题。在不考虑坝体透水性时，林凯生等[11]认为由于考虑坝基渗透性，坝基面必存在孔隙水压力，因此不需要人为定义坝体扬压力；而常晓林等[8]认为建基面的孔隙水压力并不能传递到坝体，需借助ABAQUS的Dload子程序在坝底面定义扬压力荷载。可见，对上述2个争论问题尚需进一步的探讨。

本文从理论和数值模拟2个方面出发，并结合实例分析来探讨应用ABAQUS进行混凝土重力坝渗流应力耦合分析时渗流荷载和水边界条件的施加问题，旨在为研究此类问题提供参考。

2 ABAQUS渗流应力耦合原理

2.1 有效应力方程

基于有效应力原理，总应力σ由有效应力σ′、孔隙水压力uw和孔隙气压力ua共同组成，其表达式为

σ=σ′+[χuw+(1-χ)ua]I。

(1)

式中：χ为饱和度s的函数；I为单位矩阵。当材料完全饱和时，χ=1；当材料非饱和时，χ=χ(s)。

2.2 应力平衡方程

固相材料的应力平衡方程由虚功原理表示，即在某一时刻t, 固相材料的虚功与作用在该材料上的作用力(面力和体力)产生的虚功相等[6]，即

式中：δv为虚位移场；δε为虚应变,δε=sym(∂δv/∂x)；T为单位面积的表面力；f为单位体积的体积力(不含流体重量)；n为孔隙率；V为体积；S为面积边界；ρw为流体的密度；g为重力加速度。

2.3 渗流连续性方程

流经dV的流体应满足连续性方程，使得在某时间增量内流入的流体流量等于流体体积的增加速率，其表达式[5]为

(3)

利用ABAQUS进行渗流应力耦合分析时，无需进行渗流场与应力场的反复迭代，只要按时间过程连续求解便可得到全部结果。即在上述应力平衡方程和渗流连续方程的基础上，将节点位移和孔隙水压力作为节点自由度进行空间离散，将应力平衡方程和渗流连续方程写成矩阵形式，并在同时满足位移边界和渗流边界的条件下，对每个时间步内的方程进行求解。

图1 计算模型Fig.1 Computation model

3 孔压边界和静水压力荷载施加问题

以图1所述工程问题为例来进行理论推导，并采用数值分析，比较二者结果，验证不同水力边界设置的影响。某岩基处于水下15 m，岩体的密度为2 000 kg/m3，弹性模量为27 GPa，泊松比为0.2，渗透系数为1×10-8m/s[3]，孔隙率为0.091，岩基初始饱和度为1。数值模型取计算范围为100 m×100 m。

首先，根据竖向应力平衡和孔隙静水压力平衡推导得竖向应力的计算公式为

S22=ρg(z0-z)-γws(1-n)(z0-z)=

[ρg-γws(1-n)](z0-z)=γ浮(z0-z) 。(4)

式中：S22为竖向有效应力；ρ为岩体的密度；g为重力加速度，取为10 m/s2；z0为岩基表面至计算面的深度，z0-z为岩基表面以下埋深；γw为水的重度；γ浮为岩体的浮重度。

分别采用以下2种方案计算地面以下不同深度处的竖向应力，并与理论解进行比较。

方案1：在岩基表面定义孔压边界即考虑渗透力作用，不再施加静水压力荷载。

方案2：在岩基表面同时定义孔压边界和静水压力荷载。

表1为按以上2种方案计算所得的竖向应力，由表1可知，方案1在每一个深度处计算结果与理论值相差均为0.15 MPa，恰好是地基表面静水压力值。而方案2计算结果与理论值完全相同，即考虑了式(2)中的面力项(静水压力荷载)后是正确的，说明应用ABAQUS进行渗流应力耦合分析时，定义孔压边界同时必须施加相应的静水压力荷载。另外，本算例岩石是没有左右边界的，不需要考虑静水压力；而对于有边界的物体，在水下，物体的上边界和左右边界均需定义孔压边界同时必须施加相应的水荷载。

表1 不同方案下竖向应力Table 1 Vertical stresses in different load cases

注:正值表示拉应力；负值表示压应力

关于考虑材料渗流作用后，为何还要施加表面水压力荷载，在理论上可从以下2个方面进行说明：

(1)从材料内部结构组成方面来看，不同材料的内部结构组成形式有以下3种情况：颗粒与颗粒的点接触、封闭的实体表面及二者的混合，其截面如图2所示。从宏观荷载角度上来说，对于图2(a)所示材料，水压力将表现为100%体积力的形式；对于图2(b)所示材料，水完全不能渗入，水压力将表现为100%表面水压力的形式；对于图2(c)所示材料，不仅要考虑渗透力，而且有必要考虑作用于材料表面的水压力，而岩体和混凝土恰恰属于这一类材料。需要说明的是，对于图2中的(a)和(c)所示的可渗透性材料而言，利用ABAQUS进行渗流应力耦合分析时，均需在其边界上施加静水压力荷载，并同时定义孔隙水压边界。

图3 坝体断面尺寸Fig.3 Dimensions of the cross section of dam

(2)从研究对象受力平衡方面来看，若以整个岩基(固相骨架和孔隙水)为研究对象，则外荷载包括自身重力和顶部表面水压力荷载；若以固相颗粒骨架为研究对象，根据有效应力原理可知，此时的应力平衡方程还应包含孔隙水压力项。利用ABAQUS进行渗流应力耦合计算时，在边界条件中定义的孔压并非以荷载的形式作用于骨架上，仅是用以提供孔隙水压力项，而表面水压力则是以外荷载的形式单独存在的。

4 混凝土重力坝坝底面扬压力施加问题

某重力坝，横断面见图3，其上游水位为175 m，下游无水，坝体混凝土的重度为24 kN/m3，弹性模量为25 GPa，泊松比为0.167;坝基岩体的重度为27 kN/m3，弹性模量为20 GPa，泊松比为0.25。坝体混凝土渗透系数为1×10-9m/s，初始孔隙比为0.1，假定初始饱和度为0；坝基岩体渗透系数为1×10-7m/s，初始孔隙比为0.1，假定初始饱和度为1。

图4 坝内孔压分布Fig.4 Distribution of pore pressure in the dam

4.1 坝体浸润线的确定

坝体和坝基均按可渗

材料考虑，计算得到的坝体浸润线如图4所示。从图4可以看出，坝体内浸润线明显偏高，与实际工程存在差异，原因是模型中未考虑坝内的排水设施等因素，而本模型旨在研究扬压力的施加问题，因此对下文的分析没有影响。

图5 坝内应力分布Fig.5 Distribution of stress in the dam

4.2 考虑坝基和坝体均透水情况

认为坝基和坝体均为

透水介质，该工况计算得到的应力属于有效应力(见图5)，即考虑到孔隙水压力的影响，而建基面扬压力荷载的实质就是认为坝体不透水时作用于坝体与坝基胶结面处孔隙水压力，因此不必再单独考虑坝基扬压力荷载，故以该情况作为参考基准来研究扬压力的施加问题。

4.3 考虑坝基透水，坝体不透水情况

由于坝体混凝土的渗透性较小，计算中可以将其视为不透水介质来考虑，仅认为坝基透水。为了避免建基面上共用结点的单元类型不同的问题，数值分析时通过定义坝体渗透系数为0来反映坝体的不透水性，故该工况下坝体和坝基的单元类型均为CPE8P(CPE8P是ABAQUS中的一种单元类型，即考虑孔压的8节点平面应变单元)，此时整体模型的渗流场分布如图6所示。

图6 考虑坝基的渗流作用，不考虑坝体渗流时 整体模型的渗流场Fig.6 Seepage field of the whole model in consideration of foundation seepage effect in the absence of dam body seepage

按照目前已有的理论[12-13]，应该考虑的荷载有坝基渗透力、上下游坝面静水压力及建基面扬压力，而在运用ABAQUS进行分析时，考虑坝基的渗流作用，在坝底面存在孔隙水压力，如图7所示。

图7 建基面孔隙水压力Fig.7 Pore water pressure in the dam base surface

针对图7所示孔隙水压力能否以荷载的形式传递到坝上，即此时是否应在坝底面再定义扬压力荷载的问题，分别采用以下4种方案进行分析。

方案1：考虑坝基的渗流作用，不考虑坝体渗流，也不施加坝底面扬压力。坝体内无孔隙水压作用，探讨该软件是否可自动施加坝底面的扬压力作用。

方案2：考虑坝基的渗流作用，不考虑坝体渗流，但在建基面施加扬压力。扬压力按照图7所示面力分布进行施加。

方案3：考虑坝基的渗流作用，不考虑坝体渗流，但坝体自身浸润线以下按照浮重度考虑[14]。该情况虽考虑坝体内有效应力，但不能考虑渗透力作用。

方案4：不考虑坝体和坝基的渗流作用，但在建基面施加扬压力，坝基和坝体内均无孔隙水压，相当于按总应力计算。

在以上4种分析方案中，均考虑了重力，坝体上、下游面水压力，坝基上下游表面库水压力荷载。

采用以上4种方案及考虑坝体渗流作用情况计算所得坝底面竖向应力如图8所示，坝体内部的竖向应力分布如图9所示。方案1、2和4均未考虑坝体内的孔隙水压力，因此计算得到的坝体应力为总应力，其中，方案4是不考虑坝体坝基渗透作用时，进行重力坝静力分析时所普遍采用的方法，故以方案4作为判断依据，来比较仅考虑坝基渗流作用后，ABAQUS是否可自动施加坝建基面的扬压力荷载。

图8 坝底面竖向应力Fig.8 Vertical stress in the bottom of dam

图9 坝内竖向应力分布等值线Fig.9 Distribution of vertical stress in the dam

由图8、图9知：方案1计算得到的坝内和坝底面竖向应力与方案4差异较大；方案2计算得到的坝内和坝底面的竖向应力与方案4基本相同。其原因在于，方案2不仅考虑了建基面上的孔隙水压力，而且在坝底面施加了扬压力荷载。这说明定义坝基的渗流作用后，虽然建基面上存在孔隙水压力，但是其并不能以外荷载的形式传递到坝底面，所以在利用ABAQUS进行重力坝渗流应力耦合分析时，仅考虑坝基透水，而不考虑坝体透水时，应单独定义坝体扬压力作用。

下面将针对扬压力的2种不同施加方式进行探讨，比较图9(b)和图9(a)发现，通过在坝底面施加面力来定义扬压力时，仅对坝底部分区域的应力有影响，且造成坝踵处应力集中；比较图9(c)和图9(a)发现，通过对浸润线以下坝体按浮重度考虑来定义扬压力时，虽然坝踵处也会出现应力集中，但该工况下坝内和坝底面的竖向应力与考虑坝体渗流作用时(图5)所得结果基本一致。

因此，笔者认为在利用该软件对重力坝进行渗流应力耦合分析时，若考虑坝基渗透性而不考虑坝体渗透性时应对坝基施加扬压力，但扬压力不能以面力的形式施加，宜采用浸润线以下坝体按浮重度进行计算的方式；若同时考虑坝体和坝基的渗透性时则无需再对建基面施加扬压力。

5 结 论

本文旨在结合算例探讨应用ABAQUS进行渗流应力耦合分析时模型中水荷载及边界条件施加的问题，分析结果表明：

(1)利用ABAQUS进行渗流应力耦合分析时，应在坝体和坝基表面同时定义水荷载和孔压边界条件，定义的表面水压力将以外荷载的形式作用于固相颗粒骨架上，而定义孔压边界的作用则是用以定义孔隙水压力项，以及满足渗流连续方程。而并不像有些学者认为的只要定义渗流孔压边界即考虑了渗透力的作用，从而无需再施加静水压力。

(2)不考虑坝体渗透性时，在建基面虽然存在孔隙水压力，但分析发现其并不能以荷载的形式传递到坝体上，此时需重新定义扬压力荷载，但不能以面力的形式施加，宜采用浸润线以下坝体按浮重度进行计算的方式考虑；而考虑坝体渗透性时，坝内存在孔隙水压力作用，故此时无需再单独定义扬压力。

参考文献：

[1] 杨天鸿, 唐春安, 梁正召, 等. 脆性岩石破裂过程损伤与渗流耦合数值模型研究[J]. 力学学报, 2003, 35(5): 533-541.

[2] 梁 通, 金 峰. 基于广义有效应力原理的混凝土坝分析[J]. 水力发电学报, 2007, 28(2): 47-51.

[3] 柴军瑞, 仵彦卿. 碾压混凝土坝渗流场与应力场耦合分析的数学模型[J]. 水利学报, 2000, (9): 33-35.

[4] 沈振中, 陈小虎, 徐力群. 重力坝应力-渗流相互作用的无单元耦合分析[J]. 岩土力学, 2008, 29(增): 74-79.

[5] HIBBETT D, KARLSSON B, SORENSEN P. ABAQUS/Standard: User’s Manual[K]. U.S.: Hibbitt, Karlsson & Sorenson, 1998.

[6] 陈卫忠, 伍国军, 贾善坡. ABAQUS在隧道及地下工程中的应用[M]. 北京: 中国水利水电出版社, 2010.

[7] 庄 茁, 由小川, 廖剑晖, 等.基于ABAQUS的有限元分析和应用[M]. 北京: 清华大学出版社, 2009.

[8] 常晓林, 袁 曦. 基于流-固耦合的强度折减法在重力坝抗滑稳定分析中的应用[J]. 武汉大学学报(工学版), 2012, 45(5): 545-550.

[9] 柴军瑞. 混凝土坝水荷载讨论[J]. 水电能源科学, 2000, 18(2): 18-21.

[10] 张晓咏, 戴自航. 应用ABAQUS程序进行渗流作用下边坡稳定分析[J]. 岩石力学与工程学报, 2010, 29(增): 292-294.

[11] 林凯生, 李宗利. 高渗透孔隙水压作用下混凝土损伤破坏过程数值分析[D]. 杨凌: 西北农林科技大学, 2010.

[12] 王 媛, 速宝玉. 坝基应力计算中的水荷载组合形式[J]. 勘察科学技术, 1995, (3): 3-7.

[13] 黄耀英, 沈振中, 田 斌. 地基水荷载对混凝土坝位移影响研究[J]. 水利水运工程学报, 2010, (1): 42-49.

[14] 范书立, 陈建云, 林 皋. 渗透压力对重力坝有限元分析的影响研究[J]. 岩土力学, 2007, 28(增): 575-581.