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

− 2ui, j + ui, j−1 (Δy)2
+ O(Δy)2
将上述各有限差分代入式（8－6），原动量微分方程即成如下线性体系：
Ajui, j−1 + B jui, j + C jui, j+1 = Fj
其中，
Aj
=
−
Δx 2Δy
(2vi−1, j
− vi−2, j ) −ν
Δx (Δy)2
（8－7）
同样，对 ui, j+1 、 ui, j−1 用泰勒级数展开，可得：
ui, j+1
=
ui, j
+
Δy(
∂u ∂y
)i,
j
+
(Δy)2 2!
∂ 2u ( ∂y 2 )i, j
+ O(Δy)3
ui, j−1
=
ui, j
−
Δy(
∂u ∂y
)i
,
j
+
(Δy)2 2!
∂2u ( ∂y 2 )i, j
− O(Δy)3
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 ，
(ρ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
⎡ ⎢ ⎢⎣
的切应力。
3
τs
=
μ
(
∂u ∂y
)
s
=
μ 6Δy
(18ui,2
− 9ui,3
+
2ui,4 )
+ O(Δy)3
（8－10）
此法带有一个重要的缺点，即坐标系（笛卡尔直角坐标系或极坐标系）一般不能和边 界层的外形相符（y=δ，而δ随 x 变化）。当边界层增长时，必须添加更多的网络点，或从 一开始就带有额外的网点，这样就使得在下游方向每走一步，都要核对横流向的每一步的
边界层的微分方程比积分方程给出了更为详细的物理描述，但需要把流体所受的剪切应 力和流体的传热量分别与贯穿整个边界层的速度分布和温度分布联系起来，而不像积分方程 那样，只牵涉到壁面的情况。对于层流，情况较为简单；可对于紊流，必需有一些假设来提 供必须的关系式。如常用的假设为有效粘度模型，即
τ
= μ eff
η = 0 和η = 1之间。如图 8－2 所示，其中 rI 为内边界离对称轴的距离， r 为边界内的某 一点离对称轴的距离， φ 为母线与对称轴的夹角。
由于引入了无量纲流函数η ，式（8－13）、（8－14）和（8－15）均可用 (ζ , η ) 坐标
表示成：
∂Φ + (a + bη) ∂Φ = ∂ (c ∂Φ ) + d
∂y
∂x
方程，同时可在边界层动量方程中除去横流向的速度 v ，使之简化。
由于坐标的变换，则式（8－12）中的 ∂u 、 ∂u 等项成为： ∂x ∂y
∂u = ∂u ∂ψ + ∂u ∂ζ = −ρvr ∂u + ∂u
∂x ∂ψ ∂x ∂ζ ∂x
∂ψ ∂ζ
∂u = ∂u ∂ψ + ∂u ∂ζ = ρur ∂u
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
∫ vi, j
=
−
y 0
∂udy ∂x
（8－8）
这个积分可用梯形积分计算。被积函数由一群直线段来近似，每个线段下的面积即为梯形。
故
vi, j
=
vi, j−1
−
Δy 2
⎢⎣⎡(
来自百度文库
∂u ∂x
)i
,
j
−1
+
(
∂u ∂x
)i,
j
⎤ ⎥⎦
（8－9）
其中， vi,1 = 0 。
在给定的 i 处，如 ui,1 及下游站的速度已确定后，可用一个线性插值函数来计算壁面处
∂ζ
∂η ∂η ∂η
此处，
Φ
为一充代变量，可以是
u
、
h~
或
C
* j
，系数分别定义为：
（8－16）
a ≡ rI m& I ψ E −ψ I
b ≡ rEm& E − rI m& I ψ E −ψ I
μeff
c
≡
ρur 2
(ψ
E
σ eff −ψ I )2
表 8－1 给出了和 Φ 相应的σ eff 和称为源项的 d 值。
i−2, j i−1, j i,j
对节点 (i, j) 处的速度 ui, j 和 vi, j 分别用泰勒级数展
开，可得：
ui−1, j
=
ui, j
−
Δx(
∂u ∂x
)i,
j
+
(Δx)2 2!
(
∂2u ∂x 2
)i,
j
− O(Δx)3
i−2, j−1 i−1, j−1 i, j−1
x
图 8－1 层流边界层的有限差分网点
分方程用有限差分近似转换成常微分方程，再用迭代法进行数值积分。还有一种是参数积分 法，将偏微分方程乘以自变量或他变量的权函数后，横越边界层进行积分（对 y 积分），得 出一组常微分方程。这些方程可以用矩阵表示，矩阵的微分系数为未知量，经过矩阵逆变后 用已知量表示，然后再解此方程组。
这里为了便于大家对边界层偏微分方程数值解法的了解，主要介绍最常用的数值计算 方法――斯帕丁－帕坦克(Spalding－Patankar)的方法。许多数值计算的方法现在都已编译成
表 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）的形式表示后，需将转换成有限差分方程，以便求 出数值解。
=
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）等提出了一个将紊流的能量方程和动量方程联立求解的数值解程
第八章 边界层方程的数值解
对于边界层，第五章中给出了其微分方程和积分方程。这两类方程均可用数值解法求解。 在高速数字计算机尚未普遍地得到应用之前，常采用积分方程的数值解。为了更准确地预测 边界层内的计算结果，以及在电子计算机的速度和储存量都大大增加的有利情况下，对于边 界层的微分方程也完全有能力进行数值计算。
序。
另一种称为斯帕丁－帕坦克(Spalding－Patankar)的方法，是采用普朗特混合长度的模
型，即假定
μ eff
= μ + ρlm2
∂u ∂y
（8－4）
其中， lm 是普朗特混合长度。有了这种紊流模型后，便可将边界层的微分方程化成差分方
程，然后进行数值计算。 边界层微分方程的另一种方法是横流向积分，即在一个给定地点（给定的 x），将偏微
μeff Preff
∂h~ ∂ψ
+ (μeff
− μeff Preff
)∂ ∂ψ
(
1 2
u
2
⎤ )⎥ ⎥⎦
⎪⎫ ⎬ ⎪⎭
（8－14）
其中， h~ ≡ CpT
+ 1 u 2 称为全焓 2
， Preff
为有效普朗特数。
经过这样的变换后，边界层的动量方程和能量方程具有相同的形式，即同为抛物型的
偏微分方程，因此能用相同的数值计算程序来求解。事实上，若对质量扩散方程
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
∂u ∂y
（8－1）
其中， μeff 称为有效动力粘度，
μeff = μ + ρεt
（8－2）
层流时，εt = 0 ；紊流时，必须由半经验关系式给出 μeff ，不同的模型给出的半经验关系式
也不同，因而也有不同的数值解法。
布腊德晓（Bradshaw）等提出了将有效粘度与紊流动能 K 相联系的模型：
μeff
1
了标准的计算软件，如 Ansys 软件、Fluent 软件等。 8—1 不可压缩流体层流流动动量方程的有限差分解
不可压缩流体在层流流过平板时的边界层的连续方程和动量方程（忽略体积力）为：
∂u + ∂v = 0 ∂x ∂y
（8－5）
边界条件为：
u
∂u ∂x
+v
∂u ∂y
=
−
1 ρ
dp dx
+ν
∂2u ∂y 2
对 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
∂(ρ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时，上式即简化为笛卡尔坐标系的方程。
∂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 ∂ζ
=
∂ ∂ψ
一、网格的划分
用数值法计算边界层内的动量、热量和质量交换情况时，必须把计算区域划分成网格，
并沿横流向解有限差分方程。网点之间的距离，亦即网格的尺寸应根据边界层内的速度或焓、
ui, j 值的变化，来试出边界层的外缘，因此在计算时间和储存方面都是不经济的。更好的办
法是用一个随着边界层增长而增长的坐标系。
8—2 流线坐标系
对于流体沿着平板的流动，如喷嘴或扩压器内流体的流动、沿对称物体的流动、圆管 和扁管内的流动以及各种自由射流均可建立轴对称极坐标系，使问题简化。
对于稳态、可压缩流体的流动，用轴对称极坐标系时，其连续方程和动量方程（忽略 体积力）为：
在对边界层微分方程进行数值计算时， 为了使边界层的有限差分网格会自动增长 或收敛，与定义的边界层外缘相一致，这里
特定义一个无量纲流函数η ，
等η线
yη ,ψη ,η x,ζ
E边界，η=1 I边界，η=0
η ≡ ψ −ψ I ψ E −ψ I
φ
r r1
对称轴
图 8－2 无量纲的流线坐标系
5
其中，ψ I 和ψ E 分别为所考虑的内边界 I 和外边界 E 处的流函数，通常ψ I 和ψ E 是 ζ 的函 数。这样，顺流向的变量为 ζ ，而横流向的变量则为η ，等η 线与等 ζ 正交。边界层位于
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 为有效斯密特数。
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) 个
采用轴对称坐标，仍避免不了矩形网格的缺点，为此需采用轴对称的流线坐标系。
引进流函数ψ 作为横流向的自变量，顺流向的坐标用 ζ ，ζ ≡ x 。这样，以 ( x , y )
表示的轴对称坐标系就变成了以 (ψ , ζ ) 表示的轴对称流线坐标系了。这种坐标变换称为
冯·米塞兹变换。
由于引入了流函数ψ ，则 ρur = ∂ψ ， ρvr = − ∂ψ 。显然，流函数能自动满足连续