Novel material model to predict the residual strength of a composite overwrapped pressure vessel after impact

Composite overwrapped pressure vessels are typically used to store gases such as hydrogen under high pressure. Internal damage due to an impact may reduce the strength of the vessel, through which the internal pressure at failure, the so-called burst pressure, may drop significantly. Advanced material models can be used to predict to what extend the burst pressure of a vessel decreases. However, current studies do not consider the complex mechanisms that cause a reduction of the burst pressure for vessels with a wall thickness over 20 mm (common in the automotive industry), nor do they correctly take into account the reduced stiffness during reloading. The failure mechanisms in these thick-walled vessels under impact loading localizes in shear bands that consist of fiber kinking and matrix damage. Fiber fracture due to fiber kinking is identified as the most important mechanism that reduces the burst pressure. In this paper, a new material model is proposed that adequately captures these mechanisms. The model is validated using impact loads and subsequent burst tests on pressure vessels. It is shown that the model is able to predict the burst pressure after impact well.


Introduction
Fuel cell electric vehicles hold great potential, alongside battery electric vehicles, to decarbonize our future mobility. In a fuel cell, oxygen from the air and on-board stored hydrogen are converted into electricity and water vapor. The electricity is used to power electric motors and the water vapor is a harmless exhaust product. The required hydrogen is stored in composite overwrapped pressure vessels (COPVs) under a high pressure of up to 70 MPa (700 bar). To withstand this internal pressure, vessels with automotive dimensions typically have a wall thickness exceeding 20 mm.
The use of hydrogen pressure vessels in vehicles requires a thorough understanding of any inflicted damage in the case of a crash. In [1], an experimental study was presented in which the occurring damage mechanisms in thick-walled COPVs were investigated after impact loads. After being subjected to an external contact load, the residual strength of the pressure vessel was determined using a so-called burst test. An increasing internal pressure was applied until the vessel failed (burst). The reduction of the burst pressure compared to an undamaged vessel is a measure of the induced global damage. The main conclusion of this study was that the residual burst pressure of a thick-walled pressure vessel decreases significantly when shear bands consisting of matrix damage and fiber kinking are formed. The associated fiber fracture due to fiber kinking was identified as the main cause of the measured reduction in burst pressure [1].
To improve the understanding of the inflicted damage mechanisms, an approach was developed to model a COPV with a wall thickness of over 20 mm under impact conditions [2]. It was shown that the model, which only represents interlaminar damage (delamination), was able to capture the global behavior well. Progressive intralaminar damage was not incorporated in this study. Instead, a fiber compressive failure criterion was used to determine the onset and location of failure. The failure criterion accurately predicted that the shear bands with fiber kinking initiate under an angle with respect to the indenter. The size of the damaged region in the simulations exceeded the experimental observations, which may be due to the absence of a damage evolution model which leads to damage localization.
In order to also predict a reduction of the residual strength, the modeling approach proposed in [2] should be extended to describe the evolution of intralaminar damage. The key ingredient of the damage evolution model is the loss of strength and stiffness due to fiber kinking and subsequent fiber fracture. It is known from literature that the initiation of fiber kinking does not lead to fiber fracture instanteneously [3,4]. Various authors have studied the fiber kinking mechanism in detail, of which the main results are summarized next.
Pimenta et al. [3] conducted an experimental-numerical study on the fiber kinking mechanism under uni-directional compressive loading. In this study, the following regimes were identified: (i) the elastic domain, (ii) the softening domain and (iii) the fiber failure domain. In the elastic domain, the fibers and matrix deform elastically. When reaching the peak load, the matrix starts to yield, causing a reduction of fiber support. The fibers start to rotate which causes global softening (softening domain). At large rotations, the outer fibers start to break and the fiber failure domain is initiated.
A similar trend was identified by Han et al. [4], who used a numerical model to study the fiber kinking mechanism under longitudinal compression conditions. Their results show the formation of fiber kinking after a linear-elastic regime. After the occurrence of fiber kinking, the stiffness is reduced as the fiber-matrix interface debonds and the matrix subsequently yields in shear. The peak load is reached as the interface completely debonds and matrix damage expands across the width of the bands. The fiber failure domain starts shortly after reaching the peak load, as the stress in the increasingly bending fibers reaches the fiber strength.
Nizolek et al. [5] reported a similar behavior, in which kink band formation does not immediately cause global softening. In their study, kink band propagation in Cu-Nb nanolaminate composites was investigated since the kink band formation for this composite is stable. A unidirectional compressive load was considered and digital image correlation was used to track the kink band initiation and propagation. No distinct slope change could be recognized in the stress-strain curve upon visible identification of the kink band. A sharp load drop was first identified after the kink band had grown across the specimen width.
Vogler and Kyriakides [6] studied the initiation and growth of kink bands under biaxial (combined shear and compression) loading conditions. The sample was loaded in compression up to a certain stress level. Subsequently, whilst maintaining the compressive stress, shear loading was applied until a kink band was formed. The growth of the shear band under shear loading was tracked using a camera. The shear band continued to grow under an increasing shear load. As it grew across the entire sample, the shear stress reached a limit and then decreased. The authors reported that the fibers did not break under shear-dominated fiber kink initiation and growth.
These studies summarize the complex mechanisms that take consecutively place at the micro-scale during fiber kinking. Nevertheless, these phenomena are typically not taken into account in most homogenized continuum models available in literature. Several authors presented numerical models that were used to assess the influence of an impact on the residual strength of a pressure vessel, e.g. [7][8][9][10]. However, since these studies focused on relatively thin vessels (up to 6 mm thick), where the kinking occurs due to a pure compressive stress, shear-dominated damage as observed for thick-walled composite, was not reported. Consequently, the available material models generally focused on two-dimensional failure criteria only and the role of fiber kinking was neglected or overly simplified.
The fiber kinking mechanism is included in constitutive models in commercially available Finite Element codes. In LS-DYNA, the mechanism is incorporated in the commonly used material models MAT_261 and MAT_262. These models are based on to the work of Pinho et al. [11,12] and Maimȡ et al. [13,14], respectively. Although fiber kinking is accounted for in these implementations, the degradation of the tensile properties due to fiber fracture is missing (MAT_261) or the initiation criteria for fiber kinking only relies on a two-dimensional stress state (MAT_262). As a result, these material models cannot be used to predict the residual burst pressure in thick-walled pressure vessels.
In this paper, we present a novel material model that is able to predict the residual strength of a thick-walled composite structure that is damaged by an external contact load. The model considers the complete three-dimensional stress state to initiate fiber kinking and subsequent fiber fracture.
The paper is ordered as follows. First, the computational model for the assessment of a pressure vessels is introduced. This model enables an external contact load to be applied on a pressure vessel and a burst test to be conducted subsequently. Measures are taken in order to optimize the computational efficiency whilst minimizing the loss of accuracy of the model. The details of the material model are presented in Section 3. The material and computational model are then combined in the final section, in which the results are presented. The influence of the governing material properties is shown first. Experimental results for some impact conditions are subsequently used to determine the values of these material properties. Finally, it is shown that the material model is able to predict the residual strength of this pressure vessel for different impact conditions.

Computational pressure vessel model
A numerical model is used to study the effects of an external contact loading on a composite-overwrapped pressure vessel under either quasistatic or dynamic conditions. After this contact loading, a numerical burst test is conducted in order to determine the residual strength of the vessel. The used computational model is a finite element model and uses the explicit LS-DYNA solver. The used model framework and the modeling approaches for the external contact loading and for the burst test are introduced in the following sections.

Model description
The pressure vessel is modeled as a cylinder with the length of the complete pressure vessel, see Fig. 1. The dome and boss area of the vessel are geometrically simplified and represented as a continuation of the cylindrical area.
The composite material is modeled using a combination of thick shell and cohesive elements as shown in the detail of Fig. 1. A unidirectional elasto-plastic orthotropic material model with failure is used to describe the intralaminar material behavior in the thick shell elements. The used material model is presented in Section 3. Multiple integration points are defined throughout the thickness of the thick shell elements, representing the single layers of the composite. At each of these integration points, the local lay-up of the laminate is captured by corresponding orientation of the material model.
The cohesive elements that are placed in between the thick shell elements are used to model delamination initiation and growth. This damage mechanism was identified as the first mechanism that can be recognized when a pressure vessel is loaded under external contact loading [1]. A bilinear traction-separation law with a mixed mode initiation criterion is used to describe the response of the cohesive elements.
To reduce computational costs, the cylinder is divided in multiple regions with various levels of detail: from fine underneath the indenter (shown in blue in Fig. 1) to coarse at the edges (shown in orange). The mesh size and the number of cohesive layers are different per region. In the most detailed region, the mesh size is approximately 2.5 mm and about 40 cohesive layers are in place. In the most coarse region, the mesh size is up to 25 mm and only 2 cohesive layers are used. This geometrical modeling approach was previously introduced in [2], where a complete description is provided.

Boundary conditions and loading
The simulation consists of multiple steps that are performed consecutively. First, an external contact load is applied using the indenter to locally damage the vessel wall. After reaching the maximum displacement, the indenter is moved back to its original position in order to unload the vessel wall. When the indenter is back at its original position, the numerical contact definition between the vessel and the indenter is deactivated. Finally, a numerical burst test is conducted to determine the residual burst pressure of the vessel after it has been subjected to the external contact load.
The cylinder is positioned on support blocks, which are modeled as a rigid material, see Fig. 1. An indenter (modeled as rigid) is used to apply the external contact loading. The indenter is either displacing slowly with a gradually increasing prescribed velocity (in the case of a quasistatic simulation) or moving with an initial velocity (in the case of a dynamic simulation).
The application of internal pressure induces stresses in the hoop and axial directions of a vessel. The hoop component is caused by the radial pressure in the cylindrical region of the vessel, whereas the axial component is due to the pressure at the end caps (domes and bosses). In the simulation, the radial component can be directly applied as a normal pressure on the inner surface of the cylinder, as illustrated in the crosssection view in Fig. 2.
Since the domes of the vessel are not modeled geometrically, the corresponding pressure load has to be applied at the extremities of the cylinder. Accordingly, the nodes at either end of the cylinder are tied to a single control node at their respective centers. The axial degree-offreedom of the nodes is constrained using an interpolation, i.e. the average axial displacement of the nodes at the cylinder face is equal to the control node's axial displacement. An axial force is applied to the central node, which is equal to the projected inner area of the dome multiplied by the internal pressure (πR 2 in p). The application of the axial internal pressure component is shown in Fig. 2 for one of the two faces.
The experimental burst tests were performed hydraulically using the pressure ramp described in [1]. The burst pressure is reached in approximately 15 minutes. In the simulation, the internal pressure is increased from 0 to 250 MPa in a simulation time span of only 0.1 seconds to minimize computation time. Although the loading rate is exceeding the loading rate of the experiment, the influence of dynamic effects on the burst behavior is expected to be negligible. This is confirmed by the kinetic energy being small (≪ 1%) compared to other energy measures and since the material models used are not rate-dependent.
The experimental burst pressure is the pressure at which the vessel no longer withstands the load and starts to rupture [1]. In the simulations, the pressure for which a sudden increase of dissipated internal energy is measured is taken as the burst pressure. This is later discussed in more detail in Section 4.2 and shown in Fig. 13.

Material model
Each layer in the composite-overwrapped pressure vessel can be considered as a unidirectional material consisting of fibers embedded in an epoxy matrix, see Fig. 3.
In order to describe the inhomogeneous damage behavior in composites, one needs to distinguish between transverse (matrix-dominated) and longitudinal (fiber-dominated) damage. Moreover, different damage mechanisms may take place for tensile-or compressive-dominated loading. As a result, the models takes four different damage mechanisms into account, which are next refered to as: (I) transverse tensile, (II) transverse compressive, (III) longitudinal tensile and (IV) longitudinal compressive damage. These damage mechanisms are sketched in Fig. 4. Note that these illustrations are simplified for clarity; i.e. other stress components that may contribute to these mechanisms, e.g. shearloading, are not shown.
(I) Transverse tensile loading is assumed to lead to a fracture plane in the matrix-material when a critical stress state is exceeded. This fracture plane develops perpendicular to the loading direction as illustrated in Fig. 4. Upon load reversal (from tensile to compressive), it is assumed that the crack faces close and are able to transmit stresses again. Hence, transverse tensile damage is assumed to have no effect on a subsequent transverse compressive loading. Fig. 1. Illustration of the computational model geometry. The pressure vessel is simplified as a cylinder, which is modeled using regions with various degrees of detail: Fine (blue), medium (green) and coarse (orange). The varying in-plane mesh size is shown in the overview, the varying number of cohesive layers in the cross-section. The definition of various layer orientations using integration points across the thickness is illustrated in the element view on the bottom right. Image based on illustration from [2]. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 2.
Numerical application of internal pressure. The radial component is applied as an internal surface pressure, shown in the cross-section at the right. The axial component is applied as a force which is transferred to the nodes at the cylinder ends using an interpolation constraint. Image based on illustration from [2]. For most composites, this fracture angle is approximately 53 ∘ [11,13]. The crack is assumed to be able to open under subsequent tensile loading, hence affecting loading in the transverse tensile direction after compressive damage growth. Since fracture in both transverse damage mechanisms is limited to the matrix material, the fibers are expected to remain intact.
(III) Longitudinal tensile loading is assumed to lead to fiber fracture when the tensile strength of the fibers is exceeded. This damage mechanism typically occurs abruptly and the surrounding matrix material is assumed to fail consequently. Moreover, upon load reversal, we assume that the broken fibers can no longer transmit loads.
(IV) As illustrated in Fig. 4, loading in the compressive direction of the fibers is assumed to cause fiber kinking. As reported in literature [15], the responsible mechanism is the collapse of the fibers due to damage of the matrix material. The progression of the fiber kinks consists of two stages that take place in consecutive order: (i) The fibers bend elastically and buckle, which leads to a reduction of the global stiffness and (ii) the outer fibers break as a result of the emerging bending stresses that exceed the fiber strength. When the compressive loading is interrupted during stage (i) and a longitudinal tensile load is applied subsequently, the kinked fibers straighten, regain stiffness and are able to carry the load again. However, when compressive loading enters stage (ii) and fiber fracture occurred, the broken fibers are no longer able to carry load upon load reversal. Therefore, whether or not fiber kinking affects subsequent tensile loading depends on whether the fibers fractured or not during fiber kinking. The material model should be able to account for this phenomenon.

Model framework
In the model presented here, a small strain formulation is assumed. The strain and stress components are described in a local 123-coordinate system, which is shown in Fig. 3. The 1-axis is parallel to the fiber orientation and the 23-plane is the plane of isotropy which is perpendicular to that.
The stresses and strains are expressed in Voigt notation with respect to the basis shown in Fig. 3: Additionally, a so-called effective (undamaged) stress state is defined, which we refer to as σ. This stress vector is related to the elastic component of the strain ε el by an orthotropic stiffness matrix, which describes the linear elastic material behavior in a homogenized manner, as commonly done in literature [12]: with C the stiffness matrix, which is the inverse of the compliance matrix S: which is defined as: The elastic moduli in the 11-, 22-and 33-direction of the lamina are refered to as E 11 , E 22 and E 33 , respectively. The shear moduli in the 23-, 13-and 12-directions are defined as G 23 , G 13 and G 12 . Finally, the major Poisson's ratios are described as ν 12 , ν 13 and ν 23 .
The elastic strain ε el is related to the total strain vector ε by substraction of the plastic strain ε pl : Under loading, the strain is decomposed in its elastic and plastic components. The suggestions from Pinho et al. [12] are followed here: Elasto-plastic behavior is only included for the 12-component of the lamina. The response is implemented in a similar manner and we therefore refer to [12] for a complete description. Plastic behavior is not incorporated for the other components as the plastic strains in these directions are expected to be small. A continuum damage model is used to describe the degradation of the laminate properties as a result of intralaminar damage. The damaged stress σ is related to the effective (undamaged) stress vector σ by use of the following expression: where I is the unit matrix and D is the damage matrix. The components of the damage matrix will be introduced later. In the following sections the failure initiation criteria and the damage evolution laws are first presented. The definition of the damage components in the damage matrix will be discussed thereafter.

Damage initiation
As introduced before, we distinguish between damage initiation due to longitudinal and transverse damage in the tensile and compressive directions. Different failure criteria are used to predict the initiation of these individual failure mechanisms. Here, the criteria presented by Camanho et al. [15] are used and the equations are incorporated exactly as described there. A summary of the governing equations is repeated here for completeness only. We refer to [15] for a complete overview.

Transverse failure initiation (I+II)
Camanho et al. [15] proposed the following invariant-based criterion to predict transverse failure initiation for composite materials: In here, f 22,T/C is the failure index for transverse tensile (T) or compressive (C) failure and I 1 , I 2 and I 3 are the transversely isotropic stress invariants. The prefered direction used to derive the invariants I is the 11-direction of the lamina, since the material response is invariant with respect to arbitrary rotations around this axis. Invariant I 1 describes transverse shear loading, I 2 in-plane shear loading and I 3 pressure dependent loading. The latter is used to distinguish between tensile and compressive failure. The invariants are defined as: The failure parameters α are calculated from the material strengths under transverse and shear loading. The used relations are included in A.1.

Fiber tensile failure initiation (III)
Following the suggestions from Camanho et al. [15], the maximum allowable strain criterion is used to predict fiber tensile failure: with f 11,T being the failure index for fiber tensile failure, ε 11 the tensile strain in the fiber-dominated direction and ε 11,T the fracture strain in this direction.

Fiber kinking failure initiation (IV)
As for transverse failure, Camanho et al. [15] proposed to evaluate fiber kinking failure initiation by an invariant-based criterion. An initial fiber misalignment is assumed to be present, which is represented using the dashed lines in Fig. 5. The corresponding misalignment angle is referred to as φ 0 . This initial misalignment grows under the applied loading. The resulting total fiber misalignment is illustrated using the continuous line in Fig. 5 and the corresponding misalignment angle is named φ. As the misalignment angle reaches a critical angle, the supporting matrix is damaged, resulting in collapse of the fibers. Since fiber kinking is expected to initiate due to matrix failure, a criterion resembling transverse failure is adopted [15]: with f 11,C being the failure index for fiber kinking failure, α the failure parameters and I the transversely isotropic invariants. The equations for the failure parameters are identical to those used for transverse failure and can be found in A.1. However, the prefered direction of the lamina which is used to calculate the isotropic invariants is different. In the case of fiber kinking, it is the local fiber direction in the misalignment at the moment of fiber kink initiation. The invariants for fiber kinking failure initiation are defined as [15]: I 2 =σ 2 11 sin 2 φcos 2 φ; with the misalignment angle φ being a function of different parameters, among which the compressive longitudinal strength of the lamina X C . We refer to [15] for further details on calculation of the misalignment angle.

Damage progression
When a failure criterion is reached, i.e. when f ii ≥ 1, damage is initiated. The growth of the damage variable is described using a damage evolution law. Let us consider a one-dimensional 'scalar' simplification of Equation (7) in order to illustrate the behavior: where σ is the damaged stress, σ the effective stress and D the damage variable. In this study, a linear softening behavior is adopted for transverse damage, see Fig. 6.
Initial loading is linear up to a stress σ i and a strain ε i , at which instance damage is initiated. The damaged stress (σ) starts to deviate from this point onwards as the damage variable D grows, see Equation (7). The damage evolution law is defined such that the stress linearly degrades to zero at the strain at final failure ε f : Fig. 5. Illustration of the initial and total fiber misalignment angles φ 0 and φ, based on illustrations from [15]. Fig. 6. Illustration of damage evolution. Linear softening behavior is used, which is initiated at stress σ i and strain ε i . The slope of the softening curve is driven by the energy release rate (ERR) G .
As frequently done in literature [12,14], the damage evolution law is governed by an energy release rate (ERR) G as the local nature of damage evolution is susceptible to mesh size dependency. To avoid this mesh size dependency, the ERR divided by the specific element length L is the area under the stress-strain curve, see Fig. 6. This ensures a constant energy released per unit area of crack generated. For linear softening behavior, the strain at final failure ε f is calculated by the following equation [12]: Pinho et al. [12] stated that the strain and stress components used in the above equations should be those that drive the damage variables. Accordingly, different equivalent stresses (σ eq ) and equivalent strains (ε eq ) are defined for each damage mechanism.

Transverse damage evolution (I+II)
The bilinear behavior presented in Fig. 6 is adopted for both tensile and compressive transverse damage. The ERRs are refered to as G 22,T and G 22,C for tensile and compressive failure, respectively. The equivalent strains and stresses are defined as a combination of the transverse in-plane direction (22) and shear in-plane direction (12).

Fiber tensile damage evolution (III)
For fiber tensile failure, damage growth is typically unstable and spreads throughout the ply instantly. Therefore, the damage variable D 11,T is set at a value close to one directly upon failure initiation. A damage value of one is avoided as it would reduce the stiffness to zero, which could be a source of instability issues. Moreover, upon failure initiation, the deletion flag of the integration point is set as true. When this flag is set for a specified amount of integration points within an element, the element is deleted.

Fiber kinking damage evolution (IV)
The damage evolution due to fiber kinking failure is implemented in a different manner. Following the linear behavior before damage initiation and preceeding the linear softening regime, a stress plateau is included, see Fig. 7.
This plateau represents the stress state during fiber kinking, prior to fiber fracture. Different trends are reported for the global behavior at the moment of fiber kinking without fiber fracture in the literature, varying from global softening [3] to global hardening [4,5]. Which mechanism occurs seems to depend on the support of the surrounding matrix material. Since thick wound laminates are considered here, it can be assumed that the surrounding material adequately supports the failing ply. Therefore, a rather stable behavior is expected, which is represented by the plateau.
According to literature, fiber fracture occurs at large deformations as a result of bending of the fibers [3][4][5]. As fiber fracture initiates, the linear softening regime is initiated in the model. Simultaneously, the damage variable D 11,T , representing fiber fracture, gradually grows in the presented model. This is shown in the lower graph in Fig. 7.
The amount of energy absorbed before fiber fracture initiation is controlled by the plateau ERR G Pl 11,C . The variable G Tot 11,C prescribes the total amount of energy that is dissipated, hence implicitly dictating the slope of the linear softening regime. These energies are highlighted in the two upper graphs in Fig. 7.
Three equivalent strains and one equivalent stress describe the plateau and softening behavior: ε i eq , ε p eq , ε f eq and σ i eq , see Fig. 7. The plateau is described using the following Equation, which is evaluated if ε i eq < ε eq ≤ ε p eq : The linear softening regime (ε p eq < ε eq ≤ ε f eq ) is described using following relationship: with: As stated before, the tensile damage variable D 11,T is also progressing in the linear softening regime as we assume fiber fracture to progress, see Fig. 7. The following equation is used for this purpose: The equivalent strain and stress that are used to evaluate both damage variables (D 11,C and D 11,T ) are the compressive strain and stress in the longitudinal direction:

Coupling
The described damage mechanisms and their evolution models are used to degrade the 3D stresses according to Equation (7). Accordingly, we have to introduce the damage matrix D, which is used to relate the damaged stress vector σ to the effective (undamaged) stress vector σ: As can be recognized, damage is only applied to the 11-, 22-and 12directions of the laminate. Degradation of the out-of-plane components is dominated by delamination, which is already accounted for in the model by the use of a cohesive zone model, see Section 2. The degradation of only the in-plane laminate components is common in literature, e.g. [14].
The three damage variables that are used in the damage matrix, D 11 , D 22 and D 12 , are calculated based on the damage variables that were earlier introduced for longitudinal and transverse damage in the tensile and longitudinal directions. For this purpose, load-dependent expressions are used, which are based on the underlying physical considerations.
The expression for the longitudinal damage variable D 11 is given below: The longitudinal direction is only affected by the damage that is initiated in this direction, not by transverse damage. Since fiber kinking without fiber fracture does not affect the tensile loading behavior, the damage variable D 11,C is excluded in the case of tensile loading. Fiber fracture, described by D 11,T is assumed to affect both compressive and tensile loading, and is therefore included for both loading directions.
In the transverse direction, we assume that matrix cracks, developed under tensile loading, close again under subsequent compressive loading. Therefore, the damage variable D 22,T is only included when tensile loads are applied. Matrix damage that developed under compressive loading (D 22,C ) is assumed to affect both tensile and compressive loading. Fiber fracture, represented by D 11,T , is expected to implicitly cause matrix damage and is therefore incorporated for both tensile and compressive loading. These assumptions lead to the following governing equation for the 22-direction: The in-plane shear direction, being matrix-dominated, is considered using a similar rationale as for the 22-direction. However, a different definition is used to determine the sign of the stress state, see the Equations below. The expression is based on the transverse isotropic invariants used by Camanho et al. [15] and also used to determine the sign of the stress state there.

Results
The finite element and material models are used to study the impact on the cylindrical part of a pressure vessel. The considered pressure vessels are manufactured by the company NPROXX. The length L, the outer radius R out and the wall thickness t are approximately 1750 mm, 175 mm and 25 mm, respectively. See Fig. 1 for a definition of the dimensions. About 100 plies with a thickness of 0.25 mm each compose the wall of the vessel. A multi-directional lay-up consisting of orientations between 20 ∘ and 90 ∘ forms the composite wall, with 0 ∘ being the axial direction of the vessel. The same vessel was studied in [2], in which a material characterization program and first modeling steps of these vessels were presented. The material parameters determined in [2] are used here and summarized in Tables A.3 and A.4 in A.2. In [2], the non-linear behavior under in-plane shear loading was measured using tube torsion tests. This curve is used for the elasto-plastic behavior for in-plane shear loading. Moreover, the material characterization program revealed that the constitutive behavior of this specific composite material is almost independent of the applied strain rate. Therefore, the influence of strain rate effects is not considered in this study. The tensile failure strain published by the fiber manufacturer is used (ε 11,T = 0.019).
One quasi-static test (QS) and four dynamic tests (DYN) were conducted using five samples of this pressure vessel, which were subsequently subjected to burst. Additionally, three burst tests were conducted on undamaged vessels (UND). The loading conditions and residual burst pressures are summarized in Table 1. The displacement, force and energy (external work) used throughout this section are normalized for confidentiality reasons. The results are normalized with respect to sample QS-1: The maximum displacement, maximum force and maximum energy measured for this vessel are set to 1.0 and the other data is scaled accordingly. In essence, all displacements, forces and energies are divided by the maximum displacement, force and energy measured for sample QS-1, respectively.
The external work is defined as in [1]: For the quasi-static case, the external work is the mechanical work, which is the integral of the force-displacement curve. For the dynamic case, the external work is defined as the reduction of the kinetic energy of the indenter during the impact. Therefore, the external work includes the mechanical work, but also vibrational energy and energy that is dissipated in the experimental setup.
In order to validate the simulation procedure used to determine the burst pressure (described in Section 2.2), it was conducted for an undamaged pressure vessel. The resulting numerical burst pressure of 170 MPa lies within the spread of the undamaged pressure vessels in Table 1, therefore confirming the numerical burst method and the used failure parameters. Note the experimental spread of up to 15 MPa for new samples, which can be attributed to manufacturing and material variations. The energy release rates that are used to describe the progressive damage were not determined in [2]. These parameters are therefore detemined in the first part of this section. In the second part, the external contact loading case is analyzed. In essence, the mechanical response and initiated damage in the models are compared to the experimental results. The section is concluded with a validation, in which the capability of the model to predict the residual bust pressure is shown.

Parameter identification
The residual burst pressure strongly depends on the tensile fiber damage D 11,T that is initiated during fiber kinking. The damage growth due to fiber kinking is controlled by the following two energy release rates (ERRs) of the material model: G Pl 11,C and G Tot 11,C . The fiber compressive ERR is typically determined using a compact compression test in literature [16]. However, the distinction between fiber kinking without and with fiber fracture is not commonly made. No test methods are therefore available that aim at determination of the ERR until the onset of fiber fracture, which is the parameter G Pl 11,C in the present model. Hence, these parameters were not identified in the material characterization program in [2]. Instead, the parameters are here determined using the pressure vessel tests. First, the actual influence of the two governing material parameters of the material model (G Pl 11,C and G Tot 11,C ) is shown. The combination of these material parameters that best matches the experimental results is presented afterwards.
To study the influence of the material parameters, impact simulations with varying incident energies were conducted followed by a numerical burst test. Impact tests with several energies are performed by varying the impact mass. The impact velocity was kept constant. Fig. 8 shows the resulting residual burst pressure as a function of the incident energy.
The influence of the plateau energy release rate (ERR) G Pl 11,C is studied first. Three cases are considered: One case in which the plateau is absent (indicated as G Pl 11,C = N.A.) and two different values for the ERR (G Pl 11,C = 10 N/mm and G Pl 11,C = 20 N/mm). The first case is a special implementation of the material model, in which the plateau is omitted and the softening trend is immediately started after damage initiation. Hence, both damage variables D 11,C and D 11,T immediately grow after damage initiation. The absence of the plateau causes low energy impacts to already lead to a reduction of residual strength. The steep decrease of the residual strength for low energy impacts is followed by a convergence towards a plateau at about 40 MPa for large incident energies. For the two cases in which a non-zero plateau ERR is used, a critical energy can be recognized, under which the residual strength is not or barely affected. This matches experimental observations. After exceeding the critical energy, the residual strength rapidly decreases until a plateau is again reached, which is higher than the case for which the plateau is absent. Note that the curves for an ERR of 10 and 20 N/mm are similar; the residual burst pressures at low energies are slightly higher when a higher ERR is used, and the results converge afterwards.
The influence of the total ERR G Tot 11,C is studied in the lower graph in Fig. 8. The value is varied between 25 and 100 N/mm, whilst the plateau ERR G Pl 11,C is set at 10 N/mm. In all cases, a critical incident energy can be recognized, under which the residual burst pressure is not affected. This critical energy is independent of G Tot 11,C . The small drop that follows this plateau has the same approximate magnitude for all ERRs studied. For G Tot 11,C = 100 N/mm, a second plateau can be recognized directly after the descend. The plateau is less present for G Tot 11,C = 50 N/mm, and is absent for G Tot 11,C = 25 N/mm. For large incident energies, all curves converge towards a constant burst pressure, of which the value increases as G Tot 11,C does. After establishing the influence of the two governing material parameters, the identification of the pair of values that best matches the experimental results is performed. Two tests are used to calibrate the parameters: QS-1 and DYN-2 from Table 1. The remaining tests will be used later to validate the predictive capability of the method. Simulations were conducted using various values of the two parameters G Pl 11,C and G Tot 11,C and the residual burst pressures were determined subsequently. The least squares method was used to identify the set of  parameters that best matches the experimental results, the results of which are shown in Fig. 9. The identified values for the parameters are G Pl 11,C = 25 N/mm and G Tot 11,C = 80 N/mm. The total ERR value lies within the range of values in literature, i.e. approximately between 40 and 80 N/mm, [16][17][18].
The influence of the transverse ERRs is expected to be small and no tests were therefore conducted to measure these values. The used values are based on literature [19]. The set of energy release rates is summarized in Table 2.

Analysis
The performance of the model is assessed by comparing the mechanical response and initiated damage from the simulations to experimental results. Moreover, the influence of the progressive damage model is determined by comparing the results with those from the model without intralaminar damage (as used in [2]). The effect of intralaminar damage on the global behavior was expected to be small for structures with typical automotive pressure vessel dimensions, which needs to be verified. Fig. 10 shows the experimental result and both numerical results (with and without intralaminar damage) for loading of sample QS-1 from Table 1. This confirms that the influence of intralaminar damage on the global behavior is small as both simulation results closely match. Moreover, the simulation results are in close agreement with the experimental behavior.
After external contact loading, the physical pressure vessel was studied using a computed tomography scan (CT-scan) to visualize the induced damage. Fig. 11 shows an image from a CT-scan of sample QS-1 after contact loading. Damage is mainly recognized in the regions indicated by a and c, which are not directly underneath the applied load. Only minor damage (delamination) is found under the indenter, as shown for region b. The damage observed in regions a and c mainly consists of fiber kinking and fiber fracture. Fig. 12 shows the induced damage in the simulation after unloading for the same load case and at the same position with respect to the indenter. The upper subfigure shows the compressive fiber damage index D 11,C , indicating where fiber finking is predicted. It is clear that these regions exceed the regions of visible damage in the experiment. The experimentally observed shear bands can however be clearly recognized.
The lower subfigure shows the tensile fiber damage index D 11,T , hence the locations where fiber fracture is expected to occur. Fiber fracture is observed in the shear bands mainly, hence matching the experimental observations. The distinct damage that can be observed near the outer surface in the region c in the CT-image is less pronounced in the simulation.
After external contact loading and subsequent unloading, the residual burst pressure is numerically determined using the procedure previously explained in Section 2.2. The analysis results for sample QS-1 are shown in Fig. 13.
The internal pressure is increased from 0 to 250 MPa in 0.1 seconds. As a result of the increasing internal pressure and the present damage in the vessel wall, a small amount of energy is dissipated at a simulation time of between 0.03 and 0.04 s, hence between 75 and 100 MPa of internal pressure. This can be attributed to the failure of single integration points and elements, which in this case does not yet lead to failure of the surrounding layers; the load can be transferred to the surrounding layers without causing them to fail immediately. This numerical effect can be compared to the failure of single fibers and fiber tows in the experiment before the burst pressure is reached. Note that this effect is much smoother in the experiment as in the numerical model many fibers are homogenized into one integration point. Therefore, we may expect more single failure instances in the experiment, albeit with less energy per failure event. Moreover, this effect is amplified by the fact that the fiber strength may vary up to 10% in real life. The simulation does not account for this.
The succeeding rapid increase of dissipated energy at a simulation time of 0.04 s is defined as the onset of bursting. The corresponding internal pressure of 100 MPa is the numerical burst pressure for sample QS-1, as previously listed in Fig. 9. In the simulation, some elements fail and the surrounding layers cannot carry the additional load that is transferred from these failed elements. The failure of these surrounding layers causes their neighbouring layers to fail as well, with that starting an unstable propagation of element deletion throughout the thickness of  the wall, which can be recognized from the rapid increase of dissipated energy. This effect is similar to the experimental case, in which the fibers that surround failing fibers may no longer be able to carry the additional load and fail. When this starts an unstable propagation of fiber failure throughout the thickness, the vessel can locally no longer withstand the internal pressure and bursts. Note that the experimental internal pressure is released instantly upon bursting as the pressure vessel ruptures, which is not the case in the simulation. As we can see in Fig. 13, the numerical internal pressure can be increased further as the model does not account for the fact that the vessel is no longer tight. This is an important difference between the experimental and numerical burst test and explains why we need to evaluate the dissipated energy in order to determine the numerical burst pressure.

Validation
Finally, the predictive capability of the model is studied. For that purpose, simulations are conducted for the three tests that were not used for the parameter identification (i.e. DYN-1, DYN-3 and DYN-4) and the residual burst pressures are compared to the experimental results. The numerical method introduced in Section 2.2 is used to determine the residual burst pressure. Fig. 14 shows the results.
For the test at the smallest energy (DYN-1), the residual burst pressure from the experiment is slightly underestimated. However, the experimental test has a residual burst pressure that is higher than the typical results for undamaged vessels shown in the graph. It therefore has to be assumed that the vessel was undamaged and its initial burst pressure lies at the upper tail of the distribution. Since the burst pressure of the simulation lies within the typical spread of undamaged burst pressures, the predictive capability at this energy can be confirmed.
For the test at a normalized energy of about 0.5 (DYN-3), the burst pressures of the simulation and experiment closely match. The burst pressure for the dynamic test at the highest energy (DYN-4) is overestimated by the simulation. The deviation lies within the spread of initial burst pressures that may be expected for pressure vessels, as shown for the undamaged vessels. Therefore, the overestimation at this energy can be accepted and the predictive capability of the model can be confirmed.

Conclusions
A material model has been developed to predict the residual strength of a pressure vessel after an impact. The material model describes the orthotropic elasto-plastic behavior of the carbon fiber reinforced plastic composite. Damage initiation and evolution models are defined for longitudinal and transverse damage in both the tensile and compressive directions. The focus of the model is the damage initiation criterion for longitudinal compressive (fiber kinking) damage, in which the threedimensional stress state is taken into account. When damage is initiated, the damage growth is described using damage evolution laws. The   Fig. 11. CT-images after experimental quasi-static loading of sample QS-1. Image from [2].   effect of load reversal (from tensile to compressive loading or vice versa) is addressed. More specifically, the effect of fiber fracture that occurs under fiber kinking failure is key.
To analyze the influence of the material model, it was implemented in a finite element solver and used in contact loading simulations on a pressure vessel. The influence of the two governing material parameters, the energy release rates G Pl 11,C and G Tot 11,C , on the residual burst pressure of the vessel after impact was shown first. Next, the values of these parameters were determined using two experimental results. It was shown that the initiation and growth of damage adequately matches the damage observed experimentally. Finally, the capability of the model to predict the residual strength after impact was shown by comparison of numerical and experimental results at energies that were not considered in the parameter identification process. The proposed material model specifically focuses on those failure mechanisms that are known to occur when a pressure vessel is subjected to impact. Moreover, it particularly accounts for the most important mechanism that leads to the reduction of the burst pressure after an impact (fiber fracture due to fiber kinking). In this study, it is demonstrated that the material model is valid for the application for which it was developed, which is the determination of the residual burst pressure of a pressure vessel after an impact. However, to improve its general applicability, a valuable next step would be to include additional failure mechanisms that can occur in composite materials. For that purpose, the most important mechanisms from the proposed material model could be combined with those from more complete material models, such as the LS-DYNA material models MAT_261 or MAT_262.