统一强度理论模型嵌入ABAQUS软件及在隧道工程中的应用

王俊奇，陆 峰

统一强度理论模型嵌入ABAQUS软件及在隧道工程中的应用

王俊奇1，陆 峰2

（1.华北电力大学可再生能源学院水利水电工程系，北京 102206；2.中国水利水电科学研究院 科研规划处，北京 100044）

利用非线性有限元软件ABAQUS提供的二次开发功能，实现了统一强度理论本构模型的嵌入，完成了有理论解的承受内压作用的厚壁圆筒三维算例数值测试，以及采用该模型进行隧道开挖三维数值分析。结果表明：在ABAQUS中增加统一强度理论本构模型，丰富了材料单元库，提高了计算精度和效率，而且，通过算例验证和隧道开挖模拟，说明在岩土工程中，考虑材料的中主应力效应，可以充分利用材料强度，指导工程实践，节省造价。

统一强度理论；中主应力；隧道；ABAQUS有限元

ABAQUS软件属国际上最先进的大型通用非线性有限元分析软件之一，在几何、材料和接触非线性问题方面的分析能力居世界领先水平，其中包括了岩土工程中常用的 Drucker-Prager，Mohr-Coulomb和剑桥等众多材料本构模型［1，2］。在岩土工程实践中广泛采用的Mohr-Coulomb理论认为，抗剪强度和中主应力无关。不考虑中主应力的Mohr-Coulomb准则对材料强度的反映与试验结果有一定出入［3］；考虑中主应力效应，可以提高岩石的强度20%～30%［4］。因此，发展了统一强度理论，通过取用不同参数，考虑中主应力的影响［4－6］。

虽然ABAQUS有许多本构模型可供选择，国内也开发了土工分析中常用的 Duncan-Chang模型［7，8］，但 Duncan-Chang模型实际属于非线性弹性模型，处理相对简单。为弥补ABAQUS中尚无可以考虑中主应力的统一强度理论本构模型这一不足，本文基于ABAQUS提供的二次开发用户子程序，完成了统一强度理论弹塑性本构模型的接口开发工作，与ABAQUS强大的解决大型复杂土工问题的非线性数值分析能力结合，探讨考虑中主应力对隧道开挖计算塑性区范围的影响。

1 理论基础

1.1 统一强度理论表达式

统一强度理论是西安交通大学俞茂宏教授在1961年提出的双剪屈服准则、1983年的双剪强度理_论和1987年的双剪角隅模型基础上，概括提出的工程强度理论体系，在特定条件下，可以退化为现今各种强度理论。考虑材料的拉压强度差异效应，同时考虑中主应力的影响［4－6］。

统一强度理论可以表示为主应力形式，也可表达为用应力不变量和应力角I1，J2，θ与岩土工程中常用的2个强度参数c和φ表示的形式［4－6］。

交接处的角度θb可由F＝F＇的条件求得

1.2 弹塑性增量理论的一般表达式

加载可以理解为由小屈服面F到大屈服面Φ。加载函数和加载面不仅与应力状态｛σ｝有关，而且还取决于塑性应变｛εp｝和反映加载历史的强化参数к。用Φ表示加载函数，则加载条件或加载面（后继屈服面）为 Φ（｛σ｝，｛εp｝，κ）＝0。

根据流动法则推导得出弹塑性增量理论本构关系的表达式［9－11］为

其中

2 ABAQUS用户子程序及其开发

2.1 用户子程序概况［1，2］

ABAQUS为了方便和鼓励用户开发自己研究或感兴趣的本构模型，采用FORTRAN语言接口方式，提供了若干用户子程序（User Subroutines）和在编程时可以调用的实用工具（Utility Routines）。与开发用户材料力学本构模型直接相关的子程序是UMAT，按照约定，开发者需利用UMAT子程序定义其单元材料积分点的Jacobian矩阵，即材料本构关系的刚度系数矩阵

2.2 接 口

用户材料子程序（user-defined material mechanical behavior，简称UMAT）通过与ABAQUS主求解程序的接口实现与ABAQUS的数据交流。在输入文件中，使用关键字“user material”来定义用户材料属性。

2.3 子程序结构与编程概要

根据ABAQUS的约定，用户子程序体结构应至少包括6部分，分别是：ABAQUS约定的子程序题名说明、ABAQUS定义的参数说明表、开发者定义的局部变量说明表、开发者编写的程序代码段和子程序返回与结束语句等，编程代码目的是得到式（6）的具体表达。

以下就程序设计中弹塑性增量关系的一些问题分4小节作出说明。

2.4 弹塑性状态的确定［12］

弹塑性状态确定的一般步骤如下（仅讨论各向同性硬化情况）：

（1）根据当前增量步或迭代步的位移增量计算应变增量

（2）按弹性应力应变关系计算应力增量的预测值及应力的预测值（试探性应力）：

其中tσ是上一增量步或迭代结束时的应力值。

（3）按单元内各个积分点计算本增量步或迭代结束时的t＋Δtσ，t＋Δtεp等状态量。

① 计算屈服函数值 F（t＋Δt珟σ，t珔εp），然后区分 3种情况，即

若 F（t＋Δt珟σ，t珔εp）≤0，则该积分点为弹性加载，或由塑性按弹性卸载，这时均有

若 F（t＋Δt珟σ，tεp）＞0，且 F（t珟σ，t珔εp）＜0，则该积分点为由弹性进入塑性的过渡情况，应由

来计算弹性因子m。计算m是为了确定应力到达屈服面的时刻。试探应力到达屈服面的比例因子m，可以采用线性插值和屈服函数Taylor展开修正的方法得到［9，10］，也可以采用二分的方法获得，本文采用前者。

若 F（t＋Δt珟σ，t珔εp）＞0，且 F（t珟σ，t珔εp）＝0，则该积分点为塑性继续加载，这时令m＝0。

②对于①中前两种情况，均有对应于弹塑性部分的应变增量 Δε＇，即

③ 计算弹塑性部分应力增量Δσ＇，即

一般情况下用数值积分方法进行此积分，在积分过程中可以同时得到Δ珔εp。

④ 计算本增量步或迭代结束时刻的t＋Δtσ，t＋Δtεp，即

2.5 本构关系积分（应力更新算法）

对率形式的本构方程进行积分的算法称为应力更新算法，在弹性预测阶段，塑性应变和内变量保持固定，而在塑性修正阶段，总体应变保持不变。

本文采用的是基于显式积分的切向预测径向返回的子增量法［9－12］。

2.5.1 切向预测径向返回

所谓切向预测就是将欧拉方法用于式（13），得到应力增量的预测值，即

进一步得到应力的预测值

同时还可以得到 Δt珔εp。

因为式（16）所表达的算法是显式的欧拉方法，其中的 Dep（t＋Δtσ，t珔εp）是起点切线刚度，所以 Δ珟σ是在加载曲面的切线方向。同时由于加载曲面是外凸的，因此t＋Δt珟σ总是在加载面之外。但是屈服准则要求应力t＋Δtσ只能在加载曲面之上或者之内，所以常再采用径向返回的方法以求得满足屈服条件的t＋Δtσ。具体做法是令

其中r是比例因子，它由下式

得到。

虽然经过校正以后得到的t＋Δtσ是位于屈服曲面上的，但因为假设应变增量Δε和等效塑性应变t＋Δt珔εp均保持不变，所以这样的弹塑性状态并不是完全一致的。这种不一致性随增量步长的增加而增加。为了减小由于这种不一致引起的误差，可将上述方法和子增量法相结合。

本文在径向返回过程中，做了适当调整，采用文献［9，10］的映射方法直接求得应力增量的修正值。

2.5.2 子增量法［12－14］

子增量法就是将总的弹塑性应变增量分成若干个子增量，对于每个子增量，利用上述切线预测径向返回的方法确定应力增量，将每一个子增量结束时的弹塑性状态作为下一个子增量的初始状态。子增量有助于提高应力增量的计算精度，从而加快迭代的收敛速度。

本文采用的子增量数N由等效应力确定，即

式中：σ0y为初始屈服时对应的等效应力；珚σy为按试应力得到的等效应力；σ＇y为t＋Δtσ＝tσ＋mΔ珟σ对应的等效应力。

2.6 切线刚度法［13］

隐式后向欧拉积分方法中多采用一致切线刚度矩阵。本文采用向前欧拉显式积分方法，由前面的公式（4）可见，这里形成的弹塑性刚度矩阵是标准的切线模量矩阵，即非一致切线模量矩阵。

2.7 屈服面上奇异点处理

统一强度理论在空间屈服面有交线，在π平面上并不光滑，存在奇异点，需特殊处理，一般有2种处理方式，一是Koiter提出的塑性应变增量求和法，二是 Owen等建议的回代求导法［4－6，14，15］。本文采用的方法分别如下：

（1）对于θ＝θb点产生的奇异性，用矢量和的办法处理。

（2）对于θ＝0°和θ＝60°产生的奇异，用回代求导法。

3 典型算例数值测试

3.1 模型条件

选取一个成熟的问题作为实例，进行对比和验证。长厚壁圆筒承受内压p＝160 MPa作用，内壁半径a＝0.1 m，外壁半径b＝0.2 m，如图1所示（取四分之一），材料弹性模量 E＝2.1×105MPa，泊松比ν＝0.3，单轴屈服应力 σs＝240 MPa。

图1 有限元计算网格Fig.1 Finite element calculation meshes

根据塑性力学，若采用Tresca屈服准则，弹性极限压力和塑性极限压力分别为［16］

相应于所给材料屈服条件，得pe＝90 MPa，ps＝166.35 MPa，所受内压介于弹性极限与塑性极限之间，按照理论解，圆筒处于弹塑性状态。根据与理论解的比较，验证数值试验结果的正确性。

3.2 ABAQUS程序二维和三维比较

首先比较商业软件ABAQUS就所给模型的平面和三维计算结果，采用ABAQUS中Mohr-Coulomb准则分别计算。长圆筒作为平面应变问题，划分8结点等参单元单元12个，结点51个。三维模型沿厚度方向划为四层单元，采用20结点等参块体单元12个，共有结点122个。网格轮廓见图1。

Mohr-Coulomb准则内摩擦角取 0°时退化为Tresca准则，但 ABAQUS软件中的 Mohr-Coulomb准则不允许内摩擦角取0°，计算时取一小值0.3°，剪胀角取0°，以模拟Tresca准则。计算得到二维和三维的结果完全一样，沿径向轴的位移和应力见表1。

3.3 统一强度理论UMAT程序结果比较

3.3.1 应力和位移

按照统一强度理论，当b＝0时，退化为Tresca准则。应用自编的UMAT程序，就三维模型计算得到沿径向轴的位移和应力见表1。比较表1中的数据，二者吻合很好，同时说明所编模型程序的正确。

由图可见，ABAQUS软件Mohr-Coulomb准则和b＝0时UMAT程序计算得到圆筒切向应力基本一致，径向位移在图示的结点上重合。

b增大时，圆筒周边切向应力最大值向筒内壁靠近，径向应力变化不大，位移则可以看出随b的增大而减小。

3.3.2 塑性区

图2 圆筒应力及位移（ABAQUS程序M-C准则和UMAT程序b＝0，b＝0.366，b＝1）Fig.2 Stresses and displacements of the cylinder obtained by the M-C criterion of ABAQUSprogram and b＝0，b＝0.366 and b＝1 in UMAT program respectively

表1 三维模型沿水平径向轴的位移和应力Table 1 Displacements and stresses of 3-D model along horizontal radial direction

图3 圆筒塑性区分布（等效塑性应变大于0.000 1）Fig.3 Distribution of the cylinder plastic zone with the equivalent plastic strain being more than 0.000 1

3.3.3 结果分析

将统一强度理论本构关系嵌入有限元程序，通过厚壁圆筒承受内压的算例与有关文献二维有限元分析结果对比，在ABAQUS中二次开发弹塑性程序获得成功。

4 工程应用

4.1 工程背景及材料参数

某工程隧洞开挖洞径D＝10.2 m，埋深620 m。拟以TBM法施工为主，钻爆法施工为辅。所选断面岩石为II级围岩，主要物理力学参数见表2。

表2 材料参数Table 2 Material parameters

4.2 有限元网格及边界条件

计算时依照对称原则，取结构的一半进行分析。有限元网格划分如图4。为了简化分析，自地面向下570 m的岩体部分其自重用施加等值面力来代替，约束条件，底部全约束，垂直边界法向约束。计算时不考虑分步开挖和支护，一次开挖成毛洞，考察位移和塑性区。

图4 有限元计算网格Fig.4 Finite element calculation meshes

4.3 计算结果分析

4.3.1 位 移

统一强度理论中参数b可以考虑中主应力效应，当b＝0时，退化为Mohr-Coulomb模型。本工程通过改变b值，考察中主应力效应对隧道工程的影响。将b＝0，b＝0.366和 b＝0.8时隧道开挖后的水平和竖直位移场作于图5和图6，图中可见，位移场形状一致，只是数值有别。将洞周侧壁水平方向及洞顶和洞底竖直方向的位移作于图7。可以看出，随b值的增大，3个方向的洞周位移逐渐减小。

图5 水平方向位移场Fig.5 Displacement fields in the horizontal direction

图6 垂直方向位移场Fig.6 Displacement fields in the vertical direction

图7 洞周位移比较Fig.7 Comparison for the peripheral displacements

4.3.2 塑性区

塑性区在一定程度上等同于洞室的破坏区，将b＝0，b＝0.366和 b＝0.8时隧道开挖后的塑性区作于图8，当b＝0时，即Mohr-Coulomb法则计算得到洞周塑性区出现范围为三层网格，b＝0.366时外层网格减小，b＝0.8时，两层网格，说明随b值的增大，塑性区范围在逐渐缩小。

图8 洞周塑性区范围比较（等效塑性应变大于0.0001）Fig.8 Comparison of plastic zones of the cylinder with equivalent plastic strain being more than 0.000 1

4.3.3 结果分析

在现有文献中也有采用统一强度理论研究洞室开挖的实例［4］，但仅限于二维模型，本工程采用三维分析，获得了与二维模型相似的结果，即考虑中主应力对洞室开挖破坏和位移是有不同的影响的。由工程算例可见，对应 b＝0，b＝0.366，和 b＝0.8三种条件，洞周位移随b增大而减小，塑性区范围也逐步减小。可见，随反映中间主切应力及相应面上正应力对材料破坏影响程度的系数b值的增大，可使材料强度在工程应用中得以充分发挥，对于指导工程支护措施和节约工程造价都有很强的实用上的意义。

5 结 语

研究开发工作表明：依据ABAQUS为用户提供的高起点的二次开发平台，在ABAQUS中嵌入统一强度理论本构模型是可行的。本文给出的具有理论解的典型算例结果显示了理想的计算精度。

同时，通过岩土工程中隧道开挖的工程实例计算表明，改变参数b实现中主应力效应，与不考虑中主应力的Mohr-Coulomb破坏准则比较，显示增大b值后，洞周位移及塑性区都有减小趋势，对于指导工程支护措施和节约工程造价有积极意义。

［1］ ABAQUS／Standard有限元软件入门指南［M］.庄 茁，译.北京：清华大学出版社，1998.（Instruction of ABAQUS／Standard［M］.ZHUANG Zhuo，translated.Beijing：Tsinghua University Press，1988.（in Chinese））

［2］ Hibbitte，Karlsson，SorensonINC（2002）.ABAQUS／StandardUser＇s Manual［S］.

［3］ 吴天行.土力学（第二版）［M］.冯国栋，译 .成都：成都科技大学出版社，1982.（WU Tian-xing.Soil Mechanics［M］.FENGGuo-dong，translated.Chengdu：Chengdu U-niversity of Science and Technology Press，1982.（in Chinese））

［4］ 俞茂宏 .双剪理论及其应用［M］.北京：科学出版社，1998.（YU Mao-hong.Twin Shear Theory and Application［M］.Beijing：Science Press，1998.（in Chinese））

［5］ 俞茂宏 .工程强度理论［M］.北京：高等教育出版社，1999.（YU Mao-hong.Engineering Strength Theory［M］.Beijing：Advanced Education Press，1999.（in Chinese））

［6］ 赵均海.强度理论及其工程应用［M］.北京：科学出版社，2003.（ZHAO Jun-hai.Strength Theory and Application［M］.Beijing：Advanced Education Press，2003.（in Chinese））

［7］ 徐远杰，王观琪.在ABAQUS中开发实现 Duncan-Chang本构模型［J］.岩土力学，2004，25（7）：1032－1036.（XU Yuan-jie，WANG Guan-qi.Development and implementation of Duncan－Chang constitutive model in ABAQUS［J］，Rock and Soil Mechanics，2004，25（7）：1032－1036.（in Chinese））

［8］ 张 欣，丁秀丽，李术才.ABAQUS有限元分析软件中Duncan-Chang模型的二次开发［J］.长江科学院院报，2005，22（4）：45－47.（ZHANG Xin，DING Xiu-li，LI Shu-cai.Secondary development of Duncan-Chang Model in ABAQUSsoftware［J］.Journal of Yangtze River Scientific Research Institute，2005，22（4）：45－47.（in Chinese））

［9］ 朱伯芳.有限单元法原理与应用（第二版）［M］.北京：中国水利水电出版社，1998.（Zhu Bo-fang.Finite Element Method Principle and Application［M］.Beijing：China Hydraulic and Hydroelectric Press，1998.（in Chinese））

［10］王 仁，熊祝华，黄文彬 .塑性力学基础［M］.北京：科学出版社，1982.（WANGRen，XIONGZhu-hua，HUANG Wen-bing.Elementary of Plastic Mechanics［M］.Beijing：Science Press，1982.（in Chinese））

［11］蒋友谅.非线性有限元法［M］.北京：北京工业学院出版社，1998.（JIANG You-liang.Nonlinear finite Element Method［M］.Beijing：Beijing Institute of Technology Press，1988.（in Chinese））

［12］王勖成.有限单元法［M］.北京：清华大学出版社，

2003.（WANG Xu-cheng.Finite Element Method［M］.Beijing：Tsinghua University Press，2003.（in Chinese））

［13］SMITH I M，GRIFFITHS D V.有限元方法编程，王 崧，周坚鑫，王 来，等译.北京：电子工业出版社，2003.（SMITH I M，GRIFFITHS D V.Programming the Finite Element Method［M］.Translated by WANG Song，ZHOU Jian-xin，WANG Lai，et al.Bejing：Publishing House of Electronic Industry，2003.（in Chinese））

［14］OWEN D R J，HINTON E.塑性力学有限元－理论与应用［M］.曾国平，刘 忠，徐家礼，等译.北京：兵器工业出版社，1989.（OWEN D R J，HINTON E.Finite Elements in Plasticity［M］.ZENGGuo-ping，LIU Zhong，XU Jia-li.trans.Beijing：Publishing Company of Weapon Industry，1989.（in Chinese））

［15］胡其志，周辉，杨雪强，屈服面角点奇异性的非平滑处理对塑性应变的影响［J］，岩土工程学报，2009，31（1）：66－71.（HU Qi-zhi，ZHOU Hui，YANG Xue-qiang.Effect on plastic strain through the non-smoothness management of corner singularity，Chinese Geotechnical Journal of Geotechnical Engineering，2009，31（1）：66－71.（in Chinese））

［16］王 仁，黄文彬，黄筑平 .塑性力学引论［M］.北京：北京大学出版社，1992.（WANG Ren，HUANG Wen-bin，HUANG Zhu-ping. Introduction of Plastic Mechanics［M］.Beijing：Beijing University Press，1992.（in Chinese） ）

Unified Strength Theory Constitutive Model Embedded Software ABAQUS and Its Application in Tunnel Engineering

WANG Jun-qi1，LU Feng2
（1.Department of Hydraulic and Hydropower Engineering，North China Electric Power University，Beijing 102206，China；2.Institute of Water Resources and Hydropower Research，Beijing 100044，China）

The unified strength theory constitutive model can be embedded nonlinear finite element software ABAQUSfor providing the user subroutine second development function.A numerical test using thick wall cylinder as an illustration is completed，which gives the theoretical resolution and three dimension analysis of tunnel excavation.All numerical results indicate that the approach of adding the model to the ABAQUSenriches the material element library and enhances the computation efficiency and precision，meanwhile，by the example verification，the tunnel excavation simulation demonstrates that considering influences of middle main stress in geotechnical engineering，the method can fully utilize the material strength enough，guide the engineering practice and reduce investment outlay.

unified strength theory constitutive model；middle main stress；tunnel；ABAQUSfinite element

TU311.41

A

1001－5485（2010）02－0068－07

2009-04-24

中国水利水电科学研究院开放研究基金资助项目（IWHRO2009005）

王俊奇（1971-），男，山西朔州人，博士，副教授，从事水工结构、岩土工程研究，（电话）15801452959（电子信箱）jqwangmail＠163.com。

（编辑：陈绍选）