弹性结构入水砰击载荷特性三维数值模拟研究

杨 衡，孙龙泉，龚小超，姚熊亮

（哈尔滨工程大学 船舶工程学院，哈尔滨 150001）

结构入水砰击是典型的结构与流体瞬态流固耦合作用过程，在船舶工程、航空领域等中十分常见，如船艏砰击、海上救生艇抛落、空投鱼雷入水、水上飞机降落着水问题等。入水初期砰击载荷以“瞬态”、“大幅”为特点，其砰击压力的脉宽往往在毫秒量级，瞬态冲击载荷作用下往往使结构产生高频的弹性振动，甚至对局部结构产生塑性破坏，因此，入水砰击问题越来越引起人们的重视，同时其在军事和民用上有着很现实的工程应用背景。结构自空中进入水面的过程中所承受的外部载荷大小外部载荷主要来自流体动压力即流体冲击压力。当冲击发生时，流体动压力沿物体湿表面的分布特征及其大小取决于入水速度、入水角、结构形式、流体粘性、可压缩性、结构弹性以及其它一些因素［1］，张阿漫等［2］通过实验方法得到不同环境下条件气泡脉动特性，得到刚性平面、45°斜面及舷侧曲面条件下气泡射流砰击规律，对不同入水角和不同底部形式入水砰击问题提供研究思路。由于结构物入水冲击过程涉及结构运动非线性、自由表面运动非线性、气穴效应、流固耦合效应、重力效应等，特别是伴随自由面的抬升、飞溅，入水空泡现象等，在数学处理上有很大困难，张阿漫等［3］基于边界积分开展近自由面水下爆炸气泡作用下的运动规律进行研究，对自由面运动进行较好模拟。自Von Karman［4］提出用附加质量法计算物体入水砰击载荷以来，结构的入水砰击问题得到了广泛而深入的研究，研究对象为简单二维楔形刚性体，基于势流理论进行求解。在弹性体入水砰击研究中，流体与结构的耦合效应成为研究重点，卢炽华等［1］通过非线性伯努利方程，将流体域边界元方程与结构有限元方程耦合求解，并考虑了自由表面的非线性边界条件。陈铁云等［5］采用混合欧拉－拉格朗日方法对高速船舶砰击现象的冲击载荷进行数值计算。顾慈祥等［6］计算了平头旋转壳体垂直入水时的水弹性问题。在国外对水弹性在入水砰击问题上进行分析计算的有Sharov等［7－9］。在入水问题的非线性方面，Vorus［10］提出了一种改进的wagner模型，在其计算中，考虑了自由液面条件和伯努利条件的非线性，而物面条件是利用相当平板来模拟的。Korobkin［11］发展了Logvinovich提出的模型，原始Logvinovich方法（OLM）利用相当平板理论模拟物面条件，对自由面条件进行了线性化处理，而计算压力的伯努利方程则是非线性的。Wu等［12］采用边界元法讨论了全非线性自由面条件下楔形体自由下落的模型。此外，SPH法，CIP法，VOF方法等数值方法都已经应用在入水冲击问题的求解上。本文基于边界元方法，结构边界为流固耦合交界面，计算效率较SPH、VOF方法明显提高，计算中考虑了结构运动非线性、流固耦合效应等，对结构入水问题进行了数值模拟；但基于边界元方法计算结构入水砰击问题时，无法精确模拟流体飞溅、自由液面大变形等问题。

本文的基于双渐进法［13－14］（DAA），推导了考虑流体可压缩性［15］的结构入水冲击过程非线性双渐进法（NDAA）流固耦合运动方程，研究分析弹性体结构入水冲击过程的流固耦合效应。流体域速度势采用边界元方法计算，自由面速度势满足自由表面的非线性条件。结构为弹性体，运用有限元方法获得结构的弹性响应，通过非线性伯努利方程将流体域速度势和结构弹性响应耦合求解，通过与模型试验值进行对比，验证本文方法可靠性基础上，讨论弹性体与刚性体入水砰击过程载荷特性差异，着重探讨了入水速度对弹性结构入水砰击载荷与弹性动力响应特性的影响。

1 结构入水NDAA数值计算方法

结构入水冲击过程是一个流体和结构瞬态、强非线性的相互作用过程，为了将问题简化，假设流场为各向同性、无粘、无旋但可压缩的理想流体，计算中考虑结构的重力效应，计算模型如图1所示。本文将结构视为弹塑性结构，在重力作用下以一定初速度入水，结构在外界水动力的作用下做流固耦合运动，同时结构的动力响应又反作用于流场，实现结构入水过程中流固耦合运动的数值模拟。

考虑可压缩性的流体速度势满足如下控制方程：

假设可压缩性对自由液面速度势影响可忽略，计算模型需要满足如下边界条件［16］：

自由面上SF运动学边界条件：

图1 数值模型Fig．1 Numerical model

由面上SF动力学边界条件：

在物面SB上：

在无穷远边界S∞边界条件：

初始边界条件：

式中，g为重力加速度；n是法向矢量坐标；V为结构运动速度，ω为弹性结构的振动速度。

DAA法就是从控制方程（1）的基础上经推导得到。二阶DAA基本方程［17］如下：

结构动力方程

式中：ps为流体中的散射压力，Mf为流体质量矩阵，Ωf为流体频率矩阵，Af为流体单元的面积矩阵，G为坐标转换矩阵。

DAA方程在其理论推导的过程中，均引入了一个线性假设，即

二阶DAA法是从延迟势法基础上得到，延迟势的速度势方程如下：

延迟势的压力方程

方程式（11）的求解可以得到流场散射压力ps，流场总速度势即可由得解。

对于结构入水砰击问题，以“瞬态”，“大幅”为特点的非线性运动，其非线性效应对水动力值影响比较显著，速度势的空间导数就不能忽略了。本文在二阶DAA法基础上，将速度势的空间导数引入到二阶DAA法的数值计算，从而将非线性的影响考虑到结构入水砰击的瞬态流固耦合计算当中。

在得到流场速度势的基础上，可以通过非线性伯努利方程［16］得到流场动压力Pd。

式中：V为结构运动速度。

式（12）得到的流场动压力值与二阶DAA方程得到的流场动压力值的本质区别在于式（12）中考虑了结构的速度效应项及结构的运动非线性项

求解式（12）得到考虑非线性效应的流场动压力，将该考虑非线性效应的流场动压力Pd代替结构运动方程中ps，结构运动方程变为：

联合求解式二阶DAA方程、式（12）、（13）即为非线性双渐进法（NDAA法）。

2 圆柱壳入水砰击载荷特性分析

2.1 数值验证

应用本文方法计算编制计算程序，计算结果与文献［16］某缩比模型入水过程试验数据进行对比分析。计算模型与试验模型具有相同质量分布、惯性矩、材料特性，计算初始参数如下：入水角78°，入水速度Vx＝1．341 3 m／s、Vy＝0．835 8 m／s、Vz＝－10．068 9 m／s。下图给出结构入水过程中三方向速度及砰击压力计算值与试验值对比结果。

图2 速度曲线对比Fig．2 The compare of the velocities between the numerical result and the experiment

图3 砰击载荷对比Fig．3 The compare of the slamming loads

图2是缩比模型入水过程中三个方向的速度变化曲线，V1是横向速度、V2是法向速度、V3是垂向速度，从图中可以看出，圆柱入水过程中受到的垂向载荷变化剧烈，圆柱的垂向速度变化剧烈，相比于垂向速度和法向速度，横向速度的变化曲线最为平缓。

图3给出距离模型尾部4．6 m处测点入水砰击载荷对比，从图中可以看出，应用本文计算方法得到的砰击载荷的计算值与试验值峰值、脉宽相近，说明本文方法的可靠性。5．但本文计算模型由于未考虑结构自身阻尼特性，流场的飞溅、自由液面变形引起的水动力变化，仅考虑结构自身在水动力作用下的弹性振动，计算得到的压力曲线与试验压力曲线在0．5 MPa、0．25 MPa附近的震荡形式存在一定差别。

2.2 弹性圆柱体倾斜入水模拟

入水角度为30°，垂向速度为30 m／s，水平速度为0 m／s。取圆柱底部与侧部迎水面特征节点进行分析，特征节点分布如图4所示，其中节点70为入水点。

图4 节点分布Fig．4 The distribution of notes

图5 底部节点砰击压力对比Fig．5 The compare of the slamming loads on button

采用本文计算方法计算圆柱体倾斜入水过程各特征节点砰击压力分布及总体受力。图5、6给出各节点砰击压力对比。

由图可知，各节点砰击压力呈现相同的规律，即特征点接触到流体时，会产生较大的砰击峰值，砰击压力迅速下降，形成二次砰击峰值，二次砰击峰值约为第一次砰击峰值的1／3。由于结构的弹性效应，入水初期，圆柱特征节点砰击压力呈现波动，随入水深度的增加，结构表面的压力主要为静水压力，此时，节点的压力波动情况减弱，随侵水深度的增加，压力略有增加。

比较七个特征节点的压力变化时历曲线，由图5可以看出，底部三个节点中，节点70压力峰值最大，越靠近底部中心压力峰值越小；而比较图6中各节点压力曲线，远离初始入水点的特征节点的压力峰值较大。

图6 侧部节点砰击压力对比Fig．6 The compare of the slamming loads on side

图7 水平方向合力变化曲线Fig．7 The history curve of the horizontal force

图8 垂直方向合力变化曲线Fig．8 The history curve of the vertical force

图9 俯仰力矩变化曲线Fig．9 The history curve of the moment

图10 俯仰角变化曲线Fig．10 The history curveof the pitch angle

图11 俯仰角速度变化曲线Fig．11 The history curve of the pitch angle velocity

图7、8给出圆柱入水过程水平方向和垂直方向合力时历曲线。如图6所示，水平及垂向合力存在双峰值现象，分别对应圆柱体尾部入水和头部入水两个时间点，尾入水时刻水平合力较头入水时刻略小，垂向合力较头部入水时刻减小一半。

结构入水过程中，由于流体载荷的作用不均匀，会产生绕质心转动的俯仰力矩，俯仰力矩的变化是结构入水姿态变化的主要因素。下图给出俯仰力矩变化曲线。

由图9可以看出，在圆柱头尾入水时刻，俯仰力矩也会出现两个峰值，与圆柱体所受到的合外力的变化相关。下图给出在俯仰力矩作用下，圆柱入水过程俯仰角变化及圆柱质心位置的运动轨迹。

由俯仰角曲线可以看出，在整个入水过程中俯仰角逐渐减小，在圆柱质心入水前，俯仰角速度先反向增大，质心入水后，俯仰角速度减小，俯仰角速度的变化趋势与俯仰力矩的变化曲线相同。

图12 圆柱质心运动轨迹Fig．12 The motion track of the center of mass

图12给出圆柱体入水过程质心位置变化，由图可以看出，在入水初期，圆柱质心先向水平方向正向移动，随入水深度的增加，质心向水平方向负向移动，质心运动曲线变化与圆柱水平方向合力时历曲线变化密切相关，质心运动变化率随水平方向合力变化而变化。

3 舰船入水砰击载荷特性分析

在船舶与海洋工程领域，砰击现象广泛存在，如救生艇抛落、船舶从船台下水、高速滑行艇艇底砰击、甲板上浪等，砰击的存在使船体局部产生应力集中区，船体产生局部变形甚至破坏。船舶在海上航行时，砰击可能发生在船艏、船底、船艉等位置，尤其以底部砰击最为重要，船底入水过程中受到巨大的砰击载荷，会改变船体的垂向加速度，船底局部产生高频振动，第一次世界大战期间，发生了多起因船底砰击而损坏的事故，船底砰击越来越受到人们重视。

本节在前文基础上，将双渐近法应用于大型水面舰船的砰击载荷计算中，相对于前文的回转体，在计算大型水面舰船的入水砰击流固耦合特性时，计算域变大，结构外形更加复杂，在砰击载荷作用下弹性变形更显著。同时，根据图13船体入水砰击试验，船体入水过程中无空泡产生，且本文在计算中忽略了自由液面飞溅及变形的影响，但考虑了自由液面效应、流场可压缩性、船体运动的非线性及伯努利方程的非线性。

图13 船体入水砰击Fig．13 Slamming Character on Water Entry of Hull

本节采用双渐近法对某型舰船入水砰击流固耦合运动进行数值模拟，研究船体垂直与倾斜入水砰击过程中载特性及船体变形。某型舰船计算模型如图13所示，船长为38．6 m，型深为 3．78 m。

图14 船体型线Fig．14 The shiplines

3.1 舰船垂直入水砰击

本节讨论该水面舰船以初速度（30 m／s）入水，除受到因船底砰击而产生的水动力作用外，不受其他任何外力的作用，研究该舰船的入水砰击过程中受力状态（流场动压力）、Mises应力变化及船底变形。

在该工况计算中，静止的流场在船底入水瞬间受到剧烈扰动，产生较大的砰击载荷作用于船体表面，由于船底型线不均匀性，首先接触水面的局部产生较大的砰击载荷，随入水深度的增加，应力逐步向船尾船底。下图给出船体垂直入水过程中侧部和底部Mises应力变化云图。

图15船体入水初期Mises应力变化云图，距船艏部0．177 m处首先触水，Mises应力逐渐增大，在船底，Mises应力向船艉逐渐传递，在船体侧部，Mises应力向船甲板逐渐传递。由于船体对称性，Mises应力分布呈对称分布。在整个入水过程中，船体表面产生微幅振动。

取船底三个节点进行受力分析，如图16所示，节点26为首先触水点。下图给出三个节点入水砰击过程中垂向加速度和砰击压力对比。

图15 船体Mises应力变化云图Fig．15 The nephogram of Mises stress changes

图16 特征节点分布Fig．16 The distribution of notes

由图17、18可以看出，在船底砰击过程中，越靠近船艉处，节点的垂向加速度及砰击载荷越大。当船体开始入水后，由于砰击载荷的作用范围较小，船体整体仍处于加速运动状态，砰击载荷增加。

图19，20给出船底三个节点弹性振动位移曲线及船体变形图。

由图19、20可以看出，船底有向内凹的趋势。由于船底左右底板相交的区域为尖状设计，如节点26处，并未产生明显的变形，而船底型线之间的部分及船侧部，如节点15，变性较大，局部振动明显，而尾部节点6，由于船艉板的影响，其弹性变形要小于船中部。因此，在船底局部设计过程中应对应力集中区及船底中部进行局部加强设计。

图17 垂向加速度对比Fig．17 The compare of the vertical accalarations

图18 砰击压力对比Fig．18 The compare of the slamming loads

图19 弹性振动位移对比Fig．19 The compare of the elastic vibration displacement

图20 船底变形前后对比（1∶2）Fig．20 The compare of the hull deformation（1∶2）

3.2 舰船倾斜入水砰击

本节讨论该水面舰船横向倾斜10°入水砰击，以初速度（垂向速度20 m／s、法向速度5 m／s）入水，下图给出船体倾斜入水过程中侧部和底部Mises应力变化云图。

由图21船体入水初期Mises应力变化云图可以看出，船体尾部首先触水，Mises应力逐渐增大，在船底，Mises应力呈扇形逐渐传递，在船体侧部，Mises应力向船甲板逐渐传递。

图22 特征节点布置Fig．22 The distribution of notes

取船底三个节点进行受力分析，如图22所示。图23，24给出该船倾斜入水过程中三个节点入水砰击过程中垂向加速度和砰击压力对比。

图23 垂向加速度对比Fig．23 The compare of the vertical accalarations

图24 砰击压力对比Fig．24 The compare of the slamming loads

图25 船体变形前后对比（1∶2）Fig．25 The compare of the hull deformation（1∶2）

由图23，24可以看出，在船底砰击过程中，先触水的节点的垂向加速度及砰击载荷较后触水的小。

从船底变形图中可以看出，船底倾斜入水时，船底及侧部变形较垂直入水砰击时产生较大不同，船体先触水一侧的船底与侧部产生较大变形。

4 结 论

本文基于波动理论，采用双渐近方法，建立了结构入水砰击流固耦合数值模型，计算中考虑了自由液面非线性条件、运动非线性及重力效应，对弹性圆柱体以及某型舰船入水砰击过程（以某一初速度运动、瞬态砰击载荷作用下的受迫振动以及后期自由振动）进行流固耦合数值模拟，实现对本文所开发程序的工程化应用，结果表明：

（1）结构入水过程中受到明显的砰击载荷，对于弹性体入水，砰击载荷作用后，由于结构自身的弹性，结构底部产生明显的弹性振动，砰击压力逐渐衰减。

（2）圆柱体倾斜入水时，刚弹耦合效应明显，结构入水过程中水平、垂向合力在头部入水和尾部入水时刻出现两个峰值。

（3）船体垂直入水过程中，由于船底结构形式的影响，Misses应力分布沿船底中线对称分布，船底在砰击载荷的作用下出现较大变形。船底倾斜入水时，舷侧也出现较大变形，可见，控制入水砰击的角度，对船舶入水砰击安全性至关重要。

［1］卢炽华，何友声．二维弹性结构入水冲击过程中的流固耦合效应［J］．力学学报，2010（3）：129－140．LU Chi-hua，HE You-sheng．Coupled analysis of nonlinear interaction between fluid and structure during impact［J］．Acta Mechanica Sinica，2010（3）：129－140．

［2］张阿漫，王诗平，白兆宏，等．不同环境下气泡脉动特性实验研究［J］．力学学报，2011，43（1）：71－83．ZHANG A-man，WANG Shi-ping，BAI Zhao-hong，et al．Experimental study on bubble pulse features under different circumstances［J］．Chinese Journal of Theoretical and Applied Mechanics，2011，43（1）：71－83．

［3］张阿漫，姚熊亮．近自由面水下爆炸气泡的运动规律研究［J］．物理学报，2008，57（1）：339－353．ZHANG A-man，YAO Xiong-liang．The law of the underwater explosion bubble motion near free surface［J］．Acta Physica Sinica，2008，57（1）：339－353．

［4］Von Karman T．The impact of seaplane floats during landing［R］．Washingt on DC，USA：National Advisory Committee for Aeronautics，NACA Technical Notes，1929，321．

［5］陈铁云，王刚．高速船舶结构设计中流体冲击载荷的数值计算［J］．上海交通大学学报，1998，32（11）：30－33．CHEN Tie-yun，WANG Gang． Numerical calculation of hydrodynamic impulsive loading in structural design of high speed ship［J］．Journal of Shanghai Jiaotong University，1998，32（11）：30－33．

［6］顾懋祥，程贯一，张效慈．平头旋转壳撞水水弹性效应的研究［J］．水动力学研究与进展，1991，6（1）：41－51．GU Mao-xiang， CHENG Guan-yi， ZHANG Xiao-ci． The study of hydroelastic effect on the water-entry of a revolutionary shell with a flat nose． ［J］． Journal of Hydrodynamics，1991，6（1）：41－51．

［7］Sharov V F．Ship bottom impact upon wave ［J］．Sudostroyeniye，1958，4：5－9．

［8］Meyerhoff E K．Die berechnung hzdroelastis cherstobe［J］．Schiffstechnik，1965，12：19－30．

［9］Zhao R，Faltinsen O．Water entry of two-dimensional bodies［J］．Fluid Mech，1993，（246）：593－612．

［10］Vorus W．A flat cylinder theory for vessel impact and steady planin resistance［J］．Journal of Ship Research，1996，40（2）：89－106．

［11］Korobkin A A，Scolan Y M．Three-dimensional theory of water impact：Part2．Linearized Wagner problem ［J］．Journal of Fluid Mechanics，2006，549：343－373．

［12］Wu G X，Sun H，He Y S．Numerical simulation and experimental study of water entry of a wedge in free fall motion［J］．Journal of Fluids and Strucures，2004，19：277－289．

［13］Geers T L，Residual potential and approximation methods for three dimensional fluid-structure interaction problems［J］．Acoust．Soc．Amer，1971，49：1505－1510．

［14］ Geers T L，Doubly asympotic approximation for transient motions of submerged structures［J］．Acoust．Soc．Ame，1978，64：1500－1508．

［15］王诗平，孙士丽，张阿漫，等．可压缩流场中气泡脉动数值模拟［J］．力学学报，2012，44（3）：513－519．WANG Shi-ping，SUN Shi-li， ZHANG A-man，et al．Numerical simulation of bubble dynamics in compressible fluid［J］．Acta Mechanica Sinica，2012，44（3）：513－519．

［16］孙士丽．具有自由液面效应的高频流固耦合问题［D］．哈尔滨：哈尔滨工程大学，2010．