1　数学模型

1.1　基本方程

基于N-S方程与布辛涅斯克浮力流动的假定，本文导出了湍流、水温和异质浓度的二维耦合方程。通用形式为

表1中所列各物理量均为无量纲化形式，其中vt=CK2/X为涡黏性系数；Re为雷诺数；Ri、Ric分别为温度、浓度里查逊数；Pe、Pec分别为温度、浓度皮克勒特数；Prt、Sct分别为湍流普朗特数与湍流施密特数；Ec为艾可特数；Q为热源；k为COD降解系数；c、ek、eX、C1、C2、C3、C4为模型常数；。

1.2　边界条件

图1为研究剖面示意图，断面⑳、⑲、⑰、⑮、⑭均有实测资料，本文将⑳和⑭断面的实例资料作为确定进出口边界条件的依据，⑲、⑰、⑮断面用于误差评定。

（1）进口。水温、COD浓度分布以实测值给定，速度u根据实测水面最大流速取指数分布，即u=umax（1-y/H）m，v=0，为断面平均流速，H为断面最大水深。

（2）出口。u、v、水温和COD浓度分布确定方法同进口。为简便起见，K、X的边界条件按局部坐标单向化处理[2]。

（3）水面。采用刚盖假定和对称面条件，∂u/∂y=0，v=0，∂K/∂y=0，∂x/∂y=0，∂c/∂y=0。水温的水面边界条件可以两种方式考虑：一是考虑水体与大气的热交换，取∂T/∂y=-hɑ/（dCpDZ），但这需要大量的水文气象资料[6]；二是直接按文献[4]的方法：T表=T气+Δb，按一类边界条件给定，或直接用实测水面温度值给定。

（4）库底。用壁函数处理[1]，将靠近壁面的第一个内点P置于对数分布区域，设P处平行于壁面的速度为为壁面摩阻流速；k=0.4～0.42为卡门常数；E=9.0；y+=（ywu*）/v；f为壁面切应力；d为水密度；v为运动黏性系数；yw=0.002～0.005。

假定在壁面附近湍流处于局部平衡状态（即对流与扩散可忽略，产生与耗散相平衡），则P处的湍动能与其耗散率为

温度与浓度采用绝热与无渗透假定[3]，即

1.3　坐标变换及网格生成

微分方程坐标变换的基本思想是通过求解椭圆型微分方程的边值问题，确定与计算平面内各节点相应的物理平面上的坐标。

由文献[2]知，物理平面上网格节点可由泊松方程来确定：

边界上的网格疏密控制函数为

由式（2）很容易导出混合有限分析法五点格式[5]，进而可迭代求解物理平面上的网格节点坐标。由于原始变量法求解N-S方程时常用交错网格技术，以避免出现锯齿形压力分布的不合理结果。为了直接形成u、v、P三套网络系统，网格生成时总节点取u网格节点数的两倍。各网格系统的编号方法同文献[2]，即以速度矢量箭头所指向的主节点P的编号为该速度及相应坐标的编号。

1.4　控制方程变换

由x=x（ɑ，Z），y=y（ɑ，Z）将连续方程变换为

式中：Uɑ=yZu-xZg，VZ=xɑg-yɑu，Uɑ、VZ为计算平面上的速度。

对于对流—扩散方程通式（1），∂h/∂t保持不变，其他各导数项可按文献[1]中所给公式代入。源项中若含有导数项，亦按同样方法转换，变换后方程为

式中，J=xɑyZ-xZyɑ。

设A1=（Ju-Γx）yZ/Γ+（-Jv+Γy）xZ/Γ-f，B1=（Jv-Γy）xɑ/Γ+（-Ju+Γx）yɑ/Γ-e，其中e=（yɑDx-xɑDy）/J，f=（xZDy-yZDx）/J，Dx=Txɑɑ-2UxɑZ+VxZZ，Dy=TyZZ-2UyɑZ+VyZZ。

对于u方程有

对于v方程有

对于其他各方程有

各方程源项S类似表1，但其推导方法类似，且中为略去源项S的变换后的形式。边界条件变换可参阅文献[1，2]。

方程（2）以后各表达式中变量的下标均表示对该下标的偏导数，如hZ表示∂h/∂Z，hZZ表示∂2h/∂Z2，hɑZ表示∂2h/∂ɑ∂Z。