ABAQUS在隧道及地下工程中的应用 - 万水书苑-出版资源网

知识要点：



深部岩体工程特点



ABAQUS 数值模拟功能



地下储气库的长期稳定性分析



深部引水隧洞的稳定性分析

本章导读：

本章主要研究 ABAQUS 在深部岩体工程中的应用。首先介绍 ABAQUS 在此方面的模拟

功能，然后分别以地下高应力储气库和深埋引水隧洞为工程应用对象，详细阐述 ABAQUS 在

这些方面的应用全过程，对研究类似深部岩体工程的应用者提供较好的学习和借鉴作用。

8.1 深部岩体工程简介

能源和矿产资源制约着国民经济的发展。随着浅部资源的逐渐减少和枯竭，矿物资源地

下开采的深度越来越大。人类的工程活动已经深入到地下 4000m 以下的深度，如逾千米乃

至数千米的矿山、水电工程埋深逾千米的引水隧道、核废料的深层地质处理问题、深地下防

护工程等。在我国西部水电能源开发、南水北调、西气东送等工程中普遍涉及到深部岩体地

质问题。如我国西部已建或将建的拉西瓦、二滩、天生桥等水电工程都遇到过高地应力、高

水压等问题，甚至在浅表的大坝建基面的开挖时，都会发生由于高地应力释放产生的表面岩

体成层剥落现象。西部地区的特殊地质、地形条件使得西部水电资源开发中地下深埋工程量

很大。

深部岩体所处的高地应力、高地温和高渗透水压力的特殊环境，伴随着深部岩体的力

学响应明显有别于浅部岩体的力学行为。由于深部岩体受到各种载荷作用、岩体介质作用

本身的复杂性、认知的不确定性和表现出一系列新的岩体地质力学特征，致使深部岩体的

地质结构特征、裂隙渗透特性、变形机理、破坏规模和强度特征以及高渗透水压下的流变

机理难以用传统的理论加以合理的解释，引起了国内外学者的极大关注，成为近几年来该

领域的研究热点。

由于深部岩体所处的复杂的地质环境，导致地应力高、温度高、渗透压高，加之较强的

时间效应，使深部岩体的组织结构、基本行为特征和工程响应均发生根本性变化。深部岩体表

现出明显的非线性力学特性。进入深部的工程岩体所属的力学系统不再是浅部工程岩体所属的

线性力学系统，而是非线性力学系统，传统理论、方法与技术已经部分或全部失效。另一方面，

深部岩体处于多场、多相耦合作用，地下水、瓦斯、温度均会对岩体的基本性质和工程响应带

来很大影响。因此，深部岩体非线性损伤变形特性、多场耦合作用机理等方面的研究，将为深

部岩体工程灾害的预测和控制提供理论基础，并将从新的角度和新的高度认识深部岩体的特殊

性质，以指导工程实践，创造经济和社会效益。

8.2 ABAQUS 数值模拟功能

8.2.1 损伤

ABAQUS 可 通 过 多 种 方 式 定 义 损伤 ， 诸 如 混凝 土损 伤 塑 性 模 型 、采 用 用 户子 程 序

（USDFLD、UFIELD、UMAT 等）。

混凝土损伤塑性模型适用于 ABAQUS/Standard 和 ABAQUS/Explicit 两个模块。其功能

包括：



能够模拟各种结构类型（梁、桁架、壳和实体）中混凝土和其他准脆性材料。



采用各向同性弹性损伤结合各向同性拉伸和压缩塑性理论来表征混凝土的非弹性行为。



主要用于钢筋混凝土结构分析，也能够用于素混凝土结构分析。



可采用加强筋模拟混凝土中的钢筋。



可以模拟低围压时，混凝土受到单调、循环或动载作用下的力学行为。



结合非关联多重硬化塑性和各向同性弹性损伤理论来表征材料断裂过程中发生的不

可逆损伤行为。



周期载荷反向作用时，可以控制材料的刚度恢复。



可以定义与应变率相关的性状。



在 ABAQUS/Standard 中可以结合使用粘塑性正规本构方程来提高软化区域的收敛性。



材料的弹性行为应为各向同性和线性的。

USDFLD、UFIELD 的功能有相似之处，均可以在子程序中定义场变量并输出，同时可以

定义与场变量相关的变量，诸如损伤变量等。

以 USDFLD 为例，其子程序的简要形式如下：

SUBROUTINE USDFLD(FIELD,STATEV,PNEWDT,DIRECT,T,CELENT, 1 TIME,DTIME,CMNAME,ORNAME,NFIELD,NSTATV,NOEL,NPT,LAYER, 2 KSPT,KSTEP,KINC,NDI,NSHR,COORD,JMAC,JMATYP,MATLAYO, 3 LACCFLA) C INCLUDE 'ABA_PARAM.INC' C CHARACTER*80 CMNAME,ORNAME CHARACTER*3 FLGRAY(15) DIMENSION FIELD(NFIELD),STATEV(NSTATV),DIRECT(3,3), 1 T(3,3),TIME(2) DIMENSION ARRAY(15),JARRAY(15),JMAC(*),JMATYP(*),COORD(*) CCCCCCCCCCCCCCCCCCC 定义场变量并赋予 FIELD(NFIELD) 定义与场变量相关的损伤变量,赋给状态变量 STATEV(NSTATV)

CCCCCCCCCCCCCCCCCCC RETURN END

与子程序相关的 INP 命令流格式为采用 DEPENDENCIES 语句，如：

*Elastic,DEPENDENCIES=n **n 为场变量个数 25000 , 0.25, , 0 24950.04997 , 0.25, , 0.00001 24504.96683 , 0.25, , 0.0001 20468.26883 , 0.25, , 0.001 3383.382081 , 0.25, , 0.01 1.134998244 , 0.25, , 0.05

这里定义了等效蠕变应变 CEEQ 作为场变量，场变量 CEEQ 与弹性模量 E 的关系为：

0 0.05

e

ceeq

E

E

 



（8-1）

式中，

E 为初始弹性模量。

₀

接下来采用如下语句定义状态变量的个数：

*USER DEFINED FIELD *DEPVAR

n **状态变量的个数

UMAT 是用户材料子程序（User-defined Material Mechanical Behavior，简称 UMAT）

，是

ABAQUS 提供给用户定义自己的材料属性的 FORTRAN 程序接口，使用户能使用 ABAQUS

材料库中没有定义的材料模型。UMAT 子程序具有强大的功能，不仅能定义材料的本构关系

和屈服准则，还能实现材料本构关系用于力学分析的任何过程，当然，UMAT 中必须提供材

料本构的雅可比（Jacobian）矩阵而实现增量计算。有关 UMAT 的具体介绍详见第 10 章。

8.2.2 非线性蠕变岩体本构关系

ABAQUS 有限元分析中可以定义粘弹、粘塑、粘弹塑性等各种与粘性蠕变相关的本构关

系。其中定义粘塑、粘弹塑性本构关系可采用*CREEP 或*DRUCKER PRAGER CREEP 来定义。

下面以*DRUCKER PRAGER CREEP 实现粘弹塑性本构关系为例加以介绍。

岩土材料在某些条件下将发生蠕变。尤其当蠕变时间尺度和加载率处于同一数量级时需

考虑蠕变与塑性的耦合，此时 ABAQUS 提供了一个蠕变模型来定义材料的“经典”蠕变行为，

其塑性变形服从扩展 Drucker-Prager 模型。材料的蠕变与塑性变形密切相关（此处可通过定义

蠕变流动势和试验数据），因此在材料属性中必须包含塑性及塑性硬化的定义。

1．蠕变应力的定义

采用 Drucker-Prager 屈服准则所用的蠕变势函数为双曲线函数，如果在 ABAQUS/Standard

中定义了蠕变，线性 Druker-Prager 塑性模型也采用双曲线塑性流动势。

假定存在应力点的蠕变等倾面，该等倾面具有相同的蠕变“强度”，并由等效蠕变应力确

定。当材料发生塑性变形时，需要等效蠕变面与屈服面一致，因此等效蠕变面可由屈服面等比

例缩小得到。在 p-q 平面上

，

蠕变面与屈服面平行，如图 8-1 所示。

图 8-1 p-q 面屈服面和双曲线流动势函数示意图

在定义材料蠕变特性时，等效蠕变应力

_

cr

_{可表示为：}

(

tan )

1

1

tan

3

(

tan )

1

1

tan

3

(

tan )

cr

q

p

q

p

q

p























_



_

































_



2．蠕变流动

蠕变应变率采用与塑性应变率相同的双曲线流动势函数。即：

2 2 0

(

tan )

tan

cr

G











q



p



（8-2）

式中：

. 0 pl _0, pl ₀  





 



为初始屈服应力；为偏移率，用来定义双曲线趋向于渐近线的

速率（当偏移率趋向于零时，蠕变势函数趋向于直线）；



为 p-q 平面上的剪涨角。

蠕变势函数为光滑连续的曲线，能保证蠕变流动方向唯一。在高围压作用下，该函数趋

近于线性 Druker-Prager 流动势，同时与静水压力轴正交。

当材料为非相关流动时，刚度矩阵为非对称阵，此时需要采用非对称矩阵存储和非对称

计算方法。如果摩擦角和剪涨角相差不大，且模型中非弹性变形区域受到限制时，材料刚度近

似对称矩阵会给出可以接受的收敛速度，从而不需要非对称矩阵算法。

3．蠕变定律

在 ABAQUS/Standard 中，蠕变通过等效的“单轴蠕变”定义，称为“蠕变定律”。ABAQUS

提供了三种可供用户输入相关参数来定义的蠕变定律：时间硬化蠕变定律、应变硬化蠕变定律

和 Singh-Mitchell 蠕变定律。

对于很多的实际工程问题，由于材料蠕变的函数形式往往较为复杂，因此需要采用用户

编制子程序 CREEP 来实现，采用 CREEP 子程序自行编制非线性蠕变本构模型程序与 ABAQUS

当通过单轴压缩屈服应力



_c

定义蠕变时

当通过单轴拉伸屈服应力



_t

定义蠕变时

当通过粘聚力定义蠕变时

结合计算，从而解决工程问题。

ABAQUS 可在 CREEP 子程序中定义非线性蠕变本构模型，

其参与整个计算的具体流程图

如图 8-2 所示。

ABAQUS 应变更新

( ) ( ) ( )

total total total

n n n t t t t     ABAQUS产生新的总应变增量 1 ( ) total n t  _  ABAQUS第n+1步的平衡迭代 收敛 不收敛 n=n+1，返 回第一步 减小时间 步长 ABAQUS 在平衡时间tn 时： t, , 非线性蠕变本构模型 是否进入D-P屈服 Y N ( ) ( ) ( ) e total c n n n t t t         e( ) total( ) c( ) P( ) n n n n t t t t          应力更新（塑形要迭代求解） (tn t) ( )tn ( )tn       图 8-2 CREEP 子程序计算流程框图

用户子程序 CREEP 的模板如下：

********************************************************** SUBROUTINE CREEP(DECRA,DESWA,STATEV,SERD,EC,ESW,P,QTILD, 1 TEMP,DTEMP,PREDEF,DPRED,TIME,DTIME,CMNAME,LEXIMP,LEND, 2 COORDS,NSTATV,NOEL,NPT,LAYER,KSPT,KSTEP,KINC) C INCLUDE 'ABA_PARAM.INC' C CHARACTER*80 CMNAME C DIMENSION DECRA(5),DESWA(5),STATEV(*),PREDEF(*),DPRED(*), 1 TIME(2),EC(2),ESW(2),COORDS(*) 定义 DECRA **注:最常用的是定义 DECRA(1)､DECRA(5) 当考虑材料膨胀属性时,定义 DESWA RETURN

END

**********************************************************

与用户子程序 CREEP 相关的命令流为：

*DRUCKER PRAGER CREEP,LAW=USER

8.2.3 考虑渗流的岩体本构关系

ABAQUS 能够求解多孔介质的饱和渗流、非饱和渗流及二者的混合问题（诸如坝体渗流

自由面）的计算。计算过程中可以考虑渗流重力的作用，并能够求解流体的总孔隙压力或超孔

隙压力，渗透定律可采用达西定律或更广泛的非线性定律。在流体重力载荷不可忽略或瞬态毛

细吸力较明显，即“湿化作用”不可忽略的问题中，需要求解总孔隙压力。

需要指出的是：ABAQUS 在孔隙压力解的形式上有两种，分别为总孔隙压力和超孔隙压

力，两者仅仅是重力载荷上占主导地位的情况下才有明显区别。当采用 GRAV 分布载荷类型

来定义模型的重力载荷时，ABAQUS 程序提供总孔隙压力解；其他情况（如采用分布载荷类

型 BX、BY、BZ 定义重力载荷时）提供超孔隙压力解。下面分别举例加以说明。

（1）采用 GRAV 时。

******************************** *Material, name=ROCK *Density 2500, *Elastic 1.e10, 0.3 *Permeability, specific=1.0e4 1.e-8,1.0

*Initial conditions, type=stress, geostatic Eall, 0,10, -2.0e5,0, 0.85,0.85

*Initial conditions, type=ratio Eall, 1.0

*Initial conditions, type=pore pressure Eall, 0,10, 1.0e5,0 ………. *step *Geostatic …. *Dload Eall, grav, 10, 0,-1.0 …. *end step ********************************

此处计算得到的孔压为总孔隙压力，分别须计算有效应力平衡和孔压平衡。其中，有效

应力的平衡是采用土力学公式（根据三相物理指标求的其他指标）求得有效地应力。

（2）采用 BY 或 BZ 时。

******************************** *Material, name=ROCK

*Density 2500, *Elastic 1.e10, 0.3 *Permeability, specific=1.0e4 1.e-8,1.0

*Initial conditions, type=stress, geostatic Eall, 0,10, -2.0e5,0, 0.85,0.85

*Initial conditions, type=ratio Eall, 1.0 ………. *step *Geostatic …. *Dload

Eall, BY, -3.0e4 …. *end step ********************************

此处计算得到的孔压为超孔隙压力，不考虑重力载荷的影响。ABAQUS 系统在自动平衡

地应力后，会自动计算出初始平衡孔压力，此时不需要定义初始孔压。

对于裂隙岩体固液耦合的岩体本构关系，用户可以通过定义岩体损伤、孔隙度或渗透系数

变化等实现本构关系的改变。可通过用户子程序 USDFLD 来实现，其子程序的简要形式如下：

********************************************************** SUBROUTINE USDFLD(FIELD,STATEV,PNEWDT,DIRECT,T,CELENT, 1 TIME,DTIME,CMNAME,ORNAME,NFIELD,NSTATV,NOEL,NPT,LAYER, 2 KSPT,KSTEP,KINC,NDI,NSHR,COORD,JMAC,JMATYP,MATLAYO, 3 LACCFLA) C INCLUDE 'ABA_PARAM.INC' C CHARACTER*80 CMNAME,ORNAME CHARACTER*3 FLGRAY(15) DIMENSION FIELD(NFIELD),STATEV(NSTATV),DIRECT(3,3), 1 T(3,3),TIME(2) DIMENSION ARRAY(15),JARRAY(15),JMAC(*),JMATYP(*),COORD(*) CCCCCCCCCCCCCCCCCCC 定义场变量并赋予 FIELD(NFIELD) 定义与场变量相关的损伤变量,赋给状态变量 STATEV(NSTATV) 定义渗透系数的变化量,并赋给状态变量 STATEV(NSTATV) CCCCCCCCCCCCCCCCCCC RETURN END **********************************************************

8.3 储气库的计算分析

能源是一个国家的经济命脉，一旦能源发生危机，将引起社会的严重动荡，破坏经济的

发展。在世界 GDP 前几位的国家中，中国是较晚建立国家战略储备和民间商业储备的国家之

一。地下能源储存一般放置在盐岩、非渗透性岩层以及多孔隙岩层中。目前世界上已有的地下

储气库类型主要包括枯竭油气藏、含水构造地下储气库、盐岩地下储气库、废弃矿井和水密封

岩石洞室。而盐岩具有非常低的渗透特性（渗透率小于 10

-20

m

2

）与良好的蠕变行为，能够适

应储存压力的变化；其力学性能较为稳定（损伤与损伤自我恢复）

，能够保证储存硐库的密闭

性（J.E.Quintanilha）。因此，国际上公认盐岩体是能源（石油、天然气）储存的最理想的介质。

目前全世界各地大约有五百多座地下储气库，其中有 44 座是利用含盐岩层。

就盐岩地下储库而言，目前都是根据储量和服务年限的要求，通过水溶开采形成一定规模的

溶腔群。但在某一地区形成 40～50 万 m

3

_{规模的溶腔群一般需要 5～10 年的周期，而且还需要与}

之配套的卤水处理化工厂，否则生产的卤水将对环境产生很大的污染。目前，国内外利用盐矿开

采形成的废弃溶腔作为调峰储气库尚无成功先例，利用废弃盐岩溶腔作为储气库具有建库时间短、

成本低等优点，但必须对废弃盐穴的长期变形、储库封孔的套管鞋高度等关键技术参数进行评估。

本节应用大型有限元软件 ABAQUS 卓越的非线性功能，将实验室蠕变试验所得的蠕变本

构模型植入 ABAQUS 软件，并对某废弃溶腔群储库的蠕变变形和工作性态开展数值仿真，研

究方法及成果对国内应用废弃溶腔作为储气库时的工作压力及储库压力设计具有借鉴作用。

8.3.1 流变本构模型

国内外众多学者的实验室试验结果表明：盐岩在偏应力作用下易产生蠕变，且与时间呈

高度非线性关系。盐岩的长期蠕变特性比较复杂，一般认为：盐岩的蠕变速率是作用在其上的

偏应力和温度的高阶非线性函数。盐岩在三轴压力状态下，一般经历瞬态蠕变、稳定蠕变和加

速蠕变三个阶段。由于盐岩的瞬态蠕变时间很短，因此在储气库的稳定性和长期变形时，主要

研究盐岩的稳态蠕变特性。Carter（1993）和 Chan（1997）曾提出稳态蠕变率的统一表达式：

3 1 3

(

) (

)

( )

cr

f

_c

D

H T















（8-3）

式中，

f

_c

(



₃

)

为围压的影响函数，

D

(



₁





₃

)

为偏差应力的影响函数， ( )

H T 为温度函数。

杨春和等人通过对某矿区盐岩及泥岩试样开展常温下的三轴蠕变实验研究，

得到盐岩的典

型蠕变曲线如图 8-3 所示，并提出了盐岩和泥盐岩岩的稳态蠕变率表达式：

1 3

(

)

cr

_A

n













（8-4）

式中，

_

_{ 蠕变应变率，A，n 为盐岩的材料特性参数。}

cr

8.3.2 储气库的长期稳定性分析

为科学合理地评价废弃溶腔作为储气库的可行性，必须对废弃溶腔已有的蠕变变形进行

合理的预估，并对储库储气后的盐岩蠕变变形机理及废弃溶腔群矿柱的稳定性进行数值模拟。

本节通过数值模拟研究重点讨论储库在过去 10 年中所发生的蠕变变形以及未来 10 年的储气服

务期内因工作状态变化而产生的蠕变变形和破坏范围。

0.00

0.02

0.04

0.06

0.08

0.10

0

10

20

30

40

50

60

时间/h

蠕

变

％

图 8-3 盐岩的三轴蠕变试验曲线

1．工程概况

我国某盐岩矿区的盐岩开采已具有二三十年的时间，在当地已形成一定规模的地下溶

腔群。该矿区地层平缓，构造简单，盐层分布范围大，达 60.5km

2

_{，且分布稳定，厚度大，}

一般为 100m 以上，盐层的含盐率高，夹层少且厚度较小，直接顶底板均为含钙芒硝含膏

泥岩或致密泥岩，抗压强度大，封闭性好。但该矿区水洗形成的溶腔大多数体积太小而无

大的利用价值。根据现场的声纳探测结果溶腔基本呈梨状分布，其中的废弃溶腔群依次按

1、2、3 和 4 号井命名，埋深为 900～1000m，腔体体积均为 1.5



10

5

m

3

左右，符合天然气

储存的基本要求，但储库的间距较小，一般在 20m 左右。储库区域的溶腔形态及岩层的分

布如图 8-4 所示。

1－泥岩层，2－盐层，3－泥质夹层 图 8-4 废弃溶腔岩层分布

2．数值模型

根据废弃溶腔的分布特征和埋深，数值分析模型取距离地表 525m，深 900m 的岩体进行

分析。根据岩层分布特征和溶腔形态所建立的有限元计算模型的水平切面和垂直剖面如图 8-5

所示。

（a）溶腔 1、2 的垂直纵剖面 （b）溶腔 3、4 的垂直纵剖面 （c）4 个溶腔的水平切面 图 8-5 有限元分析网格

采用 ABAQUS 数值模拟废弃溶腔作为储气库长期运行的稳定性，ABAQUS 在这方面的突

出优点是：具有杰出的非线性分析能力和广泛的模拟性能，并拥有大量不同种类的单元类型、

材料本构模型和载荷形式。计算中，盐岩体和泥岩采用 Drucker-Prager 模型，盐岩的稳态蠕变

应用上节所建议的稳态流变本构形式，并将本构模型以 FORTRAN 语言编写出子程序和

ABAQUS 主程序连接计算。

计算模型分为两部分：水平切面和纵向剖面，文件名分别为 QIE1.INP 和 QIE5.INP。

QIE1.INP：

********************************************** **节点命令流 *node,nset=nall 1, 0.000000, 840.000000 2, 200.000000, 840.000000 3, 30.310307, 840.000000 4, 57.686607, 840.000000 5, 82.412910, 840.000000 6, 104.745735, 840.000000 7, 124.916771, 840.000000 8, 143.135269, 840.000000 9, 159.590225, 840.000000 10, 174.452377, 840.000000 11, 187.875885, 840.000000 12, 0.000000, 440.000000 13, 0.000000, 809.567566 14, 0.000000, 780.567932 15, 0.000000, 752.933716 16, 0.000000, 726.600464 17, 0.000000, 701.507080 18, 0.000000, 677.595093 19, 0.000000, 654.808899 1 号井 2 号井 3 号井 4 号井 1 号井 2 号井 4 号井 3 号井

20, 0.000000, 633.095581 ………. 5944, 217.853912, 356.961243 **单元命令流 *element,type=cpe4r,elset=eall 1, 2, 11, 61, 42 2, 11, 10, 80, 61 3, 10, 9, 99, 80 4, 9, 8, 118, 99 5, 8, 7, 137, 118 6, 7, 6, 156, 137 7, 6, 5, 175, 156 8, 5, 4, 194, 175 9, 4, 3, 213, 194 10, 3, 1, 13, 213 11, 42, 61, 62, 43 12, 61, 80, 81, 62 13, 80, 99, 100, 81 14, 99, 118, 119, 100 15, 118, 137, 138, 119 16, 137, 156, 157, 138 17, 156, 175, 176, 157 18, 175, 194, 195, 176 19, 194, 213, 214, 195 20, 213, 13, 14, 214 21, 43, 62, 63, 44 22, 62, 81, 82, 63 23, 81, 100, 101, 82 24, 100, 119, 120, 101 25, 119, 138, 139, 120 26, 138, 157, 158, 139 5829, 5537, 5530, 5723, 5937 **********************************************

QIE5.INP：

********************************************** **节点命令流 1, 0.000000, 600.000000 2, 120.000000, 480.000000 3, 11.401139, 588.598877 4, 22.126673, 577.873352 5, 32.216637, 567.783386 6, 41.708691, 558.291321 7, 50.638268, 549.361755 8, 59.038700, 540.961304 9, 66.941338, 533.058655

10, 74.375687, 525.624329 11, 81.369492, 518.630493 12, 87.948860, 512.051147 13, 94.138351, 505.861633 14, 99.961067, 500.038940 15, 105.438744, 494.561249 16, 110.591827, 489.408173 17, 115.439545, 484.560455 18, 120.000000, 150.000000 19, 120.000000, 471.750000 20, 120.000000, 463.500000 ………. ………. 9883, 267.019470, 235.302170 **单元命令流 1, 1, 74, 113, 3 2, 74, 75, 128, 113 3, 75, 76, 143, 128 4, 76, 77, 158, 143 5, 77, 78, 173, 158 6, 78, 79, 188, 173 7, 79, 80, 203, 188 8, 80, 81, 218, 203 9, 81, 82, 233, 218 10, 82, 83, 248, 233 11, 83, 84, 263, 248 12, 84, 85, 278, 263 13, 85, 86, 293, 278 14, 86, 87, 308, 293 15, 87, 88, 323, 308 16, 88, 89, 338, 323 17, 89, 90, 353, 338 18, 90, 91, 368, 353 19, 91, 92, 383, 368 20, 92, 93, 398, 383 21, 93, 94, 413, 398 22, 94, 95, 428, 413 23, 95, 96, 443, 428 24, 96, 97, 458, 443 25, 97, 98, 473, 458 ………. ………. 9803, 9188, 9065, 9098, 9096 **********************************************

3．初始地应力、岩体基本力学参数和溶腔工作压力

数值分析计算范围内的岩体以泥岩、泥质夹层及盐岩层为主，根据试验结果，本次计算

所采用的岩体材料基本力学特性如表 8-1 所示。废弃溶腔的埋深在 1000m 左右，计算中假定

盐岩初始条件下处于静水压力状态。数值模拟盐矿在初始水溶法开采时，溶腔内壁的工作压力

为饱和盐水作用下的静水压力，取值 12MPa，在模拟储气阶段时，首先考虑内腔压力在 3 个

月内由 12MPa 升至 14MPa，之后在 3 个月内降至 7MPa，并在这样的工作压力下工作 9 年半，

储库的工作压力历时曲线如图 8-6 所示。根据实验室的试验结果，盐岩及含盐泥岩的流变参数

如下：

盐岩：

A=1.018×10

-8

，n=3

含盐泥岩：

A=1.018×10

-8

，n=2

表 8-1 废弃盐矿岩体基本力学参数

岩层 弹性模量 E（GPa） 泊松比 粘聚力 c（MPa） 摩擦角 抗拉强度（MPa）

泥岩 10 0.27 1.0 35o 1 盐岩 18 0.30 1.0 30o ₁ 夹层 4 0.30 0.5 30o 0.5 0 2 4 6 8 10 12 14 16

0

5

10

15

20

时间(y)

压

力

(

M

P

a

)

图 8-6 盐岩溶腔运行压力示意图

初始地应力平衡：由于模型上表面是水平的，初始地应力采用*initial conditions,type=

stress,geostatic 施加。

岩体材料力学参数：盐岩和泥岩分别定义蠕变材料特性，采用 Drucker-prager creep 定义，

夹层不定义蠕变特性。相关命令流如下：

夹层：

********************************************** *Material, name=middle *Elastic 4000., 0.3 *Density 0.0023, （10.25, 14） （10.5, 7）

*Drucker Prager 39.76, 1., 39.76 *Drucker Prager Hardening 1.0 **********************************************

泥岩：

********************************************** *Material, name=mud *Elastic, DEPENDENCIES=1 10000., 0.27,,0 10000., 0.27,,1 *Density 0.0023, *Drucker Prager 43.3, 1., 43.3

*Drucker Prager Hardening 1.96,0

*drucker prager creep,law=user *DEPVAR

2

*USER DEFINED FIELD

**********************************************

盐岩：

********************************************** *Material, name=salt *Density 0.0023, *Elastic, DEPENDENCIES=1 18000., 0.3,, 0,0 18000., 0.3,, 0.1 *Drucker Prager 39.76, 1., 39.76 *Drucker Prager Hardening 1.99,0.

*drucker prager creep,law=user *DEPVAR

2

*USER DEFINED FIELD

**********************************************

溶腔工作压力：溶腔运行压力通过*amplitude 实现。命令流：

********************************************** *amplitude,name=cc,time=total time 0, 0, 1,1.2, 2,1.2, 86402,1.2, 88562,1.4, 90722,.7, 95042,.7, 103682,.7,

112322,.7, 129602,.7,172802,.7 **********************************************

4．分析步骤

首先建立初始地应力平衡，然后水溶溶腔，接下来同采同注，给腔体内壁气压力。整个

分析过程的命令流如下：

Step 1

初始地应力平衡。

** STEP：Step-1 *Step, name=Step-1 *Geostatic 1.,1. ** ** LOADS ** ** Name：Load-1 Type：Gravity *Dload EALL, GRAV, 10., 0., -1. *dsload upsurf,p,12.42 ** ** OUTPUT REQUESTS ** *Output, field *Node Output CF, RF, U *Element Output

E, EE, FV, LE, NFORC, PE, PEEQ, PEQC, S, SDV,CE,ceeq *Output, history, variable=PRESELECT

*El Print, freq=999999 *Node Print, freq=999999 *End Step

Step 2

腔体开挖，内压 12MPa。

** STEP：Step-2 *Step, name=Step-2 *Static 1., 1., 1e-05, 1. *model change,remove exca *Dsload,amplitude=cc right, p, 10. left, p, 10. *End Step

Step 3

废弃盐穴 10 年。

*STEP,INC=3000,name=step-3 creep-1：10years *VISCO,CETOL=5E-6

0.1,311040000,1E-8 *END STEP

Step 4

开始储气 3 个月，压力由 12MPa 升至 14MPa。

*STEP,INC=3000,name=step-4 creep-2：3months

*VISCO,CETOL=5E-6 0.1,7776000,1E-8 *END STEP

Step 5

储气库压力由 14MPa 降为 7MPa。

*STEP,INC=3000,name=step-5 creep-3：3months *VISCO,CETOL=5E-6 0.1,7776000,1E-8 *END STEP

Step 6

储气库压力持续不变，运行 9 年半。

*STEP,INC=3000,name=step-6 creep-4：6months *VISCO,CETOL=5E-6 0.1,15552000,1E-8 *END STEP ** --- *STEP,INC=3000,name=step-7 creep-5：1year *VISCO,CETOL=5E-6 0.1,31104000,1E-8 *END STEP ** --- *STEP,INC=3000,name=step-8 creep-6：1year *VISCO,CETOL=5E-6 0.1,31104000,1E-8 *END STEP ** --- *STEP,INC=3000,name=step-9 creep-7：2years *VISCO,CETOL=5E-6 0.1,62208000,1E-8 *END STEP ** --- *STEP,INC=3000,name=step-10 creep-8：5years *VISCO,CETOL=5E-6 0.1,155520000,1E-8 *END STEP **************************************************

5．计算结果分析

（1）腔体蠕变变形分析。

数值计算结果表明，在溶腔形成时，溶腔围岩均有较大的变形，各溶腔围岩最大位移均

发生在腔顶周围和腔底处。以 1 号腔为例，最大位移在腔顶侧边，为 0.134m，腔顶处位移为

0.095m。在溶腔报废的 10 年间，由于内腔有恒定的 12MPa 内压作用，因此溶腔的蠕变变形不

是太大，1 号腔的最大合位移仍在腔顶侧边，为 0.140m，腔顶处合位移为 0.099m。储气后，

由于储库工作压力的波动，溶腔在经历 10 年储气后，储库盐岩的变形蠕变很大，其最大位移

转变到溶腔之间的岩柱，这是因为在溶腔之间岩柱的偏应力较大，因此该部位的蠕变变形也较

大，其中 1 号腔最大合位移为 0.803m，腔顶处移为 0.280m。表 8-2 为四个溶腔在不同时期腔

顶的位移值。

表 8-2 四溶腔腔顶位移变化表 不同时期溶腔顶位移/m 井号 溶腔期间 10 年废弃期 10 年储气期 1 号井 0.095 0.099 0.280 2 号井 0.105 0.111 0.302 3 号井 0.090 0.096 0.281 4 号井 0.083 0.088 0.272

图 8-7 至图 8-9 为四个溶腔在不同时期的位移矢量水平切面图。结果表明：储库盐岩在经

历 10 年的储气周期后，洞壁发生偏向溶腔内部的变形很大，以 1 号腔和 4 号腔之间的围岩变

形尤为明显。

图 8-7 溶腔完时位移矢量图

图 8-8 废弃 10 年后位移矢量图 图 8-9 溶腔储气 10 年时位移矢量图

（2）应力分析。

溶腔过程中围岩四周产生应力集中，在盐岩溶腔储气过程中，由于溶腔内压力发生较大

的变化，围岩的应力也相应产生变化。从溶腔开始到废弃第 10 年，由于溶腔内有饱和盐水的

静水压力作用，且压力保持在 12MPa 不变，溶腔围岩主应力变化不大。但在溶腔储气压力发

生变化后，各溶腔围岩应力均发生了较大的变化。以 1 号腔腔顶为例，溶腔完成时最大主应力

为-17.43MPa，废弃第 10 年时为-17.18 MPa，储气 10 年后为-11.91 MPa，如表 8-3 所示。

从整个流变的过程来看，由于溶腔的压力变化，溶腔围岩应力随着变化，且随着溶腔压

力的减小，溶腔围岩的稳定性变差。图 8-10 至图 8-15 为 1、2 号腔在三个不同时期的大小主

应力变化图。

表 8-3 溶腔腔顶应力变化表

1 号腔腔顶应力/MPa 2 号腔腔顶应力/MPa 3 号腔腔顶应力/MPa 4 号腔腔顶应力/MPa 不同时期 大主应力 小主应力 大主应力 小主应力 大主应力 小主应力 大主应力 小主应力 溶腔期间 -17.43 -42.73 -16.23 -38.28 -17.98 -39.00 -17.17 -36.94 10 年废弃期 -17.18 -41.15 -16.05 -37.06 -17.84 -38.37 -17.06 -36.24 10 年储气期 -11.91 -35.02 -10.95 -31.05 -12.65 -34.69 -12.31 -32.69 图 8-10 溶腔结束时 1、2 号溶腔的最大主应力图 图 8-11 溶腔废弃 10 年时 1、2 号溶腔的最大主应力图

图 8-12 溶腔储气 10 年时 1、2 号溶腔的最大主应力图

图 8-13 溶腔结束时 1、2 号溶腔的最小主应力图

图 8-15 溶腔储气 10 年时 1、2 号溶腔的最小主应力图

（3）溶腔蠕变损伤区演化特征。

数值模拟结果表明，在盐岩溶腔、溶腔报废及储气期内均发生了不同程度的蠕变现象。

蠕变主要产生在溶腔围岩，由于四个溶腔的间距较小且岩柱的偏应力水平相对其他位置都要

大一些，因此随着时间的变化，溶腔之间岩柱的蠕变特征尤其明显且有互相贯通的趋势。储

气 10 年后，溶腔之间的岩柱部位都有蠕变变形现象，蠕变最大值由溶腔围

岩转移到溶腔岩柱处，不同时期的蠕变区如图 8-16 至图 8-17 所示。由图可以看到在 1

号井和 4 号井之间的岩柱产生的蠕变很大，这将影响到储气库的长期稳定性。因此在储库压

力设计时应该充分考虑这些因素的影响。

图 8-16 10 年流变后溶腔围岩蠕变图

图 8-17 20 年流变后溶腔围岩蠕变图

随着离腔顶距离的增加，蠕变有逐渐减小的趋势，如图 8-18 所示。

0.00%

0.05%

0.10%

0.15%

0.20%

0

2

4

6

8

10

距腔顶距离(m）

粘

性

应

变

1号腔

2号腔

3号腔

4号腔

图 8-18 四个溶腔蠕变距腔顶距离变化图

另外，作为储气库封孔的重要技术参数的套管鞋的高度取决于腔顶蠕变损伤区的范围，

因为盐岩的蠕变损伤区一旦发展到腔顶上方套管鞋的位置，将影响套管的正常工作，这是工程

设计决不允许的。从本次计算的结果来看，四个溶腔中，3 号井的蠕变范围最大，接近 7.80m，

但蠕变应变不大；2 号井尽管影响范围不大，但蠕变应变最大，达到 0.18%。随着离腔顶距离

的增加，蠕变有逐渐减小的趋势。因此，在实际设计时，应充分考虑溶腔腔顶的蠕变范围和蠕

变应变值，避免将套管鞋布置在蠕变影响的范围之内。

（4）塑性区分布。

塑性区的出现是在溶腔完成之后，在溶腔腔底特别是溶腔与溶腔之间（岩柱）的区域易

出现塑性区。如图 8-19 所示，塑性区范围一般在 6～10m 之间。塑性区随盐岩的蠕变扩展范

围变化不大，而且溶腔之间塑性区也没有贯通。

图 8-19 20 年流变后溶腔塑性区分布

（5）结论。

1）由于盐岩的蠕变性，较低的溶腔压力将产生较大的围岩变形。同时由于溶腔储气压力

的变化，溶腔围岩应力随着变化，且随着溶腔压力的减小，溶腔围岩的稳定性变差。

2）对于废弃溶腔群盐穴地下储气库，盐岩最大蠕变应变发生在溶腔之间的岩柱区域，随

着时间的增长，溶腔有相互贯通的趋势。因此，在储库压力设计时，最好设计成同采同注的工

作方式，以确保储库的长期稳定性，以免带来不必要的后果。

3）盐岩溶腔和岩柱在偏应力作用下产生蠕变，储库设计中可以根据腔顶周围产生蠕变变

形的范围确定套管鞋的合理位置，由于盐岩蠕变损伤，应确保套管鞋埋设在蠕变损伤区之外。

8.4 深埋引水隧洞的稳定性分析

裂隙岩体渗流场与应力场的耦合是岩体水力学研究的热点问题。裂隙岩体的渗流场受应

力场影响很大，应力场的改变导致裂隙的张开或闭合，进而影响岩体的渗透性，从而引起岩体

渗流场的改变。而渗流场的改变将改变渗透体积力的分布，必将引起应力场的改变。

锦屏二级水电站位于四川省凉山州境内的雅砻江干流上，是雅砻江上水头最高、装机规

模最大的一座引水式地下电站。利用锦屏大河湾 150km 的天然落差，截弯取直开挖隧洞引水

发电。引水隧洞贯穿锦屏山，具有埋深大、洞线长、洞径大的特点，是锦屏二级水电站枢纽最

重要的组成部分。隧洞所在区域内地形地质条件复杂、地下水活跃，高地应力和高地下水是影

响引水隧洞围岩稳定性及衬砌结构安全性的最主要因素。任旭华、王美芹等采用考虑降雨入渗

渗流场分析的有限单元方法，对锦屏二级水电站深埋引水隧洞的外水压力进行了研究，提出富

水区深埋隧洞渗流控制“以堵为主，堵排结合”

；又对隧洞在不同渗控方案所形成的外部水环

境条件下围岩和衬砌的工作性态进行了系统的比较研究和评价，得到了一些高地下水位条件下

深埋引水隧洞的支护设计方面的结论；吴世勇等应用弹塑性有限元法分析了锦屏二级水电站引

水隧洞开挖和支护过程中围岩的应力和变形规律及特征。

目前，大多采用有限元法对引水隧洞围岩渗流场、应力场和支护衬砌的受力特征等进行

研究，且通常采用的流固耦合渗流计算将渗透系数视为常量，没有考虑渗透性因开挖扰动而引

起的动态变化特性，事实上，对于工程岩体而言，当有外加载荷作用时，不但会引起裂隙的张

开，还会引起裂隙的闭合，在裂隙张开和闭合的过程中还会使工程岩体产生弹塑性变形或蠕变

等宏观特征，目前国内外很多学者已作了这方面的研究工作：Gang Han、Maurice B. Dusseault

研究了石油井非固结砂土应力和孔隙度及渗透系数的关系；蒋中明、Dashnor HOXHA 提出了

有限元计算过程中通过调整孔隙度和渗透系数来实现热–水–力耦合分析非线性孔隙弹性的研

究方法；任长吉、黄涛研究裂隙岩体中应力场改变对地下水渗流场的影响作用机理，推导出渗

透系数 K 和给水度 μ 受应力影响的表达式；刘希亮、朱维申等对高应力作用下深部含水层 19

个试样的渗透系数作了均匀试验，研究表明：渗透系数与有效应力呈指数衰减的变化规律；杨

天鸿、唐春安等在 Biot 基本方程的基础上，增加一个反映渗透系数和孔隙变化率关系的耦合

项，提出了岩石损伤演化过程渗流－应力耦合方程。

本文从多孔介质的角度，以体积应变为切入点，建立孔隙度、渗透系数与体积应变之间

的动态演化关系，并定义了反映围岩破坏的损伤变量，建立岩体的应力渗流耦合弹塑性损伤本

构模型，以锦屏二级水电站引水隧洞长期稳定性为研究背景，分析引水隧洞围岩稳定性和衬砌

受力特征。

8.4.1 裂隙岩体应力渗流耦合本构模型

1．基于多孔介质的裂隙岩体有效应力原理

应用李培超等在论文中的观点，假定工程岩体为多孔介质，推导出基于多孔介质的有效

应力原理。该原理包含多孔介质的结构参数——孔隙度



，代替了其他有效应力公式中用的较

多的经验参数（比如常用的 biot 常数），其有效应力原理为：

(1

)

s ij ij

p

ij







 



 

（8-5）

式中，

s ij



为固体颗粒间的应力， p 为孔隙流体压力，

为孔隙度。令 (1



 

)

_ijs





_ij



为固

体骨架的有效应力（Terzaghi 有效应力），则 Terzaghi 有效应力为：

ij ij

p

ij



 





 

（8-6）

该原理较好地体现了多孔介质中的流固耦合效应。对于饱和岩体而言， p 就是孔隙水的压力。

2．裂隙岩体渗透系数动态演化模型

对于流固耦合工程岩体，孔隙度和渗透系数率 K 等参数将随岩体的应力状态不同而发生

动态变化，因此有必要建立流固耦合作用下的动态模型。

根据李培超等建立的饱和多孔介质流固耦合渗流的数学模型，可得到多孔介质孔隙度与

体积应变、温度、应力等有如下关系：

0

(1

)

1

(1

/

)

1

V s s

p K

T











 

 







（8-7）

式中，



₀

为初始孔隙度；



_V

为体积应变。其表达式为

b _{11 22 33} V b

V

V



















，

K 为多

_s

孔介质骨架固体颗粒的体积弹性压缩模量，



_s

为热膨胀系数。

若不考虑渗流工程中温度和骨架颗粒的体积变化，无扩容现象时，其孔隙度可由式（8-7）

得到

0 0

(1

)

1

1

1

V V V

















 







（8-8）

在扩容条件下，可以得到在压缩条件下孔隙度



的动态演化模型为：

0 0

(1

)

1

1

1

V V V

















 







（8-9）

根据文献[14]，可由 Kozeny-Carman 方程导出渗透系数与体积应变的关系式为：

3 0 0 0 0

(

/

)(1

)

1

1

1

V s s V

T

p K

K

K















  







_





_



_

_

（8-10）

同样，如果不考虑温度和材料骨架颗粒的体积变化，则可以得到等温渗流过程中渗透系

数的动态演化模型为：

3 0 0

1

1

1

V V

K

K













_



_



_

_

（8-11）

对于扩容后的情况，采用同样的分析方法，可以得到压缩条件下渗透系数的动态演化模

型为

3 0 0

1

1

1

V V

K

K













_



_



_

_

（8-12）

公式（8-9）～（8-12）反映岩体在应力作用下使裂隙闭合时将使孔隙度和渗透系数减小的特

点，以及当微裂隙扩展使得岩体产生扩容现象时，其孔隙度和渗透系数将相应地增大的特性。

3．裂隙岩体弹塑性应力渗流耦合损伤模型

裂隙岩体渗流对围岩的应力场有较大的影响，其渗透性与应力状态也密切相关，由于节

理裂隙的分布复杂且裂隙损伤演化依赖于加载历史，

很多学者应用自制方法或统计观点建立简

化模型研究，但他们的研究成果仅限于描述岩体弹性阶段的损伤，其损伤演化特性未能和岩体

的渗透性建立很好的关联关系。基于前述的裂隙岩体孔隙度动态演化的概念，建立岩体的损伤

演化的概念如下：

0 0 s











 



（8-13）

式中，



₀

为裂隙岩体初始孔隙度；



_s

为岩体达到强度极限时孔隙度。

当材料完全破坏时，也就是孔隙度







_s

时，损伤变量

  ；当材料没有损伤时，即

1







₀

时，损伤变量

 

0

。

将式（8-8）、式（8-9）分别代入式（8-13），可得到损伤演化方程为：

0 0 0 0 0 0 0 0

0

1

1

V V s s V V s









































_









 









_













_



弹性压密阶段

扩容前

扩容后

（8-14）

式（8-14）即为裂隙岩体损伤演化模型。由此可知，损伤变量  仅与材料的体积应变



_V

有

关，其变化规律不仅与弹性变形相关，还与塑性变形密切相关。从损伤变量的表达式可知，此

处定义的损伤变量实际上为一标量损伤。

该定义可以较好地描述岩石材料在不同应力状态下的

损伤演化规律。

根据前述有效应力原理，可得到饱和岩体的增量弹塑性损伤本构关系为：

0 0 0

d

(1

)

(d

d

)

d

d

(d

d

)

d

d

p ij ijkl kl kl w ij w ij p s ijkl kl kl w ij w ij s

E

p

p

E

p

p















 



















 



















（8-15）

式中，

0 ijkl

E

为无损伤岩体的刚度矩阵；



_kl

为总应变；



_klp

为塑性应变。

8.4.2 锦屏二级引水隧洞稳定性分析

以锦屏二级水电站引水隧洞为工程背景，采用大型有限元程序 ABAQUS 并对其进行二次

开发，将以上建立的孔隙度、渗透系数、损伤动态演化模型嵌入到程序中，通过计算，与孔隙

度、渗透系数取定值及不考虑岩体损伤的固－液耦合计算进行分析对比，研究引水隧洞在运营

期围岩的长期稳定性和衬砌结构受力特征。

1．工程概况

锦屏二级水电站引水隧洞贯穿锦屏山，长度约 16.67km，四条引水隧洞平行布置，引水隧

洞之间的中心轴线距设计为 60m，约为隧洞开挖洞径的 4.6 倍，隧洞之间的净岩体厚度为 47m，

约为隧洞开挖洞径的 3.6 倍，开挖洞径 13m，全线埋深较大，一般埋深为 1500～2000m，最大

埋深约 2525m，属于深埋长隧洞。引水隧洞所在区域地形地质条件复杂、地下水活跃，受大

气降水补给，地下水极为丰富。

引水隧洞沿线主要为三迭系中、上统的大理岩、灰岩、结晶灰岩及砂岩、板岩，其中大

部分为碳酸盐地层，其次为砂板岩，如图 8-20 所示。本文重点研究板岩地段隧洞围岩和衬砌

的水压力分布及受力特征。

引水隧洞断面型式为马蹄形，隧洞开挖和衬后断面均采用四心圆马蹄形断面，即隧洞底

部采用曲率半径较大的圆弧底，如图 8-21 所示。

K0 (m) 水平距离 1000 500 0 2000 1500 3000 2500 3500 4000 4500 （m）程 高 2569 K2 K1 1022 2132 K4 K3 3167 K6 K5 ° 4460～4469 T 2Z F25 T 1 5 F 6 F 雅 砻 江 3 3 T2Z T T 大 理 岩 6 绿泥石片岩 地层时代 T2b F x x x x z c z z z c c c 角砾状大理岩 图 例 细 砂 岩 中 砂 岩 粗 砂 岩 板 岩 角砾岩 推测地下水位线 条带状云母大理岩 断层及其编号 结晶灰岩 泥质灰岩 建筑物轮廓线 本次计算段 水　位　线 图 8-20 引水隧洞纵剖面图（部分）

砼 ° ° 总厚度60 图 8-21 引水隧洞标准断面图（长度单位：cm）

2．计算模型与参数

根据引水隧洞的布置特点，在隧洞埋深 1830m 处建立了平面应变有限元渗流计算模型。

计算范围为 640m×460m。模型的 Y 轴方向为重力方向，X 轴方向为平行于河谷方向，采用四

边形单元进行剖分，共划分 7010 个单元，7135 个节点，衬砌厚度 60cm，也采用平面四边形

单元。有限元网格模型为图 8-22，四个引水隧洞代号分别为 A、B、C、D。

图 8-22 模型网格图

边界条件为：模型左右两侧面施加 X 方向的位移约束，底面施加 Y 方向的位移约束；由

于所建模型上表面距地面 1600m，故上表面施加 44.8MPa 的压应力，上表面初始水头为 300m，

底面初始水头为 760m，左右两侧面施加沿重力方向梯度变化的水头压力。根据现场勘察和试

验成果，数值计算时围岩和结构力学参数取值如表 8-4 所示。岩体材料采用 Drucker-Prager 屈

服准则。

根据引水隧洞现场地应力的测试成果，引水隧洞初始应力场以自重应力场为主，中间主

应力平行于河谷方向（即近似垂直于隧洞轴线方向）

，其侧压力系数为 0.87。岩体充分扰动时

孔隙度取



_s

=0.01。

表 8-4 材料力学参数 材料 弹模 E/GPa 泊松比  密度 /×103kg/m3 内聚力 c/MPa 内摩擦角 /（°） 抗拉强度 t/MPa 孔隙度 /0 渗透系数 /m·s-1 岩体 20 0.25 2.8 1.5 42 2 0.003 4.63×10-9 衬砌 25 0.167 2.5 - - - 0.12 1.0×10-10

计算模型命令流：

**节点命令流 *Node,nset=nall 1, -84.0999985, 0., 0. 2, -83.5, 0., 0. 3, -85.4038086, 4.59619427, 0. 4, -83.5556107, 0.848420262, 0. 5, -83.7214813, 1.68232381, 0. 6, -83.9947815, 2.48744226, 0. 7, -84.3708344, 3.25, 0. 8, -84.8432007, 3.95694923, 0. 9, -85.8280716, 4.17192984, 0. 10, -84.1504745, 0.770104527, 0. 11, -84.3010406, 1.52703238, 0. 12, -84.5491104, 2.25783229, 0. 13, -84.8904495, 2.95000005, 0. 14, -85.3192139, 3.59169245, 0. 15, -94.1719284, 4.17192984, 0. 16, -94.5961914, 4.59619427, 0. 17, -90., 5.9000001, 0. 18, -90.7701035, 5.8495245, 0. 19, -91.5270309, 5.69896221, 0. 20, -92.2578354, 5.45088911, 0. ……… ……… 7135, 94.389473, 5.66742706, 0. **单元命令流 *Element, elset=eall,type=cpe4RP 1, 42, 53, 97, 96 2, 53, 52, 112, 97 3, 52, 51, 127, 112 4, 51, 50, 142, 127 5, 50, 49, 157, 142 6, 49, 48, 172, 157 7, 48, 47, 187, 172 8, 47, 46, 202, 187 9, 46, 45, 217, 202 10, 45, 44, 232, 217 11, 44, 43, 247, 232

12, 43, 41, 69, 247 13, 96, 97, 98, 95 14, 97, 112, 113, 98 15, 112, 127, 128, 113 16, 127, 142, 143, 128 17, 142, 157, 158, 143 18, 157, 172, 173, 158 19, 172, 187, 188, 173 20, 187, 202, 203, 188 ……… ……… 7010, 7135, 859, 548, 6883

材料力学参数：围岩体定义了力学参数随体积应变的变化关系，相关命令流如下：

*Material, name=eall **定义干密度 *Density 0.0028,0 **定义材料参数随体积应变变化 *Elastic,dep=1 20000 , 0.25,, 0.00E+00 19999.97944, 0.25,, 1.00E-07 19999.75332, 0.25,, 1.20E-06 19997.1221 , 0.25,, 1.40E-05 19958.89482, 0.25,, 2.00E-04 19876.73375, 0.25,, 6.00E-04 19671.61819, 0.25,, 1.60E-03 17964.68307, 0.25,, 1.00E-02 *Drucker Prager 54.8, 1., 0

*Drucker Prager Hardening 3.77,0. **定义渗透系数随体积应变的变化关系 *PERMEABILITY,SPECIFIC=0.01 4.63E-09 ,0.003,, 0.00E+00 4.63E-09 ,0.003,, 1.00E-07 4.63001E-09,0.003,, 1.20E-06 4.63013E-09,0.003,, 1.40E-05 4.63185E-09,0.003,, 2.00E-04 4.63556E-09,0.003,, 6.00E-04 4.64483E-09,0.003,, 1.60E-03 4.72306E-09,0.003,, 1.00E-02 ***定义场变量子程序

*USER DEFINED FIELD *depvar

4

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE USDFLD(FIELD,STATEV,PNEWDT,DIRECT,T,CELENT, 1 TIME,DTIME,CMNAME,ORNAME,NFIELD,NSTATV,NOEL,NPT,LAYER, 2 KSPT,KSTEP,KINC,NDI,NSHR,COORD,JMAC,JMATYP,MATLAYO, 3 LACCFLA) C INCLUDE 'ABA_PARAM.INC' C CHARACTER*80 CMNAME,ORNAME CHARACTER*3 FLGRAY(15) DIMENSION FIELD(NFIELD),STATEV(NSTATV),DIRECT(3,3), 1 T(3,3),TIME(2) DIMENSION ARRAY(15),JARRAY(15),JMAC(*),JMATYP(*),COORD(*) real K C if(kstep.gt.1)then CALL GETVRM('E',ARRAY,JARRAY,FLGRAY,JRCD,JMAC,JMATYP, 1 MATLAYO,LACCFLA) EPS1=ARRAY(1) EPS2=ARRAY(2) EPS3=ARRAY(3) EV=EPS1+EPS2+EPS3 if(ev.gt.0)then !D means damage D=(0.003-(0.003+EV)/(1+EV))/(0.003-0.1) !Qmeans ratio Q=(0.003+EV)/(1+EV) !K means permeability K=4.63*10**(-7.)*(1+EV/0.003)**3/(1+EV) ELSE D=(0.003-(0.003-EV)/(1-EV))/(0.003-0.1) Q=(0.003-EV)/(1-EV) K=4.63*10**(-7.)*(1+EV/0.003)**3/(1+EV) ENDIF field(1)=ev statev(1)=field(1) statev(2)=D statev(3)=Q statev(4)=K C if(kstep.eq.2.and.time(1).gt.3000.and.time(1).lt.3600)then C if(dtime.gt.30)dtime=20 C endif

if(kstep.eq.3.and.time(1).gt.3000.and.time(1).lt.3600)then if(dtime.gt.30)dtime=20 endif if(kstep.eq.4.and.time(1).gt.172000.and.time(1).lt.172800)then if(dtime.gt.60)dtime=60 endif if(kstep.eq.5.and.time(1).gt.3000.and.time(1).lt.3600)then if(dtime.gt.30)dtime=20 endif if(kstep.eq.6.and.time(1).gt.2590000.and.time(1).lt.2592000)then if(dtime.gt.60)dtime=60 endif if(kstep.eq.7.and.time(1).gt.3200.and.time(1).lt.3600)then if(dtime.gt.30)dtime=30 endif endif RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

初始条件定义如下：

饱和度：此处研究为饱和岩体。

*INITIAL CONDITIONS,TYPE=SATURATION nall,1

初始水压：定义模型上下边界的初始水压力。

***围岩体初始水压

*INITIAL CONDITIONS,TYPE=PORE PRESSURE nall, 3, 230, 7.6,-230

***衬砌初始水压

*INITIAL CONDITIONS,TYPE=PORE PRESSURE cqptn,0

初始有效应力：

*INITIAL CONDITIONS,TYPE=STRESS,geostatic eall, -41.8,230, -50.094,-230, 0.87,0.46

初始孔隙比：

***定义围岩体孔隙比 *INITIAL CONDITIONS,TYPE=RATIO nall,0.003 ***定义衬砌孔隙比 *INITIAL CONDITIONS,TYPE=RATIO cqptn,0.0008

3．计算工况与步骤

应用本文所建立的应力渗流耦合模型及损伤演化方程，研究引水隧洞在运营期围岩和衬

砌水压力分布规律及受力变形特征，并将本文的研究成果与传统的不考虑岩体渗透系数变化的

结果进行对比分析，数值模拟的工况如下：

工况一：采用本文建立的耦合模型，考虑渗透系数及损伤动态演化特征。

工况二：传统应力渗流耦合模型，考虑渗透系数为定值，不考虑损伤演化。

两种工况均按以下计算步进行：

Step 1

建立初始地应力平衡，给定水头边界。

Step 2

毛洞开挖后，洞壁水压力骤降为 0MPa。

Step 3

施作隧洞衬砌，并定义衬砌透水边界。

Step 4

引水隧洞运行，内水压力 300m 水头。

Step 5

维修期开闸放水，历时 4 小时。

分析步骤的相关命令流如下：

*********************step-1--- *step step1：GEOSTATIC *geostatic *DLOAD EALL,GRAV,10.,0.,-1. *model change,remove linern **在模型上表面定义压力，是总应力！ *dsload upSurf,p,44.8 *boundary up, 8,,3 bot,8,,7.6 *CONTROLS,ANALYSIS=DISCONTINUOUS *OUTPUT,FIELD,FREQ=1 *NODE OUTPUT POR,u *OUTPUT,FIELD *ELEMENT OUTPUT SAT,POR,s,FLVEL,peeq,sdv,e,fv *END STEP ***********************step-2--- *step,inc=3000,unsymm=yes,nlgeom=yes step2：kaiwa *SOILS,utol=10,consolidation 4,360000 *model change,remove exca,liner *boundary wyinter,8,,0 ***sflow **ss,qd,1. *end step ************************step-3--- *step,inc=3000,unsymm=yes,nlgeom=yes step3：50 liner pore pressure

*SOILS,utol=15,consolidation 4,10800 *model change,add linern *Boundary,op=new right, 1 left,1 bot,2 up, 8,,3 bot,8,,7.6 cqinter,8,,0 ***sflow **intersurf,qd,1.

*CONTROLS,PARAMETERS=FIELD,FIELD=PORE FLUID PRESSURE 0.2,1.

*end step

**************************step-4--- *step,inc=3000,unsymm=yes,nlgeom=yes

step4：3 days give water *SOILS,utol=15,consolidation 4,864000 *dsload intersurf,p,3 *Boundary,op=new right, 1 left,1 bot,2 up, 8,,3 bot,8,,7.6 cqinter,8,,3

*CONTROLS,PARAMETERS=FIELD,FIELD=PORE FLUID PRESSURE 0.2,1.

*end step

****************************step-5--- *step,inc=3000,unsymm=yes,nlgeom=yes,amplitude=ramp step5：2 hous take out water

*SOILS,utol=10,consolidation 4.,14400 *dsload intersurf,p,0 *Boundary cqinter,8,,0

*CONTROLS,PARAMETERS=FIELD,FIELD=PORE FLUID PRESSURE 0.2,1.

*end step

4．数值计算结果分析

（1）围岩应力分析。

响，以工况一为例，B 洞拱顶在施作衬砌后，最大主应力为 0.04MPa，最小主应力为-13MPa，

充水运行后，最大主应力为 0.11MPa，最小主应力为-11.0Mpa，放水后，最大主应力为 1.40MPa，

最小主应力为 1.40MPa，隧洞放水后围岩受拉，工况一受拉范围为 1.71m，工况二为 1.23m。

图 8-23 为工况一时单隧洞放水后围岩受拉应力分布图。

图 8-23 放水后围岩受拉应力分布图（工况一）

（2）衬砌受力分析。

衬砌在经历充水运行和放水的过程中，其应力变化较大。以工况一为例：衬砌在施作后，

由于围岩渗流作用，衬砌与围岩接触的外边缘受拉，最大拉应力达到 0.75MPa，位于隧洞拱脚

部位，而内边缘承受压力，最大压应力为 7.0MPa，位于隧洞拱底和拱脚部位；充水运行后，

由于内水压力的作用，在拱脚衬砌外缘局部部位产生了较大拉应力，最大达到 4.32MPa；放水

后，拱顶衬砌外边缘压应力较大，最大达到 32.9MPa。图 8-24 和图 8-25 分别为单隧洞充水运

行的应力云图。从整个衬砌受力过程来看，工况一的衬砌最大拉应力略大于工况二，图 8-26

为衬砌在工况一、二下最大拉应力变化图。

图 8-24 充水运行后衬砌最大主应力分布图（工况一）

图 8-25 充水运行后衬砌最小主应力分布图（工况一）

1

2

3

4

5

3

4

5

计算步

最

大

拉

应

力

/

M

P

a

工况一

工况二

图 8-26 两种工况下衬砌最大拉应力变化图

（3）隧洞围岩和衬砌水压力分布规律。

两种工况下的围岩水压力分布规律基本一致，但工况一的水压力变化范围较大：隧洞开

挖后，周边水压力迅速降低，尤以竖直方向降低最快，隧洞拱顶上方 14.7m 范围内有明显降

低，拱底约有 11.6m；充水运行后，由于内水压力作用，隧洞拱顶上方 18.6m 范围水压力有明

显变化，拱底约有 16.1m；放水后，拱顶上方约 10.2m，拱底下方 5.25m 范围水压力有较大变

化。表 8-5 列出了两种工况下 B 洞洞顶和洞底衬砌外水压力，由表知：隧洞在运营期最大外水

压力达到 3.62MPa，这与文献[46]中所述的通过现场监测量得最大涌水压力 3.5MPa 较为接近。

图 8-27 和图 8-28 分别为工况一隧洞充水运行和放水后围岩水压力等值线分布图。

表 8-5 B 洞洞顶和洞底衬砌外水压力

工况一 工况二

计算步

拱顶/MPa 拱底/MPa 拱顶/MPa 拱底/MPa

施作衬砌 0.82 0.81 0.77 0.74 隧洞运行 3.54 3.62 3.54 3.62 放水后 2.09 1.72 2.15 1.46 图 8-27 隧洞充水运行后围岩水压力等值线分布图（工况一）（单位：MPa） 图 8-28 隧洞放水后围岩的孔隙压力等值线分布图（工况一）（单位：MPa）

（4）位移分析。

引水隧洞内水压运行后，隧洞四周产生向岩体深处的位移，放水后围岩位移向洞内发展，

且位移变化较大。以工况一拱顶处为例，充水运行后位移为 4.66mm，但放水后位移达到

5.23mm。表 8-6 为两种工况下不同计算步时 B 洞拱顶的位移变化值，由表可知，在隧洞拱顶

处，工况一的位移均大于工况二。

表 8-6 两种工况下 B 洞拱顶的位移变化表 计算步 工况一位移/mm 工况二位移/mm 允水运行 4.66 4.37 放水后 5.23 5.01

（5）围岩塑性区。

隧洞开挖后即出现了较大的塑性破坏区，充水运行后，内水压力对围岩应力状态有一定

改善，但短时间内的放水造成了隧洞围岩更大程度的破坏，围岩产生了更大的塑性区。两种工

况下塑性区变化规律基本一致，但工况一的塑性应变较大：隧洞开挖后，塑性区范围 5.2m 左

右，工况一的最大塑性应变 8.0×10

-3

_{，工况二为 7.0×10}

-3

_{；隧洞放水后，塑性区范围 6.6m 左}

右，工况一最大塑性应变 1.1×10

-2

_{，工况二为 1.0×10}

-2

_{。图 8-29 和图 8-30 分别为工况一单隧}

洞开挖和放水后的塑性区分布图。

图 8-29 开挖后围岩塑性区分布图（工况一） 图 8-30 放水后围岩塑性区分布图（工况一）

（6）隧洞围岩损伤演化特征。

隧洞开挖后即出现损伤，隧洞充水运行后边墙损伤范围约 2.2m，竖向范围约 3.0m，最大损

伤值 0.024；放水后损伤进一步增大，边墙损伤范围为 2.5m，竖向范围约为 3.4m；由于隧洞之

间的净距较大（47m），隧洞之间损伤区并未贯通。图 8-31 为单隧洞在放水后的损伤范围云图。

图 8-31 放水后围岩损伤范围云图

（7）结论。

1）从多孔介质的角度，基于有效应力原理，以体积应变为切入点，建立了饱和岩体孔隙

度、渗透系数的动态演化固液耦合弹塑性损伤本构模型。

2）根据建立的弹塑性损伤本构模型，结合 ABAQUS 有限元程序对锦屏二级水电站引水

隧洞运营期围岩和衬砌的受力和变形特征进行计算分析，与孔隙度、渗透系数取定值及不考虑

岩体损伤的常规计算结果进行分析对比，得到以下结论：考虑渗透系数和孔隙度的动态演化损

伤计算模型，更能反映实际情况：围岩变形较小，发生了更大程度的塑性破坏，衬砌承受拉应

力略小，但所受压应力很大（尤在隧洞放水后）

；运营期衬砌的外水压力很大，围岩在边墙方

向损伤范围略小于岩体自重方向，由于隧洞之间间距较大，损伤区域并未贯通。

3）通过对实际工程的计算结果表明，本文所建的饱和岩体孔隙度、渗透系数动态演化弹

塑性损伤本构模型合理，具有较好的实用性。

8.5 本章小结

本章研究内容为 ABAQUS 在深部岩体工程中的应用。首先分别介绍了 ABAQUS 在损伤、

蠕变和渗流损伤等方面的模拟功能，然后以地下高应力储气库和深埋引水隧洞长期稳定性为工

程研究对象，详细阐述了 ABAQUS 在这些方面的分析方法和应用的全过程。对研究类似深部

岩体工程的应用者有极好的学习和借鉴作用。