基于谱单元方法的平板冲击流固耦合研究

郭 君，张阿漫，杨文山，李 佳

(哈尔滨工程大学 船舶工程学院，哈尔滨 150001)

水下爆炸载荷作用下的舰船结构毁伤一直是国内外学者研究的热点［1－3]，与之相对应的水下冲击流固耦合问题主要是求解瞬态冲击波激励下结构的响应，其研究内容涉及到结构、材料、流体力学等诸多领域。当结构遭受水下冲击波载荷作用时，可压缩瞬态波在流体中传播并撞击到结构上，在结构湿表面处会发生绕射，如果结构是弹性的，且周围水的静压力足够低，散开的稀疏波可能会导致子区域内产生气穴，我们称之为气穴现象［4]。产生气穴区域的水的体积模量比未产生气穴区域的水的体积模量要小好几个数量级，呈现非线性特性。由于结构尺寸、入射波衰减长度和气穴区域在空间尺寸上差别很大，气穴会对近自由面结构如水面舰船产生很大影响，如果一个45 kg的HBX－1炸药在水下爆炸，产生气穴区域，产生的气穴区域在水平上最长为489 m，深度上最大为9 m，其中船体尺寸为:宽9.1 m、高4.6 m、吃水 3.1 m。

为处理气穴现象，Newton［5]提出了一种非常有效的计算方案来处理气穴声学流体，该方案采用位移势作为流体域有限元方程中的主要变量，位移势是一个标量。由于传统流体有限元方程中采用的变量为位移矢量，位移矢量表达式中的变量是位移势表达式中的未知量数目的三倍，而且不能自动对流体运动实行无旋处理，因此，与位移矢量相比，位移势的选取具有相当大的优势。Newton方案的第二个特点是应用了显式时间积分方法。后来，Felippa和 DeRuntz［6]对 Newton的方法进行了扩展，对流固耦合方程进行三维表述，并对结构和流体进行有限元处理。然而，有限元方法采用的单元通常为低阶单元，在解决波的传播问题时会出现数值发散，要想得到准确的数值结果需要十分精细的网格，尤其当有气穴现象出现时，网格细化更是尤为重要，从而导致计算量很大。为此，基于Felippa和DeRuntz的三维表述和结构有限元处理，本文引进谱单元对流体进行处理，以探索一种更为有效的水下爆炸流固耦合求解方法。文中采用基于勒让德多项式的谱单元［7－8]，这是一种高阶单元，采用高阶单元可以不需要细化网格便获得收敛性。此外，谱单元还可以很好地显示波的传播特性［9－10]。谱单元方法综合了有限元方法的几何灵活性和谱方法的计算高精度性，在保证精度的同时还能提高计算速度。

1 谱单元数值模型

1.1 控制方程

在推导产生气穴的声学流体的控制方程时，首先采用双线性模型［11]来考虑流体中的气穴现象。双线性模型是一种简单合理的流体构成模型，在这种模型中，当压缩率(负体积应变)为正时，体积模量为水作为声学介质时的值;当压缩率为负时，体积模量为零。双线性状态方程如下:

式中，Ba是不存在气穴时的体积模量是流体的总压力，pv是温度T0时水的蒸汽压力，ρ0是温度T0饱和流体的密度是瞬时密度是流体质点的位置矢量，t表示时间。假定在整个讨论域内T0、ρ0和ρv是常量，定义，其中c是无气穴声学流体中的声速，S=(ρ/ρ0－1)是总的压缩率，通过这些定义可将状态方程改写为:

假设流体是无旋的，且忽略粘性的影响，从而保证了流场用势来进行表述。用位移势来表示流场的位移为，其中和分别为流体质点的总位移和平衡位移为位移势是流体质点的整体位置矢量。

根据流体的连续性方程和欧拉运动方程，并引用位移势，可得到流体的平衡方程如下:

其中peq=－ρ0(Vg+h)是平衡压力。

式(2)与式(3)便构成了流体的控制方程组，利用关系式peq=pv+ρ0(Vg+h)，并定义压缩率s=ρ0(SSeq)，位移势 ψ =ρ0Ψ，则根据式(2)和式(3)，控制方程组可写为:

式中为动态压缩率为动态位移势为流体节点坐标，peq为水的平衡压力，c为水中的声速，pv为水的汽化压力，由于pv比平衡压力小很多，故可忽略不计，上标点表示变量对时间的导数。

1.2 谱单元空间离散

与有限元离散相同，采用等参元对整个流体域进行离散，将其离散成ne个六面体单元，则每个单元的几何形状可表示为:

φ仍然采用三线性函数。

然而，与有限元离散不同的是，在自然坐标系内，将压缩率s和位移势ψ离散为:

式中:φ为由N阶多项式构成的基函数，s为每个单元内(N+1)3个节点分别对应的压缩率，ψ为每个单元内(N+1)3个节点分别对应的位移势，ξ，η，ζ为单元的自然坐标(－1≤ξ，η，ζ≤1)。

谱单元方法和有限元方法的最大区别在于φi的选择，在有限元法中，φi取为与几何变换相同的三线基函数，而在谱单元方法中，φi取为由勒让德多项式构成的一维拉格朗日插值函数［12]，这也是目前国内外应用最广泛最有效的谱单元基函数:

式中 PN为 N阶勒让德多项式，ξi为第 i个 Gauss-Lobatto-Legendre积分点，其值为下式的根:

可以证明，表达式(6)满足下列关系式:

式中δij为克罗内克δ。

用标准伽辽金法对控制方程式(4a)进行离散，用φ左乘式(4a)，在流体域上对方程积分，并应用格林第一公式可得:

式中Q，H，b分别为:

2 平板的水下爆炸流固耦合研究

2.1 计算模型

为验证所建立的数值模型的准确性，需将数值解与准确解进行对比。这里采用的计算模型为经典的Bleich-Sandler平板模型，如图4所示，平板位于半无限长流体柱上，遭受平面步指数冲击波作用，平板速度响应的解析解可由特征法得到［11]。

图1 Bleich-Sandler平板模型Fig.1 Bleich-Sandler plate model

计算模型的物理属性为:大气压力为Patm=0.101 MPa，重力加速度为g=9.81 m/s2，流体密度为 ρ=998 kg/m3，流体中声速为 c=1 450 m/s，流体长度为3.81 m，由100个等边长的六面体单元构成，平板密度为ρ=5 698 kg/m3，边长为0.038 1 m，由一个单元构成，平面步指数波的压力峰值为p0=0.71 MPa，延迟时间为τ=0.1 ms，在冲击波波阵面距离平板一个流体网格时开始积分计算［13]。

2.2 计算结果及讨论

平板中心的速度时历曲线如图2所示，实线显示的是由谱单元法得到的数值解，圆点和三角形点则分别表示的是由特征法得到的考虑气穴和不考虑气穴时的平板中心速度响应的解析解。分别计算了人工阻尼β 取作0.25、0.50、0.75 和 1.0 的四种情况，临界时间步根据2.3节中的表述进行计算，计算时间步取作临界时间步的0.75倍。

图2给出了人工阻尼β取不同值时平板中心的速度随时间变化的曲线，对计及气穴和不计气穴的两种情况进行了计算，并分别与Bleich-Sandler的解析解进行对比。从图2可以看出，由谱单元方法得到的数值结果与Bleich-Sandler解析解吻合较好，同时人工阻尼β对结构的响应幅值和响应时间几乎没有影响，图2(a)中曲线的微小振荡是由流体网格的瞬时扰动所引起，随着β的增加，振荡逐渐消失。从图2还可以看出，气穴的存在不会改变速度的峰值，但会延长速度的衰减时间，从而延长载荷对结构的作用时间，对结构产生更大的破坏。

3 谱单元方法和有限元方法的对比研究

如前面部分提到，谱单元法相比有限元法具有很多优点，这节将分别从精度和效率上对此进行讨论。

3.1 计算模型

如图3所示，计算模型为位于流体柱上的一个两自由度的质量弹簧振子，遭受平面步指数冲击波作用，图中质量块m1表示船体外板，质量块m2表示船体内部的结构和设备，两质量块由刚度系数为k的弹簧连接，m1和m2的位移分别用u1(t)和u2(t)表示，速度分别用V1(t)和V2(t)表示。

计算模型的物理属性为:大气压力为Patm=0.101 MPa，重力加速度为 g=9.81 m/s2，质量块由边长为0.3 m的正方形板表示，m1=76.9 kg，m2=384.5 kg，弹簧刚度系数为k=94 870 kg/s2，流体密度为ρ=1 026 kg/m3，流体中声速为c=1 500 m/s，流体长度为3 m，压力峰值为 p0=16.2 MPa，延迟时间为 τ=0.42 ms，冲击波波阵面距m1的距离为dinc=0.3 m时开始进行积分计算。

图2 平板中心速度随时间的变化曲线Fig.2 The flat center speed changing with time curve

3.2 计算结果及讨论

有限元方法作为一种数值方法，已经被广泛应用于水下冲击问题的求解中，并证明具有一定的精度［6]，针对3.1 节的物理模型，我们采用由高度细化的有限元网格(24 000个等长单元)产生的有限元解作为基准解，讨论谱单元方法(SEM)和有限元方法(FEM)在不同细化程度时的精度。

图3 计算模型Fig.3 Calculation Model

分别采用谱单元方法和有限元方法计算m1的速度响应，并分别与基准解对比。不同细化程度时质量块m1的速度随时间的变化曲线如图4所示。

图4给出了分别采用谱单元方法和有限元方法在不同细化程度下算得的质量块m1的速度响应曲线，分别将有限元解和谱单元解与基准解进行对比分析。误差系数见表1。

表1 有限元解及谱单元解与基准解之间的误差系数cTab.1 Finite element solution and the spectral element solution and the error c between the coefficient of benchmark solutions

图4 结构的速度响应对比曲线Fig.4 The comparison curve of the velocity response of the structure

从表1给出的误差系数可以看出，谱单元方法和有限元方法所产生的误差均随着细化阶数的升高逐渐减小，收敛于基准解。然而，在细化阶数相同(系统所具有的自由度相同)的条件下，谱单元方法得到的数值解更接近基准解，因此，与有限元方法相比，谱单元方法具有更高的精度。

从表1还可以看出，在误差相同的条件下，谱单元方法需要的自由度较少。例如:若将c=0.15设为允许误差，用谱单元方法需要不足496个自由度便能达到允许误差，而用有限元方法需要的自由度数则接近1 025。因此，谱单元方法在精度提高的同时可以大大减少计算所需要的自由度数，具有更高的计算效率。

4 结论

本文基于流体双线性构成关系考虑流体中可能发生的气穴现象，采用基于勒让德多项式的谱单元方法建立了水下爆炸流固耦合动力学模型，应用经典的平板冲击问题对模型进行了验证。分别采用谱单元方法和有限元方法对弹簧－平板模型进行水下爆炸流固耦合的求解，基于广义误差系数对谱单元方法和有限元方法进行对比研究，得到以下结论:

(1)采用谱单元方法对Bleich-Sandler平板水下爆炸流固耦合问题进行计算，得到的数值解与Bleich-Sandler的基本解吻合良好，谱单元方法能较好地应用于水下爆炸流固耦合问题的求解中。

(2)流体中产生气穴时，未产生气穴的流体区域将结构与产生气穴的流体区域分开，形成流体堆积。气穴区域闭合时流体中形成正压力冲击作用，产生二次加载现象，从而对结构产生更严重的破坏，因此在计算中应予以考虑。

(3)人工阻尼对结构响应具有光顺作用，但基本不会影响结构的响应幅值和响应时间。随着人工阻尼的增加，由流体网格瞬时扰动引起的结构响应的振荡逐渐消失。

(4)与有限元方法相比，在同样网格细化的条件下，谱单元方法具有更高的精度;在误差相同的条件下，谱单元方法所需的自由度数大大小于有限元方法所需的自由度数。因此，与有限元方法相比，谱单元方法在提高精度的同时，能大量节省计算时间。

［1] 牟金磊，朱 锡，黄晓明，等.水下爆炸气泡载荷在加筋板塑性变形中的作用［J] .振动与冲击，2010，29(5):74－77.

［2] 李海涛，朱 锡，赵小龙，等.箱形梁在水下近距非接触爆炸作用下的整体毁伤研究［J] .振动与冲击，2010，29(3):158－161.

［3] Ramajeyathilagam K，Vendhan C P，Bhujanga Rao V.Nonlinear transient dynamic response of rectangular plates under shock loading［J] . International Journal of Impact Engineering，2000，24:999 －1015.

［4] Kennard E H.Cavitation in an elastic liquid［J] .Physical Review，1943，63(5/6):172－181.

［5] Newton R E.Finite element analysis of shock-induced cavitation，1980. Preprint 80 － 110， ASCE Spring Convention.

［6] Felippa C A，DeRuntz J A.Finite element analysis of shockinduced hull cavitation［J] .Computer Methods in Applied Mechanics and Engineering，1984，44:297－337.

［7] Ronquist E M，Patera A T.A legendre spectral element method for the Stefan problem［J] .International Journal for Numerical Methods in Engineering，1987，24:2273－2299.

［8] Maday Y，Patera A T.Spectral element methods for the incompressible Navier-Stokes equations［J] .The American Society of Mechanical Engineers，New York，1989.

［9] Komatitsch D，Vilotte J P.The spectral element method:An efficient tool to simulate the seismic response of 2D and 3D geological structures［J] .Bulletin of the Seismological Society of America，1998，88(2):368－392.

［10] Mulder W A.Spurious modes in finite-element discretizations of the wave equation may not be all that bad［J] .Applied Numerical Mathematics，1999，30:425 445.

［11] Bleich H H，Sandler I S.Interaction between structures and bilinear fluids［J] .International Journal of Solids and Structures，1970，6:617－639.

［12] Ronquist E M，Patera A T.A Legendre spectral element method for the Stefan problem［J] .International Journal for Numerical Methods in Engineering，1987，24:2273－2299.