基于显式有限元方法的二维楔形体砰击压力系数预报

许蕴蕾

（海军驻上海地区舰艇设计研究军代室 上海200011）

基于显式有限元方法的二维楔形体砰击压力系数预报

许蕴蕾

（海军驻上海地区舰艇设计研究军代室 上海200011）

利用LS-DYNA软件对二维楔形体入水问题进行研究。基于显式有限元方法，选用任意拉格朗日-欧拉算法，建立了包含空气、水和楔形体的完全耦合的二维有限元模型，研究流场的射流现象与压力变化情况，预报了二维楔形体砰击压力系数，并与已公开发表的模型试验结果吻合较好，从而为后续砰击载荷计算提供可靠的方法。

LS-DYNA软件；二维楔形体；拉格朗日-欧拉算法；砰击压力系数

引 言

船舶在海上航行时，经常会遇到恶劣海况，砰击问题越来越受到学者们的关注。船舶遭受砰击的瞬时，底部受到巨大的冲击力，船体的垂向加速度会突然改变并随即出现高频振动，严重时会造成船体结构的破坏。因此对于船舶设计工作者来说，必须准确预报船体的艏部砰击压力。

砰击问题是一个包含结构、空气和水三种介质耦合的非线性定常问题。自从Von Karman和Wagner开展研究以来，针对砰击问题，人们开展了大量的理论、数值和试验研究工作。Zhao等人［1］研究了边界元方法，预报了二维楔形体和艏外飘入水砰击问题，并开展了试验对比分析工作。骆寒冰［2］采用LS-DYNA软件，对二维楔形体入水问题进行数值预报，并与Zhao等人的模型试验结果进行了比较，结果吻合较好。陈震［3］给出二维楔形体入水的基本理论，并采用有限元软件MSC-Dytran对二维楔形体的入水过程进行仿真模拟，获得砰击面上压力分布以及入水角度对砰击载荷的影响规律。

本文选用任意拉格朗日-欧拉算法，用拉格朗日网格描述结构，多物质欧拉网格描述空气和水。采用显式有限元方法、罚函数耦合算法，来模拟欧拉流体与拉格朗日固体接触面之间的耦合效应。二维楔形体模型，利用LS-DYNA软件，对斜升角为10°的刚性楔形入水问题进行仿真模拟，将数值结果与文献［4］中计算结果进行对比，验证LS-DYNA软件在流固耦合分析中的可靠性。同时，本文提出在楔形体入水问题中利用ANSYS前处理实现位移、速度曲线控制的方法，讨论了楔形体入水时液面变化情况和喷射区的特点以及流场压力变化情况，重点讨论了斜升角为30°的刚性楔形体入水问题，与Zhao等人的模型实验结果和理论计算结果进行了对比，得到了砰击压力系数。

1 ALE法

ALE法可以处理整个物体有空间的大位移并且本身有大变形等问题，ALE法最早出现于数值模拟流体动力学问题的有限差分方法中。在内部网格的划分上，吸收了欧拉网格的长处，又不完全相同，ALE法不仅使内部网格单元独立于物质实体而存在，还可以根据定义的参数在求解过程中适当调整网格位置，避免出现严重的畸变。此法可使网格与网格之间的物质流动。与欧拉网格一样，ALE网格有两层网格重叠在一起，空间网格可以在空间任意运动，其余与欧拉方法描述一样，在两层网格中进行物质输送。在结构边界运动的处理上与拉格朗日方法的特点类似，ALE法能够有效的跟踪物质结构边界的运动。

ALE法执行自动重分区时，采用拉格朗日和欧拉两种算法，包含伴随重映射的拉格朗日时间步。ALE法可以实现流体-固体耦合的动态分析，可以克服单元严重畸变引起的数值计算困难。ALE法先执行一个或几个拉格朗日时间步计算，此时材料的流动使单元网格产生变形，接着执行ALE时间步计算：

（1）变形后的物体边界条件保持不变，充分网格内部单元，网格的拓扑关系保持不变；

（2）将变形网格中的单元变量和节点速度矢量输运到重分后的新网格中，称为Advection Step。

2 二维楔形体入水问题仿真分析

2.1 建立仿真模型

图1为楔形体的砰击物理模型，楔形结构为二维对称结构，只在z轴方向运动。将静水面处取为xoy平面，在y轴方向取单位长度。β为楔形体入水角，Ωfluid代表液体欧拉域，Ωair代表空气欧拉域，ve为入水速度。

图1 楔形体砰击计算模型

楔形结构用壳单元描述，空气和水的流域用六面体欧拉单元表述。在欧拉域中，不同位置的网格大小有变化。越靠近边界处的网格划分越粗，而楔形体结构位置附近的网格划分则较细。选取楔形结构的外侧作为流固耦合作用面，楔形体的外表面法向量调整为指向流体。为了形成一个封闭的耦合面，在建模时选取对称的楔形结构的1/2部分，并将对称面处的结构单元定义为虚拟单元。楔形体只有一个自由度在z轴方向，模型只划分一个网格沿y方向，并在该方向的单元表面设置刚性墙边界条件模拟二维无限大流场。下页图2为楔形体部分网格模型。楔形体下落初速度为6.15 m/s。

图2 楔形体部分网格模型

2.2 状态方程

本文中流体的压力采用了线性多项式状态方程来描述，线性多项式状态方程表示单位体积内能的线性关系，压力值由式（1）给定：

式中：C0、C1、C2、C3、C4、C5和C6为常数。如果式中μ＜0，则C2μ2和C6μ2两项设置为0；其中μ=1/V-1，V表示相对体积。

线性多项式状态方程可以用于模拟符合γ律状态方程的气体，各项系数可设置为C0=C1=C2=C3= C6=0，C4=C5=γ-1（γ为单位热值律）。

表1 线性多项式状态方程中空气和水的参数值

2.3 砰击压力系数

基于冲量砰击理论，楔形面处的砰击力峰值可以用式（2）表示：

因此，无因次压力系数k表达式为：

式中： ρ为水的质量密度；V为楔形面在法向的速度分量；k为无因次化的压力峰值系数。

3 DYNA流固耦合分析要点

3.1 罚函数耦合

为检查每一个节点对主表面的穿透，用罚函数耦合系数追踪结构从物质-拉格朗日节点和欧拉流体间的相对位移d。罚函数的物理意义相当于为了限制从节点和主面之间的穿透，在节点和被穿透面之间放置一个法向弹簧，如果不出现穿透，就不进行任何操作。如果有穿透发生，即从节点对主物质表面的穿透，则从节点与被穿透表面之间引入一个大小与穿透深度成正比的较大界面接触力。罚函数数值的大小不仅影响穿透现象的发生，如果罚函数数值太大，还会使物质产生钢化现象，物质表面过分刚硬会影响程序计算的稳定。

3.2 刚性体模型速度位移曲线

刚性体模型在显示动力学分析中具有非常重要的意义。当定义一个刚体后，刚体内所有节点的自由度都耦合到刚性体的质量中心上，每个时间步的节点力和力矩合成为作用在刚性体上的力和力矩。然后计算刚性体的运动，再转换到节点位移。由于无论定义多少个节点，刚体仅有6个自由度，因此将有限元模型中刚硬部分用刚性体模型定义，可以大大缩减显示分析的计算时间。

定义刚体模型可约束平移参数和转动约束参数，所有刚体的加载必须在前处理ANSYS中定义的part号上而不能使用在组件上。在刚体上施加位移条件，相当于在DYNA中定义位移曲线。在前处理ANSYS中定义数组参数位移曲线是将时间间隔和对应的位移值集合在一起。它可以分成两部分，一部分为时间间隔，另一部分为与时间间隔对应的位移值，定义位移曲线之前，必须首先定义位移曲线的纵坐标和横坐标值所对应的数组，且其维数必须相同。因此，在计算楔形体以恒定速度入水时，可根据此方法控制楔形体速度保持恒定。

3.3 压力传感器

在DYNA中对刚性体模型的单元进行压力的提取，需要定义压力传感器，在输入文件中定义压力传感器位置，这些位置同样是拉格朗日算符单元所在的位置，压力传感器压力的提取是从流体单元中得到的压力。

3.4 时间步长尺度

在LS-DYNA中自动选择的时间步长是基于以下条件的：在一个单位时间步内声波不允许穿过最小的单元，时间步长也受到全系统刚度的限制，然后与时间步长相关的过大的接触刚度可能会产生非物理的能量的接触，时间步长影响了震荡压力信号，其作用效果与较高的接触刚度的作用效果一样，时间步长的改变对结果没有明显的影响。

4 计算结果

4.1 10°斜升角楔形体计算结果比较

文献［4］中二维刚性楔形体斜升角为10°，该模型采用刚性材料以5.425 m/s的速度等速入水。本文利用LS-DYNA软件进行仿真，利用前处理ANSYS中施加位移曲线的方式实现等速入水，得到如表2所示，砰击压力的理论值与数值解和对应位置本文的仿真值之间的比较结果。本文计算结果误差较小，计算结果较为可靠。图3为砰击压力时历曲线。

表2 文献值与本文仿真值之间的比较

图3 砰击压力时域结果

4.2 30°斜升角楔形体仿真模型计算结果

4.2.1 压力分布

图4 -图6显示了入水砰击过程中不同时刻的水面飞溅现象。图7 -图9显示流场中压力分布情况。在该计算模型中并没有观察到液体的渗漏现象，表明耦合模拟合理。可以看出，在射流处，砰击压力急剧减少，可以忽略不计；在外域水中，离液面升高处越远，压力越小。液面分离之前，砰击压力的最大值出现在液面升高处。在液面分离之前的入水过程，砰击压力的最大值明显降低，砰击压力的最大值出现的位置，将从楔形体底部向上端折角处移动。这是由于仿真楔形体质量较轻，速度下降较快，导致砰击压力逐渐减少，也与压力传感器得到的压力相对应，位置越靠上的压力传感器测得的压力最大值越小。

图4 入水初期

图5 入水中期

图6 入水后期

图7 0.051 s 时流场压力分布

图8 0.117 s 时流场压力分布

4.2.2 时域内砰击压力曲线

图10显示在对应于Zhao等［1］模型试验中，压力传感器位置处的仿真模型压力传感器的压力时历曲线。最大砰击压力的逐渐降低，同时也对应于图7 -图9显示的砰击压力最大值的降低。

图10 楔形体入水实验模型

图11所示为5个压力传感器测得的砰击压力时域预报结果。横坐标是入水时间，5个压力传感器位置对应于图10所示试验压力测点位置。

图11 不同位置砰击压力时历曲线

4.2.3 砰击压力系数

文献［1］中给出不同斜升角楔形体的不同理论算法算得的砰击压力系数，对于30°斜升角的砰击压力系数曲线如图12所示。

图12 砰击压力空间分布比较

本文计算结果与文献［1］ 中模型试验结果吻合较好，图中所示4个点的数据来源于前4个压力传感器。其中，文献［1］ 给出的理论计算结果是样条曲线的形式，是由于理论计算中选取计算点数目较多，将一系列计算点结果进行曲线拟合，便得到图12中的计算结果。

本文利用LS-DYNA软件进行仿真模拟，考虑到计算时间问题，选取4个点的结果进行比较，因此本文计算给出的结果是独立点，这些独立点的计算结果与理论值吻合较好。其中文献［1］ 给出的结果中，对纵坐标做了无因次化处理，Z是楔形体液面距离楔形体底部的垂向高度，V表示楔形体在入水不同时刻的垂向速度。在砰击过程中，LS-DYNA预报的砰击压力系数与试验值吻合非常好，反映显式有限元技术在预报二维刚性体砰击问题上的可靠性以及本文所采用仿真方法的合理性。

5 结 论

本文采用LS-DYNA软件，基于显式有限元方法，采用任意拉格朗日-欧拉算法和罚函数耦合算法，针对文献［1］ 中楔形体入水砰击试验结果以及理论计算结果，建立二维刚性楔形体模型，展开砰击压力系数预报工作。预报的10°楔形体砰击压力值与文献［4］中吻合较好，验证了LS-DYNA在二维楔形体入水仿真中比较可靠。

文中重点分析与Zhao等［1］楔形体入水砰击实验相同30°斜升角的楔形体入水模型，得到不同测点的砰击压力时域分布以及求得砰击压力系数，并与其理论推导结果进行比较，与模型试验以及理论推导结果吻合较好。预报的砰击压力结果可正确描述入水砰击的整个过程（包括入水初期、入水中期的砰击过程阶段）。LS-DYNA有效预报了砰击压力在楔形体表面的时间、空间分布，合理反映出砰击过程的射流以及内域和外域中的压力分布，为后续计算砰击压力提供了砰击压力系数求解的方法。

［1］ ZHAO R，FALTINSEN O M.Water entry of twodimensional bodies［J］. Journal of Fluid Mechanics，1993，246：593-612.

［2］ 骆寒冰，吴景健，王珊，等.基于显式有限元方法的二维楔形刚体入水砰击载荷并行计算预报［J］.船舶力学，2012，16（8）：907-904.

［3］ 陈震，肖熙.二维楔形体入水砰击仿真研究［J］.上海交通大学学报，2007，41（9）：1425-1428.

［4］ AQUELET N，SOULI M，OLOVSSON L. Euler-Lagrange coupling with damping effects：Application to slamming problems［J］.Computer methods in applied mechanics and engineering，2006，195（1-3）：110-132.

Slamming pressure coef fi cient prediction for water entry of 2D wedge body with explicit fi nite element method

XU Yun-lei
(Representative Of fi ce of Naval Warship Design & Research, Shanghai 200011, China)

The water entry of a two-dimension wedge body has been studied by LS-DYNA software. Based on the explicit finite element and Arbitrary Lagrangian-Eulerian (ALE) algorithm, it builds up a two-dimension coupling fi nite element model including air, water and wedge, studies jet fl ow and pressure change, and predicts the slamming pressure coef fi cient of the two-dimension wedge body. The results are consistent with the published model test results, which can provide the reliable method for calculating of the slamming pressure coef fi cient.

LS-DYNA software; 2D wedge body; Arbitrary Lagrangian-Eulerian (ALE) algorithm; slamming pressure coef fi cient

U661.44

A

1001-9855（2014）03-0042-06

2014-02-18 ；

2014-04-16

许蕴蕾（1981-），女，硕士，工程师，研究方向：船舶总体设计。