气垫船破冰机理分析及临界航速估计

1 引 言

凌汛是黄河自然灾害的主要形式之一，黄河内蒙古段处于黄河流域的最北端，纬度高、气温低（零下30℃以下）、冰层厚，是凌汛灾害的多发地段。由于黄河内蒙古段水深较浅（1～3 m），传统的破冰船在该河段航行作业较为困难。目前黄河凌汛破冰减灾主要采用飞机投弹、炮击、人工爆破等炸冰方法，危险性大、破冰范围小、破冰效率低，且需要疏散当地群众，每年耗费大量的人力物力，存在一定弊端[1]。气垫船破冰是指采用大吨位气垫船以临界航速在冰面上航行，将冰层破裂成碎冰片流向下游，避免河道中冰堆和冰坝的形成，防患凌汛于未然，是一种较好的主动防御式破冰减灾方法。气垫船的破冰范围大、效率高、危险性小，目前在美国、加拿大和俄罗斯等国已有工程应用，取得了较好的破冰效果[2-3]。

在气垫船破冰机理研究方面，Mellor提出了气垫船破冰的两种不同模式:低速模式和高速模式。在低速模式下，气垫压力使得冰层下面的水面下陷，冰层由于悬臂梁效应而断裂；在高速模式下，气垫船兴起具有冰层弯曲断裂临界尺寸的波，航速范围为5～15 m/s。高速模式能使冰层破裂成较大的冰片，尺寸比气垫船的船长大，破冰效果比低速模式更好[4]。Buck也提出了类似的两种气垫船破冰模式:水面下陷模式和兴波模式。水面下陷模式主要适用于速度低、气垫压力高的非自推进式气垫船（气垫平台）；兴波模式则为速度高、气垫压力相对较低的自推进式气垫船的主要破冰方式[5-6]。

国内对浮冰层动力学问题已开展了不少研究，如冰载荷的计算方法和实验研究[7-9]，动载荷、脉冲载荷下浮冰层的动力分析等[10-12]，但针对气垫船航行条件下浮冰层动力学问题和气垫船破冰机理的研究则较为少见。气垫船破冰是一种较为经济实用的凌汛灾害预防方法，本文采用动力学分析软件LSDYNA对气垫船破冰问题进行数值模拟，根据计算结果探讨气垫船的破冰机理及其临界航速，对加快气垫船破冰方法在黄河凌汛灾害处理中的应用具有一定意义。兴波模式的破冰效率相对较高，本文主要对气垫船兴波模式的破冰过程进行数值模拟和分析。

2 原理和方法

根据黄河宁蒙段冰凌水文资料，该河段水深1～3 m，冰层厚通常0.3～0.7 m左右，最厚达1 m，封航期3～4个月。由此确定水平分层浅水条件下，气垫船破冰数值模拟的计算模型如图1所示，由下而上依次为水层、冰层、气垫船和空气层，为图示方便，空气层未画出。本文将气垫船载荷等效为在冰面上无摩擦滑行的柔性垫，通过给柔性垫施加重力和水平速度来近似模拟气垫船对冰层的载荷作用。图1中冰层视为线弹性固体层，采用拉格朗日算法，动力学基本控制方程组如下:

图1 气垫船破冰的计算模型
Fig.1 Sketch of calculating model for ice-breaking with air cushion vehicle(ACV)

式中:λ、μ为拉梅常量；u、ρ、f分别为冰层的位移、密度和体力。根据弹性动力学的变分原理，应用有限元方法对计算域进行空间离散，得到动力学控制方程的有限元计算方程如下:

式中:M、K和C分别为质量矩阵、刚度矩阵和阻尼矩阵；U和F分别为位移和载荷向量，表达式如下:

其中:N、B和E分别为单元形函数矩阵、几何矩阵和线弹性应力应变关系矩阵；f和T¯分别为单元的体力和应力边界条件；α和β为瑞雷阻尼系数。采用显式中心差分法求解上面的有限元计算方程，显式时间积分的最小时间步长由最小单元长度和波速决定。

水层和空气层为流体层，采用任意拉格朗日—欧拉（ALE）算法，质量、动量和能量守恒方程组成的控制方程如下:

式中:ρ为流体密度；v、w分别为物质速度和对流速度；σij为应力张量；b、E分别为体力和能量。

冰层和流体层之间的流固耦合采用欧拉-拉格朗日耦合算法。冰层破裂后，水会侵蚀进入冰层内部，因此采用罚函数耦合方式中的侵蚀算法约束类型，以便今后进一步实现冰层破裂以后的流固耦合计算。

柔性垫为线弹性固体，在柔性垫和冰层之间设置接触算法。对柔性垫施加重力后，通过接触算法传递均布面载荷到冰层表面，从而近似地模拟气垫船对冰层向下的压强面载荷。

3 计算参数

在图1的计算模型中，xz平面为对称面，水底满足不可穿透边界条件。计算域长300 m，宽64 m，水层、冰层和空气层的厚度分别为1.0 m、0.8 m和0.6 m。整个模型全部采用3D实体单元进行网格划分，节点数191932，单元数158464。用于破冰的气垫船一般吨位和尺寸较大，本文取长28 m、宽16 m的气垫船进行计算分析。模拟气垫船载荷的柔性垫为线弹性材料，厚度0.2 m，对冰层传递的压强值约为3000 Pa，等效于总吨位约135吨的气垫船载荷，气垫压力约为3000 Pa。

冰层的力学特性参数和冰的成分、形成过程等因素有关，不同性质的冰层，力学特性参数有所不同。参考文献[3]，冰层的密度、弹性模量和泊松比取值如表1所示。空气和水采用线性多项式状态方程描述，空气和水的密度分别为1.25 kg/m3、1000 kg/m3。

表1 计算参数
Tab.1 Parameters of calculation

4 数值模拟结果及破冰机理分析

根据Takizawa对浮冰层在移动载荷下的下沉位移测量结果，移动载荷速度对冰层下沉位移有明显影响。移动速度过小或过大时，下沉位移均较小；移动速度适中时冰层下沉位移最大，下沉范围最小，该速度称为临界速度[13]。本文主要对气垫船以较大航速（超临界航速）和临界航速在冰面上航行的两个算例进行数值模拟。

4.1 超临界航速

气垫船航速变化如图2所示，0～3 s内气垫船不移动，施加重力载荷，3～8 s内气垫船航速线性增加到10 m/s，然后保持不变到最后。

图2 气垫船航速变化曲线
Fig.2 The curve of ACV traveling speed

图3 气垫船兴波竖直位移云图
Fig.3 The contour map of vertical displacement of ice sheet caused by ACV

图3为超临界航速下气垫船兴波竖直位移云图，位移放大200倍显示，对应计算时刻分别为3 s、5 s、10 s、11 s、15 s和20 s。3 s时刻为气垫船启动前的瞬间，主要是重力加载引起的冰面下沉（-6.4 cm）；由5 s和10 s时刻的云图可见，随气垫船航速增大，在气垫船周围冰面上兴起船波，波峰在船首前方附近，波峰线呈弧形，波谷位于船尾后方；但随着气垫船继续向前航行，气垫船逐渐追赶并超越船首的波峰，对之前所兴起的船波起相反的压制作用，如11 s时刻云图所示；到15 s时刻，气垫船将船首波峰劈分成两半，使原来位于船身附近的船波向两侧分散，冰层变形范围随之扩大；20 s时刻，船首波基本不明显，原来所兴起的船波继续向两侧分散，气垫船不再兴起峰值较大的船波。因此，气垫船在超临界航速下航行时，对船首的兴波较为不利。

4.2 临界航速

气垫船航速变化如图4所示，0～3 s内气垫船不移动，施加重力载荷，3～5.8 s内气垫船航速线性增加到5.6 m/s，然后保持不变到最后。该临界航速的取值主要通过LS-DYNA软件多次试算，以气垫船刚好不超越船首波峰为依据确定。

图4 气垫船航速变化曲线
Fig.4 The curve of ACV traveling speed

图5 气垫船兴波竖直位移云图
Fig.5 The contour map of vertical displacement of ice sheet caused by ACV

图5为临界航速下气垫船兴波竖直位移云图，可见在气垫船启动、加速和巡航的整个过程中，气垫船始终位于船首波峰后方，对船首波起持续的推波作用。船波一直分布于船身附近，冰层变形范围集中，有利于气垫船利用兴波能量进行破冰。

图6 冰面竖直位移包络曲线
Fig.6 The envelope curve of vertical displacement at the surface of ice sheet

图7 冰层上表面最大主应力包络曲线
Fig.7 The envelope curve of maximum principle stress at the surface of ice sheet

将气垫船中轴线上所有冰面节点的竖直位移随时间变化曲线在图6中列出，全部曲线的外部包络线反映了气垫船启动、加速和巡航过程中，中轴线上冰面竖直位移最大、最小值随时间的变化规律。类似地，中轴线上冰层上表面和冰层底面最大主应力值随时间的变化规律如图7和图8中包络线所示。图6中，上部包络线对应于船首波的峰值位移，下部包络线对应船尾波的谷值位移。可见，船首波的峰值位移持续增大，最大值约7 cm；船尾波的谷值位移增大后基本不变，最大值约为-10 cm。类似地，图7的包络线对应于船首波峰处应力，图8的包络线对应于船尾波谷处应力。可见，冰层上表面波峰处的最大主应力值一直增加，最大值约0.65 MPa；波谷处最大主应力则增加后略有减小，最大值0.75 MPa。冰层上表面和冰层底面最大主应力空间分布如图9～10所示，可见船尾波谷处的最大主应力值相对大一些，变形范围相对波峰处则更为集中。

图8 冰层底面最大主应力包络曲线
Fig.8 The envelope curve of maximum principle stress at the bottom of ice sheet

图9 冰层上表面最大主应力云图
Fig.9 The contour map of maximum principle stress at the surface of ice sheet

图10 冰层底面最大主应力云图
Fig.10 The contour map of maximum principle stress at the bottom of ice sheet

4.3 破冰机理分析

根据上述气垫船在超临界航速和临界航速下的数值模拟结果，初步分析气垫船破冰机理如下:

在超临界航速下，气垫船逐渐追赶并超越船首的波峰，对之前所兴起的船波起相反的压制作用。气垫船将船首波峰劈分成两半，使原来位于船身附近的船波向两侧分散，气垫船不能持续对所兴起的船波补充能量，难以兴起峰值较大的船波，对气垫船破冰不利。

当气垫船以临界航速航行时，气垫船始终位于船首波峰后方，对原来所兴起的船波起持续的推波作用，不断补充波动的能量。相应地，船首波的峰值位移持续增大；船尾波的谷值位移增大后基本维持不变。冰层上表面波峰处的最大主应力持续增加；冰层底面波谷处最大主应力增加后略有减小，但总体上变化不大。因此，当气垫船以临界航速持续推波达到一定程度时，冰面上兴起峰值较大的船波，当冰层最大主应力达到冰层抗拉强度极限时，冰层出现裂缝并扩展，从而实现气垫船破冰。

5 临界航速估计

气垫船以临界速度航行时，冰面波形曲线的数值模拟结果如图11所示。波形曲线包含不同波长的波，其中对破冰起作用的波位于气垫船船身附近，峰值位移和谷值位移均较大。该波的波长与气垫船船长、船型、冰层厚度等诸多因素有关。由图11可见，在本文的计算条件下，波峰基本位于船首前方约半个船长处，波谷基本位于船尾，波长约3倍船长左右。当气垫船的航速和该波的行进速度相同时，气垫船才能始终位于波峰后方，发挥持续的推波作用。气垫船航速大于或小于该波的行进速度，均不能起到持续推波及补充能量的作用。因此，估计气垫船破冰的临界航速，实际上就是确定船身附近表面波的行进速度。

由图1可以看出，气垫船在冰层表面兴起的波动实际上是液固分层介质中的表面波。表面波一般存在频散现象，即波长不同波速也不同。Milinazzo给出了液固分层介质中表面波相速度的表达式如下[14]:

图11 气垫船临界速度航行时的波形曲线
Fig.11 The waveform of ice-sheet when ACV traveling at critical speed

式中:g是重力加速度；H 为水深；k=κH，κ 为波数；α=D/（ ρw gH4），D=Eh3/[12 （1- ν2）]， 其中 h、E 和 ν分别为冰层的厚度、弹性模量和泊松比；γ=ηh/H，η=ρ/ρw，其中 ρ、 ρw分别是冰层和水的密度。

由κ=2π/L，L为波长，按前面估计的3倍船长取值，再将本文计算模型的参数代入（10）式，即得到船身附近表面波的波速为5.56 m/s，与数值模拟多次试算的气垫船临界航速（5.6 m/s）基本相同。另外，本文对水深为3.0 m，冰层厚度为0.9 m的破冰条件也进行了计算。按上述估计方法得出的船身附近表面波的波速为7.75 m/s，与数值模拟多次试算的气垫船临界航速（7.8 m/s）也基本相同。因此，在本文浅水厚冰层环境下，气垫船的临界航速可以按波长为3倍船长的表面波的行进速度进行估计。

6 结 论

本文采用动力学分析软件LS-DYNA对简单水平分层浅水厚冰层环境下气垫船破冰问题进行了数值模拟，数值计算结果表明:

（1）超临界航速下，气垫船逐渐追赶并超越船首的波峰，对所兴起的船波起相反的压制作用，难于在冰面上兴起峰值较大的船波，对气垫船破冰不利。

（2）当气垫船以临界航速航行时，气垫船始终位于船首波峰后方，对所兴起的船波起持续的推波作用，不断补充波动的能量。

（3）当气垫船以临界航速持续推波到一定程度时，冰面上兴起峰值较大的船波，冰层最大主应力达到冰层抗拉强度极限，冰层将产生裂缝并扩展，从而实现气垫船破冰。

（4）气垫船破冰的临界航速实际上就是船身附近表面波的行进速度，在本文浅水厚冰层环境下，气垫船的临界航速可以按波长为3倍船长的表面波的行进速度进行估计。

本文对气垫船载荷作用下浮冰层的动力学研究主要集中在冰层未发生破裂之前，浮冰层的破裂过程和破裂后的动力特性对进一步明确气垫船的破冰机理和制定气垫船破冰的实际运行模式具有一定参考价值，需要施加冰层的抗拉强度破坏准则后继续进行研究。