浮力回流

方法数值模拟^>[1]

一、引言

火力发电厂的冷却水取水水温，是影响电厂出力的重要因素；受纳水域的水温升高，会影响到水的环境。核电站冷却水有时还带有效射性物质，受纳水域除热污染外，还存在放射性物质污染问题。因此，对水电、核电建设来说，预报水域水体的温度分布与水体异物的浓度分布是十分必要的。

随着计算技术的发展，数学模拟方法日益被人们所重视，它与物理模拟一样，是环境预测技术必不可少的工具。数学模型一般分为整区模型与区域模型。目前普遍采用后者，即将受纳水域按流动的特性划分为近区与远区分别加以模拟，在每一区内按各自的特点简化控制方程。针对这两个不同区域，已有了不少数学模型>^［1］。由于近区流动的复杂性，近区模型不及远区模型成熟，不少问题尚有待解决。近区不但是污染最重的地区，直接反映着热水及其他污源排放的物理性质，是研究工程及环境的主要关注所在，而且远区模型也以近区计算成果作为其边界条件。k-ε模型用于近区相当成功>^［2］，本文的突扩部分计算结果，也证实了这一点。

区域模型虽然有利于数学处理，经济易行，但区域划分还不明确，成果衔接亦存在不少问题。因此，整区模型仍应该是发展方向。k-ε方程为建立整区模型创造了条件。

在濒临江河地区，一般是直接引取河水作为冷却水水源。图1是这种直流取水或混合取水的工程布置示意图。

这是一个典型的整区问题。20世纪60年代初，陈惠泉、许玉林对此作了势流分析、电模拟及水槽试验验证>^［3］；在数学模拟方面，麦国克作过尝试>^［2］，但未见正式发表。本文针对这一问题，用k-ε模型进行模拟，并与水槽试验成果进行对比。结果基本合理，总的趋势与试验结果定性符合，证实了这类精细模型的可用性与合理性。在计算排取水问题前，为了验证数学模型及数值方法的正确性，参照文献［4］，先用此法计算单边突扩问题，所得结果与文献［4］的试验结果符合良好，比用φ-ω方法的计算结果稍优^>[2]，详见成果分析部分。

二、数学模型

考察浅水定常出流，棱柱体渠道，选用水深平均形式的方程。流体不可压缩，不计流散、浮力及表面风应力的影响，自由面作刚盖近似，数学模型具有下述形式：

其中

以上各式中：x、y为坐标变量；U、V为相应于x、y方向的水深平均时均速度；T为水深平均时均温升（实际水温与平衡水温之差）；k、ε为水深平均形式的湍流动能及湍流动能耗散率，；c_f为经验摩擦系数，取为0.003分别为涡黏性系数与涡扩散系数；c₁、c₂、c_μ、σ_k、σ_ε、σ_T为通用常数，它们皆由基本试验确定，不依赖于特定的问题。参见文献［2］，分别取值见表1。

本模型中的k、ε两个输运方程，是从纳维—斯托克斯方程出发推导而来。解这两个方程以确定全场各点的涡黏性系数与涡扩散系数，从而充分反映各点湍流紊动的差异对动量扩散与热量（或质量）扩散的不同影响。这就是k-ε模型区别于一般常系数扩散模型的本质所在。

三、数值计算方法及程序说明

（一）统一形式的微分方程

控制方程中的每个方程都具有类似的项，为便于差分，可以表示成如下统一的形式：

或

式中：φ为任一因变量，如U、V、T、k、ε等；Γ_φ为φ变量对应的交换系数。

式（7）中〈1〉项为时间变率项，虽然考虑定常问题，但仍保留此项，为便于今后应用时间相关法；〈2〉项为对流项；〈3〉项为扩散项；〈4〉项为源项，但并非真实的物理源，而是将不能纳入前三项的一切项都归并于源项。

为清楚起见，表2列出了对应不同的φ方程式（7）中各项的具体内容。

表2中，μ_eff=ρυ_eff=ρ（ν+ν_t）=μ+μ_t，称为有效黏性系数，它是分子黏性系数与湍流黏性系数之和。

基本方程组是一多元非线性方程组，各方程间皆具有耦合性，因此只能用数值方法求解，且必须使用迭代方法。

（二）差分方程及其边界条件

1.网格系

建立一个差分方程，首先应将求解区域离散化，形成差分网格。本文采用交错网格如图2所示。

为了提高效率，网格在x、y方向皆取成非均匀，节点量随L、D、B₁、B₂有所改变，一般为10×60或12×54个。排水及取水口处应密分，一般为3～5个节点。图2中，X方向实线为U网格线，Y方向实线为V网格线，虚线及边界实线为G变量（K、ε、T、h）网格线。

2.差分方程

采用五点差分格式，把变量φ在P点上的值φ_P表示成邻节点上φ值的函数：

式中，E、W、N、S分别代表P点的东、西、北、南四点，见图3。

使用控制容积法，将统一形式的微分方程在某一单元P上积分，可得

式中：为单元体体积、单元表面积矢量；为对全部单元边界面求和（i=E，W，N，S）。

式（9）的〈1〉项为时间项，可取向后差分：

式中：P为括号内的量取t时刻的值；P^-为括号内的量取t-Δt时刻的值。

当网格不随时间而变时，不变：

式（9）的〈3〉项为源项，要线性化之后才能求解，一般将S_φ线性表示成：

式（9）的〈2〉项为总通量项，根据φ分布函数的不同取法，一般有中心格式、迎风格式、混合格式、乘方定律格式、指数格式等5种常用差分格式。本文的突扩部分计算采用迎风格式，排取水问题计算采用乘方定律格式。为简单起见，以迎风格式为例，推导式（8）中的A系数形式。

参考图3，式（9）的第〈2〉项应为：

对于迎风格式：

w面项为：

［［］］表示取两者之大。

同理，e面项为：

s面项为：

n面项为：

由上述各式，〈2〉项变为：

式（16）最后一项中，-F_w+F_e-F_s+F_n正是连续方程左边项的差分应为零，将式（10）、式（11）及式（16）代入式（9）得：

其中

式（17）即为所求的差分方程。

3.边界条件

AG：V=0或=0，T=0或=0，U、k、ε给定

FH：V=0或=0，=0（φ=U、k、ε、T）

BC：V、k、ε给定，U=0=0

DE：V、k、ε、T给定，U=0

全场水深在平均水深h₀上、下波动，但变幅Δh极小，一般｜Δh｜/h₀&lt；6%。因此，边界处的水深可以通过附近的计算值外推而得。在排、取水口及上游来流处，根据流量、平均水深及给定的速度分布型式，即可确定U₀或V₀。k₀取为×（0.05U₀）²或×（0.05V₀）²，ε₀取为，式中，c_μ如前所述，取为0.09，L为湍流长度尺度，可取各自的特征长度，如水深、水力半径等的0.03～0.05倍。

AB，CD，EF，GH：U=V=0，=0，对k、ε应作特殊处理，因这里碰到的是固体壁面问题，不能简单给定k=0，ε=0，否则将会使全场的k、ε值算小。本文采用常用的“壁函数”处理方法。设距壁面δ处一点P，平行于壁面的速度分量可用下述对数分布公式计算

且假定在壁面附近，湍流处于局部平衡状态，即认为湍流的对流与扩散可以忽略，产生与耗散相平衡。于是，p点的湍能及湍能耗散率可用下列式子计算：

ε_p=U^*3/（kδ）

式中，k=0.41，E=9.0，U^*=（τ_w/ρ）^1/2。

（三）差分方程的求解方法

如前所述，基本方程组是一相互耦合的非线性方程组。因此，求解工作应包括：①系数线性化，即在求解之前，用各变量的当前值计算出所有系数；②解除各方程之间的耦合，使每个方程都能独立求解。但是以原始变量U、V、P描述的方程式，U、V与P的耦合不易解除，本文采用h—V（或P—V）修正方法来处理这一困难。即先假定全场的水深分布，使动量方程能够独立求解。这时所得的值不是真实的，应通过连续方程来修正，从而给出新的水深分布。如此往复迭代，即能期望得到真实的解。

针对本问题的复杂性，采用线迭代四向扫描的方式进行计算，以提高稳定性，改善收敛性。

（四）程序说明

本文计算程序是在帕坦卡（Patankar）的标准层流程序的基础上修改而成。修改内容包括：①几何域、边界条件、网格划分；②改压力p为水深h；③增加湍流量k、ε，温度T的计算以及它们的壁函数处理；④增加若干子程序。

为方便起见，本程序定名为四向扫描程序。

四、计算成果及其分析

1.突扩通道的计算

为了验证所用数学模型的正确性及数值方法、计算程序的合理性与可用性，在处理河道排取水问题之前，先对如图5所示的二维单边突扩问题进行了计算。

对图5所示的通道形式，本文仅给出一种工况下的断面流速无量纲分布曲线（图6）。由图6可见，在沿程10个断面上，除进口附近回流区高度比试验值小外；其他地方都符合良好，并且比255页资料的计算结果更符合实际，特别是在上壁面附近流速最大值处更为明显。

2.排取水问题的计算

本文就宽排取水口问题作了多种工况的计算，这里仅给出图7作为示例。其中所用符号参照图4，只是x、y坐标分别应改为无量纲值x/D、y/D。Δh为相对于G点的水深，单位为mm；ΔT₀为排水口DE处的温升；Re指DE处的雷诺数；i_b为河道纵向底坡；α=Q_上/Q为河道来水量与电厂热水排放量之比，Q的单位为m³/s，ζ=L/D，ζ₁=B₁/D，ζ₂=B₂/D。

当上游来水比电厂取水为小时（即α&lt；1），排水口处的水流势必上溯到取水口，回水在全断面发生（图7）。因排取水口之间无丁坝，即使α&gt；1，在一定范围内，仍有回水存在，但回水紧贴排、取水边岸。受出口处强烈扰动掺混影响，在各种工况下，下游回流始终存在。水深的变化甚微，完全符合刚盖假定。当有极缓的正纵底坡时，水深沿程略有增加；平坡时，略有减小。当上游回水存在时，取水口的温度上升；否则，温度几乎维持平衡水温不变。这说明在近区，对流作用是主要的，而扩散影响相对较小。

图8给出了流线及无量纲k、ε等值线。这是窄排取水口问题，更符合冷却水工程实际。由于问题本身的复杂性，计算模拟困难更大。现阶段本文的计算结果定性与文献［3］的试验结果相符：排、取水口上、下游各存在一个回流区，但回流区的大小及位置等，还不能定量符合；对此种特定情形，计算的稳定性与收敛性也有待进一步改善。图8（b）、（c）分别给出了k、ε的分布。由图可见，k、ε在排、取水口处最大，说明这里的湍流最为强烈，从而湍能耗散率亦最大，符合物理规律。

五、结束语

k-ε这类精细模型用于近区小尺度问题比较成熟，现已有不少成功的算例>^［2］，本文的突扩通道，亦是典型的近区问题，计算结果与试验资料符合良好，是对k-ε方法可以成功用于近区的又一次证实；k-ε模型用于全场以建立整区模型，本文亦作过尝试（见排取水问题的计算），结果表明，这种尝试是可行的。然而，由于对远区大尺度问题，对水深较大、水流分层的问题，流散的影响与湍流的影响具有同等意义，有时它们的影响还相互抵偿>^［6］。因此，如何处理好流散效应，如何综合考虑流散与湍流的相互作用，是建立整区模型的关键所在，这也正是今后努力的方向。

k-ε方法目前仍处于发展阶段，不少问题尚不成熟。对于实际工程中复杂的边界条件，k、ε的边界处理这有待进一步改善。另外，k-ε模型增加了k、ε两个输运方程，显然要比积分模型或扩散模型更费机时。因此，目前离冷却水问题的实际应用尚有一定距离。但作为基本研究，方向无疑是正确的。对于湍流结构，k、ε的变化规律等，应进行试验研究，才能更好地了解湍流的物理本质，从而推动数学模型的发展。

Numerical Moddelling of the Intake-discharge Problem in Channels with a k-ε Turbulence Model^>[4]

Abstract

参考文献

>［1］Gerhard H.Jirka，Gerrit Abraham，Donald R.F.Harleman.An Assessment of Techniqes for Hydrothermal Prediction.Technical Report of MIT，1976.

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

>［3］陈惠泉，许玉林.河道冷却水平面问题的研究∥水利水电科学研究院论文集（第7集）.北京：水利电力出版社，1982.

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

>［5］Rastog，A.K.and Rodi，W..Prediction of heat and mass transfer in open channels.Journal of the Hydraulics Division，ASCE，No.HY3，1978.

>［6］Aspects of Modelling Horizonal Momentum Transfer in Shallow.Report on Investigations，R1150 Delft Hydraulics Laboratory，December，1981.

---