岩土工程数值分析方法.ppt

第7章 岩土工程数值分析方法,有限元法 边界元法 有限差分法 离散单元法,1 有限元法,概述基本思路：将复杂的结构看成由有限个仅在结点处联结的整体，首先对每一个单元分析其特性，建立相关物理量之间的相互联系。然后，依据单元之间的联系再将各单元组装成整体，从而获得整体特性方程，应用方程相应的解法，即可完成整个问题的分析分析过程：,结构离散化,确定单元位移模式,单元特性分析,集成总体特性,接方程求未知量,概述工程界比较流行、被广泛使用的大型有限元软件：MSC(航空航天领域，是目前规模最大的有限元分析系统)；Marc(非线性分析软件)；Adina(可进行结构、流体、热的耦合计算，具有隐式和显式两种时间积分算法，非线性运算功能强大)；ANSYS（可进行结构、流体、热、电磁场的计算，是有限元分析的通用软件包）优点：可分析几何形状及受荷条件复杂、非均质的各种实际结构；可在计算中模拟各种复杂的材料本构关系、边界条件等；前后处理技术先进。,有限单元法的理论基础虚位移原理：受给定外力的变形体处于平衡状态的充要条件是，对一切虚位移，外力所作总虚功恒等于内力总虚功,最小势能原理定义1：外力从位移状态退回到无位移的初始状态时所作的功称为外力势能定义2：形变势能和外力势能的和称为总势能最小势能原理：实际发生的位移总能使对应l的势能一阶变分为零。推导出总势能的二阶变分为正，所以实际存在的位移使变形体的总势能取极小值。,有限元法的基本方程单元位移函数,假设任意点的位移：,单元位移函数：,或：,插值函数(形函数),形函数特点：,单元应变矩阵,单元应变矩阵(几何矩阵)：,单元应力矩阵,单元应力矩阵：,单元刚度矩阵岩土体或结构体发生虚位移，单元结点的虚位移为，相应的虚应变为，则根据虚功原理有：,单元刚度矩阵：,总体刚度矩阵由于虚位移 的任意性，等式两边与其相乘的矩阵相等，则：设结构体剖分成n个单元，根据虚功有：总体刚度矩阵,等参元分析平面任意四边形单元,结点位移矩阵：,插值函数：,几何矩阵：,单元刚度矩阵：,模型范围与边界效应模型范围岩土工程涉及无限域或半无限域，但处理问题时只能对有限域进行离散化模型范围可取结构体轮廓尺寸的34倍边界效应以小变形理论为基础的有限元法中，力与变形的影响范围是无限域，因此，设定有限域，并假定模型边界的位移为零或为受力边界就会带来误差，靠近边界越近误差越大，靠近边界越远误差越小,初始地应力场与释放荷载初始地应力场自然状态的岩体处于一定的初始地应力状态，在结构荷载作用下，岩体内的应力为荷载产生的应力与初始地应力之和释放荷载由于初始地应力的存在，开挖将导致部分岩体卸荷，通常采用沿开挖面作用着与地应力等价的“开挖释放荷载”,施工建造过程的模拟开挖释放荷载,施工过程模拟,第1步开挖,第2步开挖,第3步开挖,空单元,空单元,空单元,节理及不连续面的模拟平面问题节理单元-GoodmanGoodman单元是无厚度4结点单元结点传递切向力与法向力节理应力-应变关系,节理单元本构模型,节理单元刚度矩阵(假定位移沿单元长度线性变化),考虑嵌入的节理单元模拟(考虑转动),变厚度节理单元平面六结点变厚度节理单元相当于四边形等参元,位移函数：,形函数:,具有一定厚度的单元可按四边形等参元处理当厚度很小时按等厚度或无厚度节理单元处理,多节理岩体的模拟等效连续体利用节理单元来模拟密集分布或随机分布的节理与裂隙是不适宜的，会给离散化与计算带来诸多困难与麻烦从总体上考虑节理裂隙对岩体的影响，将岩体视为等效的正交异性、各向异性或各向同性体，即等效连续体将岩体中的节理裂隙当成存在于岩体材料内的一种损伤，因如损伤力学的原理建立损伤模型来考虑节理裂隙的影响,层状岩体均质各向同性岩体受彝族节理(层理)切割形成层状节理岩体，由于节理弱面的影响使岩体具有横观各向同性的特征，可按一般横观各向同性的连续体来建立有限元模型岩层走向与纵轴平行且该纵轴为一应力主轴时,岩层走向与纵轴正交，计算平面平行于层面,岩体工程中的弹塑性问题非线性分析的基本方法分段线性增量法：将总荷载分成若干增量,荷载增量与位移增量的关系：,总位移：,总荷载：,误差修正方法,一阶自校正法：,牛顿迭代法：,有限元法的实现,模型建立(范围及参数),前处理(模型剖分),形状函数,几何矩阵,本构关系,单元刚度矩阵,总体刚度矩阵,单元结点位移,结点位移列阵,边界条件,结点位移列阵,模型应力与应变场,2 边界元法,概述边界元法是同有限元法并行发展的另一类数值方法,该方法在岩石力学中的应用自20世纪70年代以后有了较大的发展边界元法通常只须在边界上进行离散化,因而具有数据处理工作量小、占内存小、速度快等优点，但在处理多介质问题、复杂的非线性问题时效率低,边界元法有两种直接法：直接建立关于边界未知量的积分方程，通过离散化求得边界未知量，并进而求域内任一点的场函数值间接法：设定一个在域内满足支配方程但包含若干未知系数的解，在边界上强迫其满足边界条件，求得该系数，进而求得边界上及域内各点的场函数值,直接边界元法基本方程,相同结构第一状态下体积力、边界力与位移场：,相同结构第二状态下体积力、边界力与位移场：,由功的互等定理：,第一种情况外力：在无限域 上i点，沿l方向施加单位集中力 内力：在轮廓线 上，k方向的应力 位移：在 内及 上，任一点在k方向的位移第二种情况体力：在无限域 上沿k方向有分布体力表面荷载：在轮廓线 上，沿k方向荷载位移：在 内及 上，任一点在k方向的位移,由功的互等定理：,当不考虑体力时：,将边界离散成n个线段单元并假设 与 沿边界均匀分布：,直接边界元法边界支配方程：,边界应力影响系数矩阵：,边界位移影响系数矩阵：,间接边界元法基本方程(不连续应力法),外域,j,内域,作用在微段ds上的荷载为：,曲线上所有荷载在j点产生的位移：,曲线上所有荷载在j点产生的应力：,将边界离散成n个线段单元并假设 在单元内均匀分布：,曲线上所有荷载在j点产生的位移：,曲线上所有荷载在j点产生的应力：,间接边界元法边界支配方程：,边界元法求解平面问题的步骤,模型建立(范围及参数),将边界划分成单元,将原岩应力反作用在单元上,利用基本方程求解边界单元上的作用力与位移,利用开尔文基本解与功的互等定理求解内部点的应力与位移,3 有限差分法,有限差分法概述有限差分方法是将所有研究区域内的基本控制微分物理方程与边界条件近似差分方程表示，而将求解微分方程的问题变成在研究区域内特殊点上求解代数方程的问题，这些变量没有在单元内部定义。相比而言，有限元方法有一个重要的前提：应力和位移场变量应由参数控制的特征函数，以指定的模式在每一个单元内部变化。因此，有限元方法经常将单元矩阵合并为一个大的总刚度矩阵，然而，有限差分法却不这样作，而是有效地在每一步重新生成有限差分方程。,FLAC概述FLAC 是快速拉格郎日差分分析（Fast Lagrangian Analysis of Continua）的简写。这种算法可以准确地模拟材料的屈服、塑性流动、软化直至大变形，尤其在材料的弹塑性分析、大变形分析以及模拟施工过程等领域有其独到的优点。FLAC采用快速拉格朗日元法，基于显式差分来获得模型的全部运动方程的时间步长解。程序将计算模型划分为若干个不同形状的三维单元，单元之间用节点相互连接。,对某一个节点施加荷载之后，该节点的运动方程可以写成时间步长的有限差分形式。对某一个微小的时间内，作用于该点的荷载只对周围的若干节点有影响。根据单元节点的速度变化和时间，程序可以求出单元之间的相对位移，进而可以求出单元应变；根据单元材料的本构方程可以求出单元应力。随着时间的推移，这一过程将扩展到整个计算范围，直到边界。这样呈现可以追踪模型从渐进破坏直至整个破坏的全过程。,左图中表明了FLAC所包含的一般计算过程。这个过程首先调用运动方程从应力和外力导出了新的速度和位移，据速度导出应变速率，再由应变速率导出新的应力。对应于循环圈的每一个时步,值得注意的是，图中的每一个方框都根据已知值更新了网格变量，而这些已知值在方框内部操作时是保持恒定的。,差分公式有限差分离散化基础：以增量之比倒替连续导数,有限差分网格,连续函数 的泰勒展开,以0点为x，y坐标原点：,略去高阶微量：,同理：,差分公式(一阶二阶)：,差分公式(三阶四阶)：,应力函数的差分解,利用差分公式：,FLAC特点FLAC适用于多种材料模式与边界条件的非规则区域的连续问题求解；在求解过程中，FLAC 采用了离散元的动态松驰法，不需要求解大型联立方程组（刚度矩阵）；FLAC 不但可以对连续介质进行大变形分析，而且能模拟岩体沿某一软弱面产生的滑动变形；FLAC 还能在同一计算模型中针对不同的材料特性，使用相应的本构方程来比较真实地反映实际材料的动态行为，程序采用人机交互式的批命令形式执行；,FLAC优点对模拟塑性破坏和塑性流动采用的是“混合离散法“。这种方法比有限元法中通常采用的“离散集成法“更为准确、合理；即使模拟的系统是静态的，仍采用了动态运动方程，这使得FLAC3D在模拟物理上的不稳定过程不存在数值上的障碍；采用了一个“显式解”方案。因此，显式解方案对非线性的应力-应变关系的求解所花费的时间，几互与线性本构关系相同，而隐式求解方案将会花费较长的时间求解非线性问题,没有必要存储刚度矩阵。,FLAC缺点对于线性问题的求解，FLAC3D比有限元程序运行得要慢；用FLAC3D求解时间取决于最长的自然周期和最短的自然周期之比。但某些问题对模型是无效的。,FLAC3D分析步骤有限差分网格 本构特性与材料性质 边界条件与初始条件求解结果分析(标题菜单中editcopy to clipboard)完成上述工作后，可以获得模型的初始平衡状态，也就是模拟开挖前的原岩应力状态。然后，进行工程开挖或改变边界条件来进行工程的响应分析，进行一系列计算步后达到问题的解。最后对结果进行分析总结。,FLAC3D建模方法直接法。直接法是按照分析对象的几何形状利用FLAC3D内置的网格生成器建模，网格和几何模型同时生成。间接法。通过点、线、面、体，先建立对象的几何外形，再进行实体模型的分网（Meshing），以完成网格模型的建立，FLAC3D自身不具备间接法建模功能，读者可借助第三方软件与FLAC3D的接入轻松实现。,FLAC3D使用Generate zone生成基本形状网格的常用关键词,FLAC3D几何模型单元类型块体单元网格(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 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,放射状块体单元网格(Radbrick)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 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 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进行计算，所建立的模型需是一个连续的整体，否则计算结果将出现较大的误差甚至无法进行计算。对于在建立模型时，各关键点的坐标是准确无误输入且各公共面的网格数和大小均完全一致的模型，无需进行任何操作，模型即自动完成相互间的连接。此处所讨论的是公共面上网格的大小和划分的份数不一致(存在整数倍关系)或者两相邻面间存在间隙的问题。对于公共面上网格的大小和划分的份数不一致，主要采用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 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 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倍。为了保证计算结果的准确性，洞室周围的网格划分应密集一些，但也不宜过细，否则会影响计算速度。同时，单元体的边长比值要控制在一定的范围内，尽量避免比较狭长的单元体出现。为了保证网格的连续性，应保证相邻边界节点的匹配，例如，相邻网格有相同的单元体数和一致的单元体几何变化率。,FLAC3D边界条件和初始条件边界条件。在边界区域可以指定速度（位移）边界条件或应力（力）边界条件。APPLY 在模型边界上施加力，流体流动和温度条件。如：Apply szz=-1e5 sxz=-0.5e5 range z=(-.1,.1)。APPLY nstress 表示施加一个法向应力。可以使用xforce,yforce,zforce关键字将力施加到给定的节点。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，表示施加一个初始应力，并贯穿到整个网格。,FLAC3D本构模型开挖模型null 3个弹性模型：各向同性弹性；横观各向同性弹性；正交各向同性弹性。7个塑性模型：Drucker-Prager Mode德鲁克普拉格模型；Morh-Coulomb Mode摩尔库仑模型；StrainHardening/Softening Model应变硬化/软化模型；UbiquitousJoint Model遍布节理模型；Bilinear StrainHardening/Softening Model双线性应变硬化/软化模型；Modified CamClay Model 修正剑桥模型；DoubleYield Model 胡克布朗模型。,各向同性弹性模型：剪切模量；体积模量；密度,摩尔-库伦模型：剪切模量；体积模量；密度；摩擦角；粘结力；剪胀角(任选)。横观各向同性模型：剪切模量；x-模量；y-模量；密度；NUYx；NUYz。砌体节理模型：剪切模量；体积模量；密度；粘结力(整体材料)；摩擦力(整体材料)；节理粘结力；节理摩擦力；节理角；整体材料的剪胀角(任选)。德鲁克-普拉格模型：剪切模量；体积模量；密度；材料参数（Kshear；材料参数（Qdil）；材料参数（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）本构模型需要的参数没有指定时，系统使用默认值，一般为0。,FLAC3D本构模型一例子 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 ini 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模型），一般用于研究盐矿中核废料的地下储藏的有关热力学分析，也适合于软土的变形特性分析；伯格（Burger）蠕变模型和摩尔库伦模型合成的伯格蠕变粘塑性模型，它是在经典的粘弹性模型的基础上扩展而成的，它包括一个开尔文体和一个摩尔库伦体；WIPP模型和德鲁克普拉格（Drucker-Prager）模型合成的WIPP蠕变粘塑性模型，它是WIPP模型的变化形式，包括一个德鲁克普拉格（Drucker-Prager）塑性体；岩盐的本构模型，它也是WIPP模型的变化形式，它包括体积和偏量压实特性。,FLAC3D结构单元FLAC3D 可以模拟 6种结构单元（梁、锚索、桩及板壳、格栅、衬砌单元）FLAC3D 可以模拟任意形状、任意特性的结构体与岩土体的相互作用，以及力作用在结构体上或岩土体时，结构体和岩土体的力学反应。结构单元力学性态的计算模拟可以通过几何大变形和小变形两种计算模式计算，而且还设置动力状态来模拟结构与岩土体的动力反应。,梁单元需输入的参数density：density（梁单元的密度）,Emod：Youngs modulus（弹性模量）,E nu：Poissons ratio（泊松比）,pmoment：plastic moment capacity（塑性矩）,MP thexp：thermal expansion coefficient（热膨胀系数）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,梁单元的命令输入形式SEL beam（定义梁单元）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&xcarea=6e-3 xcj=0.0 xciy=200e-6 xciz=200e-6,锚杆单元需输入的参数density：density（锚杆密度）,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 ycompression：compressive yield strength（锚杆的抗压屈服强度）,Fc ytension：tensile yield strength（锚杆的拉伸屈服强度）,Ft,锚杆单元的命令输入形式SEL cable（定义锚杆单元）SEL cablesel（定义锚杆单元条件）可能的锚杆单元与已有的锚杆单元的连接是通过id命令以及id的值来确定 锚杆单元的命令一例子 sel cable id=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=25,桩单元需输入的参数density：density（桩单元的密度）,Emod：Youngs modulus（弹性模量）,Enu：Poissons ratio（泊松比）,pmoment：plastic moment capacity（塑性矩）,MP(可选的除非特别说明，假设力矩承载力是足够大的)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轴的惯性矩）,Izxcj：polar moment of inertia（极惯性矩）,J,cs_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 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,壳单元需输入的参数density：density（壳单元的密度）,isotropic：isotropic material properties:E and 和 orthotropic：orthotropic material properties(正交各向异性材料的参数)thexp：thermal-expansion coefficient(热膨胀系数)thickness：thickness(厚度),t,土工格栅单元需输入的参数density：density（单元的密度）,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：coupling 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(容许的大应变滑移值),衬砌单元需输入的参数density：density（桩单元的密度）,isotropic：isotropic material properties:E and 和 orthotropic：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 cohesion(stress units)(切向连接弹簧的残余应力),cr cs_sfric：shear coupling spring friction angle(切向连接弹簧的摩擦角),cs_sk：shear coupling spring stiffness per unit area(切向连接弹簧每单位面积的刚度),ksslide：large-strain sliding flag(default:off)(大应变滑移标记,默认为关)slide tol：large-strain sliding tolerance(容许的大应变滑移值),界面单元接触面参数的确定 虚构的为了合并节点而设置的接触面 Kn=ks=10*真实的刚性接触面 如料仓下料 c,D,Tension重要，kn,ks不重要 真实的柔性接触面 断层；水力劈裂材料 试验得到参数：对于kn,ks：岩石断层10100MPa/m(粘土);100GPa(岩石)反分析方法：通过断层中岩石的变形与原岩的变形。,建立interface方法“移来移去”(推荐)。关键要形成同一位置的两个节点(面)；建两个分开的模型；建立接触单元；通过INI*add使模型接触；注意dist的含义。“导来导去”。利用expgrid,impgrid命令进行网格导出与导入，配合DELETE命令，适于内部接触面的建立，或在已建立好的网格中加入接触面。例子（第一种方法）,建立接触面之前的几何形状,加进接触面单元,最终的几何网格形状,命令流：gen zone brick size 3 3 3 p0(0,0,0)p1(3.0,0,0)p2(0,3,0)p3(0,0,1.5)p4(3.0,3.0,0)p5(0,3.0,1.5)Lower top to complete geometryini z add-1.0 range group Topsave interface.sav,4 离散单元法,概述离散单元法是20世纪70年代初兴起的一种数值计算方法，适合节理岩体的应力分析离散元法也将模型划分成刚性单元，单元之间可以相互叠合，也可以相互分离单元之间相互作用的力可以根据力和位移的关系求出，而个别单元的运动法则完全根据该单元所受的不平衡力与不平衡力矩的大小按牛顿运动定律确定单元之间不需满足变形协调方程,离散单元法的基本方程物理方程-力和位移的关系,运动方程-牛顿第二运动定律,合力：,合力矩：,加速度：,角加速度：,加速度、速度、位移关系：,离散单元法的计算机实现力和位移的计算循环,已知力求位移,考虑阻尼的运动方程：,对运动方程进行一阶中心差分：,由上式可求得 时刻的速度与角速度：,则从 到 时刻的线位移与角位移增量为：,则 时刻的线位移与角位移为：,已知位移求力,进行坐标变换得接触点处法向与切向位移增量：,则接触点上的作用力为：,法向力,切向力,法向阻尼力,切向阻尼力,将接触点上的作用力转化到形心上：,所有接触点上的作用力转化到形心上：,UDEC简介通用离散元程序（UDEC）是一个处理不连续介质的二维离散元程序，非连续介质是通过离散的块体集合体加以表示，可模拟非连续介质承受静载或动载作用下的响应。不连续面处理为块体间的边界面，允许块体沿不连续面发生较大位移和转动。块体可以是刚体或变形体。变形块体被划分成有限个单元网格，且每一单元根据给定的“应力应变”准则，表现为线性或非线性特性。不连续面发生法向和切向的相对运动也由线性或非线性“力位移”的关系控制。允许离散块体发生有限的位移和转动，包括完全脱离，在计算过程中，自动识别新的接触面。,UDEC特性UDEC开发了人工或自动节理生成器，用以模拟产生岩体中一组或多组不连续面。地质图上所显示的节理为直线段，节理的几何粗糙度也可以通过节理材料模型加以表征。UDEC的基本公式假设为二维平面应变模型。对于动态计算，用户指定的速度或应力波可作为外部的边界条件或者内部激励直接输入到模型中。在静态分析中，包括了应力（力）和固定位移（速度为零）两种边界条件。UDEC还能够模拟通过模型中的孔隙和不连续面的流体流动。流体被处理为平行板的粘性流。程序中的结构单元可用于模拟岩体加固和工程表面支护。加固包括端部锚固、全长锚索和锚杆。表面支护模拟诸如喷射混凝土、混凝土衬砌和其他形式的隧道支护。,UDEC应用领域UDEC常用于采矿工程，已经进行了深部地下采矿洞室的静态与动态分析。洞室围岩破坏诱发的断裂、滑移是用UDEC分析研究的实例之一。UDEC还应用于地下结构和深部高辐射废料的储存研究领域。通过应用热模型，UDEC已经应用于模拟与核废料相关的热荷载效应。UDEC最初是为节理岩石边坡的稳定性分析开发的。对于块体不连续公式和运动方程（包括惯性项）采用显式时间步求解方法，便于块状岩体边坡的渐进破坏分析和大变形运动研究。UDEC程序对于模拟颗粒流动或动态分析火山喷发是不适宜的。对于该类研究，可以采用PFC2D程序。,图 UDEC模型的例子,UDEC模型：初始块体的划分UDEC模型首先生成整个计算范围的单一块体。然后，通过用地质结构特征和工程结构作为边界，切割该块体成小的块体来考虑模型特征。模型的所有块体都是通过块体质心和角点的坐标（x和y）确定。块体接触面以及变形块体的节点也通过他们的坐标位置确定。产生模型包括由端点坐标（x，y坐标）所定义的线段（splits）切割模型块体。例子 block(0,0)(0,10)(10,10)(10,0)；顺时针 Crack(0,5)(10,5)plot hold block num cont num dnum,UDEC模型：初始块体的划分1Block x1,y1 x2,y2 x3,y3；在此，（x1，y1），（x2，y2），（x3，y3）是定义块体角点的坐标对。角点必须按顺时针方向排列。round d；在此，d是圆角距离（缺省值是d=0.5）。模型中的所有圆角长度都是相同的。其目的在于避免块体悬挂在有棱角的节点上。由于块体悬挂引起应力集中。对于变形块体，最大圆角长度应当不超过块体平均棱长的1。生成地质结构的两个主要命令如下：Crack和Jset，Crack 命令用于产生块体中单一直线特征的裂缝。裂缝由端点坐标（x1，y1）和（x2，y2）所确定。Jset 命令则是自动节理组生成器。根据所给定的特征参数（即倾角、迹长、岩桥长度、间距和空间位置）产生一组裂缝。注意：UDEC需要连续断裂，刚性块体在计算过程中，或变形块体在单元划分时，没有连接形成完整块体的裂缝将被删除。,UDEC模型：初始块体的划分2例子1 round 0.1 block(0,0)(0,10)(10,10)(10,0)Crack(0,5)(10,5)crack 2.5,10 5.0,7.5 crack 5.0,7.5 7.5,10 Plot block num delete range 4.5,5.5 8,10；4.5 x 5.5 和8 y 10的范围必须包含被删除块体的形心。例子2 new round 0.01 block(0,0)(0,20)(20,20)(20,0)jset(45,0)(5,0.5)(0.5,0)(2,0)；倾角逆时针为正 jset(-10,0)(5,0.5)(0.2,0)(1.5,0),UDEC模型：初始块体的划分3通过切割UDE