基于有限体积法的控制方程离散

Chapter 10 基于有限体积法的控制方程离散

10.1 控制方程

稳态气相通用控制方程如下：

()

()()+div ρU div Γgrad S t ρφφφ∂∂=+

()rv u w

U U z r r r θU=ui+vj+wk

div ∂∂∂=∇⋅=++

∂∂∂

()()()() z r r u r v w t z r r r r S

z r r r ∂ρφ∂∂∂

ρφρφρφ∂∂∂∂θ

∂∂φ∂∂φ∂∂φ∂∂∂∂∂θ∂θ+++=Γ⎛⎫⎛⎫⎛⎫

Γ+Γ++ ⎪ ⎪ ⎪⎝⎭⎝⎭⎝⎭

上式中φ表示气相场通用变量，在三维柱坐标下可以是轴向速度u 、径向速度v 、

周向速度w 以及焓h ，φ为1时代表连续方程，Γ表示输运系数，S 为气相源项，U

为速度矢量。密度由完全气体状态方程ρRT p =确定，温度由焓温关系式(h)T f =确定。每个方程中的φ、Γ、S 如表10.1所示。

表10.1 气相方程组

§10.2 方程的离散化

§10.2.1 有限体积法

有限差分法：原理比较简单，数学推演和程序编制的工作量较小，但其灵活性较差；

有限单元法：数学推演和程序编制工作要比有限差分法复杂，因而比较费时，但其灵活性较好；

有限体积法：作为有限差分法和有限单元法的中间产物，继承了两者的优点，克服了各自的缺点，因而在流场计算中占有重要的地位。

有限体积法的基本思路是：将计算区域划分为一系列不重复的控制体积，并使每个网格点周围有一个控制体积；将待解微分方程对每一控制体积积分，得出离散方程，其中的未知数是网格点上的因变量φ。有限体积法在对控制体积积分时，必须假定φ值在网格点之间的分布。

由有限体积法导出的离散方程具有明确的物理意义，即因变量φ在有限大小的控制体积内的守恒原理，如同推导微分方程时因变量在无限小的控制体积中的守恒原理。有限体积法即使在粗网格条件下，也精确地满足

积分守恒。

图10-1 三维控制体积

()()()

[]

z r r z z r z r ,,,,,z z

r r

,,,,,i+j+k r z r r r θ

i+j+k S i+S j+S k S +S +S nb w e s n b t

nb

nb w e s n b t

nb

div div dV ndS

θθ

θθθθ==Φ=ΦΦΦ∂Φ∂Φ∂ΦΦ=++

∂∂∂Φ=Φ⋅⎡⎤=ΦΦΦ⋅⎣⎦=

ΦΦΦ⎰⎰∑

∑

()()

[][]()[][]z

z z r ,,,,,z

z z ,,,,,z

z z ,,,,,r r r

r r ,,,,,,,,,, i+0j+0k S i+S j+S k z S z

S z r S r r r r S r θnb w e s n b t nb

nb w e s n b t

nb

nb w e s n b t

nb

nb w e s n b t nb nb w e s n b t

nb

dV V V V V θθ

θθ=====∂Φ⎡⎤

=Φ⋅⎣⎦∂∂Φ=

Φ=

∆∂∂Φ=Φ∂∂Φ∂ΦΦ=Φ=+∂∂∂Φ=Φ∂∑⎰∑

∑∑∑

有限体积法的四条基本原则是： 原则1：控制体积交界面的一致性。

当一个表面为相邻的两个控制体积所共有时，在这两个控制体积的离

散方程中，通过交界面的通量表达式必须相同。 原则2：正系数。

中心节点系数a P 和相邻节点系数a nb 必须恒为正值。 原则3：源项的负坡线性化。

源项S 常依赖于因变量φ，为便于计算，通常将源项线性化，即将其表达为S S S C P =+φ，此时要求S P 必须小于或等于零。从数值计算的角度看来，可避免出现计算不稳定或结果不合理。

如径向动量方程源项中含有v 2μ

2μv, 0, r

r

C P C P S S S S --=+==

周向动量方程源项中含有ρvw w, 0, ρv C P C P S S S S -=+==- 原则4：相邻节点系数和。

当S P 为零时，中心节点系数a P 应等于各相邻节点系数之和，即满足

a a P n

b =∑。

§10.2.2方程的离散化 将通用控制方程改写为：

()∂ρφ∂∂∂ρφ∂φ∂∂∂ρφ∂φ∂∂∂θρφ∂φ∂θ t z u z r r r v r r w r S +-⎛⎝ ⎫⎭⎪+-⎛⎝ ⎫⎭⎪⎡⎣

⎢⎤

⎦⎥

+

-⎛⎝ ⎫⎭

⎪=ΓΓΓ11

（1） 首先按常用方法将源项作线性化处理，即将源项表达为： S S S c p p =+φ （2） 其次对非定常项

()

∂ρφ∂ t

进行离散化，令 ()∂ρφ∂ρφρφ t t

p p p p

=-00

∆ （3） 此处采用全隐式，假定“新”数值ρφP P ,控制着整个控制体积。上标0

表示时间步长开始时的“老”数值，不带上标的量均表示时间步长结束时的“新”数值。

接下来将方程（1）对图10-1所示控制体积积分，得：

()

ρφρφφp p p p

e w n s t b c p p t V J J J J J J S S V -⎛⎝

⎫⎭

⎪+-+-+-=+00

∆ （4） 上式中V 是控制体积的容积，可表达为

()V r r z r n s =+05. ∆∆∆θ （5）

r r n s ,代表控制体积的上、下半径，∆z 、∆r 、∆θ分别为控制体积在三个坐

标方向的长度。

我们以e ,w ,n ,s ,t ,b 来代表控制体积的6个交界面，其位置如图10-1所示。则：

J u z A J u z A J v A J v A J w r A J w r A w w w e e

e n n n s

s

s t t t b b

b =-⎛⎝ ⎫⎭⎪=-⎛

⎝ ⎫⎭⎪=-⎛⎝

⎫⎭⎪=-⎛

⎝ ⎫⎭⎪=-⎛⎝ ⎫⎭⎪=-⎛

⎝ ⎫⎭⎪ρφ∂φ∂ρφ∂φ∂ρφ∂φ∂ρφ∂φ∂ρφ∂φ∂θρφ∂φ∂θΓΓΓΓΓΓ

r r

（6）

式中A w 、A e 、A n 、A s 、A t 、A b 分别为控制体积的左、右、上、下、前、后表面积，且有

()A A r r A A r z A r z A r z w e s t b n n s s =+=== =0.5r = ;

, ;

n ∆∆∆∆∆∆∆∆θθθ;

（7）

将连续方程

()()()∂ρ∂∂ρ∂∂ρ∂∂ρ∂θ t z r +++=u v w r 0 对控制体积积分，得

()

ρρp

p e w n s t b V t

F F F F F F -+-+-+-=00∆ （8）