船舶艏部外飘砰击压力的仿真研究

赵留平

(海军驻武汉701 所军事代表室，武汉 430033)

船舶在海浪中航行，其运动过程底部不一定出水，当波浪冲击艏部外飘区域时，往往会产生较大的脉冲水动力，称之为外飘砰击［1］。外飘砰击具有砰击压力小、作用表面积大、砰击过程持续时间长等特点。确定船舶外飘压力大小，对船体结构设计至关重要，因此，对于如集装箱船等具有外飘形式的船舶来说，在设计初期，必须考虑这种砰击的影响。

对于结构物入水砰击过程，国外Von Karman［2］、Wagner［3］等诸多学者在理论上进行了深入研究，国内陈震［4］等人利用三维船体表面冲击压力计算方法与船舶耐波性预报相结合的理论手段，预报了大型集装箱船舷侧外飘砰击压力。

但由于结构体尤其是船壳等复杂结构体入水过程是一个非线性、非定常问题，物理现象十分复杂，许多砰击现象的本质特性和物理机理难以用理论来诠释，理论分析受到各种因素的限制，如物体几何外形、入水速度、非线性自由边界条件等。近年来，计算机的高速发展使得学者可以通过数值模拟来研究结构的砰击问题。Verhagen［5］、金伏生［6］等人采用有限差分法，Zhao＆Faltisen［7］等则采用边界元法，但以上大多数学者的工作主要还是以二维仿真为主。

本文应用MSC.Dytran 瞬态动力学软件，对三维船舶艏部模型进行仿真计算，考虑了大气压和重力等因素的影响，计算艏部在不同航速下的外飘砰击压力，并分析得出砰击压力的分布情况，为后续的船舶艏部线型设计奠定基础。

1 DYTRAN 基本分析理论

本文采用一般流固耦合算法，即采用拉格朗日有限元计算法求解固体模型;欧拉有限体积法求解流体模型。

1.1 拉格朗日求解法

MSC.Dytran 是瞬态动力学流—固耦合领域的高端软件，它的求解方法在时间域上均采用显式时间积分法。当前时间步是步n，显式积分方法是将运动微分方程

改写为

如果将单元质量分布在节点上，则M 为一对角阵，称为集中质量矩阵，线性方程组将成为一系列关于各个自由度的独立的一元一次方程，从而可求出节点加速度

如果假设加速度在一个时间步长内是恒定的，在时间推进上采用中心差分法

显式积分法不需要做矩阵分解，因此具有很高的计算效率。

1.2 欧拉有限体积法

初始条件已知情况下，利用控制方程［8］求解每一欧拉单元在一个时间步的密度、速度、比内能和压力。假设tn时刻各物理参数已知，对相邻单元形心处流速进行线性插值求出单元边界处流速

求出穿越单元表面的质量、动量及能量的流量

式(8 ～10)中: ΔV 表示从时刻tn～tn+1的一个时间步长内流过该单元表面的体积流量;ρ2为相邻单元密度。

1.3 流—固耦合

通过在拉格朗日模型上定义一层耦合面，建立耦合关系，作为拉格朗日部分和欧拉部分之间的传递者。拉格朗日部分直接从耦合面上接受载荷。欧拉部分将该面充当流场边界，并将欧拉单元的应力作用到耦合面上，引起拉格朗日单元的变形。主要通过3 步进行耦合计算［9］: 耦合计算步;运输步;冲量步。其中耦合步计算了耦合面与各单元相交的情况，运输步负责传递单元的质量及与质量相交的变量，冲量步加入了压力波在网格传播中的贡献。

2 计算模型

以某船艏部为例，计算船舶艏部舷侧区域的外飘砰击压力。该船的主尺度数据见表1，船艏横剖面型线见图1。

表1 主尺度数据

图1 横剖面型线

在建立三维艏部模型基础上，利用MSC. Dytran 软件进行艏部砰击问题的数值计算。在仿真计算过程中，如图2 所示，计算域尺寸的长度为艏部长度的4 倍，宽度为艏部宽度的5 倍，高度为艏部高度的7 倍，并施加无反射边界条件以保证压力波在流场中传播不至于反射回来影响砰击区域。

图2 计算域

欧拉区域采用不等密度网格划分，在离艏部模型结构较近的位置欧拉网格划分较密，远离艏部结构的区域较粗，合理布置流场区域的欧拉网格，能有效提高计算效率。

艏部模型结构用拉格朗日单元网格划分，采用刚性体材料。空气和水域场均采用欧拉单元网格划分，空气单元采用可压缩理想气体本构关系的材料填充，水单元采用无黏性、可压缩线性流体本构关系材料填充。图3 为计算模型的局部网格划分情况。

整个艏部模型封闭的外表面被定义为流固耦合面，采用一般耦合算法(General coupling)计算流体与艏部之间的相互作用。计算模型的各参数设置详见参考文献［10］。

图3 局部网格划分

3 结果与分析

船舶艏部模型在5 m 浪高海况，10.0 kn、15.0 kn 以及20.0 kn 不同航速下，根据船舶艏部线型特点，在0 ～3 站和9 ～12 m 水线船舷区域内的外飘砰击压力极值计算结果如表2 所示。

表2 外飘砰击压力极值(×105 Pa)

3.1 外飘砰击压力沿船长方向分布

图4 和图5 为5 m 浪高条件下，不同航速时船艏9 m 水线和12 m 水线处砰击压力极值沿船长方向的变化情况。

图4 9 m 水线砰击压力

图5 12 m 水线砰击压力

由图4、图5 可以看出，在5 m 浪高下，外飘砰击压力极值在同一水线高度处，沿船艏到船尾方向，砰击压力极值逐渐减小。

以上砰击压力变化趋势主要由船体表面线型特征决定，即从0 至3 号各站位处外飘角度逐渐减小，进而根据Wagner经典砰击理论可知，外飘砰击压力逐渐减小。

3.2 外飘砰击压力沿舷侧高度方向分布

图6 和图7 分别为5 m 浪高条件下，不同航速时船艏0站和3 站处砰击压力极值沿高度方向的变化情况。

图6 0 站处砰击压力

图7 3 站处砰击压力

从图6 和图7 可知，在相同浪高条件下，外飘砰击压力极值沿高度方向的变化规律一致:外飘砰击压力随着距水线高度的增加而增大。

这是由于在船艏区域，随着距水线高度的增加，舷侧线型外飘角度逐渐增大，进而由Wagner 砰击理论可知，随着高度的增加，外飘砰击压力增大。

4 结论

对于三维结构体的外飘砰击问题，理论研究有较多局限。本文通过对三维船舶艏部模型的外飘砰击压力仿真计算，得到以下结论:

1)通过仿真计算，此次分析得出外飘砰击压力分布情况满足经典Wagner 砰击理论，本文的数值仿真方法在研究三维砰击问题方面是可行的。

2)船舶艏部外飘砰击压力的大小与船体表面线型有较大关系，外飘砰击压力最大值一般出现在船艏后0 ～1 站外飘线型明显的区域。

3)对于船舶艏部线型类似集装箱船的船舶来说，沿船长方向，外飘砰击压力极值随着距船艏距离增大而降低; 在高度方向，随着距水线高度的增加而增大。

［1］戴仰山，沈进威，宋竞正.船舶波浪载荷［M］.北京:国防工业出版社，2007.1.

［2］Von Karman.The impact of seaplane floats during landing［Z］.NACA TN321，1929.

［3］Wagner. Trans. phenomena associated with impacts and sliding on liquid surfaces［Z］. Math Mech，1932，14:153 68.

［4］陈震，冯永军，肖熙.大型集装箱船舷侧外飘砰击特性研究［J］.船海工程，2011.

［5］VERHAGEN J H G.The impact of a flat plate on a water surface［J］. Journal of Ship Research，1967，11(4): 211-223.

［6］金伏生.入水冲击问题变分原理及其它［J］.应用数学和力学，1992，13(6):543-552.

［7］Zhao R，Faltinsen O. Water entry of two - dimensional bodies［J］. Journal of Fluid Mechanics，1993，246: 593-612.

［8］卞文杰，万历，吴莘馨.瞬态动力学CAE 解决方案MSC.Dytran 基础教程［M］. 北京: 北京大学出版社，2004:4-7.

［9］丁沛然，钱纯.非线性瞬态动力学分析-MSC.Dytran 理论及应用［M］.北京:科学出版社，2006:1-2.

［10］陈震，肖熙.空气垫在平底结构入水砰击中作用的仿真分析［J］.上海交通大学学报，2005，39(5):670-674.