3第二章-有限差分方法基础

TheElementsofComputationalFluidDynamics第二章有限差分方法基础§2.1有限差分方法概述§2.2导数的数值逼近方法§2.3差分格式的性质§2.4发展方程的稳定性分析§2.1有限差分方法概述以一维非定常热传导方程为例，介绍有限差分方法的概念、简单构造方法和求解过程。2.1.1基本方程和定解问题22(0)(2.1.1)uutx(,)[0,1][0,]xt求解域：(,0)()(2.1.2)(0,)(),(1,)()uxfxutatutbt初始条件：边界条件：方程(2.1.1)和初边条件(2.1.2)构成了一个适定的定解问题。有限差分方法：对于一个偏微分方程，如果把方程中的所有偏导数近似地用代数差商(AlgebraicDifferenceQuotient)代替，则可以用一组代数方程近似地替代这个偏微分方程，进而得到数值解，这种方法称为有限差分方法(FiniteDifferenceMethod)。012.1.2求解域及偏导数的离散化为了用有限差分方法求解式(2.1.1)，需要把其中的偏导数表示为代数形式，为此，首先要把自变量从连续的分布变为离散形式。这个过程称为求解域的离散化。1.空间求解域的离散化0x1x2x1MxMx把空间求解域分为M段（均匀剖分）0121,,,,,MMxxxxx网格点：=kxkx显然,=1/xM其中,，为空间步长。2.时间变量的离散化把感兴趣的时间段(t=T之前)分为N段（均匀剖分），则时间方向的求解域可以划分为01211,,,,,NNNttttt个离散时刻：=(/,)ntnttTN时间步长求解域被划分为一系列离散的时空网格点图2.1求解域的离散化3.解的离散表示目标：求出所有网格点上物理量u的近似解。(,)=(,)(0,1,,;0,1,,)knuxtukxntkMnN(,)nknkuxtu后文中，把记为。4.导数的数值逼近把方程中的偏导数项近似表示为代数形式。(,)(2.1.1)(,)(,)(2.1.3)()()(2.1.4)kntxxnntkxxkxtukxntukxntuu在网格点，方程可以表示为或000Taylor(,)(,)(,)lim(,)(,)(,)lim(,)(,)(,)lim2tttttttuuxttuxtuxttuxtuxttuxttuxttuxttuxtt按定义（利用展开式），偏导数可以写成下面几种等价形式：0limdifferencequotienttt其中，后面的项称为差商()。当足够小时，可以用差商来近似导数。(,)(,)(,)(,)(,)(,)(,)(,)(,)2tttuxttuxtuxttuxtuxttuxttuxttuxttuxtt即：1-1+1-1(,)()()(2.1.5)()()(2.1.6)()()(2.1.7)2knnnnnkktktknnnnkktktknnnnkktktkxtuuuutuuuutuuuut在网格点，有的向前差商：的向后差商：的中心差商：1111(forwarddifference)(backwarddifference)(centraldifference)nnntkkknnntkkknnntkkkuuuuuuuuuu沿时间方向的向前差分：向后差分：中心差分：ttt其中，，，分别称为时间方向前差、后差和中心差分算子。1-1+1-1()()()()()()2nnnnkkxkxknnnnkkxkxknnnnkkxkxkuxuuuuxuuuuxuuuux空间方向的一阶偏导数可以近似为的向前差商：的向后差商：的中心差商：1111nnnxkkknnnxkkknnntkkkuuuuuuuuu空间方向的向前差分、向后差分和中心差分记为xxx其中，，，分别称为空间方向前差、后差和中心差分算子。后文中，将略去差分算子的下标，简记为，，。(,)(,)(2.1.3)(2.1.3)txxukxntukxnt中的二阶偏导数应该如何近似呢？2220(,)2(,)(,)(,)limxuuxxtuxtuxxtxtxx根据数学分析中的知识，我们知道21222nnnnkkkkuuuuxx所以，二阶导数可以近似为112nnnkkkuuu称为二阶中心差分。112=()()nnnnnkkkkkuuuuu容易证明：2.1.3差分格式同一偏导数可以有不同的近似方法，不同的导数近似方法导致方程的不同的有限差分近似。1.FTCS(ForwarddifferenceinTime,CentraldifferenceinSpace)格式时间方向用前差近似，空间二阶导数用中心差分近似。11122(2.1.9)nnnnnkkkkkuuuuutx对初始条件和边界条件的离散化00()(0,1,,)(2.1.10)()(0,1,)(2.1.11)()(0,1,)(2.1.12)kknnnMnufxkMuatnubtn式(2.1.9)~(2.1.12)称为方程(2.1.1)的一个有限差分方程或有限差分格式(finitedifferencescheme)。2.BTCS(BackwarddifferenceinTime,CentraldifferenceinSpace)格式时间方向用后差近似，空间二阶导数用中心差分近似。11122(2.1.13)nnnnnkkkkkuuuuutx00()(0,1,,)()(0,1,)()(0,1,)kknnnMnufxkMuatnubtn在研究数值方法时，通常把tn时刻的物理量视为已知量，而把tn+1时刻的物理量作为待求的未知量。因此，式(2.1.13)可以改写成11111122(2.1.14)nnnnnkkkkkuuuuutx2.1.4差分方程的求解1.FTCS格式11122(2.1.9)nnnnnkkkkkuuuuutx可以改写为111(12)(2.1.15)nnnnkkkkuuuu2tx其中，可见，在FTCS格式中，某一点的数值解只依赖于前一时间步的三个点，如图2.2所示。图2.2：FTCS格式的模板点1(stencil)nku求解所涉及的典型网格点称为模板格式的点。1expliciFTCS2.1.15tn1schemenknknuu格式可以通过简单的递推关系由某一时间步的值计算出下一个时间步的值式，称为显示格。FTCS格式的求解过程01.0()()(0,1,,)nkkkknufxufxkM赋初始值令,由计算1+1112.(12)(0,1,,)nnnnnkkkkkuuuuukM内点数值解计算由计算+1+1003.()()(0,1,)nnnnnMnMuatubtnuu边界处理由，计算，1nn令4.ntT判断是否成立5.输出结果成立不成立2.BTCS格式11111122(2.1.14)nnnnnkkkkkuuuuutx可以改写为11+1+1-1-(12)+=-(2.1.16)nnnnkkkkuuuu跟FTCS格式不同，BTCS格式中同时涉及到n+1时刻的多个未知量，不能递推求解，称为隐式格式(implicitscheme)。图2.3：BTCS格式的模板点1.0nknu赋初始值(令,计算)12.nku构造求解的线性方程组3.求解线性方程组1nn令4.ntT判断是否成立5.输出结果成立不成立BTCS格式的求解过程12.1nknu构造求解时刻数值解的线性方程组11+1-1+1-(12)+=-(2.1.16)nnnnkkkkuuuu0,1,,kM列出各点差分格式的具体形式10111+1012111+1123211+1211+110:()()1:-(12)+=-2:-(12)+=-...1:-(12)+=-:=b()()nnnnnnnnnnnnnnMMMMnMnkuatkuuuukuuuukMuuuukMut边界条件边界条件写成方程组的形式1n通过求解这个线性方程组，可以得到时刻求解域上各个网格点的数值解。3.求解线性方程组系数矩阵只有主对角线和相邻的两条次对角线上有非零元素，这种形式的方程组称三对角线为方程组。可以通三对角线方程组追过赶法直接求解。：I.利用一个边界条件将三对角线方程组化为只有主对角线和相邻的一条次对角线上有非零元素的方程组；II.利用另一边界条件追赶法逐点求解。1M：是一种高效算法，计算量与未知量个数近似呈线追赶法性关系。2.1.5用时间相关方法求解定常问题考虑非定常热传导方程和定解条件22(0)(,0)()(2.1.18)(0,)(1,)uutxuxfxutaconstutbconst(2.1.18)t当时，的解与时间无关，与下面的定解问题等价。22=0(,0)()(2.1.20)(0)(1)uxuxfxuaconstubconst2.1.202.1.18显然，定常问题的解也可以通过数值求解得到。1-6(2.1.18),110nnkkuuktn在实际求解过程中，无需计算无穷多时间步，只要定解问题的数值解满足即可认为时刻的数值解为定常解。其中，是一个小的正实数，根据对定常解精度的要求事先指定，定常问如。上题的时述方法称为求解间相关方法。1.nku赋初始值(令n=0,计算)12.nku构造求解的线性方程组3.求解线性方程组1nn令-64.10判断是否成立5.输出结果成立不成立BTCS格式的求解过程FTCS格式的求解过程01.()()(0,1,,)nkkkkufxufxkM赋初始值令n=0,由计算1+1112.(12)(0,1,,)nnnnnkkkkkuuuuukM内点数值解计算由计算+1+1003.()()(0,1,)nnnnnMnMuatubtnuu边界处理由，计算，1nn令5.输出结果成立不成立-64.10判断是否成立实用中，通常采用下面的求解过程§2.2导数的数值逼近方法2.2.1精度分析在上一节，我们得到了一阶偏导数的前差、后差和中心差分近似，以及二阶导数的中心差分近似。这些近似方法逼近偏导数的程度如何呢？可以用Taylor展开式进行分析。,(,)(,)ijijuxyuixjyu记22331,,23(,)(,)(,)22331,,23(,)(,)(,)Taylor2!3!2!3!ijijijijijijijijijijuuxuxuuxxxxuuxuxuuxxxx由展开式，,1,,(,)23223(,)(,),(,)....()2!3!..()()xijijijijijijxijijuuuuTExxxuxuxTEOxxxTETruncationErrorOxuuxx：其中，称为截断向前差误差，是的量级。商一阶精度的差称是的分近似。,,-1,(,)23223(,)(,),(,)....-()2!3!xijijijijijijxijijuuuuTExxxuxuxTEOxxxuuxx