随着全球气候变暖，海冰冰厚变薄使得北极资源的开发以及航线开通成为可能。许多国家竞相储备开发北极资源和北极航道，准备随时进驻北极。因此，合理认识和利用冰缘区内波-冰作用规律，有利于加强对冰区内波浪传播特性的了解和提高极地资源利用率具有重要意义。对于浮冰与波浪的相互作用，Harms^[1-2]通过试验回归出恒密度厚度的二维浮冰漂移运动估算公式。当浮冰的横摇周期略大于波浪周期时，浮冰具有最大漂移速度。Lever等^[3]在水池试验中研究了冰山自身形状对其在波浪中运动的影响。Huang等^[4]研究发现浮冰漂移速度随着横摇周期与波浪周期的比值增大而增大。Bennetts等^[5]在试验中发现浮冰的存在阻碍了波浪传播，且随着波陡的增大，浮冰对波浪的阻碍作用增强。Yiew等^[6]发现当波陡较大时，波浪冲洗现象对浮冰的运动响应有明显抑制作用。郭春雨等^[7]试验研究了浮冰在波浪中的纵向运动，同时发现了波浪冲洗现象。波浪在冰区的传播比例随着波长和波频率的减小而增大，且波浪传播比例增大也增强了其对浮冰的破坏能力^[8]。

由于模型试验具有成本过高、存在设备干扰的缺点，数值模拟方法被越来越多的用于研究波-冰作用。其中，势流分析方法因其成本低、计算速度快等优点受到了学者的青睐。Shen等^[9]通过理论研究了反射波、阻力系数和附加质量影响浮体漂移运动的的因素，如波的反射、附加质量等。Meylan^[10]研究发现浮冰运动很大程度上依赖于浮冰形状，波浪散射主要取决于入射波方向。Bennetts等^[11]研究了非均质浮冰对波浪的散射。Smith等^[12]对均质浮冰对波浪散射现象求解进行了研究，并在此方法的基础上发展了非均质浮冰对波浪散射的求解。Meylan等^[13]通过改进莫里森方程研究大波长规则波中浮冰的纵荡运动，与试验结果具有较好的一致性。Toffoli等^[14]提出了波浪海冰相互作用的试验模型来验证冰缘区波浪衰减的理论模型。由于波-冰作用中固有的强非线性现象无法通过势流理论精确计算出，因此计算流体力学动态模拟技术的迅速发展研究波-冰作用起到了重要的作用。Bai等^[15]研究表明粘流软件比势流软件更适合于模拟浮冰运动。Huang等^[16]通过粘流软件模拟了波浪冲洗冰体表面和波浪散射现象。

在浮冰与波浪的相互作用中，浮冰在波浪作用下会发生非形变运动。与此同时，波浪场也会由于浮冰的运动而发生改变。目前关于波-冰作用的研究对于相互影响已经有了充分了解，但对波浪冲洗现象对浮冰运动的影响研究较少。波浪冲洗冰体表面现象虽有分析，但是很少对波浪冲洗现象、浮冰运动以及波浪参数之间定性定量关系的研究。本文通过粘流模型数值模拟了波浪中浮冰的运动，利用频谱分析研究了浮冰运动响应。分析了波浪冲洗现象对浮冰运动的影响，对比在有/无波浪冲洗现象工况下的浮冰运动。并研究了不同波浪参数下波浪冲洗现象对浮冰垂荡、纵摇运动的影响。

1 波-冰作用数值计算方法 1.1 流体运动控制方程

浮冰区域流体运动满足连续方程和动量方程。考虑到湍流流动，通常采用时均法将流体运动控制方程中的各物理量分解为2部分，分别为时间平均值和相对于平均值的脉动值，得到雷诺时均方程组，其表达式为：

$ \nabla \cdot \bar{\boldsymbol{v}}=0 $

(1)

$ \frac{\partial \rho \bar{\boldsymbol{v}}}{\partial t}+\nabla \cdot(\rho \overline{\boldsymbol{v} \boldsymbol{v}})=-\nabla \bar{p}+\nabla \cdot\left(\overline{{\tau}}-\rho \overline{\boldsymbol{v}^{\prime} \boldsymbol{v}}\ ^{\prime}\right)+\rho g $

(2)

式中：v是时均速度；v′是时均速度脉动值；τ是粘性应力；p是压力；$\rho \overline{\boldsymbol{v}^{\prime} \boldsymbol{v}^{\prime}} $是雷诺应力，即湍流涨落引起的时均效应，选用RKE 2L湍流模型^[17]来对控制方程进行封闭求解，该模型对旋转流动以及流动分离具有较好的模拟效果，其表达式为：

$ \begin{gathered} \frac{\partial}{\partial t}(\rho k)+\nabla \cdot(\rho k \bar{\boldsymbol{v}})=\nabla \cdot\left[\left(\mu+\frac{\mu_{t}}{\sigma_{k}}\right) \nabla k\right]+ \\ P_{k}-\rho\left(\varepsilon-\varepsilon_{0}\right)+S_{k} \end{gathered} $

(3)

$ \begin{gathered} \frac{\partial}{\partial t}(\rho \varepsilon)+\nabla \cdot(\rho \varepsilon \bar{\boldsymbol{v}})=\nabla \cdot\left[\left(\mu+\frac{\mu_{t}}{\sigma_{\varepsilon}}\right) \nabla \varepsilon\right]+ \\ \frac{1}{T_{e}} C_{\varepsilon 1} P_{\varepsilon}-C_{\varepsilon 2} f_{2} \rho\left(\frac{\varepsilon}{T_{e}}-\frac{\varepsilon_{0}}{T_{0}}\right)+S_{e} \end{gathered} $

(4)

式中：μ为动力粘度；σ_k、σ_ε、C_ε1、C_ε2为模型系数，σ_k值为1，σ_ε值为1.2;C_ε1取为0.43与$\frac{\eta}{\eta+5} $中的最大值($\eta=\frac{S k}{\varepsilon}, S=\left|\frac{1}{2}\left(\nabla \overline{\boldsymbol{v}}+\nabla \overline{\boldsymbol{v}}^{\mathrm{T}}\right)\right| $); C_ε2限值为1.9；P_k、P_ε为根据RKE 2L湍流模型变体的结果项；f₂为阻尼函数；S_k、S_ε为用户指定的源项。

1.2 数值造波理论

本文通过速度入口造波方式生成波浪，通过在入口处和出口处分别通过松弛区域消波技术和阻尼消波方法实现消波从而减少壁面反射。本文选取斯托克斯五阶波浪^[18-19]，使模拟的波浪更符合实际的规则波。在造波边界处的速度和波面瞬时升高满足以下条件：

x方向速度：

$ u_{x}=c \sum\limits_{n=1}^{5} n \lambda_{n} \cosh [n k(z+d)] \cdot \cos [n(k x-\omega t)] $

(5)

z方向速度：

$ u_{z}=c \sum\limits_{n=1}^{5} n \lambda_{n} \sinh [n k(z+d)] \cdot \sin [n(k x-\omega t)] $

(6)

波面瞬时升高：

$ \eta=\frac{1}{k} \sum\limits_{n=1}^{5} \lambda_{n} \cos [n(k x-\omega t)] $

(7)

其中各系数分别为：

$ \left\{\begin{array}{l} \lambda_{1}=\lambda \\ \lambda_{2}=\lambda^{2} B_{22}+\lambda^{4} B_{24} \\ \lambda_{3}=\lambda^{3} B_{33}+\lambda^{5} B_{35} \\ \lambda_{4}=\lambda^{4} B_{44} \\ \lambda_{5}=\lambda^{5} \mathrm{~B}_{55} \end{array}\right. $

式中：ω、d、k分别是圆频率、水深和波数；c=cosh(kd)；λ为常系数，其余各项系数参考文献[20]。

入口消波区的消波原理是通过在动量方程中添加源项q_ϕ强迫波浪在速度入口附近转换成入射波：

$ q_{\phi}=-\gamma \rho\left(\phi-\phi^{*}\right) $

(8)

式中：γ为强迫系数；ϕ为数值模拟中的动量输运方程的结果；ϕ^*为原动量方程的理论解。冰水相互作用的数值水池模型如图 1所示。

出口消波区的原理是通过在Z轴方向上添加阻尼完成消波作用，其中阻尼项为：

$ \left\{\begin{array}{l} S_{z}^{d}=\rho\left(f_{1}+f_{2}|w|\right) \frac{\mathrm{e}^{\kappa}-1}{\mathrm{e}-1} w \\ \kappa=\left(\frac{x-x_{\mathrm{sd}}}{x_{\mathrm{ed}}-x_{\mathrm{sd}}}\right)^{n_{d}} \end{array}\right. $

(9)

式中：x_sd为起始消波位置；x_ed为结束消波位置；f₁、f₂、n_d为阻尼消波的参数；w为Z方向的速度。

1.3 浮冰运动方程

在随浮冰平移和旋转的运动坐标系中，浮冰的运动方程为：

$ m\left(\frac{\mathrm{d} \boldsymbol{V}}{\mathrm{d} t}+\boldsymbol{\varOmega} \times \boldsymbol{V}\right)=\boldsymbol{f} $

(10)

$ \boldsymbol{I} \cdot \frac{\mathrm{d} \boldsymbol{\varOmega}}{\mathrm{d} t}+\boldsymbol{\varOmega} \times(\boldsymbol{I} \cdot \boldsymbol{\varOmega})=\boldsymbol{M} $

(11)

式中：m是浮冰质量；V是浮冰运动速度；Ω是浮冰旋转角速度；I是浮冰转动惯量；M是浮冰受到的合外力矩；f是浮冰受到的的合外力。

2 波-冰作用计算模型设置 2.1 波浪与计算模型验证 2.1.1 波浪验证

根据极地冰缘区出现的波浪高度和频率分布规律和波浪参数与浮冰特征尺寸的比例关系，在相同的波浪高度下，对浮冰在一系列波长(波陡偏大)波浪中的运动进行了数字模拟和分析研究，工况1~5波高恒定为0.04 m，并且入射波长为0.6、0.7、0.8、0.9、1.0 m。

液面波高探测点位置选取为圆柱所在位置，其中选取工况5的波浪来记录探测点处液面升高过程，将斯托克斯五阶波理论解与其监测到的波高时历曲线数据进行比较。对于波陡小的波浪，造波效果与其设置的波高H方向的网格尺寸Δz有重要的关联，而能否生成目标波浪决定了计算结果的准确与否。因此需要再进行收敛性分析，这里分别选取3组不同大小的自由液面处波高方向网格和3组不同时间步进行生成波浪进行分析，Δz分别为H/5、H/10及H/15。波长方向的网格尺寸Δx是波高方向网格尺寸Δz的4倍。网格与时间步设置如表 1所示。

如图 2所示，数值水池生成的波浪在频率、波高等波浪参数在网格B上的误差最小，与斯托克斯五阶波理论解吻合较好。因此，波高方向上网格尺寸Δz为H/10被选取进行数值模拟。将该套网格在3种不同时间步下进行了计算，如图 3所示，数值模拟时间步Δt为0.002 s的工况的结果收敛。因此，出于效率和资源上的综合考虑，波高方向网格尺寸Δz为H/10，时间步长Δt为0.002 s被选入进行后面的数值模拟。

2.1.2 模型验证

本文选取文献[6]中单个浮冰与规则波相互作用试验研究，提取试验中浮冰的垂荡幅值时历曲线来验证数值模拟计算结果。浮冰模型为圆柱体，半径为0.2 m，厚度为0.015 m，密度为636 kg/m³，浮冰边缘附有薄壁结构，薄壁结构高为0.05 m，厚度为0.025 m，在浮冰周围施加弹簧力以防止浮冰产生大幅度漂移。在保证无池壁反射的前提下，适当缩减了水池的长度和宽度，数值模型中水池宽2 m，长7 m，水深0.83 m。入射波频率为1.25 Hz，波幅为0.085 m。浮冰中心距离造波入口4 m处。计算模型中网格划分遵循上述网格参数设置，在浮冰周围加密以保证计算精度，网格划分图如图 4所示，其中对于k-ε模型需要满足30≤y⁺≤300，一般以接近30为佳，因此y⁺取值为30，第1层网格厚度为0.000 47 m。

由图 5可知，试验值与数值模拟中的浮冰垂荡运动基本一致，浮冰垂荡运动峰值以及到达相邻峰值所需的周期大致吻合，由于偏差较小，数值模拟结果相较于实验结果没有明显偏高或偏低趋势，因此可以将偏差视为数值计算或实验测量精度的判断依据，说明本文采用的数值模拟方法是合理和准确的。

2.2 三维模型和网格划分

如图 6所示，数值水池尺寸为4.5 m×2 m×0.8 m，水深0.5 m。数值水池的前后边界设置为对称边界条件；左侧边界、上侧边界和下侧边界设置为速度入口边界条件；右侧边界设置为压力出口边界条件。入口以及出口消波区区域大小均取为一个波长。

将浮冰视为刚体，由于实际海况上浮冰互相之间经常会受到碰撞和波浪的冲洗，浮冰边角的磨损使其形状近似于圆柱，如图 7所示。将浮冰模型进行处理便于计算，故简化为圆柱体和带一定高度薄壁的圆柱体2种形状，浮冰中心与入口处的距离是1.5 m。其中圆柱参数：直径0.5 m；厚度0.05 m；吃水0.03 m；密度600 kg/m³。为对无波浪冲洗现象进行研究，将圆柱浮冰周围增设高为0.05 m薄壁，研究表明薄壁对浮冰运动的影响可忽略不计^[21]。

本文通过重叠网格技术模拟浮冰运动，在浮体运动区域以及自由液面区域进行网格加密，如图 8为重叠网格。

3 波-冰作用计算结果及分析

波浪的参数对浮冰运动响应有着关键的影响。为了进一步研究浮冰与入射波的相互作用，通过频谱分析研究了不同波长下浮冰垂荡η₃、纵摇η₅运动响应。

3.1 波浪冲洗现象

浮冰在波浪作用下会发生大幅的摇荡运动，海水冲洗浮冰表面的现象频频出现。由数值模拟可以看出，波长增加，波浪冲洗现象会随之减弱。图 9为在不同波长条件下一个周期的浮冰表面的上浪现象。图 9(a)是波长为0.6 m时波浪冲洗状态，浮冰表面上浪严重；图 9(b)是波长为1.0 m时的波浪冲洗状态，浮冰表面未完全被水体覆盖，与图 9(a)相比，其所受的波浪冲洗现象较弱。

图 10为λ=0.6 m时的浮冰摇荡运动频谱分析图。在相同的波浪激励频率作用下，没有受到波浪冲洗的浮冰在运动幅值上要大于受到波浪冲洗的浮冰。产生上述差异的主要原因是在波浪与浮冰相互作用中存在波浪冲洗冰体现象，一方面堆积在冰体表面的水阻碍浮冰的运动；另一方面冰体表面的水体会产生碰撞继而引起水体破碎，使波浪产生一定程度的衰减，使得浮冰接受的波浪能量较少。

图 11为不同波长下的波浪冲洗现象对浮冰的垂荡和纵摇的影响。浮冰的垂荡与纵摇随着波长减小而减小。当波长较短时，波陡较大，波浪冲洗冰体表面现象较严重。由于堆积在冰体表面的水体较多以及波浪能量衰减较大，一定程度上阻碍了浮冰运动，使其运动幅度减弱。

由此可见，在给定波高情况下，波长越小，即波陡越大，波浪冲洗现象越严重。浮冰运动在波浪冲洗现象作用下受到了一定的抑制，波浪冲洗现象越强，其对浮冰的抑制作用越强。

3.2 波浪散射现象

通过模拟加薄壁的浮冰在波浪中的运动来研究波-冰作用。通过数值模拟可以看出，短波长下波-冰作用更加剧烈。浮冰的存在一定程度上阻碍了波浪的运动，波浪场存在较明显的散射现象。

图 12(a)对应波长较短的工况，波浪场产生了明显的波浪散射现象，浮冰前侧和后侧的波浪液面形状由于散射的作用发生了明显的扰动现象。图 12(b)对应波长较长的工况，波浪散射现象较弱，浮冰周围波形没有明显变化。

如图 13所示，在波浪激励频率下，图 13(a)中波长越小，浮冰的垂荡运动幅值也越小。当波长较短时，浮冰大大阻碍了波浪的传递，同时也对波浪场造成了一定散射。由于波浪的散射消耗了一定的波能，使得浮冰的运动响应和受到的波浪力也随之减弱。随着波长增加，浮冰相对波长的尺度减小，进而浮冰对波浪的作用也减小，浮冰接收的波浪能量变大，因而波浪激励频率下的浮冰运动也随之剧烈。图 13(b)中，浮冰的纵摇运动随着波长的减小愈发剧烈。然而当波长变化至0.8 m时，波-冰特征长度比值为1.6，浮冰的纵摇运动开始减弱。这是因为波浪散射的程度随着波陡的增加愈发严重，波能耗散也愈发严重，从而使浮冰受到的作用力和运动响应减弱。

综上所述，随着波陡增大，波浪的散射现象越来越明显。波浪散射现象对浮冰垂荡具有一定程度的抑制作用。研究发现波浪散射现象越强，浮冰垂荡受到的抑制作用也越强。其中，对于浮冰的纵摇运动，随着波陡的增大，纵摇运动幅值增大，当波陡增大到一定值时，浮冰纵摇运动随着波陡的减小而减小。

4 结论

1) 通过对比有冲洗/无冲洗浮冰的运动发现，波浪冲洗现象对浮冰的运动具有明显的抑制作用。

2) 随着波陡变大，波浪冲洗现象和波浪散射现象增强，波浪冲洗现象对浮冰运动的抑制作用增强。同样地，波浪散射现象对浮冰垂荡运动具有抑制作用，波长越小，浮冰垂荡运动受到的抑制作用越强。

3) 波陡越大，浮冰的纵摇运动越大，当波陡增大到一定值时，由于波浪能量耗散严重，浮冰纵摇运动幅值随波陡增大呈现出下降趋势。

本文计算方法和分析结论对于进一步研究浮冰在波浪下的运动以及波浪中浮冰与结构物的相互作用具有参考价值。