0 引　言

结构物入水砰击过程是一个复杂的强非线性问题，广泛存在于船舶、航空航天、武器装备等领域^[1]。结构物以较高速度入水时会遭受流体强烈的砰击，虽然砰击过程持续的时间较为短暂，但会给结构带来很大的瞬时压力峰值和加速度。

国内外很多学者围绕结构入水砰击过程开展了广泛研究，在公开发表的文献中，Von Karman^[2]和Wagner^[3]是较早开展入水问题研究并建立理论模型的学者。Watanabe^[4]，Howison^[5]等通过匹配渐进法对Wagner理论模型进行修正，对自由液面与物体表接触处的流场进行局部射流分析，在平板假设的基础上引入部分非线性边界条件，针对Wagner计算模型求解物面压力分布时边缘处的奇异性问题，其通过内域解和外域解匹配的方法给出了一种解决方案。Dobrovol'skaya^[6]通过研究二维楔形体零加速度速入水问题，较为完整地给出了对称楔形体匀速入水过程中受力的相似解。Zhao等^[7]运用非线性边界元算法求解了该匀速入水问题，求解精度有所提升。卢炽华^[8]运用边界元算法求解了U形带外飘船体截面的入水砰击问题，研究了典型船体剖面入水过程中水动力特性，为工程设计提供参考。Wu等^[9]进一步用边界元法分析了二维楔形体的自由入水问题，在楔形体尖点触水瞬时取相似解为流域初始解来解决速度、压力等物理量剧烈变化带来的求解发散。Xu等^[10-11]求解了对称楔形体以恒等速度斜向入水的水动力问题，后又将该方法成功拓展到楔形体的多自由度自由入水问题。Sun等^[10]用非线性边界元算法对三维物体入水过程进行研究，以圆锥体为算例给出了较为准确的计算结果。

仿真研究入水问题的常用有限元软件有Ls-dyna，Dyrtan，Fluent，Abaqus等，主流算法有拉格朗日算法、欧拉算法、ALE算法、SPH算法等。陈震^[13]用Dytran仿真求解了平板结构入水砰击问题，研究了入水过程捕捉空气形成的空气垫对平板入水过程的影响。朱显玲^[14]使用Ls-dyna求解了夹心平板结构入水载荷并据此对夹心结构参数进行优化。陈光茂^[15]使用SPH方法对楔形体入水过程进行仿真，分析了斜升角、入水角等参数对砰击载荷的影响规律。云忠^[16]使用STAR-CCM+软件对跨介质飞行器入水过程的动力学性能进行计算，为飞行器结构X设计和入水参数选择提供指导。Wang S和Soares C G用Ls-dyna分别对二维对称楔形体^[17]、二维非对称楔形体^[18]、三维椭球体结构^[19]入水过程进行了仿真求解，对物面压力分布、压力时历、结构整体砰击力、自由液面及压力场进行了有效模拟。曹正林^[20]和段芳^[21]分别用Ls-dyna和Dytran软件求解了三体船结构入水问题，验证了两软件求解负杂曲线结构入水问题的能力。

结构入水过程中涉及自由液面的排挤、扭曲、飞溅、破碎等复杂物理现象，理论和数值计算时通常需要引入一些假设对实际问题进行诸多简化。为得到砰击过程中更为真实的载荷及结构响应规律，并对理论和数值算法解加以验证，落体试验方法在研究过程中还是不可或缺的手段之一。孙辉^[22]针对6个不同角度和厚度组合的楔形体模型开展了一系列落体砰击试验，研究模型入水过程中的加速度和结构响应规律，开展了结构响应对各参数的敏感度分析。Zhao^[23]开展了斜升角30°的V型楔形体入水冲击试验，测量了模型加速度及表面砰击压力，验证了其非线性边界元算法的数值计算结果。莫立新^[24]通过开展4个不同刚度板架结构的落体砰击试验，研究了结构刚度对砰击载荷及结构动响应的影响规律。张哲^[25]在射弹入水试验中初步验证了PVDF压电薄膜测量入水压力的可行性。

本文对求解结构入水问题的基于边界元算法的数值计算方法、基于任意拉格朗日-欧拉耦合算法的数值仿真方法、自由落体入水砰击试验3种常用的研究手段进行了综合对比研究。分别用3种方法对楔形结构自由入水过程进行研究，给出了砰击压力沿物面的分布、典型点处的压力时历、物体运动特性并进行了对比分析。

1 基于边界元算法的数值计算方法 1.1 数学模型

建立如图1所示的数学模型，其中o-xz为模型坐标系，原点位于结构入水前水平面与楔形体对称线的交点处，坐标轴方向如图1所示。W为楔形体实时入水速度，是关于时间t的函数，方向定义为向下为正；s为楔形体尖点在自由液面以下的深度即入水深度，同样是关于时间t的函数；β为楔形体底升角；Ω为整个流体域；S_b为楔形体物面；S_F为自由液面；S_c为远场边界条件。

流体满足无粘、无旋、不可压缩的假设，则流体速度势 $ \phi $ 在整个流体域 $ \varOmega $ 内满足Laplace 方程：

$ {\nabla ^{\text{2}}}\phi {\text{ = 0}}。$

(1)

在楔形物面S_b上满足物面不可穿透边界条件：

$ \frac{{\partial \phi }}{{\partial n}} = - W{n_z} 。$

(2)

式中，n=(n_x，n_z)为物面的单位外法向量。

瞬时自由液面S_F上满足的欧拉形式运动学和动力学自由面边界条件分别表示如下：

$ \frac{{\partial \zeta }}{{\partial t}} = {\phi _z} - {\phi _x}{\zeta _x} ，$

(3)

$ \frac{{\partial \phi }}{{\partial t}} = - \frac{1}{2}\left( {\phi _x^2 + \phi _z^2} \right) - gz 。$

(4)

式中： $ z = \zeta (x,t) $ 为自由液面方程，g代表重力加速度。

假定在无穷远处流场不受扰动，在远场边界S_C上满足如下边界条件：

$ \frac{{\partial \phi }}{{\partial n}} \to 0，\;\;\;\;\;\;\;\sqrt {{x^2} + {z^2}} \to \infty。$

(5)

1.2 拓展坐标系

为保证计算精度和计算速度的协调统一，引入如下拓展坐标系o-αβ：

$ \varphi = \phi /sW \text{，} \alpha = x/s , \beta = z/s , s = \int_{}^{} {W{\rm{d}}t}。$

(6)

在拓展坐标系内速度势 $ \varphi $ 依然满足Laplace 方程，物面边界条件表达为：

$ \frac{{\partial \varphi }}{{\partial n}} = - W{n_\beta } 。$

(7)

式（3）和式（4）中的自由面边界条件变为：

$ \frac{{\partial s\overline \zeta }}{{\partial s}} = {\varphi _\beta } - {\varphi _\alpha }{\overline \zeta _\alpha } ，$

(8)

$ \frac{{\partial sW\varphi }}{{W\partial s}} = - \frac{1}{2}\left( {\varphi _\alpha ^2 + \varphi _\beta ^2} \right) - \frac{{gs\beta }}{{{W^2}}}。$

(9)

其中 $ \overline \zeta $ 是拓展坐标系中的自由液面方程。

远方控制面上依然满足：

$ \frac{{\partial \varphi }}{{\partial n}} \to 0，\;\;\;\;\;\;\;\sqrt {{\alpha ^2} + {\beta ^2}} \to \infty。$

(10)

1.3 控制方程离散

通过格林公式可以将拓展坐标系下的Laplace 方程转化为积分表达形式：

$ A(p)\varphi (p) = \int {\left[ {G(p,q)\frac{{\partial \varphi (q)}}{{\partial {n_q}}} - \varphi (q)\frac{{\partial G(p,q)}}{{\partial {n_q}}}} \right]} {\rm{d}}{s_q}。$

(11)

式中： $ A(p) $ 为源点p处的立体角；G代表格林函数，求解二维问题时表达式如下：

$ G(p,q) = \frac{1}{2}\ln \left[ {{{\left( {{x_p} - {x_q}} \right)}^2} + {{\left( {{z_p} - {z_q}} \right)}^2}} \right]。$

(12)

为了求解式(12)，将流体边界离散成一定数量的线性单元，假设离散为n份，认为单元内速度势线性分布，即可用线性插值函数近似表示离散单元内部的速度势φ和速度势在法线方向的方向导数φ_n。将它们的线性插值函数代入方程(12)，整理得到：

$ \left[ {{{\boldsymbol{H}}}} \right]{}_{n' \times n'} \cdot {[\varphi ]_{n' \times {\text{1}}}}{\text{ = }}\left[ {{{\boldsymbol{G}}}} \right]{}_{n' \times n'} \cdot {[{\varphi _n}]_{n' \times {\text{1}}}} 。$

(13)

式中：n为整个边界上划分节点个数， $\left[ {{{\boldsymbol{H}}}} \right]{}_{n' \times n'}$ ， $\left[ {{{\boldsymbol{G}}}} \right]{}_{n' \times n'}$ 为式(12)离散后的系数矩阵。

2 基于ALE算法的数值仿真方法 2.1 算法简介及基本控制方程

ALE算法是计算流固耦合问题时比较适用的算法之一，其综合吸收了2种算法的优势又规避了单一算法的劣势^[26]。

定义材料移动速度为v，网格移动速度为u，则w=v-u表示相对运动速度。ALE算法满足如下守恒方程^[26]：

质量守恒方程

$ \frac{{\partial \rho }}{{\partial t}} = - \rho \frac{{\partial {v_i}}}{{\partial {x_i}}} - {w_i}\frac{{\partial \rho }}{{\partial {x_i}}}，$

(14)

动量守恒方程

$ \rho \frac{{\partial {v_i}}}{{\partial t}} = \frac{{\partial {\sigma _{ij}}}}{{\partial {x_j}}} + \rho {b_i} - \rho {w_i}\frac{{\partial {v_i}}}{{\partial {x_j}}}，$

(15)

能量守恒方程

$ \rho \frac{{\partial E}}{{\partial t}} = {\sigma _{ij}}\frac{{\partial {v_i}}}{{\partial {x_j}}} + \rho {b_i}{v_i} - \rho {w_j}\frac{{\partial E}}{{\partial {x_j}}}。$

(16)

式中：ρ为流体密度；b为作用在流体上的单位体积力；E为流体比内能；σ_ij为应力张量。

2.2 软件计算流程

Ls-dyna提供了ALE算法平台，本文仿真使用MSC.PATRAN进行建模，使用LS-Prepost进行关键字修改、耦合定义及后处理，提交Ls-dyna求解器进行计算。

Ls-dyna中通过空材料加状态方程来描述水和空气2种流体材料的物理及力学特性^[27]。

理想气体用关键字*EOS_LINAR_POLYNOMIAL定义的线性多项式状态方程模拟，表达式如下：

$ P = [{C_0} + {C_1}\mu + {C_2}{\mu ^2} + {C_3}{\mu ^3}] + [{C_4} + {C_5}\mu + {C_6}{\mu ^2}] \cdot {\rho _0}e。$

(17)

式中：C₀～C₆为多项式拟合系数；ρ₀为空材料*MAT_NULL关键字里定义的材料参考密度；ρ为实时密度， $\mu = \dfrac{\rho }{{{\rho _0}}} - 1$ ；e表示流体比内能。

水用*EOS_GRUNEISEN关键字中定义的Gruneisen状态方程描述，表达式如下：

$ \begin{split}\left\{\begin{array}{l}P=\dfrac{{\rho }_{0}{C}^{2}\mu [1+(1-\dfrac{{\gamma }_{0}}{2})\mu -\dfrac{a}{2}{\mu }^{2}]}{[1-({S}_{1}-1)\mu -{S}_{2}\dfrac{{\mu }^{2}}{\mu +1}-{S}_{3}\dfrac{{\mu }^{3}}{{(\mu +1)}^{2}}}+\\ \qquad({\gamma }_{0}+a\mu ){\rho }_{0}e，\quad\quad\quad\quad\quad\,\,\,压缩状态，\\ P={\rho }_{0}{C}^{2}\mu +({\gamma }_{0}+a\mu ){\rho }_{0}e，\quad\quad \,膨胀状态。\end{array} \right. \end{split}$

(18)

式中：ρ₀为水密度；C为水中声速；S₁，S₂，S₃，γ₀为状态方程系数系数；a是对γ₀的一阶体积修正。

为与其他2种方法的结果进行对比并节省计算时间，本文计算时将入水结构定义为刚体材料（*MAT_RIGID）。结构物与流体的耦合关系通过关键字*CONSTRAINED_LAGRANGE_IN_SOLID来实现的，通过罚函数法实现流固界面力学参量的传递，建模时需注意单元法向量的统一。这种耦合定义方式建模时不需要固体网格与流体网格间共节点，这极大地降低了建模的复杂程度。

在Ls-dyna中建立楔形体入水问题数值仿真模型如图2所示，采用二维仿真方式，在楔形体长度方向只设置一层单元。入水结构物用Lagrange单元离散，赋予刚体属性；流体域离散为欧拉实体单元并划分为两部分，分别赋予空气、水属性，2种流体交界面上共节点。为在保证计算精度的前提下提升计算速度，水和空气在远离核心区的范围网格划分逐渐变粗，核心区的网格密度最终由求解前开展的网格收敛性分析确定。同时利用本文求解问题的对称性，仅建立一半模型进行仿真计算，在对称面上建立对称面边界条件。初始时刻楔形体位于水面以上较小高度处，根据自由落体计算给其以一定的初速V₀便可模拟试验中不同高度工况下的自由落体入水过程。

边界条件设置：二维仿真中约束了欧拉单元节点在图2坐标系y轴方向的运动。流体域除对称面外的其余3个面施加无反射边界条件以免反射回波影响计算结果，对称面上设置刚性墙并限制刚体的运动自由度，计算过程中考虑了重力场的影响。

3 落体砰击试验方法

落体入水试验是将试验件提升至水平面以上一定高度后释放，在重力作用下试验件获得所需的入水速度，与水发生相互作用来模拟砰击现象。本次试验在长宽深为108 m×7 m×3.5 m的水池上开展，其上设计有一套完整的落体试验装置实现试验件起吊、释放并保证所需入水姿态。本文试验的楔形体模型由某型舰船舷侧板架结构缩比设计而来，试验模型底升角为45°，长宽高为1.5 m×0.9 m×0.75 m，整个模型（含配重）总重553 kg，在模型中横剖面自上而下布置了5个压力测点，模型及测点布置局部放大图见图3。试验中通过压力传感器监测入水过程中结构底部外表面受到的砰击压力，通过拉绳式的位移传感器测量模型入水过程中的位移变化过程，取微分可得到模型入水过程中的速度。试验中每个工况至少重复3次以验证试验开展过程的稳定性。

4 结果对比与分析

分别使用3种方法对图3所示的楔形结构入水载荷和运动特性进行计算。给出入水结构与水平面接触时刻速度V₁=4.427 m/s工况下的结果对比，并将此时刻定为时间零点，P₁～P₅五个测点的具体位置见图3（b）。

将3种方法得到的各测点砰击压力信号统一时间轴绘制得到图4，可以看出Ls-dyna仿真、BEM算法和落体入水试验3种方法得到的结果一致性良好，砰击压力信号整体趋势、持续时间、压力峰值等关键参数都较为一致，具体的压力峰值差异在6.5%以内。

用3种典型方法对入水过程中的结构速度进行求解或测量，结果见图5(a)。可见在0.075 s之前三者吻合很好，0.075 s之后边界元求解结果开始偏离Ls-dyna仿真结果和试验数据。这是由于在边界元算法求解时引入了物面无限大的假设，忽略了射流越过楔形体“肩部”后真实存在的复杂的流动分离现象，图5(b)中的加速度曲线表现出了同样的规律。Ls-dyna软件为复杂结构入水问题提供了方便的建模、求解方式，相比边界元算法对所求解入水结构复杂程度的包容性较强。

图6给出了边界元算法得到的P₃点砰击压力达峰时刻的楔形体位置、自由液面形状、砰击压力物面分布情况。可见此时刻入水深度s=0.174 9 m，P₃点处于未受扰动的水平面之上，物面压力在楔形结构物尖点和入水激起的射流根部取得极值，其中第1个极值中包含了楔形体尖点的静水压力，第2个极值对应的物面尚处于水面之上，体现的全部是砰击带来的的物面载荷。第2个极值之后物面压力逐渐下降，在楔形体半宽L等于（1.4～1.6）倍入水深度s范围内迅速下降，当L=0.295 4 m（约1.7倍入水深度）时压力下降为0，即砰击载荷作用在半宽L小于1.7倍入水深度的物面范围内。

图7给出了3种方法得到的楔形体入水0.05 s（入水深度s=0.214 7 m）时刻的物面压力分布情况，整体来看3种方法求解结果一致性较好。边界元算法中通过拓展坐标变换的方法，网格划分更为灵活，可得到比较光滑的物面压力曲线；Ls-dyna仿真过程中需要在结构网格上分布设置足够多的压力监测点，在各点压力时历数据中提取关注时刻的砰击压力，进而拟合出关注时刻的物面压力分布；落体试验中则是通过在试验件上布设的点式压力传感器测量得到，要想通过试验方法准确拟合出砰击压力沿物面的分布则需要布置大量的点式压力传感器，试验成本会大幅上升。

图8给出了6个典型时刻Ls-dyna仿真方法和边界元模拟方法得到的自由液面发展情况对比，t₁～t₅分别是P₁～P₅监测点压力达到峰值的时刻，t₆为入水深度为0.35 m的时刻。从前5个时刻可见对于对砰击载荷影响较大的射流根部及以下区域，2种方法的结果较为一致，从第6个时刻可见当射流根部越过楔形体肩部位置时2种方法的仿真差异变大，原因如前所述是因为边界元算法中未考虑流场越过楔形体肩部时复杂变化，从图9试验高速摄像也可以看出，射流到达肩部之前流动比较规则，越过肩部以后，射流开始出现破碎、飞溅，这给边界元数值方法求解带来困难，相对而言Ls-dyna采用的ALE算法可以更为真实的模拟这种现象。

图10为Ls-dyna仿真得到的典型时刻流体域压力场及速度场云图。可以看出入水过程中流体域中压力较大的区域集中在物面上，流体速度较大的位置集中在射流区域。

5 结　语

本文介绍了3种求解结构入水问题的典型技术途径，通过对比分析得出以下结论：

1） Ls-dyna仿真方法和边界元数值模拟方法均可对楔形结构入水砰击问题进行有效求解，通过与试验结果对比，验证了2种方法求解入水载荷及运动特性的精度。

2） Ls-dyna仿真方法和边界元数值方法均可有效模拟出楔形体结构入水后的自由液面变化情况，Ls-dyna可对流动分离、复杂形状结构入水等数值方法较难处理的问题给出满意的求解。

3） 底升角为45°的楔形结构入水过程中砰击载荷作用范围为楔形半宽<1.7倍入水深度的区域。

4） 随着入水深度的增加，射流根部逐渐“健壮”，表现在压力上是从下至上各测点取得峰值时处于静水平面之上的高度逐渐增大，但各测点压力峰值由于结构入水速度的降低会逐渐降低。