波浪爬升是海洋工程中的热点问题^[1]。线性的势流方法只能计算线性波浪的作用^[2]。对于弱非线性入射波，二阶理论会提高预报精度^[3-4]。但仍受限于理论假设，无法准确预报强非线性入射波作用下的波浪爬升问题以及粘性效应^[5]。计算流体动力学(computational fluid dynamics，CFD)方法通过直接求解纳维-斯托克斯(Navier-Stokes，N-S)方程，没有理论假设，可以预报强非线性波浪爬升问题^[6]。但相比于势流方法，粘流方法的计算成本更高^[7]。采用势粘流单向耦合技术可以减小数值模拟计算域，降低计算成本^[8-11]。

本文参考文献[12-13]的模型试验数据，采用势粘流单向耦合技术对单直立固定圆柱的波浪爬升进行了全尺度数值模拟。

1 势粘流单向耦合技术控制方程

在实际物理流动过程中，粘性只在物体近壁面边界层起主导作用，所以在远离物体的流场中可以采用势流求解以提高计算效率。本文采用的势粘流单向耦合技术，基于STAR-CCM+求解器，把计算域分为粘流区域与势流区域，在粘流区求解不可压缩N-S方程，在势流区求解欧拉方程。在势粘流重叠区进行势粘流求解匹配，见图 1。

1.1 粘流区控制方程

本文在粘流区域采用有限体积法对控制方程进行离散，采用流体体积函数(VOF)方法捕捉气体与水之间的交界面。湍流模型选择分离涡(DES)模型^[14]：

$ \begin{array}{l}{\int\limits_{V} \frac{\partial}{\partial t}(\rho \overline{ \mathit{\boldsymbol{v}}}) \mathrm{d} V+\int\limits_A \rho \overline{ \mathit{\boldsymbol{v v}}} \mathrm{d} \mathit{\boldsymbol{a}}=} \\ {-\int\limits_{V} \nabla \overline{ \mathit{\boldsymbol{p}}} \mathrm{d} V+\int\limits_{V} \nabla \cdot\left(\tau+\tau_{t}^{\mathrm{RANS}}\right) \mathrm{d} V+\int\limits_{V} S_{\phi} \mathrm{d} V}\end{array} $

(1)

$ \begin{array}{l}{\int\limits_V \frac{\partial}{\partial t}(\rho \tilde{ \mathit{\boldsymbol{v}}}) \mathrm{d} V+\int\limits_A \rho \tilde{\mathit{\boldsymbol{v}}} \tilde{ \mathit{\boldsymbol{v}}} \mathrm{d} \mathit{\boldsymbol{a}}=} \\ {\quad-\int\limits_{V} \nabla \tilde{ \mathit{\boldsymbol{p}}} \mathrm{d} V+\int\limits_{V} \nabla \cdot\left(\tau+\tau_{t}^{\mathrm{LES}}\right) \mathrm{d} V+\int\limits_{V} S_{\phi} \mathrm{d} V}\end{array} $

(2)

式(1)、(2)分别表示分别在RANS与LES下的控制体体积为V，在表面积为A的N-S方程; a表示流体表面A的矢量。v是流体运动速度，p是压力; τ表示粘性应力张量；S_ϕ表示源项。τ_t^RANS代表RANS下的雷诺应力，根据湍流模型τ_t^RANS=f(▽· v, k, ε)来确定；k是湍动能；ε是湍流耗散率。τ_t^LES代表LES下的亚格子应力，与亚格子模型τ_t^LES= $f(\nabla \cdot \tilde{\mathit{\boldsymbol{v}}}, \Delta)$有关，Δ表示过滤尺度，与网格大小有关。

DES结合了雷诺平均(RANS)与大涡模拟(LES)。在RANS中$\phi=\overline{\phi}+\phi_{t}^{\mathrm{RANS}}$，物理量ϕ可分为时间平均值ϕ与湍流扰动量ϕ_t^RANS2个部分。在LES中$\phi=\tilde{\phi}+\phi_{\mathrm{t}}^{\mathrm{LES}}$，物理量$\tilde{\phi}$为可求解的大尺度物理量，ϕ_t^LES需要亚格子应力模型求解。

1.2 势流区控制方程

在势流区不考虑粘性的作用，可将N-S方程简化为欧拉方程：

$ \frac{\partial \boldsymbol{v}}{\partial t}+(\boldsymbol{v} \cdot \nabla) \boldsymbol{v}=-\frac{1}{\rho} \nabla p-g \boldsymbol{k} $

对势流区域流体作无旋假设之后，势流区域基于层析水波理论的二维有限水深IGN模型给出规则波。二维有限水深IGN模型的基本方程如下^[15]。

假设海底条件z=-h(x)，自由面z=ζ(x), 自由面上压力为零，即$\hat{p}(x, t)=0$。二维情况下，速度场(u, v)与流函数ψ(x, z, t)之间的关系如下：

$ u(x, z, t)=\frac{\partial \psi(x, z, t)}{\partial z} $

(3)

$ w(x, z, t)=-\frac{\partial \psi(x, z, t)}{\partial x} $

(4)

令ψ(x, z, t)在海底的值为零，ψ(x, α, t)=0。对IGN模型的速度场作近似求解：

$ \psi(x, z, t)=\sum\limits_{m=1}^{K} \psi_{m}(x, t) f_{m}(\gamma)=\sum\limits_{m=1}^{K} \psi_{m}(x, t) \gamma^{2 m-1} $

其中γ=(z+h)/(ζ+h)。

关于瞬时波面ζ(x, t)和自由水面上的速度势$\hat{\varphi}(x, t)$需要满足：

$ \zeta_{t}+\hat{\psi}_{x}=0 $

(5)

$ \hat{\varphi}_{t}=-\frac{\partial}{\partial x} \frac{\partial T}{\partial \zeta_{x}}+\frac{\partial T}{\partial \zeta}-g \zeta $

(6)

式中：$\hat{\psi}(x, t)$表示自由面上的流函数；T表示动能，由下式给出：

$ T=\frac{1}{2} \sum\limits_{m=1}^{K} \sum\limits_{n=1}^{K}\left\{(\zeta+h) A_{\mathrm{mn}} \psi_{m x} \psi_{n x}-2\left[(\zeta+h)_{x} B_{\mathrm{mn}}^{1}-\right.\right.\\ \qquad h_{x} B_{\mathrm{mn}} ] \psi_{m x} \psi_{n}+\frac{1}{\theta}\left[\left(1+h_{x}^{2}\right) C_{\mathrm{mn}}-2\left(\zeta h_{x} C_{\mathrm{mn}}^{1}+\right.\right.\\ \qquad \left(\zeta+h_{x}\right)^{2} C_{\mathrm{mn}}^{2} ] \psi_{m} \psi_{n} \} $

(7)

其中:θ=ζ+h，A_mn=∫₀¹f_m(γ)f_n(γ)dγ，B_mn¹= ∫₀¹γf_m(γ)f′_n(γ)dγ，B_mn=∫₀¹f_m(γ)f′_n(γ)dγ，C_mn= ∫₀¹f′_m(γ)f′_n(γ)dγ，C_mn¹=∫₀¹γf′_m(γ)f′_n(γ)dγ，C_mn²= ∫₀¹γ²f′_m(γ)f′_n(γ)dγ。

为了使方程闭合，给出自由水面速度势$\hat{\varphi}(x, t)$和自由水面上流函数$\hat{\psi}(x, t)$的另一组关系式：

$ f_{m}(1) \hat{\varphi}_{x}=-\frac{\partial}{\partial x} \frac{\partial T}{\partial \psi_{m, x}}+\frac{\partial T}{\partial \psi_{m}} $

1.3 势粘流重叠区控制方程

在势粘流重叠区需要对势流和粘流进行耦合匹配，通过式(8)标量ϕ的输运方程来说明这一过程：

$ \int\limits_{V} \frac{\partial}{\partial t}(\rho \phi) \mathrm{d} V+\int\limits_{A} \rho \phi \mathit{\boldsymbol{v}} \mathrm{d} \mathit{\boldsymbol{a}}=\int\limits_{A} \varGamma \nabla \phi \cdot \mathrm{d} \mathit{\boldsymbol{a}}+\int\limits_{V} S_{\phi} \mathrm{d} V $

(8)

其中：

$ S_{\phi}^{*}=S_{\phi}-\mu(x) \rho\left(\phi-\phi_{w}\right) $

(9)

式中：Γ代表流体粘性系数；S_ϕ表示源项。当物理量ϕ是流场中的速度分量时，式(8)代表动量守恒方程；当ϕ是体积分数时，式(8)代表体积项守恒方程。在重叠区中，需要对式(8)右侧的源项进行修正，消除CFD区域内的反射。式(9)是修正后的源项S_ϕ^*，其中ϕ_w是已知量，可以根据IGN得到。μ(x)表示衰减函数，μ(x)随空间x位置变化可见图 2。以DB重叠区为例，重叠区边界D上μ(x)=0，而在DB其他区域则会随着重叠区域的不同位置而变化。式(10)给出了μ(x)在DB上的计算公式。

$ \mu(x)=\left\{\begin{array}{ll}{0, } & {x<x_{D}} \\ {\mu_{\max }\left(\sin ^{2} \frac{{\rm{ \mathsf{ π} }}\left(x-x_{D}\right)}{2 L_{\mu}^{2}}\right), } & {x_{D} \leqslant x \leqslant L_{\mu}} \\ {\mu_{\max }, } & {L_{\mu}<x \leqslant x_{B}}\end{array}\right. $

(10)

式中：μ_max表示最大衰减系数，根据波周期T得到μ_max=π/T。本文中过渡区DB的长度为2D，D代表圆柱直径。

2 数值计算 2.1 试验设置

本文对模型试验做1:50.3的全尺度模拟，试验浪高仪设置见图 3。圆柱半径r=8 m，吃水24 m，坐标原点在水线处的圆心。在圆柱周围每隔45°布置浪高仪。包括圆柱内侧浪高仪WPB1~WPB5和外侧WPO1~WPO5。

图 3中黑实线代表圆柱，外侧虚线圆代表外侧浪高仪位置所在圆。内侧浪高仪与圆柱表面距离为0.206 3 m，外侧浪高仪与圆柱表面的距离等于圆柱半径。本文研究来波为规则波，波浪条件见表 1。

2.2 数值计算域设置

数值模拟的计算域长度为2(L+2D)；为减小计算量，采用对称的思想只对一半圆柱进行数值模拟，计算域宽度大小为8D；计算域深度设为L。计算域边界条件设置见图 4，对称位置采用对称条件，计算域其他边界都采用速度入口的边界条件，圆柱表面选择壁面条件。

2.3 网格设置与收敛性分析

根据STAR-CCM+造波经验，在自由液面处水平方向网格大小设为L/80，垂直方向为H/20。对圆柱壁面进行网格加密，见图 5。

圆柱壁面的加密网格在水平以及垂直方向大小相等。圆柱壁面附近不同网格大小参数以及计算域网格总数见表 2。

根据圆柱水平方向力F_x以及垂直方向受力F_z对圆柱壁面网格参数做收敛性分析，见图 6。

由图 6可知，圆柱水平力F_x在网格1时收敛，垂向力F_z在网格2时收敛。即为了预报圆柱垂向受力，圆柱壁面在垂直方向的网格大小应与自由液面处保持一致。

3 结果讨论 3.1 波面及波浪力时历曲线验证

以T=9 s，波陡H/L=1/16为例，波面以及圆柱水平、垂向波浪力的数值模拟时历曲线结果与文献中的试验值EXP、势流DIFFRACT以及粘流OpenFOAM数值结果进行对比验证。圆柱内侧浪高仪WPB2~WPB4的时历结果见图 7，外侧浪高仪WPO2~WPO4的波面时历曲线见图 8，圆柱波浪力时历曲线验证见图 9。

由图可知势流DIFFRACT求解器可以准确预报波面爬升幅值，但对于高阶项的捕捉误差较大；采用粘流OpenFOAM求解器可以准确预报波面爬升的非线性效应, 文献[13]中使用了8.29×10⁶网格。采用势粘流单向耦合技术，网格数为1.14×10⁶，并且对WPB2~WPB4和WPO2~WPO4处的波面以及圆柱波浪力数值预报结果与试验值吻合良好。

3.2 波浪绕射谐波分析

以T=9 s为例比较不同非线性下波浪爬升的绕射效应。迎浪位置WPB1处波浪爬升的高度最大^[16]。比较不同波陡下WPB1位置处在一个周期内的波面时历, 见图 10。

由图 10可知，入射波非线性越强，波浪爬升高度越大。在波陡H/L=1/10时，WPB1处的最大爬升高度可达2.0 A。

波面信号的一阶、二阶无因次化计算公式分别为η⁽¹⁾/A、η⁽²⁾r/A²。其中η⁽¹⁾、η⁽²⁾分别代表波面谐波分析后的一阶、二阶值。不同的非线性入射波下WPB1~WPB5处波面一阶、二阶谐波结果的空间分布分别可见图 11(a)、11(b)。

由图 11可知，波面谐波大小与入射波非线性大小有关。迎浪位置WPB1、WPB2处受波浪非线性影响较大。随着入射波非线性增强，一阶谐波值会变大。由于波面的二阶值与波幅的平方相比很小，二阶谐波值会随着非线性增强而变小。势流结果适用于弱非线性(H/L=1/30)波浪爬升的谐波预报。OpenFOAM和本文势粘流耦合的数值结果与试验值吻合良好，可以预报强非线性波浪爬升的谐波预报。

保持波陡H/L=1/16不变，比较在不同入射波周期时圆柱周围波面的一、二阶波面谐波值，分别见图 12和13。

由图 12可知，在相同波陡，不同波长的入射波作用下，势流和粘流对波面一阶谐波的预报结果与试验值吻合良好。由图 13知，由于圆柱对短波反射作用较强，短波下波面二阶谐波值较大。在长波T=15 s时的二阶波面谐波值接近于零。

3.3 圆柱受力分析

以T=9 s为例比较不同波陡(H/L=1/16)、不同周期(T=9 s)下圆柱波浪力，分别见图 14和图 15。

由图 14可知，相同来波周期下，波陡越大圆柱受力越大。由图 15可知，在相同的波陡下，周期大的入射波波高大，圆柱受力越大。并且F_z的平均值会随着波高变大向下偏移，对圆柱有向下“拍击”的作用。

3.4 效率比较以及圆柱绕射流场分析

本文采用的势粘流耦合方法模拟圆柱体波浪爬升问题，相比全粘流方法可以减小计算所需网格数。以T=9 s，波陡H/L=1/16为例，在相同的网格设置下，本文方法所需网格数为1.14×10⁶，而全粘流方法由于要增加计算域的长度，需要的网格数约为1.41×10⁶。

以波浪周期T=9 s，波陡H/L=1/16为例，每隔1/4周期研究圆柱周围波浪爬升过程中波面的绕射变化，见图 16。当t=T₀时, 波峰传到圆柱上游位置，由于壁面的反射作用，形成的向上游传播的圆弧形反射波在与上游入射波波峰叠加，导致圆柱上游位置出现波面爬升，见图 16(a)。当t=T₀+T/4时, 波峰到达圆柱位置，圆柱近壁面WPB1处达到最大的波面爬升高度。除了向上游传播的圆弧形反射波外，圆柱两侧会出现对称的相对扇形“凹陷”，形成向下游传播的散射波，见图 16(b)。当t=T₀+T/2时，波峰到达圆柱下游位置，入射波与向下游传播的散射波叠加，在圆柱下游位置形成较大波面爬升，见图 16(c)。当t=T₀+3T/4时，波谷到达圆柱位置，又开始形成向上游传播的圆弧形反射波，见图 16(d)。

4 结论

1) 势流方法适用于弱非线性工况(H/L=1/30)下的波浪爬升预报。势粘流耦合方法以及OpenFOAM的数值预报结果与试验吻合良好，但OpenFOAM计算需要更多的网格数。

2) 入射波的非线性越强，圆柱壁面的波浪爬升就越明显。迎浪位置WPB1、WPB2处的波面谐波值受波浪非线性影响较大。由于圆柱的反射效应，短波下波面二阶谐波值较大。而长波T=15 s时的二阶波面谐波值接近于零。

3) 圆柱收到的波浪力与来波的波高有关，波高越大，受力越大。并且随着垂向力的平均值会随着波高变大向下偏移，对圆柱有向下“拍击”的作用。

4) 在相同的网格设置下，相对于本文势粘流耦合方法，全粘流方法计算所需要的网格更多。

本文的数值模拟设置和网格参数可作为波浪爬升数值模拟的参考。后续可在此基础上对多个圆柱的波浪爬升问题进行研究。