有限体积方法

第三讲 空间离散方法—有限体积法
由于控制方程的复杂性，很难求出其解析解，一般采用数值方法对其进行求解。

采用数值求解方法，首先要对流场空间进行离散，即用一些基本体积单元对物理空间进行填充，要求这些体积单元既不能重叠，也不应有间隙，我们称这些体积单元为网格，或控制体积，填充的过程则称为网格生成。

对于二维流动，基本的网格单元有三角形和四边形网格，而对于三维流动，则基本的网格单元可由四面体、三棱柱、金字塔和六面体单元组成，图3.1即为机翼附近网格。

网格划分完成后，就可以应用相应的数值求解方法把每个网格单元中心点处的流动变量求解出来，也就完成了全部流场的计算。

有限体积法就是针对每个控制体积直接对积分形式的控制方程进行离散，从而把积分型方程近似为代数方程进行求解的方法。

图3.1 机翼附近网格
3.1 N-S 方程的半离散形式
积分形式的N-S 方程为： ∫∫Ω∂Ω
=⋅−+Ω∂∂0)(dS n F F Qd t V c r （3-1） 针对空间某一控制体I Ω，首先对时间导数项进行处理，假设守恒变量Q 在控制体积内为常数分布，即等于控制体中心点处的值I Q （也即为控制体积内守恒变量的平均值），有
∫Ω∂∂Ω=Ω∂∂t Q Qd t I （3-2） 式（3-1）变为 ∫Ω∂⋅−Ω
−=∂∂dS n F F t Q v c I r )(1 （3-3）
假设对流通量和粘性通量在控制体界面上为常值分布，且等于界面中心点（面心）处的值，则有 ⎥⎦
⎤⎢⎣⎡Δ⋅−Ω−=∂∂∑=F N m m m v c I S F F t Q 1)(1 （3-4） 对式（3-3）右端项的近似称为空间离散，而式（3-4）时间方向暂时保留连续的形式，所以称该式为半离散控制方程。

式（3-4）中的m S Δ为第m 个界面的有向面积，即该面的外法线矢量与界面面积的乘积，为一矢量，又称面积矢量。

仔细观察半离散方程可以发现：时间导数项是由单元中心点处的守恒变量值表示的，我们称其为单元中心法；式（3-4）右端项中的通量是关于界面处流动变量的函数，需由界面处的流动变量来确定，由此可看出，流动变量I Q 与流动通量m S F Δ⋅的空间存储位置不同，要想求出流动通量，需先假设流动变量在控制体积内的分布规律，这一过程称为重构，然后确定界面处的流动变量值，再求出界面处的流动通量。

这是有限体积法的主要特点之一。

根据界面处流动变量的确定方法，可将CFD 的空间离散中（主要指无粘通量的离散）的计算方法分为两大类：中心格式和迎风格式。

3.2 中心差分格式
在第m 个界面上的无粘通量为：m m c S Q F Δ⋅)(，即为界面上的守恒变量的函数。

而在中心差分格式中，界面上的守恒变量的确定方法很简单，只需取界面两边单元中心点上的守恒变量的平均值即可。

即该界面上的无粘通量为：m J I c S Q Q F Δ⋅+))(5.0(，其中I Q 和J Q 为与第m 个界面相邻的两边单元中心点上的守恒变量值。

对于粘性通量的计算，需要知道流动参数在界面面心处的空间导数，下面介绍一
种基于节点重构思想的空间导数计算方法。

图3.2 两个相邻网格单元节点图
在图3.2中，是两个相邻的四面体网格单元，网格节点为nl，n2，n3，n4，n5，网格的格心点是cl，c2，要求的是由节点n1，n2，n3所组成的公用面上的空间导数。

设空间导数x q 、y q 、z q 由下列方程组进行求解:
⎥
⎥⎥⎦⎤⎢⎢⎢⎣⎡−+−+−=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡−+−+−+−+−+−+−−−21313212213213213132132132121212)(5.0)(5.0)(5.0)(5.0)(5.0)(5.0)(5.0)(5.0n n n n n n c c z y x n n n n n n n n n n n n n n n n n n c c c c c c q q q q q q q q q q q z z z y y y x x x z z z y y y x x x z z y y x x （3-5） 在上述的方程组中，要使用到各个网格节点处的流动参数，其求解公式如下:
⎟⎟⎠
⎞⎜⎜⎝⎛⎟⎟
⎠⎞⎜⎜⎝⎛=∑∑==N i i N i i i n r r q q 11,01/ （3-6） 其中:
[]2/12,02,02,0)()()(n i n i n i i z z y y x x r −+−+−= （3-7）
i q ,,0是节点n 所在的网格格心处的流动参数，i r 是点n 到网格格心点的距离。

相比较而言，中心差分格式比较简单，程序实现上也没有多少困难。

但是由于中心差分格式是一种中性稳定格式，是非耗散的，耗散项都为零，它无法阻尼掉计算中出现的任何误差，因而在计算中会出现锯齿形的高频振荡，常常导致计算失败。

因此需要引入人工耗散项，其中二阶差分人工粘性项用于抑制解在激波附近的振荡，四阶差分人工粘性项用于阻尼高频误差，并使解趋于定态。

但是，人工粘性项系数的选择，对结果的影响是很大的，特别是在粘性流动计算上，能否选择一个合适的系数直接关系到计算结果的准确性。

因此，更受欢迎的是下面要介绍的迎风格式。

3.3 迎风格式
相对于上面提到的中心格式，迎风格式更受欢迎。

与中心格式相比，迎风格式在边界层内具有更高的计算精度及在大变形网格中适应性强等优点。

目前广泛应用的迎风格式有Roe 发展的通量差分分裂(FDS)格式和vanLeer 提出的矢通量分裂（FVS）格式。

3.3.1 Riemann 问题
所谓Riemann 问题就是求解一维Euler 方程[1]
0)(=+x t U f U （3-8） 在初值
⎩⎨⎧>≤=)0()0()0,(x U x U x U r
l （3-9）
下的解。

假设在一条均匀的直管中间有一个薄片，薄片两边的气体参数是两种不同的常数，如果突然撤去薄片，管内所发展出的流态非常复杂，这类问题在物理上称为激波管问题，它是一种典型的Riemann 问题。

激波管问题如图3.3所示，在0=t 时刻，管内速度为零，隔膜左边为高压，右边为低压，初始的压力间断将以非定常正激波的形式向右传播，同时，一个等熵非定常中心稀疏波向左传播。

图3.4（d），（e），（f）给出了撤去隔膜后1t t =时刻管内的压力、速度和密度的分布情况。

图3.3 激波管问题 图3.4 激波管流动示意图 在上述激波管问题中，为了分析简单明了，仅仅假设隔膜两边的压力不连续。

实际上，同样可以假设速度、密度等不连续，由此就导出了更多、更复杂、更一般化的间断问题，这就是Riemann 问题。

Riemann 问题可通过解析求解的方式求出其精确解。

3.3.2 Godunov 方法简述[1]
回忆一下激波管或Riemann 问题，显然，在每个网格同相邻网格的交界处，如图
3.5的)2/1(−i 和)2/1(+i 处，流动变量是间断的。

也就是说，这些地方都是激波管或Riemann 问题。

因此，可以这样进行计算：在n t =时间步，流动变量为分段常数；然
图3.5 分段界面的Riemann 问题及其分段求解的时间推进示意图
后在每个网格交界处求解Riemann 问题，接下来将得到的Riemann 问题的流动变量解在每个网格内求平均，便得到了1+=n t 时间步的流动变量分布，再对每个网格交界
面处求解Riemann 问题……如此循环，最后收敛就得到了每个网格的流动变量解，也就是全流场的数值解。

这就是Godunov 方法或Godunov 格式的基本原理。

显然，Godunov 方法在求解Riemann 问题时考虑了波的传播特性，所以它是一种迎风方法。

Godunov 方法在求解Riemann 问题时求解的是非线性一维Euler 方程，计算量非常大，从而阻碍这种方法的直接应用。

但Godunov 方法的思路却启发了后来的CFD 研究者，发展了一些Riemann 问题的近似解法以减少计算量、提高计算效率。

这些方法就成为目前CFD 的主要计算方法，如著名的Roe 格式。

3.3.3 通量差分分裂（FDS）格式
3.3.3.1 Roe 的Riemann 问题近似解[2]
考虑标量方程
0=∂∂+∂∂x
f t q （3-10） 将其线性化为
0=∂∂+∂∂x
q t q λ （3-11） 在控制体积内积分该式(化为积分形式) 0=Ω∂∂+Ω∂∂∫∫ΩΩd x f d t q （3-12） 时间导数项变量取单元内平均值i q ，第二项应用高斯定理： ∫∫Ω
∂=+∂∂Ω0dS fn t q x i i （3-13） 该式的半离散形式为（注：网格间距为x Δ，则x i Δ=Ω，界面面积为单位1） )(12/12/1−+−Δ−=∂∂i i i f f x
t q （3-14） 下面讨论如何确定界面处的通量f ，当0)(>′=q f λ时，用迎风差分，即
)(2/1)(2/1112/1i i i i i i f f f f f f −−+==+++ （3-15） 当0<λ时：
)(2/1)(2/11112/1i i i i i i f f f f f f −++==++++ （3-16） 注意到)(11i i i i q q f f −=−++λ，上面两式的统一形式为 )(2/1)(2/1112/1i i i i i q q f f f −−+=+++λ （3-17）
同理，在界面2/1−i 处的通量为 )(2/1)(2/1112/1−−−−−+=i i i i i q q f f f λ （3-18）
推广到方程组（变量为矢量），有
)(2/1)(2/1112/1i i i i i Q Q A F F F −−+=+++
)(2/1)(2/1112/1−−−−−+=i i i i i Q Q A F F F
式中Q F A ∂∂=/，是通量的Jacobian 矩阵，对角化后1−Λ=R R A ，Λ是由A 的特征值组成的对角阵，则A 定义为
1−Λ=R R A 其中Λ是A 的特征值取绝对值后组成的对角矩阵。

上述的计算是以非守恒（拟线性）方程为基础导出的，为了保持格式的守恒性，Roe 提出式（3-33）中的A 的计算必须使
)()()(l r l r Q F Q F Q Q A −=−
由于A 是守恒变量Q 的函数，Roe 的推导证明了：为了得到这样的A ，现用A ~表示，A 的计算状态需取界面左右两侧守恒变量以ρ做加权平均，即
r
l r r l l q q q ρρρρ++=~ 所以，)~(~Q A A =，式（3-33）中的A 需用A ~代替，称该加权平均为Roe 平均。

为了
保持篇幅的完整性，下面不加推导的直接给出在非结构网格中的Roe 通量差分分裂（FDS）格式的具体表达式(为了与下面高阶格式符号相一致，将下标用大写字母表示)。

采用迎风格式在界面k 的无粘通量为
[]
k L R R L k Q Q A Q F Q F F )(~)()(2
1−−+= （3-19） 上式中
541~~~)(~F F F Q Q A L R Δ+Δ+Δ=− （3-20） 其中
⎪⎪⎪⎭
⎪⎪⎪⎬⎫⎪⎪⎪⎩⎪⎪⎪⎨⎧⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡Δ−Δ+Δ+ΔΔ−ΔΔ−ΔΔ−Δ+⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎣⎡++⎟⎠⎞⎜⎝⎛Δ−Δ=ΔU U w w v v u u U n w U n v U n u w v u w v u a p U F z y x ~~~~0~2~~~~~~1~~~22221ρρ （3-21）
⎥
⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡±±±±⎟⎠⎞⎜⎝⎛Δ±Δ±=Δa U h a n w a n v a n u a U a p a U F z y x ~~~~~~~~~1~2~~~~~025,4ρ （3-22） 上面两式中，z y x n w n v n u U ~~~~++=，w n v n u n U z
y x Δ+Δ+Δ=Δ，这里符号Δ的含义为： L R Q Q Q −=Δ （3-23）
前面叙述中出现的下标L l R r ,,,解释如下：有限体积法是对每一个控制体积（单元）将积分形式的控制方程进行离散，该控制体积一般又称为左单元，与其相关的量加下标L 或l 进行表示，而与该单元相邻的控制体积，称其为右单元，与其相关的量加下标R 或r 进行表示。

3.3.3.2 高阶Godumov 格式[3]
上述计算中，假设n 时刻每个单元的守恒变量是常值分布，因此是一阶精度Godumov 格式。

为了推广到更高阶情形，用单元格心处状态重新构造每一单元中的守恒变量分布，以一维问题为例，例如假设线性分布。

也即在每个单元中假设
2/12/1),
()(+−≤≤−+=i i i i i x x x x x S Q x Q （3-24） 式中i S 是斜率，例如可选取 1
1−−−−=i i i i i x x Q Q S （3-25） 因此可找出交接面2/1+i x 的左侧状态L Q ，同样的，对1+i 单元找出2/1+i x 的右侧状态R Q 。

以此两侧状态作为Riemann 问题的初值，则对应的Godumov 格式中的通量可写为 [])(~2
12/1L R L R i Q Q A F F F −−+=+ （3-26） 推广到三维非结构网格，对应一维的斜率应由单元中心点处变量的梯度代替，则界面左侧的变量为：
⎩⎨⎧⋅∇+=⋅∇+=R J J R
L I I L r Q Q Q r Q Q Q r r （3-27） 上式中，I Q ∇为在单元I 中心点处守恒变量的梯度，L r r 为单元I 中心点到界面面心的
矢量。

对J 单元含义相同，不在重复。

单元I 中心点处守恒变量的梯度可由高斯定理计算得到
∑=Δ+Ω≈∇F N J ij ij J I I
I S n Q Q Q 1)(211v （3-28）
我们称这一重构方式为线性重构。

3.3.3.3 限制器[2] 可以证明，用式（3-24）、（3-25）的重构形式得到的格式具有二阶精度，但稳定性分析理论指出，如此构造的显式格式为无条件不稳定。

问题是经此重构后通过外插值得到的界面两侧的值可能超出原来的单元内常值分布，即可能产生新的极值。

为此，必须对式（3-24）中的斜率加以限制。

例如经过限制后的斜率可取为：
1
1)
(−−−−=i i i i i i x x Q Q r S φ （3-29） 上式中的i r 为 ⎟⎟⎠⎞⎜⎜⎝⎛−−⎟⎟
⎠⎞⎜⎜⎝⎛−−=−−++1111/i i
i i i i i i i x x Q Q x x Q Q r （3-30） )(i r φ称为限制器，这里取mod min 限制器，有
)),1min(,0max()(r r i =φ （3-31）
图3.6 限制器原理图
图3.6说明了三种典型情况，图中空心圆圈代表格心处的变量值，实心圆圈代表
加限制器后交界面处的取值。

推广到三维非结构网格，在利用式（3-27）求界面上的变量值时，应通过限制器对式中的梯度加以限制。

除了上面介绍的mod min 限制器外，还有C Osher −、VanLeer 限制器等，不做另外介绍。

至此，我们通过上面介绍的中心差分格式或迎风格式对半离散方程（3-4）的右端的通量项进行了求解，该项一般称为残值或残差，用符号i R 表示，则半离散方程（3-4）化为
Ω
−=∂∂I I R t Q （3-32） 上式右端为已知，所以该方程为关于时间的常微分方程，可以应用成熟的常微分方程求解方法对其进行求解，与其相关的内容将在下一讲讲述。

参 考 文 献
[1] 阎超.计算流体力学方法及应用[M].北京：北京航空航天大学出版社，2006.
[2] 黄明恪.计算空气动力学差分方法讲稿.南京：南京航空航天大学，1991.
[3] Blazek J. Computational Fluid Dynamics: Principles and Applications.。