水下爆炸气泡作用下船模的动态响应

水下爆炸气泡作用下船模的动态响应

这是一个非常复杂的问题，包括模拟水下爆炸负荷、爆炸气泡的动态性能、固结界面分析、水下结构的非线性动态响应等。目前,对于这类三维结构动态问题的主要研究方法是边界元/有限元结合的方法。主要采用流固耦合方程来表征耦合界面上的载荷,而不需要对水域进行建模,因此减少了大量的计算量。这种方法主要经历了平面波理论,曲面波理论,虚质量理论到双渐近(DAA)理论的发展过程。特别是由Geers所提出的DAA理论,综合并改进了前面几种方法,在高、中、低频段都有较好的精度。DAA方法经过众多学者们的改进与发展,目前是一种解决流固耦合问题的非常有效的计算方法。水下爆炸中对舰船主要产生两种载荷:一是瞬态的冲击波,其压力峰值很高但持续时间很短,通常造成水中结构的局部破坏;二是气泡脉动载荷,气泡中大约含有47%的爆炸能量,具有非常强的破坏力。这一阶段的明显特征是达到的压力峰值较冲击波的压力峰值低,而且持续时间很长,约为几百毫秒,会造成水中结构的总体响应和局部响应。对于气泡的动力学,众多学者已经进行了大量的研究工作[13,14,15,16,17,18,19]。对于气泡载荷作用下的舰船的动态响应,目前的主要研究方法都是将舰船简化为船体梁,研究船体梁的总体鞭状响应[20,21,22,23,24,25,26]。然而,对于气泡作用下三维全船结构的动态响应的细节分析,开展的研究还很少。本文致力于建立一套有限元方法与DAA方法相结合的计算程序,用以研究水下爆炸气泡作用下水面舰船的动态响应。首先,建立一个考虑迁移效应,自由面效应和气泡阻力的气泡模型和一个水面舰船结构模型,并阐述了流固耦合与结构响应相关的理论分析;然后以该船模作为算例,研究了气泡作用下船体模型的总体响应和局部响应。比较了不同位置的加速度,速度和位移时程曲线。最后详细讨论了气泡作用下船体模型动态响应的特征与机理。1理论背景1.1流固耦合方程的建立一个弹性结构受到外部激励的运动方程可以表示为:Μs⋅⋅x+Cs˙x+Κsx=F(1)Msx⋅⋅+Csx˙+Ksx=F(1)其中:Ms,Cs,Ks分别为N×N阶的质量矩阵,阻尼矩阵和刚度矩阵。N×1阶向量x,˙x,⋅⋅xx,x˙,x⋅⋅分别为结构的位移,速度和加速度。列向量F代表外力。N为结构的自由度总数量。对于浸没在流体的结构所受的外部激励可以表示为:F=-GAf(pi+ps)(2)式中:pi和ps分别为湿表面上入射流和辐射流作用下产生的节点压力向量,Af为湿表面单元的面积对角矩阵,G为结构节点力和湿表面节点力的坐标转换矩阵。当一个水中结构如水面舰船或潜艇浸没在无限声学流体介质中,结构湿表面的壳单元的控制方程可以由DAA方法来表示。对于流体中的弹性结构,DAA方法表示的结构表面流体的运动即为各正交流体边界模态的线性组合。一阶双渐近DAA方程的矩阵形式表达为:Μf˙ps+ρcAfps=ρc˙us(3)Mfp˙s+ρcAfps=ρcu˙s(3)式中:Mf为N×N阶的流体质量矩阵;ρ和c分别为流体密度和流体中声速;us是辐射流的流体法向速度向量。根据在流固耦合湿表面上,结构和流体的法向速度相等的条件,得到:GΤ˙x=ui+us(4)GTx˙=ui+us(4)式中:ui为入射流的流体法向速度向量。将式(2)代入式(1),式(4)代入式(3),得到了下面的流固耦合方程组:Μs⋅⋅x+Cs˙x+Κsx=-GAf(pi+ps)(5)Msx⋅⋅+Csx˙+Ksx=−GAf(pi+ps)(5)Μf˙ps+ρcAfps=ρc(GΤ⋅⋅x-˙ui)(6)Mfp˙s+ρcAfps=ρc(GTx⋅⋅−u˙i)(6)为了求解的方便,将式(5)和式(6)改写成对称形式。令ps=GΤ˙y‚ps=GTy˙‚这样式(5)和式(6)可以改写成:Μs⋅⋅x+Cs˙x+Κsx=-GAf(pi+GΤ˙y)(7)ΜfGΤ⋅⋅y+ρcAfGΤ˙y=ρc(GΤ⋅⋅x-˙ui)(8)将式(8)左右两边同乘以GAfM-1f,得到:GAfGΤ⋅⋅y+ρcGAfΜ-1fAfGΤ˙y=ρcGAf(GΤ⋅⋅x-˙ui)(9)这样式(7)和式(9)就可以联立进行数值耦合求解。1.2阶龙格库塔法冲击波之后,大约47%的爆炸能量仍然存留在气泡中。假设流体为无旋,不可压缩的,气泡中心在自由水面下深度d处。由于本文考虑的是中远场爆炸,假设船体对气泡没有影响,并且气泡在运动过程中保持球形。为考虑气泡的迁移效应,自由面效应和气泡阻力的影响,采用Vernon的无量纲方程组:σ=˙x‚(10a)λ=˙ζ‚(10b)˙σ=-3δ(2δ-βx)[σ2x(1-2βx3δ)-λ26x+ζxζ0-(γ-1)kx3γ+1+βx4δ2(Cdλ24x+σλ3-xζ0)](10c)˙λ=-3α[1ζ0+σλx-Cdλ24x+βx4δ2(3σ2+x˙σ)](10d)式中:x是无量纲气泡半径,ζ0为无量纲初始压头。ζ为无量纲压头,τ为无量纲时间;δ为无量纲深度,k是无量纲能量参数(对于TNT炸药,k≈0.0743(z0)1/4,z0为初始压头),γ为绝热气体参数,取为1.25,Cd为阻尼系数,取为2.5。α为气泡迁移控制系数(α只取值0或1;α=0表示不考虑气泡的迁移效应)。在本文中取α=1。β为自由面效应控制系数,同样只取值为0或1;β=0表示不考虑自由面效应。在本文中取β=1(即考虑自由面效应)。接着定义L为长度尺度因子,T为时间尺度因子。L=[3E04πρgz0]13‚Τ=1.428W13z560(11)式中:E0为爆炸的总能量,ρ为流体的密度,g为重力加速度,W为炸药量。可以得到以下无量纲参数:x=aL‚ζ=zL‚ζ0=z0L‚δ=dL‚τ=tΤ(12)其中:z为压力水头。只要给定初始条件,式(10a)～(10d)可以由四阶龙格库塔法进行求解。选定初始条件为x=x0,ζ=ζ0,σ=0,λ=0。而其中x0可以由下面能量守恒方程得到:˙x2x3+kx3(γ+1)+x3-1=0(13)方程(10a)～(10d)的解即为x(t)和ζ(t)。速度势函数可以表示为:φ=e1r+e2r2cosθ-e1r1+e2r21cosθ1(14)式中:r为所取源点到气泡中心的径向距离,r1为该源点对应的偶极到气泡中心的径向距离。θ为r方向与垂直方向的夹角,θ1为r1方向与垂直方向的夹角。e1为源强系数,e2为偶极强度系数,定义如下:e1=L3x2Τσ(15)e2=L42Τ(x3λ-x5σ4δ2)(16)˙e1=L3Τ2(2xσ2+x2˙σ)(17)˙e2=3L4σx22Τ2[λ-x2σ4δ2]+x3L42Τ2[˙λ-(2xσ+x2˙σ4δ2)](18)对于势流,速度为速度势的负梯度,即为:uΖ=-∇φ=1(X2+Y2+Ζ2)[e1Ζ(X2+Y2+Ζ2)12-e2(X2+Y2+Ζ2)12(1-3Ζ2(X2+Y2+Ζ2))](19)式中:X和Y分别为以气泡中心作为原点的直角坐标系的横坐标和纵坐标。流体加速度可以表示为:˙uΖ=∂uΖ∂t-v∂uΖ∂Ζ(20)其中:ν为气泡的垂向平均速度。式(20)中的两项可以分别表达为下面两式:∂uΖ∂t=1(X2+Y2+Ζ2)32[˙e1Ζ-˙e2(1-3Ζ2(X2+Y2+Ζ2))](21)∂uΖ∂Ζ=e1(X2+Y2+Ζ2)32(1-3Ζ2X2+Y2+Ζ2)+3Ζe2(X2+Y2+Ζ2)32(3-5Ζ2X2+Y2+Ζ2)(22)这样,可以计算每一时间步的e1,e2,˙e1,˙e2并把它们代入式(21)和式(22)中,然后由式(20),就可以求解得到流体加速度˙uγ。同样,任意点处的气泡压力p可以由伯努利方程p=-ρ∂φ/∂t-ρgw来计算,其中w为测点到自由面的垂直距离。2数值优化算法本文的程序中,首先利用有限元方法提取结构质量矩阵和刚度矩阵,然后利用一阶DAA方法来求解流体方程。一阶DAA方法实质上是在高频频域段和低频频域段分别采用平面波近似理论和虚质量近似理论进行逼近近似,中频段采用线性的过度。这样就可以在从高频频域到低频频域的区域都具有较高精度。文献通过对浸没在水中的圆形球壳的数值计算结果与精确解的比较,验证了一阶DAA方法在求解一般问题上的精度可以满足要求。而程序中的三维流体质量矩阵利用边界积分方法进行求解。在进行耦合计算时,采用分部交错迭代(Staggeredsolution)的方法,在每一个时间步内,采用预报-代入-判断-同步的方式进行迭代。结构的计算采用Wilson-θ法,考虑到计算的稳定性,流体的计算采用隐式的单步求解方法。计算程序流体图如图1所示。3船模在船底内部的垂线位移响应取一艘船模作为算例来分析其在气泡载荷作用下的动态响应。示意图如图2所示。此船模总长4.5m,型宽0.6m,型深0.45m,吃水T=0.265m,排水量390kg。船模共设7道横舱壁、两层平台,有中纵舱壁,平台2与船底外底板组成双层底。甲板、平台1、平台2的板厚4mm,其余结构的板厚3mm。采用普通钢建造,密度ρ=7850kg/m3,弹性模量E=210GPa,泊松比μ=0.3,屈服极限约为235MPa。船模的有限元模型如图3所示。船模共有2373个节点,2648个单元,其中包括661个湿表面单元。对于非接触水下爆炸,破坏力最大的情况是炸药在舰船正下方发生爆炸。所以本文中,假设药包在船中正下方爆炸,如图3所示,爆炸的位置可以由爆炸深度h来表示。取6kgTNT炸药在自由表面下10m处爆炸。为忽略冲击波效应的影响,取爆炸后0.05s作为初始时间。气泡半径,气泡压力和气泡中心深度的时间历程曲线如图4所示。从图中可以看到,气泡半径先增大到最大值,然后在t=0.26s时减小到最小半径。在这个过程中,随着气泡半径的减小,气气泡压力随着气泡半径的减小而增加。当气泡半径减小到最小值之后,气泡开始回弹,因为此时气泡内部压力已经变得非常大。随后气泡半径随着时间而增大。在回弹过程中,气泡压力迅速衰减到零以下。当气泡半径在最大值附近时,气泡中心位置上升的很缓慢,而当气泡半径要达到最小值时,气泡开始快速地向上迁移。图5给出了一系列不同时刻的船模的垂线位移响应。其中云图给出的是全船的响应分布情况,云图下方的曲线表示船底龙骨的位移响应。从图5中可以清楚地观察到船模在气泡在作用下的位移响应主要为船体的总体响应。船模按梁的一阶垂线振型进行总体鞭状运动。这主要因为气泡载荷的频率较低,与船体的低阶垂向固有频率相接近,因此激起了船体的低阶振型,使船体发生鞭状运动。船中处位移最大,图中可以看到最大位移可达到2.4cm。这样的鞭状振动使船中反复遭受作用力,严重时会对船体造成较大的总体损伤。为进一步研究船模在不同位置处的响应特性,在船模的中线面上取几个典型的位置作为测点。其中包括船底龙骨上的测点(B1,B2,B3)和主甲板上的测点(D1,D2,D3)。测点的具体位置如6所示。图7给出了不同测点的加速度,速度和位移响应的时程曲线。从图中可以观察到所有测点的加速度,速度和位移响应均开始于大约0.25s时,这是气泡载荷开始迅速增长至峰值的时刻。然后响应的幅值开始缓慢平滑的衰减。由于炸药位于船中的正下方,所以船中的测点B2和D2的加速度,速度和位移的峰值要比其他测点大很多。如B2的最大加速度值为1869.6m/s2,分别约为B1(726.8m/s2)和B3(1474.8.8m/s2)的2.57倍和1.27倍。速度和位移响应也遵循同样的规律。在同一纵向位置,主甲板和船底龙骨处的响应相差不大。如B1和D1,B2和D2,B3和D3的加速度,速度和位移的峰值都比较接近。除了船体的总体响应,船模上的一些区域也出现了明显的局部响应,如图8和图9所示。图8给出了在t=0.260s时刻的垂向位移响应。可以观察到在船中处的舷侧外板上有明显的局部大位移响应区域。在此区域取两个测点T1和T2,如图所示,此时T1的垂向位移为3.19cm,T2的垂向位移为2.73cm。两者相差17%。结合图5(a)～(c)中的位移云图,可以看到沿着船长方向,在舷侧外板上还有几个位置均出现了局部大位移响应区域,分布在船中两侧,比船中处小一些。图9给出的是t=0.271s时刻船模的垂向加速度响应。可以观察到在主甲板的舱口处,出现了明显的局部加速度集中区域。同样取测点T3和T4,如图所示,此时T3的垂向加速度为3166.12m/s2,T2的垂向加速度为1520.66m/s2。T3为T4的2.08倍。这是因为船舶甲板舱口角隅,由于形状不连续,在船体发生鞭状运动受到较大的面内载荷时,使局部的应力梯度升高,产生应力集中,严重时可能会造成塑性变形与屈服。因此在结构设计时,为了降低这种应力集中程度,应该采取加厚板、复板或形状优化设计等使舱口角隅的最大加速度极小化。通常认为,水下爆炸主要引起船体的垂直方向的响应。其实在某些位置,船体的横向响应也较为明显。图10给出的是t=0.258s时刻船模的横向位移响应。可以观察到沿着船长方向,在舷侧外板上有多处较大的局部大位移区域,尤其在艏尖舱壁位置尤为明显,在艏尖舱壁处出现了面积非常大的局部大位移响应区域,而且左右两侧的位移是相反方向的。在左右两侧分别取两个对称的测点T5和T6,它们的加速度时间历程曲线如图11所示。可见两个测点一直在做相位相反的振动,即按壳体的呼吸模态振动。而两个测点T5和T6的横向加速度峰值分别为2899.9m/s2和2784.7m/s2。它们已经和上面分析的垂向加速度在同一量级,严重时会造成船体结构的局部破坏,必须引起重视。4船