三维拉格朗日法计算原理

1 三维快速拉格朗日法的基本原理1.1 概述目前在岩土力学中常用的数值计算方法有差分方法、 有限元法、 边界元法等 几种，特别是后两种方法，随着计算机的发展其应用尤为广泛。但是，这几种方 法都是以连续介质为出发点， 而且往往囿于小变形的假定。 它们虽然也可以用来 解决由几种介质所组成的非均质的问题， 并且对于个别的断层或弱面， 也可以用 设置节理单元的办法来解决，但是用以解决富含节理和大变形的岩土力学问题， 往往所得的结果与实际的物理图景相差甚远。 于是离散单元法和拉格朗日元法就 应运而生。离散单元法是 Cundall 于上世纪 70 年代初所提出的。该法将为弱面所切割 的岩体视为复杂的块体的集合体， 允许各个块体可以平移或转动， 甚至相互分离。 拉格朗日元法则是由 Cundall 所加盟的美国 ITASCA 咨询集团于 1986 年所开发 的。该法将流体力学中跟踪流体运动的拉格朗日方法应用于解决岩体力学的问题 获得成功。三维快速拉格朗日法是一种基于三维显式有限差分法的数值分析方法， 它可 以模拟岩土或其他材料的三维力学行为。 三维快速拉格朗日分析将计算区域划分 为若干四面体单元， 每个单元在给定的边界条件下遵循指定的线性或非线性本构 关系，如果单元应力使得材料屈服或产生塑性流动， 则单元网格可以随着材料的 变形而变形，这就是所谓的拉格朗日算法， 这种算法非常适合于模拟大变形问题。 三维快速拉格朗日分析采用了显式有限差分格式来求解场的控制微分方程， 并应 用了混合单元离散模型， 可以准确地模拟材料的屈服、 塑性流动、 软化直至大变 形，尤其在材料的弹塑性分析、 大变形分析以及模拟施工过程等领域有其独到的 优点。1.2 三维快速拉格朗日分析的数学模型三维快速拉格朗日分析在求解中使用如下 3 种计算方法： (1)离散模型方法。 连续介质被离散为若干六面体单元，作用力均被集中在节点上。(2) 有限差分方法。变量关于空间和时间的一阶导数均用有限差分来近似。 (3)动态松驰方法。 由质点运动方程求解，通过阻尼使系统运动衰减至平衡状态。1.2.1 空间导数的有限差分近似快速拉格朗日分析采用混合离散方法， 将区域离散为常应变六面体单元的集 合体，又将每个六面体看作以六面体角点为角点的常应变四面体的集合体， 应力、 应变、节点不平衡力等变量均在四面体上进行计算， 六面体单元的应力、 应变取 值为其内四面体的体积加权平均。 这种方法既避免了常应变六面体单元常会遇到 的位移剪切锁死现象， 又使得四面体单元的位移模式可以充分适应一些本构的要 求，如不可压缩塑性流动等。如一四面体 ,节点编号为 1 到 4，第 n 面表示与节点 n 相对的面， 设其内一点V vi,j dVSvinjdS的速率分量为 vi,由高斯公式得：（1-1）其中 V 为四面体的体积， S 为四面体的外表面， nj 为外表面的单位法向向量分量。 对于常应变单元， vi 为线性分布， nj 在每个面上为常量。对式 （1-1）积分得：图 1-1 四面体单元的面和节点4 (f )Vvi,jvinj(f)S( f)(1-2)f1式中，上标（ f）指面 f 的相关变量值，vi指 i速度分量的均值。若速度呈线性变化，则：vi(f)(l) vi1-3)1-4)1-5)l 1,l f上标 l 指节点 l 的值。将上式代入式( 1-2)，有：1 4 41 l (f) ( f) Vvi,jvinj S3l1 f 1,f l在式（ 1-1）中，若 vi=1，应用高斯法则可得：4n(jf)S( f)0f1所以，式（ 1-4）两边同除以 V，则有：141-6)1 l (l) (l ) vi,j vi nj S3V l 1而应变速率张量则可由下式表示：应变速率张量的分量形式为：4ij1 l (l) l (l) (l ) (vin jvjni )S6V l 1(1-7)1.2.2 节点运动方程一定时域内，静力平衡问题可通过以下的平衡方程求解得到：ij ,jBi0（1-8）式中： 为介质密度， Bibi dvi ，bi 为介质单位质量的体积力。i i dt i根据虚功原理，作用于单个四面体上的节点力 fl（ l=（1,4）与四面体应力和等效体力相平衡。引入节点虚速度vl（它在四面体中产生线性速度场v 和常应变速率 ），则节点力 Fl 和体力B 产生的外力功功率等于内部应力ij 产生的内力功功率。外力功功率可表示为：4Evin finvi BidVn 1 V1-9)而内力功功率：I ij ij dV（1-10）由式（ 1-7），对常应变速率的四面体有：14I 1vil ijnlj vlj ijnilS（l）6l11-11）应力张量是对称张量，定义矢量 Tl ：Tilijnlj S(l )iij j（1-12）则：1 4 l l IvilTil3 l 1（1-13）式（ 1-8）代入式（ 1-9），有：4n n b I Evi f i E En1（1-14）Eb和 EI分别为体力 bi和惯性力所作的外力功功率。 若四面体内体力 bi为常数，则有：Eb bi vidV（1-15）E I vi dvi dVV i dt（1-16）根据有限差分近似， 速度场在四面体内线性变化。 为描述它， 引进一个参考坐标系（它的坐标原点则四面体的中心上） ，则有：4vivinN nn1（1-17）式中 Nn（n=1 ，4）为一线性函数：N n c0n c1nx1 c2n x2 c3n x3（1-18）其中， c0n,c1n,c2n,c3n （n=1，4）为下述方程的解：Nn x1j ,x2j ,x3jnj（1-19）式中， nj 是克罗内克尔增量（ Kronecker delta）。通过中心点的定义，所有形如 xj dV的积分均为 0，将式（ 1-18）、式（ 1-17）代入式（ 1-14）得：Eb4bivinc0nV1-20)n1由克雷姆定律，解式 （1-19）得：nc01-21)将上式代入式（ 1-20），有：E bvin biV41-22)n1同理，将式（ 1-17）代入式（ 1-16）得到：n1n dvidt dV1-23)将式（ 1-22）和（ 1-23）代入式（ 1-14）：4nEvin fin1nbiV4V Nn ddvti dV1-24)对任何虚速度，外虚功率 E 等于内虚功率I：finTi nbiVNndvi3 4 V dt dV1-25)在四面体范围内，加速度场空间变化是微小的，则有：V Nn ddvti dVddvtiV NndV1-26)为不变量，则上式可写为：V Nn ddvti dV4V ddvti1-27)用假想的节点质量 mn 代替上式中的质量则，式（ 1-25）可写为：TinbiV mn4ndvidt1-28)对于等效体系，可以建立平衡状态， 要求在每个节点上静态等效载荷之和为 零。可以写出全部节点上牛顿定律表达式：FilMldvi ldtl 1,nn(1-29)式中， nn介质中的所有的节点总数，节点质量定义为：1-30)M l m l不平衡力F 定义为：lTibiV l l1-31)Fi liiPi li34i当介质达到平衡时，不平衡力等于 01.2.3 增量形式的本构方程 快速拉格朗日分析中，假定时间 t 内速度为常数，增量形式的本构方程可 表示为：ijHi*j(,ij ijt)(1-32)式中， ij 称为共转( co-rotational)应力增量， Hi*j 为一给定的函数。共转co-rotational)应力速率张量 ij 等于给定参考系的介质内一点应力的偏导数和以瞬时角速度 的转动，数学表达式为： ijd ijdtwij ij ijwkj(1-33)式中，w 为转动速率张量。利用有限差分方程，可以得到转动速率张量的分量形式：41 l (l) vi n j 6Vvj ni(l)S(l)l1(1-34)式中符合同前。1.2.4 时间导数的有限差分近似由本构方程(式( 1-32)和变形速率与节点速率之间的关系(式( 1-7)， 式( 1-26)可表示为一般的差分方程：dvi 1lFil t, vi 1 ,vi 2,vi 3 , vi p,k l 1,nn(1-35)dtM li i ii i n式中， 是指在计算过程中全局节点 l 节点速度值的子集(式( 1-29)。 在时间间隔 t 中实际节点的速度假定是线性变化的， 式(1-35)左边导数用中心有 限差分估算。vil (t2t)vil(t2t)M1lFi lt, vi1 ,vi2,vi3,vipl ,k(1-36)类似地，节点的位置也用中心有限差分进行迭代：xi l (t t) xi l (t) tvi l (t t)(1-37)因此，节点位移也有如下关系：ui l (t t) ui l (t) tvi l (t t)(1-38)21.2.5 阻尼力 为使运动方程获得静态或准静态 (非惯性) 解，快速拉格朗日分析的静力分 析中，在式 (1-29)中加入非粘性阻尼力。则式 (1-29)变为：Fil li M l dvil 1,ni i dt(1-39)式中： lil 为阻尼力， liFilsignvil ， 为阻尼系数，其默认值为 0 8。1, ify0sign( y)1, ify0 (1-40)0 ify01.3 FLAC 3D 简介 由以上原理可以看出，无论是动态问题，还是静态问题，三维快速拉格朗日 分析均由运动方程用显式方法进行求解， 这使得它很容易模拟动态问题， 如振动、 失稳、大变形等。对显式法来说非线性本构关系与线性本构关系并无算法上的差 别，对于已知的应变增量，可很方便地求出应力增量，并得到不平衡力，就同实 际中的物理过程一样， 可以跟踪系统的演化过程。 此外，显式法不形成刚度矩阵， 每一步计算所需计算机内存很小，使用较少的计算机内存就可以模拟大量的单 元，特别适于在微机上操作。在求解大变形过程中，因每一时步变形很小，可采 用小变形本构关系， 只需将各时步的变形叠加， 即得到了大变形。 这就避免了通 常大变形问题中推导大变形本构关系及其应用中所遇到的麻烦， 也使它的求解过 程与小变形问题一样。根据前述原理，美国 Itasca Consulting Group开发了三维快速拉格朗日分析程 序 FLAC 一 3D ，该程序能较好地模拟地质材料在达到强度极限或屈服极限时发 生的破坏或塑性流动的力学行为， 特别适用于分析渐进破坏和失稳以及模拟大变 形。它主要有如下一些特点 :(l) 应用范围广泛，可以模拟复杂的岩土工程或力学问题。 FLAC 3D 包含了 10 种弹塑性材料本构模型，有静力、动力、蠕变、渗流、温度五种计算模式，各种 模式间可以互相藕合， 以模拟各种复杂的工程力学行为。 FLAC - 3D 可以模拟多 种结构形式，如岩体、土体或其他材料实体，梁、锚元、桩、壳以及人工结构如 支护、衬砌、锚索、岩栓、土工织物、摩擦桩、板桩等，另外， FLAC 3D 设有界 面单元，可以模拟节理、断层或虚拟的物理边界等；(2) FLAC 3D 具有强大的内嵌程序语言 FISH，使得用户可以定义新的变量或函 数，以适应用户的特殊需要。例如，利用 FISH，用户自己设计 FLA C 3D内部没 有的特殊单元形态； 用户可以在数值试验中进行伺服控制； 可以指定特殊的边界 条件，自动进行参数分析； 可以获得计算过程中节点、 单元参数， 如坐标、位移、 速度、材料参数、应力、应夺、不平衡力等；(3) FLAC 3D具有强大的前后处理功能。 FLAC 3D具有强大的自动三维网格生成 器，内部定义了多种基本单元形态， 可以生成非常复杂的三维网格。 在计算过程 中用户可以用高分辨率的彩色或灰度图或数据文件输出结果， 以对结果进行实时 分析，图形可以表示网格、结构以及有关变量的等值线图、矢量图、曲线圈等， 可以给出计算域的任意截面上的变量等值线图和矢量图。FLAC 3D 具有如下缺陷：(1)对于线性问题， FLAC 3D 要比相应的有限元花费更多的计算时间， FLAC 3D在模拟非线性问题、大变形问题或动态问题时更有效(2)FLAC 3D 的收敛速度取决于系统的最大固有周期与最小固有周期的比值， 这使得它对某些问题的模拟效率非常低， 如单元尺寸或材料弹性模量相差很大的 情况。1.4 本构模型FLAC 3D 包含了 10种弹塑性材料本构模型， 以下简单介绍报告中涉及到的 4 类本构模型。1.4.1 空单元模型( Null Model ) 空单元材料用于描述从模型中删除或开挖掉的部分。在模拟的后续阶段，空 单元材料可以被转变成不同的材料模型。 用这种方式， 例如，能够模拟回填开挖。 在空单元区域内的所有应力被自动设置成零。(1-41)1.4.2 Mohr-Coulomb 模型( Mohr-Coulomb Model )1.4.2.1 组合的破坏准则和流动法则在这个模型的破坏准则是张拉剪切组合的 Mohr-Coulomb 准则。这个准则可 以用图 (1-2)解释。用 Mohr-Coulomb 破坏准则描绘从点 A 到点 B 破坏包络线f s 0 ，即：f s 1 3N 2c N(1-42)用式 f t 0 张拉破坏准则描绘从 B 点到点 C 的包络线：t是张拉强度，且有：(1-43)式中， 是摩擦角， C 是粘聚力，1 s in ()1 s in ()(1-44) 张拉强度不超过 3 值。最大值由下式给定：1-45)tcmaxt an分别用两个定义剪切塑性流动和张拉塑性流动的函数 gs和 gt 描述势函数。函数 gs 有如下形式：g 1 3N1-46)式中， 是膨胀角。N 1 sin( )N1 sin( )(1-47)函数 gt 符合相关流动法则，写成：tg3(1-48)图 1-2 Mohr-Coulomb 破坏准则用式 h( 1, 3) 0 将流动法则写成统一的形式：h3PPt aP ( 1P)(1-50)式中， aP和 P 是由下式定义的常数：1-51)tN 2c N1-52)1.4.3 遍布节理模型( Ubiquitous-Joint Model)这个模型描述了 Mohr-Coulomb 模型中弱面的力学行为。弱面的方向给定，弱 面 的 破 坏 准 则 由 组 合 的 Mohr-Coulomb 张 拉 剪 切 破 坏 包 络 线 组 成 。 与Mohr-Coulomb 模型一样，弱面方向上的破坏准则在剪切破坏时采用非相关联的流动法则，在拉张破坏时采用相关联的流动法则。非弱面方向的破坏准则与Mohr-Coulomb 模型一致，如前所述，以下简单说明以下弱面上的破坏准则。在以下各式中， 各应力分量皆是弱面局部坐标系下的应力分量， 局部坐标系规定如下：x轴与弱面的法向一致， y 轴水平， z轴符合右手螺旋法则。FLAC-3D 模型 中弱 面的 破坏准则 是应 力 ( 33, )表示的组 合的Mohr-Coulomb 张拉剪切破坏准则，如图 (1-3) 。从点 A 到点 B 定义的包络线 f s( 33 , ) 0：33 tan j cj(1-52)图 1-3 遍布节理模型的破坏准则t c j jmax tan j用式 f t 0 张拉破坏准则描绘 从 B 点到点 C 的包络线：tft33 j(1-53)式中，j，Cj，jt 分别是弱面 的摩擦角、 粘聚力和张拉强度， 对 于非零摩擦角的弱面， 张拉强度的 最大值为：(1-54)势函数由 gs和 gt 两个函数组成，这两个函数分别定义剪切和张拉塑性流动。 函数 gs 有如下形式：gs33 t an j(1-55)式中，j 是弱面膨胀角。流动法则用统一式表示为：P P th j a j （ 33 j ）（1-56）式中， jP 和 aPj 是由下式定义的常数：Ptj cj t an j j（1-57） a Pj1 t an j2 t a n j（1-58）1.4.4 修正剑桥模型（ Modified Cam-Clay Model ） 土的有效应力弹塑性本构关系， 即应力增量 与应变增量 的关系可用下式 表达： =Cep （159）式中 Cep 为弹塑性矩阵。（1-60）式中Ce为弹性矩阵； f 为屈服函数； g 为势函数。 FLAC 3D 中，采用相关联 流动法则，故取 g= f，A为表示硬化规律的参数。图 1-4为修正剑桥模型，它采 用的是塑性体积应变的硬化规律，相关联的流动法则（即 g =f）和具有椭圆轨迹 的屈服函数，如式 （1-61）所示 :22f (q, p) q2 M 2p( p pc)（1-61）式中， p和 q定义为：1p 3 ii（1-62）q 3J 2（1-63）M 为材料常数。 剑桥理论采用体应变 硬化模式，其屈服函数可 表示为：图 1-4 修正剑桥模型示意图f ( p,q) H ( eP)(1-64)PH ( eP) e xpeP e0 e k ln p ( 1-66) p0相应增量为 dp和 dq时有：ffdf dp dqpq( e p)d( ep )(1-67)p p f 而 d( ep ) (1 e) vp (1 e)d p1-68)故有A df 1 e f f kp图 1-5 等向压缩与膨胀曲线 （图中 v 为比容）1-69)式中 d 为非负值的函数； e 为孔隙比；为等向压缩曲线的斜率； 为等向回弹曲线的斜率，如图 1-5 所示故有：1-70)f 1 q2 11 2 2p M 2 p 21-71)f 2q2 q M p