天然河流、海湾的边界曲折、地形复杂，采用坐标变换是解决问题的途径之一。目前多数N-S方程的坐标变换中，流程全部采用逆变分量，这样就增加了方程的复杂程度。于是忽略掉方程中的非正交项，利用正交变换下的方程进行数值求解[7,8]。对于具有复杂边界的海湾及弯曲的河流，坐标变换中很难保证每个点都正交，特别是边界附近。水位变化是水力计算中难点之一，在目前的紊流数学模型中，多简单的利用“冻结法”，这样做将失去对边界出流动模拟的准确性。

本文研究中，采用正交曲线坐标变换生成数值网格，而数值计算中采用非正交曲线坐标下的k-ε双方程紊流模型，这样可以自动修正网格生成中的非正交项。流速除对流项中采用逆变分量，在其余各项中均采用原始分量，这样使得方程书写简单，有利于将各方程写成通用形式，编写的程序变得更规范。作者受Jian Ye同位网格[9]的启发，对普通交错网格做了修改，即采用B型交错网格，使得u,v,k,ε的计算布置在一个节点上，有利于节省计算程序代码，使程序书写更加规范。引入动边界扫描技术，结合紊流模型的壁面函数法，使壁面随着真实边界而变化。数值求解时，采用控制体积法离散方程，运用SIMPLEC算法，使计算的流场更符合实际流场。

1 数值网格
本文对计算区域用Laplace方程实施坐标变换，生成正交的贴体网格，控制方程：

方程组(1)可采用SOR方法求解，然后用正交曲线网格边界正交化处理的三次样条插值边界滑动法[10]，处理其边界处的网格，以提高其正交性。与其它方法相比，此方法虽在一定程度上提高了边界处网格的正交性，但难以保证完全正交，如果再把曲线坐标系下N-S方程中的非正交项去掉，由斜交网格引起的计算误差是不容忽视的。为了修正此项误差，本文采用非正交曲线坐标下的平面二维k-ε双方程紊流模型，模型能自动修正网格生成中的非正交项。

2 数学模型
笛卡儿坐标下深度平均的k-ε双方程紊流模型的通用微分方程：

方程(2)转换到曲线坐标(ξ、η)下，仅在对流项中使用流速的逆变分量，而在其它项中使用原始变量，这样既简化了方程，又使所有方程仍可写为曲线坐标下的通用方程，模型的微分方程可写为如下通用形式：

式中：Φ为所求问题的因变量；U和V分别为直角坐标下流速u和v的逆变分量，仅在对流项中出现；ΓΦ为扩散系数；SΦ为源项。当Φ表示某一特定量时，ΓΦ和SΦ对此特定量有特定的意义和表达式，这时方程(3)亦赋予特定的意义，模型的控制方程组如表1所示。

式中：u和v分别为笛卡儿坐标下x和y方向的深度平均的流速分量；k和ε分别为深度平均的紊动动能和紊动动能耗散率；H为水深；zs为水位；zb为河床高程；ν为分子运动粘性系数；νt为漩涡运动粘性系数；τbx和τby分别为x和y方向上的底摩擦力μeff为有效粘性系数；Cμ、Cε1、Cε2、σk和σε为经验常数，其中Cμ为0.09，Cε1为1.44，Cε2为1.92，σk为1.0，σε为1.3。

3 模型离散

如图1所示，本文采用B型交错网格，即将u,v,k，ε及各标量布置在网格中心，将水深相关项H、z?s、z?b项布置在网格的左下方节点上。这样布置有以下两个优点：①有利于边界条件的引入；②在编写程序u,v,k,ε各方程共用一套系数，可节省程序代码，使程序简单明了。
本文用控制容积法，采用混合格式来导出离散方程。其边界条件为：
在进口边界，所有边界条件都按本征条件给出，即：u=u0,v=0,k=k0,ε=ε0；
在出口边界，给定出口处的水位zs，令v=0，其余变量都取法向导数为零，即：

在水陆边界，在计算河流、河口、海湾的非恒定流场时，水位不是恒定不变的，而且每一时间步长的水位可能都不一样，这就带来了水陆边界的移动，如果每一时间步长都进行坐标变换而生成新的网格是非常不经济的。于是本文在前人“冻结法”的基础上，提出了适合k-ε双方程紊流模型的“移动边界的壁面函数法”。如图2、图3，水位变化后，使边界处的网格干出，原来的岸边界12、34变成了新的岸边界PQ、ST，而壁面函数仍然布置在12、34岸边界处；水域内岛屿由于水位的变化而出露，其岛屿内部按“冻结法”来处理，其边界ABCD周围也需用壁面函数法来处理。本文采用“水位扫描法”来判断新的水陆边界：①从边界12向34扫描，如水深不为零，则令其为起始点；②从边界34向12扫描，如水深不为零，则令其为终点；③然后扫描水域内部，确定岛屿的边界ABCD。之后对水陆边界实施壁面函数法。

在壁面边界，k-ε双方程紊流模型是高Reynolds数模型，只适用于充分发展的紊流，而在近壁处，粘性效应起主导作用。另外，由于近壁处物理量变化很快，需要非常密的网格，才能真实地模拟实际流动。因此本文采用壁面函数法[11]，即用半分析的方法得到的解来近似由壁面到紊流核心区之间的流速、紊动动能、紊动动能耗散率的分布规律，将壁面的影响(如壁面应力)附加到差分方程中(应先将有关的边界系数置为零)[12]。

4 SIMPLEC 算法求解模型

SIMPLE 算法的基本步骤如下：(1)计算坐标变换的相关系数；(2)根据进、出口边界的水位，给定全场水位，计算全场初始水深H*；(3)解动量方程求un,vn;(4)解k-ε紊流模型求k、ε；(5)解水深校正方程求H＇，并修正水深 H=H*+αhH' ；(6)修正流速u、v及有效粘性系数μeff，然后重复(2)～(5)步直至收敛。

模型求解中，为了利于非线性迭代的收敛，计算中采用亚松弛技术。计算中采用ADI技术和TDMA算法，以连续方程的误差小于给定的值作为判断收敛的依据。在恒定流计算中，以相邻时间步之间全场各点流速误差的最大绝对值小于给定的值作为达到定常状态的依据。

5 模型应用

为了使本文数学模型具有代表意义，本文采用天然连续弯道河流——美国Colorado洲Fall River部分河段的实测水流资料进行验证。该河段长120m，平均河宽8.5m，其河床地形资料如图4所示。计算流量为4.1m３/s，上游水位2.85m，下游水位2.61m。用坐标变换的方法将物理平面转换到计算平面。如图5所示，整个河段可划分为91×20个网格，网格长度平均0.4m，宽度平均0.09m。

本模型计算所选用的糙率平均为0.62，并在计算中自动调整[13]。本文属于恒定流计算，所采用的精度为10-5，计算经过1500步收敛，得全河段流场图，如图6所示。
为了进一步对本模型进行验证，作者选用具有实测资料的6个断面的流速和水位与计算值进行比较。由于本河段属于冲积性多沙河流，计算中糙率以及河底高程很难调整到与实际情况完 全相符，使得计算中存在一定的误差，而且在一些断面上个别点的计算结果与实测值存在着较大的离散，但总体看来计算结果能够满足工程需要。流速比较见图7，水位比较见图8。

6 结论

本文建立了非正交曲线坐标下的k-ε双方程紊流模型，可以对具有不规则边界的水域，如海湾、天然连续弯道水流的进行较好的模拟。模型中采用了一系列针对天然水流的处理技巧：特殊的交错网格、与“水拉扫描法”结合的壁面函数法的动边界处理、保留坐标变换后模型方程中的非正交项等等，均提高了模型的实用性，显示了本模型向水库、河流、海湾等水域的污染预报模型发展的前景。