FLAC3D理论基础_屏幕版

FLAC3D理论基础这部分阐述的是FLAC3D的有关理论。FLAC3D很大一部分是二维FLAC的扩展，而显式有限差分法是FLAC和FLAC3D的共同的理论基础，有关这一部分，可参考FLAC用户手册。尽管如此，二维和三维的方程还是有一些明显的不同，特别是在数学模型的扩展上。这里主要讨论三维模型在FLAC3D中的实现方法。1.三维显示差分模型模型的构成FLAC3D是显式有限差分程序，可以模拟连续三维介质达到平衡状态或稳定塑性流动时的力学行为。这种力学行为，可以通过建立特定的数学模型和特定的数字模拟方法来实现。下面就来阐述这两方面的有关内容。1.1数学模型介质的力学特征可通过一般的力学关系（如应变的定义、运动方程等）和理想介质的本构方程进行推导。所得到的数学表达式是一系列的偏微分方程及相关变量如：静力学中应力和动力学中的应变速率、速度等。对于特定的具有几何特征和特殊性质的介质，这些方程和变量在给定的边界条件和初始条件下，可以求解。尽管FLAC3D主要是研究处于极限平衡状态下的介质变形及应力状态，但它的模型里可以包含有运动方程是它的一大特色。在进行数字模拟过程中，由于惯性物体将达到稳定状态或平衡状态。1.1.1符号约定在FLAC3D的拉格朗日公式中，用矢量（iiivux，，），dtdvi（其中i=1，3）来分别表示介质中点的空间位置、位移、速度和加速度。作为一种符号约定，据上下文的不同，斜体字可以矢量和张量。如：符号ia表示笛卡儿坐标系下矢量][a的i分量；Aij表示张量[A]的（i,j）分量。还有，i,表示对xi的偏导数（其中可以是标量，也可以矢量或张量的分量。规定：拉力和张力为正。爱因斯坦的求和约定只适用于i,j,k（i,j,k=1，2，3）1.1.2应力给定点的应力状态可用一个对称的应力张量ji,来表示。由柯西定理，若一个面的单位法矢量为[n]，则它的拖曳矢量[t]：jijint（1）1.1.3应变速率与转动速率假定介质颗粒以速度[v]运动，则在无穷小的时间内，发生无穷小应变dtvi，相应的应变张量可写为：ijjiijvv,,21（2）式中是对空间位置矢量的偏导数。在下面的论述中，第一应变速率变量是表征单元体积的膨胀率的量。介质的形变率除了张量ij外，还有刚体位移[v]以及转动速率：jkijkiwe21（3）式中：eijk为符号函数，[w]为转动速率张量，其分量定义为：ijjiijvvw,,21（4）1.1.4运动方程与平衡方程由运动定律的连续形式的柯西运动方程：dtdvbiijij,（5）式中：ρ介质单位体积的质量，[b]作用于单位质量上的体力，d[v]/dt是速度对时间的偏导数。这些定理控制了介质单元体在力的作用下运动状态。当介质处于静力平衡时，加速度d[v]/dt为0，由上式得平衡状态下的偏微分方程：0,ijijb（6）1.1.5边界条件及初始条件边界条件包括施加于边界上的牵引力和（或）速度（由位移而产生的）。另外，还应考虑体力和初始应力状态。1.1.5本构方程运动方程式（5）和应变速率的定义式（2）共有9个方程和15个未知量，这15未知量包括6个应力速率分量、6个应变速率分量及3个速度矢量分量。考虑特定的材料性质，有另外6个相联系的本构方程，通常用下式表示：kHijijijij,,][（7）式中：ij][为共转（co-rotational）应力速率张量，[H]为一给定的函数，k是与加载路径有关的参数。共转（co-rotational）应力速率张量，等于给定参考系的介质内一点的应力的偏导数和以瞬时角速度的转动，数学表达式如下：kjijijijijijwwdtd][（8）式中，dtd][是应力][对时间的偏导数，[w]为转动速率张量。1.2数学模拟FLAC3D有以下三种逼近方法：（1）有限差分法。假定在变量在空间和时间内线性变化，用变量对空间和时间的一阶导数来近似等于它的有限差分值。（2）离散单元法。将连续介质离散为等效块体集合体，所有的力（施加的作用力与相互作用力）作用在三维网格的节点上。（3）动力学解法。运用运动方程求解所研究系统达到平衡状态时的参量。利用以上逼近方法，连续介质的运动定律可变为节点上的牛顿定律的离散形式。从而可通过显式有限差分法来求解一般的差分方程。等价介质空间偏导数在由速度定义的应变速率用到。因此，为了定义速度变量和相应的空间间距，介质需离散为常应变速率四面体单元，它的顶点即为网格的节点（图1）。图1四面体单元的面和节点1.2.1空间微分与有限差分的近似作为基础，下面由节点运动方程推导四面体的应变速率张量各分量的有限差分公式。四面体的节点1到4，面n指与节点n相对的面（见图1）。由四面体高斯离散理论，可知：dSnvdVvSjiVji,（9）上式，左边和右边分别为对四面体体积和面积积分，[n]为面的单位外法矢量。对常应变速率四面体，速度场是线性变化的，而每一个面的外法矢量[n]也为一常值。因此，式（9）积分后有：)()()(41,ffjffijiSnvVv（10）式中，上标（f）指面f的相关变量值，iv指i速度分量的均值。若速度呈线性变化，则：4,1)()(31flllifivv（11）上标l指节点l的值。将上式代入式（10），有：)(4,1)(41,31flfffjllijiSnvVv（12）在式（9）中，若vi=1，则由离散理论可得：0)(41)(fffjSn（13）所以，式（12）两边同除以V，则有：)(41)(,31llljlijiSnvVv（14）而应变速率张量则可由下式表示：)(41)()(,)(61llliljljlijiSnvnvV（15）1.2.2运动定律的节点方程对静力学问题，一般用虚功原理来推导运动定律的节点方程。节点习惯术语也与解相应的平衡方程的所用的相类似（式（6））。一定的时间t内，我们认为静力平衡问题一般可通过以的平衡方程进行求解：0,ijijB（16）式中，体力的定义见式（5），其中：dtdvbBiii（17）在这里所用的有限差分逼近格式中，假设介质为受体力[B]的常应变四面体的连续集合体。根据虚功原理，作用于单个四面体上的节点力[f]n（n=（1,4））与所受到的四面体应力和等效体力达到一种静态平衡。若节点的虚速度][v（它在四面体中产生线性速度场][v和常应变速率][），则由此产生的节点力和体力等外力功功率等于内部应力产生的内力功功率。根据1.1.1节的的符号约定（上标表示变量的在该节点上的值）和爱因斯坦的求和约定，外力功功率可表示为：dVBvfvEinVinini41（18）而内力功功率：dVIijVij（19）由式（15），对常应变速率的四面体有：)(4161llliijljljijliSnvnvI（20）应力张量是对称张量，定义矢量Tl：)(lljijliSnT（21）则：4131lliliTvI（22）式（17）代入式（18），有：IbninniEEfvE41（23）Eb和EI分别体力ib和惯性力所作的外力功功率。若四面体内体力ib为常数，则有：ViibdVvbE（24）EI为：ViiIdVdtdvvE（25）由前述有限差分逼近假定，四面体单元内的速度场符合线性变化。因此，我们引进一个参考坐标系（它的坐标原点则四面体的中心上），则有：41nnniiNvv（26）式中Nn（n=1，4）为一线性函数：3322110xcxcxccNnnnnn（27）其中，nnnncccc3210,,,（n=1，4）为下述方程的解：njjjjnxxxN321,,（28）式中，nj是克罗内克尔增量（Kroneckerdelta）。根据重心的定义，所有形如VjdVx'的积分为0，由式（27），再将式（26）代入式（23）得：410nnniibVcvbE（29）克莱姆（Cramer）解式（28），再由重心的性质，有：410nc（30）将上式代入式（29），有：414niVnibbvE（31）同理，将式（26）代入式（25）得到：41nVinniIdVdtdvNvE（32）将式（31）和（32）代入式（23）：414nVinininidVdtdvNVbfvE（33）在等效问题中，处于静态平衡的四面体，在任意的虚速度下，内力功功率（式（22））等于外力功功率所以，由上述方程式，有：VinininidVdtdvNVbTf43（34）若四面体内加速度在均值附近只有微小的空间变化，则有：VnniVindVNdtdvdVdtdvN（35）四面体内，为一常数，由重心的性质（式（27）和（30）），则上式可写为：niVindtdvVdVdtdvN4（36）考虑平衡方程数学上的对称性，用虚拟节点质量mn代替上式中的质量4V：ninVindtdvmdVdtdvN（37）则，式（34）可写为：ninininidtdvmVbTf43（38）由等效系统的平衡条件，各个节点上，等效静力-[f]与有何载荷集中力引起的作用在四面体内和节点上的外力[P]的和应等于0。为了说明这种情况，我们引入一些符号，上标l表示在所有节点中标号为l的节点上的变量，[[.]]l表示作用在所有与包含节点l的四面体上的变量的和。根据上述约定，可将节点牛顿定律表达为：nlillinldtdvMF,1（39）式中，nn介质中的所有的节点总数，节点质量定义为：llmM（40）不平衡力[F]l定义为：liliiliPVbTF43（41）当介质达到平衡时，不平衡力等于0。1.2.3显式有限差分逼近的时间导数考虑本构方程（式（7））和变形速率与节点速率之间的关系（式（15）），式（35）可表示为一般的差分方程：nlpiiiilillinlkvvvvtFMdtdv,1,,....,,,,1321（42）式中，{}l是指节点l速度的一个子集（式（39））。在FLAC3D中，是利用显式差分方程进行求解的。在这种逼近方法中，假定节点速度随时间间隔t线性变化，式（42）左边的速度对时间的偏导数用中心差分格式，因此，随与位移和力有关的半时间步长进行迭代。节点的速度利用下式进行计算：kvvvvtFMttvttvlpiiiilillili,,....,,,,1)2()2(321（43）类似地，节点的位置也用中心差分进行迭代：)2()()(tttvtxttxlilili（44）由上述可以看出，中心差分可避免一阶误差，具有二阶精确性。因此，节点位移也有如下关系：)2()()(tttvtuttulilili（45）其中，0)0(liu。1.2.4本构方程的增量形式FLAC3D中，时间t内为常数。本构方程（式（7））增量表示为：),(*tHijijijij（46）在t内，若位移和位移梯度很小，在看认为：ttijij（47）式中ij为t时刻时的结构的应变增量值。应力增量ij可通过ij进行计算：Cijijij（48）式中，Cij为修正应力，定义如下：twwkjkjkjijCij)(（49）由转动张量分量用（式（4））和有限差分方程（式（14））得：)(41)()(61llliljljliijSnvnvVw（50）1.2.5大应变和小应变模式上述的有关大应变的变形模式的数学