三维拉格朗日法计算原理

1三维快速拉格朗日法的基本原理1.1概述目前在岩土力学中常用的数值计算方法有差分方法、有限元法、边界元法等几种，特别是后两种方法，随着计算机的发展其应用尤为广泛。但是，这几种方法都是以连续介质为出发点，而且往往囿于小变形的假定。它们虽然也可以用来解决由几种介质所组成的非均质的问题，并且对于个别的断层或弱面，也可以用设置节理单元的办法来解决，但是用以解决富含节理和大变形的岩土力学问题，往往所得的结果与实际的物理图景相差甚远。于是离散单元法和拉格朗日元法就应运而生。离散单元法是Cundall于上世纪70年代初所提出的。该法将为弱面所切割的岩体视为复杂的块体的集合体，允许各个块体可以平移或转动，甚至相互分离。拉格朗日元法则是由Cundall所加盟的美国ITASCA咨询集团于1986年所开发的。该法将流体力学中跟踪流体运动的拉格朗日方法应用于解决岩体力学的问题获得成功。三维快速拉格朗日法是一种基于三维显式有限差分法的数值分析方法，它可以模拟岩土或其他材料的三维力学行为。三维快速拉格朗日分析将计算区域划分为若干四面体单元，每个单元在给定的边界条件下遵循指定的线性或非线性本构关系，如果单元应力使得材料屈服或产生塑性流动，则单元网格可以随着材料的变形而变形，这就是所谓的拉格朗日算法，这种算法非常适合于模拟大变形问题。三维快速拉格朗日分析采用了显式有限差分格式来求解场的控制微分方程，并应用了混合单元离散模型，可以准确地模拟材料的屈服、塑性流动、软化直至大变形，尤其在材料的弹塑性分析、大变形分析以及模拟施工过程等领域有其独到的优点。1.2三维快速拉格朗日分析的数学模型三维快速拉格朗日分析在求解中使用如下3种计算方法：(1)离散模型方法。连续介质被离散为若干六面体单元，作用力均被集中在节点上。(2)有限差分方法。变量关于空间和时间的一阶导数均用有限差分来近似。(3)动态松驰方法。由质点运动方程求解，通过阻尼使系统运动衰减至平衡状态。1.2.1空间导数的有限差分近似快速拉格朗日分析采用混合离散方法，将区域离散为常应变六面体单元的集合体，又将每个六面体看作以六面体角点为角点的常应变四面体的集合体，应力、应变、节点不平衡力等变量均在四面体上进行计算，六面体单元的应力、应变取值为其内四面体的体积加权平均。这种方法既避免了常应变六面体单元常会遇到的位移剪切锁死现象，又使得四面体单元的位移模式可以充分适应一些本构的要求，如不可压缩塑性流动等。如一四面体,节点编号为1到4，第n面表示与节点n相对的面，设其内一点的速率分量为vi,由高斯公式得：dSnvdVvSjiVji,(1-1)其中V为四面体的体积，S为四面体的外表面，nj为外表面的单位法向向量分量。对于常应变单元，vi为线性分布，nj在每个面上为常量。对式(1-1)积分得：)()()(41,ffjffijiSnvVv（1-2）式中，上标（f）指面f的相关变量值，iv指i速度分量的均值。若速度呈线性变化，则：4,1)()(31flllifivv（1-3）上标l指节点l的值。将上式代入式（1-2），有：)(4,1)(41,31flfffjllijiSnvVv（1-4）在式（1-1）中，若vi=1，应用高斯法则可得：0)(41)(fffjSn（1-5）所以，式（1-4）两边同除以V，则有：图1-1四面体单元的面和节点)(41)(,31llljlijiSnvVv（1-6）而应变速率张量则可由下式表示：应变速率张量的分量形式为：)(41)()()(61llliljljliijSnvnvV(1-7)1.2.2节点运动方程一定时域内，静力平衡问题可通过以下的平衡方程求解得到：0,ijijB(1-8)式中：ρ为介质密度，dtdvbBiii，bi为介质单位质量的体积力。根据虚功原理，作用于单个四面体上的节点力fl（l=（1,4））与四面体应力和等效体力相平衡。引入节点虚速度δvl（它在四面体中产生线性速度场δv和常应变速率δξ），则节点力Fl和体力B产生的外力功功率等于内部应力σij产生的内力功功率。外力功功率可表示为：dVBvfvEinVinini41（1-9）而内力功功率：dVIijVij（1-10）由式（1-7），对常应变速率的四面体有：)(4161llliijljljijliSnvnvI（1-11）应力张量是对称张量，定义矢量Tl：)(lljijliSnT（1-12）则：4131lliliTvI（1-13）式（1-8）代入式（1-9），有：IbninniEEfvE41（1-14）Eb和EI分别为体力ib和惯性力所作的外力功功率。若四面体内体力ib为常数，则有：ViibdVvbE（1-15）ViiIdVdtdvvE（1-16）根据有限差分近似，速度场在四面体内线性变化。为描述它，引进一个参考坐标系（它的坐标原点则四面体的中心上），则有：41nnniiNvv（1-17）式中Nn（n=1，4）为一线性函数：3322110xcxcxccNnnnnn（1-18）其中，nnnncccc3210,,,（n=1，4）为下述方程的解：njjjjnxxxN321,,（1-19）式中，nj是克罗内克尔增量（Kroneckerdelta）。通过中心点的定义，所有形如VjdVx'的积分均为0，将式（1-18）、式（1-17）代入式（1-14）得：410nnniibVcvbE（1-20）由克雷姆定律，解式(1-19)得：410nc（1-21）将上式代入式（1-20），有：414niVnibbvE（1-22）同理，将式（1-17）代入式（1-16）得到：41nVinniIdVdtdvNvE（1-23）将式（1-22）和（1-23）代入式（1-14）：414nVinininidVdtdvNVbfvE（1-24）对任何虚速度，外虚功率E等于内虚功率I：VinininidVdtdvNVbTf43（1-25）在四面体范围内，加速度场空间变化是微小的，则有：VnniVindVNdtdvdVdtdvN（1-26）为不变量，则上式可写为：niVindtdvVdVdtdvN4（1-27）用假想的节点质量mn代替上式中的质量4V：则，式（1-25）可写为：ninininidtdvmVbTf43（1-28）对于等效体系，可以建立平衡状态，要求在每个节点上静态等效载荷之和为零。可以写出全部节点上牛顿定律表达式：nlillinldtdvMF,1（1-29）式中，nn介质中的所有的节点总数，节点质量定义为：llmM（1-30）不平衡力[F]l定义为：liliiliPVbTF43（1-31）当介质达到平衡时，不平衡力等于0。1.2.3增量形式的本构方程快速拉格朗日分析中，假定时间t内速度为常数，增量形式的本构方程可表示为：),(*tHijijijij(1-32)式中，ij称为共转（co-rotational）应力增量，*ijH为一给定的函数。共转（co-rotational）应力速率张量ij等于给定参考系的介质内一点应力的偏导数和以瞬时角速度的转动，数学表达式为：kjijijijijijwwdtd][(1-33)式中，w为转动速率张量。利用有限差分方程，可以得到转动速率张量的分量形式：)(41)()(61llliljljliijSnvnvVw(1-34)式中符合同前。1.2.4时间导数的有限差分近似由本构方程（式（1-32））和变形速率与节点速率之间的关系（式（1-7）），式（1-26）可表示为一般的差分方程：nlpiiiilillinlkvvvvtFMdtdv,1,,....,,,,1321（1-35）式中，{}l是指在计算过程中全局节点l节点速度值的子集（式（1-29））。在时间间隔t中实际节点的速度假定是线性变化的，式(1-35)左边导数用中心有限差分估算。kvvvvtFMttvttvlpiiiilillili,,....,,,,1)2()2(321（1-36）类似地，节点的位置也用中心有限差分进行迭代：)2()()(tttvtxttxlilili（1-37）因此，节点位移也有如下关系：)2()()(tttvtuttulilili（1-38）1.2.5阻尼力为使运动方程获得静态或准静态（非惯性）解，快速拉格朗日分析的静力分析中，在式(1-29)中加入非粘性阻尼力。则式(1-29)变为：nldtdvMFlillili,1(1-39)式中：lil为阻尼力，lililivsignF，为阻尼系数，其默认值为08。000,10,1)(yifyifyifysign(1-40)1.3FLAC3D简介由以上原理可以看出，无论是动态问题，还是静态问题，三维快速拉格朗日分析均由运动方程用显式方法进行求解，这使得它很容易模拟动态问题，如振动、失稳、大变形等。对显式法来说非线性本构关系与线性本构关系并无算法上的差别，对于已知的应变增量，可很方便地求出应力增量，并得到不平衡力，就同实际中的物理过程一样，可以跟踪系统的演化过程。此外，显式法不形成刚度矩阵，每一步计算所需计算机内存很小，使用较少的计算机内存就可以模拟大量的单元，特别适于在微机上操作。在求解大变形过程中，因每一时步变形很小，可采用小变形本构关系，只需将各时步的变形叠加，即得到了大变形。这就避免了通常大变形问题中推导大变形本构关系及其应用中所遇到的麻烦，也使它的求解过程与小变形问题一样。根据前述原理，美国ItascaConsultingGroup开发了三维快速拉格朗日分析程序FLAC一3D，该程序能较好地模拟地质材料在达到强度极限或屈服极限时发生的破坏或塑性流动的力学行为，特别适用于分析渐进破坏和失稳以及模拟大变形。它主要有如下一些特点:(l)应用范围广泛，可以模拟复杂的岩土工程或力学问题。FLAC3D包含了10种弹塑性材料本构模型，有静力、动力、蠕变、渗流、温度五种计算模式，各种模式间可以互相藕合，以模拟各种复杂的工程力学行为。FLAC-3D可以模拟多种结构形式，如岩体、土体或其他材料实体，梁、锚元、桩、壳以及人工结构如支护、衬砌、锚索、岩栓、土工织物、摩擦桩、板桩等，另外，FLAC3D设有界面单元，可以模拟节理、断层或虚拟的物理边界等；(2)FLAC3D具有强大的内嵌程序语言FISH，使得用户可以定义新的变量或函数，以适应用户的特殊需要。例如，利用FISH，用户自己设计FLAC3D内部没有的特殊单元形态；用户可以在数值试验中进行伺服控制；可以指定特殊的边界条件，自动进行参数分析；可以获得计算过程中节点、单元参数，如坐标、位移、速度、材料参数、应力、应夺、不平衡力等；(3)FLAC3D具有强大的前后处理功能。FLAC3D具有强大的自动三维网格生成器，内部定义了多种基本单元形态，可以生成非常复杂的三维网格。在计算过程中用户可以用高分辨率的彩色或灰度图或数据文件输出结果，以对结果进行实时分析，图形可以表示网格、结构以及有关变量的等值线图、矢量图、曲线圈等，可以给出计算域的任意截面上的变量等值线图和矢量图。FLAC3D具有如下缺陷：(1)对于线性问题，FLAC3D要比相应的有限元花费更多的计算时间，FLAC3D在模拟非线性问题、大变形问题或动态问题时更有效。(2)FLAC3D的收敛速度取决于系统的最大固有周期与最小固有周期的比值，这使