1. Introduction
Because of its abundant supply, low pollution, and high calorific value, hydrogen energy has emerged as a significant answer to the energy problem and to achieving the “double carbon” target. Hydrogen-powered vehicles have emerged as the future vehicle trend [
1,
2,
3]. A composite-overwrapped pressure vessel is a fiber-wrapped composite construction with many layers. These tanks are typically used for high-pressure hydrogen gas storage [
4]. Composite materials used in composite-overwrapped pressure vessels are preferred over traditional materials due to their higher specific strength, stiffness, and fatigue resistance. However, the interlaminar fracture strength of composite materials is poor, and the incidence and extension of delamination significantly lower a structure’s strength and stiffness, posing a threat to the structure’s integrity and safety [
5].
The cohesion zone model (CZM) has been widely used in the simulation of delamination failure of composite structures in recent years [
6]. The main simulation methods to describe the delamination failure of composite structures are the cohesive element method and the tiebreak algorithm based on the cohesive contact method, etc.
Song, L. A., et al., used a progressive damage model based on the Puck criterion to model matrix fracture and fiber fracture (in-plane failure) in composite pressure vessels. Cohesive elements were inserted between pressure vessel layers to model interlayer failure. The simulation results were consistent with the experimental results, and the key parameters for the progressive damage analysis of composite pressure vessels were determined [
7]. Liao, B., et al., embedded bilinear cohesive elements in a composite pressure vessel model and investigated the evolution of interlaminar damage under low-velocity impact loading evolution. More severe delamination damage occurred in composite pressure vessels with aluminum lining compared to EPDM lining [
8]. Weerts et al., prelaminated a hydrogen storage tank and simulated delamination failure of the column section under quasistatic loading conditions by inserting zero-thickness cohesive elements. Interlaminar failure was usually the first visible mechanism for thick rings, whereas fiber fracture prevailed for thin rings [
9]. Aleksandr Cherniaev et al., used the tiebreak contact algorithm to model interlaminar damage behavior and analyzed the effect of nonphysical parameters of the interlaminar material model on the prediction of compressive damage loads in damaged composite cylinders. The peak stress parameter NLFS of the interlaminar model was able to accurately predict the damage loads observed in physical experiments when the parameter SLIMC, associated with the continuous damage criterion, was between 0.6 and 1 [
10]. Jinyang Zheng et al., used the tiebreak contact algorithm to simulate interlaminar failure behavior. The amplitude and phase inconsistency between adjacent layers of composite materials during vibration under explosive loading was the main cause of delamination failure [
11]. Montes De Oca Valle developed a finite element model of a specimen with less complexity than a whole cylinder. An explicit solution was used to predict the delamination damage between the composites and the separation of the metal from the composite by inserting cohesive elements. Finally, the metal–composite interface properties were obtained and calibrated through experimental tests [
12]. Drew E. Sommer et al., conducted an explicit finite element analysis of the damage response of a broken carbon fiber composite tube impacted by a falling hammer. A continuous damage mechanic (CDM) model was used to simulate interlaminar damage to the composite, and a tiebreak contact algorithm was used to simulate delamination damage to the composite. The interaction between delamination fracture energy and friction could greatly affect the results [
13]. These researchers have thoroughly investigated the external construction of composite columnar containers but have not evaluated the influence of different lying angles on the inside.
George Edward Street et al., investigated the effect of preload on large composite structures. When the loading was not parallel to the major fiber direction, the effect of preload forces could be severe. Matrix cracking and delamination could be more severe under low-velocity impact loading due to the effect of preload. The relationship between the loading main fiber direction and angle was analyzed. When the two were not parallel, the preloading had a greater impact [
14]. Naseer H. Farhood et al., predicted the first ply failure pressure based on finite element simulations of Tsai-Wu and maximum stress failure criteria. The effects of ply stacking order and orientation angle on the burst pressure were investigated with a constant ply thickness PV. However, the layering method was simple and did not separately consider the effect of different angular proportions [
15]. Zhengyun Hu et al., designed twelve stacking schemes for Type IV hydrogen storage vessels and evaluated the effect of different stacking sequences on the ultimate strength. The effect of the lamination order on the rupture pressure of the hydrogen vessel was about 15%. The laying sequence of the separated hoop and helical layers could increase the bursting pressure [
16]. Salamat-Talab et al., prepared and tested laminated composites with 0//θ (θ = 0, 15, 45, 60, and 90) interfaces, respectively, and obtained Mode II laminate fracture toughness and cohesive shear traction–separation models. The experimental results showed that the crack extension behavior in the specimens became more stable with the increase in the interfacial fiber angle [
17]. At the moment, research on the effect of varied laying orientations on composite outsourcing pressure vessels is insufficient. Few researchers have investigated the influence of a single laying method on structural strength, particularly the influence of a laying ratio of plus or minus 45 degrees on the structural strength of a thick-walled composite cylinder.
In this paper, a method of simulating the delamination of Model I and Model II composites using the tiebreak contact algorithm under different mesh sizes is verified by modifying the interlaminar strength. Additionally, under the verification of this method, a quasistatic loading simulation is carried out on the cylinder part of a thick-walled composite-overwrapped pressure vessel with a reasonable mesh size, and the location of the first significant delamination phenomenon of the thick-walled composite-overwrapped pressure vessel is verified. By increasing the proportion of spiral-wound wires (plus or minus 45 degrees) in the polar-wound layup, the effect of different plus or minus 45 degrees layup ratios on the structural strength of the compression ring is analyzed. In comparison, the fibers laid alternately at plus or minus 45 degrees can provide better structural strength in the circumferential and radial directions.
2. Problem Description
In a cohesive model, the amount of stress separation in the cohesive model artificially reduces the overall stiffness of the material due to the presence of the initial elastic phase. It is well-suited to describe the delamination phenomenon between layers in composite materials.
Although a CZM has requirements on the model mesh size, for larger simulation structures, larger mesh sizes can be used appropriately to save computation time when similar results can be achieved by correcting the interlayer strength. The tiebreak contact algorithm differs from the cohesive element method in that the surface-based cohesive behavior is a contact property rather than a material property. The tiebreak contact algorithm does not involve the usage of building finite elements, but rather exploits the nodes between the upper and lower sublayer elements to mimic interlayer bonding, as illustrated in
Figure 1 [
18]. In solving interface problems, the advantages of cohesive behavior over cohesive elements are as follows: no extra cohesive elements are required, the number of participating elements is reduced, and the results are relatively easier to converge.
To represent delamination, the two contact interfaces between the TSHELL element layers are modeled using *CONTACT_AUTOMATIC_ONE_WAY_SURFACE_TO_SURFACE_TIEBREAK_{OPTION} with OPTION = 9. OPTION = 9 is based on the fracture model in the viscous material model *MAT_COHESIVE_MIXED_MODE. The principle is shown in Equation (1):
where
is the normal stress,
is the shear stress, NFLS is the normal interlaminar strength, and SFLS is the shear interlaminar strength [
19]. In the load–displacement curve, this option is feasible for obtaining similar results with coarse meshes as with fine meshes. However, coarser meshes require lower peak traction forces [
20]. As the load keeps increasing, the stress σ between the layers of the specimen continues to increase; when the stress σ reaches the failure criterion, the contact point starts to break down, and the interface starts to separate. After that, as the interface separation distance δ decreases linearly, the interface separates further until it reaches the interface limit separation distance, the stress σ decreases to 0, the contact point is completely destroyed, and the interface separates completely, marking the beginning of crack generation. With the complete destruction of contact between multiple groups of nodes, the crack starts to expand. The energy release rate at the interface separation is, theoretically, equal to the area of the region below the σ–δ curve, as shown in
Figure 2.
In order to obtain accurate calculation results, the model needs to be divided into a fine grid, with the length of the layered process area to be divided into a certain number of elements to meet the requirements of Equation (2):
where
is the length of the cohesive zone,
is the length of a single element on the crack extension surface, and
is the number of elements to be divided. Turon et al., showed that it was not less than 3. That is, the size of the elements in the cracking direction was not larger than 0.5 mm to ensure the accuracy of the calculation [
21]. This severely limits the application of the tiebreak algorithm in large-scale simulations. Rice et al., proposed a series of cohesive models where the definition of the length of the stratified process area could all be expressed uniformly in Equation (3) [
22]:
where
is the Young’s modulus of elasticity of the material;
is the critical energy release rate;
is the interlaminar strength; and
M is a constant less than 1, where the value of M has not been determined. The most commonly used values are Rice [
22] and Hillerborg [
23], which are considered to be 0.88 or 1. Pedro PC et al.’s, study showed that the accuracy of the simulation results was mainly affected by the energy release rate
, and the correction of the interlaminar strength could be ensured by the constant energy release rate
(as shown in
Figure 2). The simulation accuracy could be guaranteed when the two areas were equal. It can be seen from Equation (2) that, when the interlaminar strength value is reduced, the length of the cohesive zone becomes larger. Under the condition that the number of grids in the cohesion zone is sufficient, a larger grid size can be used to quickly and accurately calculate the damage and failure load of the adhesive layer, thereby improving the calculation efficiency. This provides the possibility to simulate composite delamination using the tiebreak algorithm under meshes of different sizes.
3. Simulation of Model I and Model II Delamination Failure of Composite Materials
In order to model delamination failure appropriately, two delamination modes were considered to be relevant within the vessel model: BCD (opening delamination) and ENF (shear forward delamination) [
24]. The tiebreak contact algorithm determined the failure equation of the interface from two main types of interlaminar fracture modes, model I and model II, and could accurately simulate the crack generation and expansion process. Since the test methods for pure model I and model II fracture toughness have become mature, the fracture toughness (
,
) in the model I and model II modes is usually selected as the basic parameter to describe the failure response of mixed-mode composite structures, as shown in
Figure 3.
The specimen material parameters were selected based on the literature [
25]. The specimen length was 100 mm, the width was 20.0 mm, the height was 2 h, each h thickness was 1.55 mm, and the initial delamination length was 35 mm. In the finite element model calculation, the model elements model was selected as TSHELL elements for the model I and model II delamination failure simulations, with element lengths of 0.4 mm, 0.5 mm, 1.0 mm, 1.5 mm, and 2.0 mm. The total numbers of elements were 100, 400, 64,000, 8160, 1742, and 1020, as shown in
Figure 4.
3.1. DCB Delamination Failure Simulations
In the model Ⅰ DCB test, the energy release rate versus time is shown in
Figure 5, which shows that
played a major role, and the effect of
on delamination was negligible. The calculation results of the displacement of the load with the load application point for different element sizes are given in
Figure 6.
The calculation results for element lengths of 0.4 mm and 0.5 mm were overlapped. An element length of less than 0.5 mm did not affect the calculation accuracy. The calculation results of the element length of 1.00 mm was very close to the calculation results of 0.4 mm and 0.5 mm.
Figure 6a shows the uncorrected DCB load–displacement curves based on the tiebreak contact algorithm, from which it can be seen that the slopes of the load–displacement curves remained in good agreement for different grid sizes for the model I DCB test. When
and
were kept constant, the greater the mesh size, the greater the stress required for interface separation and the more noticeable the ensuing oscillation in delamination damage.
Figure 6b depicts the load–displacement curve based on the tiebreak contact correction.
It can be seen that the different mesh densities that could be obtained after the conversion of interlayer strength could be very approximate. However, with increases in mesh size and length, the oscillation in delamination damage was still obvious. For an element length of 1.0 mm, the load–displacement curve based on the tiebreak contact algorithm agreed well with the load–displacement curves of 0.4 mm and 0.5 mm. The load–displacement curve for the element length of 1.0 mm ensured the accuracy of the calculation and had better calculation efficiency compared with the element lengths of 0.4 mm and 0.5 mm,
Figure 7.
3.2. ENF Delamination Failure Simulations
The energy release rate versus time in the model II ENF test is shown in
Figure 8, which shows that
played a major role, and the effect of
on delamination was negligible. The results of load and displacement calculations for the loading of members with different element lengths are given in
Figure 9. Element lengths less than 0.5 mm did not affect the calculation accuracy. The calculation results for the element length of 1.00 mm were also very close to the calculation results of 0.4 mm and 0.5 mm.
Figure 9a shows the uncorrected ENF load–displacement curve based on the tiebreak contact algorithm, from which it can be seen that the slope of the load–displacement curve also maintained a good agreement for different grid sizes for the model II ENF test. When
and
were kept constant, the larger the mesh size, the larger the load required for interface separation. Compared with the model I DCB test, the oscillation in delamination damage by increasing mesh size was not obvious. The load–displacement curves based on the tiebreak contact correction in
Figure 9b was obtained through relatively reducing interlaminar strength. It can be seen that, for the model II ENF test, the results could still be obtained very approximately after the conversion of interlayer strength for different grid densities. The comparison between the two figures shows that the load–displacement curve based on the tiebreak contact correction with an element length of 1.0 mm still had a good agreement with the calculation results of 0.5 mm while improving the calculation efficiency.
In both the model I DCB test and the model II ENF test, the error increased with the increase in element length. However, the results for the element length of 1.0 mm showed good similarity with element lengths of 0.5 mm and 0.4 mm. This is in agreement with the findings of Camanho et al., in the literature [
26].
5. Conclusions
(1) In order to model delamination failure appropriately, a method of simulating the delamination of Model I and Model II composites using the tiebreak contact algorithm under different mesh sizes was verified by modifying the interlaminar strength. If the interlaminar strength σ remained constant as the displacement δ at delamination failure increased, the area of the energy release rate for delamination damage increased, and the load required for delamination damage increased with the corresponding increase. The oscillation caused by the delamination damage became more visible as well. If the interlaminar strength σ was reduced at this time, the area of the delamination damage’s energy release rate remained constant. The interlaminar strengths of various meshes were converted to produce very close results. Because of the fewer elements, the model with element lengths of 1.0 mm had higher computational efficiency than the models with element lengths of 0.4 mm and 0.5 mm, and the computational results were still in good agreement with those of the models with element lengths of 0.4 mm and 0.5 mm.
(2) When the thick-walled composite-overwrapped pressure vessels were quasistatically loaded, there was almost no in-plane failure, and the fibers were in the elastic–plastic stage. The elastic–plastic process was verified by selecting a progressive damage model with quasistatic nonlinear tensile shear of the material. The slope of the simulated numerical section was basically similar to that of the experiment, and the stress increased uniformly with strain, which could well-reflect the elastic phase of the tensile process. In the plastic phase, the stress could be kept constant at about 720 MPa by setting the stiffness retention factor. The progressive damage model could respond to the elastic–plastic phase of nonlinear shear of composite materials ideally.
(3) A 25 mm composite-overwrapped pressure vessel cylindrical model was established, and a quasistatic loading simulation was carried out. The first delamination failure occurred approximately in the center. In comparison, the fibers laid alternately at ±45° could provide better structural strength in the circumferential and radial directions. For axial strength, the [(45°)5 /0°10]5 combined layup model was higher than single polar winding [0°]100 and spiral winding [45°]50s. The ±45° layup had a direct effect on the flexural properties of the composite. The proportion of fibers at 50% had a higher modulus value than a pure 0° ply, and the structure could be stronger than a 0° ply in terms of load capacity after adding ±45° ply, but too large a ratio was detrimental to the improvement of structural properties.
This simulation study analyzed the difference in structural strength between the combined plus and minus 45-degree laying method and single laying methods, and specific ratio optimization will be the focus of future research.