三维拉格朗日法计算原理Word文档格式.docx
《三维拉格朗日法计算原理Word文档格式.docx》由会员分享，可在线阅读，更多相关《三维拉格朗日法计算原理Word文档格式.docx（14页珍藏版）》请在冰豆网上搜索。
对于常应变单元，vi为线性分布，nj在每个面上为常量。
对式（1-1）积
分得：
图1-1四面体单元的面和节点
4（f）
Vvi,jvinj（f）S（f）（1-2）
f1
性变化，则：
1-4）
1-5）
l1,lf
上标l指节点l的值。
将上式代入式（1-2），有：
144
1l（f）（f）Vvi,jvinjS
3l1f1,fl
在式（1-1）中，若vi=1，应用高斯法则可得：
4
n（jf）S（f）0
所以，式（1-4）两边同除以V，则有：
14
1-6）
1l（l）（l）vi,jvinjS
3Vl1
而应变速率张量则可由下式表示：
应变速率张量的分量形式为：
ij
1l（l）l（l）（l）（vinjvjni）S
6Vl1
（1-7）
1.2.2节点运动方程
一定时域内，静力平衡问题可通过以下的平衡方程求解得到：
ij,jBi0
（1-8）
式中：
ρ为介质密度，Bibidvi，bi为介质单位质量的体积力。
iidti
根据虚功原理，作用于单个四面体上的节点力fl（l=（1,4））与四面体应力
和等效体力相平衡。
引入节点虚速度
δvl（它在四面体中产生线性速度场
δv和
常应变速率δξ），则节点力Fl和体力
B产生的外力功功率等于内部应力
σij产生
的内力功功率。
外力功功率可表示为：
EvinfinviBidV
n1V
1-9）
而内力功功率：
IijijdV
（1-10）
由式（1-7），对常应变速率的四面体有：
I1vilijnljvljijnilS（l）
6l1
1-11）
应力张量是对称张量，定义矢量Tl：
TilijnljS（l）
iijj
（1-12）
则：
14llIvilTil
3l1
（1-13）
式（1-8）代入式（1-9），有：
nnbIEvifiEE
n1
（1-14）
Eb和EI分别为体力bi和惯性力所作的外力功功率。
若四面体内体力bi为常数，
则有：
EbbividV
（1-15）
EIvidvidV
Vidt
（1-16）
根据有限差分近似，速度场在四面体内线性变化。
为描述它，引进一个参考
坐标系（它的坐标原点则四面体的中心上），则有：
vivinNn
（1-17）
式中Nn（n=1，4）为一线性函数：
Nnc0nc1nx1c2nx2c3nx3
（1-18）
其中，c0n,c1n,c2n,c3n（n=1，4）为下述方程的解：
Nnx1j,x2j,x3jnj
（1-19）
式中，nj是克罗内克尔增量（Kroneckerdelta）。
通过中心点的定义，所有
形如x'
jdV的积分均为0，将式（1-18）、式（1-17）代入式（1-14）得：
由克雷姆定律，解式（1-19）得：
将上式代入式（1-20），有：
同理，将式（1-17）代入式（1-16）得到：
将式（1-22）和（1-23）代入式（1-14）：
在四面体范围内，加速度场空间变化是微小的，则有：
为不变量，则上式可写为：
用假想的节点质量mn代替上式中的质量
则，式（1-25）可写为：
1-28）
对于等效体系，可以建立平衡状态，要求在每个节点上静态等效载荷之和为零。
可以写出全部节点上牛顿定律表达式：
（1-29）
式中，nn介质中的所有的节点总数，节点质量定义为：
1-30）
Mlml
不平衡力[F]<
l>
定义为：
lTibiVll
1-31）
FiliiPil
i34i
当介质达到平衡时，不平衡力等于0
1.2.3增量形式的本构方程快速拉格朗日分析中，假定时间t内速度为常数，增量形式的本构方程可表示为：
（1-32）
式中，ij称为共转（co-rotational）应力增量，Hi*j为一给定的函数。
共转
co-rotational）应力速率张量ij等于给定参考系的介质内一点应力的偏导数和
以瞬时角速度的转动，数学表达式为：
（1-33）
式中，w为转动速率张量。
利用有限差分方程，可以得到转动速率张量的分量形式：
（1-34）
式中符合同前。
1.2.4时间导数的有限差分近似
由本构方程（式（1-32））和变形速率与节点速率之间的关系（式（1-7）），式（1-26）可表示为一般的差分方程：
dvi1lFilt,vi1,vi2,vi3,,vip,kl1,nn（1-35）
dtMliiiiin
式中，{}<
是指在计算过程中全局节点l节点速度值的子集（式（1-29））。
在时间间隔t中实际节点的速度假定是线性变化的，式（1-35）左边导数用中心有限差分估算。
vil（t2t）vil（t2t）M1lFilt,vi1,vi2,vi3,,vipl,k（1-36）
类似地，节点的位置也用中心有限差分进行迭代：
xil（tt）xil（t）tvil（tt）（1-37）
因此，节点位移也有如下关系：
uil（tt）uil（t）tvil（tt）（1-38）
2
1.2.5阻尼力为使运动方程获得静态或准静态（非惯性）解，快速拉格朗日分析的静力分析中，在式（1-29）中加入非粘性阻尼力。
则式（1-29）变为：
FilliMldvil1,n
iidt
（1-39）
lil为阻尼力，liFilsignvil，为阻尼系数，其默认值为08。
1,ify0
sign（y）1,ify0（1-40）
0ify0
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的收敛速度取决于系统的最大固有周期与最小固有周期的比值，这使得它对某些问题的模拟效率非常低，如单元尺寸或材料弹性模量相差很大的情况。
1.4本构模型
FLAC3D包含了10种弹塑性材料本构模型，以下简单介绍报告中涉及到的4类本构模型。
1.4.1空单元模型（NullModel）空单元材料用于描述从模型中删除或开挖掉的部分。
在模拟的后续阶段，空单元材料可以被转变成不同的材料模型。
用这种方式，例如，能够模拟回填开挖。
在空单元区域内的所有应力被自动设置成零。
（1-41）
1.4.2Mohr-Coulomb模型（Mohr-CoulombModel）
1.4.2.1组合的破坏准则和流动法则
在这个模型的破坏准则是张拉剪切组合的Mohr-Coulomb准则。
这个准则可以用图（1-2）解释。
用Mohr-Coulomb破坏准则描绘从点A到点B破坏包络线
fs0，即：
fs13N2cN
（1-42）
用式ft0张拉破坏准则描绘从B点到点C的包络线：
σt是张拉强度，且有：
（1-43）
式中，φ是摩擦角，C是粘聚力，
1sin（）
（1-44）张拉强度不超过σ3值。
最大值由下式给定：
tc
max
tan
数gs有如下形式：
g13N
1-46）
式中，ψ是膨胀角。
法则写成统一的形式：
（1-50）
式中，aP和P是由下式定义的常数：
1-51）
tN2cN
1-52）
1.4.3遍布节理模型（Ubiquitous-JointModel）
这个模型描述了Mohr-Coulomb模型中弱面的力学行为。
弱面的方向给定，
弱面的破坏准则由组合的Mohr-Coulomb张拉剪切破坏包络线组成。
与
Mohr-Coulomb模型一样，弱面方向上的破坏准则在剪切破坏时采用非相关联的
流动法则，在拉张破坏时采用相关联的流动法则。
非弱面方向的破坏准则与
Mohr-Coulomb模型一致，如前所述，以下简单说明以下弱面上的破坏准则。
在以下各式中，各应力分量皆是弱面局部坐标系下的应力分量，局部坐标系
规定如下：
x'
轴与弱面的法向一致，y'
轴水平，z'
轴符合右手螺旋法则。
FLAC-3D模型中弱面的破坏准则是应力（3'
3'
）表示的组合的
Mohr-Coulomb张拉剪切破坏准则，如图（1-3）。
从点A到点B定义的包络线fs（3'
）0：
用式ft0张拉破坏准则描绘从B点到点C的包络线：
t
ft3'
j
（1-53）
式中，φj，Cj，σjt分别是弱面的摩擦角、粘聚力和张拉强度，对于非零摩擦角的弱面，张拉强度的最大值为：
（1-54）
势函数由gs和gt两个函数组成，这两个函数分别定义剪切和张拉塑性流动。
函数gs有如下形式：
gs3'
tanj
（1-55）
式中，ψj是弱面膨胀角。
流动法则用统一式表示为：
PPt
hjaj（3'
j）
（1-56）
式中，jP和aPj是由下式定义的常数：
Pt
jcjtanjj
（1-57）aPj1tanj2tanj
（1-58）
1.4.4修正剑桥模型（ModifiedCam-ClayModel）土的有效应力弹塑性本构关系，即应力增量'
与应变增量的关系可用下式表达：
{'
}=[Cep]{}
（1－59）
式中[Cep]为弹塑性矩阵。
（1-60）
式中[Ce]为弹性矩阵；
f为屈服函数；
g为势函数。
FLAC3D中，采用相关联流动法则，故取g=f，A为表示硬化规律的参数。
图1-4为修正剑桥模型，它采用的是塑性体积应变的硬化规律，相关联的流动法则（即g=f）和具有椭圆轨迹的屈服函数，如式（1-61）所示:
22
f（q,p）q2M2p（ppc）
ePe0eklnp（1-66）p0
相应增量为dp和dq时有：
向回弹曲线的斜率，如图1-5所示
故有：
1-70）
f1q21
122
pM2p2
1-71）
f2q
2qMp
升级会员 