FLAC3D基本原理

2.5 三维显示有限差分基本方程

当FLAC3D达到平衡或是稳定的塑性流动时，它通过显示有限差分来模拟三维连续介质的力学行为。监控的力学响应主要是通过特殊的数学模型和数值计算过程得到。接下来介绍这两方面。

2.5.1 数学模型描述

介质的力学行为主要来源于一般原理（应变定义、运动规律），和理想材料的本构关系。这个数学结果表达式通常是一些偏微分方程，涉及到力学（应力）和运动学（应变率、速度）变量。这些偏微分方程联合个别的几何关系、材料参数，以及给定的边界条件和初始条件就可以求解。

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

（1） 符号约定

3在FLAC3D中采用拉格朗日算法，介质中的一个点，通过矢量xi，ui，vi和dvidt，i1，来

定义一个点的坐标，位移，速度和加速的。记号ai表示矢量的第i个分量，在笛卡尔坐

a标系中；

Aij表示张量A的第（i，j）个分量。a，i表示变量对xi的偏导数。（变量a可以使

标量，矢量和张量）

默认结构受拉为正，变形伸长为正。爱因斯坦求和记号只针对下标，i，j，k（i，j，k=1，2，3）。

（2） 应力

来表示。任意斜面上的应力矢量t介质中一已知点的应力状态是通过对称应力张量可以通过柯西公式得到（拉为正），如下：

tiijnjij （2.37）

n表示任意斜面上的单位法向矢量

（3） 应变率和转动率

假设介质的离子以张量运动。在一个无限短时间dt内，介质产生一个无限小的应变

v为vidt，相关的应变率张量可以写成如下：

1vi,jvj,i2 （2.38）

ij第一应变率张量不变量描述了体积单元的的膨胀程度。张量

ij中没有包含变形率，由

于速度矢量的平移和角速度的转动，一个体积单元会产生一个瞬间的刚体位移，如下：

1ieijkjk2 （2.39）

eijk表示置换符号，矢量表示转动率张量，定义如下：

ij1vi,jvj,i2 （2.40）

（4） 运动平衡方程

采用连续介质的动量原理和柯西公式，平衡方程如下：

dvidt （2.41）

ij,jbi为介质的密度，b表示单位体力，dvdt表示速度矢量对时间的导数。当力开始施

加到介质上，平衡方程贯穿在整个的数学模型中，以及单元介质的运动中。在静力计算过程中，

dvdt为0，公式2.41简化为如下偏微分方程：

ij,jbi0 （2.42）

（5） 边界条件和初始条件

应力边界条件主要是通过公式2.37来表示，位移边界主要是通过指定边界的速度分量为0来实现。初始条件中体力也是可以施加的，需要指定初始应力状态。

（6） 本构方程

公式2.41和公式2.38中包含有9个方程，15个未知量，这15个未知量是6个应力分量和6个应变分量，以及3个速度分量。通过本构方程可以提供额外的6个方程，这个6个方程描述介质的应力和应变之间的关系。一般定义如下：

ˆijHijij,ij, （2.43）

ˆijH为共轭应力张量，为一已知函数，为考虑加载历史变量，

ˆijdijdtikkjikkj （2.44）

这里

ddt表示应力矢量对时间的实导数，表示转动率张量。

2.5.2 数值方程

FLAC3D通过以下三步骤来求解：

（1）有限差分逼近（变量的一阶导数、时间导数用有限差分来逼近，假定变量在很短的空间内和时间间隔内线性变化）

（2）离散逼近（将连续介质离散为与之相当的网格，在这个网格中，所有的力包括外力和内力，都作用在单元节点的三个方向上）

（3）动态求解方法（在平衡方程中引入惯性定律，使得系统慢慢达到平衡）

连续介质的运动定律通过以上三步骤，转化为离散单元节点上的牛顿定律。一般普通的微分方程可以通过时间显示差分求解。

介质的偏微分平衡方程中涉及到的空间导数，出现在以速度来定义的应变率中。为了进一步的定义速度变量和相关的空间间隔，连续介质被划分为常应变率的四面体单元，这个四面体单元的顶点就是网格单元的顶点。如图2.2所示。

图2.2 四面体单元

（1） 空间导数的有限差分形式

四面体的应变率张量分量的有限差分方程是通过节点平衡方程得到。四面体共有1-4个节点，节点n所对的面，称之为面n。

通过高斯散度定律可得如下方程：

vVi,jdVvinjdsS （2.45）

通过高斯散度定律将四面体上的体积分转化为四个面上的面积分。假设四面体是常应变率，所以速度矢量是线性变化的，每个面的法向方向也是常量，公式2.45简化为：

Vvi,jvif14(f)nj(f)S(f) （2.46）

上标（f）表示面f，vi表示变量vi的平均，假设速度线性变化，可以得出：

vi(f)14lvi3l1,lf （2.47）

这里的上标l指的是节点。将公式2.47带入2.46可得：

14l4Vvi,jvin(jf)S(f)3l1f1,fl （2.48）

假设在公式2.45中vi=constant，我们得到散度定理如下：

nf14(f)jS(f)0 （2.49）

利用时2.49可以将2.48化简为：

1vi,j3Vvnlil14(l)jS(l) （2.50）

因此应变率张量就可以表示为：

1ij6Vvnlil14(l)jvljni(l)S(l) （2.51）

（2） 节点平衡运动方程

节点运动平衡方程通过虚功原理来推导，在任一个瞬时，得到一个等价静态问题。通过在节点惯性平衡方程中引入虚功原理，来求解整个结构。

固定时间t，我们通过平衡方程来研究静力等价状态问题

ij,jBi0 （2.52）

体力的定义如公式2.41

dvBibiidt （2.53）

在有限差分的框架中，介质由一些连续的承受体力f态平衡方程中，作用在单个四面体上的节点力nB的常应变四面体单元代替。在静

,n1,4，以及四面体应力和体力都是通过

虚功方程推导出来。假设四面体有一个虚拟的速度化的速度场

vvn（它在四面体内会产生一个线性变

n和一个常应变率

ijBf），我们假定外力虚功由节点力和体力产生，内

力虚功由虚速度下的应力产生。

外力功率为：

EvinfinviBidVn1V4 （2.54）

内力虚功率为：

IijijdVV （2.55）

利用公式2.51，公式2.55可以写成常应变率的形式：

14Ivilijn(jl)vljijni(l)S(l)6l1 （2.56）

由于应力张量是对称张量，现在定义一个矢量T如下:

lTilijn(jl)S(l) （2.57）

由式2.57，式2.56可写成如下：

14IvilTil3l1 （2.58）

将公式2.53带入公式2.54中，可得：

EvinfinEbEIn14 （2.59）

Eb为体力bi对外力虚功率的贡献，EI为惯性力对虚功率的贡献。对于一个四面体常

体力：

EbbividVV （2.60）

EIviVdvidVdt （2.61）

根据以前的假设，在四面体内速度场是线性变化的。为了描述它，现在定义一个新的局部坐标系，坐标系原点在四面体形心，坐标轴为x1,x2,x3。四面体内任一点速度可以用四



面体4个顶点的速度差值得到，如下：

vivinNnn14 （2.62）

Nn,n1,4为线性函数

nnnc2c3Nnc0c1nx1x2x3 （2.63）

nnnc0,c1n,c2,c3,n1,4,为常数，可以通过求解以下方程得出：

j,x2j,x3jnjNnx1 （2.64）

nj为克罗内克符号，通过定义质心，这项VxjdV0，将公式2.62和2.63带入2.60

中得：

nEbivinc0Vbn14 （2.65）

由式2.64和形心的性质可知

nc014，带入式2.65中，得：

bVi4Evinbn14 （2.66）

将式2.62带入2.61中，得：

4EvinNnIn1VdvidVdt （2.67）

最后将式2.66和2.67带入2.59中可得：

bVdvEvinfiniNnidVV4dt （2.68） n14在四面体的静力平衡方程中，在任何虚速度下，外力虚功率等于内力虚功率：由式2.58和2.68可知：

TinbVdvfiiNnidVV34dt （2.69）

n假定在四面体内，加速度是不变的，则式2.69中最后一项可改写为：

dvindviNdVdtVdtnVNndV （2.70）

nxjdV0c0V又因为为常数，，14，则式2.7为：

dviVNdVVdt4ndvidt （2.71）

nn定义V4为虚拟的节点质量m，节点质量可由式2.72确定，节点质量可以保证在平

衡的过程中数值计算的稳定。式2.71改写为：

ndvindviNdVmdtVdt （2.72）

n因此式2.69改写为：

TinbVdvfiimni34dt （2.73）

nn至此等价系统的平衡方程已建立，如式2.73，在每个节点，总的静态平衡力括所有四面体的贡献、节点贡献荷载pf，包

以及集中力为0。为了描述以上规律，我们约定如

下记号：如果一变量带有上标，就表示节点变量，表示该节点在全局节点编号中的

.编号。l表示与节点l相连的所有四面体对节点l的贡献之和。采用以上记号，节点牛

顿运动定律可以写成如下形式：

FildvMlidtl l1,nn （2.74）

nn为介质的节点总数，节点质量Ml定义如下：

Mlml （2.75）

F不平衡力定义如下：

lFilTbVii43lPil （2.76）

当系统达到平衡时，不平衡力渐渐趋于0。

（3） 时间导数的显示有限差分

考虑本构方程2.43和变形率和节点速度之间的关系2.51，公式2.74可以写成如下形式：

dvidtl1l123Ft,v,v,v,iiiilM,vipl, l1,nn （2.77）

l{}记号表示计算节点l速度的速度子集。在FLAC3D中求解2.77一般采用对时间t

的显示有限差分。在这个方法中假定节点速度在时间间隔t内线性变化。公式2.77中左边的导数采用中心差分来代替。中心差分时，当计算位移、力时只存取一半时间步的速度。节点速度通过以下迭代公式计算：

ttt)vil(t)lFilt,vi1,vi2,vi3,22Mvil(t,vipl, （2.78）

同理节点坐标更新也是通过中心差分的形式得到：

t)2 （2.79）

xil(tt)xil(t)tvil(t从式2.78和2.79中可知，一阶误差项不存在，因此上述迭代方法是二阶精度的节点位移迭代公式如下：

uil(tt)uil(t)tvil(tt)2 （2.80）

uil(0)0

（4） 本构方程的增量形式

在FLAC3D里面假定在时间间隔t，速度为常量，因此本构方程2.43可以写成如下形式：

*ˆijHijij,ijt （2.81）

ˆij为共轭应力增量，

*Hij为一已知函数。

假定位移在时间t内，线性变化：

ijtij （2.82）

ij为时间t时的应变改变量。

ˆij应力增量通过计算得到，如下：

Cˆijijij （2.83）

Cij为一应力修正，定义如下：

Cijikkjikkjt （2.84）

转动率张量由公式2.40或是2.51计算

（5） 大小应变模式

以上描述的微分方程都是描述大变形，涉及到大位移，以及位移的线性变化和转动。在FLAC3D里面称之为大变形模式。

ijij在实际应用中，转动足够小，以至于分量与整体相比很小，可以忽略不计。张

I量可以用来代替，公式2.83中的应力修正项可以忽略。当然在小位移和位移变化中，

应变率张量中涉及的空间导数，可以通过初始值估计。节点坐标也不再由2.79式更新。

在FLAC3D里面小变形假定小位移、位移梯度和转动，节点坐标不再更新，以及不考虑应力转动修正。

（6） 保持数值稳定的力学时间步

公式2.78能求解出正确解答的前提是数值计算的稳定。通过一个理想的质点-弹簧系统，可以得出一些观点。质点-弹簧系统的平衡方程如下（矩阵符号表示）：

dv*PKuMdt （2.85）

PKM这里大括号表示节点矢量，表示外力，表示弹簧系统的刚度矩阵，表示分

*布质量矩阵。如果我们能解释公式2.74中的不平衡力作为外力施加到弹簧系统上，以及公

式2.85中弹簧的反力，那么理想材料的近似就是可取的。

在采用有限差分来分析质点-弹簧振荡系统时，时间步不能超过临界时间步，而这临界时间步与整个系统的最小固有周期有关。因此在有限差分中必须提供时间的一个上限，以此来保证数值计算的稳定。

推导临界时间步的关系，需要一些固有周期方面的知识。然而在实际中，全局的特征值分析是不可行的，也不经常采用这种方法来得到固有周期。在FLAC3D里面采用在稳定系统中的局部扰动方法来实现一个目标。在整个数值计算过程中，整个系统都是采用一不变的单位时间步t1。

假定公式2.74中右边节点质量为变量，并且参与整个局部平衡过程。

图3 质点弹簧系统1

首先考虑一维自由度弹簧系统，见图3.质点在给定初始位移情况下的运动方程如下：

d2xkxm2dt （2.86）

k为弹簧的劲度，m节点质量。这个方程的二阶有限差分临界时间步如下：

Tt （2.87）

这里T为系统周期。

T2mk （2.88）

图4 质点弹簧系统2

现在假定无穷个质点和弹簧见图4，由于结构对称，这个组合系统的力学行为可以采用图4所示简化方法分析。图4中（a）所示系统可以简化为一刚度矩阵为4k的单一弹簧系统。从式2.87和2.88中可推出极限稳定临界公式如下：

mkt2 （2.89）

假定t1，如果节点质量大于或是等于弹簧劲度，则系统式稳定的。在局部分析中，公式2.89的有效性，通过解释m作为节点l的节点质量贡献m，以及k和相关的节点刚

l

度贡献k扩展到四面体。节点质量贡献通过无穷的准则来提供一个时间步的上限来导出。节点刚度通过局部劲度矩阵的简单对角化得到。

l四面体在节点l的内力贡献为：

Til3，这个力可以假定弹簧的恢复力

llkijuj（式2.85）。

假定在时间间隔dt内，可得如下公式：

dTilllkijvjdt3 （2.90）

将公式2.57带入2.90中可得：

dij3lln(jl)S(l)kijvjdt （2.91）

假定在节点l的q方向上有一单位速度，其他方向上的速度为0，则从2.91式可得：

dqj3kqqdtn(jl)S(l) （2.92）

式2.92中重复指标q并不累计求和。引入小应变的增量虎克定律来描述应力应变之间的关系，式2.92可得：

qq1dt3kqqdt(l)(l)nqS （2.93）

式2.93中1K43G，K为体积模量，G为剪切模量。由于速度已经确定，则由式2.51可得：

qq1(l)(l)nqS3V （2.94）

将式2.94带入2.93中得：

1kqqnS9V （2.95）

(l)q(l)2现在定义节点刚度矩阵的上限为：

klmax(k11,k22,k33) （2.96）

从式2.89可得四面体的节点质量贡献为：

1ml(l)(l)maxnS,i1,3i9V （2.97）

2节点质量是为了保证数值计算的稳定。

（7） 力学阻尼

平衡方程需要提供阻尼，以得到静态解或是准静态解（非惯性）。在FLAC3D中局部非粘性阻尼为默认的静力计算阻尼。在动力分析中也是可以应用局部阻尼的。

在FLAC3D中求解均匀的稳定的运动时需要选择一种阻尼。这种情况发生在蠕变计算和桩的极限承载力计算中。这个阻尼称之为组合阻尼。组合阻尼在消减动能方面比局部阻尼效率要高。

在动力计算中一般采用瑞利阻尼和人工粘性阻尼。

局部阻尼：

在平衡运动方程2.74中施加局部阻尼力，则平衡方程变为：

lFilf(i)dvMlidtl l1,nn （2.98）

lf(i)为阻尼力，通过普遍的平衡力

Fil和速度

vil定义如下：

lllf(Fsignvi)ii1,y0;sign(y)1,y0;0,y0 （2.99）

这个阻尼力通过阻尼常数来控制，默认为0.8。

局部阻尼具有如下特点：

1、 只有存在加速度时才会产生阻尼力，稳态运动是不会产生阻尼力。

2、 阻尼常数为一无量纲量。

3、 阻尼力是与频率无关的，在一个集合体中，不同的固有周期，采用相同的阻尼常数。

局部阻尼与滞后阻尼有点相似，每次循环结束时能量损失都是独立的，或是不同的。

由式2.99可知，阻尼力是与与总合成力成比例的。而不像粘性阻尼是与速度的大小成比例的。当

Flli和

vi同号时，2.98可写为：

dvliFli(1)dtMl 当

Flli和

vi异号时，2.98可写为：

dvliFli(1)dtMl 以上两式可以写成表观质量的形式如下：

dvllliFidviFldtm或是dtimMllmM式2.102中

1,m1 2.100）

2.101）

2.102）

（ （ （

图5 单自由度弹簧系统的运动

单自由度的弹簧系统的运动图3，当给定一个初始位移，在t=0，x=a见图5所示。像这样的运动，在每个循环中，当速度为0时，表观质量增加两次；当速度到达最大值时，其减小两次。阻尼力在每个循环中消减动能两次，在两个速度的峰值，部分质量被删减。注意：加速度在整个过程中是连续的，当质量被移除时，加速度瞬时值为0。

每个循环过程中移除的能量是通过放弃B点的动能或是：

21W(i)2m(i)m(i)v(il)2 （2.103）

瞬间移除的平均动能为：

1l2W(i)m(i)m(i)v(i)4 （2.104）

每个循环中最大储存能量的损失率称为单位能耗，在FLAC3D中，单位能耗可以采用阻尼常数来表述。对一个单自由度系统，或是单模态的振动，动能的峰值与储存能量的峰

值一样，因此单位能耗可由下式表达：

4m(i)m(i)W(i)W(i)m(i)m(i)4 （2.105）

临界阻尼比D可以在小的阻尼中，建立与单位能耗之间的关系：

DW(i)W(i)4 （2.106）

的默认值为0.8，所以D的默认值为0.25

公式2.105可以通过单个单元或是一列三个单元在突然施加重力荷载的情况下，系统振动到平衡来证实。两个实例都表明WW值都很接近，说明局部阻尼式与频率无关的。然而对一个多模态同时作用下的系统，局部阻尼对每个模态都采用相同的阻尼力。

混合阻尼

式2.99中考虑的阻尼力，仅仅考虑当速度分量的符号改变时才改变。在实际，有一重要的匀速运动（与需要衰减的振动大小相比），速度可能没有0点，因此能量并不消耗。通过一个单自由度系统承受谐波振动，可以推导出一可选方程（低效率）来考虑上述问题。这个系统的不平衡力是与速度的负值成比例。例如：

vasint （2.107）

d2v2asint2dt （2.108）

因此dFdt是与速度v成比例的。在一个振荡系统中，dFdt的平均值为0，即使速度v的平均值不为0。在FLAC3D里面采用一种交替形式的阻尼，称之为混合阻尼来实现这个目的。阻尼力由速度和dFdt两项来构成，如下：

fl(i)Fil2dFillsign-signvidt （2.109）

这种阻尼适合于一振荡耗散运动中，有一个刚体运动的系统。然而混合阻尼在能量消散方面要比基于速度的局部阻尼效率低。因此局部阻尼依然是首先的阻尼。

2.5.3 网格离散

在三维常应变率单元中，四面体有不会产生沙漏变形的优点（由节点速度产生的变形形式，没有应变率，因此没有节点力增量）。当在塑性变形的框架下，这些单元不能够提供足够的变形模式。在某些特殊情况下，当某些重要的本构方程，要求体积没有变化时，这时他们不能个别的变形。在这些情况，单元响应比理论预期值要显得过渡坚硬。在FLAC3D中，采用一种混合离散技术来处理上述问题。

混合离散技术的原理是通过适当调整四面体第一应变率张量不变量，使得单元体积具有更大的灵活性。在这个方法中，单元中最粗糙的离散是采用四面体单元来离散，单元中某个四面体的第一应变率张量不变量做为整个单元中所有四面体的体积平均值。这个方法见图6。在草图的特殊变形模式中，个体的常应变率将会产生一个体积改变，而这与不可压缩的塑性流动理论不符。在这个例子中，集合四面体（长方体）的体积仍然为常数，采

用混合离散技术能使每个四面体都能反映单元这个性质，因此将会与理论预见值比较相符。

图6 混合离散技术的单元变形

在FLAC3D中，一个单元是由nt四面体组成，如图7所示nt=5。假定一个单元，假定四面体l的应变率单元式首先考虑的，然后将他们分解为偏量和静水压力项：

l3ijijllij （2.110）

ll式2.110中应变率张量偏量，为第一应变率张量不变量。

llii （2.111）

图7 8节点单元两种离散方法（每种五个四面体）

第一应变率不变量，由每个四面体对体积的平均值得到：

zVV （2.112）

ntkkk1ntkk1V这里为四面体k的体积。因此得应变率张量的计算式如下：

kijijllz3ij （2.113）

当屈服发生时，膨胀本构模型将会导致平均正应力的变化。同理，由应变率增量张量导出的第一应力不变量张量也需要在整个单元上进行体积平均。在这过程中，四面体l的应力张量也是分解为偏张量和静水压力项：

lllijsijij （2.114）

s为应变率张量，为平均正应力。

llllii13 （2.115）

单元的第一应力不变量，通过在所有四面体中进行体积平均得到：

zVV （2.116）

ntkkk1ntkk1由上2.116，式2.114可写为：

llijsijzij （2.117）

在FLAC3D中，离散的过程在初始网格形成时：当单元被定义时,在内部就把其离散化为四面体。例如一个8节点单元，可以离散为两种不同的5个四面体的组合体，见图7。节点应力的计算（基于应变率和应力的计算）可以使用其中一种方法，或是两种方法的一个组合。两种离散方法的优点在于，对称单元在对称荷载作用下，产生对应的响应。在这中情况下，就采用这两种方法结合的混合离散方法。节点力由这两种方法计算出来进行平均得到。

2.5.4 数值计算流程

2.5.4.1网格离散

在FLAC3D中一般由用户将模拟区域进行网格离散，每个单元将会由程序自动离散为

一些列四面体单元。用户能够选择是采用一种离散方法还是采用混合离散方法。一般在应力和位移变化比较剧烈的地方推荐采用混合离散方法。程序默认采用混合离散的方法。

2.5.4.2初始条件很边界条件

求解问题的边界条件包括面力、集中荷载以及位移边界。另外体力和初始应力条件也是需要施加的。程序开始时，所有的应力和节点速度都为0，然后初始应力开始施加。集中荷载指定在表面的节点上，位移边界条件是由节点的速度来精确控制的。体力和面力在内部被转化成一些列等效的节点荷载。以上构成了数值计算的初始状态。

2.5.4.3主要的计算步骤

FLAC3D采用显示的时间历程有限差分，对每一时间步，主要的计算过程如下：

（1） 由节点速度计算新的应变率；

（2） 通过应变率，采用本构方程计算新的上一步的应力；

（3） 由应力和力，通过运动方程导出新的节点速度和位移。

每个时间步都重复上述三个步骤，在求解的过程中监测最大不平衡力。这个力要不趋于0，表示整个系统达到平衡状态，要不等于一非0常数，表明部分或是整个系统达到一个材料的稳定塑性流动。为了分析计算结果，计算过程可以随时打断。

2.5.4.4应变率计算

从一个已知的速度场，开始计算应变率，对单元中的每个四面体来说，使用以下有限差分方程来计算：

1ij6Vvnlil14(l)jvljni(l)S(l) （2.118）

采用混合离散技术以后，通过式2.112和2.113计算新的对角应变率张量。对一个四面体l来说，公式如下：

ijijllz3ij （2.119）

z为单元的第一应变率不变量平均值，计算公式如下：

zVV （2.120）

ntkkk1ntkk1nt为单元中四面体的组合数，V为四面体k的体积。

k2.5.4.5应力计算

在一个单元中，通过增量本构方程

*Hij来计算每个四面体的应力增量。

Cˆijijij （2.121）

*ˆHijij,ijij （2.122）

tij6Vvnlil14(l)jvljni(l)S(l) （2.123）

在小变形模式中，不考虑应力修正项

Cij，在大变形中，需要修正：

Cijikkjikkjt （2.124）

1ij6Vvnlil14(l)jvljni(l)S(l) （2.125）

新的应力通过附加应力增量得到。采用混合离散技术后，通过式2.116和2.117计算新的对角化的应力张量。

llijsijzij （2.126）

zkkk1Vntk1Vntk （2.127）

2.5.4.6节点质量计算

四面体节点l的质量贡献公式如下：

1mmaxni(l)S(l),i1,39V （2.128）

l2注意所有的时间步都是统一的。全局节点质量M2.128求和所得。

l由所有与之共点的四面体采用式

Mlml （2.129）

在小变形模式中，节点质量值在循环开始前只计算一次，在大变形中，没10步计算一次。