Dynamic damage analysis of carbon fiber reinforced polymer composite pressure vessels

Carbon Fiber Reinforced Polymers (CFRPs) are widely utilized in various industries due to their exceptional strength-to-weight and stiffness-to-weight ratios, making them ideal for a range of demanding applications [1]. One of the critical uses of CFRPs is the construction of confinement vessels, which play a vital role in preventing catastrophic events and containing explosion debris [2]. Consequently, understanding the behavior of CFRP pressure vessels under dynamic loading conditions is essential for ensuring their safety and reliability. However, experimental testing of such parts under representative deformation rates is often expensive and resource-intensive, leading to a growing trend toward reducing destructive testing through the development of high-fidelity Finite Element Analysis (FEA) models [3]. These models offer numerous advantages, including cost-effectiveness by minimizing material waste and energy consumption, improved safety, and accelerated research and development processes. The principal requirement for the use of such models for design is the confirmation of their validity through rational experimental-computational correlation (RECC) [4].
The behavior of CFRP pressure vessels is characterized by significant complexity, influenced by such factors as the number of layers, fiber orientation, and loading rate [5-7]. Experimental investigations, including the use of the split-Hopkinson pressure bar technique, have demonstrated that strain rate plays a critical role in determining the post-impact residual strength and fracture resistance of these materials [8]. Specifically, increased strain rates could enhance the elastic modulus and tensile strength of CFRP composite materials, due to the influence of strain rate on the deformation and damage characteristics of composites [9], [10]. Such high-strain-rate loading conditions are commonly observed in impact scenarios,
particularly within the automotive and aerospace industries, where pressure vessels are often subjected to extreme dynamic loads. To accurately predict the material behavior under these conditions, constitutive material models have been developed [7], [11]. Among these, the Johnson-Cook [12] and Zerilli-Armstrong [13] models have been extensively employed to analyze the dynamic response of CFRP composite tubes under axial impact loading.
Impact-induced damage in composite materials can be categorized into in-plane (intralaminar) damage and interlaminar damage, the latter often referred to as spallation. Spallation occurs when material layers delaminate due to high-impact loading. While previous studies have attempted to model impact damage in composite materials [14], inconsistencies remain in accurately representing intralaminar damage, fiber failure, and matrix cracking. Recent applications of advanced composite-reinforced structures in structural analysis, such as in [15-18] highlight the importance of accurate simulation methods with integrated micromechanical based modeling for predicting damage, free vibration and thermomechanical behavior under severe dynamic loads. Additionally, Fedorenko et al. [19] investigated the behavior of filament-wound T700/LY113 carbon/epoxy composites under shock loading using the electric wire explosion method used in [20], revealing the significant influence of strain rate on the post-impact strength. This study [19] presented valuable experimental data, but did not provide a comprehensive analysis of the intralaminar or interlaminar damage mechanisms.
The present study bridges this gap and builds upon previous research by investigating the effects of high strain rates induced by internal blast loading on CFRP pressure vessels, with a particular focus on the circumferential strain evolution and the concurrent development of intralaminar and interlaminar damage exploiting both Hashin and VCCT criteria, respectively. For modeling CFRP structures under blast and dynamic loading, it is commonly accepted to combine the use of the Virtual Crack Closure Technique (VCCT) for interlaminar delamination and the Hashin criteria for intralaminar damage. Nevertheless, each method has certain presumptions and restrictions. Although it is not necessarily rate-dependent and may need to be extended for high strain-rate applications, Hashin damage model is effective for predicting intralaminar damage because it is physically based and differentiates between fiber and matrix failure modes. While the Tsai-Wu criterion provides a more general but less physically descriptive interaction-based failure prediction, other criteria, like the Puck model, are frequently better suited for matrix-dominated failures or situations involving strong fiber-matrix interfaces. Cohesive zone models (CZM) provide more flexibility for complex crack paths by simulating both crack initiation and growth, whereas VCCT is reliable in capturing crack propagation for delamination modeling. However, fracture toughness, which is known to be strain-rate sensitive, must be accurately characterized for both VCCT and CZM. Static formulations are unable to adequately capture additional inertial effects that may impact delamination growth under blast and highvelocity impact loading. Consequently, the primary failure mechanisms must be taken into account when choosing a model: Puck may perform better in matrix-dominated situations, Hashin is usually effective for fiber/matrix failures, and CZM offers benefits over VCCT when nonlinear crack evolution is substantial [21]. Although its application at high strain rates requires careful calibration and awareness of inherent limitations, the combined Hashin-VCCT framework is chosen because it has been extensively validated for CFRPs under dynamic conditions.
The work of Fedorenko et al. [19] is utilized as the source of information regarding the material properties, experimental configuration, modeling approach, and key relationships. Major additional developments are introduced to enhance the accuracy of numerical simulations. The main purpose of the present paper is to validate complex models based on the RECC approach. The novelty of this work lies in the successful concurrent utilization of the two damage criteria: Hashin's failure criterion for intralaminar damage, and the Virtual Crack Closure Technique (VCCT) for modeling interlaminar delamination in a CEL framework. The results demonstrate improved agreement with experimental data compared to previous studies, particularly in modeling damage patterns observed in explosion shock pulse experiments. This research contributes to the ongoing effort to formulate and validate accurate predictive models for CFRP structures subjected to extreme dynamic loads.

The numerical model used in this study was developed based on the experimental configuration described in [19]. All modeling and simulations were conducted using Abaqus/Explicit FEA software (version 2021) on a highperformance computer [22]. The pressure vessel (p.v.) exploding-wire sudden internal pressurization model setup consists of a carbon fiber reinforced polymer (CFRP) composite p.v. shell, a PMMA solid tube insert, and a Eulerian domain representing the vapor generated from the explosion of the copper wire. The sketch of the model with characteristic dimensions as well as the numerical assembly is illustrated in Fig. 1. The copper wire positioned inside the tube insert is represented by the Eulerian domain that represents the vapor resulting from the wire explosion. In Eulerian analysis, nodes
remain fixed, and material flows through non-deforming elements, whereas in Lagrangian analysis, nodes move with the material. In the real experiment, the copper wire is subjected to a sudden discharge voltage to stimulate an internal explosion. Together with the PMMA cylindrical insert this assembly ensures the passage of the pressure pulse to the CFRP shell, and replicates the important real-life scenario of pressure vessel suffering abrupt internal pressurization due to the explosion of contained reactive chemicals or due to the externally triggered incendiary event. CFRP structures investigated in this study featured two different winding arrangements, namely, with uniaxial fiber reinforcement in the circumferential orientation ( ${\left[0^{\circ }\right]}_{10}$${\left[0^{\circ }\right]}_{10}$[0^(@)]_(10)\left[0^{\circ}\right]_{10}${\left[0^{\circ }\right]}_{10}$ ) and an angle-ply orientation with two directions, ${[+{45}^{\circ }/-{45}^{\circ }]}_5$${\left[+{45}^{\circ }/-{45}^{\circ }\right]}_5$[+45^(@)//-45^(@)]_(5)\left[+45^{\circ} /-45^{\circ}\right]_{5}${[+{45}^{\circ }/-{45}^{\circ }]}_5$. Both samples of ${\left[0^{\circ }\right]}_{10}$${\left[0^{\circ }\right]}_{10}$[0^(@)]_(10)\left[0^{\circ}\right]_{10}${\left[0^{\circ }\right]}_{10}$ and ${[+{45}^{\circ }/-{45}^{\circ }]}_5$${\left[+{45}^{\circ }/-{45}^{\circ }\right]}_5$[+45^(@)//-45^(@)]_(5)\left[+45^{\circ} /-45^{\circ}\right]_{5}${[+{45}^{\circ }/-{45}^{\circ }]}_5$ orientations had the overall thickness of 1 mm . The thickness of each ply was 0.1 mm .

Since a single homogeneous set cannot accurately represent the structure of the composite shell due to the variation in mechanical properties between individual layers, and because we are concerned about failure that happens by delamination due to weakened interfaces between plies; each layer must be modeled with its own set of elements [23]. In this study, the composite shell layup was modeled using 10 separate rings, each with its respective orientation prescribed, as illustrated in Fig. 2. The initial mesh (coarse) for each composite ring consists of 4000 reduced continuum shell elements in the Lagrangian formulation (designated as SC8R in Abaqus, Simulia Inc.), with the sectional circumference discretized into 80 elements and the element dimension along the $z$$z$zz$z$-axis (longitudinal axis of the cylinder) is 0.6 mm . By giving each ply a local material coordinate system, the finite element model explicitly took into account the various fiber orientations, guaranteeing that the orthotropic properties matched the actual fiber directions. Two layup configurations were examined after modeling the composite shell as a stack of ten different plies: a unidirectional ${\left[0^{\circ }\right]}_{10}$${\left[0^{\circ }\right]}_{10}$[0^(@)]_(10)\left[0^{\circ}\right]_{10}${\left[0^{\circ }\right]}_{10}$ stack, where all plies were oriented along the vessel hoop (circumferential) direction; and an angle-ply ${[+{45}^{\circ }/-{45}^{\circ }]}_5$${\left[+{45}^{\circ }/-{45}^{\circ }\right]}_5$[+45^(@)//-45^(@)]_(5)\left[+45^{\circ} /-45^{\circ}\right]_{5}${[+{45}^{\circ }/-{45}^{\circ }]}_5$ stack, where the plies alternated between $+{45}^{\circ }$$+{45}^{\circ }$+45^(@)+45^{\circ}$+{45}^{\circ }$ and $-{45}^{\circ }$$-{45}^{\circ }$-45^(@)-45^{\circ}$-{45}^{\circ }$ with respect to the hoop direction.
The "orientation" keyword in Abaqus/Explicit was used to assign a local material orientation to each ply. The ply fiber direction and the fiber axis ( 1 -direction) were lined up, and the rotation plane was established by defining a reference axis in the global cylindrical coordinate system, usually the vessel axial direction. The fiber axis was then rotated to its correct orientation using the designated rotation angle $(0^{\circ },+{45}^{\circ }$$\left(0^{\circ },+{45}^{\circ }\right.$(0^(@),+45^(@):}\left(0^{\circ},+45^{\circ}\right.$(0^{\circ },+{45}^{\circ }$, or $-{45}^{\circ })$$\left.-{45}^{\circ }\right)${:-45^(@))\left.-45^{\circ}\right)$-{45}^{\circ })$. This process made sure that the Hashin damage criteria, which rely on stresses in the local coordinate system, and the stiffness matrix (with maximum stiffness ${\mathrm{E}}_1$${\mathrm{E}}_1$E_(1)\mathrm{E}_{1}${\mathrm{E}}_1$ along the fiber axis) were accurately assessed for every ply. This method is shown in Fig. 3, where the local directions are indicated by arrows: red for the radial direction ( 3 -axis), yellow for the in-plane transverse direction ( 2 -axis), and blue for the fiber direction ( 1 -axis). The model represents the anisotropic response of each ply by explicitly assigning material orientations, avoiding homogenization assumptions. Because of the accuracy in defining fiber orientation, this method is crucial for precisely forecasting both intralaminar failure and interlaminar delamination.

The initial mesh (coarse) of the PMMA section was composed of 80000 reduced solid elements (designated as C3D8R in the notation of Abaqus, Simulia Inc.) as shown in Fig. 4(a). In "Region 1", comprising the surrounding of the central hole, the mesh size evolves, from in out, from 0.2 mm to 0.3 mm , while in "Region 2", it evolves from 0.45 mm to 0.7 mm along the same direction. The PMMA insert experiences large deformations, wave propagation, and material failure (modeled with JH-2). The C3D8R element is robust and efficient for such large-strain, explicit dynamics simulations involving complex constitutive models and element deletion. The inner channel section representing the fixed space which will be occupied by the copper vapor during its expansion is referred to as the Eulerian part. Beyond its volume, the vapor pressure is deemed insignificant for further examination. As shown in Fig. 4(b), the Eulerian part is modeled using 1440000 cuboid elements, each with a side length of 0.1 mm , in the Eulerian formulation (denoted as EC3D8R in Abaqus notation). The EC3D8R element is the standard first-order brick element for Eulerian analyses in Abaqus. Moreover, the Eulerian formulation is necessary to handle the extreme deformation and flow of the vaporized wire products without the mesh distortion problems that would cripple a Lagrangian mesh.

The reason of representing the p.v. with continuum shell elements (SC8R) is that they provide a computationally effective way to simulate three-dimensional systems on a local level, such as a multilayered laminated composite structure. They combine the three-dimensional geometry of solid elements with the planar nature of conventional shells, making them suitable for global analyses rather than local failure predictions. These elements allow mechanical and/or thermal loading for both static and dynamic procedures, with only displacement degrees of freedom. They also use linear interpolation and are suitable for nonlinear geometric analysis due to their ability to accommodate large rotations, finite membrane deformation, and the consequences of thickness change and transverse shear deformation. Continuum shell elements utilize first-order layer-wise composite theory, allowing for a refined radial response and more accurate transverse shear stress and force prediction by stacking them, unlike conventional shells [22].
A mesh study analysis is essential in FEA simulations, so failure stresses are compared based on changing the mesh size. So, medium and fine mesh configurations are created for the 2 structures, the ${[+{45}^{\circ }/-{45}^{\circ }]}_5$${\left[+{45}^{\circ }/-{45}^{\circ }\right]}_5$[+45^(@)//-45^(@)]_(5)\left[+45^{\circ} /-45^{\circ}\right]_{5}${[+{45}^{\circ }/-{45}^{\circ }]}_5$, and the ${\left[0^{\circ }\right]}_{10}$${\left[0^{\circ }\right]}_{10}$[0^(@)]_(10)\left[0^{\circ}\right]_{10}${\left[0^{\circ }\right]}_{10}$ orientated structures. The fine mesh of each composite ring consists of 49600 continuum shell elements; the circumference consists of 124 elements; and the along-z-axis element dimension is 0.075 mm . The along-z-axis element dimension and number of elements along the circumference for both the PMMA and composite rings are kept equal to ensure better bond separation during spall analysis. Regarding the PMMA insert, , the element length of the refined PMMA mesh was set to 0.08 mm in the radial direction in zone 1 to ensure accurate representation. Zone 2 grew in element size beyond this zone, from 0.22 mm to 0.4 mm , for elements at the outer surface of the cylinder. Furthermore, a stacking direction is introduced to the model to be from bottom to top, which is logical regarding the internal blast loading inside the vessel. The Eulerian mesh is kept the same because it was made as fine as possible from the beginning of the simulation to ensure having accurate results.
The intermediate mesh of each composite ring consists of 20000 continuum shell elements; the circumference consists of 100 elements; and the along-z-axis element dimension is 0.150 mm . In zone 1, the element length of the refined PMMA mesh was set to 0.150 mm in the radial direction to ensure accurate representation. Zone 2 grew in element size beyond this zone, from 0.23 mm to 0.44 mm , for elements at the outer surface of the cylinder. The main mesh parameters of both the p.v. and the PMMA insert are summarized in Tab. 1.

The PMMA material properties, the CFRP elastic properties and strength under quasi-static loading conditions for the simulation were taken from [19] and are summarized in Tabs. 2, 3 and 4, respectively. The elastic engineering constants in Tab. 3 represent the longitudinal elasticity modulus ( ${\mathrm{E}}_1$${\mathrm{E}}_1$E_(1)\mathrm{E}_{1}${\mathrm{E}}_1$ ), the transverse elasticity moduli ( ${\mathrm{E}}_2$${\mathrm{E}}_2$E_(2)\mathrm{E}_{2}${\mathrm{E}}_2$ and ${\mathrm{E}}_3$${\mathrm{E}}_3$E_(3)\mathrm{E}_{3}${\mathrm{E}}_3$ ), Poisson's ratios ( $\nu$$\nu$nu\nu$\nu$ ) and shear moduli ( G ) in all the respective local directions. The in-plane shear modulus ${\mathrm{G}}_{12}(3200\mathrm{MPa})$${\mathrm{G}}_{12}(3200\mathrm{MPa})$G_(12)(3200MPa)\mathrm{G}_{12}(3200 \mathrm{MPa})${\mathrm{G}}_{12}(3200\mathrm{MPa})$ was given the same value as the transverse shear moduli ${\mathrm{G}}_{13}$${\mathrm{G}}_{13}$G_(13)\mathrm{G}_{13}${\mathrm{G}}_{13}$ and ${\mathrm{G}}_{23}$${\mathrm{G}}_{23}$G_(23)\mathrm{G}_{23}${\mathrm{G}}_{23}$. When precise radial shear test data is not available, macro-scale composite modeling frequently uses this simplification, which makes the assumption that the material is transversely isotropic with regard to shear behavior [24]. It is justified for this study because the model is less sensitive to the exact value of ${\mathrm{G}}_{23}$${\mathrm{G}}_{23}$G_(23)\mathrm{G}_{23}${\mathrm{G}}_{23}$ because the main failure modes of interest: intralaminar damage and delamination are primarily driven by in-plane stresses and VCCT Mode-I fracture.

Table 2: PMMA parameters used for the model [19].

For CFRP, the intralaminar fracture energies that represent Hashin were initially calibrated based on [25]. This calibration was followed by extensive trial and error simulations to match the experimental data and damage profiles in [19], and the used values in the FE simulations are given in Tab. 5.

The response of the PMMA insert was modeled using the Johnson-Holmquist (JH-2) constitutive formulation [19], which accounts for high strain-rate effects, damage evolution, and brittle fracture. This implementation captures the volumetric, deviatoric, and damage-dependent behavior of PMMA subjected to shock and impact loading. The volumetric response was represented using a polynomial EOS , consistent with the $\mathrm{JH}-2$$\mathrm{JH}-2$JH-2\mathrm{JH}-2$\mathrm{JH}-2$ model where the pressure $(P)$$(P)$(P)(P)$(P)$ is defined as:

for intact material (damage parameter $(D)=0$$(D)=0$(D)=0(D)=0$(D)=0$ ), and:

for damaged states $(0\le D\le 1)$$(0\le D\le 1)$(0 <= D <= 1)(0 \leq D \leq 1)$(0\le D\le 1)$, where $K_1$$K_1$K_(1)K_{1}$K_1$ represents the bulk modulus while $K_2$$K_2$K_(2)K_{2}$K_2$ and $K_3$$K_3$K_(3)K_{3}$K_3$ are EOS constants that are derived from PMMA Hugoniot data [19]. $\mu$$\mu$mu\mu$\mu$ is the volumetric strain; ${\varrho }_0$${\varrho }_0$ϱ_(0)\varrho_{0}${\varrho }_0$ and $\varrho$$\varrho$ϱ\varrho$\varrho$ are the reference and current density values; and $\mathrm{\Delta }P$$\mathrm{\Delta }P$Delta P\Delta P$\mathrm{\Delta }P$ is an additional pressure term accounting for bulking and dilatation effects during fracture:

where $\beta$$\beta$beta\beta$\beta$ is the elastic energy loss parameter and $\mathrm{\Delta }U$$\mathrm{\Delta }U$Delta U\Delta U$\mathrm{\Delta }U$ is the energy increment.
Regarding strength and damage, the normalized stresses ${\sigma }_i^{\ast },{\sigma }^{\ast }$${\sigma }_i^{\ast },{\sigma }^{\ast }$sigma_(i)^(**),sigma^(**)\sigma_{i}^{*}, \sigma^{*}${\sigma }_i^{\ast },{\sigma }^{\ast }$, and ${\sigma }_f^{\ast }$${\sigma }_f^{\ast }$sigma_(f)^(**)\sigma_{f}^{*}${\sigma }_f^{\ast }$ out of the von Mises equivalent stresses ${\sigma }_i,\sigma$${\sigma }_i,\sigma$sigma_(i),sigma\sigma_{i}, \sigma${\sigma }_i,\sigma$, and ${\sigma }_f$${\sigma }_f$sigma_(f)\sigma_{f}${\sigma }_f$ express the strength of intact, damaged, and fractured PMMA, respectively, as:

where $A,B,C,M,N$$A,B,C,M,N$A,B,C,M,NA, B, C, M, N$A,B,C,M,N$ are material constants, $T^{\ast }$$T^{\ast }$T^(**)T^{*}$T^{\ast }$ is the maximum tensile hydrostatic pressure that the material can endure, $p^{\ast }$$p^{\ast }$p^(**)p^{*}$p^{\ast }$ is the normalized pressure, and ${\epsilon }^{\cdot \ast }$${\epsilon }^{\cdot \ast }$epsi^(***)\varepsilon^{\cdot *}${\epsilon }^{\cdot \ast }$ is the normalized of equivalent strain rate through dividing the strain rate by a reference value equal to $1{\mathrm{s}}^{-1}$$1{\mathrm{s}}^{-1}$1s^(-1)1 \mathrm{~s}^{-1}$1{\mathrm{s}}^{-1}$. Normalization is performed with respect to the Hugoniot elastic limit (HEL):

where ${\sigma }_{HEL},p_{HEL}$${\sigma }_{HEL},p_{HEL}$sigma_(HEL),p_(HEL)\sigma_{H E L}, p_{H E L}${\sigma }_{HEL},p_{HEL}$, and $s_{HEL}$$s_{HEL}$s_(HEL)s_{H E L}$s_{HEL}$, respectively, are the equivalent stress, volumetric pressure and deviatoric stress at HEL while $P$$P$PP$P$ is the actual pressure and $T$$T$TT$T$ is the tensile strength.
Damage accumulation was calculated incrementally as:

where $\mathrm{\Delta }{\epsilon }_p$$\mathrm{\Delta }{\epsilon }_p$Deltaepsi_(p)\Delta \varepsilon_{p}$\mathrm{\Delta }{\epsilon }_p$ is the plastic strain increment, and ${\epsilon }_{p}f$${\epsilon }_{p}f$epsi_(p)f\varepsilon_{p} f${\epsilon }_{p}f$ is the strain to fracture expressed as:

where $D_1$$D_1$D_(1)D_{1}$D_1$ and $D_2$$D_2$D_(2)D_{2}$D_2$ are PMMA fracture constants and when $D$$D$DD$D$ is equal to 1 , the material is already fractured.

The experimental testing evidence in [19] shows that both types of damage; interlaminar and intralaminar, occur at a high strain rate during blast loading, so a faithful model should contain them in numerical representation. Using continuum shell elements, the delamination damage was modelled using the Virtual Crack Closure Technique (VCCT), while the intralaminar damage was modeled according to the Hashin failure criterion. To avoid mesh distortion problems, a Coupled EulerianLagrangian (CEL) analysis was applied, in which the Copper vapor at the Eulerian part region was simulated using Eulerian mesh while the PMMA was represented according to the Lagrangian formulation. The application of the CEL analysis allows for efficient contact conditions between the Eulerian and Lagrangian parts, allowing for the efficient transfer of pressure between the PMMA/vapor phases, especially after the fully fractured PMMA elements are excluded from the analysis.
For VCCT-based interlaminar damage, a general self-contact interaction with default interaction properties was defined at the initial step for the tube mesh which consists of the PMMA tube and the composite thin cylinder. In Step 1, this general contact was modified to include the 9 contact surface pairs of the 10 composite layers, ensuring that each ring surface interacted with the opposite surface of the adjacent ring. An interaction property incorporating both cohesive behavior and VCCT debonding criteria was defined for these nine new contact pairs. Potential crack surfaces are bonded in Abaqus/Explicit using the general contact's definition. The capability is applied through a pure main-secondary formulation. To explicitly identify the crack tips, the predefined crack surfaces are assumed to be initially partially bonded. Defining the initially bonded node set, a contact clearance is first defined, then assigned to each of the 9 general contact pairs that each of them consists of two single-sided surfaces. This is how the initial crack is defined. The section that is not bonded functions like a typical contact surface. It is assumed that the nodes in the node set are initially bonded in every direction.
The specified node set is part of the secondary surface, which has to extend beyond it. If not, the crack cannot spread, and the surfaces cannot not debond. The bonded set is shown in Fig. 5, illustrating the initial created crack. This initial clearance value is assigned to each of the 9 contact pairs. The definition of the crack propagation capability was finished by defining a cohesive behavior, based on fracture, to the surface interaction. It was assigned to each pair of initially partially bonded surfaces to initiate the crack propagation. Between these two surfaces, crack propagation takes place when the fracture criterion is satisfied. The of the bonds' elastic behavior is also described by cohesive behavior.

When the fracture criterion parameter $(f)$$(f)$(f)(f)$(f)$ of VCCT approaches 1.0 within a specified tolerance, the crack-tip node starts debonding as according to:

where $f_{LL}=1$$f_{LL}=1$f_(LL)=1f_{L L}=1$f_{LL}=1$ is $f$$f$ff$f$ upper limit and $f_{UL}=1+f_{\text{tol}}$$f_{UL}=1+f_{\text{tol }}$f_(UL)=1+f_("tol ")f_{U L}=1+f_{\text {tol }}$f_{UL}=1+f_{\text{tol }}$ is the upper limit, where $f_{\text{tol}}$$f_{\text{tol }}$f_("tol ")f_{\text {tol }}$f_{\text{tol }}$ is the tolerance parameter whose default value is taken as 0.20 , the default in Abaqus. In the general case when three modes of fracture are considered, $f$$f$ff$f$ is defined as:

where the nodal equivalent strain energy-release rate is $g_{\text{equiv}}$$g_{\text{equiv }}$g_("equiv ")g_{\text {equiv }}$g_{\text{equiv }}$, and the critical equivalent strain energy release rate is $g_{\text{C,equiv}}$$g_{\text{C,equiv }}$g_("C,equiv ")g_{\text {C,equiv }}$g_{\text{C,equiv }}$ which is determined by the interface bond strength and the user-specified mode-mix criterion. In this study $g_{G,\text{equiv}}$$g_{G,\text{ equiv }}$g_(G," equiv ")g_{G, \text { equiv }}$g_{G,\text{ equiv }}$ is evaluated by the formula of Benzeggagh and Kenane (BK) [26], as in Eqns. (10) and (11), where $g_I,g_{II}$$g_I,g_{II}$g_(I),g_(II)g_{I}, g_{I I}$g_I,g_{II}$ and $g_{III}$$g_{III}$g_(III)g_{I I I}$g_{III}$ are the energy-release rate values of the three modes of fracture, while $g_{IC},g_{IIC}$$g_{IC},g_{IIC}$g_(IC),g_(IIC)g_{I C}, g_{I I C}$g_{IC},g_{IIC}$, and $g_{IIIC}$$g_{IIIC}$g_(IIIC)g_{I I I C}$g_{IIIC}$ are their critical values, correspondingly. $g_{IC}=950\mathrm{J}/{\mathrm{m}}^2$$g_{IC}=950\mathrm{J}/{\mathrm{m}}^2$g_(IC)=950J//m^(2)g_{I C}=950 \mathrm{~J} / \mathrm{m}^{2}$g_{IC}=950\mathrm{J}/{\mathrm{m}}^2$ and $g_{\text{IIC}}=g_{\text{IIIC}}=1250\mathrm{J}/{\mathrm{m}}^2$$g_{\text{IIC }}=g_{\text{IIIC }}=1250\mathrm{J}/{\mathrm{m}}^2$g_("IIC ")=g_("IIIC ")=1250J//m^(2)g_{\text {IIC }}=g_{\text {IIIC }}=1250 \mathrm{~J} / \mathrm{m}^{2}$g_{\text{IIC }}=g_{\text{IIIC }}=1250\mathrm{J}/{\mathrm{m}}^2$ are the chosen critical energy-release rates according to [27] because of the similarity in used materials.

In intralaminar damage analysis based on Hashin's damage criteria, the fibers in fiber-reinforced materials are usually thought to be parallel, and degradation at a point is referred to as damage initiation. Material stiffness deterioration is a crucial sign of damage in the analysis of composite materials reinforced with fibers. Elastic-brittle behavior is dominant in unidirectional fiber-reinforced composite materials. Damage to such materials occurs without significant plastic deformation. As a result, when simulating the behavior of these materials, plasticity is often overlooked.
The initiation failure criteria of fiber tension, fiber compression, matrix tension and matrix compression, respectively, have the following general forms [22]:

where $X^t$$X^t$X^(t)X^{t}$X^t$ is the longitudinal tensile strength; $X^c$$X^c$X^(c)X^{c}$X^c$ is the longitudinal compressive strength; $Y^t$$Y^t$Y^(t)Y^{t}$Y^t$ is the transverse tensile strength; $Y^c$$Y^c$Y^(c)Y^{c}$Y^c$ is the transverse compressive strength; $S^L$$S^L$S^(L)S^{L}$S^L$ is the longitudinal shear strength; $S^T$$S^T$S^(T)S^{T}$S^T$ is the transverse shear strength; $\alpha$$\alpha$alpha\alpha$\alpha$ is the coefficient of shear stress contribution to fiber tensile-initiation criterion; and ${\sigma }_{11}{}^{\wedge }{\sigma }_{22}{}^{\wedge }{}_{12}$${\sigma }_{11}{}^{\wedge }{\sigma }_{22}{}^{\wedge }{}_{12}$sigma_(11)^sigma_(22)^_(12)\sigma_{11}{ }^{\wedge} \sigma_{22}{ }^{\wedge}{ }_{12}${\sigma }_{11}{}^{\wedge }{\sigma }_{22}{}^{\wedge }{}_{12}$ are the effective components of stress tensor $({\sigma }^{\wedge }=M^{\ast }\sigma )$$\left({\sigma }^{\wedge }=M^{\ast }\sigma \right)$(sigma^=M^(**)sigma)\left(\sigma^{\wedge}=M^{*} \sigma\right)$({\sigma }^{\wedge }=M^{\ast }\sigma )$, where $\sigma$$\sigma$sigma\sigma$\sigma$ is the true stress, and $M$$M$MM$M$ is the damage operator. $M$$M$MM$M$ is determined by Eqn. (13):

where $d_f,d_m$$d_f,d_m$d_(f),d_(m)d_{f}, d_{m}$d_f,d_m$ and $d_s$$d_s$d_(s)d_{s}$d_s$ are the fiber, matrix, and shear damage parameters; while the superscripts $t$$t$tt$t$ and $c$$c$cc$c$ refer to tensile and compressive effects, respectively. The values of $X^t,X^c,Y^t,Y^c,S^L$$X^t,X^c,Y^t,Y^c,S^L$X^(t),X^(c),Y^(t),Y^(c),S^(L)X^{t}, X^{c}, Y^{t}, Y^{c}, S^{L}$X^t,X^c,Y^t,Y^c,S^L$ and $S^T$$S^T$S^(T)S^{T}$S^T$ used in the current FE study are attached in Tab. 4. Regarding intralaminar damage propagation, material properties must be specified in a local coordinate system that is userdefined, where the fiber direction is aligned with the local 1 -direction. Due to matrix microcracking, the nonlinear behavior of the matrix-which includes both stiffness and plasticity degradation-dominates the shear response for bidirectional fiber-reinforced composites. Orthogonality between the fiber directions is assumed. With the local 1 - and 2 -direction aligned with the fiber directions, material properties must be specified in a user-defined local coordinate system. Damaged elasticity characterizes the material response along the fiber orientations. In addition to distinguishing between tensile and compressive fiber failure modes, the model takes into account varying initial (undamaged) stiffness in tension and compression. To define the undamaged response, therefore, bilaminar elasticity is preferred.
The material stiffness is reduced in accordance with the given damage evolution law for each specific damage initiation criterion once it is met. There is no degradation of the material stiffness without assigning a law damage evolution. The rate at which the material stiffness degrades after meeting the appropriate initiation criteria is described by the damage evolution law. The stress tensor of the material at any point during the analysis is related to the elastic strain ( ${\epsilon }^{el}$${\epsilon }^{el}$epsi^(el)\varepsilon^{e l}${\epsilon }^{el}$ ) by:

where $C_d$$C_d$C_(d)C_{d}$C_d$ is the matrix-damaged elasticity. Hashin damage evolution in this FE study is energy based, and the values of fracture energy are attached in Tab. 5. According to Abaqus, the deterioration of the stiffness of the matrix and fiber components can be described as a function of 3 scalar damage parameters, which correspond to the current states of shear, matrix, and fiber damage. It is necessary to use the Hashin's damage model with plane-stress elements, such as the continuum shell elements used for the composite p.v. rings throughout this study.

The physics of wire explosion is more difficult to model than that of the application of EOS for well-studied chemical explosives because it involves a complex combination of phase transitions [28]. Previous research has proposed the governing equations for wire explosion, which include a circuit equation combined with wide-range EOS, magnetohydrodynamic equations, and electrical conductivity for copper [29]. However, in this study, a simplified model is employed for the sake of computational efficiency and the viability of fracture analysis of composites.
It is assumed that a uniform column of vaporized products forms at the moment of explosion, and that this column expands to produce an equivalent cylindrical shock wave. One can estimate the initial conditions of this state by consulting existing literature. For example, the critical density ( $\varrho$$\varrho$ϱ\varrho$\varrho$ ), boundary velocity, and plasma pressure of a $500\mu \mathrm{m}$$500\mu \mathrm{m}$500 mum500 \mu \mathrm{~m}$500\mu \mathrm{m}$ copper wire at the beginning of explosion in water with a charging voltage of 30 kV were estimated to be $2600\mathrm{kg}/{\mathrm{m}}^3,2000\mathrm{m}/\mathrm{s}$$2600\mathrm{kg}/{\mathrm{m}}^3,2000\mathrm{m}/\mathrm{s}$2600kg//m^(3),2000m//s2600 \mathrm{~kg} / \mathrm{m}^{3}, 2000 \mathrm{~m} / \mathrm{s}$2600\mathrm{kg}/{\mathrm{m}}^3,2000\mathrm{m}/\mathrm{s}$, and 10 GPa , respectively [30]. In this case, at the time of explosion, the melted wire radius is roughly 0.46 mm . Furthermore, [31] reports an experimentally measured shock wave velocity of approximately $2000\mathrm{m}/\mathrm{s}$$2000\mathrm{m}/\mathrm{s}$2000m//s2000 \mathrm{~m} / \mathrm{s}$2000\mathrm{m}/\mathrm{s}$ in the ambient air.
The specific energy deposited into the wire ( $E_d$$E_d$E_(d)E_{d}$E_d$ ) determines the shock wave peak pressure and decay time. For example, $E_d=19.44\times {10}^3\mathrm{kJ}/\mathrm{kg}$$E_d=19.44\times {10}^3\mathrm{kJ}/\mathrm{kg}$E_(d)=19.44 xx10^(3)kJ//kgE_{d} =19.44 \times 10^{3} \mathrm{~kJ} / \mathrm{kg}$E_d=19.44\times {10}^3\mathrm{kJ}/\mathrm{kg}$ was measured in [32] for a copper wire with a $25\mu \mathrm{m}$$25\mu \mathrm{m}$25 mum25 \mu \mathrm{~m}$25\mu \mathrm{m}$ diameter and a 20 kV discharge voltage at the air breakdown point. In [20], the deposited energy at various charging voltages is shown for wire diameters up to $240\mu \mathrm{m}$$240\mu \mathrm{m}$240 mum240 \mu \mathrm{~m}$240\mu \mathrm{m}$. By extending these findings to a wire with a diameter of $500\mu \mathrm{m}$$500\mu \mathrm{m}$500 mum500 \mu \mathrm{~m}$500\mu \mathrm{m}$ and a 20 kV discharge, $E_d=20\times {10}^3\mathrm{kJ}/\mathrm{kg}$$E_d=20\times {10}^3\mathrm{kJ}/\mathrm{kg}$E_(d)=20 xx10^(3)kJ//kgE_{d}=20 \times 10^{3} \mathrm{~kJ} / \mathrm{kg}$E_d=20\times {10}^3\mathrm{kJ}/\mathrm{kg}$.
Under the adiabatic assumption, the ideal gas EOS is used to model the vaporized wire:

where $\gamma$$\gamma$gamma\gamma$\gamma$ denotes the isentropic index, while $E_m(\mathrm{kJ}/\mathrm{kg})$$E_m(\mathrm{kJ}/\mathrm{kg})$E_(m)(kJ//kg)E_{m}(\mathrm{~kJ} / \mathrm{kg})$E_m(\mathrm{kJ}/\mathrm{kg})$ is the internal specific energy related to $E_d$$E_d$E_(d)E_{d}$E_d$ owing to the transformation of electric energy. The ideal gas EOS assumes the explosion products behave as a perfect gas with constant specific heats. This simplification ignores the complex phase transitions (solid-liquid-vapor-plasma) and the real gas effects of the copper vapor plasma at extremely high pressures and temperatures. The primary implication is that the initial peak pressure and the detailed temporal decay of the pressure pulse may not be perfectly captured. However, for the purpose of this study, which is to simulate the mechanical response of the composite structure to a shock pulse, the key requirements are to accurately replicate the impulse (integral of pressure over time) and the general pulse shape that loads the PMMA wall.
The published data in [19],[20],[29-32] is used to determine the EOS parameters of the exploded wire: the initial value of $\varrho$$\varrho$ϱ\varrho$\varrho$ is equal to $2600\mathrm{kg}/{\mathrm{m}}^3,E_d=20\times {10}^3\mathrm{kJ}/\mathrm{kg}$$2600\mathrm{kg}/{\mathrm{m}}^3,E_d=20\times {10}^3\mathrm{kJ}/\mathrm{kg}$2600kg//m^(3),E_(d)=20 xx10^(3)kJ//kg2600 \mathrm{~kg} / \mathrm{m}^{3}, E_{d}=20 \times 10^{3} \mathrm{~kJ} / \mathrm{kg}$2600\mathrm{kg}/{\mathrm{m}}^3,E_d=20\times {10}^3\mathrm{kJ}/\mathrm{kg}$, and the isentropic index ( $\gamma$$\gamma$gamma\gamma$\gamma$ ) was set to 1.4. The specific energy of copper is equal to $50\mathrm{kJ}/\mathrm{kg}$$50\mathrm{kJ}/\mathrm{kg}$50kJ//kg50 \mathrm{~kJ} / \mathrm{kg}$50\mathrm{kJ}/\mathrm{kg}$. In [19], further explosion experiments on pure PMMA samples were carried out to find the best fit between the numerical simulations and the experimental free-surface velocity history in order to establish the proper Equation of State (EOS), calibrating the parameters. Two discharge voltage scenarios were chosen based on earlier research on unidirectional composite cylinders: $20\mathrm{kV}(1.2\mathrm{kJ})$$20\mathrm{kV}(1.2\mathrm{kJ})$20kV(1.2kJ)20 \mathrm{kV}(1.2 \mathrm{~kJ})$20\mathrm{kV}(1.2\mathrm{kJ})$ and $25\mathrm{kV}(1.875\mathrm{kJ})$$25\mathrm{kV}(1.875\mathrm{kJ})$25kV(1.875kJ)25 \mathrm{kV}(1.875 \mathrm{~kJ})$25\mathrm{kV}(1.875\mathrm{kJ})$ [19]. In the EOS model, the internal energy $E_m$$E_m$E_(m)E_{m}$E_m$ was connected to the deposited energy $E_d$$E_d$E_(d)E_{d}$E_d$ with a correction factor that took conversion losses into account. Trial-and-error simulations were carried out using an iterative methodology, altering the pressure-internal energy relationship to correspond with experimental findings. The pressure history at a fixed Eulerian coordinate on the channel boundary confirmed consistency with previously reported experimental data. The $25\mathrm{\%}$$25\mathrm{\%}$25%25 \%$25\mathrm{\%}$ increase in discharge voltage ( 5 kV ) resulted in a $\sim 56\mathrm{\%}$$\sim 56\mathrm{\%}$∼56%\sim 56 \%$\sim 56\mathrm{\%}$ rise in stored energy, with a much greater energy conversion efficiency at higher voltages. Pressure differences between the two situations decreased over time, and the peak pressure in PMMA close to the channel was approximately 1.6 GPa for 20 kV , which was in good agreement with previous experimental observations. This process ensures that the resulting pressure pulse applied to the Lagrangian structure (PMMA/composite) in the current study FE simulations is empirically validated against physical experiments, mitigating the inaccuracies of the simplified EOS.
The discrete field of the initial Eulerian volume-fraction of the copper wire during the explosion was defined as a uniform volume-fraction-type initial condition to facilitate the CEL interaction for modeling the shock wave. To track material flow, Abaqus/Explicit employs the volume-of-fluid method, which determines the Eulerian Volume Fraction (EVF) for each element. An EVF of one indicates complete material presence, while zero signifies its absence. When the volume fraction is less than one, the remaining portion of the element is filled with "void" material, which has no mass or strength. In this
study, the void material represents air, with negligible ambient pressure of only 0.1 MPa . The initial volume fraction value of 1 was assigned to the Eulerian elements located within the central hole of the PMMA cylinder, occupying a 0.92 mm diameter cylinder representing the copper wire region. Because the Eulerian elements are cubic, the volume fraction values ranged from 0.283 to 1 depending on whether the whole element is inside the wire region completely or partially.
CEL formulation governing the interaction between the Eulerian (vapor) and Lagrangian (PMMA and composite) domains is modeled in Abaqus/Explicit by an initial-condition-defined pressure pulse acting on the PMMA insert at the region Eulerian part inside the PMMA-composite region. Along with defining the Eulerian volume fraction, the vapor applies a pressure load P as a surface traction at the fluid-structure interface created by this conventional CEL method.
On the other hand, the Eulerian material flow is constrained by the motion of the Lagrangian boundaries, which inhibits penetration. This approach replicates the experimental transfer of energy via a shock wave through the PMMA to the composite vessel by accurately confining the expanding vapor within the channel and ensuring realistic transfer of the timedependent pressure pulse from the expanding vapor to the structure.

A dynamic explicit step of ${10}^{-5}$${10}^{-5}$10^(-5)10^{-5}${10}^{-5}$ seconds, incorporating non-linear geometry calculations, was used to model the internal blast loading. The output requests included stress parameters, Hashin's failure criteria, and specific general contact parameters such as bond status (BDSTAT), energy release rate (ENRRT) and effective energy release rate (EFENRRTR).
Numerical viscosity is included to alleviate numerical stability issues related to explicit integration schemes. To dampen high frequencies in Abaqus, the linear bulk viscosity pressure is introduced and formulated as follows [17]: