摘要：作者针对近岸波浪与结构物相互作用问题提出了一种耦合数值方法，即用时均化的二维雷诺平均的Navier stokes方程-流体体积法模型表达内域流动，用一维Boussinesq方程表达外域流动，通过速度、压力和波面匹配边界条件实现两种数值模型的同步求解。耦合模型中的二维子模型能够较好地表达结构物附近流动的细部特征，包括漩涡结构；一维子模型的计算效率很高，可通过延长其计算域以达到有效地避免二次反射波的影响。所建立的耦合数值模型被证实可应用于幕墙式消浪结构防波性能的研究。
关键词：幕墙式消浪结构 耦合数值方法 VOF方法 Boussinesq方程
以往的研究成果表明，迎浪面开孔的沉箱直立堤可以有效地减小反射波，但消浪室的宽度(即开孔前墙和不透水后墙之间的宽度)一般应达到当地波长的四分之一[1]。如果入射波为涌浪或者其他类型的长波，这意味着理想的消浪室宽度在实际工程上可能无法实现。最近，日本学者提出了一种能有效消减直立堤前反射波的新型——幕墙式消浪[2](curtain-walled dissipater)，其断面如图1所示。设在直立墙前的垂直屏障称为幕墙。幕墙至直立墙的距离B为消浪室的宽度。
图1 幕墙式消浪
幕墙吃水深度用c表示。由于入射波引起的消浪室内水体振荡运动和幕墙下面的涡旋运动使得波能大量耗散，从而实现消减反射波浪的目的。这种新型消浪的主要优点是能够有效地减小消浪室的设计宽度。
数值模拟是揭示幕墙式消浪水力学性能和消浪机理的有效手段。由于在幕墙和直立墙处产生的反射波在造波边界处可能形成二次反射，通常的方法需要在二次反射波传播到物之前停止计算。这意味着计算域的长度必须足够大。然而，在一个很长的立面二维计算域上全部采用粘流波浪数值模型做精细模拟一方面计算工作量很大，另一方面必要性也不充分。为此，本文提出了一个耦合求解策略，即将二维RANS-VOF模型与一维Boussinesq方程模型耦合起来解决问题。在耦合模型中，一维子模型由于其计算效率很高，可以考虑足够长的计算域；二维子模型则能够较好地反映流场的细部，包括粘性对流动的影响。
1 耦合模型的原理
如图2所示，耦合模型是将整个计算域Ω划分为Ω1和Ω2两个子域。这两个子域通过一条公共的重叠带衔接起来。Ω1为包含幕墙和直立墙的近场，流动以二维紊流运动方程，即Reynolds方程(RANS)为控制方程，采用标准k-ε紊流模型封闭，并在近壁区应用壁面函数理论[6]；自由水面的描述采用定义流体体积函数的方法[5]。
图2 耦合模型区域划分
Ω2域内流动的控制方程采用色散性改进的Boussinesq方程[4]的一维形式，经差分离散后得到系数矩阵为三对角矩阵形式的代数方程组，采用追赶法快速求解。在耦合模型中，两个子模型RANS和Boussinesq各自独立求解，耦合的实现体现在重叠带上流动信息的匹配。
为了便于耦合处理，Boussinesq方程和RANS方程均采用交错网格进行差分离散。其中，Boussinesq方程的离散参考了Madsen和Sorensen所用的格式[4]。RANS动量方程中时间项的离散格式为向前差分，粘性项的离散格式为二阶中心差分。为消除数值粘性的影响，动量方程中对流项的离散格式采用了三阶迎风差分格式[3]。差分方程的求解采用了SOLA-VOF方法[5]。其基本思想是：首先用前一时刻的流场计算结果代入动量方程的显式差分格式，求出当前时刻流场的近似值；再通过对压力厨行迭代修正，使得连续方程在一定的精度条件下得以满足，对表面单元要求满足自由表面的动力学边界条件，即通过线性插值确定表面单元中心处的压力值；在完成压力迭代后，再对速度进行校正，然后用校正后的速度值代入k-ε方程相应的差分格式求解紊动动能和紊动动能耗散率；最后，应用施主与受主单元模型计算当前时刻的流体体积函数，确定流体自由表面的位置。
由于动量方程、紊动动能方程和紊动动能耗散率方程对近壁区网格细密程度的要求不同，耗散率方程的要求最严，动量方程和动能方程的要求基本一致，为了既保证解的精度而又不致使网格划分太密，本文在固壁区附近采用了壁面函数方法[6]进行处理。即在壁面附近引入以下关系
其中β0是常数，与壁面粗糙度有关，本文取β0=0.0005；L是特征长度，计算中取为近壁区网格中心到壁面的距离。
合理地设置匹配边界条件，使得内域流动和外域流动在匹配边界处光滑而连续地过渡，是保证子模型耦合的关键。本文的做法如图3所示。匹配边界ΓB是Ω2域的出流边界，Boussinesq模型执行每一时间步的计算之前需要预知该边界上的速度和波面值，由于边界ΓB同时又在Ω1域的计算节点上，于是ΓB边界的速度条件[]ΓB可利用Ω1域得到的流场信息表达如下，
图3 匹配边界附近的差分网格单元
(1)
式中：u为水平速度，F为流体体积函数[5]，δy为垂向的网格步长，j为垂向网格节点编号，Jmax表示垂向网格节点的最大编号。
给ΓB边界的波面赋予匹配条件时数值试验表明，需要利用连续方程反映的水位流量关系给出匹配条件其效果好于直接给定水位过程条件。这是因为通过水位和流量的相互调整，计算域内的反射波可得以减弱。因而，在实际计算中ΓB边界的波面可表示为
利用Boussinesq方程推导过程中引入的近似展开关系，在波面函数η和(深度平均)速度导的结果，ΓV边界上水平速度u\，垂向速度v和压力p的匹配条件可按下面表达式给出：
耦合模型同步求解过程可简单概括如下。首先，考虑Ω2域左边界处的入射波条件，在Ω2域内执行Boussinesq模型，当Ω2域中接近匹配边界的节点上的水平速度值第一次达到10-3m/s量级时，开始在Ω1域内执行RANS-VOF模型。在某一时间步n，执行Boussinesq模型所需要的匹配边界ΓB处的速度和波面条件按式(1)和式(2)给出；Boussinesq模型在当前时间步的计算完成后，随即利用式(3)～(5)计算出匹配边界ΓV处的波面、速度和压力边界条件，并启动RANS-VOF模型；RANS-VOF模型在当前时间步的计算完成后，即按照式(1)和(2)计算出匹配边界ΓB所在位置的速度和波面匹配条件。然后进入下一时间步的计算。
2 计算结果
为了验证计算结果，参考已知的模型试验条件，取耦合模型为等水深(d=0.42m)，计算域全长18.36m.其中，子区域Ω2和Ω1的长度分别为12.58m和5.78m.Ω1域中幕墙厚度q=0.042m，幕墙下端为一45°尖角。幕墙吃水深度按相对吃水c/d=0.3～0.7计算了多个工况。幕墙和直立墙均作为全反射边界处理。消浪室宽度B=0.29m取为固定值。通过改变入射波的周期，消浪室的相对宽度在B/L=0.08～0.16(L为波长)的范围内变化。计算网格的划分情况如下。Ω2域按等步长(Δx=0.042m)划分为300个单元格；Ω1域的垂向也按等步长(Δz=0.021m)划分，但其水平方向是按变步长划分的，共划分了150个单元格。Ω1域中单元编号i=1～120的各单元为Δxi=0.042m，i=121～125的各单元为Δxi=0.035m，i=126～130的各单元为Δxi=0.028m，i=131～150的各单元为Δxi=0.021m.幕墙附近和消浪室区域的网格划分得比较细的目的是为了更有效地分辨幕墙下面涡旋流动的特征。Ω2域左边界入射波的波高H=0.06m.应用合田方法[8]测定反射率系数，在幕墙前大约2倍波长以外处设置了适当间隔的两个浪高仪(如图2所示)记录波面变化。图4为执行耦合模型得到的瞬时(t=19.972s)波面(图中点划线为匹配边界位置).图5显示了波浪反射系数Cr随相对消浪室宽度B/L变化的关系。当幕墙吃水深度c/d=0.5一定时，Cr随B/L变化，在B/L=0.12附近出现了最小值。本文的数值结果同模型试验的结果十分接近，说明1D/2D耦合模型是有足够精度的。
计算过程中观察到了消浪室内水体做整体升降振荡的现象。取消浪室平均水位变化的幅值(HT)与原入射波高(H)之比描述波浪在消浪室中被激励的情况。图6显示了波浪在消浪室中水面波动被放大的比率随消浪室相对宽度B/L变化的情况。可以看出，大约从B/L=0.12开始，消浪室中波浪的放大率随B/L的增大而迅速减小的特点十分明显。这说明，在给定幕墙吃水深度(c=0.12m)和消浪室宽度(B=0.29m)的情况下，波长(L)相对较大的波浪经垂直幕墙透射到消浪室后引起的水体振荡幅度更大。图7显示了幕墙反射系数Cr随其相对吃水深度c/d变化的关系。当消浪室相对宽度B/L=0.12一定时，Cr随c/d变化，在c/d=0.5附近出现了最小值。
H=0.06m，B/L=0.12，c/d=0.5
图4 整个计算域瞬时波面 d=0.42m，c=0.21m，H=0.06m
图5 Cr-B/L关系
d=0.42m，c=0.21m，H=0.06m
图6 消浪室中波浪放大率 B/L=0.12，H=0.06m
图7 Cr-c/d关系
对于幕墙反射系数为最小的情形(B/L=0.12，c/d=0.5)，图8显示了消浪室内的平均水位变化。从计算中可以看出，消浪室内水体的运动呈整体活塞式的运动。幕墙末端附近的复杂流态对自由水面并没有大的影响。图9显示了幕墙附近涡旋流动的演化过程。可以看出，消浪室内常驻着一个逆时针方向的漩涡，其尺度明显受控于消浪室宽度。在消浪室水位由高向低变化的过程中，幕墙开口处外侧有一个顺时针方向的漩涡产生；随着消浪室内水位的回升，该漩涡逐渐进入消浪室。
图8 消浪室内平均水位变化
(H=0.06m，B/L=0.12，c/d=0.5)
受消浪室内常驻漩涡的影响，进入消浪室内的漩涡有尺度减小的趋势。消浪室水位再一次由高向低变化时，幕墙开口处外侧将有新的顺时针方向的漩涡产生。如此反复，涡漩的频繁产生和移动导致了波浪能量的大量耗散，从而实现了消减反射波浪的目的。
3 结语
以上针对近岸波浪与物相互作用问题提出了一种耦合数值方法，即用二维RANS-VOF模型表达内域流动，一维Boussinesq方程模型表达外域流动，通过速度、压力和波面匹配边界条件实现两个子模型的同步求解。耦合模型中二维子模型能够较好地表达物附近流动的细部特征，包括漩涡；一维子模型的计算效率很高，使得大域问题的求解成为可能。所建立的耦合方法能有效地用于研究幕墙式消浪的水力性能。计算结果表明：(1)消浪室相对宽度B/L和幕墙相对吃水深度c/d都是影响垂直幕墙波浪反射系数Cr的重要参数。当c/d=0.5一定时，在B/L=0.12附近Cr为最小；当B/L=0.12一定时，在c/d=0.5附近Cr为最小。也即，当消浪室宽度接近0.12倍波长和幕墙吃水接近0.5倍水深时，波浪反射系数达到最小。这一结果对于工程上减小消浪室的宽度和确定合适的幕墙吃水深度是有重要参考价值的。
H=0.06m，B/L=0.12，c/d=0.5，t0=19.376s
图9 幕墙附近局部流场
(2)幕墙式消浪消浪机理的重要方向是幕墙下面(特别是消浪室内)涡旋的发生。强烈的涡旋流动将波浪能量大量耗散，从而使得消减反射波浪的目的得以实现。(3)计算获得的Cr-B/L关系等结果与实验结果符合良好，说明耦合模型具有良好的计算精度。