FLAC3D基本原理

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

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

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

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

（1） 符号约定

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

dvidt，i?1，3来定义一个点的坐标，位移，速度和加速的。记号ai表示矢量?a?的第i个分量，在笛卡尔坐标系中；Aij表示张量?A?的第（i，j）个分量。a，i表示变量对xi的偏导数。（变量a可以使标量，矢量和张量）

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

（2） 应力

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

ti??ijnj

（2.37）

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

（3） 应变率和转动率

假设介质的离子以张量?v?运动。在一个无限短时间dt内，介质产生一个无限小的应变为vidt，相关的应变率张量可以写成如下：

?ij?12?vi,j?vj,i?

（2.38）

第一应变率张量不变量描述了体积单元的的膨胀程度。张量?ij中没有包含变形率，由于速度矢量的平移和角速度的转动，一个体积单元会产生一个瞬间的刚体位移，如下：

?i??eijk12eijk?jk （2.39）

表示置换符号，矢量???表示转动率张量，定义如下：

?ij?12?vi,j?vj,i?

（2.40）

（4） 运动平衡方程

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

?ij,j??bi???dvidt （2.41）

为介质的密度，?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个方程描述介质的应力和应变之间的关系。一般定义如下：

?Hij??ij,?ij,?? （2.43）

????ij????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??Svinjds （2.45）

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

4Vvi,j??vf?1(f)inj(f)S(f) （2.46）

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

vi(f)?134?l?1,l?fvil （2.47）

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

Vvi,j?14li4v??3l?1njS(f)(f) （2.48）

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

4?nf?1(f)jS(f)?0 （2.49）

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

vi,j??13V4?vnl?1li(l)jS(l) （2.50）

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

?ij??16V??vnlil?14(l)j?vjnil(l)?S(l) （2.51）

（2） 节点平衡运动方程

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

?ij,j??Bi?0 （2.52）

体力的定义如公式2.41

dvi??Bi???bi??dt?? （2.53）

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

n力和体力都是通过虚功方程推导出来。假设四面体有一个虚拟的速度??v?（它在四面体内会产生一个线性变化的速度场??v?和一个常应变率????），我们假定外力虚功由节点力?f?和体力?B?产生，内力虚功由虚速度下的应力?ij产生。 外力功率为：

n4E???vn?1nifi?n?V?viBidV （2.54）

内力虚功率为：

I??4V??ij?ijdV （2.55）

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

I??1??v??6lil?1ijnj??vj?ijni(l)l(l)?S(l) （2.56）

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

Ti??ijnjSl(l)(l) （2.57）

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

I??134??viTi （2.58）

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

4E???vn?1nifi?E?EnbI （2.59）

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

面体常体力：

E??bi??vidV （2.60）

VbE?????viVIdvidtdV （2.61）

根据以前的假设，在四面体内速度场是线性变化的。为了描述它，现在定义一个

?,x3?。四面体内任一新的局部坐标系，坐标系原点在四面体形心，坐标轴为x1?,x2点速度可以用四面体4个顶点的速度差值得到，如下：

4?vi?N,n?1,4为线性函数

n?n?1?viN （2.62）

nnNn??c3x3? ?c0?c1x1??c2x2nnnn （2.63）

c0,c1,c2,c3,n?1,4nnnn,为常数，可以通过求解以下方程得出：

百度搜索“77cn”或“免费范文网”即可找到本站免费阅读全部范文。收藏本站方便下次阅读，免费范文网，提供经典小说综合文库FLAC3D基本原理在线全文阅读。