畸形波(又称凶波或大波)是指海上偶然出现的具有超大波高的波浪，其概念由Draper^[1]首次提出，几十年来不断有畸形波造成船舶失事和海洋结构物破坏的报道。畸形波出现在除北冰洋以外的三大洋中的不同海域，北海^[2-3]、南非东南海域^[4]和日本海域^[5]等都有观测到畸形波的记录。畸形波的形式与产生条件多样，其可发生在不同水深(深水、有限水深和浅水)，可伴随显著海流或无明显海流，可能为孤立的大波或一组大波，但其共同特征为超高的波高(强非线性)、出现的偶然性及存在的短暂性。目前常从波高角度定义畸形波，常用为Kjeldsen^[6]及Klinting等^[2]的定义。Kjeldsen定义畸形波为波高超过有意波高2倍的波浪。Klinting等在Kjeldsen定义的基础上考虑了畸形波相邻波浪的影响。

畸形波结构物相互作用研究多在实验室或数值水池中开展^[7-12]。由于研究建立在畸形波机理及模拟的基础上，故相关研究仍不成熟。其中数值方法多基于商业势流软件或CFD求解器^[13-14]，目标畸形波多基于线性叠加模型。前者可方便快捷地得到预报结果，然而往往无法准确考虑畸形波的强非线性及期间伴随的上浪砰击效应。后者可模拟上浪砰击及过程中的湍流效应，但是计算耗时，且多数研究不考虑结构的变形(即流固耦合)。无论采用势流方法亦或CFD方法，正确模拟强非线性畸形波都是前提也是难点。在畸形波上浪砰击过程中，巨大的砰击压力使局部结构在很短时间内改变运动状态，进而影响流场的状态和砰击压力。因此的瞬态砰击作用过程中有必要考虑流固耦合效应。

本文构建求解N-S方程的二维数值波浪水槽，基于速度入口和守恒消波法模拟某北海实测畸形波^[3]。采用一种基于SIMPLE法的隐式算法求解畸形波对垂直弹性板的上浪砰击问题，并分析不同条件对结构响应的影响。

1 数值波浪水槽

数值水槽利用交错网格离散求解域，采用有限差分法求解N-S方程，使用部分单元参数处理复杂流场边界，采用体积分数(volume of fluid, VOF)法重构自由液面。数值波浪水槽的控制方程为

连续性方程:

$ \nabla \left( {\lambda \mathit{\boldsymbol{u}}} \right) = 0 $

(1)

动量方程:

$ \frac{{\partial \mathit{\boldsymbol{u}}}}{{\partial t}} + \nabla \mathit{\boldsymbol{uu}} = - \frac{1}{\rho }\nabla p + \mathit{\boldsymbol{f + }}\frac{\gamma }{\rho }{\nabla ^2}\mathit{\boldsymbol{u}} $

(2)

体积分数输运方程:

$ \frac{{\partial \lambda F}}{{\partial t}} + \nabla \left( {\lambda \mathit{\boldsymbol{u}}F} \right) = 0 $

(3)

式中：u为流场速度，p为流场压力，ρ为流体密度，γ为动力粘性系数，λ为部分单元参数，f为体积力，F为体积分数，▽为散度(梯度)算子。利用SIMPLE算法求解速度压力耦合方程，使用VOF-Youngs法^[15]重构自由液面，压力梯度项和速度耗散项采用中心差分法离散，通过一阶迎风差分与中心差分的混合格式离散扩散项，其中迎风差分的权重设为0.4。忽略表面张力，则自由液面条件为

$ \frac{{\partial {\mathit{\boldsymbol{u}}_n}}}{{\partial \tau }} + \frac{{\partial {\mathit{\boldsymbol{u}}_\tau }}}{{\partial n}} = 0 $

(4)

$ - p + 2\mu \frac{{\partial {\mathit{\boldsymbol{u}}_n}}}{{\partial n}} = - {p_0} $

(5)

式中：u_n和u_τ分别为自由液面处速度矢量的法向和切向分量，p₀为自由液面上方大气压力。水底采用无滑移条件。采用速度入口方式造波。造波边界的速度、压力和自由液面抬升根据相应模型给出的理论值设定。

消波方法的好坏直接影响数值水槽模拟波浪的质量。Mayer等^[16]提出了松弛消波法，其优点为无需改变控制方程，可以吸收多种类型的波浪。松弛法的表达式为

$ \left\{ \begin{array}{l} {u_1} = C{u_0}\\ {v_1} = C{v_0}\\ {\eta _1} = C{\eta _0} \end{array} \right. $

(6)

式中：下标0表示消波前的物理量，下标1为消波后的物理量，C为空间分布的消波权重且仅与水平位置有关。在消波区始端C取1，末端C取0。实际应用中C的函数形式多种多样，本文取

$ C = 1 - \exp \left[ { - 5{{\left( {1 - {x^ * }} \right)}^n}} \right] $

(7)

式中：x^*=(x-x_a)/(x_b-x_a)，x_a和x_b分别为消波区始端和末端的水平坐标位置。Hu等^[17]提出了一种守恒消波法以适应变化的流场状态，提高消波效率，其中参数C表示如下

$ \int_{{x_{\rm{a}}}}^{{x_{\rm{b}}}} {\left[ {CJ + \left( {C - 1} \right)\eta } \right]{\rm{d}}x} = 0 $

(8)

式中：$J = \Delta t\int_0^{{Y_{\rm{H}}}} {\frac{{\partial (\lambda u)}}{{\partial x}}} {\rm{d}}y$，Y_H为水底到水槽顶端的垂向距离(即水槽高度)。

2 畸形波模拟

采用波列叠加模型模拟类似北海实测的确定性畸形波。波列叠加模型将畸形波视为随机波列和汇聚波浪按一定能量比例的叠加，其表达式为

$ \begin{array}{*{20}{c}} {\eta = \sum\limits_{i = 1}^{{N_1}} {{a_{1i}}\cos \left( {{k_{1i}}x - {\omega _{1i}}t + {\theta _{1i}}} \right)} + }\\ {\sum\limits_{i = 1}^{{N_2}} {{a_{2i}}\cos \left[ {{k_{2i}}\left( {x - {x_c}} \right) - {\omega _{2i}}\left( {t - {t_c}} \right)} \right]} } \end{array} $

(9)

式中：第一个求和项为随机波列，第二个求和项为汇聚波列；下标1和2分别代表随机波列和汇聚波列，a_1i和a_2ii个组成波的波幅；k_i、ω_i和θ_ii个组成波的波数、圆频率和随机相位；常数x_c和t_c分别为汇聚位置和汇聚时刻；N₁和N₂为相应波列的组成波数目。组成波波幅常通过海浪谱(如P-M谱)按比例分配得到，即

$ \left\{ \begin{array}{l} {a_{1i}} = \sqrt {2{p_1}S\left( {{\omega _{1i}}} \right)\Delta {\omega _1}} \\ {a_{2i}} = \sqrt {2{p_2}S\left( {{\omega _{2i}}} \right)\Delta {\omega _2}} \end{array} \right. $

(10)

式中辽阔p₁、p₂影响了畸形波波高及模拟效率。本文的目标谱为P-M谱，表达式为S(ω)=Aexp(-B/ω⁴)/ω⁵，其中A和B为常数。目标谱的选取参考北海实测畸形波^[3]，观测位置水深130 m，有意波高为5.65 m，畸形波波高为18.04 m，波峰高度13.9 m。能量比例取为p₁=0.9与p₂=0.1。汇聚波列组成波数目为N₁=250与N₂=500。数值水槽长680 m，高150 m。畸形波最大波高出现位置距造波端120 m，出现时刻t=432 s。图 1为北海实测畸形波液面时历曲线，图 2为数值水槽模拟得到的液面时历，可见两者在畸形波出现时刻较接近。图 3对比了实测谱和数值模拟得到的波浪谱，两者吻合较好。表 1对比了数值模拟和观测数据得到的畸形波参数，两者吻合良好。表 1中，α=H_j/H_s、β₁=H_j/H_j-1，β₂=H_j/H_j+1、μ=η_j/H_j，其中H_j、H_j-1和H_j+1分别为畸形波，畸形波前方相邻波及畸形波后方相邻波的波高，η_j为畸形波波峰液面抬升。以上对比表明模拟生成的畸形波是正确合理的，可用于进一步研究上浪砰击。图 4给出了结构前端的畸形波轮廓，可见模拟得到的畸形波前端陡峭，后端平坦，最大波峰处出现波浪翻卷，具有强非线性。

3 流固耦合算法求解

通过向数值波浪水池求解器中添加结构求解模块和耦合处理算法实现流固耦合模拟。结构空间离散采用有限元法。由于在二维数值波浪水槽中平板可简化为板条梁模型，因此采用Euler梁单元进行结构离散。每个单元具有两个节点，每个节点具有2个自由度，即线位移u_i和角位移θ_i。欧拉梁满足∂u_i/∂y=θ_i(其中y轴垂直于梁截面)，因此单元具有三阶精度。动力分析中采用集中质量矩阵和比例阻尼矩阵，结构控制方程为

$ \mathit{\boldsymbol{M\ddot U}} + \mathit{\boldsymbol{C\dot U}} + \mathit{\boldsymbol{KU}} = \mathit{\boldsymbol{F}} $

(11)

式中：U为节点位移，$\boldsymbol{\dot U}$和$\boldsymbol{\ddot U}$分别为节点速度和节点加速度，M、C和K分别为整体质量、阻尼和刚度矩阵，F为节点力。比例阻尼矩阵为C=P₁M+P₂K，其中，P₁=2×10^-5，P₂=8×10^-5。

采用Houbolt格式求解方程(11)。Houbolt是一种无条件稳定格式，相应加速度和速度的离散格式如下

$ {{\mathit{\boldsymbol{\ddot U}}}_{t + \Delta t}} = \frac{1}{{\Delta {t^2}}}\left( {2{\mathit{\boldsymbol{U}}_{t + \Delta t}} - 5{\mathit{\boldsymbol{U}}_t} + 4{\mathit{\boldsymbol{U}}_{t + \Delta t}} - {\mathit{\boldsymbol{U}}_{t - \Delta t}}} \right) $

(12)

$ {{\mathit{\boldsymbol{\dot U}}}_{t + \Delta t}} = \frac{1}{{6\Delta t}}\left( {11{\mathit{\boldsymbol{U}}_{t + \Delta t}} - 18{\mathit{\boldsymbol{U}}_t} + 9{\mathit{\boldsymbol{U}}_{t - \Delta t}} - 2{\mathit{\boldsymbol{U}}_{t - 2\Delta t}}} \right) $

(13)

将式(12)、(13)代入式(11)，得到离散的结构方程为

$ \begin{array}{*{20}{c}} {\left( {2{a_1}\mathit{\boldsymbol{M}} + 11{a_2}\mathit{\boldsymbol{C}} + \mathit{\boldsymbol{K}}} \right){\mathit{\boldsymbol{U}}_{t + \Delta t}} = }\\ {{\mathit{\boldsymbol{F}}_{t + \Delta t}} + \left( {5{a_1}\mathit{\boldsymbol{M}} + 18{a_2}\mathit{\boldsymbol{C}}} \right){\mathit{\boldsymbol{U}}_t} - }\\ {\left( {4{a_1}\mathit{\boldsymbol{M}} + 9{a_2}\mathit{\boldsymbol{C}}} \right){\mathit{\boldsymbol{U}}_{t - \Delta t}} + \left( {{a_1}\mathit{\boldsymbol{M}} + 2{a_2}\mathit{\boldsymbol{C}}} \right){\mathit{\boldsymbol{U}}_{t - 2\Delta t}}} \end{array} $

(14)

其中，a₁=1/Δt²，a₂=1/6Δt。已知F_t+Δt则求解方程(14)可得t+Δt该时刻节点位移，进而利用式(12)、(13)得到节点加速度和节点速度。同时，作为流场的边界条件，结构的运动也会改变流场的状态，即U_t+Δt和${\boldsymbol{\dot U}_{t + \Delta t}}$，进而会影响到节点力。这说明式(11)或(14)与流场求解是耦合的。本文的求解策略为：1)考虑流场压力对节点力F_t+Δt的影响，求解相应的节点位移、速度与加速度；2)考虑节点位移与速度对流场边界条件的影响，计算节点力；3)考虑策略1)和2)的相互作用进行耦合迭代求解。Hu等^[18]讨论了该流固耦合算法的推导、验证与收敛性分析等，这里仅给出流固耦合求解器的流程。

如图 5所示，“数据输入”模块读取计算必须的流体与结构信息，以及计算控制参数。其中流体运动粘性系数为0.000 001 145 m²/s，流体密度为1 000 kg/m³，重力加速度为9.8 m/s²，板长为22 m，弹性模量为2.1×10¹¹ Pa，Poisson比为0.3，结构密度为7 860 kg/m³，横截面积为0.2 m²，横截面惯性矩为0.000 666 7 m⁴。“初始化”模块划分计算网格并设置网格属性。“确定时步”模块根据稳定性条件计算最大允许时步。“设置造波边界”模块基于式(9)设置造波端的边界条件实现造波。在“临时速度场计算”完成后进行n+1时步的压力场迭代(“压力场计算”模块)。根据更新后的压力场进行“速度场矫正”，同时将压力场转化为结构节点力进行“结构速度计算”。若流场速度与结构速度在耦合面处相等，则迭代收敛，否则继续调整压力场。“自由面速度赋值”模块根据收敛后的速度场设置自由液面边界条件，接着采用施主-受主网格技术进行“流体体积计算”。根据新时步的流体体积进行“自由液面更新”。重复上面过程直至计算结束。

4 上浪砰击计算结果

模拟总网格数约50 000，并在水平与垂直板附近进行局部加密。模拟时步应足够小以捕捉到结构振动的不同模态和脉冲压力，这里取为0.000 5 s。垂直弹性板(简化为板条梁模型)等分为44个Euler单元，并从底至顶依次编号。模拟过程中，畸形波从甲板左侧涌入，向右传播并最终砰击垂直板。整个上浪砰击过程持续时间约10 s，不考虑相邻波浪的连续上浪砰击。

图 6为上浪砰击过程中不同阶段的液面形态。从图 6中可以看出，在到达垂直板之前，甲板上的波浪已经溃塌并伴随破碎现象，高速水体到达垂直板时对板底部造成较大的砰击载荷。由于板的阻碍作用，上浪水体开始沿垂直板爬升，此时，由于流体的反作用力，垂直板底端前部的水平甲板受到较大局部载荷。当爬升至最高点后，流体开始下落，并从甲板左端流入“大海”，此过程伴随着大量的波浪翻卷、破碎以及气泡现象。

为研究不同边界条件下的垂直板动力响应，这里选择三种支撑形式：两端简支、两端固支和底端固支顶端自由。图 7为以上三种支撑形式的板节点位移，其中节点12、23和34分别距底部6、11.5和17 m，图中正向表示与上浪水体传播方向相同，波浪砰击开始于约302 s。可以看出，垂直板的节点位移曲线并不规则，且包括多阶模态的振动，这与上浪砰击现象的复杂性与瞬态冲击载荷有关。为研究激发的振动模态，图 8给出了三种支撑形式板节点位移的FFT曲线，注意这里选择节点4以便捕捉到振动的所有模态。

表 2为不同支撑形式板的干频率以及通过FFT变换得到的湿频率。从表 2中可以看出砰击载荷激发了三阶模态。在图 8中可发现某频率附近出现双峰值，在表 3中体现为某干频率附近有2个湿频率与之对应。双峰值的产生是由于在上浪砰击的过程中结构呈现两种运动。以图 7(a)为例，在437~441 s时间范围内，由于流体与结构接触充分，附加质量效应明显，结构振动周期较长；而441 s后由于流体下落回流，流体与结构接触面积变小，附加质量效应减弱，结构振动频率相应增大，此时结构振动频率一般与干频率接近。对于不同支撑形式，其FFT变换峰值对应频率不同，但规律大体相似。

表 3给出了三种支撑形式下不同高度位置处的最大砰击压力，并与刚性板进行对比。可以看出，考虑结构水弹性效应后，流场的砰击压力有所减小。在垂直板底部，简支边界对最大砰击压力的减小幅度较大，固支与悬臂边界对于砰击压力的减少效果小于简支边界，这是由于固支和悬臂边界在根部处的约束更强，更趋近于刚性边界。同时可以发现，固支和悬臂边界的垂直板底部砰击压力相同，这是由于两种支撑在底部接近，而顶端约束区别对于底部局部砰击压力的影响可以忽略。

5 结论

1) 上浪砰击整体过程中垂直板振动频谱有双峰值现象发生。

2) 由于受到瞬态的砰击载荷，结构的动力响应中往往包含高阶模态。

3) 上浪后期水弹性附加质量效应影响较小，结构振动频率与干频率相当。

4) 受到水弹性效应的影响，与刚性结构相比，弹性结构表面的局部砰击压力有所降低，且结构的刚度越低，砰击压力的减小约明显。