0.   引　言

内爆（implosion）是中空壳体在外界静水压力作用下发生屈曲，导致流场发生水锤效应，流体动能转化为冲击波能量的过程。内爆机制受到结构承载能力和缓冲吸能两方面的共同影响，材料属性和结构形式是影响内爆变形模式的重要因素。

在过去20年中，许多学者对圆柱壳内爆过程进行了研究^[1-2]，分析了结构在不同变形阶段与周围流体压力的关系。纪冲等^[3]将钢质圆柱壳置于TNT药柱产生的爆炸场中进行冲击实验，并基于数值计算确定了壳壁发生破裂的临界装药距离。Sun等^[4]通过模拟金属圆柱壳受到静水压力和水下冲击波载荷的内爆过程，揭示了初始静水压力和流固耦合对金属圆柱壳体动力稳定性的影响，初始静水压力的存在降低了结构刚度。

在复合材料中空结构水下内爆方面，近年来有不少学者进行了研究。Pinto等^[5]开展了碳纤维增强复合材料（carbon fiber reinforced plastic，CFRP）圆柱壳的内爆实验，相比于单向纤维管，编织纤维管内爆会产生更高的冲击波载荷，其失效特征主要为纵向和环向脆性裂纹。Pinto等^[6]对比了玻璃纤维/聚酯复合材料圆管和碳纤维/环氧复合材料圆管的结构损伤和失效传递的内爆过程，前者内爆产生的冲击波压力约为后者的一半。Moon等^[7]研究发现纤维角度变化会显著影响CFRP管的内爆变形模式。CFRP的抗冲击性能已经获得学术界的广泛认可^[8-10]。目前，学界已开展了许多关于CFRP的抗冲击实验和研究，但对CFRP圆柱壳在静水压力和冲击载荷共同作用下内爆问题的相关研究很少。

本文拟探讨CFRP圆柱壳在静水压力和冲击载荷作用下的失效模式与动态响应特征，通过数值模拟，研究水下爆炸载荷作用下复合材料圆柱壳的内爆问题，以及对圆柱壳的静力屈曲和冲击屈曲、失效和断裂机制、结构动态响应、流固耦合过程等，并考虑结构参数的影响，分析圆柱壳的长度直径比（以下简称“长径比”）、纤维层数和冲击块速度对失效模式的影响。

1.   数值模拟方法验证

1.1   等效水下冲击载荷

真实的水下爆炸试验存在危险性高、对场地要求高、费用高昂以及可重复性低等缺点，因此通常设计水下爆炸模拟装置来进行实验测试。水下爆炸模拟装置实验的基本流程为：弹片撞击活塞，使活塞中产生应力波并向后传播；由于活塞的后方是密封的圆柱形水域，根据应力波理论，应力波传到活塞与水的交界面时，会向水中发射冲击波（此时水中也会向活塞内部反射稀疏波），在水中传播的冲击波向后传播直到作用于靶板。透射波与水下爆炸冲击波衰减的理论公式在形式上是统一的，可以通过改变弹片和活塞的材料、接触横截面积、弹片质量或弹片初速度等参数，控制压力波的压力峰值和衰减时间，从而模拟不同强度的冲击波。本文重点关注水下爆炸冲击波载荷作用下CFRP圆柱壳的动态响应，为了较准确地模拟冲击波载荷的产生和传播过程以及流固耦合作用效果，采用水下爆炸冲击波模拟装置的方式加载。

水下爆炸冲击波压力$p\left(t\right)$随时间呈指数衰减的经验公式可以表示为^[9]

式中：${p}_{0}$，$\theta$分别为峰值压力和衰变时间，两者均取决于爆炸材料的质量、类型以及距离。基于经典 Taylor分析，较远流场处的峰值压力可以通过改变飞片的初始速度进行调节，冲击波的衰减速率可以通过改变飞片和活塞的质量进行调节。一维水下冲击载荷的脉冲强度可以表示为

式中：${c}_{\mathrm{w}}$和${\,\rho }_{\mathrm{w}}$分别为水中的声速和水的密度；${v}_{0}$为冲击块速度；${m}_{\mathrm{f}}$为冲击块的质量。

1.2   计算模型

本文采用非线性有限元软件ABAQUS显式求解器中基于流体体积法（(volume of fluid，VOF）的耦合欧拉−拉格朗日方法（coupled Euler Lagrange，CEL）。CEL有限元方法是欧拉有限元方法与拉格朗日有限元方法相互结合的一种耦合有限元算法，该算法解决了分析物质产生的大变形流动问题，其中，结构部分用拉格朗日单元模拟，流体部分用欧拉单元模拟，计算过程中欧拉单元空间网格保持不变，物质在网格间运动。本文对Farhat等^[11]的实验工况进行模拟，并结合文献[11]的实验参数验证CEL有限元方法对水下爆炸冲击问题处理的可靠性。

本文模拟参照文献[11]中金属圆柱壳水下爆炸实验工况，其内爆模型如图1所示。圆柱壳直径$D=38.2\;\mathrm{m}\mathrm{m}$，长度$L=76.2\;\mathrm{m}\mathrm{m}$，厚度0.711 mm；水域采用与圆柱壳相似的圆柱体，水域直径${D}_{\mathrm{w}}=200\;\mathrm{m}\mathrm{m}$，长度${L}_{\mathrm{w}}=400\;\mathrm{m}\mathrm{m}$。圆柱壳材料设置为双线性弹塑性固体，弹性模量为69.6 GPa，屈服应力为292 MPa，屈服后模量为674 MPa，泊松比为0.3，密度为2779 kg/ m³。金属圆管采用壳单元，计算流体域采用欧拉单元。金属圆管表面的单元尺寸为1 mm，水域网格和圆柱壳网格共结点，其尺寸向外逐渐变大，水域外部边界单元尺寸为10 mm。模型有79494个单元。计算初始，只有圆管外有水，圆管内为空气。

图  1  整体计算模型

Figure  1.  Volume type global computing model

下载: 全尺寸图片 幻灯片

本文采用U_s− U_p状态方程来模拟不可压缩黏性流体的层状流动，并通过N-S方程模拟水在冲击波载荷作用下的运动过程。水介质状态方程的具体形式为^[12]

式中：P_n为水介质内的压力；${P}_{\mathrm{H}}$为Hügoniot压力；${\varGamma }_{0}$为材料常数；$\eta$为名义体积压缩应变，其表达式为$\eta =1-\rho/{\rho }_{0}$，其中，ρ为压缩状态下的水介质密度，${\rho }_{0}$为参考密度； ${E}_{\mathrm{m}}$为单位质量内能。表1为水介质材料参数，表中 C₀为 U_s− U_p曲线的截距，S为 U_s− U_p曲线斜率的系数， μ为动力黏度。

表  1  水介质材料参数

Table  1.  Water material parameters

参数	数值
密度$/({\rm{g} }\cdot {\rm{c}\rm{m} }^{-3})$	1
$C_{0}/(\rm{m}\rm{m}\cdot {\rm{s}}^{-1})$	1 480 000
$S$	0
${{\varGamma }}_{0}$	0
$\mu /(\rm{P}\rm{a}\cdot \rm{s})$	1.0×10−3

下载: 导出CSV 

| 显示表格

本文采用理想气体状态方程描述空气，其状态方程为

式中：p为压强；n为物质的量；R为比例系数；T 为温度；V为体积。理想气体状态方程具体参数如表2所示。

表  2  空气介质材料参数

Table  2.  Air material parameters

参数	数值
密度$/(\mathrm{g}\cdot {\mathrm{c}\mathrm{m} }^{-3})$	$1.293\times {10}^{-3}$
气体常数	287 000 000
环境压力/MPa	0.101 3
$\mu /(\mathrm{P}\mathrm{a}\cdot \mathrm{s})$	$8.25\times {10}^{-5}$

下载: 导出CSV 

| 显示表格

在圆柱壳跨中处均匀分布了如图2(a)所示的6个测点，以便输出压力时程曲线。这些测点的选取位置与文献[11]中6个中心传感器晶体的位置一致，如图2(b)所示。流场的初始压力与文献[11]保持一致，即流场内使圆柱壳发生内爆的临界静水压力${P_{{\text{CO}}}} = 4.46\;{\text{MPa}}$。在欧拉域外边界加载压力，压力与预定义场压力相等。

图  2  测点及压力传感器位置^[11]

Figure  2.  Positions of measuring point and pressure sensor ^[11]

下载: 全尺寸图片 幻灯片

1.3   变形模式和测点压力验证

图3显示了金属圆管实验的高速摄影和数值模拟结果，结合左右两边金属圆管的变形和损伤情况，可以发现金属圆管的模拟变形模式与实验结果吻合良好。第4模态坍塌清晰可见，随着时间的推移，金属圆管的变形越来越严重。在静水载荷的作用下，圆管的损伤区域由跨中向两端扩散。

图  3  金属圆管实验^[11]与数值模拟结果对比

Figure  3.  Comparison of experiment^[11] and numerical simulation results of metal tube

下载: 全尺寸图片 幻灯片

图4所示为实验中传感器测得的压力时程曲线与数值模拟得到的压力时程曲线的对比结果，为了便于比较实验和数值模拟结果，移动时间轴，使二者从同一时刻开始加载。图中${P}_{\mathrm{d}\mathrm{y}\mathrm{n}}$为测点处的动压载荷，即传感器测得的实际压力减去水中初始静水压力的相对压力。由于静水压力达到圆柱壳的临界载荷，圆柱壳发生屈曲，圆柱壳周围流场压力下降，流场中水粒子获得了流向圆柱壳壁面的加速度。由于流场中粒子运动的惯性作用，导致圆柱壳发生壁面接触时，流场中粒子与圆柱壳壁面碰撞，产生脉冲载荷（图4中B点处），并以压缩波的形式向外传播。

图  4  实验^[11]和数值模拟的压力信号对比

Figure  4.  Comparison of pressure signals between experiment ^[11] and numerical simulation

下载: 全尺寸图片 幻灯片

2.   数值模拟及结果分析

2.1   计算模型与网格收敛性

计算模型如图1所示，前半部分是为了模拟水下爆炸冲击波，后半部分的水域里放置了CFRP薄壁圆柱壳。通过罚函数近似方法完成拉格朗日界面和欧拉界面的耦合计算，而不需要设置公用的单元节点或单元面^[13]。在预定义场内确定流场的初始压力，使圆柱壳发生内爆的临界静水压力${P}_{\mathrm{C}\mathrm{O}}=1.5\;\mathrm{M}\mathrm{P}\mathrm{a}$。然后，在欧拉网格的外边界加载压力，其大小与预定义场压力大小相等，方向相反。冲击块和活塞的材料为钢，密度为7890 $\mathrm{k}\mathrm{g}/{\mathrm{m}}^{3}$，杨氏模量为206 GPa，泊松比为0.3，质量分别为0.79和0.95 kg。冲击块的初始速度为10 m/s。考虑冲击块和活塞的质量、冲击块速度以及加载面积对冲击波压力峰值与衰减时间常数的影响，式(2)和式(3)可进一步表示为^[14]：

式中：$k$为影响系数；${D}_{1}$为冲击波加载区域直径；${m}_{\mathrm{p}}$为活塞质量；a，b，c，d，e为控制系数。

以长径比L/D=5的CFRP圆柱壳模型为例，计算网格收敛性，结果如图5所示。图中，α 为网格尺寸与CFRP圆柱壳单层碳纤维布厚度之比。由图可见，网格尺寸越小，CFRP圆柱壳的临界静水压力P_CO逐渐收敛。本文CFRP圆柱壳单层碳纤维布厚度为0.25 mm，当$\alpha$≤8，即网格尺寸不大于2 mm时，网格尺寸对临界静水压力的影响微弱。因此，将α设置为8，水域（欧拉域）共有31148个网格，L/D=5的CFRP圆柱壳共有60531个网格。

图  5  长径比为5的CFRP圆柱壳模型网格收敛性

Figure  5.  Mesh convergence of CFRP cylindrical shell model with a length-to-diameter ratio of 5

下载: 全尺寸图片 幻灯片

计算模型的水介质和空气介质材料参数与上文验证实验模型中保持一致。CFRP板由多层碳纤维预浸料和黏合材料组合而成，每层碳纤维布采用壳单元模拟，两层板之间采用黏性材料模拟，计及Quads损伤及其损伤演化。内聚力单元的损伤起始由名义应力相关的二次函数决定：

式中：t_n为沿x轴的轴力；t_s和t_t分别为沿y轴和z轴的剪力；${t}_{\mathrm{n}}^{0}，{t}_{\mathrm{s}}^{0}，{t}_{\mathrm{t}}^{0}$分别为沿3个方向的最大名义应力；“< >”为麦考利算子。当式（8）的值达到1，内聚力单元开始损伤，然后基于Benzeggagh–Kenane (BK)准则损伤演化：

式中：${G}^{\mathrm{C}}$ 和 $G$ 分别为临界断裂能和外力功；${G}_{\mathrm{n}}，{G}_{\mathrm{s}}\mathrm{和}{G}_{\mathrm{t}}$为沿3个方向的外力功。CFRP板的内聚力单元力学参数引自文献[15]。对于CFRP编织材料，纤维方向平面内应力−应变关系假定是正交线弹性。CFRP材料参数和损伤演化参数也引自文献[15]，具体材料参数如表3所示，单层碳纤维编织布厚度为 0.25 mm，树脂含量为 30%，面密度为$270\;\mathrm{k}\mathrm{g}/{\mathrm{m}}^{2}$。

表  3  CFRP板的本构模型参数^[15]

Table  3.  Constitutive model parameters of CFRP plate^[15]

参数	数值
密度$/(\mathrm{g}\cdot {\mathrm{c}\mathrm{m} }^{-3})$	4
杨氏模量E1$/$GPa	42.7
杨氏模量E2$/$GPa	42.7
剪切模量$/$GPa	4.4
泊松比${\nu }_{1}$	0.05
泊松比${\nu }_{2}$	0.05

下载: 导出CSV 

| 显示表格

2.2   模拟结果及分析

图6为圆柱壳跨中测点的动压时程曲线。根据 CFRP圆柱壳的失效过程，可以将其水下内爆过程分为3个阶段：屈曲阶段、壁面接触阶段、失效扩展阶段。在0.5 ms附近出现一个较大的压力峰值，该峰值由冲击波传播产生。在3.4 ms时，压力从A点开始平稳下降，这是由于周围流体沿着凹陷的圆柱壳壁面加速流动，流体压力势能转换为流体动能，导致局部动压下降。此时，圆柱壳进入屈曲阶段，但未发生明显损伤，初始屈曲是从跨中处开始，方向由冲击波传递过来的一侧逐渐向对侧凹陷。在4.1 ms左右，压力突然下降至B点，圆柱壳继续向内凹陷，变形区域从跨中向两端扩展开来。在达到最小值后，压力进入短暂的平稳阶段，维持大概0.5 ms。之后，圆柱壳发生壁面接触，壁面接触突然阻止水流运动，流体动能快速衰减，产生冲击波，所以在C点出现了一个尖锐的压力峰值，此时的圆柱壳结构发生大规模破坏。随着复合材料的屈曲和失效效应逐渐扩展，水流发生多次撞击，产生多个正负交替的峰值。图7为CFRP圆柱壳的位移云图，右边是跨中至端部一半圆柱壳的变形图，对应图6中的C点。

图  6  CFRP圆柱壳跨中测点的动压（P_CO=1.5 MPa）

Figure  6.  Dynamic pressure of the measuring points in the span of CFRP cylindrical shell（P_CO=1.5 MPa）

下载: 全尺寸图片 幻灯片

图  7  CFRP圆柱壳位移云图

Figure  7.  Deformation contours of CFRP cylindrical shell

下载: 全尺寸图片 幻灯片

3.   参数对失效模式的影响

3.1   长径比

在静水加载条件下，将水、钢、CFRP等各材料参数设置一致，圆柱壳直径不变，仅改变圆柱壳的长度，在L/D 分别为5和8的情况下，进行数值模拟。两圆柱壳测点的净水内爆压力−时间曲线如图8所示。当测点处压力P开始下降，圆柱壳开始发生屈曲。线性加载结果表明，L/D=5的圆柱壳在P=3.6 MPa开始屈曲，而L/D=8的圆柱壳在P=2.9 MPa开始屈曲。这表明L/D越大，圆柱壳发生屈曲的临界静水压力越小。其原因是，当圆柱壳直径保持一定，随着圆柱壳长度的降低，结构刚度逐渐增大，因此承载能力增大。随后，圆柱壳发生壁面接触，壁面接触停止后，曲线出现压力峰值。对比可见，L/D=5的圆柱壳产生的压力脉冲振幅比L/D=8的高，但持续时间短，部分原因是L/D小的圆柱壳刚度较大，坍塌压力较高。

图  8  静水内爆压力−时间曲线对比

Figure  8.  Comparison of hydrostatic implosion pressure versus time curves

下载: 全尺寸图片 幻灯片

在静水载荷和冲击波载荷共同作用下，圆柱壳的应力云图如图9所示，图形显示为沿中剖面取一半长度的圆柱壳轴测图。L/D=5的圆柱壳横剖面呈工字形，L/D=8的圆柱壳横剖面呈T字形。这是由于L/D越大，圆柱壳的轴向越容易发生失稳，根据经典的弹性塑性屈曲理论和最小势能原理，金属圆柱壳的L/D会影响其屈曲模式，CFRP圆柱壳的L/D也会影响其失效模式。CFRP材料不同于具有良好延展性的金属材料，无法抵抗过大的塑性变形，因此CFRP材料容易发生不规则的失效变形，其规律在于圆柱壳L/D越大，圆柱壳的壁面接触程度越大。

图  9  圆柱壳应力云图

Figure  9.  Stress contours of cylindrical shell

下载: 全尺寸图片 幻灯片

将2组圆柱壳测点的压力减去各自的临界静水压力P_CO得到相应的动压载荷P_dyn，2种工况下的动压变化曲线如图10所示。首先，由于受到冲击波的压力脉冲，2条曲线的A点均出现一个较大的压力峰值，且峰值接近。然后，圆柱壳开始屈曲，流体沿壁面流动，局部压力开始下降，L/D=5的圆柱壳附近测点的最小负压（C点）明显小于L/D=8的圆柱壳对应测点的最小负压（B点）。接着，圆柱壳发生壁面接触，水无法继续向内流动，此时出现压力峰值，且L/D=5的圆柱壳的压力峰值（E点）要大于L/D=8的圆柱壳的压力峰值（D点）。

图  10  长径比为5和8的圆柱壳动压变化

Figure  10.  Dynamic pressure of cylindrical shells with length-to-diameter ratios of 5 and 8

下载: 全尺寸图片 幻灯片

3.2   纤维层数

保持圆柱壳的直径与长度均不变，改变纤维层数。单层CFRP的厚度是0.25 mm，所以纤维每减少1层，圆柱壳厚度相应减小0.25 mm。纤维层数减少会降低圆柱壳的临界静水压力P_CO，且纤维层数对临界静水压力为非线性影响，即随着纤维层数的增大，圆柱壳结构的静水承载能力增长速率增加。参数设置如表4所示。

表  4  圆柱壳纤维层数的参数设置

Table  4.  Parameter setting with different fiber layers of cylindrical shell

纤维层数	长度L/mm	直径D/mm	厚度/mm	临界静水压力PCO/MPa	冲击块速度v0/(m·s−1)
4	292	36.5	1.00	1.5	10
3	292	36.5	0.75	1.2	10
2	292	36.5	0.50	1.0	10

下载: 导出CSV 

| 显示表格

数值模拟得到的测点动压载荷时程曲线如图11（a）所示，随着圆柱壳纤维层数的减少，压力开始下降的时间提前，同时内爆压力峰值P_im相应地增大，内爆压力峰值P_im与内爆时间随层数的变化更加直观地反映在图11（b）中。这说明圆柱壳结构纤维层数越少，越不稳定，越容易发生内爆，且内爆产生的压力峰值越大。同样，纤维层数对圆柱壳内爆产生的峰值压力为非线性影响。即随着纤维层数的增大，圆柱壳结构的抗冲击能力增长速率增加。

图  11  纤维层数对内爆的影响

Figure  11.  Effect of fiber layers on implosion

下载: 全尺寸图片 幻灯片

3.3   冲击块速度

采用L/D=8的圆柱壳为模拟对象，设置冲击块速度为10，20和30 m/s，探究冲击块速度对圆柱壳失效模式的影响，模拟得到的变形结果如图12所示。图中，左侧为完整的圆柱壳位移云图，右侧为沿中剖面取一半长度的圆柱壳轴测图。

图  12  冲击块速度对圆柱壳变形模式的影响

Figure  12.  Influence of impactor velocity on deformation mode of cylindrical shell

下载: 全尺寸图片 幻灯片

对比不同冲击块速度对应的圆柱壳变形模式，可以发现增大冲击块速度会造成圆柱壳产生更加严重的损伤。在折叠模式方面，第1组和第3组为T形折叠，第2组为C形折叠。在冲击波速度变化的情况下，圆柱壳的失效模式具有随机性，冲击波速度越大，圆柱壳的壁面接触和失效扩展越显著，发生的基体断裂更多，且裂纹在圆柱壳长度方向上有明显增大趋势。这是由于冲击块速度较小时，圆柱壳的整体弯曲小，载荷主要由圆柱壳一侧承担，容易产生周向裂纹；冲击块速度较大时，圆柱壳受到的弯矩增大，整体变形程度越大。无法承担过大的应变，导致圆柱壳的裂纹沿基体方向迅速扩展。

图13为各组的动压模拟结果。对比发现，增大冲击块速度能够提高冲击波载荷能量，压力曲线的第1个峰值会升高。而且，随着冲击块速度的增大，圆柱壳会提前进入屈曲阶段，即压力下降的时间会提前，可见圆柱壳将更快发生损伤，产生更高的内爆冲击波载荷。

图  13  冲击块速度对动压的影响

Figure  13.  Influence of impactor velocity on dynamic pressure

下载: 全尺寸图片 幻灯片

4.   结　论

本文以水下爆炸冲击波载荷下深水设备的变形毁伤和防护为研究背景，利用数值模拟方法对水下爆炸载荷作用下CFRP圆柱壳的内爆问题展开了研究，主要包括圆柱壳的静力屈曲和冲击屈曲、失效和断裂机制、结构动态响应、流固耦合过程，并讨论了圆柱壳的长径比、纤维层数和冲击块速度对其失效模式的影响。

1） 采用CEL法，利用ABAQUS软件构建水下爆炸载荷下金属结构内爆实验的数值模型；根据金属圆管的变形模式和测点压力，与文献中的实验仿真结果进行对比，验证了数值方法的有效性。

2） 建立水下爆炸载荷下CFRP圆柱壳内爆实验的数值模型，通过网格收敛性分析得到了合适的网格尺寸。计算结果显示，CFRP圆柱壳水下内爆可以分为3个阶段：屈曲阶段、壁面接触阶段、失效扩展阶段。利用水域中四等分位置的3个压力测点得到冲击波时程曲线，观察到水下冲击波压力呈指数型衰减，与实际情况相符。

3） 改变圆柱壳的长径比、纤维层数、冲击块速度等结构参数，分析水下爆炸载荷作用下CFRP圆柱壳结构的动态响应及变形模式。结果发现，减小CFRP圆柱壳长径比能提高结构的抗冲击能力，且影响其失效模式。随着纤维层数的增加，圆柱壳结构的静水承载能力和抗冲击能力增长速率增加；增加冲击块速度，圆柱壳的壁面接触和失效扩展越显著，发生的基体断裂更多，且裂纹在圆柱壳长度方向上有明显增大趋势。