单边突扩通道中kε方法数值模拟^>[1]

一、引言

如图1所示的单边突扩通道形式，是在燃烧室设计、电厂冷却水研究等工程问题上经常遇到的流动现象。由于流动具有分离，湍流结构复杂，不宜采用传统的近区积分模型。因为这种积分模型需要假定射流边缘上的掺入率及流速、温度的相似剖面，经验性较大；在确切判断半经验公式的系数上，遇到不少困难。近年来，人们开始应用湍流模型来处理这类问题，它可以算得所论区域的湍流场、掺入率及速度、温度剖面而无需经验系数的事先假定。这种方法已由单方程模型发展到双方程模型、多方程模型，是国际上迅速发展的一种新的理论和方法。在上述模型中，k-ε双方程湍流模型应用最为广泛。对浮射流、淹没射流、尾流、回流、次生流等问题，无论定常或非定常、无论二维或三维，都已有若干成功的算例。结果令人鼓舞，预示了这类模型强大的生命力。

二、数学模型

不可压缩流体作平面二维定常运动，当忽略质量力的影响时，其控制方程可以表示为：

称为有效黏性系数，由分子黏性系数及湍流黏性系数所组成。其中k、ε分别称为湍流动能及湍流动能耗散率，它们由以下两输运方程确定：

为湍流动能的产生项。

不考虑散热损失，能量方程可以表示成一无源的对流扩散方程：

方程（1）～（8）组成一封闭系统。若要考虑密度ρ随时间T的变化，例如温差异重流问题，则还应加入状态方程ρ=ρ（T）联立求解。

数学模型中的六个通用常数，系由基本实验确定，根据文献［2］，它们的值见表1。

三、数值计算

差分方程的推导及“壁函数”的应用在文献［3］中已有详细叙述，本文拟就P-V迭代计算中出现的速度修正方程及压力修正方程作一推导，并对CHAMPION程序的初值设定问题提出新的方法。

1.P-V迭代

本文采用斯波尔丁学派创立的SIMPLE算法，以压力P和速度V为基本变量进行P-V迭代计算。即先假定全场压力分布，使运动方程能独立求解。这时解出的速度值自然不满足连续方程，应通过速度修正方程来修正，将修正后的速度值代入压力修正方程又给出新的压力值。如此往复迭代，即能期望得到真实的解。

2.速度修正方程

参见图2，U分量差分方程为：

式中：b为除压力梯度外的源项。

V分量差分方程为：

式（9）与式（10）的A_i是分别对U、V单元而言的。

以初设压力P^*代入上两式得：

令

式（9）-式（9′）及式（10）-式（10′）得：

为了简化计算，通常略去（参见文献［4］），故得：

式中

3.压力修正方程

连续方程（1）的差分形式为：

参见图2，知：

将式（11）、式（12）、式（14）代入式（13），整理即得关于压力修正P′的方程：

其中

式（15）即为所求的压力修正方程，它具有统一形式，可以与其他方程一样求解。

式（16）的B是压力修正方程的源项，它具有连续方程的形式。若B=0，说明U^*、V^*满足连续方程，显然它们又满足运动方程，因此必然就是精确解。通常称B为质量源，以全场各点B值的绝对值之和趋于零作为迭代收敛的判据。

4.边界条件

参见图1，各变量的边界条件可表述为：

AB面：

k=k₀=1.5×（0.05×U₀）²

DE面：

其他各面：

U=V=0

T=T_W

式中：T_W为壁面温度。至于k、ε两个变量，则采用“壁函数”方法处理>^［3］。

5.迭代方法及初值设定

迭代方法一般有点迭代与线迭代两种。点迭代计算较稳定，但太费机时；本文采用线迭代，在每一条线上使方程的系数矩阵成为三对角矩阵形式，可用通用追赶法程序求解。线迭代不如点迭代稳定，为此，应用低松弛法提高稳定性和收敛性。

变量全场初值是否合理，直接影响到计算收敛问题。因为完全不合理的初值，可能一开始就使问题超出收敛半径的范围之外，计算无法进行。对初值设定问题，通常的办法是用粗网格或令交换系数为常数进行粗略计算，将所得结果作为初值。这种做法，不仅同样存在不收敛的危险，同时还使数据的输入输出变得很复杂，增加许多工作量。本文在迭代计算的前若干步，控制变量不超过其可能最大值，保证计算可以进行。经过这若干步的计算，变量初值很快趋于合理。此时取消其控制，任其自由收敛。此法对网格过分不均匀或算法上有本质错误的问题不起任何作用，即使每一步都加以控制，其余源及质量源也永远不会变小，不可能得到收敛解。

收敛判据以相对质量源的绝对值小于10^-2为准，在M-160机上，需CPU时间约6min。

四、计算结果及其分析

本文对标准CHAMPION程序进行改进，将原来的单向扫描改为来回双向扫描。这对回流问题，因为反向扫描可以迅速地将下游影响传递到区域各点，大大提高了计算的稳定性与收敛性。另外，去掉原来的逐线总体连续修正，仅保留出口处一条线不变，计算反而容易收敛。

定义进口雷诺数Re为：

式中：U_m为进口平均流速；ν为水的运动黏性系数；H为入口高度。

为了便于与计算结果^>[2]>[3]及实验资料>^［5］对比，分别取Re=5322、3441、1803进行计算。图3～图5是各个工况下本文的计算结果、文献［1］的计算结果及文献［5］的实验结果。图中各量系下列无量纲数：

1.流线图及时均速度分布

由图3（a）、图4（a）、图5（a）可知，本文计算所得回流区长度在台阶下游7h～9h之间，详见表2。

由表2可以看出，本文计算结果稍偏大。这一方面是由于数学模型的常数c₁对L／h值有一定影响，另一方面是实验资料的精度亦不太高。文献［5］的L／h值是由10个断面的X方向时均速度分布内插而得；文献［2］的结果是由雷诺数Re为5900和8340的L／h值线性外推得来。但事实上，当Re高达一定数值时，L/h值与Re无关，线性外推不合理。因此，文献［2］中，Re=5322时的L/h值显然应比7.0为大。

在图3（b）、图4（b）、图5（b）中，分别给出了所取三个雷诺数情况下，十个横断面上的时均速度分布及文献［1］、［5］的结果。由图可知，除进口附近由于存在三维分离区而使回流区高度较实验值小外，其他的地方均符合良好，并且比文献［1］的计算结果更符合实验资料，特别是在上壁面附近流速最大值处更为明显。

2.湍流动能k、湍流动能耗散率ε及涡黏性系数μ_t分布

图3（c）、图4（c）、图5（c）分别给出了三个雷诺数下的k剖面曲线，图3（d）、图4（d）、图5（d）分别给出了ε剖面曲线，图3（e）给出了Re=5322时的μ_t剖面曲线。剖面图中同时还给出了文献［1］的计算结果。由图可见，吻合情况良好，只是ε在壁面附近符合稍差，这显然是对ε的壁函数的不同处理而引起的。

3.相对压力及温度分布

图5（e）给出了Re=1803时的相对压力剖面曲线，图4（e）给出了Re=3441时的等温线。这里所指的压力，是相对于进口处的压力。由图可见，相对压力是在回流区最小。雷诺数越大（亦即进口速度越大），相对压力越易出现负值，即越易引起通道内部出现低压区，这与高速水流容易引起气蚀的现象是吻合的。对温度的计算，由于本文所选通道尺寸为长×高=20×2cm²，区域范围很小，若按一般绝热要求取作为壁面条件，则计算的温度分布几乎维护进口温度不变。作为数学模型及程序的验证，本文取T_W=0（即壁面处温度等于自然水温）作为边界条件来计算温度。由图4（e）可见，温度亦是在回流区最小。

Numerical Simulation of a Single Side Sudden Enlargement Channel with a k-ε Turbulence Model

Abstract

参考文献

>［1］John F.Ripken，Norio Hayakawa.Cavitation in High-head Conduit Control Dissipators.ASCE，HY₁，1972：239-256.

［2］Rodi，W..Turbulent Models and Their Application in Hydraulics.A state of the Art Review，SFB 80/T/127，May，1978.

>［3］王能家.河道排取水问题的k-ε方法数值模拟∥水利水电科学研究院论文集（第26集）.北京：水利电力出版社，1986.

［4］Suhas V.Patankar.Nemerical Heat Transfer and Fluid Flow，McGraw Hill，New York，1980.

>［5］沈熊，廖湖生.应用频移激光测速系统测量高湍流度回流流动.力学学报，1982（5）：505-511.

---