有限差分方法

1有限差分法2†有限差分法是数值求解微分问题的一种重要工具，在很多领域（如传热分析、流动分析、扩散分析等）取得了显著的成就；†有限差分法在材料成形领域的应用较为普遍，与有限差分法一起成为材料成形计算机模拟技术的主要两种数值分析方法。目前材料加工中的传热分析（如铸造成形过程的传热凝固、塑性成形中的传热、焊接成形中的热量传递等）、流动分析（如铸件充型过程，焊接熔池的产生、移动，激光熔敷中的动量传递等）都可用有限差分法进行模拟分析。†特别是在流动场分析方面，与有限元相比，有限差分法有独特的优势。†被认为是有限差分法的弱项——应力分析，也取得了长足进步。一些基于差分法的材料加工领域的应力分析软件纷纷推出。†使得流动、传热、应力统一于差分方式下。3第一节差分原理及逼近误差/差分原理1．差分原理设有x的解析函数y=f(x)，从微分学知道函数y对x的导数xxfxxfxydxdyxxΔ−Δ+=ΔΔ=→Δ→Δ)()(limlim00(1-1)dxdy是函数对自变量的导数，又称微商；dy、dx分别是函数及自变量的微分。yΔ、xΔ分别称为函数及自变量的差分，xyΔΔ为函数对自变量的差商。4第一节差分原理及逼近误差/差分原理向前差分)()(xfxxfy−Δ+=Δ)()(xxfxfyΔ−−=Δ)21()21(xxfxxfyΔ−−Δ+=Δ向后差分中心差分x（1-2）（1-3）（1-4）Δ〉0yyyyyy5第一节差分原理及逼近误差/差分原理上面谈的是一阶导数，对应的称为一阶差分。对一阶差分再作一阶差分，所得到的称为二阶差分，记为。y2Δ以向前差分为例，有)()(2)2()]()([)]()2([)()()]()([)(2xfxxfxxfxfxxfxxfxxfxfxxfxfxxfyy+Δ+−Δ+=−Δ+−Δ+−Δ+=Δ−Δ+Δ=−Δ+Δ=ΔΔ=Δ（1-5）6第一节差分原理及逼近误差/差分原理依此类推，任何阶差分都可由其低一阶的差分再作一阶差分得到。例如n阶前差分为)]}()(([{)]}([{)]([)(21xfxxfyyyynnn−Δ+ΔΔΔ=ΔΔΔΔ=ΔΔΔ=ΔΔ=Δ−−ΛΛΛΛ（1-6）7第一节差分原理及逼近误差/差分原理函数的差分与自变量的差分之比，即为函数对自变量的差商。一阶向前差商为xxfxxfxyΔ−Δ+=ΔΔ)()(一阶向后差商为xxxfxfxyΔΔ−−=ΔΔ)()(（1-7）（1-8）8第一节差分原理及逼近误差/差分原理一阶中心差商为xxxfxxfxyΔΔ−−Δ+=ΔΔ)21()21(或xxxfxxfxyΔΔ−−Δ+=ΔΔ2)()(（1-9）（1-10）9第一节差分原理及逼近误差/差分原理二阶差商多取中心式，即222)()()(2)(xxxfxfxxfxyΔΔ−+−Δ+=ΔΔ当然，在某些情况下也可取向前或向后的二阶差商。（1-11）10第一节差分原理及逼近误差/差分原理以上是一元函数的差分与差商。多元函数f(x,y,…)的差分与差商也可以类推。如一阶向前差商为,),,(),,(xyxfyxxfxfΔ−Δ+=ΔΔΛΛΛΛΛΛ,),,(),,(yyxfyyxfyfΔ−Δ+=ΔΔ（1-12）（1-13）11第一节差分原理及逼近误差/逼近误差差商与导数之间的误差表明差商逼近导数的程度，称为逼近误差。由函数的Taylor展开，可以得到逼近误差相对于自变量差分（增量）的量级，称为用差商代替导数的精度，简称为差商的精度。))(()(!4)()(!3)()(!2)()()()(5432xOxfxxfxxfxxfxxfxxfIVΔ+Δ+′′′⋅Δ+′′⋅Δ+′⋅Δ+=Δ+（1-14）)()())(()(!4)()(!3)(!2)()()()(432xOxfxOxxfxxfxxfxfxxfxxfIVΔ+′=Δ+Δ+Δ′′′+Δ′′+′=Δ−Δ+∴（1-15）一阶向前差商2．逼近误差一阶向前差商具有一阶精度O()表示与()中的量有相同量级的量。将函数f(x+Δx)在x的Δx邻域作Taylor展开12第一节差分原理及逼近误差/逼近误差)()()()(),)(()(!4)()(!3)()(!2)()()()(5432xOxfxxxfxfxOxfxxfxxfxxfxxfxxfIVΔ+′=ΔΔ−−∴Δ+Δ+′′′Δ−′′Δ+′⋅Δ−=Δ−（1-16）一阶向后差商：一阶向后差商也具有一阶精度。()()nOxΔxΔ中的指数n作为精度的阶数是个小量，因此阶数越大则表明精度越高。xΔ13第一节差分原理及逼近误差/逼近误差将)(xxfΔ+与)(xxfΔ−的Taylor展开式相减可得))(()(2)()(2xOxfxxxfxxfΔ+′=ΔΔ−−Δ+（1-17）可见一阶中心差商具有二阶精度。一阶中心差商))(()(2)()(2xOxfxxxfxxfΔ+′=ΔΔ−−Δ+一阶中心差商14第一节差分原理及逼近误差/逼近误差))(()()()(2)(22xOxfxxxfxfxxfΔ+′′=ΔΔ−+−Δ+这说明二阶中心差商的精度也为二阶（1-18）将)(xxfΔ+与)(xxfΔ−的Taylor展开式相加可得二阶中心差商15第一节差分原理及逼近误差/非均匀步长2−Δix在有些情况下要求自变量的增量本身是变化的，如图2-1中的1−Δix1+ΔΔiixx和，是不相等的，相应的差分和差商就是不Ox2−ix1−ixix1+ix2+ix2−Δix1−ΔixixΔ2+Δix图1-1非均匀步长差分3．非均匀步长一阶向后差商11)()(−−ΔΔ−−iiiixxxfxf一阶中心差商11)()(−−Δ+ΔΔ−−Δ+iiiiiixxxxfxxf（1-22）（1-23）不等距的。基于不等距离（一维）、不等规格（二维、三维）网格：变网格技术可以在保证计算精度的前提下有效地提高计算速度。16第一节差分原理及逼近误差/非均匀步长图1-2均匀和非均匀网格实例117第一节差分原理及逼近误差/非均匀步长图1-3均匀和非均匀网格实例218第二节差分方程0=∂∂+∂∂xtζαζ差分相应于微分，差商相应于导数。差分和差商是用有限形式表示的，而微分和导数则是以极限形式表示的。如果将微分方程中的导数用相应的差商近似代替，就可得到有限形式的差分方程。现以对流方程为例，列出对应的差分方程。（2-1）：物质的量；为运动速度。ζααα为常数，一维常系数对流方程不为常数，一维变系数对流方程19Λ,2,1,0,0=Δ+=ixixxiΛ,2,1,0,=Δ=ntntn图2-1差分网格划分矩形网格用差商近似代替导数时，首先选定Δx和Δt，称为步长。然后在x-t坐标平面用平行于坐标的两族直线：通常空间步长Δx取为相等的，而时间步长Δt与Δx以及α有关，当Δx和α为常数时，Δt也取常数。步长定律第二节差分方程1txα⋅Δ=Δ直线t=tn称为第n层。网格交叉点称为节点。20若时间导数用一阶向前差商近似代替，即ttnininiΔ−≈⎟⎠⎞⎜⎝⎛∂∂+ζζζ1空间导数用一阶中心差商近似代替，即xxnininiΔ−≈⎟⎠⎞⎜⎝⎛∂∂−+211ζζζ则在),(nitx点的对流方程就可近似地写作02111=Δ−+Δ−−++xtninininiζζαζζ（2-3）（2-4）（2-2）第二节差分方程21第二节差分方程按照前面关于逼近误差的分析知道，用时间向前差商代替时间导数时的误差为,)(tOΔ用空间中心差商代替空间导数时的误差为，))((2xOΔ因而对流方程与对应的差分方程之间也存在一个误差，它是))(,())(()(22xtOxOtORniΔΔ=Δ+Δ=这也可由Taylor展开得到。因为))(,()(!31212),(),(),(),(223322xtOxtxtxtttxtxxtxxttxttxnininininininininiΔΔ+⎟⎠⎞⎜⎝⎛∂∂+∂∂=⎥⎥⎦⎤⎢⎢⎣⎡+Δ⎟⎟⎠⎞⎜⎜⎝⎛∂∂+⎟⎠⎞⎜⎝⎛∂∂++Δ⎟⎟⎠⎞⎜⎜⎝⎛∂∂+⎟⎠⎞⎜⎝⎛∂∂=ΔΔ−−Δ++Δ−Δ+ζαζζζαζζζζαζζΛΛ（2-5）（2-6）22第二节差分方程按照前面关于逼近误差的分析知道，这种差分方程近似代替微分方程所引起的误差，称为截断误差。这里的误差量级相当于Δt的一次式、Δx的二次式。一般情况下，则说对Δt精度一阶，对Δx精度二阶。23第二节差分方程一个与时间相关的物理问题，应用微分方程表示时，还必须给定初始条件，从而形成一个完整的初值问题。对流方程的初值问题为⎪⎩⎪⎨⎧==∂∂+∂∂)()0,(0xxxtζζζαζ这里)(xζ为某已知函数。同样，差分方程也必须有初始条件：⎪⎩⎪⎨⎧==Δ−+Δ−−++)(020111iininininixxtζζζζαζζ初始条件是一种定解条件。如果是初边值问题，定解条件中还应有适当的边界条件。差分方程和其定解条件一起，称为相应微分方程定解问题的差分格式。（2-7）（2-8）24第二节差分方程⎪⎩⎪⎨⎧=−ΔΔ−=−++)()(20111iininininixxtζζζζαζζFTCS格式（2-9）⎪⎩⎪⎨⎧=−ΔΔ−=++)()(011iininininixxtζζζζαζζFTFS格式（2-10）⎪⎩⎪⎨⎧=−ΔΔ−=−+)()(011iininininixxtζζζζαζζFTBS格式ForwardTimeCentralSpaceForwardTimeForwardSpaceForwardTimeBackwardSpace（2-11）25显式有限差分计算模板⎪⎩⎪⎨⎧=−ΔΔ−=−++)()(20111iininininixxtζζζζαζζ26第二节差分方程(a)FTCS(b)FTFS(c)FTBS图2-2差分格式27第二节差分方程FTCS格式的截断误差为))(,(2xtORniΔΔ=FTFS和FTBS格式的截断误差为),(xtORniΔΔ=3种格式对tΔ都有一阶精度。（2-12）（2-13）28例：考虑长度为1的均匀直杆，其表面是绝热的，而且杆截面足够细，可以把截面上的所有点的温度看成是相同的。X轴取为沿杆轴方向，x=0，x=1对应杆的端点，则杆内温度分布T(x,t)随时间变化由下面的传热方程来描述：()()()22,000,1001,100TxTtTTttxTβ∂∂=∂⎧⎪⎪⎪∂=⎨⎪=⎪⎪=⎩(0x1)计算时考虑210,0.1,0.5xtβ−=Δ=Δ=29时间导数用一阶向前差商近似代替：1nnniiiTTTtt+−∂⎛⎞≈⎜⎟∂Δ⎝⎠空间导数用二阶中心差商近似代替：211222nnnniiiiTTTTxx+−⎛⎞−+∂≈⎜⎟∂Δ⎝⎠1112(2)nnnnniiiiitTTTTTxβ++−Δ=−−+Δ210,0.1,0.5xtβ−=Δ=Δ=取那么则()11112nnniiiTTT+−+=+300.00.10.20.30.40.50.60.70.80.91.00.01000000000001000.5100500000000501001.010050250000025501001.510062.52512.500012.52562.51002.010062.537.512.56.2506.2512.537.562.51002.510068.837.521.96.256.256.2521.937.568.81003.010068.845.321.914.16.2514.121.945.368.8100计算结果31(1)请推导出y=f(x)的二阶向前差分，向后差分以及中心差分。(2)用中心差商数值微分公式计算在时的一阶导数，另用向前差商计算二阶导数值，取步长。(3)请推导出对流方程FTCS格式，并指出其截断误差的精度。（4）试推导出一维扩散微分方程（考虑D为常数）差分方程，并指出其截断误差。⎪⎩⎪⎨⎧==∂∂+∂∂)()0,(0xxxtζζζαζ作业()fxx=2x='(2)f''(2)f0.5h=220CCDtx∂∂−=∂∂