ABAQUS在深部岩体工程中的应用

第8章ABAQUS在深部岩体工程中的应用知识要点：;深部岩体工程特点;ABAQUS数值模拟功能;地下储气库的长期稳定性分析;深部引水隧洞的稳定性分析本章导读：本章主要研究ABAQUS在深部岩体工程中的应用。首先介绍了ABAQUS在此方面的模拟功能，然后分别以地下高应力储气库和深埋引水隧洞为工程应用对象，详细阐述了ABAQUS在这些方面的应用全过程。对研究类似深部岩体工程的应用者提供较好的学习和借鉴作用。8.1、深部岩体工程简介能源和矿产资源制约着国民经济的发展。随着浅部资源的逐渐减少和枯竭，矿物资源地下开采的深度越来越大。人类的工程活动已经深入到地下4000m以下的深度，如逾千米乃至数千米的矿山、水电工程埋深逾千米的引水隧道、核废料的深层地质处理问题、深地下防护工程等。在我国西部水电能源开发、南水北调、西气东送等工程中普遍涉及到深部岩体地质问题等等。如我国西部已建或将建的拉西瓦、二滩、天生桥等水电工程都遇到过高地应力、高水压等问题，甚至在浅表的大坝建基面的开挖时，都会发生由于高地应力释放产生的表面岩体成层剥落现象。西部地区的特殊地质、地形条件使得西部水电资源开发中地下深埋工程量很大。深部岩体所处的高地应力、高地温和高渗透水压力的特殊环境，伴随着深部岩体的力学响应明显有别于浅部岩体的力学行为。由于深部岩体受到各种荷载作用、岩体介质作用本身的复杂性、认知的不确定性和表现出一系列新的岩体地质力学特征，致使深部岩体的地质结构特征、裂隙渗透特性、变形机理、破坏规模和强度特征以及高渗透水压下的流变机理难于用传统的理论加以合理的解释，引起了国内外学者的极大关注，成为近几年来该领域的研究热点。由于深部岩体所处的复杂的地质环境，导致地应力高、温度高、渗透压高，加之较强的时间效应，使深部岩体的组织结构，基本行为特征和工程响应均发生根本性变化。深部岩体表现出明显的非线性力学特性。进入深部的工程岩体所属的力学系统不再是浅部工程岩体所属的线形力学系统，而是非线性力学系统，传统理论、方法与技术己经部分或全部失效。另一方面，深部岩体处于多场、多相耦合作用，地下水、瓦斯、温度均会对岩体的基本性质和工程响应带来很大影响。因此，研究深部岩体非线性损伤变形特性、多场耦合作用机理等方面的研究，将为深部岩体工程灾害的预测和控制提供理论基础，并将从新的角度和新的高度认识深部岩体的特殊性质，以指导工程实践，创造经济和社会效益。8.2ABAQUS数值模拟功能8.2.1损伤ABAQUS可通过多种方式定义损伤，诸如混凝土损伤塑性模型、采用用户子程序（USDFLD,UFIELD，UMAT等）。混凝土损伤塑性模型：适用于ABAQUS/Standard和ABAQUS/Expllicit两个模块。其功能包括：（1）能够模拟各种结构类型（梁、桁架、壳和实体）中混凝土和其他准脆性材料；（2）采用各向同性弹性损伤结合各向同性拉伸和压缩塑性理论来表征混凝土的非弹性行为；（3）主要用于钢筋混凝土结构分析，也能够用于素混凝土结构分析；（4）可采用加强筋模拟混凝土中的钢筋；（5）可以模拟低围压时，混凝土受到单调、循环或动载作用下的力学行为；（6）结合非关联多重硬化塑性和各向同性弹性损伤理论来表征材料断裂过程中发生的不可逆损伤行为；（7）周期荷载反向作用时，可以控制材料的刚度恢复；（8）可以定义与应变率相关的性状；（9）在ABAQUS/Standard中可以结合使用粘塑性正规本构方程来提高软化区域的收敛性；（10）材料的弹性行为应为各向同性和线性的。USDFLD,UFIELD的功能有相似之处，均可以在子程序中定义场变量并输出，同时还可以定义与场变量相关的变量，诸如损伤变量等。以USDFLD为例，其子程序的简要形式如下：SUBROUTINEUSDFLD(FIELD,STATEV,PNEWDT,DIRECT,T,CELENT,1TIME,DTIME,CMNAME,ORNAME,NFIELD,NSTATV,NOEL,NPT,LAYER,2KSPT,KSTEP,KINC,NDI,NSHR,COORD,JMAC,JMATYP,MATLAYO,3LACCFLA)CINCLUDE'ABA_PARAM.INC'CCHARACTER*80CMNAME,ORNAMECHARACTER*3FLGRAY(15)DIMENSIONFIELD(NFIELD),STATEV(NSTATV),DIRECT(3,3),1T(3,3),TIME(2)DIMENSIONARRAY(15),JARRAY(15),JMAC(*),JMATYP(*),COORD(*)CCCCCCCCCCCCCCCCCCC定义场变量并赋予FIELD(NFIELD)定义与场变量相关的损伤变量，赋给状态变量STATEV(NSTATV)CCCCCCCCCCCCCCCCCCCRETURNEND与子程序相关的INP命令流格式为采用DEPENDENCIES语句，如：*Elastic,DEPENDENCIES=n**n为场变量个数25000,0.25,,024950.04997,0.25,,0.0000124504.96683,0.25,,0.000120468.26883,0.25,,0.0013383.382081,0.25,,0.011.134998244,0.25,,0.05这里定义了等效蠕变应变CEEQ作为场变量，场变量CEEQ与弹性模量E的关系为：00.05ceeqEEe−⋅=（8-1）上式中：0E为初始弹性模量。接下来采用如下语句定义状态变量的个数：*USERDEFINEDFIELD*DEPVARn**状态变量的个数UMAT是用户材料子程序（User-definedMaterialMechanicalBehavior，简称UMAT）是ABAQUS提供给用户定义自己的材料属性的Fortran程序接口，使用户能使用ABAQUS材料库中没有定义的材料模型。UMAT子程序具有强大的功能，不仅能定义材料的本构关系和屈服准则，还能实现材料本构关系用于力学分析的任何过程，当然，UMAT中必须提供材料本构的雅可比（Jacobian）矩阵而实现增量计算。有关UMAT的具体介绍详见第章。8.2.2非线性蠕变岩体本构关系ABAQUS有限元分析中可以定义粘弹、粘塑、粘弹塑性等各种与粘性蠕变相关的本构关系。其中定义粘塑、粘弹塑性本构关系可采用*CREEP或*DRUCKERPRAGERCREEP来定义。下面以*DRUCKERPRAGERCREEP实现粘弹塑性本构关系为例加以介绍。Drucker-Prager蠕变模型岩土材料在某些条件下将发生蠕变。尤其当蠕变时间尺度和加载率处于同一数量级时需考虑蠕变与塑性的耦合，此时ABAQUS提供了一个蠕变模型来定义材料的“经典”蠕变行为，其塑性变形服从扩展Drucker-Prager模型。材料的蠕变与塑性变形密切相关（此处可通过定义蠕变流动势和试验数据），因此在材料属性中必须包含塑性及塑性硬化的定义。1.蠕变应力的定义采用Drucker-Prager屈服准则所用的蠕变势函数为双曲线函数，如果在ABAQUS/Standard中定义了蠕变，线性Druker-Prager塑性模型也采用双曲线塑性流动势。假定存在应力点的蠕变等倾面，该等倾面具有相同的蠕变“强度”，并由等效蠕变应力确定。当材料发生塑性变形时，需要等效蠕变面与屈服面一致，因此等效蠕变面可由屈服面等比例缩小得到。在p-q平面上，蠕变面与屈服面平行。如图8-1。在定义材料蠕变特性时，等效蠕变应力crσ可表示为：(tan)1(1tan)3(tan)1(1tan)3(tan)crqpqpqpβββσββ−⎧⎪−⎪⎪⎪−=⎨⎪+⎪⎪⎪−⎩图8-1p-q面屈服面和双曲线流动势函数示意图2.蠕变流动蠕变应变率采用与塑性应变率相同的双曲线流动势函数。即：220(tan)tancrGqpσψψ=∈+−（8-2）式中：.00,0plplεεσσ===为初始屈服应力；∈为偏移率，用来定义双曲线趋向于渐近线的速率（当偏移率趋向于零时，蠕变势函数趋向于直线）；ψ为p-q平面上的剪涨角。蠕变势函数为光滑连续的曲线，能保证蠕变流动方向唯一。在高围压作用下，该函数趋近于线性Druker-Prager流动势，同时与静水压力轴正交。当材料为非相关流动时，刚度矩阵为非对称阵，此时需要采用非对称矩阵存储和非对称计算方法。如果摩擦角和剪涨角相差不大，且模型中非弹性变形区域受到限制时，材料刚度近似对称矩阵会给出可以接受的收敛速度，从而不需要非对称矩阵算法。3.蠕变定律在ABAQUS/Standard中，蠕变通过等效的“单轴蠕变”定义被称为“蠕变定律”。ABAQUS提供了三种可供用户输入相关参数来定义的蠕变定律：时间硬化蠕变定律、应变硬化蠕变定律和Singh-Mitchell蠕变定律。对于很多的实际工程问题，由于材料蠕变的函数形式往往较为复杂，因此需要采用用户编制子程序CREEP来实现，采用CREEP子程序自行编制非线性蠕变本构模型程序与ABAQUS结合计算从而解决工程问题。ABAQUS可在CREEP子程序中定义非线性蠕变本构模型，其参与整个计算的具体流程图如图8-2：当通过单轴压缩屈服应力定义cσ定义蠕变时当通过单轴拉伸屈服应力定义tσ定义蠕变时当通过粘聚力定义蠕变时ABAQUS应变更新()()()totaltotaltotalnnnttttεεε+Δ=+ΔABAQUS产生新的总应变增量1()totalntε+ΔABAQUS第n+1步的平衡迭代收敛不收敛N=N+1,返回第一步减小时间步长ABAQUS在平衡时间tn时：Δt,ε,Δε非线性蠕变本构模型是否进入D-P屈服YN()()()etotalcnnntttεεεΔ=Δ−Δ()()()()etotalcPnnnnttttεεεεΔ=Δ−Δ−Δ应力更新（塑形要迭代求解）()()()nnnttttσσσ+Δ=+Δ图8-2CREEP子程序计算流程框图用户子程序CREEP的模版如下：**********************************************************SUBROUTINECREEP(DECRA,DESWA,STATEV,SERD,EC,ESW,P,QTILD,1TEMP,DTEMP,PREDEF,DPRED,TIME,DTIME,CMNAME,LEXIMP,LEND,2COORDS,NSTATV,NOEL,NPT,LAYER,KSPT,KSTEP,KINC)CINCLUDE'ABA_PARAM.INC'CCHARACTER*80CMNAMECDIMENSIONDECRA(5),DESWA(5),STATEV(*),PREDEF(*),DPRED(*),1TIME(2),EC(2),ESW(2),COORDS(*)定义DECRA**注：昀常用的是定义DECRA(1)、DECRA(5)当考虑材料膨胀属性时，定义DESWARETURNEND**********************************************************与用户子程序CREEP相关的命令流为：*DRUCKERPRAGERCREEP,LAW=USER8.2.3考虑渗流的岩体本构关系ABAQUS能够求解多孔介质的饱和渗流，非饱和渗流及二者的混合问题（诸如坝体渗流自由面）的计算。计算过程中可以考虑渗流重力的作用，并能够求解流体的总孔隙压力或超孔隙压力，渗透定律可采用达西定律或更广泛的非线性定律。在流体重力荷载不可忽略或较明显瞬态毛细吸力较明显，即“湿化作用“不可忽略的问题中，需要求解总孔隙压力。需要指出的是：ABAQUS在孔隙压力解的形式上有两种，分别为总孔隙压力和超孔隙压力，两者仅仅是重力荷载上占主导地位的情况下才有明显区别。当采用GRAV分布荷载类型来定义模型的重力荷载时，ABAQUS程序提供总孔隙压力解；其他情况（如采用分布荷载类型BX、BY、BZ定义重力荷载时）提供超孔隙压力解。下面分别举例加以说明：（1）采用Grav时************