基于两相湍流珊瑚礁地形的波浪数学模型

邹学锋 朱良生 张善举

(华南理工大学土木与交通学院，广东广州510640)

珊瑚礁广泛分布于太平洋和印度洋的亚热带地区，其中也分布于我国南海地区。珊瑚礁的特点是礁坪水深很浅，常伴有泻湖、凹地，并且礁前斜坡陡峭。近年来，为了开发南海的自然资源和提高南海的经济价值，我国在南海诸岛上建立了多个重要港口，它们多建立在珊瑚礁地形上，按地形分有泻湖港、礁缘港、岛缘港、连岛港和人工岛港等，我国在南沙群岛上施工的几处较大的人工岛都包含港口［1］。在岛礁区周围、岛礁之间都已经形成了一些习惯性的航道。波浪是港口工程里常见的水动力现象，当波浪传播到礁面时，波浪在斜坡上发生浅水变形，波高增大，影响周围水域内船舶的航行。另外，波浪在礁缘或者礁坪上破碎，发生增水，对礁坪后方的港口建筑物产生威胁。因此，研究岛礁的水动力特征，对港口的建设和船舶的航行都有着十分重要的意义。

波浪破碎一直是港口工程领域内的重要课题。一方面，波浪破碎会给海床、海岸建筑物施加额外的作用力，带动泥沙输运，改变着海岸线的面貌;另一方面，波浪破碎后引起的增水和波高变化对堤岸和港口的设计与建造有着主要的影响。波浪破碎的地方称为破波区［2］。波浪在破碎前的形态发生剧烈变化，并可能在破波区发生多次破碎，最终大部分能量转化为湍流动能。

Galvin［3］将波浪破碎分为3种类型:崩破波、卷破波、激破波。Battjes［4］用 Iribarren数 ξ0来判断波浪破碎类型，ξ0＜0.5时为崩破波，0.5＜ξ0＜3.3时为卷破波，ξ0＞3.3时为激破波。在波浪破碎过程中，波浪开始形成涡旋，湍流的产生是其中的控制性因素之一。关于崩破波和卷破波的实验和数学模型一直以来都在不断发展。任冰等［5］通过物理实验模拟了规则波在陡峭珊瑚礁地形上的传播和破碎变化过程;李绍武等［6］通过物理实验研究了不规则波况下岸礁礁坪上的增水;诸裕良等［7］通过物理实验模拟了复合坡度下珊瑚礁地形上的波浪破碎。另一方面，数学模型能提供波浪破碎时的具体水动力学变量值，有些在实验中不易获取，如流速和湍流动能分布等。在模拟波浪破碎的数学模型中，主要有Bousniessq类模型，基于纳维埃-斯托克斯方程的非静压模型和VOF模型，Bousniessq方程是垂向积分方程，模拟波浪破碎时需要在动量方程中加入一个耗散项来将破碎过程参数化，但是参数不易选取［8］。而非静压模型不能处理破波区的掺气作用和水体分离，VOF模型相较于前两者能较好地模拟考虑空气作用时的波浪破碎，以及波浪破碎时水气界面上的湍流分布。Zhang等［9］使用改进的Boussinesq方程模拟了极陡珊瑚礁地形上的波浪传播;房克照等［10］将非静压模型分解为静压步和非静压步两个过程混合求解，处理了海岸区域常见的波浪破碎和水-陆动边界问题。在以往的波浪破碎模型中一些学者采用了单相流模型［2，8，10-11］，但是也有些学者认为单相流模型无法考虑波浪破碎中卷入空气的能量耗散，会造成湍流动能整体偏大［12-14］。Xie［15］在两相流中采用了 k-ε 模型模拟了崩破波和卷破波;Jacobsen等［16］在两相流中使用k-ω模拟了波浪破碎。k-ω SST模拟存在反向压力梯度的流动场有较好的表现［17］。不过，在以往的波浪破碎模型中，使用标准的湍流模型会造成波浪在水面处的过度衰减，本文采用k-ω SST湍流模型，并按照Brown等［18］的方法建立了两相湍流模型，模拟了波浪在珊瑚礁水域的波高、水位和湍流动能的变化规律，数学模型使用Yao等［19］的物模结果进行验证。

1 数值模型

1.1 控制方程

本文采用OpenFOAM中interFOAM的两相流求解器进行模拟，并使用Higuera等［20］的速度边界造波和主动消波边界来实现入口的造波和出口消波。控制方程如下:

式中:U为速度矢量;t为时间;p*为动压力;ρ为计算单元体积密度;g为重力加速度矢量;X=(x，y，z)，为位置矢量;σ为表面张力系数;κ=，为界面的平均曲度;α为体积分数函数;μeff为有效动力黏度，μeff=μ +ρνt，μ 为分子动力黏度，νt为湍流运动黏度，由所选的湍流模型决定。流体体积方程为

其中，α为计算单元中水的体积分数，是个有界值，即0≤α≤1。

在OpenFOAM里没用压缩差分格式计算自由表面，而是加入了人为压缩项·［Urα (1-α)］，使得空气与水的交界面保持尖锐，因此修改后的α方程变为

式中，Ur=min［cαU，max(U)］，与交界面垂直，cα为自定义压缩系数，决定了空气和水体交界面的压缩程度，增大cα的数值可以加大交界面的压缩性，使交界面更加尖锐，但同时也会产生更大的速度梯度，不利于数值稳定［13］，这里取默认值cα=1。体积分数α方程的有界性通过MULES(多维通用显式求解限制器)来完成，它通过离散散度项中的限制因子来确保α最后的值位于0到1之间。

1.2 湍流模型

在这里，使用标准k-ω SST模型模拟了波浪在珊瑚礁上的破碎，并根据Brown等［18］的建议将单元格密度显式地加入到方程中，建立了两相湍流模型。

标准k-ω SST模型的方程如下:

式中:k为湍流动能;ν为分子运动黏度;ω为湍流内能转化率;Pk为湍流动能源项，其定义见本文2.1节;b1=，S为流体应变率，F2为混合函数［17-18］;常系数β*=0.09;σω、σk、β、γ是混合变量，由混合表达式 (7)计算，默认值见表1，

F1为混合函数，详细定义参见文献 ［17-18］。

表1 混合函数的默认值Table 1 Default values for blender functions

进行显式密度后的k-ω SST模型的方程如下:

1.3 边界条件

左侧边界采用速度边界入口，湍流动能入口边界采用Lin等［8］的方法，即k=(cI)2，其中，c和I分别为波形传播速度和湍流强度。计算区域内部的湍流动能初始值与入口处相同，入口边界湍流内能转化率，其中，这里采用Lin等［8］的取值0.0025，底部边界的k和ω使用壁面函数［17］。底部边界的速度设置为Dirichlet边界条件，压力和体积分数为Neumann边界条件。顶部的速度、压力和体积分数都为Dirichlet-Neumann混合边界条件。边界条件的具体设置详见表2。

1.3.1 造波边界

这里使用Higuera等［20-21］的速度边界造波法，在入口处使用速度边界造波，并能主动吸收从岛礁反射回来的波浪。在OpenFOAM中，为了方便造波，假定重力方向为Z轴负方向，同一个造波边界面的最低点处在同一水平线上，但不同的造波边界面可以在不同水平线上。通过不同的波浪理论，能够计算出任意坐标位置上自由水面的高度和质点速度，称为Dirichlet边界条件，而压力采用Open-FOAM里的NeuMann边界条件来计算。1.3.2 消波边界

表2 OpenFOAM中边界条件的设置Table 2 Settings of boundary conditions in OpenFOAM

模型在入口边界和出口边界使用主动消波。主动消波是物模与数模中其中一个重要的消波方法。在宽阔的海域里，波浪反射的地方距离研究区域甚远。但在物模与数模中受到空间的限制，如有限长度的波浪水槽和受计算机性能限制的有限网格数的计算单元体反射的波浪会影响研究区域的实验结果。以前普遍采用的是在波浪出口边界设置海绵层、多孔介质等被动消波区域，或者是人为增加一个消波阻力。被动消波存在缺点，如对于长波无法在消波层上消除，因为它不破碎，因此仍然会发生反射。在数模中，学者多使用基于雷诺应力平均的纳维埃-斯托克斯方程的多孔介质和松弛区域的消波方法［22-23］。另一种方法是主动消波［24］，这种方法最初是在推波板造波的物模实验中实现，该造波机通过测量波高来调整移动速度，以此产生需要的波浪，同时防止反射至造波处的波浪二次反射影响造波。这种主动消波的方法也适用于只用来消波的出口边界。在数模中，由于波浪边界条件固定，所以这种方法通过加入校正速度边界的方法来实现。在二维波浪数模中，为了消除反射波，边界必须产生一个与入射波的速度相同，但方向相反、指向计算区域的波浪，这也是主动消波法的实现原理。本次数模在入口、出口处设置主动消波。

2 数模设置及数值结果分析

本文采用Yao等［19］的实验进行数模的验证。实验在长36m、宽0.55m、高0.6m的波浪水槽中进行，其中水深0.45m，在水槽一端设置平移式推波板造波，在距离推波板16.35m处设置1∶6的斜坡，斜坡由PVC材料制作，斜坡长2.1m，斜坡后接水平平台，水平平台长9.9 m、距离水槽底部0.35m。在波浪水槽的两端均设置多孔介质材料以减少波浪反射的影响，具体设置详见文献［19］。

本文的数模计算区域示意图如图1所示。二维数值波浪水槽长25m、高0.75m。选定与浪高仪相同位置的12个测点的水体体积分数计算波高，计算区域的坐标原点选在造波边界入口与底部边界的交点处。12个测点的坐标值如表3所示。Open-FOAM有两种网格划分工具，分别为blockMesh和snappyHexMesh。blockMesh用于创建由8个顶点、6个边界面组成的“块”。“块”的生成需要定义8个顶点的坐标位置和边界面的边界条件，这里使用blockMesh生成结构矩形背景网格。snappyHexMesh使用三角形构成的STL文件切割blockMesh生成的六面体。这里使用snappyHexMesh生成斜坡和礁坪面，并移除斜坡以下的非计算区域的网格。整个计算区域横向网格单元长度为0.02m，纵向网格长度为0.01m，在水面处附近进行横向加密，加密后的横向单元长度为0.01m，纵向单元长度为0.005m。斜坡位置上的网格示意图如图1所示。坐标原点选取在测点的位置见图1。

图1 数值模型计算区域和浪高仪位置Fig.1 Computational domains of numerical model and the position of wave gauges

表3 浪高仪的坐标位置Table 3 Coordinates of wave gauges

本文研究的规则波的入射波高为0.095 m，周期为1.25s，为斯托克斯二阶波［25］，计算得Iribarren数为2.36，因此属于卷破波。模型计算耗时60s。

2.1 湍流动能与雷诺应力

在雷诺应力模型中，湍流动能表示为

式中，U′为脉动流速。在k-ω SST模型中，湍流动力黏度表示为

湍流动能的源项表示为

对于波陡较大的波浪 (波高与波长的比值H/L＞0.07)，在水面处容易产生非常大的湍流动能，这是因为空气与水有很大的密度比 (ρa/ρw=1∶1000)。

2.1.1 深水区湍流动能

图2分别是标准k-ω SST模型和两相k-ω SST模型下4个波长范围内的湍流动能分布云图，图中的左端是入口边界，黑色线表示是水面。标准k-ω SST模型在波浪破碎前段的水面处存在着湍流动能过度偏大的现象，这会引起波浪在传播过程中波能过度消耗，使得传播到破碎带前的波高预估不足;两相k-ω SST模型在水面处限制了湍流动能发生不合理的增长，因而也降低了湍流动力黏度，避免在使用湍流模型后波高发生衰减。

2.1.2 破波带湍流动能

图3分别是标准k-ω SST模型和两相k-ω SST模型下斜坡上的湍流分布云图，其中t/T为无量纲化时间。以41.5s开始为例，两种模型在破碎前的自由表面湍流动能分布有很大的差异。在波浪破碎位置上，两相k-ω SST模型的波浪破碎前的湍流动能很小，破碎后的湍流动能从破波点开始急剧增大，在波峰前端位置形成具有较大湍流动能的区域。破碎后的湍流动能继续向礁岸传播，并在礁坪上的传播过程中衰减。此外，在波浪破碎前，两种模型在右端礁坪区域的湍流动能都比较大，这表明在下一个波浪发生破碎时，前一个波浪产生的湍流动能还未完全耗散。

图2 深水区湍流动能云图Fig.2 Turbulence energy contours in the deep water area

2.1.3 雷诺应力

其中I为单位矩阵。因此，在雷诺应力模型中，计算雷诺应力改为计算湍流动能和湍流运动黏度，并由式 (12)确定雷诺应力。雷诺应力的大小=。图4为雷诺应力的分布云图，从图中可以看出，波浪在发生卷破后水面附近的雷诺应力增大，雷诺应力的空间分布类似于湍流动能，主要集中在破碎位置的波峰处。

2.2 波高与平均水位

图5为波高和平均水位的数模模拟结果和实验值的对比。图中的RMSE代表均方根误差，可按下式计算:

式中:N为测量点总数;Qic为第i个测量点的变量模拟值;Qi

图3 标准k-ω SST模型和两相k-ω SST模型下斜坡上的湍流动能云图 (仅水体)Fig.3 Turbulence energy contours on the slope with standard kω SST and multiphase k-ω SST models(Only for water phase)

m为第i个测量点的变量实验值。图5(a)给出了波浪在标准k-ω SST模型和两相k-ω SST模型下的数模结果与文献 ［19］的实验结果沿珊瑚礁方向传播的波高分布。其中，RMSE文献、RMSE标准、RMSE两相分别为由文献 ［19］、标准k-ω SST模型、两相k-ω模型获得的均方根误差。从图中可以看出，两个模型在破碎前段波高都与实验值相差不大，波高的沿程起伏变化是由于礁前陡坡对波浪有反射作用，实际波面是由入射波和反射波叠加而成，因而在深水区形成了不完全驻波。值得注意的是，相较于两相k-ω SST模型，标准k-ω SST模型在礁前发生了明显的衰减，这是考虑了湍流模型后，湍流动能在水面处产生了不合理的源项引起的。当波浪到达斜坡顶端附近时，波高达到最大，波浪发生破碎，波能急剧耗散，波高急剧变小。由于底部是光滑面，底摩阻引起的波能损失较小，而在礁坪上，波浪破碎后的波高由当地水深决定。Yao等［19］的数模中模拟湍流采用的是混合长度模型，从图5(a)中可以看出，该模型在破碎前的波高预测值较实验值和k-ω SST模型都偏大，这是由于混合长度模型低估了破碎前的湍流强度，使得湍流运动黏度νt偏小［26］。图5(b)为波浪在标准k-ω SST模型和两相k-ω SST模型下的数模结果与实验结果沿珊瑚礁方向传播的平均水位分布。平均水位通过对自由表面的历时线取时间平均求得。以静水位为纵轴零点，正值表示增水，反之为减水。从图中可以看出，两种模型在深水区的模拟结果相差不大，在礁前都发生了减水，平均水位在礁坪边缘处达到了最低值，两相k-ω SST模型的礁前减水更明显一些。数值模型的水位回升比实验值缓慢一些，这也是由波高衰减得比较缓慢引起的。而在礁坪区域，两相k-ω SST模型的增水比标准k-ω SST模型要更大一些。

图4 标准k-ω SST模型和两相k-ω SST模型下斜坡上的雷诺应力云图 (仅水体)Fig.4 Reynolds stress contours on the slope with standard k-ω SST and multiphase k-ω SST models(Only for water phase)

图5 数值模型与实验值的波高和平均水位对比Fig.5 Comparisons of wave height and mean water level between numerical and experimental results

2.3 波浪传播变形

图6 比较了两种数学模型下，在具有代表性的浪高仪位置 (#1、#3、#5、#7、#9、#11)处的自由表面的时间序列。从整体上看，两种模型在连续6个周期内的波形基本保持一致，这说明k-ω SST模型在模拟波浪破碎的过程时，不仅在破碎前有很好的稳定性，在破碎后的区域也能够保持波形稳定。浪高仪#1位于破碎前段，表现为波峰和波谷对称的规则波，两种数学模型与实验值基本相符。当波浪传播到浪高仪#3位置时，底部为礁前斜坡，浅化作用使得单个周期内的波峰和波谷开始出现不对称特征，具体表现为波峰向前倾斜、波谷变长。值得注意的是，标准k-ω SST模型的波高明显偏低，而且出现了相位侧移。浪高仪#5位于破碎区内，波高急剧变小，波峰和波谷严重不对称，实验结果出现明显的锯齿状，周期间的波形差异较大，而两个数值模型在破波区内的波能耗散比较缓慢，与实验值相比都偏大，这表明数值模型在破碎带的能量耗散得比实验值慢。此外，本文采用体积平均法计算波高，这在破碎带外计算波高确实可行，但在破碎区，由于流体体积法 (VOF)计算出来的波面附近含有强烈的破碎和掺气作用，水面处的一些单元格水体发生飞溅，空气也卷入到了水体中，这一区域的波高计算方法仍然有待探讨。在浪高仪#7、#9处，波浪破碎停止，在礁坪区域重新生成了波高较小的段波。由于礁坪区域水深变浅，受非线性作用影响而呈现显著的非线性特征，具体表现为这两处的波峰变短抬升且向前倾，波谷变长且更加平坦，不对称特征非常明显。

图6 浪高仪#1、#3、#5、#7、#9、#11处波面高度的时间序列Fig.6 Time-series of the surface elevations at six wave gauges#1，#3，#5，#7，#9，#11

2.4 流速分布

图7 是两种湍流模型在浪高仪#7处一个周期内不同时刻的流速垂向分布图。以波浪开始发生破碎的时刻为零点，时间间隔为0.1 s。可以看出，在浪高仪#7处底部存在着明显的下层逆流，且两相k-ω SST模型的逆向流速要更大一些。此外，水平流速梯度在破碎前后的方向相反，底部水体的流速很小，并沿着水面的方向流速逐渐增大。另外，k-ω SST模型下的最大流速略大于标准k-ω SST模型，位置也更接近礁坪面。

图7 标准k-ω SST模型与两相k-ω SST的水平速度垂向分布比较Fig.7 Comparisons of vertical distributions of horizontal velocity between standard k-ω SST model and multiphase k-ω SST model

3 结论

本文使用OpenFOAM中的标准k-ω SST模型和两相k-ω SST模型模拟了波浪在珊瑚礁地形上的破碎过程，并对数值模型的波高、水位、流速进行了分析。总体上看，数模结果与实验结果大体一致，主要结论如下:

(1)两相k-ω SST湍流模型相比标准k-ω SST模型对波高、水位的模拟更加准确。密度修正后的模型在波浪破碎前的层流区域能有效减轻湍流动能的过度产生，使得波浪在传播过程中不会发生不合理的明显衰减，也使得对破碎前的波高预测更加准确。此外，两相k-ω SST湍流模型在礁坪上的增水更加明显，更接近实验值。

(2)在破碎区域，在自由表面附近产生了明显的湍流动能且主要集中在波峰处，这说明在波浪破碎区域应当使用湍流模型，且选择适当的湍流模型对模拟结果精度十分重要。

(3)使用两相流的模拟结果表明，在波浪的破碎过程中，不管是上层空气还是底部水体，都存在着明显的逆向梯度流，使用k-ω SST模型计算得到的波面具有较好的稳定性。