岩土工程数值分析方法

1、第7章 岩土工程数值分析方法 有限元法 边界元法 有限差分法 离散单元法 1 有限元法 概述基本思路：将复杂的结构看成由有限个仅在结点处联结的整体，首先对每一个单元分析其特性，建立相关物理量之间的相互联系。然后，依据单元之间的联系再将各单元组装成整体，从而获得整体特性方程，应用方程相应的解法，即可完成整个问题的分析分析过程：结构离散化确定单元位移模式单元特性分析集成总体特性接方程求未知量 概述工程界比较流行、被广泛使用的大型有限元软件：MSC(航空航天领域，是目前规模最大的有限元分析系统)；Marc(非线性分析软件)；Adina(可进行结构、流体、热的耦合计算，具有隐式和显式两种时间积分算法，

2、非线性运算功能强大)；ANSYS（可进行结构、流体、热、电磁场的计算，是有限元分析的通用软件包）优点：可分析几何形状及受荷条件复杂、非均质的各种实际结构；可在计算中模拟各种复杂的材料本构关系、边界条件等；前后处理技术先进。 有限单元法的理论基础虚位移原理：受给定外力的变形体处于平衡状态的充要条件是，对一切虚位移，外力所作总虚功恒等于内力总虚功ATidxdyWieWWaSTSATbeldsFldxdyFW最小势能原理定义1：外力从位移状态退回到无位移的初始状态时所作的功称为外力势能定义2：形变势能和外力势能的和称为总势能最小势能原理：实际发生的位移总能使对应l的势能一阶变分为零。推导出总势能的二

3、阶变分为正，所以实际存在的位移使变形体的总势能取极小值。)(*aSTSATbpldsFldxdyFE)(21aSTSATbATpldsFldxdyFdxdyE0pE 有限元法的基本方程单元位移函数),( | ),( | ) 1 ,(1111vuyxi),( | ),( | )3 ,(3333vuyxk),( | ),( | )2 ,(2222vuyxjyaxaayxvyaxaayxu654321),(),(假设任意点的位移：单元位移函数：332211332211),(),(),(),(),(),(),(),(vyxNvyxNvyxNyxvuyxNuyxNuyxNyxu Nvu或：3213210

4、00000NNNNNNN Tvuvuvu332211插值函数(形函数)321321000000NNNNNNN)()(),(312312133221yxyxyxyxyxyxycxbayxNiiiijiiyxNiyxNijiiii,3 , 2 , 1, 0),(3 , 2 , 1, 1),(311),(iiyxN形函数特点：单元应变矩阵 BNLvuxyyxxyyx00 xNyNxNyNxNyNyNyNyNxNxNxNNLB332211321321000000单元应变矩阵(几何矩阵)： 321332211321321000000BBBxNyNxNyNxNyNyNyNyNxNxNxNB 3 , 2 ,

5、 1,002100ibccbxNyNyNxNBiiiiiiiii单元应力矩阵 SBDDxyyx 321321SSSBBBDBDS单元应力矩阵： iiBDS 单元刚度矩阵岩土体或结构体发生虚位移，单元结点的虚位移为 ，相应的虚应变为 ，则根据虚功原理有： * tdAtdsPNtdAFNTAATTATTnn* * eeTeTeATTeTkFtdABDBFn* tBDBtdABDBkTATen单元刚度矩阵：总体刚度矩阵由于虚位移 的任意性，等式两边与其相乘的矩阵相等，则：设结构体剖分成n个单元，根据虚功有：总体刚度矩阵 * eekF PUKkFnieeTnieT1*1* 荷载列阵结点位移列阵总体刚度

6、矩阵PUK等参元分析平面任意四边形单元 Tmmvuvuvu2211结点位移矩阵：01234mmNNNNNNN0000002121插值函数：)1)(1 (41),()1)(1 (41),(43NN)1)(1 (41),()1)(1 (41),(21NN mixNyNyNxNBiiiii, 3 , 2 , 100几何矩阵： mBBBB21 iiiiNNJyNxN1单元刚度矩阵： mmpllpiTiTTAeWWJtBDBddJtBDBtdABDBkn111111)( 模型范围与边界效应模型范围岩土工程涉及无限域或半无限域，但处理问题时只能对有限域进行离散化模型范围可取结构体轮廓尺寸的34倍边界效应以

7、小变形理论为基础的有限元法中，力与变形的影响范围是无限域，因此，设定有限域，并假定模型边界的位移为零或为受力边界就会带来误差，靠近边界越近误差越大，靠近边界越远误差越小 初始地应力场与释放荷载初始地应力场自然状态的岩体处于一定的初始地应力状态，在结构荷载作用下，岩体内的应力为荷载产生的应力与初始地应力之和释放荷载由于初始地应力的存在，开挖将导致部分岩体卸荷，通常采用沿开挖面作用着与地应力等价的“开挖释放荷载” 施工建造过程的模拟开挖释放荷载1ii1ihvihPivP空单元施工过程模拟Hh第1步开挖第2步开挖第3步开挖Hv空单元空单元空单元 节理及不连续面的模拟平面问题节理单元-GoodmanG

8、oodman单元是无厚度4结点单元结点传递切向力与法向力节理应力-应变关系nsnsnsuuKK00 xynsns节理单元本构模型nFnutannksFsutannsckmaxsF00节理单元刚度矩阵(假定位移沿单元长度线性变化) nnnnssssnnnsssnnssnseKKKKKKKKKKKKKKKKKKKKtlk200020200202202002020200202026称）（对考虑嵌入的节理单元模拟(考虑转动) nnnnssssnnnsssnnssnseKKKKKKKKKKKKKKKKKKKKtlk200020200202202002020200202024称）（对 nsnsnsuuKK

9、KM0000000nKLK341变厚度节理单元平面六结点变厚度节理单元相当于四边形等参元2134566161iiiiiivNuNvu位移函数：41)1)(1 (41iNiiii6 , 5)1)(1 (212iNii形函数:具有一定厚度的单元可按四边形等参元处理当厚度很小时按等厚度或无厚度节理单元处理 1111ddJBDBdVBDBkTVTensnsnsuuKK00 多节理岩体的模拟等效连续体利用节理单元来模拟密集分布或随机分布的节理与裂隙是不适宜的，会给离散化与计算带来诸多困难与麻烦从总体上考虑节理裂隙对岩体的影响，将岩体视为等效的正交异性、各向异性或各向同性体，即等效连续体将岩体中的节理裂隙

10、当成存在于岩体材料内的一种损伤，因如损伤力学的原理建立损伤模型来考虑节理裂隙的影响层状岩体均质各向同性岩体受彝族节理(层理)切割形成层状节理岩体，由于节理弱面的影响使岩体具有横观各向同性的特征，可按一般横观各向同性的连续体来建立有限元模型岩层走向与纵轴平行且该纵轴为一应力主轴时xyyxxyyxmGnvvvvvnvM0001)1 (0)1 (122121222)1 (2/)21)(1/(1111222111vEGEEnnvvvEMyxz岩层走向与纵轴正交，计算平面平行于层面xyyxxyyxnvvvnvvnvvnvM2210001012212222222122)1 (2/)21)(1/(11112

11、22111vEGEEnnvvvEMyxz 岩体工程中的弹塑性问题非线性分析的基本方法分段线性增量法：将总荷载分成若干增量 iiiPUK1荷载增量与位移增量的关系： niiUUU10总位移： niiPPP10总荷载：误差修正方法nPPU1P0P0U误差 11iiiiPPUK一阶自校正法：iiiUUU11UPU0P0U误差11iAiiPPUK牛顿迭代法：2UAPAU1P 有限元法的实现模型建立(范围及参数)前处理(模型剖分)形状函数几何矩阵本构关系单元刚度矩阵总体刚度矩阵单元结点位移结点位移列阵边界条件结点位移列阵模型应力与应变场2 边界元法 概述边界元法是同有限元法并行发展的另一类数值方法,该方

12、法在岩石力学中的应用自20世纪70年代以后有了较大的发展边界元法通常只须在边界上进行离散化,因而具有数据处理工作量小、占内存小、速度快等优点，但在处理多介质问题、复杂的非线性问题时效率低边界元法有两种直接法：直接建立关于边界未知量的积分方程，通过离散化求得边界未知量，并进而求域内任一点的场函数值间接法：设定一个在域内满足支配方程但包含若干未知系数的解，在边界上强迫其满足边界条件，求得该系数，进而求得边界上及域内各点的场函数值域边界元 直接边界元法基本方程*,ubtubt,相同结构第一状态下体积力、边界力与位移场：相同结构第二状态下体积力、边界力与位移场：由功的互等定理：dbudstuudbud

13、st*第一种情况外力：在无限域 上i点，沿l方向施加单位集中力 内力：在轮廓线 上，k方向的应力 位移：在 内及 上，任一点在k方向的位移第二种情况体力：在无限域 上沿k方向有分布体力表面荷载：在轮廓线 上，沿k方向荷载位移：在 内及 上，任一点在k方向的位移il*klP*klukbkPku由功的互等定理：dbudsPududsuPkklkklkilkkl*函数，具有如下性质：被称为DeltaDiracil 内点在外点在外点在ixfixfidxfiiil,21, 0上点在区域，上点在边界，式中：iiCuCduiilikil121)(2),( 1),(2),( 1(*方向方向方向方向yxlyxk

14、dsPudsuPuCdbudsPududsuPkklkklilikklkklkilkkl当不考虑体力时：将边界离散成n个线段单元并假设 与 沿边界均匀分布： dsuUdsPPIPUuPuICPdsuudsPuCTijTijnjjijnjjijiinjjkklnjjkklilijj)()(*111*1*；二阶单位矩阵；式中：kPku直接边界元法边界支配方程： PUuH 边界位移列阵u 边界应力列阵P nnnnnPPPPccH1111100边界应力影响系数矩阵：边界位移影响系数矩阵： nnnnUUUUU1111 间接边界元法基本方程(不连续应力法)外域),(yxpj内域dsyx上分布的荷载沿闭合边

15、界曲线),(yxp dsPdsPPFFFyxyx作用在微段ds上的荷载为：曲线上所有荷载在j点产生的位移： dsPuuj*曲线上所有荷载在j点产生的应力： dsPSj*将边界离散成n个线段单元并假设 在单元内均匀分布： niijiniLijijPUPdsudsPuui11* niijiniLijijPSPdsSdsPSi11*曲线上所有荷载在j点产生的位移：曲线上所有荷载在j点产生的应力：间接边界元法边界支配方程： 应力影响系数矩阵量边界元表面应力总体向量边界元表面位移总体向式中： 1 SuKSUu边界元法求解平面问题的步骤模型建立(范围及参数)将边界划分成单元将原岩应力反作用在单元上利用基本

16、方程求解边界单元上的作用力与位移利用开尔文基本解与功的互等定理求解内部点的应力与位移3 有限差分法有限差分法概述有限差分方法是将所有研究区域内的基本控制微分物理方程与边界条件近似差分方程表示，而将求解微分方程的问题变成在研究区域内特殊点上求解代数方程的问题，这些变量没有在单元内部定义。相比而言，有限元方法有一个重要的前提：应力和位移场变量应由参数控制的特征函数，以指定的模式在每一个单元内部变化。因此，有限元方法经常将单元矩阵合并为一个大的总刚度矩阵，然而，有限差分法却不这样作，而是有效地在每一步重新生成有限差分方程。 FLAC概述FLAC 是快速拉格郎日差分分析（Fast Lagrangian

17、 Analysis of Continua）的简写。这种算法可以准确地模拟材料的屈服、塑性流动、软化直至大变形，尤其在材料的弹塑性分析、大变形分析以及模拟施工过程等领域有其独到的优点。FLAC采用快速拉格朗日元法，基于显式差分来获得模型的全部运动方程的时间步长解。程序将计算模型划分为若干个不同形状的三维单元，单元之间用节点相互连接。对某一个节点施加荷载之后，该节点的运动方程可以写成时间步长的有限差分形式。对某一个微小的时间内，作用于该点的荷载只对周围的若干节点有影响。根据单元节点的速度变化和时间，程序可以求出单元之间的相对位移，进而可以求出单元应变；根据单元材料的本构方程可以求出单元应力。随着

18、时间的推移，这一过程将扩展到整个计算范围，直到边界。这样呈现可以追踪模型从渐进破坏直至整个破坏的全过程。 平衡方程 （运动方程） 应力、应变关系 （本构关系）新的应力或力新的速率和位移初始化平衡 平衡判断mechanic atioprecision） 终止该部分计算YesNo 图 FLAC基本显式计算循环 左图中表明了FLAC所包含的一般计算过程。这个过程首先调用运动方程从应力和外力导出了新的速度和位移，据速度导出应变速率，再由应变速率导出新的应力。对应于循环圈的每一个时步,值得注意的是，图中的每一个方框都根据已知值更新了网格变量，而这些已知值在方框内部操作时是保持恒定的。 差分公式有限差分离

19、散化基础：以增量之比倒替连续导数xuxudxdux0limyxxxi1i1iuxu有限差分网格连续函数 的泰勒展开yxu,3033202200! 31! 21xxuxxuxxuuu30332022003! 31! 21xxuxxuxxuuu30332022001! 31! 21xxuxxuxxuuu以0点为x，y坐标原点：203102231022xuuuxuxuuxu略去高阶微量：同理：204202242022yuuuyuyuuyuxuuxu2310差分公式(一阶二阶)：20310222xuuuxuyuuyu242020420222yuuuyu 75860241uuuuyxyxu)22(211

20、1931032uuuuxxu差分公式(三阶四阶)：11930140444641uuuuuxxu876543210220224222241uuuuuuuuuyxyxu)22(21121042032uuuuyyu121042040444461uuuuuyyu应力函数的差分解xyxyxyyx22222为应力函数，则：设 8675220200312022004220220412121xyxyxxyyxyyx利用差分公式： FLAC特点FLAC适用于多种材料模式与边界条件的非规则区域的连续问题求解； 在求解过程中，FLAC 采用了离散元的动态松驰法，不需要求解大型联立方程组（刚度矩阵）； FLAC 不但

21、可以对连续介质进行大变形分析，而且能模拟岩体沿某一软弱面产生的滑动变形； FLAC 还能在同一计算模型中针对不同的材料特性，使用相应的本构方程来比较真实地反映实际材料的动态行为，程序采用人机交互式的批命令形式执行； FLAC优点对模拟塑性破坏和塑性流动采用的是“混合离散法“。这种方法比有限元法中通常采用的“离散集成法“更为准确、合理； 即使模拟的系统是静态的，仍采用了动态运动方程，这使得FLAC3D在模拟物理上的不稳定过程不存在数值上的障碍；采用了一个“显式解”方案。因此，显式解方案对非线性的应力-应变关系的求解所花费的时间，几互与线性本构关系相同，而隐式求解方案将会花费较长的时间求解非线性问

22、题,没有必要存储刚度矩阵 。 FLAC缺点对于线性问题的求解，FLAC3D比有限元程序运行得要慢 ； 用FLAC3D求解时间取决于最长的自然周期和最短的自然周期之比。但某些问题对模型是无效的。 FLAC3D分析步骤有限差分网格 本构特性与材料性质 边界条件与初始条件求解结果分析(标题菜单中editcopy to clipboard) 完成上述工作后，可以获得模型的初始平衡状态，也就是模拟开挖前的原岩应力状态。然后，进行工程开挖或改变边界条件来进行工程的响应分析，进行一系列计算步后达到问题的解。最后对结果进行分析总结。 FLAC3D建模方法直接法。直接法是按照分析对象的几何形状利用FLAC3D内

23、置的网格生成器建模，网格和几何模型同时生成。 间接法。通过点、线、面、体，先建立对象的几何外形，再进行实体模型的分网（Meshing），以完成网格模型的建立，FLAC3D自身不具备间接法建模功能，读者可借助第三方软件与FLAC3D的接入轻松实现。 FLAC3D使用Generate zone生成基本形状网格的常用关键词关键词含义关键词含义dimension内部区域的尺寸p0p16各种形状网格的参考点edge网格的边长ratio单元尺寸大小比率fill用网格填充内部区域size网格在每个坐标方向上的单元数目dip对称面与水平面的夹角dd对称面的法线在xy面的投影与y轴的正向夹角 FLAC3D几何模

24、型单元类型块体单元网格(Brick) 例如一长为8m，宽为6m， 高为8m的长方体，建模：gen zone brick p0 (0,0,0)& p1(6,0,0) p2(0,8,0) p3(0,0,8)& size 6 8 8退化的块体单元网格(Dbrick)gen zone dbrick p0(0,0,0)& p1(10,0,0)p2(0,10,0)& p3(0,0,8) size 5 5 8楔形单元网格(Wedge)gen zone wedge p0(0,0,0)& p1(5,0,0) p2(0,6,0)& p3(0,0,3) size 10

25、 12 6棱锥体单元网格(Pyramid)gen zone pyramid p0(0,0,0)& p1(10,0,0) p2(0,8,0) & p3(0,0,6) size 10 8 6四面体单元网格(Tetrahedron)gen zone tet p0(0,0,0)& p1(10,0,0) p2(0,8,0) & p3(0,0,6) size 10 8 6圆柱体单元网格(Cylinder)gen zone cyl p0 0 0 0& p1 1 0 0 p2 0 4 0& p3 0 0 1 size 4 4 6 放射状块体单元网格(Radbri

26、ck)gen zone radbrick p0 (0,0,0)& p1 (10,0,0) p2 (0,10,0)& p3 (0,0,10) size 3,5,5,7 & ratio 1,1,1,1.5 dim 1 4 2放射状矩形隧洞网格(Radtunnel)gen zone radtun p0 (0,0,0)& p1 (10,0,0) p2 (0,10,0)& p3 (0,0,10) size 3,5,10,7 & ratio 1,1,1,1.5 dim 2 1 2 1 fill放射状圆柱形隧洞网格 (Radcylinder) gen zone

27、 radcyl p0(0,0,0)& p1(10,0,0) p2(0,20,0)& p3(0,0,12)size 5 10 6 12& dim 2 1 2 1圆柱形壳体单元网格(Cshell)gen zon cshell p0 0 0 0& p1 6.0 0 0 p2 0 10 0& p3 0 0 5.0 size 3 10 8 3& dim 5.6 4.6 5.6 4.6 圆柱形交叉形隧洞单元网格(Cylint)gen zon cylint p0 0 0 0& p1 7.5 0 0 p2 0 7.5 0& p3 0 0 4.75

28、 dim 2 2 2 2 2 2 2& size 6 6 6 6 6 ratio 1 1 1 1 矩形交叉形隧洞单元网格(Tunint)gen zon tunint p0 0 0 0& p1 7.5 0 0 p2 0 7.5 0& p3 0 0 4.75 dim 2 2 2 2 2 2 2& size 6 6 6 6 6 ratio 1 1 1 1 网格单元间的连接 采用FLAC3D进行计算，所建立的模型需是一个连续的整体，否则计算结果将出现较大的误差甚至无法进行计算。 对于在建立模型时，各关键点的坐标是准确无误输入且各公共面的网格数和大小均完全一致的模型，无需

29、进行任何操作，模型即自动完成相互间的连接。此处所讨论的是公共面上网格的大小和划分的份数不一致(存在整数倍关系)或者两相邻面间存在间隙的问题。对于公共面上网格的大小和划分的份数不一致，主要采用attach命令来进行连接，而对两相邻面间存在间隙的问题，则采用gen merge 命令来进行连接的操作。 Gen zone brick size 4 4 4 Gen zone brick size 4 4 4 p0 4.1 0 0 p1 8 0 0 p2 4.1 4 0 p3 4.1 0 4 Gen merge 0.1尺寸和份数要求为？gen zone brick size 4 4 4 p0 0,0,0

30、p1 4,0,0 p2 0,4,0 p3 0,0,2gen zone brick size 8 8 4 p0 0,0,2 p1 4,0,2 p2 0,4,2 p3 0,0,4attach face range z 1.9 2.1model elasprop bulk 8e9 shear 5e9fix z range z -.1 .1fix x range x -.1 .1fix x range x 3.9 4.1fix y range y -.1 .1fix y range y 3.9 4.1apply szz -1e6 range z 3.9 4.1 x 0,2 y 0,2solvesave

31、 att.sav一建模例子gen zone radcyl p0 0,0,0 p1 100,0,0 p2 0,200,0 p3 0,0,100 size 5 10 6 12 dim 4 4 4 4 ratio 1 1 1 1.2 gen zone radtun p0 0,0,0 p1 0,0,-100 p2 0,200,0 p3 100,0,0 size 5 10 5 12 dim 4 4 4 4 ratio 1 1 1 1.2 gen zone reflect dip 90 dd 90注意：*为了避免边界效应的影响，模型的外边界尺寸最小应为应为地下开挖洞室直径的35倍。*为了保证计算结果的准确

32、性，洞室周围的网格划分应密集一些，但也不宜过细，否则会影响计算速度。同时，单元体的边长比值要控制在一定的范围内，尽量避免比较狭长的单元体出现。*为了保证网格的连续性，应保证相邻边界节点的匹配，例如，相邻网格有相同的单元体数和一致的单元体几何变化率。 FLAC3D边界条件和初始条件边界条件。在边界区域可以指定速度（位移）边界条件或应力（力）边界条件。 APPLY 在模型边界上施加力，流体流动和温度条件。 如：Apply szz=-1e5 sxz=-0.5e5 range z=(-.1,.1) 。APPLY nstress 表示施加一个法向应力。可以使用xforce, yforce, zforce

33、关键字将力施加到给定的节点。APPLY sxx-1e6 gradient 0,0,1e5 range z -100,0，可以施加渐变力。 FIX 是保持网格节点指定参数（速度、压力和温度）的值不变。如fix x range x -0.1 0.1。也可以给出初始应力条件，包括重力荷载以及地下水位线。所有的条件都充许指定变化梯度。 Initial sxx=-5e6 syy=-1e7 szz=-5e6，表示施加一个初始应力，并贯穿到整个网格。 FLAC3DFLAC3D本构模型本构模型开挖模型开挖模型null 3个弹性模型：各向同性弹性；横观各向同性弹性；个弹性模型：各向同性弹性；横观各向同性弹性；正

34、交各向同性弹性。正交各向同性弹性。7个塑性模型：个塑性模型：Drucker-Prager Mode德鲁克德鲁克普普拉格模型；拉格模型；Morh-Coulomb Mode摩尔摩尔库仑模型；库仑模型；StrainHardening/Softening Model应变硬化应变硬化/软化模型；软化模型；UbiquitousJoint Model遍布节理遍布节理模型；模型；Bilinear StrainHardening/Softening Model双线性应变硬化双线性应变硬化/软化模型；软化模型；Modified CamClay Model 修正剑桥模型；修正剑桥模型；DoubleYield Mod

35、el 胡克布朗模型。胡克布朗模型。模型模型材料特性材料特性实际应用实际应用空模形空模形 Null空空孔洞，挖开，后续施工材料（如回填）孔洞，挖开，后续施工材料（如回填）各向同性弹性模型各向同性弹性模型 elastic均匀各向同性的线性本构关系均匀各向同性的线性本构关系低于强度极限的人工材料（如钢筋）；安全系低于强度极限的人工材料（如钢筋）；安全系数计算数计算正交各向同性弹性正交各向同性弹性 orthotropic正交各向同性材料正交各向同性材料不超过强度极限的柱状岩体不超过强度极限的柱状岩体横观各向同性弹性横观各向同性弹性anisotropic横观各向同性弹性（即板岩）横观各向同性弹性（即板岩

36、）不超过强度极限的层压材料不超过强度极限的层压材料德鲁克德鲁克普拉格模型普拉格模型drucker极限分析，低摩擦角的软粘土极限分析，低摩擦角的软粘土与隐式有限元程序相比的常用模型与隐式有限元程序相比的常用模型摩尔摩尔库仑模型库仑模型 mohr松散或胶结的粒状材料：土，松散或胶结的粒状材料：土，岩石，混凝土岩石，混凝土岩土力学通常用的模型（边坡稳定性分析，地岩土力学通常用的模型（边坡稳定性分析，地下开挖）下开挖）应变硬化应变硬化/软化摩尔软化摩尔库仑模库仑模型型 ss存在非线性硬化或软化的粒状存在非线性硬化或软化的粒状材材料料破坏后研究（失稳过程，立柱屈服，顶板崩落）破坏后研究（失稳过程，立柱屈

37、服，顶板崩落）遍布解理模型遍布解理模型 ubiquitous具有强度各向异性的层状材料具有强度各向异性的层状材料（即板岩）（即板岩）松散沉积地层中的开挖松散沉积地层中的开挖双线性应变硬化双线性应变硬化/软化遍布解软化遍布解理模型理模型subiquitous具有非线性材料硬化或软化的具有非线性材料硬化或软化的层层状材料状材料层状材料破坏后研究层状材料破坏后研究双屈服面塑性模型双屈服面塑性模型doubleyield轻胶结的粒状材料，在压力作轻胶结的粒状材料，在压力作用用下导致永久体积减小下导致永久体积减小修正剑桥模型修正剑桥模型 canrclay变形和抗剪强度是体变的函数变形和抗剪强度是体变的函数

38、粘土粘土胡克胡克布朗模型布朗模型各向同性的岩石材料各向同性的岩石材料岩石岩石各向同性弹性模型：各向同性弹性模型：剪切模量；体积模量；密度,摩尔摩尔-库伦模型：库伦模型：剪切模量；体积模量；密度；摩擦角；粘结力；剪胀角 (任选)。横观各向同性模型：横观各向同性模型：剪切模量；x-模量；y-模量；密度； NUYx；NUYz 。砌体节理模型砌体节理模型：剪切模量；体积模量；密度；粘结力 (整体材料)；摩擦力 (整体材料)；节理粘结力；节理摩擦力；节理角；整体材料的剪胀角 (任选) 。德鲁克德鲁克-普拉格模型普拉格模型：剪切模量；体积模量；密度；材料参数（Kshear；材料参数（Qdil）；材料参数（

39、Qvol）；抗拉强度。 本构模型的选择本构模型的选择 Model range 如：如：model mohr；model null range x=2,4 y=2,6 z=5,10 Prop prop bulk=1e8 shear=0.3e8 fric=35 prop coh=1e10 tens=1e10 注意注意:（1）材料的本构模型必须先定义，以便绘图）材料的本构模型必须先定义，以便绘图或显示材料参数。（或显示材料参数。（2）如果材料参数的关键字与本构）如果材料参数的关键字与本构模型不协调，则弹出警告信息。（模型不协调，则弹出警告信息。（3）本构模型需要的）本构模型需要的参数没有指定时，系统

40、使用默认值，一般为参数没有指定时，系统使用默认值，一般为0。 FLAC3DFLAC3D本构模型一例子本构模型一例子 gen zone cyl p0 0 0 0 p1 1 0 0 p2 0 2 0 p3 0 0 1 size 4 5 4 gen zone reflect norm 1,0,0 gen zone reflect norm 0,0,1 model mohr prop bulk 1.19e10 shear 1.1e10 prop coh 2.72e5 fric 44 ten 2e5 fix x y z range y -.1 .1 fix x y z range y 1.9 2.1 i

41、ni yvel 1e-7 range y -.1 .1 ini yvel -1e-7 range y 1.9 2.1 hist gp ydisp 0,0,0 hist zone syy 0,1,0 hist zone syy 1,1,0 step 3000蠕变模型蠕变模型 经典的粘弹性模型，是经典的粘弹性模型，是Maxwell体的经典表达式；体的经典表达式； 二分量幂定律，可用于采矿业；二分量幂定律，可用于采矿业；用于核废料隔离研究的参考蠕变公式（用于核废料隔离研究的参考蠕变公式（WIPP模型），一般用模型），一般用于研究盐矿中核废料的地下储藏的有关热力学分析，也适合于研究盐矿中核废料的地下储

42、藏的有关热力学分析，也适合于软土的变形特性分析于软土的变形特性分析 ； 伯格（伯格（Burger）蠕变模型和摩尔库伦模型合成的伯格蠕变）蠕变模型和摩尔库伦模型合成的伯格蠕变粘塑性模型，它是在经典的粘弹性模型的基础上扩展而成的，粘塑性模型，它是在经典的粘弹性模型的基础上扩展而成的，它包括一个开尔文体和一个摩尔库伦体；它包括一个开尔文体和一个摩尔库伦体； WIPP模型和德鲁克普拉格（模型和德鲁克普拉格（Drucker-Prager）模型合成）模型合成的的WIPP蠕变粘塑性模型，它是蠕变粘塑性模型，它是WIPP模型的变化形式，包括模型的变化形式，包括一个德鲁克普拉格（一个德鲁克普拉格（Drucker

43、-Prager）塑性体；）塑性体； 岩盐的本构模型，它也是岩盐的本构模型，它也是WIPP模型的变化形式，它包括体积模型的变化形式，它包括体积和偏量压实特性。和偏量压实特性。 FLAC3DFLAC3D结构单元结构单元FLAC3D 可以模拟可以模拟 6种结构单元（梁、锚索、桩及种结构单元（梁、锚索、桩及板壳、格栅、衬砌单元）板壳、格栅、衬砌单元） FLAC3D 可以模拟任意形状、任意特性的结构体与可以模拟任意形状、任意特性的结构体与岩土体的相互作用，以及力作用在结构体上或岩土岩土体的相互作用，以及力作用在结构体上或岩土体时，结构体和岩土体的力学反应。体时，结构体和岩土体的力学反应。 结构单元力学性

44、态的计算模拟可以通过几何大变形结构单元力学性态的计算模拟可以通过几何大变形和小变形两种计算模式计算，而且还设置动力状态和小变形两种计算模式计算，而且还设置动力状态来模拟结构与岩土体的动力反应。来模拟结构与岩土体的动力反应。梁单元梁单元需输入的参数需输入的参数density：density（梁单元的密度）（梁单元的密度）, Emod：Youngs modulus（弹性模量）（弹性模量）, E nu：Poissons ratio（泊松比）（泊松比）, pmoment： plastic moment capacity（塑性矩）（塑性矩）, MP thexp：thermal expansion coe

45、fficient（热膨胀系数）（热膨胀系数）xcarea：cross-sectional area（横截面积）（横截面积）, A xciy：second moment with respect to beamSEL y-axis（截面对（截面对Y轴的惯性矩）轴的惯性矩）, Iy xciz：second moment with respect to beamSEL z-axis（截面对（截面对Z轴的惯性矩）轴的惯性矩）, Iz xcj：polar moment of inertia（极惯性矩）（极惯性矩）, J t梁单元的命令输入形式梁单元的命令输入形式SEL beam （定义梁单元）（定义梁单

46、元）SEL beamsel （定义梁单元条件）（定义梁单元条件） 可能的梁单元与已有的梁单元的连接是通过可能的梁单元与已有的梁单元的连接是通过id命令以及命令以及id的值来的值来确定确定 梁单元的命令一例子梁单元的命令一例子 sel beam id=1 begin=(0,0,0) end=(3,0,0) nseg=3 sel beam id=1 begin=(3,0,0) end=(6,0,0) nseg=4 sel beam id=1 begin=(6,0,0) end=(9,0,0) nseg=3 sel beam id=1 prop emod=2e11 nu=0.30 & xca

47、rea=6e-3 xcj=0.0 xciy=200e-6 xciz=200e-6selbeamselnodesnd1nd2selbeambeginx,y,zendx,y,zselbeamkeywordapplykeywordydistvaluezdistvaluepropertykeywordvalue(1)densityvalue(2)emodvalue(3)nuvalue(4)pmomentvalue(5)thexpvalue(6)xcareavalue(7)xciyvalue(8)xcizvalue(9)xcjvalue锚杆单元锚杆单元需输入的参数需输入的参数density：densit

48、y（锚杆密度）（锚杆密度）, Emod：Youngs modulus（弹性模量）（弹性模量）, E gr_coh：grout cohesive strength（水泥浆的粘聚力）（水泥浆的粘聚力）, cg gr_fric ：grout friction angle（水泥浆的内摩擦角）（水泥浆的内摩擦角）, g gr_k：grout stiffness（水泥浆的刚度）（水泥浆的刚度）, kg gr_per：grout exposed perimeter（锚杆的外周长）（锚杆的外周长）, pg xcarea：cross-sectional area（锚杆的横截面积）（锚杆的横截面积）, A yco

49、mpression：compressive yield strength（锚杆的抗压屈服（锚杆的抗压屈服强度）强度）, Fc ytension：tensile yield strength（锚杆的拉伸屈服强度）（锚杆的拉伸屈服强度）, Ft 锚杆单元的命令输入形式锚杆单元的命令输入形式SEL cable （定义锚杆单元）（定义锚杆单元）SEL cablesel （定义锚杆单元条件）（定义锚杆单元条件） 可能的锚杆单元与已有的锚杆单元的连接是通过可能的锚杆单元与已有的锚杆单元的连接是通过id命令以及命令以及id的值的值来确定来确定 锚杆单元的命令一例子锚杆单元的命令一例子 sel cable i

50、d=1 begin=(0.1, 0.5, 0.1) end=(11.9, 0.5, 0.1) & nseg=13 sel cable prop xcarea=2e-3 emod=200e9 yTens=1e20 & gr_k=1e10 gr_coh=1e20 sel cable prop gr_per=0.314 gr_fric=25selkeywordcableselnodel_idnode2_idcablebeginx1,y1,z1endx2,y2,z2pretensionvaluepropertykeyword(1)densityvalue(2)emodevalue(3)

51、xcareavalue(4)gr_cohvalue(5)gr_fricvalue(6)gr_kvalue(7)gr_pervalue(8)ycompressionvalue(9)ytensionvalue桩单元桩单元需输入的参数需输入的参数density：density（桩单元的密度）（桩单元的密度）,Emod：Youngs modulus（弹性模量）（弹性模量）, Enu：Poissons ratio（泊松比）（泊松比）, pmoment： plastic moment capacity（塑性矩）（塑性矩）, MP (可选的可选的除非特除非特别说明，假设力矩承载力是足够大的别说明，假设力矩承

52、载力是足够大的) thexp：thermal expansion coefficient（热膨胀系数）（热膨胀系数） xcarea：cross-sectional area（桩的横截面积）（桩的横截面积）, A perimeter：exposed perimeter（桩的外周长）（桩的外周长）, pxciy：second moment with respect to beamSEL y-axis（截面对（截面对Y轴轴的惯性矩）的惯性矩）, Iy xciz：second moment with respect to beamSEL z-axis（截面对（截面对Z轴的轴的惯性矩）惯性矩）, Izx

53、cj：polar moment of inertia（极惯性矩）（极惯性矩）, Jcs_scoh：shear coupling spring cohesion per unit length（切向耦合弹簧粘聚强度）（切向耦合弹簧粘聚强度）cs_sfric：shear coupling spring friction angle（切向耦合弹簧（切向耦合弹簧内摩擦角）内摩擦角）, s cs_sk：shear coupling spring stiffness per unit length（切向（切向耦合弹簧刚度）耦合弹簧刚度）, kscs_ncoh：normal coupling spring

54、cohesion per unit length（法向耦合弹簧粘聚强度）（法向耦合弹簧粘聚强度）, cn cs_nfric：normal coupling spring friction angle（法向耦合弹（法向耦合弹簧内摩擦角）簧内摩擦角）, n cs_ngap：normal coupling spring gap-use flag（桩和土交界（桩和土交界面之间的法向间隔）面之间的法向间隔）, g cs_nk：normal coupling spring stiffness per unit length（法（法向耦合弹簧刚度）向耦合弹簧刚度）, kn selpileselnodesnd

55、1nd2selpilebeginx,y,zendx,y,zselpilekeywordapplykeywordydistzdistpropertykeyword(1)densityvalue(2)emodvalue(3)nuvalue(4)pmomentvalue(5)thexpvalue(6)xcareavalue(7)perimetervalue(8)xciyvalue(9)xcizvalue(10)xcjvalue(11)cs_scohvalue(12)cs_sfricvalue(13)cs_skvalue(14)cs_ncohvalue(15)cs_nfricvalue(16)cs_n

56、kvalue(17)cs_ngapon/off壳单元壳单元需输入的参数需输入的参数density：density（壳单元的密度）（壳单元的密度）,isotropic：isotropic material properties: E and 和和 orthotropic：orthotropic material properties(正交各向正交各向异性材料的参数异性材料的参数)thexp：thermal-expansion coefficient(热膨胀系数热膨胀系数) thickness：thickness(厚度厚度), t 土工格栅单元土工格栅单元需输入的参数需输入的参数density：d

57、ensity（单元的密度）（单元的密度）,isotropic：isotropic material properties: E and 和和 orthotropic：orthotropic material properties(正交各向异性材料的参数正交各向异性材料的参数)thexp：thermal-expansion coefficient(热膨胀系数热膨胀系数) thickness：thickness(厚度厚度), t cs_scoh：coupling spring cohesion (stress units)(连接弹簧的粘聚力连接弹簧的粘聚力), c cs_sfric：couplin

58、g spring friction angle(连接弹簧的线的摩擦角连接弹簧的线的摩擦角) cs_sk：coupling spring stiffness per unit area(连接弹簧的每单位面积连接弹簧的每单位面积刚度刚度), k slide：large-strain sliding flag (default: off)(大应变滑移标记大应变滑移标记,默认为默认为关关) slide tol：large-strain sliding tolerance(容许的大应变滑移值容许的大应变滑移值) selgeogridselnodesnd1nd2nd3selgeogridbeginx,y,

59、zendx,y,zselgeogridkeywordapplypressurepropertykeyword(1)densityvalue(2a)isotropicvalue(2b)orthotropicvalue(3)thexpvalue(4)thicknessvaluecs_scohvaluecs_sfricvaluecs_skvalueslideon/offslide_tolvalue衬砌单元衬砌单元需输入的参数需输入的参数density：density（桩单元的密度）（桩单元的密度）,isotropic：isotropic material properties: E and 和和 o

60、rthotropic：orthotropic material properties(正交各向异性材料的参数正交各向异性材料的参数)thexp：thermal-expansion coefficient(热膨胀系数热膨胀系数) thickness：thickness(厚度厚度), t cs_ncut：normal coupling spring tensile strength (法向连接弹簧的拉伸强度法向连接弹簧的拉伸强度) cs_nk：normal coupling spring stiffness per unit area(法向连接弹簧单位面积上刚度法向连接弹簧单位面积上刚度)cs_scoh：shear coupling spring cohesion (stress units)(切向连接弹簧的粘聚力切向连接弹簧的粘聚力:应力应力单位单位), c cs_scohres：shear coupling spring residual