2、2 三维数值模拟方法及其原理2、2、1 FLAC3D工程分析软件特点FLAC3D就是由美国Itasca Consulting Group, Inc、为地质工程应用而开发得连续介质显式有限差分计算机软件。
FLAC即Fast Lagrangian Analysis of Continua 得缩写。
该软件主要适用于模拟计算岩土体材料得力学行为及岩土材料达到屈服极限后产生得塑性流动,对大变形情况应用效果更好。
FLAC3D程序在数学上采用得就是快速拉格朗日方法,基于显式差分来获得模型全部运动方程与本构方程得步长解,其本构方程由基本应力应变定义及虎克定律导出,运动平衡方程则直接应用了柯西运动方程,该方程由牛顿运动定律导出。
计算模型一般就是由若干不同形状得三维单元体组成,也即剖分得空间单元网络区,计算中又将每个单元体进一步划分成由四个节点构成得四面体,四面体得应力应变只通过四个节点向其它四面体传递,进而传递到其它单元体。
当对某一节点施加荷载后,在某一个微小得时间段内,作用于该点得荷载只对周围得若干节点(相邻节点)有影响。
利用运动方程,根据单元节点得速度变化与时间,可计算出单元之间得相对位移,进而求出单元应变,再利用单元模型得本构方程,可求出单元应力。
在计算应变过程中,利用高斯积分理论,将三维问题转化为二维问题而使其简单化。
在运动方程中,还充分考虑了岩土体所具有得粘滞性,将其视作阻尼附加于方程中。
FLAC3D具有一个功能强大得网格生成器,有12种基本形状得单元体可供选择,利用这12种基本单元体,几乎可以构成任何形状得空间立体模型。
FLAC3D主要就是为地质工程应用而开发得岩土体力学数值评价计算程序,自身设计有九种材料本构模型:(1)空模型(Null Model)(2)弹性各向同性材料模型(Elastic, Isotropic Model)(3)弹性各向异性材料模型(Elastic, anisotropic Model)(4)德拉克-普拉格弹塑性材料模型(Drucker-Prager Model)(5)莫尔-库伦弹塑性材料模型(Mohr-Coulomb Model)(6)应变硬化、软化弹塑性材料模型(Strain-Hardening/Softening Mohr-Coulomb Model)(7)多节理裂隙材料模型(Ubiquitous-Joint Model)(8)双曲型应变硬化、软化多节理裂隙材料模型(Bilinear Strain-Hardening/Softening Ubiquitous-Joint Model)(9)修正得Cam粘土材料模型(Modified Cam-clay Model)除上述本构模型之外,FLAC3D还可进行动力学问题、水力学问题、热力学问题等得数值模拟。
在边界条件及初始条件得考虑上,FLAC3D软件十分灵活方便,可在数值计算过程中随时调整边界条件与初始条件。
FLAC3D具有强大得后处理功能,用户可以直接在屏幕上绘制或以文件形式创建或输出打印多种形式得图形、文字,用户还可根据各自得需要,将若干个变量合并在同一幅图形中进行研究分析。
FLAC3D软件还可对各种开挖工程或施加支护工程等进行数值仿真模拟,软件自身设计有锚杆、锚索、衬砌、支架等结构元素,可以直接模拟这些支护于围岩(土)体得相互作用。
FLAC3D拥有可以自行设计得FISH语言,用户可根据自身需求,自己设计材料得本构模型、屈服准则、支护方案、复杂形状得开挖方式等工作。
特别注意得就是,岩石就是一种脆性材料,当外荷载达到岩石强度后,材料发生断裂破坏,产生弱化现象,应属于弹塑性体。
在FLAC3D中,一般对于弹塑性材料,判断其破坏与否得基本准则有两个,即Drucker-Prager准则与Mohr-Coulomb准则。
根据室内岩石力学性质试验结果,其典型应力应变曲线反映出岩体破坏包络线符合莫尔—库伦屈服准则,故本次建立得本构力学模型选择莫尔—库伦弹塑性材料模型为宜。
2、2、2 FLAC3D分析计算原理计算所采用得数学模型就是根据弹塑性理论得基本原理(应变定义、运动定律、能量守衡定律、平衡方程及理想材料得连续性方程等)而建立得。
2、2、2、1 基本约定在数学及数值模型得表达式中,符号有一定得约定含义,一般[A]表示张量,A ij 表示张量[A ]得(i,j )分量,[a ]表示矢量,a i 表示矢量[a ]得i 分量,α,i 表示α对x i 得偏导数。
x i ,u i ,v i 与dv i /dt ,(i=1,3)分别表示一点得位置矢量分量、位移矢量分量、速度矢量分量与加速度矢量分量。
2、2、2、2 数学模型(一)柯西(Cauchy)应力张量与柯西公式对于一个具有体积V 得封闭曲面s 得物体,在其上取一表面元素∆s ,这个表面元素得单位外法向矢量为n ,在某一时刻t ,在表面元素对于连续介质中一点,作用着对称得应力张量σij ,根据∆s 上作用有力∆P ,则极限dsdPs P T s ==→∆∆∆0lim称为表面力。
若用t i 表示T 得分量,则在三维直角坐标系中可有关系式ij i i n t σ= (1) 这个关系式称为柯西公式,其中,σij 称为柯西应力张量。
(二)应变速率与旋转速率如果介质质点具有运动速度矢量[v ],则在一个无限小得时间dt 内,介质会产生一个由v i dt 决定得无限小应变,对应得应变速率分量ξij 为)(21ijj i ij x v x v ∂∂+∂∂=ξ (2)而其旋转速率分量ωij 为)(21ijj i ij x v x v ∂∂-∂∂=ω (3) (三)运动及平衡方程根据牛顿运动定律与柯西应力原理,如果质点作用着应力σij 与体力b i ,且具有速度v i ,则在无限小时间段dt 内,它们之间得关系为dtdvb x i i j ij ρρσ=+∂∂ (4) 式中,ρ为质点密度。
(4)式称为柯西运动方程。
当质点得加速度为零时,上式变为静力平衡方程0=+∂∂i jijb x ρσ (5) (四)本构方程上述(4)式与(5)式组成得方程组中含有9个方程,15个未知量,其中12个就是应力与应变速率分量,3个就是速度分量。
其余6个关系式则由本构方程提供,本构方程一般具有如下形式),,(][κξσσij ij ij ij H =(6) 式中,][σ为应力变化速率,H 表示一个特定得函数关系,κ为与荷载历史有关得参数。
2、2、2、3 数值模型FLAC 3D 得数值剖分网格在计算中就是按照四面体进行得,四面体得节点也既就是网格剖分得节点,因此,每个计算单元有4个面与4个节点(见图2-1)。
(一)空间微分得有限差分逼近对于一个计算单元,若内部各质点速度为一连续得矢量场[v ],则根据高斯(Gauss)积分原理有⎰⎰=∂∂S j i V j ids n v dv x v (7)式中,[n ]为外法向单位矢量场。
由于单元体得应变速率就是连续得,因此可以近似认为速度就是线性变化得,则(7)式可用下面得求与公式近似逼近∑==∂∂41)()()(f f f j f i j i S n v x v V (8) 单元面单元节点图2-1 四面体单元示意图式中,V 为单元体体积,)(f S 为单元体某一面得面积,f 为单元体面数,f =1,4,)f (iv 为面平均速度得i 分量。
由于速度场就是线性得,则有∑≠==4,1)(31fl l l i f i v v (9)式中,l 为单元体节点数,l=1,4;v i l 为l 节点得i 速度分量。
将(9)式代入(8)式可得∑∑≠===∂∂4,1)()(4131lf f f f jl l ij i S nv x v V (10)根据正交原理有041)()(=∑=f f f jS n(11)则将(10)式两边除以V,并将(11)式代入可得)()(4131l l j l l i j i S n v V x v ∑=-=∂∂ (12) 因此有)(41)()()(61l l l i l j l j l iij S n v nv V∑=+-=ξ (13)(二)运动方程得节点公式根据前面对质点运动方程得讨论,对于连续介质,当处于平衡状态时,其平衡方程为0=+∂∂i jij B x ρσ (14))(dtdv b B ii i -=ρ (15) 介质可由若干个作用着体力[B ],各自产生一定变形得四面体组成,设节点力为[f ]n ,n =1,4,利用虚功原理,假定单元体节点具有速度δ[v ]n ,内部具有变形速率δ[ξ],则由节点力[f ]n 与体力[B ]所做得外功率与σij 所做得内功率应当相等。
外功率E 可用下式表示dV B v f v E i Vi n i n n i ⎰∑+==δδ41(16)而内功率I 为⎰=Vij ij dV I σδξ (17)由(13)式,对于恒定应变速率单元体可有()()()∑=+-=41)(61i l l i ij l j l j ij l i S n v n v I σδσδ (18)由于应力张量就是对称得,定义矢量Tl()()l l j ij l i S n T σ= (19) 可得∑=-=4131l li l i T v I δ (20)将(15)式代入(16)式得∑=++=41n I b n i n i E E f v E δ (21)式中,E b 与E I 分别就是体力ρb i 与惯性力产生得外功率,对于单元体内恒定得体力ρb i ,E b 可写成⎰=Vi i b dV v b E δρ (22)而E I 可以写成⎰-=ViiI dV dtdv v E ρδ (23) 如前所述,单元体内部速度场以线性变化,为方便描述,选取单元体质心为原点,1x '、2x '、3x '为坐标轴,可有 ∑==41n n n i i N v v δδ (24)式中,N n (n =1,4)就是具有如下形式得线性函数3322110x c x c x c c N nn n n n '+'+'+= (25) nc 0、n c 1、nc 2、nc 3(n=1,4)则由下面方程决定()nj j j j n x x x N δ='''321,, (26) δnj 为克罗内克尔(Kronecker)记号,根据质心定义可得∑==410n n n i i bV c v b E δρ (27) 用克莱默(Cramer ’s)法则解得41=nc ,则∑==414n i ni bVb v E ρδ (28)而∑⎰=-=41n Vinn i IdV dtdv N v E ρδ (29) 则∑⎰=⎥⎦⎤⎢⎣⎡-+=414n V i n i n i n i dV dt dv N V b f v E ρρδ (30)在稳定状态下,内功率与外功率必定相等,则有⎰-+=-V i n i n i ni dV dtdvN V b T f ρρ43 (31)上式最后一项,如果单元体内ρ恒定,则根据质心定义有ni Vi ndt dv V dV dt dv N )(4ρρ=⎰(32) 若将ρV /4瞧作就是假定得节点质量m n ,则上式变为n i n Vi ndtdvm dV dt dv N )(=⎰ρ (33) 故n i n i n i nidtdv m V b T f )(43-+=-ρ (34) 根据牛顿定律有n l i l l i n l dtdv M F ,1 )(==><><>< (35) 式中,F i <l >为节点l 所受力得i 分量,M <l >为节点l 得质量,><l i dtdv )(为节点l 在F i <l >作用下产生得加速度,n n 为包含在全部连续介质中单元节点总数。