入水砰击问题的三维数值模拟

入水砰击问题的三维数值模拟

0船舶理解和试验研究结构物进入水中的问题存在于许多情况下，例如船体沉降、救援、飞机（水上飞机、宇宙返回等）掉落、武器（深水爆炸、鱼雷等）爆炸等。当砰击发生时,结构物在非常短的时间内受到了较大的冲击力,有时这种冲击力甚至会破坏物体结构和内部仪器,因此,确定结构物入水时受到的冲击力大小,对结构设计至关重要。世界上几个主要船级社,比如挪威船级社、英国劳氏船级社、中国船级社等均对船舶砰击进行了研究并提出了计算方法。VonKarman、Wagner等诸多学者在理论上相继进行了深入研究,但由于结构物入水过程是一个涉及固-气-液三相耦合的非线性、非定常问题,物理现象十分复杂,许多砰击现象的本质特性和物理机理难以用理论分析来诠释。同时,数值模拟的结果要依靠试验来检验,因而试验研究就显得尤为重要。Chuang等人率先完成了平板、楔形体、圆锥体、船模的入水试验,分析了砰击压力与相关因素的关系,提出了砰击压力的回归经验公式以及其他重要结论;Moghisi等人完成了球体入水试验。但理论与试验研究均受到各种因素限制,如物体几何外形、入水速度、非线性自由表面边界条件、试验成本、尺缩效应、测量仪器准确度等。近年来,计算机的高速发展使得人们可以通过数值模拟来研究结构物的入水问题。采用有限差分法的主要有Verhagen、Koehler、金伏生、顾懋祥等人;采用边界元法的主要有Geers、Zhao和Faltinsen以及叶取源;采用有限元法的主要有Marcal、陈震等人。以上大多学者的工作主要还是以二维仿真为主。本文基于商用软件Fluent对三维平板与圆锥体的砰击压力进行了数值计算,并与黄震球、ChuangSL的试验结果作了对比。1平板的黄震球响应特性本文对平板的两种入水方式进行了计算与比较:(1)自由落体,通过6DOF选项实现;通过编写用户自定义函数(UDF)赋予平板与黄震球试验模型一致的质量940kg;为保证平板基本只有垂向运动,3个惯性矩被赋予了一个大数。(2)等速入水,通过UDF给平板赋予入水速度v,平板匀速下降。v由下式计算:v=√2gh(1)v=2gh−−−√(1)式中:g为重力加速度;h为落高,即平板下表面距初始静水面高度。1.1数学的物理模式1.1.1流体运动方程连续性方程是质量守恒定律在流体运动中的表述:∂ρ∂t+∇⋅(ρ→v)=Sm(2)∂ρ∂t+∇⋅(ρv⃗)=Sm(2)式中:ρ为密度;t为时间;→vv⃗为速度;Sm为源项。运动方程是动量守恒定律在流体运动中的表述:∂∂t(ρ→v)+∇⋅(ρ→v→v)=-∇p+∇⋅(ˉˉτ)+ρ→g+→F(3)∂∂t(ρv⃗)+∇⋅(ρv⃗v⃗)=−∇p+∇⋅(τ¯¯¯)+ρg⃗+F⃗(3)式中:p为压力;ˉˉττ¯¯¯为应力张量;ρ→gg⃗为重力体积力;→F为外部体积力。1.1.2离散方程的定义Fluent软件基于有限体积法进行计算。有限体积法于20世纪60、70年代发展起来,其基本原理为:将计算域划分网格,并使每个网格点周围有一个互不重复的控制体积,将待解的微分方程(控制方程)对每个控制体积积分,从而得到一组离散方程。其中的未知数是网格节点上的因变量。子域法加离散,即是有限体积法的基本思路。1.1.3基于mac法的网格计算采用流体体积法(VOF法)追踪自由液面。该方法由Hirt和Nichols等人基于MAC法得到,其原理为计算网格单元中流体体积和网格单元体积比函数F来确定自由面的形状、位置。若某时刻网格单元中F=1时,则该网格全部被指定相流体占据;若F=0,则网格没有指定相流体;若0<F<1,则网格包含有多相流体交界面。1.2各区域网格及网格划分建立平板仿真模型如图1所示。坐标原点在平板下表面几何中心以下0.3m处,计算域上表面在平板上表面以上0.2m处,平板尺寸为1m×0.8m×0.5m。网格在靠近平板的区域划分较密,其余区域逐步变疏,各区域尺寸及网格数量见表1。动网格:在Fluent软件中,运动物体的模拟通过动网格技术来实现。本文采用的网格运动方式为Layering,即整体网格随平板一起下落,在计算域上表面边界产生新的网格,计算域下表面边界湮灭旧的网格。这样既能满足平板入水时附近有足够多的网格数量来保证计算的精确性,又能减少整体网格的数量加快计算速度。计算域流体属性:各算例根据物体落高不同将计算域划分为空气和水上下两部分。边界条件定义如下:平板表面、计算域四周外表面及下表面定义为壁面边界条件(WALL);计算域上表面定义为压力出口边界条件(PRESSURE_OUTLET);其余定义为内部边界条件(INTERIOR)。1.3平板测试结果平板落高分别为H=0.3m、0.4m、0.5m,对应的入水速度分别为2.426m/s、2.801m/s、3.123m/s。平板下表面上的压力测点分布如图2所示:沿长边方向中心线上布置11个测点,间隔90mm;沿短边方向中心线上布置5个测点,间隔85mm。图3给出了平板以3.123m/s等速入水时各测点压力的时间历程。图4给出了不同落高时各测点压力峰值的计算结果和试验结果(文献中,落高H=0.4m时的W2测点,落高H=0.5m的E2测点,这两个压力峰值数据有缺失),由于对称关系,下面只给出模型一边的计算结果。可以看到:(1)平板各处的砰击压力分布不均匀,其峰值自中心向边缘逐渐减小(试验的压力最大值不在中心,偏向了附近位置,这可能是由于试验模型不是完全的刚体,受到冲击后变形影响了流体载荷的作用);(2)各测点压力到达峰值的时刻较为接近,持续时间也基本相同;(3)比较等速入水和自由落体入水的仿真计算,两者得到的砰击压力峰值基本一致,这是由于结构物质量较大,在入水达到砰击压力峰值之前不会被压力明显地减速从而减小峰值。2三维圆入口本文对四个具有不同底升角的圆锥进行了自由落体入水仿真。2.1建立网格模型分别建立底升角为0°、1°、3°、6°四个圆锥仿真模型。以3°底升角为例,为方便划分结构化网格,取整个计算域为圆柱体,坐标原点在圆锥顶点处,计算域上表面在圆锥上表面以上0.01m处,圆锥上表面半径为0.2032m。模型示意图见图5,各计算域尺寸及网格数量见表2。赋予圆锥质量为28.58kg,与ChuangSL试验模型一致。其余设置同平板入水。2.2覆合物入水过程中面料的动态每个圆锥分别从距静水面高H=0.3810m、0.3048m、0.2286m、0.1524m、0.0762m处落下,对应的入水速度分别为2.734m/s、2.445m/s、2.118m/s、1.729m/s、1.223m/s。图6给出了圆锥入水过程中液面变化情况。图7给出了各圆锥从不同高度下落时受到的表面最大砰击压力峰值,以及ChuangSL根据试验得到的回归结果。可以看到:(1)比较0°底升角圆锥和平板入水,在大致相同的入水速度下,两者表面受到的最大砰击压力峰值基本一致且均在砰击面的几何中心处。这说明砰击面为平面时,其几何外形对最大砰击压力峰值影响不大;(2)对于三维问题,砰击产生的最大压力峰值从0°底升角开始随角度增大而迅速增大,在大约1°底升角达到最大值,然后又随底升角的增大而逐渐减小。而在二维问题(楔形体入水)中,砰击压力峰值的最大值要在3°底升角左右产生。这可能是由于随着底升角的增大,空气在圆锥体下要比楔形体下逃逸得更快,从而使结构物底面空气垫效应迅速减弱。3不同几何形状对最大bm对于三维结构物入水问题,理论研究和试验研究有较多局限。本文基于Fluent软件平台对平板和圆锥体进行了三维数值模拟。结果表明:(1)砰击面为平面时,其几何形