《计算流体力学2011-ch3 [兼容模式].pdf》由会员分享，可在线阅读，更多相关《计算流体力学2011-ch3 [兼容模式].pdf（15页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、1Chapter 3 有限差分方法有限差分方法Finite Difference Method(FDM)Finite Difference Method(FDM)3.1.3.1.引言引言3.1.3.1.引言引言3.1.1.有限差分基本慨念3.1.1.有限差分基本慨念离散的方式表示连续的流体：用在空间和时间有限网格点上流体的状况近似地描述连续流体的状况。FDM 是所有近似方法中最简单的一种，在FDM中：离散的方式表示连续的流体：用在空间和时间有限网格点上流体的状况近似地描述连续流体的状况。FDM 是所有近似方法中最简单的一种，在FDM中：微分?差分微分方程?代数方程组求解微分方程?求解代数方程组
2、微分?差分微分方程?代数方程组求解微分方程?求解代数方程组将时间和空间划分成网格后，本章中运用下面的基本符号时间步数（Time level）n:n表示时间步n的t=时间间隔=tn+1 tn.n-1:past;n:present;n+1:futuret=n t，n=时间步数=0,1,2,3,.,NT=N t=最后时间T N t 最后时间空间位置:下标i,j,k,（for x,y,and z）.x，y,z:x，y,z向网格距坐标：xi=ixx，y,z和t 经常为常数Note:离散化后会丢失信息 the greater the number of points,the more accurate w
3、ill be the representation?See Figure.3.1.2.3.1.2.数值算法中需要考虑的几个定量性质数值算法中需要考虑的几个定量性质3.1.2.3.1.2.数值算法中需要考虑的几个定量性质数值算法中需要考虑的几个定量性质流体控制方程(PDEs)具有的一些物理性质希望在离散化后仍然保持。比如1)守恒性守恒性 控制方程经常表现出一些物理量的守恒性质(在一个封闭的空间物理量的积分不随时间改变).质量守恒E.g.,质量守恒在一个封闭的箱体积分离散化后是否保持守恒？答案：不一定。比如，还可以将方程写成考虑网格点上 和 v 是分离的(staggered grids)，一个小网
4、格带内的质量变化是通过边界的流出（入）引起的，为了计算边界处须将 平均到点上为了计算边界处v，须将平均到u,v点上。2如果用另外的计算网格(non-staggered grids):为了计算网格边界处通量又要采用不同的平均方法,得到的守恒结果将会不同。因此，为了保证守恒性要仔细设计网格和算法。因此，为了保证守恒性要仔细设计网格和算法。2)非负性非负性 一些正的物理量(mass,energy,water vapor)不能变成负的.数值计算过程中不一定能够保证保证非负性的计算方案称为正定格式线性内插方案不会产生新的极值避免含有外插性质的格式3)可逆性3)可逆性 对输送过程应该可逆性，但是扩散过程没
5、有。实际很难达到：计算误差不可避免。4)精度对输送过程应该可逆性，但是扩散过程没有。实际很难达到：计算误差不可避免。4)精度 影响精度的因素：截断误差；计算机精度；影响精度的因素：截断误差；计算机精度；数值解的稳定性数值解的稳定性。数值解的稳定性数值解的稳定性。参考文献：参考文献：1、王斌，季仲贞：大气科学中数值新方法及其运用，科学出 版社，、王斌，季仲贞：大气科学中数值新方法及其运用，科学出 版社，20052、廖洞贤：大气数值模式的设计，气象出版社，、廖洞贤：大气数值模式的设计，气象出版社，19993.2.3.2.有限差分表达式的产生方法有限差分表达式的产生方法3.2.3.2.有限差分表达式
6、的产生方法有限差分表达式的产生方法两种常用方法:两种常用方法:1.Taylor 级数展开 （最常用）；1.Taylor 级数展开（最常用）；2.多项式拟合（插值）最一般性的方法，Taylor 级数展开可以看为其子方法。3.2.1.3.2.1.3.2.1.3.2.1.Taylor Taylor 级数展开方法级数展开方法Taylor Taylor 级数展开方法级数展开方法给出 u(x0),对于u(x0+x)可以写出 Taylor 级数展开=+=+33322200000.!3)(!2)()()(xxxxxxxuxxuxxuxxuxxu给出有限差分表达式是反过来：对有限的x 给出u/x 的近似表达式=
7、00!)(nxxnnnxunx记:ui+1=u(x0+x).由 Taylor 级数表达式得到u/x：右端第一项是 u(x)在点(x0)和(x0+x)之间的平均斜率 可以作为该点斜率的近似1)(.!3)(!2)()()(0003322200+=xxxxxxxuxxuxxxuxxuxu率,可以作为该点斜率的近似：比较偏导数定义)()()(000 xoxxuxxuxuxx+=xxuxxuxuxxx+=)()(lim00003类似可以用(x0 x)的值代替得到(1)称为向前差分（forward difference）(2)称为向后差分（backward difference）2)(.!3)(!2)()
8、()(0003322200+=xxxxxxxuxxuxxxxuxuxu如果有方程u/t+cu/x=0,c 0,可以用 u/t+c(ui-ui-1)/x=0 称为上游差分格式或迎风格式（Upwind Difference）；也可以是u/t+c(ui+1-ui)/x=0（下游格式）。对这个问题上游差分格式比下游格式好,原因是和信号的传递方向一致。从特征线看，没有信息从+x 方向传递过来。如果 c 2x=xi+1xi-1,得到xuuxuii=+2)(11自己验算。0 ,)2(3321122=+=+xuxuuuxuiii5气象.1990.016(005).11-14,42 函数F(x)一阶导数的5点对
9、称差分格式的一般形式为（7）s 是待定系数。用傅里叶级数拟合F(x)：F(x)A exp(imx)（8）波数2/L对于特定波长 Lkh可以选取 使得差分格式的hFFshFFsdxdFjjjjj4/)(2/)(1()(2211+=波数 m=2/L。对于特定波长 Lk=kh，可以选取s使得差分格式的T.E为零。由（8）式得 dF/dx=imF。将（8）代入（7）得到(dF/dx)j=Rk(imFj)，这里让R=1,则T.E0)/4/()/4sin()/2sin()1(2kksksRk+=)/2sin(2)/4/sin()/2sin()/2(2kkkks=解得：问题：k该取多大为好？Sk=-0.46
10、50这一组方案是精度最高的一种，称为最优五点差分格式。这一组方案是精度最高的一种，称为最优五点差分格式。3.2.3.3.2.3.插值方法和特征线插值方法和特征线3.2.3.3.2.3.插值方法和特征线插值方法和特征线本质上讲，插值方法是一切数值方法的核心，下面就从插值方从插值方法角度来构造差分格式法角度来构造差分格式例如1-D 平流方程:（c 0 and constant)0=+xuctu画一个x-t平面的图注意：沿特征线dx/dt=c,du=0。目标:决定)etc,(111ninininiuuufu+=根据特征线上的守恒性质：一般根据特征线上的守恒性质：一般 x*不会落在计算格点上，可通过周
11、围点插值得到。选用不同的插值方式可得到不同精度的格式。如果用最简单的线性内插可以得到不会落在计算格点上，可通过周围点插值得到。选用不同的插值方式可得到不同精度的格式。如果用最简单的线性内插可以得到(see figure),(*1nnnitxuuu=+现在，一个时间步长现在，一个时间步长t内以速度内以速度c移动距离为移动距离为x,即即x=c t.代入上式得到代入上式得到niniuxxuxxu)1(1*+=011=+xuuctuunininini6得到一个上游差分格式，以后会看到它有很强的衰减作用，因为线性内插会削弱峰值。如果 c=constant and t=x/c，得到的解没有误差。（难于碰到
12、）如果将点i+1也包含在插值公式内得到xtcuuuuuuuninininininini=+=+/where)2(2)(2112111称为 Lax-Wendroff scheme（“三条腿”格式）注意：右边最后一项是一个扩散项的差分表达式，它对保证计算稳定是必要的。如果没有这项，成为时间向前差，空间中央差格式（FTCS)（forward-in-time and centered-in space），以后可以看到是不稳定的。课堂练习：分析上式的截断误差课堂练习：分析上式的截断误差。稳定条件和对时间步长的限制稳定条件和对时间步长的限制从上面的插值计算可以看到，保证插值是内插而非外插的条件是它限制时间
13、步长：t 0 and constant)FTCS scheme:12uuuuunnnn+运用Taylor级数展开方法决定它的截断误差它的截断误差.0 when x 0 and t 0，i.e.,the scheme is consistent.2112xuuutuuiiiii+=+),(22222txotut+=一个反例:Dufort-Frankel 格式：时间和空间都是中央差（2阶精度）,计算截断误差：2111111)()(2xuuuutuunininininini+=+H.R.T 6)()()(12)(3322222442+=tuttuxtxux“空心”菱形格式空心”菱形格式可以看出除第二
14、项外，t,x0 时 0。如果则第二项也满足一致性要求，因此必须要求 t 比 x更快0。否则，假设有0lim0,=xttx=xttx0,lim)(于是方程变成于是方程变成222uuu=+2220,0 limtutx=变成一个完全不同的方程:带阻尼的波动方程带阻尼的波动方程！22 xtt+3.3.2.3.3.2.收敛性收敛性收敛性收敛性仅一般讨论，定义：一般情况收敛性不好证明，特别是非线性方程),(lim0,ninitxtxuu=般情况收敛性不好证明，特别是非线性方程。Laxs 定理对于说明线性方程的收敛性很有用，也常常推广到非线性方程。首先以扩散方程为例来证明收敛性。1-D 扩散方程的收敛性证明
15、扩散方程的收敛性证明微分方程 u/t=2u/x2，其中 x 0，L，大于0常数I.C：B.C.：u(0,t)=u(L,t)=0 是线性的、适定的初值问题。)()/sin()0,(1xfLxkaxukk=a、先求出解析解：、先求出解析解：由于是线性的，只探讨如下形式的一个单波解:)/sin()(),(LxktAtxukk=代入方程得到振幅:由 I.C.,Ak(0)=ak,于是得出单波解析解：并记B、考虑考虑(FTCS)格式格式)/(exp)0()(2tLxAtAkk=)/sin()/(exp),(2LxktLxatxukk=),(),(txutxukk=（9）B、考虑考虑(FTCS)格式格式目标
16、:证明),(lim0,ninitxtxuu=第一步求数值解：数值解：JiLxktAtxuinkJkni,.2,1,0for )/sin()(),(0=8将离散的Fourier 级数表达式代入方程，先让 I.C.离散化为总格点数为 J+1。通过离散Fourier变换得到系数:kaJiLxkaxuxfikJkii,.2,1,0for )/sin()0,()(0=注意1注意1:L=J x，随x 0,J,离散的Fourier 级数变成连续的，并且kkaa JkLxkxfallJlk,.2,1,0for )/sin()(J20=注意2注意2:在网格空间所能表现的Fourier 级数的波数是格点数(J+1
17、)（自由度）的函数。由上面的表达式中，不难得到波长和 波数k的关系：kLLkLxk/22)/()/sin(=即2L/k最长的波=,相当于波数k=0，下一个最长的波=2L (k=1)，最短的波(k=J):=2L/J=2Jx/J=2x。结论结论:波长为2x的波是能被分辨的最短的波（至少要3个格点的值才能表示一个波）。2x格距的波在数值计算中是描写最差的（越光滑误差越小），也常常有一些特别的性质。C、考虑其中一个波解、考虑其中一个波解(波数为k,满足B.C.)：n是时间步数，要求满足I.C.)/sin()()(LxknAuikkni=kkaA)0(=21112xuuutuuninininini+=+
18、代入FDE记我们有 )/sin(Lxksii=211)(2)()()1(xsssnAtsnAsnAiiikikik+=+ikknisnAu)()(=（10）由于xi=ix,xi+1=(i+1)x,得代入（10）式，利用三角函数关系得到/)(sin1Lxxksii+=+)2(sin41)()1(2xknAnAkk=+2)(where,xt=)2/(sin41)(2LxkkM如果，记那末我们有对n=0,1,2,n写出：)()()1(nAkMn Akk=+kkkkkkakMAkMAakMAkMA)()1()()2()()0()()1(2=得到差分解knkkakMnAkMnA)(.)1()()(=)/
19、sin()()()(LxkakMsnAuiknikkni=这里M(k)是一个振幅因子，如果|M(k)|1，解不会随n 而增长，这是我们期望的（微分方程的解随时间衰减的），要求|M(k)|1，意味着如果取最大可能值i2()1得到1)2/(sin41 2Lxk如果取最大可能值sin2()=1,得到?11-41 得到?1/2.（对所有k）根据的定义,这个条件成为以后我们看到它是保持稳定性的条件。22xt9现在检验它的收敛性：或者，如果我们记)/(sin)/(sin41)(2220.0.limlimLxkLxkxtauiniktxknitx=)/(sin4)(22Lxkxxf)/(sin)(1)(20
20、0limlimLxkxtfauinktxknitx=因为所以0.0.txtx1sinlim0=yyyaLkLxkLxkLkxfxx=22200)/(2/)2/sin()/(lim)(lim可以证明，如果f(x)是一个以实数（比如x)为自变量的复函数,则有（证明见后面）ttnextfatntx=其中,)(1 lim0.此外，kkxaa=lim0于是)2/sin()/(exp)(lim20,LxktLkauikknitx=收敛性得到证明 与单波解析解（9）比较atntxextf=)(1 lim0.公式补充证明公式补充证明3.3.3.3.3.3.数值收敛性数值收敛性3.3.3.3.3.3.数值收敛性
21、数值收敛性要从理论上论证收敛性很困难，因为解析解一般得不到，然而我们可以通过数值方法找出一个差分格式的收敛性。基本思想是不断提高网格的分辨率，看看误差如何随分辨率变化:x 0，0？以及减少得有多快？（可能会花大量计算时间，当x很小时）。一般我们用L2模或者 均方根误差RMS(root mean square error)来表述误差：这里ui是真值或者是“收敛的”数值解（如果没有真值）。nuuLiixi/)(2,2=Fletcher book section Fletcher book section on numerical convergence 10此处此处 s 即为上面定义的即为上面定义
22、的，按稳定性条件要小于等于，按稳定性条件要小于等于0.50f(x,3.3.4.3.3.4.稳定性分析稳定性分析稳定性稳定性分析稳定性分析稳定性：对于 一个稳定的数值格式，初始时刻的误差不会随时间无界增长。几种判定方法几种判定方法:能量方法 von Neumann 方法 矩阵方法(对方程组)直接扰动法(不讨论)能量法能量法比von Neumann 方法用的少。好处：可以用于非线性问题，不要求周期边条件。证明的关键：证明一些非负量，比如对所有的时间n 有界。2)(niu以平流方程u/t+cu/x=0的upstream-forward scheme 来说明0)(11=+xuuctuunininini
23、(11)()1(2)(1 ()(:)1(/2121222111ninininiiniinininiuuuuuuuuxtc+=+=求和两边平方，对所有各点让212 )()(B.C.niniuu=假定周期性22222122121)()()1(2)1()(,(11),0)-(1 )()()()()(Schwarz niininiiniiiniinininiiiiuuuuuuuu=+=+得到右端所有项系数是正的如果和对于两个向量不等式应用VUVUV,U11方案稳定!条件(1-)0 给出=ct/x 1正如前面指出的，t内信号传播的距离不超过x。von Neumann 方法方法其实我们早已用过这一方法确定
24、1-D扩散方程FTCS格式的稳定性（收敛性一节中）。我们发现un+1=M(t)n+1u0M(t)=振幅因子。如果|M(t)|1，得到un+1 un解不会随时间增长，稳定！von Neumann 方法：将解按照Fourier 级数展开，找出振幅因子，看看它满足 1 的条件。假定1.方程是线性的，并且是常系数；2.解是周期性的。（要按照Fourier 级数展开）:U是振幅（可以是复数），=R+ii是频率.(12)(exp),(,tmzlykxiUtzyxumlkmlk+=实际上R给出波的传播速度，i给出增长（衰减）率如果i0,解随时间指数增长.t-itti-itiRiiReeee=+)(差分方程得
25、：将上面表达代入以及的一个波考虑波数为这里格式扩散方程，titkxikxiktkxiknititkxiknitkxikninininininixkeeUeeUueeUueUukxtuuuuu+=+=+=0)2/(sin41 :)/()2(:FTCS D-1 1 Example2)()(11)(1)(2111。振幅因子，相当于前面的是实际的振幅因子这里对于非零解，要求tininitititikeuuexkexkexkeeU+=+/M ,)2/(sin410)2/(sin410)2/(sin41122实际上=1/2时候,对于2x波(k=22 xx，sin 2(k x/2)=1)(振幅因子1)每一步
26、解在和之间摆动 这是不真实的因此般会)!(1/2 1)2/(sin41112同样的结果稳定要求xk解在1 和+1 之间摆动,这是不真实的。因此，一般会要求1/4。在稳定条件下，扩散项会起到衰减短波的作用，然而不满足稳定条件时，他自己造成不稳定。3.3.5.3.3.5.隐式方法隐式方法3.3.5.3.3.5.隐式方法隐式方法前面讨论的差分方程可以归结为称为显式格式：下一步每个格点的值只依赖于当前时间或者过去时间的值，因而对每一点的方程可以独立求解。隐式格式：,.),(11+=nnniuufu隐式格式：它类似一个隐式方程（如x=sin(x)）可以想像隐式格式求解更困难，解一个联立的代数方程组。,.
27、),(,.),(.,111111+=nininininiuufuuug12Example.一维扩散方程F.D.方程:=0,FTCS scheme，显式的=1/2,隐式的（Crank-Nicolson scheme）.)x)/2(:(Note13 )1(/)(21111+=+=iiixxnixxnixxniniuuuuuKtuu符号）（对于其他的,一般是隐式的.可以证明如果=1/2,时间空间都是2阶精度.否则时间是一阶精度。如果 =0,变成空间4阶精度.运用von Neumann方法分析(13)的稳定性.看几种情况=(13)(14)j .)14()(代入将表示格点号。现在用下标这里tixikjnkxikjtniktkxiknjeeUeeUeUu)2/(sin41)2/(sin)1(4122xkxk+=看几种情况:绝对稳定（无条件稳定）.成立所有的对定）（和前面相同，条件稳显示格式 1)2/(sin21)2/(sin21Nicolson)-(Crank 1/2 :II Case2/1)2/(sin41),0(:I Case222+=xkxkxkCase III:=1,相对于右边项全隐式，时间是向前差。然而对时间只是一阶精度(非中央差)。,1)2/(sin2112也是绝对稳定的对所有+=xk一般情况：0 0为常数为常数。