第四章-有限差分法及软件

第四章有限差分法及软件第四章第1页本章主要内容一有限差分法的理论基础二平面问题有限差分数值原理与方法三三维问题有限差分数值原理与方法四FLAC软件介绍五有限元法与有限差分的比较第四章第2页一、有限差分的理论基础xyh0123456789101112h弹性力学中的差分法是建立有限差分方程的理论基础在弹性体上用相隔等间距h而平行于坐标轴的两组平行线划分成网格。),(yxff弹性体内某一个连续函数，它可能是某一个应力分量或位移分量，也可能是应力函数、温度、渗流等等。图4-1有限差分网格第四章第3页这个函数，在平行于x轴的一根格线上，例如在3-0-1上xyh0123456789101112h它只随x坐标的变化而改变。在邻近结点0处，函数f可以展开为泰勒级数：400443003320022000)(!41)(!31)(!21)(xxxfxxxfxxxfxxxfff（4-1）第四章第4页在结点3及结点1，x分别等于x0-h及x0+h，即：x-x0分别等于-h和h。将其代入上式，有：xyh0123456789101112h0444033302220032462xfhxfhxfhxfhff（4-2）0444033302220012462xfhxfhxfhxfhff（4-3）第四章第5页假定h是充分小的，因而可以不计它的三次幂及更高次幂的各项，则(4-2)式及(4-3)式简化为：02220032xfhxfhff02220012xfhxfhff（4-4）（4-5）第四章第6页联立求解(4-4)式及(4-5)式，得到差分公式：hffxf231020310222hfffxf（4-6）（4-7）同理可以求得y方向差分公式：hffyf242020420222hfffyf（4-8）（4-9）第四章第7页公式(4-6)至(4-9)是基本差分公式，通过这些公式可以推导出其它的差分公式。例如，利用(4-6)和(4-8)式，可以导出混合二阶导数的差分公式：075862002)()(41ffffhyfxyxf（4-10）用同样的方法，由公式(4-7)及(4-9)可以导出四阶导数的差分公式。有限差分法不仅仅局限矩形网格，Wilkins(1964)提出了推导任何形状单元的有限差分方程的方法。与有限元法类似，有限差分方法单元边界可以是任何形状、任何单元可以具有不同的性质和值的大小。第四章第8页本章主要内容一有限差分法的理论基础二平面问题有限差分数值原理与方法三三维问题有限差分数值原理与方法四FLAC软件介绍五有限元法与有限差分的比较第四章第9页1.平面问题有限差分基本原理二平面问题有限差分数值原理与方法对于平面问题，将具体的计算对象用四边形单元划分成有限差分网格，每个单元可以再划成两个常应变三角形单元。三角形单元的有限差分公式用高斯发散量定理的广义形式推导得出：dAxffdsnAisi（4-11）s——绕闭合面积边界积分；ni——对应表面s的单位法向量；f——标量、矢量或张量；xi——位置矢量；ds——微量弧长；A——对整个面积A积分。abababui(b)ui(a)ΔSFini(1)ni(2)s(1)s(2)图4-2四边形单元划分成两个常应变三角形单元第四章第10页在面积A上，定义f的梯度平均值为：AiidAxfAxf1（4-12）将(4-11)式代入上式，得：SiifdsnAxf1（4-13）对一个三角形子单元，(4-13)式的有限差分形式为：SiisnfAxf1（4-14）式中，Δs是三角形的边长，求和是对该三角形的三个边进行。f的值取该边的平均值。第四章第11页平面问题有限差分法基于物体运动与平衡的基本规律。最简单的例子是物体质量为m、加速度为dū/dt与施加力F的关系，这种关系随时间而变化。牛顿定律描述的运动方程为：Fdtudm.（4-15）当几个力同时作用于该物体时，如果加速度趋于零，即：F=0(对所有作用力求和)，(4-15)式也表示该系统处于静力平衡状态。对于连续固体，(4-15)式可写成如下广义形式：ijijgxtu.（4-16）——物体的质量密度；t——时间；xi——坐标矢量分量；gi——重力加速度(体力)分量；ij——应力张量分量。该式中，下标i表示笛卡尔坐标系中的分量，复标喻为求和。第四章第12页利用(4-14)式，将f替换成单元每边平均速度矢量，这样，单元的应变速率ēij可以用结点速度的形式表述：snuuAxujSbiaiji)(21)(.)(..（4-17）][21...ijjiijxuxue（4-18）式中：(a)和(b)是三角形边界上两个连续的结点。注意：如果结点间的速度按线性变化，(4-17)式平均值与精确积分是一致的。通过(4-17)式和(4-18)式，可以求出应变张量的所有分量。第四章第13页根据力学本构定律，可以由应变速率张量获得新的应力张量：),,(:.keMijijij（4-19）式中，M()——表示本构定律的函数形式；k——历史参数，取决于特殊本构关系；:=——表示“由…替换”。通常，非线性本构定律以增量形式出现，因为在应力和应变之间没有单一的对应关系。当已知单元旧的应力张量和应变速率(应变增量)时，可以通过(4-19)式确定新的应力张量。例如，各向同性线弹性材料本构定律为：teGeGKijkkijijij}2)32({:..（4-20）第四章第14页式中，δij——Kronecker记号；t——时间步；G，K——分别是剪切模量和体积模量。在一个时步内，单元的有限转动对单元应力张量有一定的影响。对于固定参照系，此转动使应力分量有如下变化：tkjikkjikijij)(:（4-21）}{21..ijjiijxuxu式中，（4-22）第四章第15页在大变形计算过程中，先通过(4-21)式进行应力校正，然后利用(4-20)式(或本构定律(4-19)式)计算当前时步的应力。计算出单元应力后，可以确定作用到每个结点上的等价力。在每个三角形子单元中的应力如同在三角形边上的作用力，每个作用力等价于作用在相应边端点上的两个相等的力。每个角点受到两个力的作用，分别来自各相邻的边(图4-2)。因此：)(21)2()2()1()1(SnSnFjjiji（4-23）由于每个四边形单元有两组两个三角形，在每组中对每个角点处相遇的三角形结点力求和，然后将来自这两组的力进行平均，得到作用在该四边形结点上的力。在每个结点处，对所有围绕该结点四边形的力求和Fi，得到作用于该结点的纯粹结点力矢量。该矢量包括所有施加的载荷作用以及重力引起的体力Fi(g)。abababui(b)ui(a)ΔSFini(1)ni(2)s(1)s(2)第四章第16页gigimgF)(（4-24）式中，mg是聚在结点处的重力质量，定义为联结该结点的所有三角形质量和的三分之一。如果四边形区域不存在（如空单元），则忽略对Fi的作用；如果物体处于平衡状态，或处于稳定的流动（如塑性流动）状态，在该结点处的Fi将视为零。否则，根据牛顿第二定律的有限差分形式，该结点将被加速：mtFuutittitti)()2/(.)(.（4-25）第四章第17页式中，上标表示确定相应变量的时刻。对大变形问题，将(4-25)式再次积分，可确定出新的结点坐标：tuxxttititti)2/(.)(.)(.（4-26）注意到(4-25)式和(4-26)式都是在时段中间，所以对中间差分公式的一阶误差项消失。速度产生的时刻，与结点位移和结点力在时间上错开半个时步。第四章第18页2.显式有限差分算法——时间递步法期望对问题能找出一个静态解，然而在有限差分公式中包含有运动的动力方程。这样，可以保证在被模拟的物理系统本身是非稳定的情况下，有限差分数值计算仍有稳定解。对于非线性材料，物理不稳定的可能性总是存在的，例如：顶板岩层的断裂、煤柱的突然垮塌等。在现实中，系统的某些应变能转变为动能，并从力源向周围扩散。有限差分方法可以直接模拟这个过程，因为惯性项包括在其中——动能产生与耗散。相反，不含有惯性项的算法必须采取某些数值手段来处理物理不稳定。尽管这种做法可有效防止数值解的不稳定，但所取的“路径”可能并不真实。第四章第19页图4-3是显式有限差分计算流程图。计算过程首先调用运动方程，由初始应力和边界力计算出新的速度和位移。然后，由速度计算出应变率，进而获得新的应力或力。每个循环为一个时步，图4-3中的每个图框是通过那些固定的已知值，对所有单元和结点变量进行计算更新。平衡方程（运动方程）应力/应变关系（本构方程）新的应力或力新的速度和位移图4-3有限差分计算流程图第四章第20页例如，从已计算出的一组速度，计算出每个单元的新的应力。该组速度被假设为“冻结”在框图中，即：新计算出的应力不影响这些速度。这样做似乎不尽合理，因为如果应力发生某些变化，将对相邻单元产生影响并使它们的速度发生改变。然而，如果我们选取的时步非常小，乃至在此时步间隔内实际信息不能从一个单元传递到另一个单元(事实上，所有材料都有传播信息的某种最大速度)。因为每个循环只占一个时步，对“冻结”速度的假设得到验证——相邻单元在计算过程中的确互不影响。当然，经过几个循环后，扰动可能传播到若干单元，正如现实中产生的传播一样。显式算法的核心概念是计算“波速”总是超前于实际波速。所以，在计算过程中的方程总是处在已知值为固定的状态。这样，尽管本构关系具有高度非线性，显式有限差分数值法从单元应变计算应力过程中无需迭代过程，这比通常用于有限元程序中的隐式算法有着明显的优越性，因为隐式有限元在一个解算步中，单元的变量信息彼此沟通,在获得相对平衡状态前，需要若干迭代循环。显式算法的缺点是时步很小，这就意味着要有大量的时步。因此，对于病态系统——高度非线性问题、大变形、物理不稳定等，显式算法是最好的。而在模拟线性、小变形问题时，效率不高。第四章第21页由于显式有限差分法无需形成总体刚度矩阵，可在每个时步通过更新结点坐标的方式，将位移增量加到结点坐标上，以材料网格的移动和变形模拟大变形。这种处理方式称之“拉格朗日算法”，即：在每步计算过程中，本构方程仍是小变形理论模式，但在经过许多步计算后，网格移动和变形结果等价于大变形模式。用运动方程求解静力问题，还必须采取机械衰减方法来获得非惯性静态或准静态解，通常采用动力松弛法，在概念上等价于在每个结点上联结一个固定的“粘性活塞”，施加的衰减力大小与结点速度成正比。前已述及，显式算法的稳定是有条件的：“计算波速”必须大于变量信息传播的最大速度。因此，时步的选取必须小于某个临界时步。若用单元尺寸为Δx的网格划分弹性体，满足稳定解算条件的时步Δt为：Cxt（4-27）第四章第22页式中，C是波传播的最大速度，典型的是P波Cp：3/4GKCp（4-28）对于单个质量—弹簧单元，稳定解的条件是：kmt2（4-29）式中，m是质量，k是弹簧刚度。在一般系统中，包含有各种材料和质量—弹簧联结成的任意网络，临界时步与系统的最小自然周期Tmin有关：minTt（4-30）第四章第23页下面，通过一个简单例子，说明显式有限差分法在解题过程