FLAC3D基本原理

2.5 三维显示有限差分基本方程当FLAC3D 达到平衡或是稳定的塑性流动时，它通过显示有限差分来模拟三维连续介质的力学行为。

监控的力学响应主要是通过特殊的数学模型和数值计算过程得到。

接下来介绍这两方面。

2.5.1 数学模型描述介质的力学行为主要来源于一般原理（应变定义、运动规律），和理想材料的本构关系。

这个数学结果表达式通常是一些偏微分方程，涉及到力学（应力）和运动学（应变率、速度）变量。

这些偏微分方程联合个别的几何关系、材料参数，以及给定的边界条件和初始条件就可以求解。

虽然FLAC3D 在平衡状态附近，主要关注介质的应力状态和变形，但是必须要注意到该数学模型中的运动方程。

（1） 符号约定在FLAC 3D 中采用拉格朗日算法，介质中的一个点，通过矢量i i i x u v ，，和13i dv dt i =，，来定义一个点的坐标，位移，速度和加速的。

记号i a 表示矢量[]a 的第i 个分量，在笛卡尔坐标系中；ij A 表示张量[]A 的第（i ，j ）个分量。

i a ，表示变量对i x 的偏导数。

（变量a 可以使标量，矢量和张量）默认结构受拉为正，变形伸长为正。

爱因斯坦求和记号只针对下标，i ，j ，k （i ，j ，k =1，2，3）。

（2） 应力介质中一已知点的应力状态是通过对称应力张量ij σ来表示。

任意斜面上的应力矢量[]t 可以通过柯西公式得到（拉为正），如下：i ij j t n σ= （2.37）[]n 表示任意斜面上的单位法向矢量（3） 应变率和转动率假设介质的离子以张量[]v 运动。

在一个无限短时间dt 内，介质产生一个无限小的应变为i v dt ，相关的应变率张量可以写成如下：(),,12ij i j j i v v ξ=+ （2.38） 第一应变率张量不变量描述了体积单元的的膨胀程度。

张量ij ξ中没有包含变形率，由于速度矢量的平移和角速度的转动，一个体积单元会产生一个瞬间的刚体位移，如下：12i ijk jk e ωΩ=- （2.39）ijk e 表示置换符号，矢量[]ω表示转动率张量，定义如下：(),,12ij i j j i v v ω=- （2.40） （4） 运动平衡方程采用连续介质的动量原理和柯西公式，平衡方程如下：,i ij j i dvb dtσρρ+= （2.41）ρ为介质的密度，[]b 表示单位体力，[]d v dt 表示速度矢量对时间的导数。

当力开始施加到介质上，平衡方程贯穿在整个的数学模型中，以及单元介质的运动中。

在静力计算过程中，[]d v dt 为0，公式2.41简化为如下偏微分方程：,0ij j i b σρ+= （2.42）（5） 边界条件和初始条件应力边界条件主要是通过公式2.37来表示，位移边界主要是通过指定边界的速度分量为0来实现。

初始条件中体力也是可以施加的，需要指定初始应力状态。

（6） 本构方程公式2.41和公式2.38中包含有9个方程，15个未知量，这15个未知量是6个应力分量和6个应变分量，以及3个速度分量。

通过本构方程可以提供额外的6个方程，这个6个方程描述介质的应力和应变之间的关系。

一般定义如下：[]()ˆ,,ij ij ij ij H σσξκ= （2.43） []ˆij σ为共轭应力张量，[]H 为一已知函数，κ为考虑加载历史变量， []ˆij ik kj ik kj ij d dtσσωσσω=-+ （2.44）这里[]d dt σ表示应力矢量[]σ对时间的实导数，[]ω表示转动率张量。

2.5.2 数值方程FLAC3D 通过以下三步骤来求解：（1）有限差分逼近（变量的一阶导数、时间导数用有限差分来逼近，假定变量在很短的空间内和时间间隔内线性变化）（2）离散逼近（将连续介质离散为与之相当的网格，在这个网格中，所有的力包括外力和内力，都作用在单元节点的三个方向上）（3）动态求解方法（在平衡方程中引入惯性定律，使得系统慢慢达到平衡） 连续介质的运动定律通过以上三步骤，转化为离散单元节点上的牛顿定律。

一般普通的微分方程可以通过时间显示差分求解。

介质的偏微分平衡方程中涉及到的空间导数，出现在以速度来定义的应变率中。

为了进一步的定义速度变量和相关的空间间隔，连续介质被划分为常应变率的四面体单元，这个四面体单元的顶点就是网格单元的顶点。

如图2.2所示。

图2.2 四面体单元（1） 空间导数的有限差分形式四面体的应变率张量分量的有限差分方程是通过节点平衡方程得到。

四面体共有1-4个节点，节点n 所对的面，称之为面n 。

通过高斯散度定律可得如下方程：,i j i j VSv dV v n ds =⎰⎰ （2.45）通过高斯散度定律将四面体上的体积分转化为四个面上的面积分。

假设四面体是常应变率，所以速度矢量是线性变化的，每个面的法向方向也是常量，公式2.45简化为：4()()(),1f f f i j ij f Vv v n S ==∑ （2.46）上标（f ）表示面f ，i v 表示变量i v 的平均，假设速度线性变化，可以得出：4()1,13f lii l l fv v =≠=∑ （2.47） 这里的上标l 指的是节点。

将公式2.47带入2.46可得：44()(),11,13l f f i j i j l f f lVv v n S ==≠=∑∑ （2.48）假设在公式2.45中i v =constant ，我们得到散度定理如下：4()()10f f jf nS ==∑ （2.49）利用时2.49可以将2.48化简为：4()(),113l l l i j i jl v v nS V==-∑ （2.50）因此应变率张量就可以表示为：()4()()()116l l l l l ij ijj i l v nv n S Vξ==-+∑ （2.51）（2） 节点平衡运动方程节点运动平衡方程通过虚功原理来推导，在任一个瞬时，得到一个等价静态问题。

通过在节点惯性平衡方程中引入虚功原理，来求解整个结构。

固定时间t ，我们通过平衡方程来研究静力等价状态问题,0ij j i B σρ+= （2.52）体力的定义如公式2.41ii i dv B b dt ρ⎛⎫=- ⎪⎝⎭（2.53） 在有限差分的框架中，介质由一些连续的承受体力[]B 的常应变四面体单元代替。

在静态平衡方程中，作用在单个四面体上的节点力[],1,4nf n =，以及四面体应力和体力都是通过虚功方程推导出来。

假设四面体有一个虚拟的速度[]nv δ（它在四面体内会产生一个线性变化的速度场[]v δ和一个常应变率[]δξ），我们假定外力虚功由节点力[]nf 和体力[]B 产生，内力虚功由虚速度下的应力ij σ产生。

外力功率为：41n n i i i i Vn E v f v B dV δδ==+∑⎰ （2.54）内力虚功率为：ij ij VI dV δξσ=⎰ （2.55）利用公式2.51，公式2.55可以写成常应变率的形式：()4()()()116l l l l l i ij j j ij i l I v n v n S δσδσ==-+∑ （2.56）由于应力张量是对称张量，现在定义一个矢量l T 如下:()()l l l i ij j T n S σ= （2.57）由式2.57，式2.56可写成如下：4113l l i i l I v T δ==-∑ （2.58）将公式2.53带入公式2.54中，可得：41n n b I i i n E v f E E δ==++∑ （2.59）b E 为体力i b ρ对外力虚功率的贡献，I E 为惯性力对虚功率的贡献。

对于一个四面体常体力：b i i VE b v dV ρδ=⎰ （2.60）I iiVdv E v dV dtρδ=-⎰ （2.61） 根据以前的假设，在四面体内速度场是线性变化的。

为了描述它，现在定义一个新的局部坐标系，坐标系原点在四面体形心，坐标轴为123,,x x x '''。

四面体内任一点速度可以用四面体4个顶点的速度差值得到，如下：41n n i i n v v N δδ==∑ （2.62）,1,4n N n =为线性函数0112233n n n n nN c c x c x c x '''=+++ （2.63） 0123,,,,1,4n n n nc c c c n =,为常数，可以通过求解以下方程得出：()123,,n j j j nj N x x x δ'''= （2.64） nj δ为克罗内克符号，通过定义质心，这项0j Vx dV '=⎰，将公式2.62和2.63带入2.60中得：401bn ni i n E b v c V ρδ==∑ （2.65） 由式2.64和形心的性质可知014nc =，带入式2.65中，得： 414bnii n bV E v ρδ==∑ （2.66）将式2.62带入2.61中，得：41In nii Vn dv E v N dV dtδρ==-∑⎰ （2.67） 最后将式2.66和2.67带入2.59中可得：414n n n i i i i V n bV dv E v f N dV dt ρδρ=⎡⎤=+-⎢⎥⎣⎦∑⎰ （2.68）在四面体的静力平衡方程中，在任何虚速度下，外力虚功率等于内力虚功率：由式2.58和2.68可知：34n nn i i i i V T bV dvf N dV dtρρ-=+-⎰ （2.69）假定在四面体内，加速度是不变的，则式2.69中最后一项可改写为：nn n ii VVdv dv N dV N dV dt dt ρρ⎛⎫= ⎪⎝⎭⎰⎰（2.70）又因为ρ为常数，0j Vx dV '=⎰，014nc =，则式2.7为： 4nn ii Vdv dv V N dV dt dt ρρ⎛⎫= ⎪⎝⎭⎰ （2.71） 定义4V ρ为虚拟的节点质量n m ，节点质量可由式2.72确定，节点质量可以保证在平衡的过程中数值计算的稳定。

式2.71改写为：nn n i i Vdv dv N dV m dt dt ρ⎛⎫= ⎪⎝⎭⎰ （2.72） 因此式2.69改写为：34nn nn i i i i T bV dv f m dt ρ⎛⎫-=+- ⎪⎝⎭（2.73）至此等价系统的平衡方程已建立，如式2.73，在每个节点，总的静态平衡力[]f -，包括所有四面体的贡献、节点贡献荷载[]p 以及集中力为0。

为了描述以上规律，我们约定如下记号：如果一变量带有<l >上标，就表示节点变量，<l >表示该节点在全局节点编号中的编号。