深部软岩工程大变形力学分析设计系统何满潮

第26卷第5期岩石力学与工程学报Vol.26No.52007年5月ChineseJournalofRockMechanicsandEngineeringMay，2007收稿日期：2006–09–26；修回日期：2006–12–18基金项目：国家重点基础研究发展(973)项目(2006CB202200)；国家自然科学基金重大项目(50490273)；国家基金委创新群体基金资助项目(50221402)；教育部科学技术研究重大项目(10405)作者简介：何满潮(1956–)，男，博士，1981年毕业于长春地质学院工程地质专业，现任教授、博士生导师，主要从事软岩、边坡及工程地质、地热工程等方面的教学与研究工作。E-mail：hemanchao@163.com深部软岩工程大变形力学设计系统何满潮1，陈新1，梁国平2，钱华山2，周永发2，庄小燕2(1.中国矿业大学力学与建筑工程学院，北京100083；2.中国科学院数学与系统科学研究院飞箭软件有限公司，北京100098)摘要：采用飞箭软件公司的有限元程序自动生成系统作为软件开发平台，进行“深部软岩工程大变形力学分析设计系统”的合作开发工作。该软件系统具有如下特色：(1)以和分解有限变形力学分析模型为主体，包括极分解力学分析模型供用户进行对比分析；(2)强调非线性大变形力学设计特色：力学对策设计、过程设计和参数设计；(3)包括多种岩土材料本构模型和工程材料单元，可以进行边坡工程、地基工程、地下洞室工程等大变形力学分析；(4)具有我国自主知识产权。目前该软件的弹性、弹塑性大变形二维计算程序和界面已经开发完成。通过对典型算例的弹性数值分析，验证大变形程序的正确性。该软件系统的问世，可为深部软岩巷道等岩体大变形力学问题提供具有分析和设计功能的数值分析软件系统。关键词：岩石力学；和分解；大变形；有限元软件系统；深部软岩工程中图分类号：TU452文献标识码：A文章编号：1000–6915(2007)05–0934–10SOFTWARESYSTEMFORLARGEDEFORMATIONMECHANICALANALYSISOFSOFTROCKENGINEERINGATGREATDEPTHHEManchao1，CHENXin1，LIANGGuoping2，QIANHuashan2，ZHOUYongfa2，ZHUANGXiaoyan2(1.SchoolofMechanicsandCivilEngineering，ChinaUniversityofMiningandTechnology，Beijing100083，China；2.BeijingFEGENSoftwareCo.，Ltd.，AcademyofMathematicsandSystemsScience，ChineseAcademyofSciences，Beijing100098，China)Abstract：Uptonow，thereisnofiniteelementsoftwareincludingnonlinearlargedeformationtheorybasedonS-Rdecompositiontheorem，inwhichthedeformationgradientFisdecomposedintotheadditionofthedeformationandtherotationbyaco-movingcoordinatesystemmethod.Inthiswork，suchasoftwaresystemwawdeveloped，whichisnamedasASoftwareonLargeDeformationAnalysisofSoftRockEngineeringatGreatDepth.Twoworkgroups，workedtogetheronthesoftwarebyusingfiniteelementprogramgenerator(FEPG)developedbyFEGENSoftwareCompany.Thesoftwarehasthefollowingfeatures：(1)itincludeslargedeformationtheorybasedonS-Rdecompositiontheoremfirsttimeaswellasonebasedonpolardecompositiontheoremforcomparison，thelatterisincludedbymostcommercialfiniteelementsoftware；(2)threedesignmethodsfornonlinearmechanicalproblem，i.e.countermeasuredesign，proceduredesignandparameterdesign，arefulfilledhereforlargedeformationanalyses；(3)itincludesgenerallyusedelastoplasticconstitutivemodelsforgeomaterialsandelementtypesforsupports，andthereforeisfitforlargedeformationmechanicalanalysesforslope，foundationandundergroundengineering；and(4)ithasacopyrightofourowncountry.Atpresent，the第26卷第5期何满潮，等.深部软岩工程大变形力学分析设计系统•935•interfaceofthesoftwarehasaccomplishedandprogramsofelasticandelastoplasticconstitutivemodelsbasedontwolargedeformationtheoriesintwo-dimensionhavebeenfulfilled，testedandimplementedintotheinterface.Representativenumericalexamplesconfirmtherationalityoftheprograms.Thesoftwarewillbeausefulnumericaltoolforlargedeformationanalysesandnonlinearmechanicaldesignforsoftrockengineeringatgreatdepth.Keywords：rockmechanics；S-Rdecomposition；largedeformation；finiteelementsoftwaresystem；softrockengineeringatgreatdepth1引言随着开采深度的增加，岩爆、瓦斯突出、流变、底板突水等非线性动力学灾害现象日趋增多，严重影响了深部资源的安全高效开采。深部开采工程中产生的岩石力学问题已经成为目前国内外采矿及岩石力学界研究的焦点。深部软岩巷道工程在巷道开挖后，巷道周围岩体的变形一般都较大，出现巷道冒顶、底膨和侧胀等现象都是围岩大变形的结果[1～4]。为了合理地进行岩体变形分析，除了要很好地建立符合岩体变形的物理、力学特性本构模型外，还必须采用能正确描述大变形的非线性几何理论，才能得到较为合理的结果。随着非线性几何场论的研究[5～9]不断取得进展，目前形成了两大有限变形理论：(1)采用固定坐标系描述方法，以Green非线性应变作为应变度量和Finger极分解定理得到的转动张量为转动度量的经典有限变形理论；(2)采用拖带坐标系描述方法，基于应变和转动的和分解定理(additivedecompositiontheorem，简称S-R分解定理)的有限变形理论。Green非线性应变张量克服了大转动时Cauchy应变张量在物体发生有限转动时对变形度量的失真，给出了大变形的确定性度量，因而被许多非线性力学著作采用。而Green应变张量没有相应的转动张量与之匹配，对全面研究大变形、大转动问题是一个重大缺陷。极分解定理定义的转动张量虽然在一定程度上弥补了Green应变张量没有相匹配的转动张量的缺点，然而它包含左、右极分解2个分解式，变形和转动的分解有先后之分，分解式不唯一，分解得到的伸长张量不适合作为变形的度量。针对经典有限变形理论的这些问题，陈至达教授经过多年研究，在拖带系描述法、张量分析和微分变换群等数学理论的基础上，于1979年提出了应变和转动的和分解定理，完善了Stokes在1945年提出的固体位移场分解理论，克服了以往大变形理论的缺陷，构成了非线性连续体力学的新体系。十几年来，以拖带系的和分解定理为基础，许多学者[10～14]研究了拖带系中的应变、转动、应力、应力客观率的度量和变分原理，形成了基于和分解的有限变形理论。该理论的应用涉及岩石、混凝土、金属等刚体的大变形、非线性效应、断裂力学、细观力学、生物力学、分形几何、边坡稳定、软岩巷道支护、地质构造大变形、实验力学的激光频谱分析法以及物性方程等诸多领域，取得了令人瞩目的成果。在深部软岩工程中，由于涉及到物理非线性、几何非线性和接触边界非线性等力学问题，因此其理论解的求解在数学上遇到非常大的困难，需要借助于有限元、有限差分和离散元等数值方法和软件。目前，商业有限元软件的有限变形力学模型中，无一例外的都是采用Green非线性应变和Finger极分解定理来度量应变和转动，没有提供基于和分解的有限变形力学模型，这是很不全面的。为此，重大项目首席科学家何满潮教授组织了中国矿业大学(北京)力学与建筑工程学院和中国科学院数学与系统科学研究院飞箭软件公司2个单位的研究力量，策划和实施了“深部软岩工程大变形力学分析设计系统”的合作开发工作。该软件系统以软岩大变形理论和非线性力学设计方法[14]为理论基础，通过有限元程序自动生成平台FEPG5.2，建立深部三相(气、液、固)和三场(应力场、渗流场、温度场)耦合作用的数值模拟系统，应用于理论科研与工程实践，这对研究深部资源开采面临的岩石力学问题具有重大现实意义。•936•岩石力学与工程学报2007年2和分解有限变形力学基本方程2.1位形变换的和分解定理和分解有限变形力学采用嵌含在可变形体中的拖带系描述方法，在变形体的运动过程中，变形体中各质点的拖带坐标不变，而度量拖带坐标的基矢量则不断发生变化。考察变形体的运动，设在初始平衡状态0t时刻，变形体中某质点在固定参考系中的空间位置为P，矢径为r，拖带坐标为ix，基矢量为i0g。在变形后的某个平衡状态t时刻，该点运动到P′的空间位置，矢径变为R(R为正交变换，表现为点集的转动)，基矢量变为ig。变形前、后质点的协变基矢量分别定义为⎪⎪⎭⎪⎪⎬⎫∆∆=∂∂=∆∆=∂∂=→∆→∆ixiiixiixxxxiiRRgrrg000limlim(1)变形前、后质点的矢径与拖带坐标分别存在如下微分关系：⎪⎪⎭⎪⎪⎬⎫=∂∂==∂∂=iiiiiiiixxxxxxdddddd0gRRgrr(2)质点的位移矢量为rRu−=(3)变形前、后质点的协变基矢量之间存在如下变换关系：ijjijijjiiuFF|0+==δ，gg(4)式中：jiδ为Kronecker符号；iju|为质点位移在初始拖带系中的协变导数，可表示为ljilijijuxuuΓ0|+∂∂=(5)式中：jilΓ0为初始拖带系的Christoffel符号。因此，变形体的运动变换可以用质点局部嵌含坐标系的基矢量变换(式(4))来表示。和分解定理指出，若变换系数F在形变体点集域内是单值连续的，处处具有一阶导数，则此运动变换可以唯一分解为正交与对称2个子变换的直和，即RSF+=(6)式中：S为对称变换，表现为点集的形变。变换系数F的S-R分解分量形式为⎪⎪⎭⎪⎪⎬⎫−++=−−+=+=⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅kjikijjijikjikTjijiijijijijLLLRLLuuSRSF)cos1(sin)cos1()||(21ΘΘδΘ(7)式中：Θ，ijL⋅⋅分别为平均整旋角和转轴方位张量分量，可分别定义如下：⎪⎭⎪⎬⎫=−±=⋅⋅)/(sin)2/arcsin(21ΘωωωΘjiijjiijL(8)式中：ijω为位移梯度的反对称部分，可表示为2/)||(Tjijiijuu−=ω(9)位移梯度的物理分量为jjiijijigguu0|=(10)应变张量、转角与转角方位的物理分量为⎪⎪⎪⎪⎪⎭⎪⎪⎪⎪⎪⎬⎫−==⎟⎠⎞⎜⎝⎛−±=−−+=⋅⋅⋅⋅⋅⋅)||(21sin121arcsin)cos1()||(2121TjijiijijijjiijkjikTjijiijuuLLLuuSωωΘωωΘΘ(11)2.2广义Euler运动学在节2.1考察变形体的位形变换时，采用了初始拖带系基矢量ig为基准；对于非线性大变形运动，采用瞬时拖带系为参考系可更为准确地进行描述。考察变形体从t时刻到tt∆+时刻的位形变换，拖带系的基矢量变换为：iiiggg∆+=0。基矢量的变化率为jijitiVttggg||lim0=∆∆=∂∂→∆(12)式中：ijV||为变形体中一点的速度在瞬时拖带系中的协变导数，可表示为jkikijijVxVVΓ+∂∂=||(13)第26卷第5期何满潮，等.深部软岩工程大变形力学分析设计系统•937•式中：jkiΓ为t时刻瞬时拖带系的Christoffel符号。陈至达将Euler运动学公式推广到了变形体的运动，得到了变形体运动的速度梯度和分解式：jijiijRSV+=||(14)式中：jiS，jiR分别为该点的应变速率和局部转动速度，可分别表示为⎪⎪⎭⎪⎪⎬⎫−==+=)||||(21)||||(21TjijiijijTjijiijVVLRVVSΘ(15)2.3应力、应力速率、平衡方程在瞬时拖带系中，定义)(jg为单位化的逆变基矢量，)(ig为与)(jg满足对偶关系jijδ=⋅)()(ggi的协变基矢量：iiiijjjjgg)()(gggg==，(16)Euler应力张量可以在瞬时拖带系中分解为)()(..jiijjiijggggσσ==σ(17)式中：ij.σ为Euler应力张量的物理分量，可表示为iijjijijgg..σσ=(18)应力的物质导数为kkijjkikijijt||||dd...νσνσσσ−−=(19)应力的物质导数不是客观的时间导数，与刚体的转动有关，去掉转动项得到客观应力速率ij⋅∆σ为kjikijijRt..ddσσσ+=⋅∆(20)即kkijkjjkijijSS...σσσσ+−=∆(21)平衡方程是建立在变形后的位形上的，可表示为0||.=+jiijfσ(22)式中：iij||.σ为应力在瞬时拖带系中的协变导数，且jkkijkjiikiijiij....||σΓσΓσσ−+=，；jf为体积力。2.4物性关系当变形体内没有体矩时，应力是对称张量，记为ijσ。对几何非线性程度不高的弹性变形体，可采用如下全量弹性应力–应变关系：kllijkijSE..⋅⋅=σ(23)式中：lijkE..⋅⋅为弹性刚度矩阵。对各向同性材料，lijkE..⋅⋅可表示为)(..ljikjkillkijlijkggEδδµδλδ++=⋅⋅(24)式中：λ，µ均为材料的Lame常数，可由实验测定。此时的物性关系为ijijkkijSSµδλσ2.+=(25)对几何非线性程度较高的弹性体或弹塑性体，变形较大，微元体将发生畸变。全量形式的本构定律没有考虑位形改变的影响，严格来说是不精确的，必须采用如下应力客观率形式本构关系：kllijkijSE..⋅⋅∆=σ(26)式中：lijkE..⋅⋅为材料的切线刚度矩阵。利用式(21)可得到应力率形式的本构关系：kkijkjikkllijkijSSSEσσσ−+=⋅⋅..(27)在实际计算中，应力、应变要采用物理分量。例如，对各向同性弹性变形体的平面应变问题，增量应力–应变关系为⎪⎪⎭⎪⎪⎬⎫∆++∆−=∆∆+∆++∆−=∆∆+∆−+∆+=∆1211111212122122112222211222111111)2()2()()()2(SSSSSSSSσµσσσµλσλσσσλµλσ(28)2.5变分原理在变形体发生大变形时，外力与位移呈非线性关系，且边界随时间变动，合理的能量原理应以瞬时位形的功率形式表达，下面应用势能率原理来建立增量虚位移的变分方程。在瞬时拖带系中，势能率原理为：对一切满足速度、角速度与应变速率谐调条件和速度边界条件的所有可能形变状态，满足平衡方程以及力边界条件的真实应力状态，将使其能量泛函J的变化率取驻值：0ddd=δ−δ−δ=δ∫∫∫aVTVfSJiSiiijiijσΩΩΩΩσ(29)式中：iT为表面力。•938•岩石力学与工程学报2007年如果表面力、体积力不随时间变化，应力与荷载相关但不为时间函数，势能驻值可以表示为如下的全量形式的变分原理：0ddd=δ−δ−δ=δ∫∫∫auTufSJiSiiijiijσΩΩΩΩσ(30)严格地说，只有率形式的变分原理是普遍适用的。但在小位移条件下，全量变分原理仍有相当的准确性。3基于FEPG有限变形力学程序开发在深部软岩工程大变形力学分析设计系统中，所有的有限元程序都是利用其自动生成系统(FEPG)开发的。有限元程序自动生成系统(FEPG)由北京飞箭软件有限公司的创始人之一、中国科学院数学与系统研究院博士生导师梁国平研究员历经十余年的潜心研究而独创，该系统的主要设计思想是采用元件化的程序设计方法和有限元语言，根据有限元方法统一的数学原理及其内在规律，以类似于数学公式推理的方式，由有限元问题的微分方程表达式及其求解算法自动产生有限元程序。该系统适用于各种领域的各种工程与科学计算问题，突破了目前通用有限元程序只适用于特定领域和特定问题的限制。用户只需输入有限元方法所需的各种表达式和公式，即可自动产生所需的全部有限元程序，免去了大量繁琐的有限元编程劳动，并保证了程序的正确性和统一性。在系统的有限元程序中，极分解增量分析模块采用了更新拉格朗日法[15]，和分解增量分析模块采用了更新拖带坐标法[13]，位移或位移增量的求解采用牛顿迭代法。4系统的功能深部软岩工程大变形力学分析设计系统的最终目标是建立考虑深部三相(气、液、固)和三场(应力场、渗流场、温度场)耦合作用的非线性力学数值模拟系统，来研究深部资源开采面临的非线性岩石力学问题。该系统将分阶段、分层次开发完成。第1阶段，该系统称为“深部工程软岩大变形力学分析设计系统LDEAS1.0版”。LDEAS1.0分4个功能独立的子系统：前处理子系统、参数设置子系统、计算子系统和后处理子系统。通过界面将这些子系统组合起来。系统的主流程和主界面分别如图1，2所示。图1系统的主流程图Fig.1Mainflowchartofsoftwaresystem图2系统的主界面Fig.2Maininterfaceofsoftware该版本具有二维固体应力场的非线性力学数值模拟分析功能，其有限变形力学分析模型包括极分解、和分解2个模块，计算模型包括弹性全量、弹性增量、弹塑性增量的平面应力以及平面应变问题共计16个分析计算程序。单元模型有实体单元、Goodman节理单元、殷有泉弹塑性节理单元和刚架、锚杆、桁架等结构单元，可以考虑桁架、锚杆、锚索支护和开挖以及建造施工过程。该系统采用力学对策设计、过程设计和参数设计的非线性大变形力学设计方法。5典型算例本节采用板的有限转动与伸长、板的简单剪切、采场顶板垮落大变形和梁的弯曲大变形4个典型算例系统中各程序的正确性。在增量分析中，将固体问题的加载途径划分成若干段平衡状态：)()1()1()0(ΝΝΩΩΩΩ，，，，−。运行计算保存退出建立新项目开始参数设置(线荷载、工程材料、单元材料、地应力等信息)GID后处理(位移、应力、应变等计算结果的图形显示)GID前处理(几何建模、设定材料号和边界条件、单元剖分)第26卷第5期何满潮，等.深部软岩工程大变形力学分析设计系统•939•其中，)0(Ω为初始状态，)(nΩ为其中任意一个中间状态，)(ΝΩ为最后的平衡状态。对第1(，=nn)2N，，个增量步，平衡状态由)1(−nΩ变换到)(nΩ。5.1算例1：板的有限转动与伸长图3给出了板的有限转动与伸长变换。2cm×2cm的正方形板由0123→01′2′3′，整体转动45°，同时均等扩张为23cm×23cm的正方形。材料参数：E=1.0×1010Pa，ν=0.3。图3算例1中板的有限转动与伸长变换Fig.3Rotationandexpansionofaplateinexample1取初始拖带系与固定直线直角坐标系一致，变形体在平衡状态)(nΩ的固定坐标为⎭⎬⎫++=−+=)0(2)0(1)(2)0(2)0(1)(1)](5.01[)(5.1)(5.1)](5.01[xnKxnKxxnKxnKxnn(31)其中，NnnK/)(=采用和分解弹性平面应力增量程序进行计算，通过理论分析，可以得到第n个增量步的板应变增量计算公式为⎪⎭⎪⎬⎫=∆=∆∆=∆=∆−0)()(21.)(12.)1()(22.)(11.nnnnnSSKKSSα(32)其中，)5.21(251)(2KKKK+++=α(33)板的应力增量计算公式为⎪⎭⎪⎬⎫∆=∆=∆∆⎟⎠⎞⎜⎝⎛−−=∆=∆−−KSEnnnnnnn)1(12.)(21.)(12.)(11.)1(11.)(22.)(11.1σσσσνσσ(34)平衡状态)(nΩ板的应变、应力计算公式分别为⎪⎭⎪⎬⎫∆+=∆+=−−)(.)1(.)(.)(.)1(.)(.nijnijnijnijnijnijSSSσσσ(35)取10=N个增量步计算，网格为2020×，各增量步板的应变和应力解析解分别如表1，2所示。表1算例1中各增量步板应变解析解Table1Analyticalsolutiontostrainsofplateatsomestepsinexample1n)(11.nS)(22.nS)(12.nS20.1166670.116667040.2755570.275557060.4412430.441243080.5981660.5981660100.7417590.7417590表2算例1中各增量步板应力解析解Table2Analyticalsolutiontostressesofplateatsomestepsinexample1Pan)(11.nσ)(22.nσ)(12.nσ21.619048×1091.619048×109043.551792×1093.551792×109065.256592×1095.256592×109086.617906×1096.617906×1090107.679448×1097.679448×1090本算例中板的有限转动与伸长变换的最终应变和应力数值计算结果分别如图4，5所示。由此可知，(a))10(11.S和)10(22.S(b))10(12.S图4算例1中板的最终应变数值计算结果Fig.4Numericalresultsoffinalstrainsofplateinexample1(a))10(11.σ和)10(22.σ(b))10(12.σ图5算例1中板的最终应力数值计算结果(单位：Pa)Fig.5Numericalresultsoffinalstressesofplateinexample1(unit：Pa)g20g10g1g23′2′1′x1x1x2x2•940•岩石力学与工程学报2007年计算结果与理论分析的解析解完全吻合。5.2算例2：板的简单剪切图6给出了板的简单剪切变换。2cm×2cm的正方形ABCD发生45°的剪切平滑，变形至A′B′C′D′。图6算例2中板的简单剪切变换Fig.6Simpleshearofplateinexample2取初始拖带系与固定直线直角坐标系一致，变形体在平衡状态)(nΩ的固定坐标为⎪⎭⎪⎬⎫=+=)0(2)(2)0(2)0(1)(1)(xxnKxxxnn(36)采用和分解弹性平面应力增量程序进行计算，通过理论分析，可以得到第n个增量步的板的应变增量计算公式为⎪⎭⎪⎬⎫∆=∆=∆=∆=∆2/0)(21.)(12.)(22.)(11.KSSSSnnnn(37)板的应力增量计算公式为⎪⎭⎪⎬⎫∆+=∆=∆∆=∆=∆−−KGKnnnnnn)2/(2)1(11.)(21.)(12.)1(12.)(22.)(11.σσσσσσ(38)取10=N个增量步计算，网格为2020×，各增量步板的应变和应力解析解分别如表3，4所示。表3算例2中各增量步板的应变解析解Table3Analyticalsolutiontostrainsofplateatsomestepsinexample2n)(11.nS)(22.nS)(12.nS2000.104000.206000.308000.4010000.50表4算例2中各增量步板应力解析解Table4Analyticalsolutiontostressesofplateatsomestepinexample2Pan)(11.nσ)(22.nσ)(12.nσ21.92308×1071.92308×1077.69231×10841.15433×1081.15433×1081.54231×10962.89183×1082.89183×1082.32694×10985.41830×1085.41830×1083.13090×109108.75506×1088.75506×1083.96214×109本算例中板的简单剪切变形的最终应变、应力数值计算结果分别如图7，8所示。由此可知，计算结果与理论分析的解析解完全吻合。(a))10(11.S和)10(22.S(b))10(12.S图7算例2中板的最终应变数值计算结果Fig.7Numericalresultsoffinalstrainsofplateinexample2(a))10(11.σ和)10(22.σ(b))10(12.σ图8算例2中板的最终应力数值计算结果(单位：Pa)Fig.8Numericalresultsoffinalstressesofplateinexample2(unit：Pa)g2γg20ADC(－1，1)(－1，－1)(1，－1)B(1，1)A′D′C′x2x2B′x10g1，g1第26卷第5期何满潮，等.深部软岩工程大变形力学分析设计系统•941•5.3算例3：采场顶板垮落大变形采场顶板垮落问题可以简化为端部作用有集中荷载的板的弯曲变形问题。如图9所示的悬臂板，其端部作用有集中力P，尺寸为：L=25m，H=2m。材料参数为：10100.1×=EPa，3.0=ν。图9算例3中端部作用有集中力的悬臂板的弯曲大变形Fig.9Bendingofacantileverplatewithaconcentratedforceatthefreeendinexample3采用和分解弹性平面应力增量程序进行计算，网格为110×，取增量步为90，每步增量荷载为6108.0×=∆PN。本算例中板的最终应变和应力数值计算结果分别如图10，11所示。由此可知，板的固定端处，剪应变、水平方向正应变较大。由于板的固定端弯矩很大，所以水平方向正应力很大，呈线性分布，在板的上下边界达到最大。(a))90(11.S(b))90(22.S(c))90(12.S图10算例3中板的最终应变计算结果Fig.10Calculatingresultsoffinalstrainsofplateinexample35.4算例4：梁的弯曲大变形如图12所示的悬臂梁，其端部承受的弯矩为M，梁的尺寸为：L=25m，H=2m。材料参数为：(a))90(11.σ(b))90(22.σ(c))90(12.σ图11算例3中板的最终应力计算结果(单位：Pa)Fig.11Calculatingresultsoffinalstressesofplateinexample3(unit：Pa)图12算例4中端部作用有弯矩的悬臂梁弯曲大变形Fig.12Bendingofacantileverbeamwithamomentatthefreeendinexample4E=1.0×1010Pa，ν=0.3。采用逐级更新坐标法的梁单元程序，计算在弯矩M作用下梁的弯曲变形。网格为200×1，取13个增量步分析，每步的弯距增量为=∆M1.0×107N·m。图13给出了算例4中梁的各增量步变形图。图13算例4中梁的各增量步变形图Fig.13Deformationofbeamateverystepinexample4LHPAB321045678910111213•942•岩石力学与工程学报2007年从图13中可以看出，随着弯矩的增加，梁逐渐弯曲直至最终变成一个圆形。6工程实例本文分别采用软件的大变形和小变形弹性增量平面应变计算模块，对柳海运输大巷分步开挖的软岩变形进行数值模拟。该矿软岩的蒙脱石含量为90%～96%，其大巷底臌大变形破坏照片见图14。巷道断面为直墙半圆拱型，圆顶半径为1.910m，直墙高和宽均为1.810m，计算区域的高和宽均为30m。分5步开挖，其有限元计算模型和网格见图15，其中，1～10表示材料号，且材料号1～5分别对应于岩体1～5，材料号6～10表示分5步开挖的大巷岩体材料。边界条件为：底边和侧边简支，顶边施加上部围岩自重作用的应力12MPa(相当于埋深500m)。运输大巷岩体计算参数见表5，通过对岩体3，4降低弹性模量来考虑软岩遇水软化效应。图14柳海运输大巷的底臌大变形破坏照片Fig.14ApictureforlargedeformationofatransporttunnelinLiuhaiMine图15柳海运输大巷的有限元计算模型和网格Fig.15FiniteelementmodelandmeshesoftransporttunnelinLiuhaiMine两种计算模型所得柳海运输大巷的有限元网格大、小变形模型变形图分别如图16，17所示；两种表5岩体计算参数Table5Calculatingparametersofrockmasses岩体E/(104MPa)容重/(kN·m－3)泊松比ν岩体11.00016.620.15～0.20岩体21.00028.570.15～0.20岩体30.05012.800.25岩体4，50.00522.940.25图16柳海运输大巷的有限元网格大变形模型变形图Fig.16DeformedmeshesoftransporttunnelinlargedeformationmodelofLiuhaiMine图17柳海运输大巷的有限元网格小变形模型变形图Fig.17DeformedmeshesoftransporttunnelinsmalldeformationmodelofLiuhaiMine计算模型分析所得巷道变形量对比情况见表6。由此可以看出，大变形模型能更好地模拟巷道由于底板软岩遇水软化引起的底膨大变形。表6大变形和小变形计算模型分析所得巷道变形量对比Table6Comparisonofdeformationoftunnelbytwomodelswithlargedeformationandsmalldeformationrespectively计算模型底膨变形量/m两帮移近量/m大变形1.118320.011404小变形0.091470.04321661724589103第26卷第5期何满潮，等.深部软岩工程大变形力学分析设计系统•943•7结论深部软岩工程大变形力学分析设计系统是旨在建立考虑深部三相(气、液、固)和三场(应力场、渗流场、温度场)耦合作用的非线性力学数值模拟系统，应用于研究深部资源开采面临的非线性岩石力学问题。系统具有如下的4个特色：(1)以和分解有限变形力学分析模型为主体，包括极分解力学分析模型供用户进行对比分析；(2)强调非线性大变形的力学设计特色：力学对策设计、过程设计和参数设计；(3)包括多种岩土材料本构模型和工程材料单元，可以进行边坡工程、地基工程、地下洞室工程等岩土工程的大变形力学分析；(4)具有自主知识产权。本文简要地介绍了现阶段推出的软件版本LDEAS1.0的和分解有限变形力学基本方程、软件功能和开发平台。典型算例的数值分析结果说明了该软件系统程序的正确性。对柳海软岩巷道开挖过程分别采用大变形和小变形弹性增量模块进行对比分析，计算结果表明，大变形模型分析能更好地模拟巷道底臌大变形。参考文献(References)：[1]赵长海，周小兵，贺建国，等.极软岩隧洞的设计与施工[J].岩石力学与工程学报，2006，25(增1)：3034–3039.(ZHAOChanghai，ZHOUXiaobing，HEJianguo，etal.Designandconstructionoftunnelinextremelysoftrock[J].ChineseJournalofRockMechanicsandEngineering，2006，25(Supp.1)：3034–3039.(inChinese))[2]郭志飚，胡永光，任爱武，等.深部膨胀性软岩巷道修复技术研究[J].采矿与安全工程学报，2006，23(3)：316–319.(GUOZhibiao，HUYongguang，RENAiwu，etal.Roadwayrepairtechnologyinswellsoftrockwithhighdepthandhighstress[J].JournalofMiningandSafetyEngineering，2006，23(3)：316–319.(inChinese))[3]万援朝.二次支护原理在深井软岩硐室支护中的实践[J].煤炭科学技术，2006，34(9)：5–7.(WANYuanchao.Practiceofsecondsupporttheoryappliedtodeepminesoftrocksectionchambersupport[J].CoalScienceandTechnology，2006，34(9)：5–7.(inChinese))[4]柴肇云，康天合，李义宝.物化型软岩微结构单元特征及其胀缩性研究[J].岩石力学与工程学报，2006，25(6)：1265–1269.(CHAIZhaoyun，KANGTianhe，LIYibao.Studyonmicrostructureunitcharacteranditsswell-shrinkpropertyforphysico-chemicalsoftrock[J].ChineseJournalofRockMechanicsandEngineering，2006，25(6)：1265–1269.(inChinese))[5]TRUESDELLC，NOLLW.Nonlinearfieldtheoriesofmechanics[M].[S.l.]：[s.n.]，1965.[6]BIOTMA.Mechanicsofincrementaldeformations[M].NewYork：JohnWileyandSons，Inc.，1965.[7]郭仲衡.非线性弹性理论[M].北京：科学出版社，1980.(GUOZhongheng.Nonlinearelasticitymechanics[M].Beijing：SciencePress，1980.(inChinese))[8]陈至达.连续介质有限变形力学几何场论[J].力学学报，1979，11(2)：107–117.(CHENZhida.Geometrictheoryoffinitedeformationmechanicsforcontinuum[J].ActaMechanicaSinica，1979，11(2)：107–117.(inChinese))[9]陈至达.理性力学[M].重庆：重庆大学出版社，2000.(CHENZhida.Rationalmechanics[M].Chongqing：ChongqingUniversityPress，2000.(inChinese))[10]何满潮.露天矿高边坡工程[M].北京：煤炭工业出版社，1991.(HEManchao.Highslopeengineeringofopenpitmine[M].Beijing：ChinaCoalIndustryPublishingHouse，1991.(inChinese))[11]WANGC，CHENZD.Microrotationanalysisofmaterialcrackingandtoughness[J].Int.J.Fracture，1992，54(4)：359–369.[12]QINZ，CHENZD.LargedeformationanalysisofshellswithfiniteelementmethodbasedontheS-Rdecompositiontheorem[J].ComputerandStructures，1988，30(4)：957–961.[13]李平.非线性大变形有限元分析的更新拖带坐标法及其应用[博士学位][D].北京：中国矿业大学，1991.(LIPing.Theupdatedco-movingcoordinateformulationforthenonlinearlargedeformationfiniteelementanalysisandapplication[Ph.D.Thesis][D].Beijing：ChinaUniversityofMiningandTechnology，1991.(inChinese))[14]何满潮.软岩巷道工程概论[M].徐州：中国矿业大学出版，1993.(HEManchao.Generaltheoryofsoftrocktunnelengineering[M].Xuzhou：ChinaUniversityofMiningandTechnologyPress，1993.(inChinese))[15]KYUICHIROWS.Variationalmethodsinelasticityandplasiticity[M].2nded.[S.l.]：PergamonPress，1975.