差分方法基础

第二讲 有限差分法基本原理
一般的流体控制方程都是非线性的偏微分方程。

在绝大多数情况下，这些偏微分方程无法得到精确解；而CFD 就是通过采用各种计算方法得到这些偏微分方程的数值解，或称近似解。

当然这些近似解应该满足一定的精度。

目前，主要采用的CFD 方法是有限差分法和有限体积法。

本讲主要介绍有限差分法，它也是下一讲中的有限体积法的基础[1]。

有限差分法求解流动控制方程的基本过程是：首先将求解区域划分为差分网格，用有限个网格点代替连续的求解域，将待求解的流动变量（如密度、速度等）存储在各网格点上，并将偏微分方程中的微分项用相应的差商代替，从而将偏微分方程转化为代数形式的差分方程，得到含有离散点上的有限个未知变量的差分方程组。

求出该差分方程组的解，也就得到了网格点上流动变量的数值解。

2.1 差分和逼近误差
由于通常数字计算机只能执行算术运算和逻辑运算，因此就需要一种用算术运算来处理函数微分运算的数值方法。

而有限差分法就是用离散网格点上的函数值来近似导数的一种方法。

设有x 的解析函数)(x f y =，从微分学知道函数y 对x 的导数为 x
x f x x f x y dx dy x x ∆-∆+=∆∆=→∆→∆)()(lim lim 00 （2-1） dy 、dx 分别是函数及自变量的微分，dx dy /是函数对自变量的导数，又称微商。

相应地，上式中的x ∆、y ∆分别称为自变量及函数的差分，x y ∆∆/为函数对自变量的差商。

在导数的定义中x ∆是以任意方式逼近于零的，因而x ∆是可正可负的。

在差分方法中，x ∆总是取某一小的正数。

这样一来，与微分对应的差分可以有三种形式：
向前差分 )()(x f x x f y -∆+=∆
向后差分 )()(x x f x f y ∆--=∆
中心差分 )2
1()21(x x f x x f y ∆--∆+=∆
上面谈的是一阶导数，对应的称为一阶差分。

对一阶差分再作一阶差分，就得到二阶差分，记为y 2∆。

以前向差分为例，有
[][][]
)
()(2)2()()()()2()()()()()
(2x f x x f x x f x f x x f x x f x x f x f x x f x f x x f y y +∆+-∆+=-∆+-∆+-∆+=∆-∆+∆=-∆+∆=∆∆=∆ （2-2）
依次类推，任何阶差分都可以由低一阶再作一阶差分得到。

函数的差分与自变量的差分之比，即为函数对自变量的差商。

如一阶向前差商为
x
x f x x f x y ∆-∆+=∆∆)()( 一阶向后差商为
x
x x f x f x y ∆∆--=∆∆)()( 一阶中心差商为
x
x x f x x f x y ∆∆--∆+=∆∆)21()21( 或
x
x x f x x f x y ∆∆--∆+=∆∆2)()( 二阶差商多取中心格式，即
222)
()()(2)(x x x f x f x x f x y ∆∆-+-∆+=∆∆
图2.1 差商与导数的关系
差商与导数的关系可见图2.1。

由导数（微商）和差商的定义知道，当自变量的差分（增量）趋近于零时，就可以由差商得到导数。

因此在数值计算中常用差商近似代替导数。

差商与导数之间的误差表明差商逼近导数的程度，称为逼近误差。

由函数的Taylor 展开，可以得到逼近误差相对于自变量差分（增量）的量级，称为用差商代替导数的精度，简称为差商的精度。

现以一阶向前差商为例来分析其精度。

将函数)(x x f ∆+在x 的x ∆邻域作Taylor 展开：
))(()(!
3)()(!2)()()()(43
2x O x f x x f x x f x x f x x f ∆+'''⋅∆+''⋅∆+'⋅∆+=∆+ 将上式代入一阶向前差商表达式中，有
)
()())(()(!
3)(!2)()()()(32x O x f x O x x f x x f x f x x f x x f ∆+'=∆+∆'''+∆''+'=∆-∆+ 这里符号)(O 表示与括号中的量有相同的量级。

上式表明一阶向前差商的逼近误差与自变量的增量为同一量级。

把)(n x O ∆中x ∆的指数作为精度的阶数。

这里1=n ，故一阶向前差商具有一阶精度。

由于x ∆是个小量，因此阶数越大精度越高。

采用同样的办法可知一阶向后差商也具有一阶精度。

对于一阶中心差商，将函数)(x x f ∆+与)(x x f ∆-在x 的x ∆邻域作Taylor 展开并代入一阶中心差商的表达式中，有
))(()(2)()(2x O x f x
x x f x x f ∆+'=∆∆--∆+ （2-3） 可见一阶中心差商具有二阶精度。

同样，二阶中心差商的精度也为二阶。

2.2 差分方程、截断误差和相容性
从上节所述可知，差分相应于微分，差商相应于导数。

只不过差分和差商是用有限形式表示的，而微分和导数则是以极限形式表示的。

如果将微分方程中的导数用相
图2.2 网格划分
应的差商近似代替，就可以得到有限形式的差分方程。

现以对流方程（2-4）为例，列出相对应的差分方程。

0=∂∂+∂∂x
u a t u （2-4） 用差商近似代替导数时，首先要选定x ∆和t ∆，称为步长。

然后在t x -坐标平面上用平行于坐标轴的两族直线：
,......
2,1,0,,......2,1,0,
0=∆==∆+=n t n t i x i x x n i 划分出矩形网格，如图2.2所示。

这里x ∆和t ∆取常数。

直线n t t =称为第n 层，网格交叉点称为结点。

网格点划定后，就可针对某一结点，例如图2.2中的结点),(n i t x ，用差商近似代替导数。

现用n i )(表示括号内函数在),(n i t x 点的值，则对流方程在该点为 0=⎪⎭⎫ ⎝⎛∂∂+⎪⎭⎫ ⎝⎛∂∂n
i
n i x u a t u （2-5） 如果时间导数用一阶向前差商近似代替： t u u t u n i n i n
i
∆-≈⎪⎭⎫ ⎝⎛∂∂+1 空间导数用一阶中心差商近似代替： x u u x u n i n i n
i
∆-≈⎪⎭⎫ ⎝⎛∂∂-+211 则对流方程在),(n i t x 点对应的差分方程为 02111=∆-+∆--++x u u t u u n i n i n i n i α （2-6） 按照前面关于逼近误差的分析知道，用时间向前差商代替时间导数的误差为)(t O ∆，用空间中心差商代替空间导数时的误差为))((2x O ∆，因而对流方程与对应的差分方程之间也存在一个误差，这一误差可由Taylor 展开确定，即 ))(,(...)(!31......212),(),()(),(223322x t O x u a t
u x x u x u a t t u t u x
t x x u t x x u a t t x u t t x u n
i n i n i n i n i n i n i n i n i ∆∆+⎪⎭⎫ ⎝⎛∂∂+∂∂=⎥⎥⎦
⎤⎢⎢⎣⎡+∆⎪⎪⎭⎫ ⎝⎛∂∂+⎪⎭⎫ ⎝⎛∂∂++∆⎪⎪⎭⎫ ⎝⎛∂∂+⎪⎭⎫ ⎝⎛∂∂=∆∆--∆++∆--∆+ （2-7） 这种用差分方程近似代替微分方程所引起的误差，称为截断误差。

这里误差量级相当
于t ∆的一次式、x ∆的二次式。

一个与时间相关的物理问题，应用微分方程表示时，还必须给定初始条件，从而形成一个完整的初值问题。

对流方程的初值问题为 ⎪⎩⎪⎨⎧==∂∂+∂∂)
()0,(0x u x u x u a t u （2-8） 这里)(x u 为某已知函数。

同样，差分方程也必须有初始条件： ⎪⎩
⎪⎨⎧==∆-+∆--++)(020111i i n i n i n i n i x u u x u u t u u α （2-9） 初始条件是一种定解条件，差分方程和其定解条件一起，称为相应微分方程定解问题的差分格式。

将式（2-9）中第1+n 时间层的量放在等号左边，将其余时间层的量放在等号右边，有 ⎪⎩
⎪⎨⎧=-∆∆-=-++)()(20111i i n i n i n i n i x u u u u x t a u u （2-10） 称其为FTCS 格式（时间前差、空间中差）。

若时间和空间都用向前差分，则得 ⎪⎩
⎪⎨⎧==∆-+∆-++)(0011i i n i n i n i n i x u u x u u t u u α （2-11） 同样，将第1+n 时间层的量放在等号左边，将其余时间层的量放在等号右边，有 ⎪⎩
⎪⎨⎧=-∆∆-=++)()(011i i n i n i n i n i x u u u u x t a u u （2-12） 该格式称为FTFS 格式。

若时间采用向前差分、空间采用向后差分，则得到FTBS 格式：
⎪⎩
⎪⎨⎧=-∆∆-=-+)()(011i i n i n i n i n i x u u u u x t a u u （2-13） 观察这三种差分格式，可以看出若知道第n 时间层的u ，则可以由一个差分式子直接算出第1+n 时间层的u ，称这类格式为显式格式。

差分方程的相容性：如果当x ∆、0→∆t 时，此差分方程的截断误差的某种范数也趋近于零，则表明从截断误差的角度来看，此差分方程是能用来逼近微分方程的，通常称这样的差分方程和相应的微分方程相容。

2.3 收敛性和稳定性
当步长趋于零时，要求差分格式的解趋于微分方程的解，称这种是否趋于微分方程定解问题的解的情况为差分格式的收敛性。

在有限差分法的具体运算中，计算误差总是不可避免的，如舍入误差，以及这种误差的传播、积累。

如果这一误差对以后的影响越来越小，或是这个误差保持在某个限度内，那么就称这个差分格式在给定的条件下稳定。

根据理论分析可以知道，上面介绍的几种差分格式是条件稳定的。

2.4 差分格式介绍
2.4.1 迎风格式
前面已经指出，微分问题
⎪⎩⎪⎨⎧==∂∂+∂∂)
()0,(0x u x u x u a t u 的FTBS 格式，在0>a 和1≤∆∆x t a 的条件下稳定，而FTFS 格式在0<a 和x
t a ∆∆≤-1的条件下稳定。

这里，当a 的符号改变时，为了使差分格式稳定，空间差分的方向也作了相应的变化。

由于a 是与速度对应的量，a 的正负表示速度方向的不同，即表示流（风）向：a 为正，看作风向沿着正方向吹；a 为负，则风朝着负方向吹。

而迎着方向往上游取空间差分，所得到的差分格式才可能是稳定的。

因为对流方程0>a 时的波形传播方向沿x 轴正向，上游的量经过一段时间要传播到下游，1+n 时刻i 站的量要受到上游站n 时刻量的影响，故只可以迎着风向取空间差分，而不可以顺着风向取空间差分。

这种格式是迎着风向往上游作差分所得到的，称为迎风格式。

上述FTBS 格式和FTFS 格式都必须在迎风时有条件稳定。

2.4.2 隐式格式
前面介绍的显式格式往往是有条件稳定的，甚至完全不稳定。

如FTCS 是完全不
稳定的，FTBS 格式是条件稳定的。

对于FTBS 格式，在0>a 和1≤∆∆x
t a 的条件下稳定，即要求a
x t ∆≤∆，当要求空间步长x ∆很小时，时间步长也必须取的很小，才能保证格式稳定，而t ∆取得小，计算工作量就大大增加，经济上也不合算。

而本节将要介绍的隐式格式常常是无条件稳定的，因此在许多情况下受到重视并被广泛应用。

隐式格式相当于从),(t t x n i ∆+点出发，用时间的向后差分把第1+n 时间层的量与已知时间层的量联系起来。

现以对流方程为例，从),(t t x n i ∆+点出发取BTCS 差分可得
0211111=∆-+∆-+-+++x
u u t u u n i n i n i n i α （2-14） 或改写为
n i n i n i n i u u x
t a u u x t a -=∆∆--∆∆++++-1111122 （2-15） 由于该方程含有三个第1+n 时间层上的函数值，即一个方程含有三个未知量，必须解联立方程才能得到第1+n 时间层上的未知量，故称该格式为隐式格式。

可以证明，用于对流方程的隐式格式是完全稳定的。

由于完全稳定，时间步长可以取得大些，从这一点来说，工作量减少了。

但隐式格式要解代数联立方程组，在每一时间步长内工作量有所增加。

2.5 耗散与色散
现以对流方程为例，采用时间向前差分、空间向后差分：
011=∆-+∆--+x
u u t u u n i n i n i n i α （2-16） 利用Taylor 展开，得到：
......)123(6)()1(21332222+∂∂--∆+∂∂-∆=∂∂+∂∂x
u r r x a x u r x a x u a t u （2-17） 上式就是差分方程（2-16）实际所模拟的微分方程，与原对流方程相比，多了二次导数项和三次导数项。

一般说：截断误差中含偶次导数项时，将引起耗散。

在流体力学方程中，二次导数项是与粘性项相关的。

不同的是，在这里该项是差分方程数值离散的结果，因此纯粹的数值引发，没有物理意义。

因此，在CFD 方法中，类似的项被称为数值耗散。

相似的，该项的系数，即)1(2
1r x a -∆，其作用类似于物理粘性，因而被称为人工粘性。

虽然数值耗散损害了计算精度，但却改善了计算的稳定性。

因此，很多不稳定的计算格式或显式、或隐含地添加了人工粘性后，就变得稳定了。

同样，在截断误差中还含奇次导数项，该项将引起数值色散。

最后，总结一下有限差分法的优缺点。

优点：有限差分法只须构造偏导数的离散方法，这使得它比较容易推广到高阶精度。

缺点，只能在直角网格上进行差分离散，需要物理空间到计算空间的坐标变换，使控制方程变得更为复杂，同时当进行复杂外形流动计算时带来网格划分的不便。

参 考 文 献
[1] 顾尔祚.有限差分法基础.上海：上海交通大学出版社，1988.。