Progressive Damage Analysis (PDA) of Carbon Fiber Plates with Out-of-Plane
Fold under Pressure

The out-of-plane fold is a common defect of composite materials during the manufacturing process and will greatly affect the compressive strength as well as the service life. Making it of great importance to investigate the influence of out-of-plane defects to the compressive strength of laminate plates of composite materials, and to understand the patterns of defect evolution. Therefore, the strip method is applied in this article to create out-of-plane defects with different aspect ratios in laminated plates of composite materials, and a compressive performance test is conducted to quantify the influence of out-of-plane defects. The result shows that the compressive strength becomes weaker with a greater aspect ratio. When the highest aspect ratio is set to 0.12 in the experiment, the compressive strength reduces by 36.1%. Then we establish a 3-D progressive damage model based on continuum mechanics, and write it into the UMAT subroutine together with the 3-D Hashin criteria and the non-linear degradation criteria of materials. 3-D solid modeling is performed for the samples with an out-of-plane fold based on ABAQUS, and progressive damage analysis is conducted to acquire the inplane evolution process of initial failure strength with different laminates. The experimental results agree well with the simulation results.


Introduction
The carbon-fiber reinforced resin compound is widely applied in many fields of engineering due to its great mechanical performance including high specific strength and high specific stiffness, corrosion resistance and designability, such as automobile industry, aerospace industry, marine ships industry and medical field [1][2][3]. Tawfi analyzed the free vibration of composite laminates by using Neural Network-Based Second Order Reliability Method [4]. Hong studied the long-term durability of composite laminates [5]. However, the process of forming composite materials is very complicated. Especially in the process of manufacturing composite components of large dimensions and great thickness, many defects may frequently occur, such as corrugation, pore, fold, dry spot, and delamination, which may greatly decrease relevant mechanical performance and usage of composite materials [6]. Especially out-of-plane fold defects in the forming process could greatly reduce the compressive strength of laminated plates of composite materials and the service life of composite components. Therefore, it is currently a hot topic to effectively investigate the influence of out-of-plane folds, and to evaluate the strength of composite laminated plates with out-of-plane folds. Domestic and international scholars have made relevant researches.
The current study of mechanical performance can be divided into two aspects in general. On the one aspect, experiments are conducted to quantify the influence of an out-of-plane fold to the strength of composite laminate plates. The relevant representative researches include: Wu fabricated laminated plates with out-of-plane folds of different aspect ratios by applying high temperature curing pressurization. Then conducting the mechanical performance test, to research the affect of out-of-plane defects on the mechanical performance of composite laminate plates [7]. Allison and Evans fabricated fiberglass laminated plates with local laminate corrugation. And conducted a series of three-point bending tests to describe how the fiber corrugation influences the strength of fiberglass composite materials [8]. Adams fabricated T300/P1700/carbon/polysulfone composite laminate plates with fiber corrugation, and fabricated laminated plates with isolated corrugation by a three-step method, having studied the bearing capability of laminated plates of different waviness, with results showing that the strength of composite laminate plates would decrease more with more adopted waves [9]. Bradley produced out-of-plane folds in composite laminated plates with a secondary forming method, and conducted compressive performance test of the samples, with results that the compressive strength would decrease more when adopting 3 layers of waves than 1 layer [10]. Dodwell established a 1-D model of out-of-plane folds in the external radius curing process, with predicted fold wavelength and critical fold conditions closely matching the fold defects observed in the crystal display [11]. On the other aspect, the failure mechanism and damage behaviors are investigated by combining the theory with computer simulation software including ABAQUS and ANSYS, of which representative studies include: Zhou conducting research on damage evolution of 2-D plain weave under single axial and double axial load by applying multi-scale modeling [12]. Hsiao et al. [13] established the analytical model of laminated plates with fibers of undulating form, and qualitatively investigated the declining value of the corrugated fiber stiffness and strength under compressive load. Takeda established the constitutive relation of composite laminated plates with out-ofplane wave defects based on a two-step assumption method to predict the effective elastic modulus of 3-D composite laminate plates [14]. Zeng et al. [15] constructed composite laminated plate models with fibers in sine form and cosine form, and acquired the initial damage strength of the carbon fiber/epoxy resin composite laminated plate by applying this model based on the Tsai-Wu criteria as the failure criteria. Wu et al. [16] established the composite laminated plate model with corrugation by applying the finite element analysis software, predicted the equivalent stiffness of composite laminated plates with corrugation defects under compressive load, and analyzed the influence of relevant corrugation parameters to the mechanical performance of composite laminated plates.
In summary, there are many researches on the strength of out-of-plane folds on laminates, but little on the damage evolution of composite laminates with out-of-plane folds. In this paper, composite laminates with symmetrical plies are first produced in the experiment, and then out-of-plane wrinkles are introduced. The effect of out-of-plane wrinkles on the strength of composite laminates is studied. Meanwhile, a 3-D continuous damage model was established based on continuum mechanics. A 3-D model of the plates was produced using the ABAQUS numerical simulation software, and by applying the established 3-D continuous damage model and the 3-D Hashin failure criteria adapted with secondary strain in the UMAT subprogram, the progressive damage analysis acquired the damage evolution process.

Experiment
The raw material was fabricated with the T300/C4211 carbon fiber prepreg provided by Shandong AVIC Taida Composite Material Co., Ltd., by autoclave processing. The prepreg cloth has a thickness of 0.2 mm, and the prepreg fiber has a volume content of 65%. The laminated plate is comprised of 546 CMES, 2020, vol.124, no.2 , totaling 16 layers with a thickness of 3.2 mm after formation. The fold area is at the center of the top 2 layers 0°, the direction of the laminate and the fold position are shown in the following Fig. 1:

Geometrical Characteristics and Description of Samples
The composite laminated plates with out-of-plane fold in this article have similar form in the fold area with a cosine curve in appearance. The aspect ratio of the fold is used in this article to describe the degree of out-of-plane fold. The appearance characteristics and geometrical parameters of the fold are shown in the following Fig. 2: The aspect ratio j is applied in this article to describe the degree of fold with the following expression: where h is the height of the out-of-plane fold, and b is the range.

Sample Fabrication
The autoclave process in the experiment to fabricate test samples improves the stripe method provided by Wang to apply out-of-plane folds to the composite laminate plates, which is to fabricate the plates by inserting prepreg strip in the normal composite material laminates [17]. The characteristics of the formed out-of-plane folds are controlled by changing the inserted layers of the prepreg stripe. The inserted prepreg stripe in this experiment has a dimension of 10 mm × 15 mm and a laminate direction of 90°, and is inserted under the top 2 layers 0°. The insertion of 90°prepreg has two advantages [18]: On the one hand, that the prepreg direction is perpendicular to the loading direction would minimize the influence on the overall load-displacement; on the other hand, it is easy to define in simulation. The selection of laminate position in the prepreg relies on the fact that the out-of-plane folds mostly appear in the upper laminates, and the fold would not be clear if it were inserted in the medium layers.
Four groups of laminate plates are fabricated in this experiment by inserting 0, 2, 4 and 6 layers of prepreg strip respectively including control group. The sample information is shown in Tab. 1:  The strip is cut into rectangles of 30 cm × 35 cm and overlaid according to the design. Then cured in the autoclave. The laminated plates are cut into test samples of 140 mm × 15 mm after curing. The compressive performance of the sample is tested with the ZY-44 microcomputer controlled electronic universal testing machine according to the ASTM-D6641 [19] standard. The detailed process is shown in Fig. 3.

Experiment Results
Each group of molding samples are cut into five standard samples, and compressive performance tested according to the standard. The five data are then imported to the ORIGIN software for processing, with the stress-strain curve of each sample revealing the average value of the five test data. The calculated stress-strain curve of fold samples with different aspect ratios is shown in Fig. 4:  Compare the compressive strength of the fold with different aspect ratios with the flat plates to acquire the strength and decreasing ratio. The result is shown in Tab. 2.
Seen from Tab. 2, the critical failure strength of the laminated plate is 558.3 MPa without out-of-plane fold. With the increasing of the aspect ratio of the out-of-plane fold, the compressive strength of the sample decreases. When the aspect ratio reaches its maximum value of 0.12, the compressive strength decreases at the most by 36.1%. By drawing a smooth curve with the percentage decline according to the aspect ratio, the result is shown in Fig. 5: From Figs. 4 and 5, we can see that the decline of the compressive strength of the composite laminate plates relates to the aspect ratio of the fold. We can also find out some other influences of the out-of-plane fold from the fractured structure. The fractured structure is shown in comparison Fig. 6: About the two fractured structure samples in Fig. 6, the upper one is the control group without fold, and the lower one is a test sample with the aspect ratio of 0.12. Without fold we see the fracture occurs at the fixed end area of the sample, and some may occur at the load end together with slippage and debonding of some stiffeners. According to the design and theoretical assumptions, the mechanical properties of the whole test region of the plate should be uniform, so theoretically failure may occur at any place of the test area. The occurance at the fixed end might result from the stress concentration at the fixed end, which may be related to the load method and test standard. However, the fractured structure of the sample with an outof-plane fold should be close to the fold area under the same standard. This should be due to the stress concentration produced by the fiber bending at the fold. Therefore, we can conclude that the out-of-plane fold would not only decrease the strength of the laminate plate, but also affect its fracture position.

Progressive Damage Analysis
By experiment we have concluded that an out-of-plane fold would decrease the compressive strength as well as affect the fracture position. To further investigate the damage evolution behavior under compressive load, we apply the method of progressive damage analysis to conduct numerical simulation analysis. By studying the constitutive equation of the composite material and by combining the theory of continuum mechanics. We apply the constitutive equation of the composite material and the 3-D Hashin failure criteria as well as relevant material stiffness degradation criteria in the UMAT subprogram, so as to study the damage evolution behavior of different laminates of the composite laminated plates with out-of-plane folds.

Constitutive Equation of the Composite Material
Firstly we need to understand the constitutive equation of the composite material, which is as follows [20]: where s is the tensor of stress, ε is the tensor of strain, expressed as a matrix: e ¼ e 11 ; e 22 ; e 33 ; e 12 ; e 13 ; e 23 ½ T r ¼ r 11 ; r 22 ; r 33 ; r 12 ; r 13 ; r 23 ½ T C is the stiffness matrix of the material: The constitutive equation of the composite material is a mathematical model that reflects the macroscopic property of composite materials, and describes the damage to the composite material with the deformation value. In continuum mechanics, the simplest way to describe the material damage is to describe its internal damage deformation, quantizing the microscopic damage distribution with a proper tensor field. This analysis is focused on the compressive fracture of the composite material, so the tensor field is applied to describe the direction and density of the micro-fractured crack in materials. The micro Figure 6: Fracture of specimen with j ¼ 0 and j ¼ 0:12 crack of the composite material after damage would decrease the cross sectional area of the material from A to A. The stress would change due to the change of the cross sectional area, so s is applied to represent the stress after damage, which qualitatively should be greater than that before damage. The damage variable factor d is determined as the following equation: The stress r after damage is defined as: The material model in this article is the orthotropic composite laminated plate. So the orthotropic damage tensor matrix D is applied to promote calculation precision, which is shown as follows: where d i represents the decline proportion of the effective bearing area of the i-th main directional plane of the material. The first main direction of the material is the same as the fiber, while the second as the crosswise direction of the slice plane, and the third as the laminate overlap direction. The typical value of the damage tensor should be within the range of 0 ≤ di ≤ 1, where di = 0 shows that the plane of the i-th main direction of the material has no damage, and di = 1 indicates that the plane of the i-th main direction of the material has fractured structure. The stresss by applying the damage tensor matrix D can be defined as: where M (D) can be expressed as: where: where the strain energy of the material without damage by applying the effective stress in the equivalence criteria of strain energy could be expressed as: The strain energy of the material with damage can be expressed as: The strain occurs after the material is loaded, so the strain can be calculated as: The stiffness matrix of the composite material after damage is: where the non-zero component of C d is given as:

Practical Criterion of Damage
In the development of composite materials, there are many failure principles that judge the internal failure of composite materials, among which the commonly used are 3-D Hashin failure, Chang-Chang failure and Cai-Wu failure criteria [21][22][23]. This article applies the 3-D Hashin failure criteria. The 3-D Hashin failure criteria has many editions in which many scholars have made improvement. The 3-D Hashin failure criteria applied in this article is the practical criteria after secondary strain modification. Compared with regular stress failure criteria, the failure criteria applied in this article is relatively simple and easy to compile in the subprogram. In this article, F f , F m , F z represent fiber failure, substrate failure and interlayer delamination failure respectively, and the integration point only leads to compressive failure. The specific performance is shown as follows: The interlayer fractured structure is not taken into consideration at present due to axial compression. The strain ε of above all directions is expressed as:

Rule of Damage Evolution
When the material in a composite material structure gets damaged, the characteristics of the materials should be adjusted according to the degradation criteria, and the corresponding material property in finite element model should be modified according to specific degradation criteria. Slight has proposed a rather comprehensive degradation model in this aspect [24]. Among them there are instantaneous unloading model, linear degradation model and non-linear degradation model. The instantaneous unloading model can by applied to approximately estimate the critical failure strength of the material, but not to precisely simulate the process of material failure. The linear degradation model is generally applied for isotropic linear materials, but since the composite material is anisotropic, a non-linear degradation model is applied in this simulation for calculation to precisely simulate the failure behavior of the composite material. The nonlinear degradation model is shown below: where L c is the characteristic length of the unit, and G c,i is the critical fracture energy release rate of the material at the i-th direction. F f ; ðk ¼ f ; m; z) is the damage failure vector.

Finite Element Modeling
The mechanical property data of the carbon fiber cloth to fabricate the laminate plate is provided by Shandong Chonhunteda Composite Co., Ltd. With detailed parameters as shown in Tab. 3: During sample fabrication, we saw that the out-of-plane fold has the shape of a cosine ellipse as shown in Fig. 1, and the characteristics of the curve get closer with more layers inserted. Therefore, the sample with the aspect ratio of 0.12 is selected for the modeling and analysis, resulting in the sample characteristics close to a cosine curve. The top point of the fold and the two points at the edge of the scope are selected for the modeling of the fold area. The discrete coordinate system is applied for the overlaying of the composite laminates of the fold area. The test area takes mesh refinement due to the existence of the curve area by applying neutral axis drawing. The mesh grid has a dimension of 0.4 mm × 0.4 mm under the cloth with the boundary condition of both fixed ends, and the modeling is conducted to the clamp with the analytical rigid body. The central test area has a total of 13420 units. The boundary condition is set to have the both ends of the sample fixed, and the compressive test is simulated by applying displacement to the clamp. The universal frictional contact is set between the clamp and the composite material model to form constraints. Wherein above the boundary condition, the bottom two clamps are set with normal contact, and the upper two clamps have a displacement of 0.01 mm separately to simulate the tightening of the clamp. The left two clamps are fully constrained to simulate the fixation of the sample, and the right two clamps are set with a displacement of -0.6 mm to simulate the compressive test. The finite element model is shown in Fig. 7:

Optimization of Modeling
Since the 3-D Hashin failure criteria after secondary strain modification is applied, the strain energy would affect the calculation results. The calculated energy may decrease as the mesh grid becomes smaller, so the dimension of the mesh grid may affect the precision of the calculation result, especially for the existence of the fold curve in this model. To optimize the mesh grid, we chose to refine the mesh grid. So the fracture energy is applied in the variables of damage evolution to reduce the influence of the mesh grid dimension to the final result. In the implicit analysis of ABAQUS, the decline and contact on the unit stiffness will lead to non-convergence of the analysis. To ensure the convergence of the analysis, this article applies the viscidity regularization method in the subprogram to guarantee the tangent stiffness matrix of the damaged unit could stay positive definite in relatively small increment. The expression of the damage variable after applying viscidity regularization method is: where Δt is the time increment, d v j;old is the viscidity damage variable of the former step, and h is the viscidity coefficient, which is usually 10 -5 ∼10 -3 according to the research of Chen et al. [25], and is set to 0.00016 in this article.

Simulation Result and Damage Evolution Analysis
The computed initial failure strength is 335.41 MPa, and by comparing the simulation result with the test result of 356.2 MPa, And the relative error of the trend slope between the simulated stress-strain curve and the experiment result is less than 5%, and the error less than 10% is acceptable. The simulation environment is more ideal than the experiment, and the aspect ratio of the folded part in the simulation is the same as that in the experiment, but the folded part in the experiment is not as smooth as that in the simulation, and the error less than 10% is acceptable. The comparison results are shown in Fig. 8: At the same time, we can get the failure stress cloud diagram in the simulation, and define the output variables in the subroutine. Considering that the direction of the stress is axial compression, the interlayer damage is not obvious enough, so only fiber failure and matrix failure are considered here. The failure stress cloud diagram of the finite element model is shown in Fig. 9: In Fig. 10 we can also observe the failure form and fracture area. The fracture structure in the numerical simulation is located at the center of the fold area, and the main fracture structure has the form of fiber compressive fracture and substrate compressive fracture, which can be seen in the fracture sample in Fig. 11. The main fracture position and failure strength of the sample suit basically well with the experiment result. The ABAQUS analysis has the advantage that we can acquire the initial failure strength, and can meanwhile observe the damage evolution point of each different laminate, of which the analysis is shown in the following Figs. 11-15: In the above failure evolution of different laminates, SDV1 represents the fiber damage, SDV2 for the substrate damage, and Step Time for the starting time of the analysis. Figs. 11-15 shows the damage evolution nephogram of different layers obtained by ABAQUS. We can get the failure time of different layers, as shown in Tab. 4.
In the above Fig. 11 we can see that the 0°bottom layer firstly has fiber compressive failure with failure position at the center of the laminate plate. With the vertical expansion of its cross section, the substrate damage then occurs. In Fig. 12 we can see that the substrate damage of the 90°laminate precedes the fiber damage, followed by the vertical expansion. In Fig. 14 we can see that the failure of the 0°fold area starts last, and the damage position is different from the others, located at the bending area of the fiber. The failure time of each layer is shown in Tab. 4. The damage may also occur in the center, which should be interlayer expansion damage, and meanwhile may fit with the above idea that the fiber fold would lead to stress concentration. The fracture form and position as well as the initial failure strength of the laminated plates in real situations also match with the simulation on the whole.   1. The influence of the aspect ratio of out of plane folds on the compressive strength of composite laminates was studied by experiments. The results show that the out of plane fold will seriously affect the compressive strength of composite laminates, and with the increase of aspect ratio, the compressive strength of composite laminates is lower. In the experiment, when the aspect ratio is 0.12, the compressive strength of the laminate is reduced by 36.1%. 2. The ABAQUS software was used to model the j ¼ 0:12 specimen, and the damage evolution analysis was carried out with UMAT subroutine. The damage evolution nephogram of each layer is obtained. The failure modes and damage evolution of composite laminates with out of plane folds are explored. By comparing the failure strength and failure location between the simulation results and the experimental results, it is found that the two results are in good agreement. It also shows that the damage evolution nephogram obtained in this paper has certain reference for the same type of research.
Funding Statement: The author(s) received no specific funding for this study.

Conflicts of Interest:
The authors declare that they have no conflicts of interest to report regarding the present study.