第八章 边界层方程的数值解

u
∂C
* j
∂x
+
v
∂C
* j
∂y
=
D
∂
2C
* j
∂y 2
（ C 为质量浓度， D 为质扩散系数）。进行变换，也是变成抛物
型的偏微分方程，也能用该数值计算程序来求解，如下式：
∂C
* j
=
∂
(ρur 2 μeff
∂C
* j
+
Rj
)
∂ζ ∂ψ
Seeff ∂ψ ρu
（8－15）
其中， R j 为由于化学反应等原因组份 j 的产生率， Seeff 为有效斯密特数。
ui−2, j
= ui, j
−
2Δx(
∂u ∂x
)i
,
j
+
(2Δx)2 2!
(
∂2u ∂x 2
)i
,
j
− O(Δx)3
由此可得， ui, j = 2ui−1, j − ui−2, j + O(Δx)2 类 似 地 ， 对 于 vi, j ， 可 得
vi, j = 2vi−1, j − vi−2, j + O(Δx)2
(ρur 2μeff
∂u ∂ψ
)−
1 ρu
dp dζ
（8－13）
这就是以轴对称流线坐标系表示的动量方程。
同样，可将能量方程 ρc p (u
∂t ∂x
+ v ∂t ) ∂y
=
λ
∂2t ∂y 2
+
μ( ∂u )2 也改写成轴对称流线坐标系 ∂y
中的形式：
∂h~ ∂ζ
=∂ ∂ψ
⎪⎧ ⎨ ⎪⎩
ρur
2
⎡ ⎢ ⎢⎣
在对边界层微分方程进行数值计算时， 为了使边界层的有限差分网格会自动增长 或收敛，与定义的边界层外缘相一致，这里
特定义一个无量纲流函数η ，
等η线
yη ,ψη ,η x,ζ
E边界，η=1 I边界，η=0
η ≡ ψ −ψ I ψ E −ψ I
φ
r r1
对称轴
图 8－2 无量纲的流线坐标系
5
其中，ψ I 和ψ E 分别为所考虑的内边界 I 和外边界 E 处的流函数，通常ψ I 和ψ E 是 ζ 的函 数。这样，顺流向的变量为 ζ ，而横流向的变量则为η ，等η 线与等 ζ 正交。边界层位于
y = 0，u = 0，v = 0
（8－6）
y = δ ， u = u∞ x 和 y 分别是顺流向和横流向的坐标，u 和 v 分别是顺流向和横流向的分速度。如用矩
形网格，则层流边界层的有限差分网点如图 8－1 所示。其中， Δx 和 Δy 分别是顺流向和横
流向的步长，均取为常数。
y
i−2, j+1 i−1, j+1 i, j+1
=
lk
ρK
1 2
（8－3）
其中，lk
是距壁面距离的函数，K
是流体脉动的平均动能
(K
=
1 2
ρu′2
+
1 2
ρv′2
+
1 2
ρw′2 )
。
由于式（8－3）中包含有紊流动能 K ，因此必须联立求解紊流的能量方程和动量方程，因
此，布腊德晓（Bradshaw）等提出了一个将紊流的能量方程和动量方程联Leabharlann 求解的数值解程一、网格的划分
用数值法计算边界层内的动量、热量和质量交换情况时，必须把计算区域划分成网格，
并沿横流向解有限差分方程。网点之间的距离，亦即网格的尺寸应根据边界层内的速度或焓、
对 y 的导数用中心差分公式，对 x 的导数用向后差分公式，则得：
2
∂u ( ∂x )i, j
=
3ui, j
− 4ui−1, j 2Δx
+ ui−2, j
− O(Δx)2
∂u ( ∂y )i, j
=
ui, j+1 − ui, j−1 2Δy
+ O(Δy)2
及
(
∂2u ∂y 2
)i,
j
=
ui, j+1
μeff Preff
∂h~ ∂ψ
+ (μeff
− μeff Preff
)∂ ∂ψ
(
1 2
u
2
⎤ )⎥ ⎥⎦
⎪⎫ ⎬ ⎪⎭
（8－14）
其中， h~ ≡ CpT
+ 1 u 2 称为全焓 2
， Preff
为有效普朗特数。
经过这样的变换后，边界层的动量方程和能量方程具有相同的形式，即同为抛物型的
偏微分方程，因此能用相同的数值计算程序来求解。事实上，若对质量扩散方程
采用轴对称坐标，仍避免不了矩形网格的缺点，为此需采用轴对称的流线坐标系。
引进流函数ψ 作为横流向的自变量，顺流向的坐标用 ζ ，ζ ≡ x 。这样，以 ( x , y )
表示的轴对称坐标系就变成了以 (ψ , ζ ) 表示的轴对称流线坐标系了。这种坐标变换称为
冯·米塞兹变换。
由于引入了流函数ψ ，则 ρur = ∂ψ ， ρvr = − ∂ψ 。显然，流函数能自动满足连续
∂y
∂x
方程，同时可在边界层动量方程中除去横流向的速度 v ，使之简化。
由于坐标的变换，则式（8－12）中的 ∂u 、 ∂u 等项成为： ∂x ∂y
∂u = ∂u ∂ψ + ∂u ∂ζ = −ρvr ∂u + ∂u
∂x ∂ψ ∂x ∂ζ ∂x
∂ψ ∂ζ
∂u = ∂u ∂ψ + ∂u ∂ζ = ρur ∂u
∂u ∂y
（8－1）
其中， μeff 称为有效动力粘度，
μeff = μ + ρεt
（8－2）
层流时，εt = 0 ；紊流时，必须由半经验关系式给出 μeff ，不同的模型给出的半经验关系式
也不同，因而也有不同的数值解法。
布腊德晓（Bradshaw）等提出了将有效粘度与紊流动能 K 相联系的模型：
μeff
ui, j 值的变化，来试出边界层的外缘，因此在计算时间和储存方面都是不经济的。更好的办
法是用一个随着边界层增长而增长的坐标系。
8—2 流线坐标系
对于流体沿着平板的流动，如喷嘴或扩压器内流体的流动、沿对称物体的流动、圆管 和扁管内的流动以及各种自由射流均可建立轴对称极坐标系，使问题简化。
对于稳态、可压缩流体的流动，用轴对称极坐标系时，其连续方程和动量方程（忽略 体积力）为：
ui, j ( j = 2,3,L, N −1) 的线性代数方程组。因系数 Aj 、 B j 、 C j 、 Fj 中只包含已知量（即
预先算出的或初始的
u
和
v
值，以及已规定的
dp dx
值），ui,
j
的各未知值即可由形成的矩阵式
计算得出，边界条件为 ui,1 = 0 ， ui,N = u∞ 。得出 ui, j 后，可由连续方程求 vi, j ，
∂y ∂ψ ∂y ∂ζ ∂y
∂ψ
4
dp = dp （忽略横流向压力变化，即忽略 dp ）
dx dζ
dψ
1 r
∂(rτ ∂y
)
=
ρu
∂(rτ ∂ψ
)
经过转换，式（8－12）就成为：
∂u = ∂(rτ ) − 1 dp ∂ζ ∂ψ ρu dζ
当把τ
=
μ eff
∂u 代入后，即得： ∂y
∂u ∂ζ
=
∂ ∂ψ
∂(ρur) + ∂(ρvr) = 0
∂x
∂y
（8－11）
ρ(u ∂u + v ∂u ) = − dp + 1 ∂(rτ ) ∂x ∂y dx r ∂y
（8－12）
其中， x 和 y 仍分别是顺流向和横流向的距离， u 和 v 仍分别是顺流向和横流向的分速度。 r 为离对称轴的半径， r = 1时，上式即简化为笛卡尔坐标系的方程。
第八章 边界层方程的数值解
对于边界层，第五章中给出了其微分方程和积分方程。这两类方程均可用数值解法求解。 在高速数字计算机尚未普遍地得到应用之前，常采用积分方程的数值解。为了更准确地预测 边界层内的计算结果，以及在电子计算机的速度和储存量都大大增加的有利情况下，对于边 界层的微分方程也完全有能力进行数值计算。
表 8－1 和 Φ 相应的σ eff 和称为源项的 d 值
他变量 Φ
σ eff
d
u
1
1− 1 dp ρu dζ
h~
Preff
∂ ∂η
⎡ ρur 2
⎢ ⎢⎣
(ψ
E
−ψ
I
)2
(μeff
−
μeff Preff
)∂ ∂η
(
1
u
2
⎤ )⎥
2 ⎥⎦
C
* j
Seeff
Rj ρu
8—3 边界层方程的有限差分形式
6
将基本的偏微分方程用式（8－16）的形式表示后，需将转换成有限差分方程，以便求 出数值解。
Bj
=
3 2
(2ui−1,
j
− ui−2, j ) + 2ν
Δx (Δy)2
Cj
=
−
Δx 2Δy
(2vi−1, j
− vi−2, j ) −ν
Δx (Δy)2
Fj
=
1 2
(2ui−1, j
− ui−2, j )(4ui−1, j
− ui−2, j ) −
Δx ρ
dp ( dx )i
对于横贯边界层的 N 个网络点（y 方向除 y=0 和 y=δ以外）而言，就有 (N − 2) 个