不同湍流模型模拟近岸潮流差异性比较

郝嘉凌，刘镇，胡春也

（河海大学港口海岸与近海工程学院，江苏 南京 210098）

0 引言

近岸水流在大部分时间段为单向水流，潮流在水流运动中起着主导作用，潮流的边界层如明渠单向水流的边界层一样可以得到充分的发展，使边界层所在水深范围很大，甚至可以扩展到整个水深[1]。对于该区域流场，虽然在底床附近仍是低雷诺数的湍流，但是随着垂向尺度增大，逐渐发展的涡旋使湍流变成高雷诺数的状态。因此，简单地将计算管流或者水深很小的明渠流的传统湍流模型应用到该区域，模拟得到的结果是不精确甚至是不适用的。

在第二届国际阻力预测会议中发现，湍流模型的选取是解决湍流问题最重要的影响因素，对计算结果的影响占到15%[2]。雷诺平均模拟方法（RANS），以雷诺平均运动方程与脉动运动方程为基础，通过引入相关假设建立起封闭模型，避免了对湍流运动的直接模拟，降低了大量的计算成本，同时又能保证较高的精确度，是目前最流行的湍流模拟方法[3]。雷诺平均模拟方法数值模型众多，其中k-ε 模型和Mellor-Yamada 模型在海洋水动力模型中应用最为广泛，但Mellor-Yamada模型会低估充分发展湍流的垂向混合，导致水流和湍流耗散之间的相位滞后偏小[4]。针对充分发展的高雷诺数湍流，以Standard k-ε 模型、RNG k-ε模型、Realizable k-ε 模型应用最为广泛。另外，对于床层附近的低雷诺数区域湍流，不同壁面函数的选取对模拟结果有很大影响，针对这3 种k-ε模型，标准壁面函数能表现出对网格和流动良好的适应性[5]。

潮流的一些重要特征量能反映流场的相关特性，除了能较为直接地描述水流特性的流速剖面和床层剪切应力外，湍动黏度的数值和湍动能的衰减可以反映水流的流速梯度及扩散程度[6]，雷诺应力是水流质点脉动产生的且能与平均速度梯度建立关系[7]。本文采用3 种k-ε 模型，用标准壁面函数法解决在底床附近低雷诺数湍流的问题，对潮流进行了一维垂向数值模拟，通过对模拟结果的分析，得出计算不同的湍流特征值时适用的模型及其优缺点，为进一步研究沿岸泥沙的迁移扩散、侵蚀沉积过程、悬浮浓度分布等提供参考[8]。

1 控制方程和湍流模型

近岸水流的一维垂向的基本方程为连续方程和动量方程[9]。

连续方程为：

动量方程为：

其中：

式中：u为时均水平速度；t为时间；ρ为流体密度；w为时均垂向流速；z为垂向坐标，在床底z=0，在水表面z=h+A，h为平均水深，A为潮流振幅；P为作用在流体微元体上的压力；μt为湍流运动黏度；τij为雷诺应力；分别为x（平行潮流流速方向）、z方向上的速度时均值。

1.1 Standard k-ε 模型

Launder 和Spaldingt 提出了基于湍动能k和湍动耗散率ε的Standard k-ε 模型，对于不可压流体，与之相应的输运方程分别为：

式中：湍动黏度μt=ρCuk2/ε，Cu为黏度系数，取Cu=0.09；μ为流体动力黏度；P为湍流动能的产生项，P=μt（∂u/∂z）；σk、σε分别为湍动能k、湍动耗散ε对应的普朗特数，取σk=1.0，σε=1.3；C1、C2为模型常数，取C1=1.44，C2=1.92。

1.2 RNG k-ε 模型

RNG k-ε 模型的湍动能k和湍动耗散率ε分别对应的输运方程如下：

式中：μeff为修正湍动黏度；Cu为黏度系数，取Cu=0.085；、η为模型系数；Eij为时均应变率；模型常数η0=4.377，β=0.012，αk=αε=1.39，C1ε=1.42，C2ε=1.68。

1.3 Realizable k-ε 模型

Realizable k-ε 模型的湍动能k和湍动耗散率ε分别对应的输运方程如下：

其中：

式中：模型常数σk=1.0，σε=1.2，=1.9；、A0、As为模型系数，A0=4.0；φ、W、U*为通用变量；为时均转动速率张量，对于无旋流场为0；ωk为角速度。

2 3 种k-ε 模型在近岸潮流计算中的应用

2.1 计算区域

英国Menai 海峡连通Caernarfon 湾和北爱尔兰海，全长约20 km，平均宽度约为800 m。海峡受以半日潮M2 为主的潮流影响，潮流振幅为2.25 m，大潮时，海峡西南入口的水流速度可达2.5 m/s[10]。当水流流向为东北方向时，海平面处于上升状态，该过程为涨潮（定义为x正方向）。当水流流向为西南方向时为落潮（定义为x负方向），水道见图1。

图1 测点位置图Fig.1 Location map of measuring points

采用1 200 kHz 的ADCP 布置在海峡截面中部的底床上，对该水道的垂向流速进行了测量。测点平均水深12.5 m 左右，床层以上1 m 以及水面下1 m 区域是没有数据的，垂向数据点间隔距离为0.5 m。本文模拟时长远大于实测时长，对一个潮周期的模拟值和实测数据进行比较分析。

2.2 数值方法

Menai 海峡受强潮影响，含沙量较小，淡水输入量小，且湍流发展充分，无分层现象[10-11]。因此针对计算区域的水流，不考虑盐度分层及悬沙分层的影响，对水流进行一维垂向数值模拟。

在以上3 个模型中，对k 方程和ε 方程求解时都采用理论上无条件稳定的Crank-Nicolson 时间积分方案。由于潮流水深随时间在变化，网格数设置为nz-1不变，网格尺寸Δz为变化值，为了避免造成稳定性问题，选定足够小的时间步长Δt=0.5 s。变量u，ρ，μt，k和ε设置在网格中心，只有w设置在网格底部。网格设置如图2 所示。

图2 垂向一维水动力模型中设置的网格和尺寸Fig.2 Mesh and dimensions set in vertical one-dimensional hydrodynamic model

模型采用无滑移不可入固底边界条件，即z=0 处u=0，在无风情况下的自由表面z=h+η处，流速梯度du/dz=0，us=ks=εs=0。为了计算精确，将垂向水流分为黏性层、过渡层和对数层，由于黏性层近似层流，k-ε 模型对于该区域是不适用的。对黏性层的k和ε值进行线性处理，在黏性层与过渡层交界，采用标准壁面函数令该处的湍动能，湍动耗散率，其中：u*为摩阻流速；κ为卡门常数取为0.41；Zb为黏性底层厚度，调试后可取Zb=0.1 m[8，12]。

2.3 计算结果与分析

2.3.1 流速剖面

模型对恒定流情况进行了模拟，选取Gonzalez等人用1 200 kHz ADCP 在Chicago Sanitary and Ship Canal 测量的稳定明渠流数据进行验证[13]。比较3 种湍流模型在恒定流和潮流情况下计算流速剖面的差异，见图3、图4。

图3 恒定流条件下3 种模型计算的流速剖面Fig.3 Velocity profiles calculated by three models under steady flow conditions

图4 潮流条件下3 种模型计算的流速剖面Fig.4 Velocity profiles calculated by three models under tidal current conditions

从图3 模拟结果与实测数据可以看出，对于恒定流而言，3 种模型计算的流速剖面结果很相似，与实测数据都很接近。然而近岸水流以潮流场为主，且涨、落潮峰值和潮流时域特性对港工建筑、工程整治措施影响显著[14-15]。潮流条件下3种模型计算的流速剖面结果有较为明显的差异，结合表1 中流速方差值可知，RNG k-ε 模型和Realizable k-ε 模型相比于Standard k-ε 模型对潮流流速剖面模拟更为精确，尤其是在落急时刻Standard k-ε 模型误差较大。对于恒定流和潮流模拟，模拟结果与实测数据的误差随着水深增加而变大，且模拟值较实测值偏大，这可能是由于模型中未考虑风应力条件导致的。

表1 潮流流速剖面与剪切应力的模拟结果与实测值的方差Table 1 Square deviation of simulated and measured values of tidal current velocity profile and shear stress

确定垂向流速为对数模式u=，通过对实测数据的拟合可得到摩阻流速u*，水流底部剪切应力，因此流速剖面和底部剪切应力可以相互验证。对一个周期内的潮流底部剪切应力进行了模拟，结果见图5。结合表1 中潮流底部剪切应力的方差值可以看出，在涨潮过程，3种模型的模拟结果与实测值都很接近，而落潮过程的模拟结果在相位和数值大小上与实测值有一定的误差。整个涨落潮过程，Realizable k-ε 模型模拟最为精确，RNG k-ε 模型次之。

图5 潮流底部剪切应力Fig.5 Shear stress at the bottom of the tidal current

2.3.2 雷诺应力

在涡流模型方法中，不直接处理雷诺应力项，而是引入湍动黏度，把湍流应力表示成湍动黏度的函数[7]。涡黏假定建立了雷诺应力与平均速度梯度的关系，即：

式中：δij为“Kronecker delta”符号（当i=j时，δij=1；当i≠j时，δij=0）

雷诺应力模拟结果见图6。从图6 可以看出，RNG k-ε 模型和Realizable k-ε 模型能对涨潮情况下的潮流雷诺应力精准模拟，Standard k-ε 模型模拟结果误差较大。对于雷诺应力分布比较分散的落潮情况，3 种模型模拟的结果与实测值都有较大的误差。在潮流底部，模拟值急剧减小，是由于模型对于垂向上进行了分层，从对数层计算到黏性底层时，水流接近层流，流速梯度骤减，这与实际的湍流情况是相符的。

图6 雷诺应力模拟结果Fig.6 Reynolds stress simulation results

2.3.3 湍动能与湍动黏度的无量纲数

湍动能与湍动黏度的无量纲数模拟结果见图7、图8。湍动能的无量纲数、湍动黏度的无量纲数N′=μt/κu*h，其中，湍动能的产生项P=。从图7 可以看出，对于湍动能的无量纲数的模拟，在涨急时刻，RNG k-ε 模型的模拟结果要明显优于其他两种；在落急时刻，3种模型模拟的结果与实测数据趋势一致，RNG k-ε模型模拟结果较优，但仍有一定误差。从图8 可知，对于湍动黏度的无量纲数的模拟，在涨急时刻，Realizable k-ε 模型模拟结果与实测数据最接近，其他2 种模型的计算值相对实测值较小；在落急时刻，3 种模型模拟结果与实测数据趋势一致，湍动黏度的分布规律符合一维垂向湍动黏度的标准抛物线形[12]，但存在一定误差，表现为在底部和表面位置模拟值偏大，而在中部位置模拟值偏小。

图7 湍动能的无量纲数模拟结果Fig.7 Dimensionless number simulation results of turbulent kinetic energy

图8 湍动黏度的无量纲数模拟结果Fig.8 Dimensionless number simulation results of turbulent viscosity

3 结语

根据以上模拟结果与分析，可得如下结论：

1）3 种模型在涨急时刻的特征值模拟结果比在落急时刻的精确。

2）对于恒定流流速剖面的模拟，3 种模型的模拟值几乎没有差异且均适用；对于潮流流速剖面模拟，RNG k-ε 模型和Realizable k-ε 模型模拟都比较精确。

3）在涨潮过程中，对于底部剪切应力、雷诺应力、湍动黏度的无量纲数，选用Realizable k-ε模型计算最精确，RNG k-ε 模型模拟湍动能的无量纲数最精确；对于落潮过程，3 种模型能正确地模拟特征值的分布规律，但有一定误差，原因有两方面：①实测切应力数据垂向分布发生局部突变时，模型没有考虑此效应，落潮切应力垂向分布发生突变主要受涨落潮不对称和地形变化影响；②由于该区域海峡呈收缩状，涨潮水流收束，落潮水流分散，模型只计算了x方向特征值而忽略横向水流的影响。