Method for predicting residual strength of composite laminated plate containing delamination damage by using finite element method
Through the finite element method and the three-dimensional Hashin failure criteria, the problem that the two-dimensional Hashin criterion in ABAQUS software is solved, and the accurate prediction of the residual strength of composite laminates and the detailed description of the material characteristics is achieved.
Patent Information
- Application Number
- CN202510291307.X
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- Filing Date
- 2025-03-12
- Publication Date
- 2025-06-27
AI Technical Summary
The two-dimensional Hashin criterion in the existing ABAQUS software is difficult to fully describe complex failure modes such as three-dimensional stress state and stratified damage, and it is impossible to effectively predict the residual strength of composite laminates.
Using the finite element method, the composite laminate model is constructed in ABAQUS, the three-dimensional Hashin failure criteria are set, the expansion process of layered damage is simulated, and the grid division and boundary conditions are performed are carried out to calculate the residual strength of the composite laminate.
It can accurately describe the complex stress states and multiple failure modes of composite laminates in three-dimensional space, simulate damage states that cannot be predicted by two-dimensional models, evaluate the residual strength of composite laminates, and describe the anisotropy and damage characteristics of the material.
Smart Images

Figure CN120217773A_ABST
Abstract
Description
Technical Field
[0001] The present invention relates to a method for predicting the remaining strength of a composite laminate. Background Art
[0002] Carbon fiber composite materials have the characteristics of high strength, light weight, corrosion resistance, and excellent insulation performance, and are easy to be processed into various complex shapes to meet diverse design requirements. They are widely used in many fields such as aerospace, automotive manufacturing, building structures, sports goods, and energy, infrastructure, ocean, military, etc., providing high-performance, durable, and innovative solutions for various industries. The development of modern industry has also continuously advanced the design and manufacturing technologies of composite structures, and the usage amount of composite materials shows an increasing trend year by year.
[0003] A composite laminate is composed of unidirectional composite laminates, and the laminate design is used to achieve the optimization and complementarity of performance. However, if the composite laminate structure is not properly operated during use, or is subjected to factors such as low-energy impact during assembly, delamination damage is easily induced. Such damage is often located inside the laminate, resulting in a decrease in the stiffness of the component. Under the action of external loads, internal damage is easily further expanded, leading to phenomena such as laminate buckling. These phenomena often cause the structure to fail prematurely before reaching the failure strength of the material itself, seriously weakening the stability and safety of the structure.
[0004] Although the two-dimensional Hashin criterion in the ABAQUS software considers the anisotropy of materials to a certain extent, it has limitations in predicting the remaining strength of composite laminates with delamination damage and is difficult to comprehensively describe complex failure modes such as three-dimensional stress states and delamination damage. Summary of the Invention
[0005] In order to solve the problem that the two-dimensional Hashin criterion in the existing ABAQUS software is difficult to comprehensively describe complex failure modes such as three-dimensional stress states and delamination damage, the present invention proposes a method for predicting the remaining strength of a composite laminate with delamination damage by using the finite element method.
[0006] The method for predicting the remaining strength of a composite laminate with delamination damage by using the finite element method according to the present invention is carried out according to the following steps:
[0007] 1. Determine the location of the delamination damage inside the composite laminate test piece. The upper sub-plate is above the location of the delamination damage, and the lower sub-plate is below the location of the delamination damage. Determine the geometric dimensions of the upper sub-plate and the lower sub-plate, and construct the upper sub-plate model and the lower sub-plate model in abaqus; then divide the plies of the upper sub-plate model and the lower sub-plate model.
[0008] II. Construct a cohesive layer model in Abaqus, and cut holes in the cohesive layer model according to the shape and spatial position of the delamination damage in the composite laminate test piece to simulate the damage in the composite laminate test piece;
[0009] III. In Abaqus, make the upper sub-plate model, the lower sub-plate model and the cohesive layer model contact each other to form a cuboid as a whole, and obtain the composite laminate model;
[0010] IV. Establish a composite laminate failure criterion in Abaqus; then set the constitutive relationship of the cohesive layer model, set the performance parameters of the materials of the upper sub-plate model and the lower sub-plate model, and set the performance parameters of the cohesive layer model; judge the initiation of damage in the cohesive layer model and define the delamination damage propagation process;
[0011] V. Set the friction contact between the upper sub-plate model and the lower sub-plate model in the area where the delamination damage of the composite laminate model exists; use hexahedral elements to mesh the upper sub-plate model and the lower sub-plate model, and mesh the cohesive layer model;
[0012] VI. Determine the boundary conditions and apply loads, calculate the composite laminate model, obtain the calculation results after calculation and extract the load-displacement curve to obtain the residual strength of the composite laminate.
[0013] The beneficial effects of the present invention are as follows:
[0014] The present invention uses the three-dimensional Hashin failure criterion to describe the complex stress state and various failure modes of the composite laminate in three-dimensional space, can simulate damage states that cannot be predicted by two-dimensional models such as out-of-plane shear damage, can be used for laminated plates with delamination damage of any shape, quantity, size and internal thickness position, can be used to evaluate the residual strength of composite laminates containing internal delamination, and can describe in detail the anisotropy of materials and the damage mechanical behavior and damage characteristics, which is convenient for designing any reasonable evaluation of the structural strength of the laminate, and provides a reference method for the delamination propagation simulation and residual strength evaluation of damaged composite parts. Description of the Drawings
[0015] Figure 1 The cohesive layer model constructed in step II of Example 1;
[0016] Figure 2 The composite laminate model constructed in step III of Example 1;
[0017] Figure 3 The damage schematic diagram of the cohesive layer model after calculating the composite laminate model in step VI of Example 1;
[0018] Figure 4 It is a damage schematic diagram of the lower sub - plate model after calculating the composite laminate model in Step Six of Embodiment 1;
[0019] Figure 5 It is the load - displacement curve extracted in Step Six of Embodiment 1. Specific Embodiments
[0020] The technical solution of the present invention is not limited to the specific embodiments listed below, but also includes any reasonable combination between specific embodiments.
[0021] Specific Embodiment One: The method for predicting the residual strength of a composite laminate with delamination damage using the finite element method in this embodiment is carried out according to the following steps:
[0022] 1. Determine the location of the delamination damage inside the composite laminate test piece. The part above the delamination damage location is the upper sub - plate, and the part below the delamination damage location is the lower sub - plate. Determine the geometric dimensions of the upper sub - plate and the lower sub - plate, and construct the upper sub - plate model and the lower sub - plate model in abaqus; then divide the plies of the upper sub - plate model and the lower sub - plate model.
[0023] 2. Construct a cohesive layer model in abaqus, and cut holes in the cohesive layer model according to the shape and spatial position of the delamination damage in the composite laminate test piece to simulate the damage in the composite laminate test piece.
[0024] 3. In abaqus, make the upper sub - plate model, the lower sub - plate model and the cohesive layer model contact each other to form a cuboid as a whole, and obtain the composite laminate model.
[0025] 4. Establish a composite laminate failure criterion in abaqus; then set the constitutive relationship of the cohesive layer model, set the performance parameters of the materials of the upper sub - plate model and the lower sub - plate model, and set the performance parameters of the cohesive layer model; judge the initiation of damage in the cohesive layer model and define the delamination damage propagation process.
[0026] 5. Set the contact between the upper sub - plate model and the lower sub - plate model as frictional contact in the area where the delamination damage exists in the composite laminate model; use hexahedral elements to mesh the upper sub - plate model and the lower sub - plate model, and mesh the cohesive layer model.
[0027] 6. Determine the boundary conditions and apply loads, calculate the composite laminate model, obtain the calculation results after calculation and extract the load - displacement curve to obtain the residual strength of the composite laminate.
[0028] The beneficial effects of this embodiment are:
[0029] In this embodiment, the three-dimensional Hashin failure criterion is adopted to describe the complex stress state and various failure modes of the composite laminate in three-dimensional space. It can simulate the damage states that cannot be predicted by two-dimensional models such as out-of-plane shear damage. It can be used for laminated plates with delamination damage of any shape, quantity, size, and internal thickness position. It can be used to evaluate the remaining strength of composite laminates containing internal delaminations, and to describe in detail the anisotropy of materials, as well as the damage mechanical behavior and damage characteristics, which is convenient for designing any reasonable method to evaluate the structural strength of laminates, providing a reference method for simulating the delamination propagation and evaluating the remaining strength of damaged composite parts.
[0030] Specific Embodiment 2: The difference between this embodiment and Specific Embodiment 1 is that: the composite laminate test piece described in Step 1 is composed of multiple unidirectional composite layers.
[0031] Specific Embodiment 3: The difference between this embodiment and Specific Embodiment 1 or 2 is that: the ply thickness and ply quantity described in Step 1 are the same as those of the unidirectional composite layers inside the composite laminate test piece.
[0032] Specific Embodiment 4: The difference between this embodiment and any one of Specific Embodiments 1 to 3 is that: both the upper sub-plate model and the lower sub-plate model described in Step 1 are cuboids.
[0033] Specific Embodiment 5: The difference between this embodiment and any one of Specific Embodiments 1 to 4 is that: the thicknesses of both the upper sub-plate model and the lower sub-plate model described in Step 1 are the same as those of the upper sub-plate and the lower sub-plate in the composite laminate.
[0034] Specific Embodiment 6: The difference between this embodiment and any one of Specific Embodiments 1 to 5 is that: the cohesive layer model described in Step 2 is a cuboid.
[0035] Specific Embodiment 7: The difference between this embodiment and any one of Specific Embodiments 1 to 6 is that: the thickness of the cohesive layer model in Step 2 is 5% of the thickness of the unidirectional composite layer in the composite laminate test piece.
[0036] Specific Embodiment 8: The difference between this embodiment and any one of Specific Embodiments 1 to 7 is that: in Step 5, the upper sub-plate model and the lower sub-plate model are meshed using C3D8R elements.
[0037] Specific Embodiment 9: The difference between this embodiment and any one of Specific Embodiments 1 to 8 is that: in Step 5, the cohesive layer model is meshed using COH3D8.
[0038] Specific Embodiment 10: The difference between this embodiment and any one of Specific Embodiments 1 to 9 is that: in Step 4, the three-dimensional Hashin failure criterion formula based on stress for the composite laminate failure criterion is established in abaqus;
[0039] σ 11 When it is greater than 0, it is the fiber tensile mode, and the three-dimensional Hashin failure criterion formula is:
[0040]
[0041] σ 11 When it is less than or equal to 0, it is the fiber compression mode, and the three-dimensional Hashin failure criterion formula is:
[0042]
[0043] σ 22 +σ 33 When it is greater than 0, it is the matrix tensile mode, and the three-dimensional Hashin failure criterion formula is:
[0044]
[0045] σ 22 +σ 33 When it is less than 0, it is the matrix compression mode, and the three-dimensional Hashin failure criterion formula is:
[0046]
[0047] The three-dimensional Hashin failure criterion formula for the in-plane shear mode is:
[0048]
[0049] The three-dimensional Hashin failure criterion formula for the out-of-plane shear mode is:
[0050]
[0051] In Formulas 1 to 6: σ 11 is the principal stress in the x direction of the unidirectional composite layer; σ 22 is the principal stress in the y direction of the unidirectional composite layer; σ 33 is the principal stress in the z direction of the unidirectional composite layer; σ 12 is the shear stress in the xy plane of the unidirectional composite layer; σ 13 is the shear stress in the xz plane of the unidirectional composite layer; σ 23 is the shear stress in the yz plane of the unidirectional composite layer; X T is the longitudinal tensile strength of the unidirectional composite layer, X C is the longitudinal compressive strength of the unidirectional composite layer; S 12 is the shear strength in the xy plane of the unidirectional composite layer; S 13 is the shear strength in the xz plane of the unidirectional composite layer; S 23 is the shear strength in the yz plane of the unidirectional composite layer;
[0052] In Step 4, the constitutive relation of the cohesive layer model is set as follows:
[0053] τ = (1 - d)Kλ Formula 7
[0054] In Formula 7: τ is the normal cohesive force of the cohesive layer model, K is the initial penalty stiffness of the cohesive layer model, d is the damage coefficient related to the displacement of the cohesive layer model, and λ represents the normal relative displacement of the cohesive layer model;
[0055] In Step 4, the quadratic nominal stress criterion is adopted to judge the initiation of damage in the cohesive layer model:
[0056]
[0057] In Formula 8: σ z is the normal tensile stress of the cohesive layer model perpendicular to the xy plane; τ xz is the shear stress in the xz direction of the cohesive layer model, and τ yz is the shear stress in the yz direction of the cohesive layer model; T is the interlayer tensile strength of the cohesive layer model; S is the interlayer shear strength of the cohesive layer model;
[0058] In Step 4, the power exponent criterion is adopted to define the process of delamination damage propagation:
[0059]
[0060] In Formula 9: α, β, and γ are power factors respectively, with values of 2; G I is the energy release rate of type I crack propagation; G II is the energy release rate of type II crack propagation; G III is the energy release rate of type III crack propagation; G IC is the fracture toughness of type I; G IIC is the fracture toughness of type II; G IIIC is the fracture toughness of type III.
[0061] Example 1
[0062] The method for predicting the residual strength of a composite laminate containing delamination damage using the finite element method in this example is carried out according to the following steps:
[0063] 1. Determine the location of delamination damage inside the composite laminate test piece. The upper sub - plate is above the delamination damage location, and the lower sub - plate is below the delamination damage location. Determine the geometric dimensions of the upper sub - plate and the lower sub - plate, and construct the upper sub - plate model and the lower sub - plate model in abaqus; then divide the plies for the upper sub - plate model and the lower sub - plate model;
[0064] The composite laminate test specimen is composed of multiple unidirectional composite layers; the size of the composite laminate test specimen is 150×100×3.2 mm, the thickness of the unidirectional composite layer is 0.2 mm, the unidirectional composite layer is a T700 carbon fiber reinforced epoxy composite material, with a total of 16 layers, and the delamination damage is located at 1 / 4 of the laminate. The upper sub - plate contains 4 unidirectional composite layers, and the lower sub - plate contains 12 unidirectional composite layers;
[0065] The ply thickness and ply number are the same as those of the unidirectional composite layers inside the composite laminate test specimen;
[0066] Both the upper sub - plate model and the lower sub - plate model are cuboids;
[0067] The thicknesses of both the upper sub - plate model and the lower sub - plate model are the same as those of the upper sub - plate and the lower sub - plate in the composite laminate;
[0068] Second, construct a cohesive layer model in abaqus, and cut holes in the cohesive layer model according to the shape and spatial position of the delamination damage in the composite laminate test specimen to simulate the damage in the composite laminate test specimen; the hole is an ellipse with a major axis of 40 mm and a minor axis of 20 mm; Figure 1 The cohesive layer model constructed for Step 2;
[0069] The cohesive layer model is a cuboid, and the thickness of the cohesive layer model is 5% of the thickness of the unidirectional composite layer in the composite laminate test specimen;
[0070] Third, in abaqus, make the upper sub - plate model, the lower sub - plate model and the cohesive layer model contact each other to form a cuboid as a whole, and obtain the composite laminate model; Figure 2 The composite laminate model constructed for Step 3;
[0071] Fourth, establish a composite laminate failure criterion in abaqus; then set the constitutive relationship of the cohesive layer model, set the performance parameters of the materials of the upper sub - plate model and the lower sub - plate model, and set the performance parameters of the cohesive layer model; judge the initiation of damage in the cohesive layer model and define the process of delamination damage propagation; the performance parameters of the materials are shown in Table 1; the performance parameters of the cohesive layer model are shown in Table 2;
[0072] Table 1
[0073]
[0074] In Table 1, E1, E2, and E3 are the elastic moduli in the x, y, and z directions of the unidirectional composite layer respectively, G12, G13, G 23 are the shear moduli in the xy, xz, and yz planes of the unidirectional composite layer respectively, ν12 and ν 13 and ν 23 are the Poisson's ratios between the x - direction and the y - direction, the x - direction and the z - direction, and the y - direction and the z - direction of the unidirectional composite material layer, respectively.
[0075] Table 2
[0076]
[0077] V. Set the contact between the upper sub - plate model and the lower sub - plate model as frictional contact in the area where delamination damage exists in the composite laminate model; use hexahedral elements to mesh the upper sub - plate model and the lower sub - plate model, and mesh the cohesive layer model;
[0078] In step V, use C3D8R elements to mesh the upper sub - plate model and the lower sub - plate model;
[0079] In step V, use COH3D8 to mesh the cohesive layer model;
[0080] VI. Determine the boundary conditions and apply loads, calculate the composite laminate model, obtain the calculation results after calculation and extract the load - displacement curve to obtain the residual strength of the composite laminate. In this embodiment, refer to GB / T 21239 - 2022 "Test Method for Compressive Properties of Fiber - Reinforced Plastics Laminates After Impact", and determine the boundary conditions according to the load application method of the composite test piece; in this example, the applied load is a compressive load, and the applied boundary conditions are fixed at both sides and simply supported at the other two sides.
[0081] In step IV, establish the stress - based three - dimensional Hashin failure criterion formula for the composite laminate in abaqus;
[0082] σ 11 When σ > 0, it is the fiber tensile mode, and the three - dimensional Hashin failure criterion formula is:
[0083]
[0084] σ 11 When σ ≤ 0, it is the fiber compression mode, and the three - dimensional Hashin failure criterion formula is:
[0085]
[0086] σ 22 +σ 33 When σ + σ > 0, it is the matrix tensile mode, and the three - dimensional Hashin failure criterion formula is:
[0087]
[0088] σ 22+σ 33 When it is less than 0, it is the matrix compression mode, and the three-dimensional Hashin failure criterion formula is:
[0089]
[0090] The three-dimensional Hashin failure criterion formula for the in-plane shear mode is:
[0091]
[0092] The three-dimensional Hashin failure criterion formula for the out-of-plane shear mode is:
[0093]
[0094] In Formulas 1 to 6: σ 11 is the principal stress in the x-direction of the unidirectional composite layer; σ 22 is the principal stress in the y-direction of the unidirectional composite layer; σ 33 is the principal stress in the z-direction of the unidirectional composite layer; σ 12 is the shear stress in the xy plane of the unidirectional composite layer; σ 13 is the shear stress in the xz plane of the unidirectional composite layer; σ 23 is the shear stress in the yz plane of the unidirectional composite layer; X T is the longitudinal tensile strength of the unidirectional composite layer, X C is the longitudinal compressive strength of the unidirectional composite layer; S 12 is the shear strength of the xy plane of the unidirectional composite layer; S 13 is the shear strength of the xz plane of the unidirectional composite layer; S 23 is the shear strength of the yz plane of the unidirectional composite layer;
[0095] In Step 4, the constitutive relationship of the cohesive layer model is set as:
[0096] τ=(1 - d)Kλ Formula 7
[0097] In Formula 7: τ is the normal cohesive force of the cohesive layer model, K is the initial penalty stiffness of the cohesive layer model, d is the damage coefficient related to the displacement of the cohesive layer model, and λ represents the normal relative displacement of the cohesive layer model;
[0098] In Step 4, the quadratic nominal stress criterion is used to judge the initiation of damage in the cohesive layer model:
[0099]
[0100] In Formula 8: σ z is the normal tensile stress perpendicular to the xy plane of the cohesive layer model; τ xz is the shear stress in the xz direction of the cohesive layer model, τyz is the shear stress in the yz direction of the cohesive layer model; T is the interlaminar tensile strength of the cohesive layer model; S is the interlaminar shear strength of the cohesive layer model;
[0101] In step 4, the power-law criterion is adopted to define the delamination damage propagation process:
[0102]
[0103] In Equation 9: α, β, and γ are power factors, with a value of 2; G I is the energy release rate for mode I crack propagation; G II is the energy release rate for mode II crack propagation; G III is the energy release rate for mode III crack propagation; G IC is the fracture toughness of mode I; G IIC is the fracture toughness of mode II; G IIIC is the fracture toughness of mode III.
[0104] Figure 3 is the damage schematic diagram of the cohesive layer model after the calculation of the composite laminate model; Figure 4 is the damage schematic diagram of the lower subplate model after the calculation of the composite laminate model; The finite element model is simulated and analyzed. Figure 5 is the extracted load-displacement curve, and the residual strength of the composite laminate is obtained as 90.447 kN.
Claims
1. A method for predicting the residual strength of a composite laminate with delamination damage using a finite element method, characterized in that: The method for predicting the residual strength of composite laminates with delamination damage using the finite element method is carried out in the following steps:
1. Determine the location of the delamination damage inside the composite laminate test piece. The upper sub-plate is above the location of the delamination damage, and the lower sub-plate is below the location of the delamination damage. Determine the geometric dimensions of the upper sub-plate and the lower sub-plate, and build the upper sub-plate model and the lower sub-plate model in abaqus; then divide the upper sub-plate model and the lower sub-plate model into layers; Second, a cohesive layer model is constructed in Abaqus. According to the shape and spatial position of the delamination damage in the composite laminate test piece, holes are cut in the cohesive layer model to simulate the damage in the composite laminate test piece.
3. In abaqus, the upper sub-plate model, the lower sub-plate model and the cohesive layer model are contacted with each other to form a rectangular whole, so as to obtain a composite laminate model; Fourth, establish the failure criterion of composite laminates in abaqus; then set the constitutive relationship of the cohesive layer model, set the material performance parameters of the upper sub-plate model and the lower sub-plate model, and set the performance parameters of the cohesive layer model; determine the damage initiation of the cohesive layer model and define the delamination damage extension process; 5. In the area where delamination damage exists in the composite laminate model, the upper sub-plate model and the lower sub-plate model are set to be in friction contact; the upper sub-plate model and the lower sub-plate model are meshed using hexahedral units, and the cohesive layer model is meshed; 6. Determine the boundary conditions and apply loads, calculate the composite laminate model, obtain the calculation results and extract the load-displacement curve after calculation, and obtain the residual strength of the composite laminate.
2. The method for predicting the residual strength of a composite laminate containing delamination damage by using the finite element method according to claim 1, characterized in that: Step 1: The composite laminate test piece is composed of multiple unidirectional composite material layers.
3. The method for predicting the residual strength of a composite laminate containing delamination damage by using the finite element method according to claim 1, characterized in that: The ply thickness and ply quantity in step 1 are the same as those of the unidirectional composite material layer inside the composite laminate test piece.
4. The method for predicting the residual strength of a composite laminate containing delamination damage by using the finite element method according to claim 1, characterized in that: The upper sub-plate model and the lower sub-plate model described in step 1 are both rectangular parallelepipeds.
5. The method for predicting the residual strength of a composite laminate containing delamination damage by using the finite element method according to claim 1, characterized in that: The thickness of the upper sub-plate model and the lower sub-plate model in step 1 are the same as the upper sub-plate and the lower sub-plate in the composite laminate.
6. The method for predicting the residual strength of a composite laminate containing delamination damage by using the finite element method according to claim 1, characterized in that: The cohesive layer model in step 2 is a cuboid.
7. The method for predicting the residual strength of a composite laminate containing delamination damage by using the finite element method according to claim 1, characterized in that: Step 2: The thickness of the cohesive layer model is 5% of the thickness of the unidirectional composite material layer in the composite laminate test piece.
8. The method for predicting the residual strength of composite laminates containing delamination damage by using the finite element method according to claim 1, characterized in that: In step five, the upper sub-plate model and the lower sub-plate model of C3D8R unit are used for meshing.
9. The method for predicting the residual strength of composite laminates containing delamination damage by using the finite element method according to claim 1, characterized in that: In step five, COH3D8 is used to mesh the cohesive layer model.
10. The method for predicting the residual strength of composite laminates containing delamination damage by using the finite element method according to claim 1, characterized in that: In step 4, a stress-based three-dimensional Hashin failure criterion formula for composite laminate failure criterion is established in abaqus; σ 11 When >0, it is the fiber tensile mode, and the three-dimensional Hashin failure criterion formula is: σ 11 When ≤0, it is the fiber compression mode, and the three-dimensional Hashin failure criterion formula is: σ 22 +σ 33 When >0, it is the matrix tensile mode, and the three-dimensional Hashin failure criterion formula is: σ 22 +σ 33 When <0, it is the matrix compression mode, and the three-dimensional Hashin failure criterion formula is: The three-dimensional Hashin failure criterion formula for the in-plane shear mode is: The three-dimensional Hashin failure criterion formula for the out-of-plane shear mode is: In Formula 1 to Formula 6: σ 11 is the principal stress in the x direction of the unidirectional composite material layer; σ 22 is the principal stress in the y direction of the unidirectional composite material layer; σ 33 is the z-direction principal stress of the unidirectional composite material layer; σ 12 is the shear stress in the xy plane of the unidirectional composite material layer; σ 13 is the shear stress in the xz plane of the unidirectional composite material layer; σ 23 is the yz plane shear stress of the unidirectional composite material layer; X T is the longitudinal tensile strength of the unidirectional composite material layer, X C is the longitudinal compressive strength of the unidirectional composite material layer; S 12 is the shear strength of the unidirectional composite layer in the xy plane; S 13 is the shear strength of the unidirectional composite layer in the xz plane; S 23 is the shear strength of the unidirectional composite layer in the yz plane; In step 4, the constitutive relation of the cohesive layer model is set as: τ=(1-d)Kλ Formula 7 In formula 7: τ is the normal cohesion of the cohesive layer model, K is the initial penalty stiffness of the cohesive layer model, d is the damage coefficient related to the displacement of the cohesive layer model, and λ represents the normal relative displacement of the cohesive layer model; In step 4, the quadratic nominal stress criterion is used to determine the damage initiation of the cohesive layer model: In formula 8: σ z is the normal tensile stress of the cohesive layer model perpendicular to the xy plane; τ xz is the shear stress in the xz direction of the cohesive layer model, τ yz is the shear stress in the yz direction of the cohesive layer model; T is the interlaminar tensile strength of the cohesive layer model; S is the interlaminar shear strength of the cohesive layer model; In step 4, the power exponential criterion is used to define the layered damage propagation process: In formula 9, α, β, and γ are power factors, each with a value of 2; G I is the energy release rate of mode I crack growth; G II G is the energy release rate of mode II crack growth; III G is the energy release rate of mode III crack growth; IC is the type I fracture toughness; G IIC is the type II fracture toughness; G IIIC It is the type III fracture toughness.