第26卷 第5期
岩石力学与工程学报 Vol.26 No.5
2007年5月 Chinese Journal of Rock Mechanics and Engineering May，2007
深部软岩工程大变形力学分析设计系统
何满潮1，陈 新1，梁国平2，钱华山2，周永发2，庄小燕2
(1. 中国矿业大学 力学与建筑工程学院，北京 100083；
2. 中国科学院 数学与系统科学研究院飞箭软件有限公司，北京 100098)
摘要：采用飞箭软件公司的有限元程序自动生成系统作为软件开发平台，进行“深部软岩工程大变形力学分析设计系统”的合作开发工作。该软件系统具有如下特色：(1) 以和分解有限变形力学分析模型为主体，包括极分解力学分析模型供用户进行对比分析；(2) 强调非线性大变形力学设计特色：力学对策设计、过程设计和参数设计；(3) 包括多种岩土材料本构模型和工程材料单元，可以进行边坡工程、地基工程、地下洞室工程等大变形力学分析；(4) 具有我国自主知识产权。目前该软件的弹性、弹塑性大变形二维计算程序和界面已经开发完成。通过对典型算例的弹性数值分析，验证大变形程序的正确性。该软件系统的问世，可为深部软岩巷道等岩体大变形力学问题提供具有分析和设计功能的数值分析软件系统。
关键词：岩石力学；和分解；大变形；有限元软件系统；深部软岩工程
中图分类号：TU 452 文献标识码：A 文章编号：1000–6915(2007)05–0934–10
SOFTWARE SYSTEM FOR LARGE DEFORMATION MECHANICAL ANALYSIS OF SOFT ROCK ENGINEERING AT GREAT DEPTH
HE Manchao1，CHEN Xin1，LIANG Guoping2，QIAN Huashan2，ZHOU Yongfa2，ZHUANG Xiaoyan2
(1. School of Mechanics and Civil Engineering，China University of Mining and Technology，Beijing 100083，China；2. Beijing FEGEN Software Co.，Ltd.，Academy of Mathematics and Systems Science，Chinese Academy of Sciences，Beijing 100098，China)
Abstract：Up to now，there is no finite element software including nonlinear large deformation theory based on S-R decomposition theorem，in which the deformation gradient F is decomposed into the addition of the deformation and the rotation by a co-moving coordinate system method. In this work，such a software system waw developed，which is named as A Software on Large Deformation Analysis of Soft Rock Engineering at Great Depth. Two work groups，worked together on the software by using finite element program generator(FEPG) developed by FEGEN Software Company. The software has the following features：(1) it includes large deformation theory based on S-R decomposition theorem first time as well as one based on polar decomposition theorem for comparison，the latter is included by most commercial finite element software；(2) three design methods for nonlinear mechanical problem，i.e. countermeasure design，procedure design and parameter design，are fulfilled here for large deformation analyses；(3) it includes generally used elastoplastic constitutive models for geomaterials and element types for supports，and therefore is fit for large deformation mechanical analyses for slope，foundation and underground engineering；and (4) it has a copyright of our own country. At present，the
收稿日期：2006–09–26；修回日期：2006–12–18
基金项目：国家重点基础研究发展计划(973)项目(2006CB202200)；国家自然科学基金重大项目(50490273)；国家基金委创新群体基金资助项目(50221402)；教育部科学技术研究重大项目(10405)
作者简介：何满潮(1956–)，男，博士，1981年毕业于长春地质学院工程地质专业，现任教授、博士生导师，主要从事软岩、边坡及工程地质、地热
第26卷 第5期 何满潮，等. 深部软岩工程大变形力学分析设计系统 • 935 •
interface of the software has accomplished and programs of elastic and elastoplastic constitutive models based on two large deformation theories in two-dimension have been fulfilled，tested and implemented into the interface. Representative numerical examples confirm the rationality of the programs. The software will be a useful numerical tool for large deformation analyses and nonlinear mechanical design for soft rock engineering at great depth.
Key words：rock mechanics；S-R decomposition；large deformation；finite element software system；soft rock engineering at great depth
授经过多年研究，在拖带系描述法、张量分析和微分变换群等数学理论的基础上，于1979年提出了应变和转动的和分解定理，完善了Stokes在1945年提出的固体位移场分解理论，克服了以往大变形理论的缺陷，构成了非线性连续体力学的新体系。十几年来，以拖带系的和分解定理为基础，许多学者[10
～14]
1 引 言
随着开采深度的增加，岩爆、瓦斯突出、流变、底板突水等非线性动力学灾害现象日趋增多，严重影响了深部资源的安全高效开采。深部开采工程中产生的岩石力学问题已经成为目前国内外采矿及岩石力学界研究的焦点。
深部软岩巷道工程在巷道开挖后，巷道周围岩体的变形一般都较大，出现巷道冒顶、底膨和侧胀等现象都是围岩大变形的结果[1
～4]
研究了拖带系中的应变、转动、应力、应力
客观率的度量和变分原理，形成了基于和分解的有限变形理论。该理论的应用涉及岩石、混凝土、金属等刚体的大变形、非线性效应、断裂力学、细观力学、生物力学、分形几何、边坡稳定、软岩巷道支护、地质构造大变形、实验力学的激光频谱分析法以及物性方程等诸多领域，取得了令人瞩目的成果。
在深部软岩工程中，由于涉及到物理非线性、
～9]
。为了合理地进
行岩体变形分析，除了要很好地建立符合岩体变形的物理、力学特性本构模型外，还必须采用能正确描述大变形的非线性几何理论，才能得到较为合理的结果。
随着非线性几何场论的研究[5
不断取得进
展，目前形成了两大有限变形理论：(1) 采用固定坐标系描述方法，以Green非线性应变作为应变度量和Finger极分解定理得到的转动张量为转动度量的经典有限变形理论；(2) 采用拖带坐标系描述方法，基于应变和转动的和分解定理(additive decomposition theorem，简称S-R分解定理)的有限变形理论。
Green非线性应变张量克服了大转动时Cauchy应变张量在物体发生有限转动时对变形度量的失真，给出了大变形的确定性度量，因而被许多非线性力学著作采用。而Green应变张量没有相应的转动张量与之匹配，对全面研究大变形、大转动问题是一个重大缺陷。极分解定理定义的转动张量虽然在一定程度上弥补了Green应变张量没有相匹配的转动张量的缺点，然而它包含左、右极分解2个分解式，变形和转动的分解有先后之分，分解式不唯一，分解得到的伸长张量不适合作为变形的度量。
针对经典有限变形理论的这些问题，陈至达教
几何非线性和接触边界非线性等力学问题，因此其理论解的求解在数学上遇到非常大的困难，需要借助于有限元、有限差分和离散元等数值方法和软件。
目前，商业有限元软件的有限变形力学模型中，无一例外的都是采用Green非线性应变和Finger极分解定理来度量应变和转动，没有提供基于和分解的有限变形力学模型，这是很不全面的。为此，重大项目首席科学家何满潮教授组织了中国矿业大学(北京)力学与建筑工程学院和中国科学院数学与系统科学研究院飞箭软件公司2个单位的研究力量，策划和实施了“深部软岩工程大变形力学分析设计系统”的合作开发工作。该软件系统以软岩大变形理论和非线性力学设计方法[14]为理论基础，通过有限元程序自动生成平台FEPG5.2，建立深部三相(气、液、固)和三场(应力场、渗流场、温度场)耦合作用的数值模拟系统，应用于理论科研与工程实践，这对研究深部资源开采面临的岩石力学问题具有重大现实意义。
• 936 • 岩石力学与工程学报 2007年
2 和分解有限变形力学基本方程
2.1 位形变换的和分解定理
和分解有限变形力学采用嵌含在可变形体中的拖带系描述方法，在变形体的运动过程中，变形体中各质点的拖带坐标不变，而度量拖带坐标的基矢量则不断发生变化。
考察变形体的运动，设在初始平衡状态t0时刻，变形体中某质点在固定参考系中的空间位置为P，矢径为r，拖带坐标为xi，基矢量为gi。在变形后的某个平衡状态t时刻，该点运动到P′的空间位置，矢径变为R(R为正交变换，表现为点集的转动)，基矢量变为gi。变形前、后质点的协变基矢量分别定义为
∂r∆r⎫
lim=
⎪∂xi∆xi→0∆xi⎪
⎬ (1)
∂R∆R⎪gi=i=lim
∆xi→0∆xi⎪∂x⎭gi=
式中：S为对称变换，表现为点集的形变。
变换系数F的S-R分解分量形式为
⎫
⎪⎪1
Sij=(ui|j+ui|Tj)−(1−cosΘ)L⋅ki⋅L⋅kj⋅⎬ (7)
2⎪
⎪Rij=δij+L⋅ij⋅sinΘ+(1−cosΘ)L⋅ki⋅L⋅kj⋅⎭式中：Θ，L⋅ij⋅分别为平均整旋角和转轴方位张量分量，可分别定义如下：
Fji=Sij+Rij
Θ=±arcsin(−ωijωij/2)2⎫⎪
L=ωi/(sinΘ)
⋅i
j⋅
j
⎬ (8) ⎪⎭
式中：ωij为位移梯度的反对称部分，可表示为
ωij=(ui|j−ui|Tj)/2 (9)
位移梯度的物理分量为
ui|j=ui
j
giijj (10)
变形前、后质点的矢径与拖带坐标分别存在如下微分关系：
∂ri0⎫dr=idx=gidxi⎪
∂x⎪
⎬ (2)
∂Ri
dR=idx=gidxi⎪
⎪∂x⎭
应变张量、转角与转角方位的物理分量为 ⎫1
Sij=(ui|j+ui|Tj)−(1−cosΘ)L⋅ki⋅L⋅kj⋅⎪
2
⎪
⎪⎛1 i j⎞
⎪Θ=±arcsin⎜−ωjωi⎟
2⎪⎝⎠ ⎬ (11)
⎪1 i
L⋅ij⋅=ωj⎪
sinΘ⎪
⎪ i1 i iT
ωj=(u|j−u|j)⎪
2⎭2.2 广义Euler运动学公式
质点的位移矢量为
u=R−r (3)
变形前、后质点的协变基矢量之间存在如下变换关系：
gi=Fijgj，Fij=δij+uj|i (4)
在节2.1考察变形体的位形变换时，采用了初始拖带系基矢量gi为基准；对于非线性大变形运动，采用瞬时拖带系为参考系可更为准确地进行描述。
考察变形体从t时刻到t+∆t时刻的位形变换，拖带系的基矢量变换为：gi=gi+∆gi。
基矢量的变化率为
∂gi∆g
=limi=Vj||igj (12) ∂t∆t→0∆t
式中：δij为Kronecker符号；uj|i为质点位移在初始拖带系中的协变导数，可表示为
∂uj0jl
u|i=i+Γilu (5)
∂x
j
式中：0Γilj为初始拖带系的Christoffel符号。
因此，变形体的运动变换可以用质点局部嵌含坐标系的基矢量变换(式(4))来表示。
和分解定理指出，若变换系数F在形变体点集域内是单值连续的，处处具有一阶导数，则此运动变换可以唯一分解为正交与对称2个子变换的直和，即
式中：Vj||i为变形体中一点的速度在瞬时拖带系中的协变导数，可表示为
∂Vj
V||i=+VkΓkij (13) i
j
第26卷 第5期 何满潮，等. 深部软岩工程大变形力学分析设计系统 • 937 •
式中：Γkij为t时刻瞬时拖带系的Christoffel符号。
陈至达将Euler运动学公式推广到了变形体的运动，得到了变形体运动的速度梯度和分解式：
为σij。对几何非线性程度不高的弹性变形体，可采用如下全量弹性应力–应变关系：
..lk
σij=E⋅ijk⋅Sl (23)
j+R j (14) Vj||i=Sii
..l
式中：E⋅ijk⋅为弹性刚度矩阵。
..l
对各向同性材料，E⋅ijk⋅可表示为
..lililil
E⋅ijk⋅=λδjδk+µ(ggjk+δkδj) (24)
j，R j分别为该点的应变速率和局部转动式中：Sii
速度，可分别表示为
⎫ i=1(Vi||+Vi||T)Sjjj⎪⎪2
⎬ (15)
1 i=LiΘ =(Vi||−Vi||T)⎪Rjjjj⎪2⎭
2.3 应力、应力速率、平衡方程
在瞬时拖带系中，定义g(j)为单位化的逆变基矢量，g(i)为与g(j)满足对偶关系g(i)⋅g(j)=δij的协变基矢量：
gj=gjjg(j)，gi=g(i)
gii (16)
式中：λ，µ均为材料的Lame常数，可由实验测定。此时的物性关系为
σij=λSkkδ.ij+2µSij (25)
对几何非线性程度较高的弹性体或弹塑性体，变形较大，微元体将发生畸变。全量形式的本构定律没有考虑位形改变的影响，严格来说是不精确的，必须采用如下应力客观率形式本构关系：
∆i
..lk
σj=⋅ijk⋅Sl (26)
Euler应力张量可以在瞬时拖带系中分解为
..l
式中：⋅ijk⋅为材料的切线刚度矩阵。
σ=σ.ijgigj=σ.jg(i)g(j) (17)
i
利用式(21)可得到应力率形式的本构关系：
..l ki ki k
ij=⋅ijkσ⋅Sl+σkSj−σjSk (27)
式中：σ.ij为Euler应力张量的物理分量，可表示为
在实际计算中，应力、应变要采用物理分量。例如，对各向同性弹性变形体的平面应变问题，增量应力–应变关系为
σ.ij=σ.ijgjjgii (18)
应力的物质导数为
dσ.ijdt
ij−σ.ikνk||j−σ.ijνk||k (19) =σ
应力的物质导数不是客观的时间导数，与刚体的转动有关，去掉转动项得到客观应力速率σ⋅j为
∆i
1 1 2 1 2⎫
)∆S2+σ2∆S1∆σ1=(λ+2µ)∆S11+(λ−σ1
2 2 1⎪ 2 2 1⎪
∆σ2=(λ−σ2)∆S1+(λ+2µ)∆S2+σ1∆S2⎬ (28)
⎪ 1 1 1 1 1
∆σ2=−σ2∆S1+(2µ+σ1)∆S2⎪⎭2.5 变分原理
在变形体发生大变形时，外力与位移呈非线性关系，且边界随时间变动，合理的能量原理应以瞬时位形的功率形式表达，下面应用势能率原理来建立增量虚位移的变分方程。
在瞬时拖带系中，势能率原理为：对一切满足速度、角速度与应变速率谐调条件和速度边界条件的所有可能形变状态，满足平衡方程以及力边界条件的真实应力状态，将使其能量泛函J的变化率取驻值：
ii jdΩ− =σiδSδJfδVdΩ−TδVda=0 jiii∫∫∫
Ω
Ω
Sσ
即
σ⋅j=
∆i
dσ
i
.j
dt
k (20) +σ.ikRj
k+σiS k ij−σ.kjSσ.j=σj.jk (21)
∆i
平衡方程是建立在变形后的位形上的，可表示为
σ.ij||i+fj=0 (22)
式中：σ.ij||i为应力在瞬时拖带系中的协变导数，且
i
σ.ij||i=σ.ij，i+Γikσ.kj−Γijkσ.kj；fj为体积力。
2.4 物性关系
(29)
式中：为表面力。
• 938 • 岩石力学与工程学报 2007年
如果表面力、体积力不随时间变化，应力与荷载相关但不为时间函数，势能驻值可以表示为如下的全量形式的变分原理：
δJ=∫σijδSijdΩ−∫fiδuidΩ−∫Tiδuida=0 (30)
Ω
Ω
Sσ
和后处理子系统。通过界面将这些子系统组合起来。
系统的主流程和主界面分别如图1，2所示。
严格地说，只有率形式的变分原理是普遍适用的。但在小位移条件下，全量变分原理仍有相当的准确性。
3 基于FEPG有限变形力学程序开发
在深部软岩工程大变形力学分析设计系统中，所有的有限元程序都是利用其自动生成系统
(FEPG)开发的。
有限元程序自动生成系统(FEPG)由北京飞箭软件有限公司的创始人之一、中国科学院数学与系统研究院博士生导师梁国平研究员历经十余年的潜心研究而独创，该系统的主要设计思想是采用元件化的程序设计方法和有限元语言，根据有限元方法统一的数学原理及其内在规律，以类似于数学公式推理的方式，由有限元问题的微分方程表达式及其求解算法自动产生有限元程序。该系统适用于各种领域的各种工程与科学计算问题，突破了目前通用有限元程序只适用于特定领域和特定问题的限制。用户只需输入有限元方法所需的各种表达式和公式，即可自动产生所需的全部有限元程序，免去了大量繁琐的有限元编程劳动，并保证了程序的正确性和统一性。
在系统的有限元程序中，极分解增量分析模块采用了更新拉格朗日法[15]，和分解增量分析模块采用了更新拖带坐标法[13]，位移或位移增量的求解采用牛顿迭代法。
图1 系统的主流程图
Fig.1 Main flowchart of software system
图2 系统的主界面 Fig.2 Main interface of software
4 系统的功能
深部软岩工程大变形力学分析设计系统的最终目标是建立考虑深部三相(气、液、固)和三场(应力场、渗流场、温度场)耦合作用的非线性力学数值模拟系统，来研究深部资源开采面临的非线性岩石力学问题。
该系统将分阶段、分层次开发完成。第1阶段，该系统称为“深部工程软岩大变形力学分析设计系统LDEAS1.0版”。LDEAS1.0分4个功能独立的子系该版本具有二维固体应力场的非线性力学数值模拟分析功能，其有限变形力学分析模型包括极分解、和分解2个模块，计算模型包括弹性全量、弹性增量、弹塑性增量的平面应力以及平面应变问题共计16个分析计算程序。单元模型有实体单元、Goodman节理单元、殷有泉弹塑性节理单元和刚架、锚杆、桁架等结构单元，可以考虑桁架、锚杆、锚索支护和开挖以及建造施工过程。该系统采用力学对策设计、过程设计和参数设计的非线性大变形力学设计方法。
5 典型算例
本节采用板的有限转动与伸长、板的简单剪切、采场顶板垮落大变形和梁的弯曲大变形4个典型算例说明系统中各程序的正确性。
在增量分析中，将固体问题的加载途径划分成(0)(1)(Ν−1)(Ν)
第26卷 第5期 何满潮，等. 深部软岩工程大变形力学分析设计系统 • 939 •
其中，Ω(0)为初始状态，Ω(n)为其中任意一个中间状态，Ω(Ν)为最后的平衡状态。对第n(n=1， 2， "，N)个增量步，平衡状态由Ω(n−1)变换到Ω(n)。 5.1 算例1：板的有限转动与伸长
图3给出了板的有限转动与伸长变换。2 cm×2 cm的正方形板由0123→01′2′3′，整体转动45°，同时均等扩张为32cm×32cm的正方形。材料参数：E = 1.0×1010 Pa，ν = 0.3。
x22′
x2
1′
g2
g2
.ij(n)=.ij(n−1)+∆.ij(n)⎫⎪
⎬ (35) i(n)i(n−1)i(n)
.j=.j+∆.j⎪⎭
取N=10个增量步计算，网格为20×20，各增量步板的应变和应力解析解分别如表1，2所示。
表1 算例1中各增量步板应变解析解
Table 1 Analytical solution to strains of plate at some
steps in example 1
n 2 4
(n). 1
2(n)
.2
(n)
.1 2
0.116 667 0.275 557 0.441 243 0.598 166 0.741 759
0.116 667 0.275 557 0.441 243 0.598 166 0.741 759
0 0 0 0 0
x1
3′
g1g1
6 8 10
x1
表2 算例1中各增量步板应力解析解
Table 2 Analytical solution to stresses of plate at some steps
in example 1 Pa
n 2 4 6 810
(n)
. 1
(n)
.2 2
(n)
. 2
图3 算例1中板的有限转动与伸长变换 Fig.3 Rotation and expansion of a plate in example 1
1.619 048×109 1.619 048×109 3.551 792×109 3.551 792×109 5.256 592×109 5.256 592×109 6.617 906×109 6.617 906×109 448×109 7.679 448×109 7.679
0 0 0 0 0
取初始拖带系与固定直线直角坐标系一致，变
形体在平衡状态Ω(n)的固定坐标为
(0)
⎫x1(n)=[1+0.5K(n)]x1(0)−1.5K(n)x2
(31) (n)(0)(0)⎬x2=1.5K(n)x1+[1+0.5K(n)]x2⎭
本算例中板的有限转动与伸长变换的最终应变
其中，
和应力数值计算结果分别如图4，5所示。由此可知，
K(n)=n/N
采用和分解弹性平面应力增量程序进行计算，通过理论分析，可以得到第n个增量步的板应变增量计算公式为
(n)2(n)
∆.=∆.2=α(K(n−1))∆K⎫1⎪
(32) ⎬ (n)2(n)
⎪∆.2=∆.1=0⎭
(10)2(10)(10)
(a) .和.2 (b) .1 12
图4 算例1中板的最终应变数值计算结果
Fig.4 Numerical results of final strains of plate in example 1
其中，
α(K)=
1+5K
(33) 2
2(1+K+2.5K)
板的应力增量计算公式为
⎛E(n)2(n)(n−1)⎞1(n)⎫∆S∆.1=∆=−⎜⎟1.2.1.1⎪⎝1−ν⎠⎬ (34)
⎪(n)2(n)(n−1)
∆.1=∆.1=.∆K22⎭
(a)
(10)
.1
(10)
和.22
(10)
.2
(b)
图5 算例1中板的最终应力数值计算结果(单位：Pa) Fig.5 Numerical results of final stresses of plate in example
1(unit：Pa)
平衡状态(n)
• 940 • 岩石力学与工程学报 2007年
计算结果与理论分析的解析解完全吻合。 5.2 算例2：板的简单剪切
图6给出了板的简单剪切变换。2 cm×2 cm的正的剪切平滑，变形至A′B′C′D′。
方形ABCD发生45°
表4 算例2中各增量步板应力解析解 Table 4 Analytical solution to stresses of plate at some
step in example 2 Pa
n
(n)
. 1
(n)
.2 2
(n)
. 2
2
′
4 6
1.923 08×107 1.923 08×107 7.692 31×108 1.154 33×108 1.154 33×108 1.542 31×109 83×108 2.326 94×109 2.891 83×108 2.891
5.418 30×108 5.418 30×108 3.130 90×109 8.755 06×108 8.755 06×108 3.962 14×109
810
图6 算例2中板的简单剪切变换 Fig.6 Simple shear of plate in example 2
本算例中板的简单剪切变形的最终应变、应力数值计算结果分别如图7，8所示。由此可知，计算结果与理论分析的解析解完全吻合。
取初始拖带系与固定直线直角坐标系一致，变形体在平衡状态Ω(n)的固定坐标为
(0)
x1(n)=x1(0)+x2K(n)⎫⎪
⎬ (36) (n)(0)
⎪x2=x2⎭
(10)2(10)
(a) S.和S.2 1
采用和分解弹性平面应力增量程序进行计算，通过理论分析，可以得到第n个增量步的板的应变增量计算公式为
(n)2(n)
∆.1=∆.2=01(n)2(n)∆.1=∆.12
⎫⎪
(37) ⎬
=∆K/2⎪⎭
(10)
(b) . 2
板的应力增量计算公式为
(n)2(n)(n−1)
∆.=∆.2=∆K.12(n)2(n)
∆.=∆.12
图7 算例2中板的最终应变数值计算结果
Fig.7 Numerical results of final strains of plate in example 2
⎫⎪
⎬ (38) (n−1)
=(G+.1/2)∆K⎪⎭
取N=10个增量步计算，网格为20×20，各增量步板的应变和应力解析解分别如表3，4所示。
表3 算例2中各增量步板的应变解析解
Table 3 Analytical solution to strains of plate at some steps
in example 2
n 2 4 6 8
(n)
.1 1
2(n)
.2
(n)
. 2
(10)(10)
(a) .和.2 12
0 0 0 0
0 0 0 0
0.10 0.20 0.30 0.40
(10)
(b) . 2
图8 算例2中板的最终应力数值计算结果(单位：Pa) Fig.8 Numerical results of final stresses of plate in example
2(unit：Pa)
10 0 0 0.50
第26卷 第5期 何满潮，等. 深部软岩工程大变形力学分析设计系统 • 941 •
5.3 算例3：采场顶板垮落大变形
采场顶板垮落问题可以简化为端部作用有集中荷载的板的弯曲变形问题。如图9所示的悬臂板，其端部作用有集中力P，尺寸为：L = 25 m，H= 2 m。材料参数为：E=1.0×10
10
(90)
(a) . 1
Pa，ν=0.3。
(90)(b) .2 2
图9 算例3中端部作用有集中力的悬臂板的弯曲大变形 Fig.9 Bending of a cantilever plate with a concentrated force
at the free end in example 3
采用和分解弹性平面应力增量程序进行计算，网格为10×1，取增量步为90，每步增量荷载为∆P=0.8×106 N。本算例中板的最终应变和应力数
值计算结果分别如图10，11所示。由此可知，板的固定端处，剪应变、水平方向正应变较大。由于板的固定端弯矩很大，所以水平方向正应力很大，呈线性分布，在板的上下边界达到最大。
(90)
(c) . 2
图11 算例3中板的最终应力计算结果(单位：Pa) Fig.11 Calculating results of final stresses of plate in example
3(unit：Pa)
(90)
(a) . 1
图12 算例4中端部作用有弯矩的悬臂梁弯曲大变形 Fig.12 Bending of a cantilever beam with a moment at the
free end in example 4
E = 1.0×1010 Pa，ν = 0.3。
采用逐级更新坐标法的梁单元程序，计算在弯
(b)
2(90)
.2
矩M作用下梁的弯曲变形。网格为200×1，取13个增量步分析，每步的弯距增量为∆M=1.0×107
N·m。图13给出了算例4中梁的各增量步变形图。
5
6
7
4
3
2
(90)
(c) . 2
891011
1213
图10 算例3中板的最终应变计算结果
Fig.10 Calculating results of final strains of plate in example 3
1
5.4 算例4：梁的弯曲大变形
如图12所示的悬臂梁，其端部承受的弯矩为M，
图13 算例4中梁的各增量步变形图
Fig.13 Deformation of beam at every step in example 4
梁的尺寸为：L = 25 m，H = 2 m。材料参数为：
• 942 • 岩石力学与工程学报 2007年
从图13中可以看出，随着弯矩的增加，梁逐渐弯曲直至最终变成一个圆形。
表5 岩体计算参数
Table 5 Calculating parameters of rock masses
岩体 岩体1 岩体2 岩体3 岩体4，5
E/(104 MPa)
容重/(kN·m3)
－
6 工程实例
本文分别采用软件的大变形和小变形弹性增量平面应变计算模块，对柳海运输大巷分步开挖的软岩变形进行数值模拟。该矿软岩的蒙脱石含量为
泊松比ν
1.000 16.62 0.15～0.20 1.000 28.57 0.15～0.20
0.050 12.80 0.25 0.005 22.94 0.25
90%～96%，其大巷底臌大变形破坏照片见图14。巷道断面为直墙半圆拱型，圆顶半径为1.910 m，直墙高和宽均为1.810 m，计算区域的高和宽均为
30 m。分5步开挖，其有限元计算模型和网格见 图15，其中，1～10表示材料号，且材料号1～5分别对应于岩体1～5，材料号6～10表示分5步开挖的大巷岩体材料。边界条件为：底边和侧边简支，顶边施加上部围岩自重作用的应力12 MPa(相当于埋深500 m)。运输大巷岩体计算参数见表5，通过对岩体3，4降低弹性模量来考虑软岩遇水软化效应。
图16 柳海运输大巷的有限元网格大变形模型变形图 Fig.16 Deformed meshes of transport tunnel in large
deformation model of Liuhai Mine
图14 柳海运输大巷的底臌大变形破坏照片 Fig.14 A picture for large deformation of a transport tunnel
in Liuhai Mine
图17 柳海运输大巷的有限元网格小变形模型变形图 Fig.17 Deformed meshes of transport tunnel in small
deformation model of Liuhai Mine
6
1724589103
计算模型分析所得巷道变形量对比情况见表6。由此可以看出，大变形模型能更好地模拟巷道由于底板软岩遇水软化引起的底膨大变形。
表6 大变形和小变形计算模型分析所得巷道变形量对比
图15 柳海运输大巷的有限元计算模型和网格
Table 6 Comparison of deformation of tunnel by two models
with large deformation and small deformation respectively
计算模型 大变形 小变形
底膨变形量/m 1.118 32 0.091 47
两帮移近量/m 0.011 404 0.043 216
Fig.15 Finite element model and meshes of transport tunnel in
Liuhai Mine
两种计算模型所得柳海运输大巷的有限元网格
第26卷 第5期 何满潮，等. 深部软岩工程大变形力学分析设计系统 • 943 •
7 结 论
深部软岩工程大变形力学分析设计系统是旨在建立考虑深部三相(气、液、固)和三场(应力场、渗流场、温度场)耦合作用的非线性力学数值模拟系统，应用于研究深部资源开采面临的非线性岩石力学问题。系统具有如下的4个特色：(1) 以和分解有限变形力学分析模型为主体，包括极分解力学分析模型供用户进行对比分析；(2) 强调非线性大变形的力学设计特色：力学对策设计、过程设计和参数设计；(3) 包括多种岩土材料本构模型和工程材料单元，可以进行边坡工程、地基工程、地下洞室工程等岩土工程的大变形力学分析；(4) 具有自主知识产权。
本文简要地介绍了现阶段推出的软件版本
[4] 柴肇云，康天合，李义宝. 物化型软岩微结构单元特征及其胀缩性
研究[J]. 岩石力学与工程学报，2006，25(6)：1 265–1 269.(CHAI Zhaoyun，KANG Tianhe，LI Yibao. Study on microstructure unit character and its swell-shrink property for physico-chemical soft rock[J]. Chinese Journal of Rock Mechanics and Engineering，2006，25(6)：1 265–1 269.(in Chinese))
[5] TRUESDELL C，NOLL W. Nonlinear field theories of mechanics[M].
[S. l. ]：[s. n. ]，1965.
[6] BIOT M A. Mechanics of incremental deformations[M]. New York：
John Wiley and Sons，Inc.，1965.
[7] 郭仲衡. 非线性弹性理论[M]. 北京：科学出版社，1980.(GUO
Zhongheng. Nonlinear elasticity mechanics[M]. Beijing：Science Press，1980.(in Chinese))
[8] 陈至达. 连续介质有限变形力学几何场论[J]. 力学学报，1979，
11(2)：107–117.(CHEN Zhida. Geometric theory of finite deformation mechanics for continuum[J]. Acta Mechanica Sinica，1979，11(2)：107–117.(in Chinese))
[9] 陈至达. 理性力学[M]. 重庆：重庆大学出版社，2000.(CHEN Zhida.
Rational mechanics[M]. Chongqing：Chongqing University Press，2000.(in Chinese))
[10] 何满潮. 露天矿高边坡工程[M]. 北京：煤炭工业出版社，1991.(HE
Manchao. High slope engineering of open pit mine[M]. Beijing：China Coal Industry Publishing House，1991.(in Chinese))
[11] WANG C，CHEN Z D. Micro rotation analysis of material cracking
and toughness[J]. Int. J. Fracture，1992，54(4)：359–369.
[12] QIN Z，CHEN Z D. Large deformation analysis of shells with finite
element method based on the S-R decomposition theorem[J]. Computer and Structures，1988，30(4)：957–961.
[13] 李 平. 非线性大变形有限元分析的更新拖带坐标法及其应用[博
士学位论文][D]. 北京：中国矿业大学，1991.(LI Ping. The updated co-moving coordinate formulation for the nonlinear large deformation finite element analysis and application[Ph. D. Thesis][D]. Beijing：China University of Mining and Technology，1991.(in Chinese)) [14] 何满潮. 软岩巷道工程概论[M]. 徐州：中国矿业大学出版，1993.
(HE Manchao. General theory of soft rock tunnel engineering[M]. Xuzhou：China University of Mining and Technology Press，1993.(in Chinese))
[15] KYUICHIRO W S. Variational methods in elasticity and plasiticity[M].
2nd ed.[S. l. ]：Pergamon Press，1975.
LDEAS1.0的和分解有限变形力学基本方程、软件功能和开发平台。典型算例的数值分析结果说明了该软件系统程序的正确性。对柳海软岩巷道开挖过程分别采用大变形和小变形弹性增量模块进行对比分析，计算结果表明，大变形模型分析能更好地模拟巷道底臌大变形。 参考文献(References)：
[1] 赵长海，周小兵，贺建国，等. 极软岩隧洞的设计与施工[J]. 岩石
力学与工程学报，2006，25(增1)：3 034–3 039.(ZHAO Changhai，ZHOU Xiaobing，HE Jianguo，et al. Design and construction of tunnel in extremely soft rock[J]. Chinese Journal of Rock Mechanics and Engineering，2006，25(Supp.1)：3 034–3 039.(in Chinese)) [2] 郭志飚，胡永光，任爱武，等. 深部膨胀性软岩巷道修复技术研究[J].
采矿与安全工程学报，2006，23(3)：316–319.(GUO Zhibiao，HU Yongguang，REN Aiwu，et al. Roadway repair technology in swell soft rock with high depth and high stress[J]. Journal of Mining and Safety Engineering，2006，23(3)：316–319.(in Chinese)) [3]
万援朝. 二次支护原理在深井软岩硐室支护中的实践[J]. 煤炭科学技术，2006，34(9)：5–7.(WAN Yuanchao. Practice of second support theory applied to deep mine soft rock section chamber
support[J]. Coal Science and Technology，2006，34(9)：5–7.(in Chinese))
第26卷 第5期
岩石力学与工程学报 Vol.26 No.5
2007年5月 Chinese Journal of Rock Mechanics and Engineering May，2007
深部软岩工程大变形力学分析设计系统
何满潮1，陈 新1，梁国平2，钱华山2，周永发2，庄小燕2
(1. 中国矿业大学 力学与建筑工程学院，北京 100083；
2. 中国科学院 数学与系统科学研究院飞箭软件有限公司，北京 100098)
摘要：采用飞箭软件公司的有限元程序自动生成系统作为软件开发平台，进行“深部软岩工程大变形力学分析设计系统”的合作开发工作。该软件系统具有如下特色：(1) 以和分解有限变形力学分析模型为主体，包括极分解力学分析模型供用户进行对比分析；(2) 强调非线性大变形力学设计特色：力学对策设计、过程设计和参数设计；(3) 包括多种岩土材料本构模型和工程材料单元，可以进行边坡工程、地基工程、地下洞室工程等大变形力学分析；(4) 具有我国自主知识产权。目前该软件的弹性、弹塑性大变形二维计算程序和界面已经开发完成。通过对典型算例的弹性数值分析，验证大变形程序的正确性。该软件系统的问世，可为深部软岩巷道等岩体大变形力学问题提供具有分析和设计功能的数值分析软件系统。
关键词：岩石力学；和分解；大变形；有限元软件系统；深部软岩工程
中图分类号：TU 452 文献标识码：A 文章编号：1000–6915(2007)05–0934–10
SOFTWARE SYSTEM FOR LARGE DEFORMATION MECHANICAL ANALYSIS OF SOFT ROCK ENGINEERING AT GREAT DEPTH
HE Manchao1，CHEN Xin1，LIANG Guoping2，QIAN Huashan2，ZHOU Yongfa2，ZHUANG Xiaoyan2
(1. School of Mechanics and Civil Engineering，China University of Mining and Technology，Beijing 100083，China；2. Beijing FEGEN Software Co.，Ltd.，Academy of Mathematics and Systems Science，Chinese Academy of Sciences，Beijing 100098，China)
Abstract：Up to now，there is no finite element software including nonlinear large deformation theory based on S-R decomposition theorem，in which the deformation gradient F is decomposed into the addition of the deformation and the rotation by a co-moving coordinate system method. In this work，such a software system waw developed，which is named as A Software on Large Deformation Analysis of Soft Rock Engineering at Great Depth. Two work groups，worked together on the software by using finite element program generator(FEPG) developed by FEGEN Software Company. The software has the following features：(1) it includes large deformation theory based on S-R decomposition theorem first time as well as one based on polar decomposition theorem for comparison，the latter is included by most commercial finite element software；(2) three design methods for nonlinear mechanical problem，i.e. countermeasure design，procedure design and parameter design，are fulfilled here for large deformation analyses；(3) it includes generally used elastoplastic constitutive models for geomaterials and element types for supports，and therefore is fit for large deformation mechanical analyses for slope，foundation and underground engineering；and (4) it has a copyright of our own country. At present，the
收稿日期：2006–09–26；修回日期：2006–12–18
基金项目：国家重点基础研究发展计划(973)项目(2006CB202200)；国家自然科学基金重大项目(50490273)；国家基金委创新群体基金资助项目(50221402)；教育部科学技术研究重大项目(10405)
作者简介：何满潮(1956–)，男，博士，1981年毕业于长春地质学院工程地质专业，现任教授、博士生导师，主要从事软岩、边坡及工程地质、地热
第26卷 第5期 何满潮，等. 深部软岩工程大变形力学分析设计系统 • 935 •
interface of the software has accomplished and programs of elastic and elastoplastic constitutive models based on two large deformation theories in two-dimension have been fulfilled，tested and implemented into the interface. Representative numerical examples confirm the rationality of the programs. The software will be a useful numerical tool for large deformation analyses and nonlinear mechanical design for soft rock engineering at great depth.
Key words：rock mechanics；S-R decomposition；large deformation；finite element software system；soft rock engineering at great depth
授经过多年研究，在拖带系描述法、张量分析和微分变换群等数学理论的基础上，于1979年提出了应变和转动的和分解定理，完善了Stokes在1945年提出的固体位移场分解理论，克服了以往大变形理论的缺陷，构成了非线性连续体力学的新体系。十几年来，以拖带系的和分解定理为基础，许多学者[10
～14]
1 引 言
随着开采深度的增加，岩爆、瓦斯突出、流变、底板突水等非线性动力学灾害现象日趋增多，严重影响了深部资源的安全高效开采。深部开采工程中产生的岩石力学问题已经成为目前国内外采矿及岩石力学界研究的焦点。
深部软岩巷道工程在巷道开挖后，巷道周围岩体的变形一般都较大，出现巷道冒顶、底膨和侧胀等现象都是围岩大变形的结果[1
～4]
研究了拖带系中的应变、转动、应力、应力
客观率的度量和变分原理，形成了基于和分解的有限变形理论。该理论的应用涉及岩石、混凝土、金属等刚体的大变形、非线性效应、断裂力学、细观力学、生物力学、分形几何、边坡稳定、软岩巷道支护、地质构造大变形、实验力学的激光频谱分析法以及物性方程等诸多领域，取得了令人瞩目的成果。
在深部软岩工程中，由于涉及到物理非线性、
～9]
。为了合理地进
行岩体变形分析，除了要很好地建立符合岩体变形的物理、力学特性本构模型外，还必须采用能正确描述大变形的非线性几何理论，才能得到较为合理的结果。
随着非线性几何场论的研究[5
不断取得进
展，目前形成了两大有限变形理论：(1) 采用固定坐标系描述方法，以Green非线性应变作为应变度量和Finger极分解定理得到的转动张量为转动度量的经典有限变形理论；(2) 采用拖带坐标系描述方法，基于应变和转动的和分解定理(additive decomposition theorem，简称S-R分解定理)的有限变形理论。
Green非线性应变张量克服了大转动时Cauchy应变张量在物体发生有限转动时对变形度量的失真，给出了大变形的确定性度量，因而被许多非线性力学著作采用。而Green应变张量没有相应的转动张量与之匹配，对全面研究大变形、大转动问题是一个重大缺陷。极分解定理定义的转动张量虽然在一定程度上弥补了Green应变张量没有相匹配的转动张量的缺点，然而它包含左、右极分解2个分解式，变形和转动的分解有先后之分，分解式不唯一，分解得到的伸长张量不适合作为变形的度量。
针对经典有限变形理论的这些问题，陈至达教
几何非线性和接触边界非线性等力学问题，因此其理论解的求解在数学上遇到非常大的困难，需要借助于有限元、有限差分和离散元等数值方法和软件。
目前，商业有限元软件的有限变形力学模型中，无一例外的都是采用Green非线性应变和Finger极分解定理来度量应变和转动，没有提供基于和分解的有限变形力学模型，这是很不全面的。为此，重大项目首席科学家何满潮教授组织了中国矿业大学(北京)力学与建筑工程学院和中国科学院数学与系统科学研究院飞箭软件公司2个单位的研究力量，策划和实施了“深部软岩工程大变形力学分析设计系统”的合作开发工作。该软件系统以软岩大变形理论和非线性力学设计方法[14]为理论基础，通过有限元程序自动生成平台FEPG5.2，建立深部三相(气、液、固)和三场(应力场、渗流场、温度场)耦合作用的数值模拟系统，应用于理论科研与工程实践，这对研究深部资源开采面临的岩石力学问题具有重大现实意义。
• 936 • 岩石力学与工程学报 2007年
2 和分解有限变形力学基本方程
2.1 位形变换的和分解定理
和分解有限变形力学采用嵌含在可变形体中的拖带系描述方法，在变形体的运动过程中，变形体中各质点的拖带坐标不变，而度量拖带坐标的基矢量则不断发生变化。
考察变形体的运动，设在初始平衡状态t0时刻，变形体中某质点在固定参考系中的空间位置为P，矢径为r，拖带坐标为xi，基矢量为gi。在变形后的某个平衡状态t时刻，该点运动到P′的空间位置，矢径变为R(R为正交变换，表现为点集的转动)，基矢量变为gi。变形前、后质点的协变基矢量分别定义为
∂r∆r⎫
lim=
⎪∂xi∆xi→0∆xi⎪
⎬ (1)
∂R∆R⎪gi=i=lim
∆xi→0∆xi⎪∂x⎭gi=
式中：S为对称变换，表现为点集的形变。
变换系数F的S-R分解分量形式为
⎫
⎪⎪1
Sij=(ui|j+ui|Tj)−(1−cosΘ)L⋅ki⋅L⋅kj⋅⎬ (7)
2⎪
⎪Rij=δij+L⋅ij⋅sinΘ+(1−cosΘ)L⋅ki⋅L⋅kj⋅⎭式中：Θ，L⋅ij⋅分别为平均整旋角和转轴方位张量分量，可分别定义如下：
Fji=Sij+Rij
Θ=±arcsin(−ωijωij/2)2⎫⎪
L=ωi/(sinΘ)
⋅i
j⋅
j
⎬ (8) ⎪⎭
式中：ωij为位移梯度的反对称部分，可表示为
ωij=(ui|j−ui|Tj)/2 (9)
位移梯度的物理分量为
ui|j=ui
j
giijj (10)
变形前、后质点的矢径与拖带坐标分别存在如下微分关系：
∂ri0⎫dr=idx=gidxi⎪
∂x⎪
⎬ (2)
∂Ri
dR=idx=gidxi⎪
⎪∂x⎭
应变张量、转角与转角方位的物理分量为 ⎫1
Sij=(ui|j+ui|Tj)−(1−cosΘ)L⋅ki⋅L⋅kj⋅⎪
2
⎪
⎪⎛1 i j⎞
⎪Θ=±arcsin⎜−ωjωi⎟
2⎪⎝⎠ ⎬ (11)
⎪1 i
L⋅ij⋅=ωj⎪
sinΘ⎪
⎪ i1 i iT
ωj=(u|j−u|j)⎪
2⎭2.2 广义Euler运动学公式
质点的位移矢量为
u=R−r (3)
变形前、后质点的协变基矢量之间存在如下变换关系：
gi=Fijgj，Fij=δij+uj|i (4)
在节2.1考察变形体的位形变换时，采用了初始拖带系基矢量gi为基准；对于非线性大变形运动，采用瞬时拖带系为参考系可更为准确地进行描述。
考察变形体从t时刻到t+∆t时刻的位形变换，拖带系的基矢量变换为：gi=gi+∆gi。
基矢量的变化率为
∂gi∆g
=limi=Vj||igj (12) ∂t∆t→0∆t
式中：δij为Kronecker符号；uj|i为质点位移在初始拖带系中的协变导数，可表示为
∂uj0jl
u|i=i+Γilu (5)
∂x
j
式中：0Γilj为初始拖带系的Christoffel符号。
因此，变形体的运动变换可以用质点局部嵌含坐标系的基矢量变换(式(4))来表示。
和分解定理指出，若变换系数F在形变体点集域内是单值连续的，处处具有一阶导数，则此运动变换可以唯一分解为正交与对称2个子变换的直和，即
式中：Vj||i为变形体中一点的速度在瞬时拖带系中的协变导数，可表示为
∂Vj
V||i=+VkΓkij (13) i
j
第26卷 第5期 何满潮，等. 深部软岩工程大变形力学分析设计系统 • 937 •
式中：Γkij为t时刻瞬时拖带系的Christoffel符号。
陈至达将Euler运动学公式推广到了变形体的运动，得到了变形体运动的速度梯度和分解式：
为σij。对几何非线性程度不高的弹性变形体，可采用如下全量弹性应力–应变关系：
..lk
σij=E⋅ijk⋅Sl (23)
j+R j (14) Vj||i=Sii
..l
式中：E⋅ijk⋅为弹性刚度矩阵。
..l
对各向同性材料，E⋅ijk⋅可表示为
..lililil
E⋅ijk⋅=λδjδk+µ(ggjk+δkδj) (24)
j，R j分别为该点的应变速率和局部转动式中：Sii
速度，可分别表示为
⎫ i=1(Vi||+Vi||T)Sjjj⎪⎪2
⎬ (15)
1 i=LiΘ =(Vi||−Vi||T)⎪Rjjjj⎪2⎭
2.3 应力、应力速率、平衡方程
在瞬时拖带系中，定义g(j)为单位化的逆变基矢量，g(i)为与g(j)满足对偶关系g(i)⋅g(j)=δij的协变基矢量：
gj=gjjg(j)，gi=g(i)
gii (16)
式中：λ，µ均为材料的Lame常数，可由实验测定。此时的物性关系为
σij=λSkkδ.ij+2µSij (25)
对几何非线性程度较高的弹性体或弹塑性体，变形较大，微元体将发生畸变。全量形式的本构定律没有考虑位形改变的影响，严格来说是不精确的，必须采用如下应力客观率形式本构关系：
∆i
..lk
σj=⋅ijk⋅Sl (26)
Euler应力张量可以在瞬时拖带系中分解为
..l
式中：⋅ijk⋅为材料的切线刚度矩阵。
σ=σ.ijgigj=σ.jg(i)g(j) (17)
i
利用式(21)可得到应力率形式的本构关系：
..l ki ki k
ij=⋅ijkσ⋅Sl+σkSj−σjSk (27)
式中：σ.ij为Euler应力张量的物理分量，可表示为
在实际计算中，应力、应变要采用物理分量。例如，对各向同性弹性变形体的平面应变问题，增量应力–应变关系为
σ.ij=σ.ijgjjgii (18)
应力的物质导数为
dσ.ijdt
ij−σ.ikνk||j−σ.ijνk||k (19) =σ
应力的物质导数不是客观的时间导数，与刚体的转动有关，去掉转动项得到客观应力速率σ⋅j为
∆i
1 1 2 1 2⎫
)∆S2+σ2∆S1∆σ1=(λ+2µ)∆S11+(λ−σ1
2 2 1⎪ 2 2 1⎪
∆σ2=(λ−σ2)∆S1+(λ+2µ)∆S2+σ1∆S2⎬ (28)
⎪ 1 1 1 1 1
∆σ2=−σ2∆S1+(2µ+σ1)∆S2⎪⎭2.5 变分原理
在变形体发生大变形时，外力与位移呈非线性关系，且边界随时间变动，合理的能量原理应以瞬时位形的功率形式表达，下面应用势能率原理来建立增量虚位移的变分方程。
在瞬时拖带系中，势能率原理为：对一切满足速度、角速度与应变速率谐调条件和速度边界条件的所有可能形变状态，满足平衡方程以及力边界条件的真实应力状态，将使其能量泛函J的变化率取驻值：
ii jdΩ− =σiδSδJfδVdΩ−TδVda=0 jiii∫∫∫
Ω
Ω
Sσ
即
σ⋅j=
∆i
dσ
i
.j
dt
k (20) +σ.ikRj
k+σiS k ij−σ.kjSσ.j=σj.jk (21)
∆i
平衡方程是建立在变形后的位形上的，可表示为
σ.ij||i+fj=0 (22)
式中：σ.ij||i为应力在瞬时拖带系中的协变导数，且
i
σ.ij||i=σ.ij，i+Γikσ.kj−Γijkσ.kj；fj为体积力。
2.4 物性关系
(29)
式中：为表面力。
• 938 • 岩石力学与工程学报 2007年
如果表面力、体积力不随时间变化，应力与荷载相关但不为时间函数，势能驻值可以表示为如下的全量形式的变分原理：
δJ=∫σijδSijdΩ−∫fiδuidΩ−∫Tiδuida=0 (30)
Ω
Ω
Sσ
和后处理子系统。通过界面将这些子系统组合起来。
系统的主流程和主界面分别如图1，2所示。
严格地说，只有率形式的变分原理是普遍适用的。但在小位移条件下，全量变分原理仍有相当的准确性。
3 基于FEPG有限变形力学程序开发
在深部软岩工程大变形力学分析设计系统中，所有的有限元程序都是利用其自动生成系统
(FEPG)开发的。
有限元程序自动生成系统(FEPG)由北京飞箭软件有限公司的创始人之一、中国科学院数学与系统研究院博士生导师梁国平研究员历经十余年的潜心研究而独创，该系统的主要设计思想是采用元件化的程序设计方法和有限元语言，根据有限元方法统一的数学原理及其内在规律，以类似于数学公式推理的方式，由有限元问题的微分方程表达式及其求解算法自动产生有限元程序。该系统适用于各种领域的各种工程与科学计算问题，突破了目前通用有限元程序只适用于特定领域和特定问题的限制。用户只需输入有限元方法所需的各种表达式和公式，即可自动产生所需的全部有限元程序，免去了大量繁琐的有限元编程劳动，并保证了程序的正确性和统一性。
在系统的有限元程序中，极分解增量分析模块采用了更新拉格朗日法[15]，和分解增量分析模块采用了更新拖带坐标法[13]，位移或位移增量的求解采用牛顿迭代法。
图1 系统的主流程图
Fig.1 Main flowchart of software system
图2 系统的主界面 Fig.2 Main interface of software
4 系统的功能
深部软岩工程大变形力学分析设计系统的最终目标是建立考虑深部三相(气、液、固)和三场(应力场、渗流场、温度场)耦合作用的非线性力学数值模拟系统，来研究深部资源开采面临的非线性岩石力学问题。
该系统将分阶段、分层次开发完成。第1阶段，该系统称为“深部工程软岩大变形力学分析设计系统LDEAS1.0版”。LDEAS1.0分4个功能独立的子系该版本具有二维固体应力场的非线性力学数值模拟分析功能，其有限变形力学分析模型包括极分解、和分解2个模块，计算模型包括弹性全量、弹性增量、弹塑性增量的平面应力以及平面应变问题共计16个分析计算程序。单元模型有实体单元、Goodman节理单元、殷有泉弹塑性节理单元和刚架、锚杆、桁架等结构单元，可以考虑桁架、锚杆、锚索支护和开挖以及建造施工过程。该系统采用力学对策设计、过程设计和参数设计的非线性大变形力学设计方法。
5 典型算例
本节采用板的有限转动与伸长、板的简单剪切、采场顶板垮落大变形和梁的弯曲大变形4个典型算例说明系统中各程序的正确性。
在增量分析中，将固体问题的加载途径划分成(0)(1)(Ν−1)(Ν)
第26卷 第5期 何满潮，等. 深部软岩工程大变形力学分析设计系统 • 939 •
其中，Ω(0)为初始状态，Ω(n)为其中任意一个中间状态，Ω(Ν)为最后的平衡状态。对第n(n=1， 2， "，N)个增量步，平衡状态由Ω(n−1)变换到Ω(n)。 5.1 算例1：板的有限转动与伸长
图3给出了板的有限转动与伸长变换。2 cm×2 cm的正方形板由0123→01′2′3′，整体转动45°，同时均等扩张为32cm×32cm的正方形。材料参数：E = 1.0×1010 Pa，ν = 0.3。
x22′
x2
1′
g2
g2
.ij(n)=.ij(n−1)+∆.ij(n)⎫⎪
⎬ (35) i(n)i(n−1)i(n)
.j=.j+∆.j⎪⎭
取N=10个增量步计算，网格为20×20，各增量步板的应变和应力解析解分别如表1，2所示。
表1 算例1中各增量步板应变解析解
Table 1 Analytical solution to strains of plate at some
steps in example 1
n 2 4
(n). 1
2(n)
.2
(n)
.1 2
0.116 667 0.275 557 0.441 243 0.598 166 0.741 759
0.116 667 0.275 557 0.441 243 0.598 166 0.741 759
0 0 0 0 0
x1
3′
g1g1
6 8 10
x1
表2 算例1中各增量步板应力解析解
Table 2 Analytical solution to stresses of plate at some steps
in example 1 Pa
n 2 4 6 810
(n)
. 1
(n)
.2 2
(n)
. 2
图3 算例1中板的有限转动与伸长变换 Fig.3 Rotation and expansion of a plate in example 1
1.619 048×109 1.619 048×109 3.551 792×109 3.551 792×109 5.256 592×109 5.256 592×109 6.617 906×109 6.617 906×109 448×109 7.679 448×109 7.679
0 0 0 0 0
取初始拖带系与固定直线直角坐标系一致，变
形体在平衡状态Ω(n)的固定坐标为
(0)
⎫x1(n)=[1+0.5K(n)]x1(0)−1.5K(n)x2
(31) (n)(0)(0)⎬x2=1.5K(n)x1+[1+0.5K(n)]x2⎭
本算例中板的有限转动与伸长变换的最终应变
其中，
和应力数值计算结果分别如图4，5所示。由此可知，
K(n)=n/N
采用和分解弹性平面应力增量程序进行计算，通过理论分析，可以得到第n个增量步的板应变增量计算公式为
(n)2(n)
∆.=∆.2=α(K(n−1))∆K⎫1⎪
(32) ⎬ (n)2(n)
⎪∆.2=∆.1=0⎭
(10)2(10)(10)
(a) .和.2 (b) .1 12
图4 算例1中板的最终应变数值计算结果
Fig.4 Numerical results of final strains of plate in example 1
其中，
α(K)=
1+5K
(33) 2
2(1+K+2.5K)
板的应力增量计算公式为
⎛E(n)2(n)(n−1)⎞1(n)⎫∆S∆.1=∆=−⎜⎟1.2.1.1⎪⎝1−ν⎠⎬ (34)
⎪(n)2(n)(n−1)
∆.1=∆.1=.∆K22⎭
(a)
(10)
.1
(10)
和.22
(10)
.2
(b)
图5 算例1中板的最终应力数值计算结果(单位：Pa) Fig.5 Numerical results of final stresses of plate in example
1(unit：Pa)
平衡状态(n)
• 940 • 岩石力学与工程学报 2007年
计算结果与理论分析的解析解完全吻合。 5.2 算例2：板的简单剪切
图6给出了板的简单剪切变换。2 cm×2 cm的正的剪切平滑，变形至A′B′C′D′。
方形ABCD发生45°
表4 算例2中各增量步板应力解析解 Table 4 Analytical solution to stresses of plate at some
step in example 2 Pa
n
(n)
. 1
(n)
.2 2
(n)
. 2
2
′
4 6
1.923 08×107 1.923 08×107 7.692 31×108 1.154 33×108 1.154 33×108 1.542 31×109 83×108 2.326 94×109 2.891 83×108 2.891
5.418 30×108 5.418 30×108 3.130 90×109 8.755 06×108 8.755 06×108 3.962 14×109
810
图6 算例2中板的简单剪切变换 Fig.6 Simple shear of plate in example 2
本算例中板的简单剪切变形的最终应变、应力数值计算结果分别如图7，8所示。由此可知，计算结果与理论分析的解析解完全吻合。
取初始拖带系与固定直线直角坐标系一致，变形体在平衡状态Ω(n)的固定坐标为
(0)
x1(n)=x1(0)+x2K(n)⎫⎪
⎬ (36) (n)(0)
⎪x2=x2⎭
(10)2(10)
(a) S.和S.2 1
采用和分解弹性平面应力增量程序进行计算，通过理论分析，可以得到第n个增量步的板的应变增量计算公式为
(n)2(n)
∆.1=∆.2=01(n)2(n)∆.1=∆.12
⎫⎪
(37) ⎬
=∆K/2⎪⎭
(10)
(b) . 2
板的应力增量计算公式为
(n)2(n)(n−1)
∆.=∆.2=∆K.12(n)2(n)
∆.=∆.12
图7 算例2中板的最终应变数值计算结果
Fig.7 Numerical results of final strains of plate in example 2
⎫⎪
⎬ (38) (n−1)
=(G+.1/2)∆K⎪⎭
取N=10个增量步计算，网格为20×20，各增量步板的应变和应力解析解分别如表3，4所示。
表3 算例2中各增量步板的应变解析解
Table 3 Analytical solution to strains of plate at some steps
in example 2
n 2 4 6 8
(n)
.1 1
2(n)
.2
(n)
. 2
(10)(10)
(a) .和.2 12
0 0 0 0
0 0 0 0
0.10 0.20 0.30 0.40
(10)
(b) . 2
图8 算例2中板的最终应力数值计算结果(单位：Pa) Fig.8 Numerical results of final stresses of plate in example
2(unit：Pa)
10 0 0 0.50
第26卷 第5期 何满潮，等. 深部软岩工程大变形力学分析设计系统 • 941 •
5.3 算例3：采场顶板垮落大变形
采场顶板垮落问题可以简化为端部作用有集中荷载的板的弯曲变形问题。如图9所示的悬臂板，其端部作用有集中力P，尺寸为：L = 25 m，H= 2 m。材料参数为：E=1.0×10
10
(90)
(a) . 1
Pa，ν=0.3。
(90)(b) .2 2
图9 算例3中端部作用有集中力的悬臂板的弯曲大变形 Fig.9 Bending of a cantilever plate with a concentrated force
at the free end in example 3
采用和分解弹性平面应力增量程序进行计算，网格为10×1，取增量步为90，每步增量荷载为∆P=0.8×106 N。本算例中板的最终应变和应力数
值计算结果分别如图10，11所示。由此可知，板的固定端处，剪应变、水平方向正应变较大。由于板的固定端弯矩很大，所以水平方向正应力很大，呈线性分布，在板的上下边界达到最大。
(90)
(c) . 2
图11 算例3中板的最终应力计算结果(单位：Pa) Fig.11 Calculating results of final stresses of plate in example
3(unit：Pa)
(90)
(a) . 1
图12 算例4中端部作用有弯矩的悬臂梁弯曲大变形 Fig.12 Bending of a cantilever beam with a moment at the
free end in example 4
E = 1.0×1010 Pa，ν = 0.3。
采用逐级更新坐标法的梁单元程序，计算在弯
(b)
2(90)
.2
矩M作用下梁的弯曲变形。网格为200×1，取13个增量步分析，每步的弯距增量为∆M=1.0×107
N·m。图13给出了算例4中梁的各增量步变形图。
5
6
7
4
3
2
(90)
(c) . 2
891011
1213
图10 算例3中板的最终应变计算结果
Fig.10 Calculating results of final strains of plate in example 3
1
5.4 算例4：梁的弯曲大变形
如图12所示的悬臂梁，其端部承受的弯矩为M，
图13 算例4中梁的各增量步变形图
Fig.13 Deformation of beam at every step in example 4
梁的尺寸为：L = 25 m，H = 2 m。材料参数为：
• 942 • 岩石力学与工程学报 2007年
从图13中可以看出，随着弯矩的增加，梁逐渐弯曲直至最终变成一个圆形。
表5 岩体计算参数
Table 5 Calculating parameters of rock masses
岩体 岩体1 岩体2 岩体3 岩体4，5
E/(104 MPa)
容重/(kN·m3)
－
6 工程实例
本文分别采用软件的大变形和小变形弹性增量平面应变计算模块，对柳海运输大巷分步开挖的软岩变形进行数值模拟。该矿软岩的蒙脱石含量为
泊松比ν
1.000 16.62 0.15～0.20 1.000 28.57 0.15～0.20
0.050 12.80 0.25 0.005 22.94 0.25
90%～96%，其大巷底臌大变形破坏照片见图14。巷道断面为直墙半圆拱型，圆顶半径为1.910 m，直墙高和宽均为1.810 m，计算区域的高和宽均为
30 m。分5步开挖，其有限元计算模型和网格见 图15，其中，1～10表示材料号，且材料号1～5分别对应于岩体1～5，材料号6～10表示分5步开挖的大巷岩体材料。边界条件为：底边和侧边简支，顶边施加上部围岩自重作用的应力12 MPa(相当于埋深500 m)。运输大巷岩体计算参数见表5，通过对岩体3，4降低弹性模量来考虑软岩遇水软化效应。
图16 柳海运输大巷的有限元网格大变形模型变形图 Fig.16 Deformed meshes of transport tunnel in large
deformation model of Liuhai Mine
图14 柳海运输大巷的底臌大变形破坏照片 Fig.14 A picture for large deformation of a transport tunnel
in Liuhai Mine
图17 柳海运输大巷的有限元网格小变形模型变形图 Fig.17 Deformed meshes of transport tunnel in small
deformation model of Liuhai Mine
6
1724589103
计算模型分析所得巷道变形量对比情况见表6。由此可以看出，大变形模型能更好地模拟巷道由于底板软岩遇水软化引起的底膨大变形。
表6 大变形和小变形计算模型分析所得巷道变形量对比
图15 柳海运输大巷的有限元计算模型和网格
Table 6 Comparison of deformation of tunnel by two models
with large deformation and small deformation respectively
计算模型 大变形 小变形
底膨变形量/m 1.118 32 0.091 47
两帮移近量/m 0.011 404 0.043 216
Fig.15 Finite element model and meshes of transport tunnel in
Liuhai Mine
两种计算模型所得柳海运输大巷的有限元网格
第26卷 第5期 何满潮，等. 深部软岩工程大变形力学分析设计系统 • 943 •
7 结 论
深部软岩工程大变形力学分析设计系统是旨在建立考虑深部三相(气、液、固)和三场(应力场、渗流场、温度场)耦合作用的非线性力学数值模拟系统，应用于研究深部资源开采面临的非线性岩石力学问题。系统具有如下的4个特色：(1) 以和分解有限变形力学分析模型为主体，包括极分解力学分析模型供用户进行对比分析；(2) 强调非线性大变形的力学设计特色：力学对策设计、过程设计和参数设计；(3) 包括多种岩土材料本构模型和工程材料单元，可以进行边坡工程、地基工程、地下洞室工程等岩土工程的大变形力学分析；(4) 具有自主知识产权。
本文简要地介绍了现阶段推出的软件版本
[4] 柴肇云，康天合，李义宝. 物化型软岩微结构单元特征及其胀缩性
研究[J]. 岩石力学与工程学报，2006，25(6)：1 265–1 269.(CHAI Zhaoyun，KANG Tianhe，LI Yibao. Study on microstructure unit character and its swell-shrink property for physico-chemical soft rock[J]. Chinese Journal of Rock Mechanics and Engineering，2006，25(6)：1 265–1 269.(in Chinese))
[5] TRUESDELL C，NOLL W. Nonlinear field theories of mechanics[M].
[S. l. ]：[s. n. ]，1965.
[6] BIOT M A. Mechanics of incremental deformations[M]. New York：
John Wiley and Sons，Inc.，1965.
[7] 郭仲衡. 非线性弹性理论[M]. 北京：科学出版社，1980.(GUO
Zhongheng. Nonlinear elasticity mechanics[M]. Beijing：Science Press，1980.(in Chinese))
[8] 陈至达. 连续介质有限变形力学几何场论[J]. 力学学报，1979，
11(2)：107–117.(CHEN Zhida. Geometric theory of finite deformation mechanics for continuum[J]. Acta Mechanica Sinica，1979，11(2)：107–117.(in Chinese))
[9] 陈至达. 理性力学[M]. 重庆：重庆大学出版社，2000.(CHEN Zhida.
Rational mechanics[M]. Chongqing：Chongqing University Press，2000.(in Chinese))
[10] 何满潮. 露天矿高边坡工程[M]. 北京：煤炭工业出版社，1991.(HE
Manchao. High slope engineering of open pit mine[M]. Beijing：China Coal Industry Publishing House，1991.(in Chinese))
[11] WANG C，CHEN Z D. Micro rotation analysis of material cracking
and toughness[J]. Int. J. Fracture，1992，54(4)：359–369.
[12] QIN Z，CHEN Z D. Large deformation analysis of shells with finite
element method based on the S-R decomposition theorem[J]. Computer and Structures，1988，30(4)：957–961.
[13] 李 平. 非线性大变形有限元分析的更新拖带坐标法及其应用[博
士学位论文][D]. 北京：中国矿业大学，1991.(LI Ping. The updated co-moving coordinate formulation for the nonlinear large deformation finite element analysis and application[Ph. D. Thesis][D]. Beijing：China University of Mining and Technology，1991.(in Chinese)) [14] 何满潮. 软岩巷道工程概论[M]. 徐州：中国矿业大学出版，1993.
(HE Manchao. General theory of soft rock tunnel engineering[M]. Xuzhou：China University of Mining and Technology Press，1993.(in Chinese))
[15] KYUICHIRO W S. Variational methods in elasticity and plasiticity[M].
2nd ed.[S. l. ]：Pergamon Press，1975.
LDEAS1.0的和分解有限变形力学基本方程、软件功能和开发平台。典型算例的数值分析结果说明了该软件系统程序的正确性。对柳海软岩巷道开挖过程分别采用大变形和小变形弹性增量模块进行对比分析，计算结果表明，大变形模型分析能更好地模拟巷道底臌大变形。 参考文献(References)：
[1] 赵长海，周小兵，贺建国，等. 极软岩隧洞的设计与施工[J]. 岩石
力学与工程学报，2006，25(增1)：3 034–3 039.(ZHAO Changhai，ZHOU Xiaobing，HE Jianguo，et al. Design and construction of tunnel in extremely soft rock[J]. Chinese Journal of Rock Mechanics and Engineering，2006，25(Supp.1)：3 034–3 039.(in Chinese)) [2] 郭志飚，胡永光，任爱武，等. 深部膨胀性软岩巷道修复技术研究[J].
采矿与安全工程学报，2006，23(3)：316–319.(GUO Zhibiao，HU Yongguang，REN Aiwu，et al. Roadway repair technology in swell soft rock with high depth and high stress[J]. Journal of Mining and Safety Engineering，2006，23(3)：316–319.(in Chinese)) [3]
万援朝. 二次支护原理在深井软岩硐室支护中的实践[J]. 煤炭科学技术，2006，34(9)：5–7.(WAN Yuanchao. Practice of second support theory applied to deep mine soft rock section chamber
support[J]. Coal Science and Technology，2006，34(9)：5–7.(in Chinese))