计算力学第13课 差分原理

E-mail:，heu.ComputationalMechanics

(13)

本次课的主要问题：差分原理应力函数的差分解基本内容差分原理

设有x的解析函数y=f(x)，从微分学知道函数y对x的导数为式中

dy/dx是函数对自变量的导数，又称微商Dy、Dx分别称为函数及自变量的差分Dy/Dx为函数对自变量的差商。差分原理常微分方程数值解工程上很多数学表述都可以归结为常微分的定解问题，很多偏微分方程，也可以化为常微分方程问题来近似求解。对于一个常微分方程通常会有无穷个解，如因此，要加上一个定解条件，如差分原理常微分方程的解是一个函数，但是，计算机没有办法对函数进行运算，因此，常微分方程的数值解并不是求函数的近似值，而是求解函数在某些节点上的近似值。这种方法，称为数值离散方法。求的是未知函数在一系列离散点上的值的近似。得到高精度方法的一个直接想法是利用Taylor展开。假设式中的f(x,y)充分光滑，将y(xi+1)在xi点作Taylor展开，若取右端不同的有限项作为y(xi+1)的近似值，就可得到计算y(xi+1)的各种不同截断误差的数值公式。差分原理例如：取前两项可得到截断误差为O[(x)2]的公式若取前三项，可得到截断误差为O[(x)3]的公式类似地，若取前p+1项作为y(xi+1)的近似值，便得到差分原理欧拉(Euler)法寻求通过某种方法去获得解y(x)在点xi上的近似值yi，即为实现这一目标，Euler方法首先将微分算子离散化即推广这就是显式的Euler公式，又称差分格式。显然与y(xi+1)的Taylor展开式的前两项完全相同，即局部截断误差为O[(x)2]。Euler公式可改写成差分原理改进的欧拉法对方程y=f(x,y)，两边由xi到xi+1积分，并利用梯形公式，有：故有公式：此即改进的欧拉法。差分原理同理，改进Euler公式可改写成局部截断误差为O[(x)3]。上述两组公式在形式上共同点：都是用f(x,y)在某些点上值的线性组合得出y(xi+1)的近似值yi+1，且增加计算的次数，可提高截断误差的阶。欧拉法：每步计算一次f(x,y)的值，为一阶方法。改进欧拉法：需计算两次f(x,y)的值，为二阶方法。差分原理龙格-库塔(Runge-Kutta)法于是可考虑用函数f(x,y)在若干点上的函数值的线性组合来构造近似公式，构造时要求近似公式在(xi,yi)处的Taylor展开式与解y(x)在xi处的Taylor展开式的前面几项重合，从而使近似公式达到所需要的阶数。既避免求高阶导数，又提高了计算方法精度的阶数。或者说，在[xi,xi+1]这一步内多计算几个点的斜率值，然后将其进行加权平均作为平均斜率，则可构造出更高精度的计算格式。这就是龙格-库塔(Runge-Kutta)法的基本思想。差分原理一般龙格-库塔方法的形式为称为p阶龙格-库塔方法。其中ai，bij，ci为待定参数，要求上式yi+1在点(xi，yi)处作Taylor展开，通过相同项的系数确定参数。一阶龙格-库塔公式，也就是欧拉法二阶龙格-库塔公式，也就是改进的欧拉法差分原理

差分方法是一类重要的数值解法。

寻求一系列离散节点x1<x2<…<xn<…上的近似解y1，y2，…，yn，…。

h=xn+1-xn称为步长。

初值问题差分方法的特点：

步进式——求解过程顺着节点排列的次序一步一步地向前推进。

描述这种算法，只要给出从已知信息yn，yn-1，yn-2

，…计算yn+1的递推公式——差分格式。

求解的核心——消除导数，离散化方法。差分原理有限差分法(finitedifferencemethod)微分方程和积分方程数值解的方法。基本思想把弹性力学的基本方程和边界条件(一般均为微分方程)近似地改用差分方程(代数方程)来表示，把求解微分方程的问题改换成为求解代数方程的问题，属于数学上的近似。(有限单元法是以有限个单元的集合体来代替连续体，属于物理上的近似)。有限差分法的主要内容包括：根据问题的特点将定解区域作网格剖分把原微分方程离散化为差分方程组解差分方程组此外为了保证计算过程的可行和计算结果的正确，还需从理论上分析差分方程组的性态，包括解的唯一性、存在性和差分格式的相容性、收敛性和稳定性有限差分法简介构造差分的方法有多种形式，目前主要采用的是泰勒级数展开方法。其基本的差分表达式主要有三种形式：一阶向前差分一阶向后差分一阶中心差分和二阶中心差分其中前两种格式为一阶计算精度，后两种格式为二阶计算精度。通过对时间和空间这几种不同差分格式的组合，可以组合成不同的差分计算格式。向前差分向后差分中心差分上面谈的是一阶导数，对应的称为一阶差分。有限差分法简介对一阶差分再作一阶差分，得到二阶差分，记为D2y。以向前差分为例，有依此类推，任何阶差分都可由其低一阶的差分再作一阶差分得到。例如n阶前差分为差分原理以一维的情况为例向前差分一阶差分二阶差分向后差分一阶差分二阶差分中间差分一阶差分二阶差分差分原理yi-2yi-1yiyi+1yi+2yi-2yi-1yiyi+1yi+2yi-2yi-1yiyi+1yi+2yi-2yi-1yiyi+1yi+2yi-2yi-1yiyi+1yi+2yi-2yi-1yiyi+1yi+2yi-1/2yi+1/2函数的差分与自变量的差分之比，即为函数对自变量的差商。

一阶向前差商为一阶向后差商为一阶中心差商为或差分原理yi-2yi-1yiyi+1yi+2yi-1/2yi+1/2二阶差商多取中心式，即

对于多元函数f(x,y,…)的差分与差商也可以类推。如一阶向前差商为差分原理截断误差差商与导数之间的误差表明差商逼近导数的程度，称为截断误差。由函数的Taylor展开，可以得到截断误差相对于自变量差分(增量)的量级，称为用差商代替导数的精度，简称为差商的精度。一阶向后差商也具有一阶精度。差分原理将f(x+Dx)和f(x-Dx)的Taylor展开式相减可得可见一阶中心差商具有二阶精度。可见二阶中心差商的精度也为二阶精度。差分原理非均匀步长在有些情况下要求自变量的增量本身是变化的，相应的差分和差商就是不等距的。一阶向后差商一阶中心差商差分原理Ox弹性体上，用相隔等间距h而平行于坐标轴的两组平行线织成正方形网格，Dx=Dy=h，如图。设f=f(x，y)为弹性体内的某一个连续函数。该函数在平行于x轴的一根网线上，如在3-0-1上，它只随x坐标的改变而变化。在邻近节点0处，函数f可展为泰勒级数如下：应力函数的差分解只考虑离开节点0充分近的那些节点，即(x-x0)充分小。于是可不计(x-x0)的三次及更高次幂的各项，则上式简写为：在节点3，x=x0-h，在节点1，

x=x0+h，代入上式得：联立上两式，解得差分公式应力函数的差分解同理，在网线4-0-2上可得到差分公式差分公式中分母为2h的是以相隔2h的两节点处的函数值来表示中间节点处的一阶导数值，可称为中点导数公式。以相邻三节点处的函数值来表示一个端点处的一阶导数值，可称为端点导数公式。中点导数公式与端点导数公式相比，精度较高。因为前者反映了节点两边的函数变化，而后者却只反映了节点一边的函数变化。因此，尽可能应用前者，而只有在无法应用前者时才不得不应用后者。应力函数的差分解以上4式是基本差分公式，从而可导出其它的差分公式如下：应力函数的差分解当不计体力时，将弹性力学平面问题归结为在给定边界条件下求解双调和方程的问题。用差分法解平面问题，就应先将双调和方程变换为差分方程，而后求解。一旦求得弹性体全部节点的j值后，就可按应力分量差分公式(对节点0)算得弹性体各节点的应力。应力函数的差分解可见，用差分法解平面问题，共有两大任务：

建立差分方程

将(D.1)~(D.3)代入双调和方程对于弹性体边界以内的每一节点，都可以建立这样一个差分方程。应力函数的差分解

联立求解这些线性代数方程，就能求得各内节点处的值

一般建立和求解差分方程，在数学上不会遇到很大困难。但是，当对于边界内一行的(距边界为h的)节点，建立的差分方程还将涉及边界上各节点处的j值，并包含边界外一行的虚节点处的j值。为了求得边界上各节点处的j值，须要应用应力边界条件，即：应力函数的差分解由右图可见于是，上式可改写为由此得应力函数的差分解关于边界上任一点处的j/y，j/x值，可将上式从A点到B点对s积分得到：应力函数的差分解由高等数学可知将此式也从A点到B点沿s进行积分，就得到边界上任一点B处的j值。为此利用分部积分法，得：整理得：可见，已知即可根据面力分量及导数求得应力函数的差分解由分析可知，把应力函数加上一个线性函数，并不影响应力。因此，可设想把应力函数加上a+bx+cy，然后调整a，b，c三个数值，使得这样应力函数的差分解从图看出，前面3式右边的积分式分别表示A与B之间的x方向的面力之和、A与B之间的y方向的面力之和、A与B之间的面力对于B点的矩。至此，解决了计算边界上各节点的值的问题。应力函数的差分解

边界外一行虚节点处j的值，则可用边界上节点处的j/y或j/x值和边界内一行相应节点处j的值来表示。例如，对于图中的虚节点14，因为有所以有当求出全部节点上的j值以后，我们就可按应力分量的差分公式计算应力分量。应力函数的差分解用差分法解弹性平面问题时，可按下列步骤进行：1.在边界上任意选定一个节点作为基点A，取然后由面力的矩及面力之和算出边界上所有各节点处j的值，以及所必需的一些j/y及j/x值，即垂直于边界方向的导数值。2.应用公式(D.6)，将边界外一行虚节点处的j值用边界内的相应节点处的j值来表示。3.对边界内的各节点建立差分方程(D.5)，联立求解这些节点处的值。应力函数的差分解4.按照公式(D.6)，算出边界外一行的各虚节点处的j值。5.按照公式(D.4)计算应力的分量。应力函数的差分解说明以上是针对单连体导出的结果。如果一部分边界是曲线的，或是不与坐标轴正交，则边界附近将出现不规则的内节点。对于这样的节点，差分方程(D.4)必须加以修正。应力函数的差分解差分方程的相容性、收敛性、稳定性相容性对于一个微分方程建立的各种差分格式，为了有实用意义，一个基本要求是它们能够任意逼近微分方程。导数与其差分近似式之间存在截断误差。因此，差分方程的解并不是严格的，而是近似地满足原来的偏微分方程。但是，当时间步长Dt和空间步长Dx都趋近于零时，差分方程的截断误差也趋近于零，差分方程的极限形式就是原偏微分方程。这时，认为差分方程与偏微分方程是相容的，这种相容性表示差分方程“收敛”于原偏微分方程。是收敛性的必要条件。有限差分法收敛性一个差分格式是否有用，最终要看差分方程的精确解能否任意逼近微分方程的解。差分方程的解，即当步长Dt→0

、Dx→0时收敛于原偏微分方程的解。有限差分法稳定性差分格式的计算过程是逐步推进的，在计算第n+1步的近似值时要用到第n步的近似值，直到与初始值有关。

前面各步若有舍入误差，必然影响到后面各步的值，如果误差的影响越来越大，以致差分格式的精确解的面貌完全被掩盖，这种格式是不稳定的。

相反如果误差的传播是可以控制的，或者说若这种误差积累保持有界，就认为格式是稳定的，只有在这种情形，差分格式在实际计算中的近似解才可能任意逼近差分方程的精确解。有限差分法优点差分法是解微分方程边值问题和弹性力学问题的有效方法；差分法简便易行；对于某些结构，为了更精确地分析局部的应力状态，可以用差分法进行分析；缺点对于曲线边界和斜边界等产生的不等间距网格的处理，比较麻烦和易于出错；比较适用于求解二维问题或平面问题；比较适用于等间距网格，对于应力变化较为剧烈时，需要采用二次网格进行计算。有限差分法有限差分法项目有限差分法有限元法试函数局部近似局部近似程序难易程度很好好程序灵活性好很好精确性差好计算效率好好适宜的方程各类型椭圆型主要优点经济、程序简单灵活性好主要缺点较难推广到高阶不经济以位势问题为例，例如对于各向同性二维稳态导热方程，如果物体内部热源密度为0，方程可写成如下形式：如方程和