FLAC3D基本原理

2.5三维显示有限差分基本方程当FLAC3D达到平衡或是稳定的塑性流动时，它通过显示有限差分来模拟三维连续介质的力学行为。监控的力学响应主要是通过特殊的数学模型和数值计算过程得到。接下来介绍这两方面。2.5.1数学模型描述介质的力学行为主要来源于一般原理（应变定义、运动规律），和理想材料的本构关系。这个数学结果表达式通常是一些偏微分方程，涉及到力学（应力）和运动学（应变率、速度）变量。这些偏微分方程联合个别的几何关系、材料参数，以及给定的边界条件和初始条件就可以求解。虽然FLAC3D在平衡状态附近，主要关注介质的应力状态和变形，但是必须要注意到该数学模型中的运动方程。（1）符号约定在FLAC3D中采用拉格朗日算法，介质中的一个点，通过矢量iiixuv，，和13idvdti，，来定义一个点的坐标，位移，速度和加速的。记号ia表示矢量a的第i个分量，在笛卡尔坐标系中；ijA表示张量A的第（i，j）个分量。ia，表示变量对ix的偏导数。（变量a可以使标量，矢量和张量）默认结构受拉为正，变形伸长为正。爱因斯坦求和记号只针对下标，i，j，k（i，j，k=1，2，3）。（2）应力介质中一已知点的应力状态是通过对称应力张量ij来表示。任意斜面上的应力矢量t可以通过柯西公式得到（拉为正），如下：iijjtn（2.37）n表示任意斜面上的单位法向矢量（3）应变率和转动率假设介质的离子以张量v运动。在一个无限短时间dt内，介质产生一个无限小的应变为ivdt，相关的应变率张量可以写成如下：,,12ijijjivv（2.38）第一应变率张量不变量描述了体积单元的的膨胀程度。张量ij中没有包含变形率，由于速度矢量的平移和角速度的转动，一个体积单元会产生一个瞬间的刚体位移，如下：12iijkjke（2.39）ijke表示置换符号，矢量表示转动率张量，定义如下：,,12ijijjivv（2.40）（4）运动平衡方程采用连续介质的动量原理和柯西公式，平衡方程如下：,iijjidvbdt（2.41）为介质的密度，b表示单位体力，dvdt表示速度矢量对时间的导数。当力开始施加到介质上，平衡方程贯穿在整个的数学模型中，以及单元介质的运动中。在静力计算过程中，dvdt为0，公式2.41简化为如下偏微分方程：,0ijjib（2.42）（5）边界条件和初始条件应力边界条件主要是通过公式2.37来表示，位移边界主要是通过指定边界的速度分量为0来实现。初始条件中体力也是可以施加的，需要指定初始应力状态。（6）本构方程公式2.41和公式2.38中包含有9个方程，15个未知量，这15个未知量是6个应力分量和6个应变分量，以及3个速度分量。通过本构方程可以提供额外的6个方程，这个6个方程描述介质的应力和应变之间的关系。一般定义如下：ˆ,,ijijijijH（2.43）ˆij为共轭应力张量，H为一已知函数，为考虑加载历史变量，ˆijikkjikkjijddt（2.44）这里ddt表示应力矢量对时间的实导数，表示转动率张量。2.5.2数值方程FLAC3D通过以下三步骤来求解：（1）有限差分逼近（变量的一阶导数、时间导数用有限差分来逼近，假定变量在很短的空间内和时间间隔内线性变化）（2）离散逼近（将连续介质离散为与之相当的网格，在这个网格中，所有的力包括外力和内力，都作用在单元节点的三个方向上）（3）动态求解方法（在平衡方程中引入惯性定律，使得系统慢慢达到平衡）连续介质的运动定律通过以上三步骤，转化为离散单元节点上的牛顿定律。一般普通的微分方程可以通过时间显示差分求解。介质的偏微分平衡方程中涉及到的空间导数，出现在以速度来定义的应变率中。为了进一步的定义速度变量和相关的空间间隔，连续介质被划分为常应变率的四面体单元，这个四面体单元的顶点就是网格单元的顶点。如图2.2所示。图2.2四面体单元（1）空间导数的有限差分形式四面体的应变率张量分量的有限差分方程是通过节点平衡方程得到。四面体共有1-4个节点，节点n所对的面，称之为面n。通过高斯散度定律可得如下方程：,ijijVSvdVvnds（2.45）通过高斯散度定律将四面体上的体积分转化为四个面上的面积分。假设四面体是常应变率，所以速度矢量是线性变化的，每个面的法向方向也是常量，公式2.45简化为：4()()(),1fffijijfVvvnS（2.46）上标（f）表示面f，iv表示变量iv的平均，假设速度线性变化，可以得出：4()1,13fliillfvv（2.47）这里的上标l指的是节点。将公式2.47带入2.46可得：44()(),11,13lffijijlfflVvvnS（2.48）假设在公式2.45中iv=constant，我们得到散度定理如下：4()()10ffjfnS（2.49）利用时2.49可以将2.48化简为：4()(),113lllijijlvvnSV（2.50）因此应变率张量就可以表示为：4()()()116lllllijijjilvnvnSV（2.51）（2）节点平衡运动方程节点运动平衡方程通过虚功原理来推导，在任一个瞬时，得到一个等价静态问题。通过在节点惯性平衡方程中引入虚功原理，来求解整个结构。固定时间t，我们通过平衡方程来研究静力等价状态问题,0ijjiB（2.52）体力的定义如公式2.41iiidvBbdt（2.53）在有限差分的框架中，介质由一些连续的承受体力B的常应变四面体单元代替。在静态平衡方程中，作用在单个四面体上的节点力,1,4nfn，以及四面体应力和体力都是通过虚功方程推导出来。假设四面体有一个虚拟的速度nv（它在四面体内会产生一个线性变化的速度场v和一个常应变率），我们假定外力虚功由节点力nf和体力B产生，内力虚功由虚速度下的应力ij产生。外力功率为：41nniiiiVnEvfvBdV（2.54）内力虚功率为：ijijVIdV（2.55）利用公式2.51，公式2.55可以写成常应变率的形式：4()()()116llllliijjjijilIvnvnS（2.56）由于应力张量是对称张量，现在定义一个矢量lT如下:()()llliijjTnS（2.57）由式2.57，式2.56可写成如下：4113lliilIvT（2.58）将公式2.53带入公式2.54中，可得：41nnbIiinEvfEE（2.59）bE为体力ib对外力虚功率的贡献，IE为惯性力对虚功率的贡献。对于一个四面体常体力：biiVEbvdV（2.60）IiiVdvEvdVdt（2.61）根据以前的假设，在四面体内速度场是线性变化的。为了描述它，现在定义一个新的局部坐标系，坐标系原点在四面体形心，坐标轴为123,,xxx。四面体内任一点速度可以用四面体4个顶点的速度差值得到，如下：41nniinvvN（2.62）,1,4nNn为线性函数0112233nnnnnNccxcxcx（2.63）0123,,,,1,4nnnnccccn,为常数，可以通过求解以下方程得出：123,,njjjnjNxxx（2.64）nj为克罗内克符号，通过定义质心，这项0jVxdV，将公式2.62和2.63带入2.60中得：401bnniinEbvcV（2.65）由式2.64和形心的性质可知014nc，带入式2.65中，得：414bniinbVEv（2.66）将式2.62带入2.61中，得：41InniiVndvEvNdVdt（2.67）最后将式2.66和2.67带入2.59中可得：414nnniiiiVnbVdvEvfNdVdt（2.68）在四面体的静力平衡方程中，在任何虚速度下，外力虚功率等于内力虚功率：由式2.58和2.68可知：34nnniiiiVTbVdvfNdVdt（2.69）假定在四面体内，加速度是不变的，则式2.69中最后一项可改写为：nnniiVVdvdvNdVNdVdtdt（2.70）又因为为常数，0jVxdV，014nc，则式2.7为：4nniiVdvdvVNdVdtdt（2.71）定义4V为虚拟的节点质量nm，节点质量可由式2.72确定，节点质量可以保证在平衡的过程中数值计算的稳定。式2.71改写为：nnniiVdvdvNdVmdtdt（2.72）因此式2.69改写为：34nnnniiiiTbVdvfmdt（2.73）至此等价系统的平衡方程已建立，如式2.73，在每个节点，总的静态平衡力f，包括所有四面体的贡献、节点贡献荷载p以及集中力为0。为了描述以上规律，我们约定如下记号：如果一变量带有l上标，就表示节点变量，l表示该节点在全局节点编号中的编号。.l表示与节点l相连的所有四面体对节点l的贡献之和。采用以上记号，节点牛顿运动定律可以写成如下形式：1,llliindvFMlndt（2.74）nn为介质的节点总数，节点质量lM定义如下：llMm（2.75）不平衡力lF定义如下：34llliiiiTbVFP（2.76）当系统达到平衡时，不平衡力渐渐趋于0。（3）时间导数的显示有限差分考虑本构方程2.43和变形率和节点速度之间的关系2.51，公式2.74可以写成如下形式：1231,,,,,,1,lllpiiiiiinldvFtvvvvlndtM（2.77）记号{}l表示计算节点l速度的速度子集。在FLAC3D中求解2.77一般采用对时间t的显示有限差分。在这个方法中假定节点速度在时间间隔t内线性变化。公式2.77中左边的导数采用中心差分来代替。中心差分时，当计算位移、力时只存取一半时间步的速度。节点速度通过以下迭代公式计算：123()(),,,,,,22llllpiiiiiiiltttvtvtFtvvvvM（2.78）同理节点坐标更新也是通过中心差分的形式得到：()()()2llliiitxttxttvt（2.79）从式2.78和2.79中可知，一阶误差项不存在，因此上述迭代方法是二阶精度的节点位移迭代公式如下：()()()2llliiituttuttvt（2.80）(0)0liu（4）本构方程的增量形式在FLAC3D里面假定在时间间隔t，速度为常量，因此本构方程2.43可以写成如下形式：*ˆ,ijijijijHt（2.81）ˆij为共轭应力增量，*ijH为一已知函数。假定位移在时间t内，线性变化：ijijt（2.82）ij为时间t时的应变改变量。应力增量通过ˆij计算得到，如下：ˆCijijij（2.83）Cij为一应力修正，定义如下：Cijikkjikkjt（2.84）转动率张量由公式2.40或是2.51计算（5）大小应变模式以上描述的微分方程都是描述大变形，涉及到大位移，以及位移的线性变化和转动。在FLAC3D里面称之为大变形模式。在实际应用中，转动足够小，以至于分量ijij与整体相比很小，可以忽略不计。张量可以用I来代替，公式2.83中的应力修正项可以忽略。当然在小位移和位移变化中，应变率张量中涉及的空间导数，可以通过初