Identification of Multiple Mechanical Properties of Laminates from a Single Compressive Test

In-plane elastic and interlaminar properties of composite laminates are commonly obtained through separate experiments. In this paper, a simultaneous identification method for both properties using a single experiment is proposed. The mechanical properties of laminates were treated as uncertainties and Bayesian inference was employed with measured strain-load curves in compression tests of laminates with embedded delamination. The strain–load curves were separated into two stages: the pre-delamination stage and the post-delamination stage. Sensitivity analysis was carried out to determine the critical properties at different stages, in order to alleviate the ill-posed problem in inference. Results showed that the in-plane Young’s modulus and shear modulus in elastic properties are dominant in the pre-delamination stage, and the interlaminar strength and type I fracture toughness in interlaminar properties are dominant in the post-delamination stage. Five times of property identification were carried out; the maximum coefficient of variation of identified properties was less than 1.11%, and the maximum error between the mean values of the identified properties and the ones from standard experiments was less than 5.44%. The proposed method can reduce time and cost in obtaining multiple mechanical properties of laminates.


Introduction
Carbon fiber reinforced composite has been widely used in the aerospace field, such as for the tail and fuselage of aircraft [1], due to its high specific strength, specific stiffness, and corrosion resistance [2,3]. Delamination is an important damage mechanism for this kind of material [4,5]. It is of great significance to obtain the mechanical properties of the laminate with delamination efficiently for evaluating the state of the laminate and predicting the remaining life [6][7][8]. The ability to obtain a single mechanical property, such as the elastic property and the interlaminar property, has formed a variety of standards. However, considering the cost of economy and time, how to obtain multiple properties through a single experiment has been of wide concern.
Molimard [9] identified multiple in-plain mechanical properties by using the tensile test of perforated thin plate, moiré technology and Levenberg-Marquardt (L-M) optimization method. Lecompte [10] carried out a biaxial loading test on glass fiber reinforced materials, obtained field strain information by an optical measurement method, and effectively identified multiple material performance properties based on a mixed numericalexperiment method. Lee [11] identified the elastic properties of flexibly supported rectangular laminated composite sandwich plates by using measured natural frequencies.
Zhuo [12] proposed an inverse method to identify the mechanical properties of fiber metal laminates on the basis of the measured and calculated frequency response functions. Michopoulos [13][14][15][16] used the energy method to construct a general material constitutive model, and utilized the independently developed 6-DoF loading system and complex sample design to realize properties identification based on combined use of MATLAB and Ansys. The maximum difference between the obtained identification results and the standard test result was no more than 3.5%. Chen [17] proposed an inverse method for identification properties of variable stiffness composite laminates based on the approximate Bayesian computation, which can avoid the calculation of complex likelihood. Bouhala [18] proposed a method to identify the interface properties of unidirectional carbon/epoxy composite, based on the results of double cantilever beam and the corresponding a counterpart extended finite element method cohesive zone model. Su [19] developed an automatic identification of the interfacial cohesive properties between fiber-reinforced polymers and concrete based on a machine learning-based artificial neural network. Based on the fullfield strain information obtained by digital image correlation and global sensitivity analysis method, Alfano [20] determined the main interface properties to be identified and the selection of the most suitable data area for identifying the properties. This provides guidance for the selection of appropriate observation data and experimental design.
Although much research has been done in identification of multiple properties of laminates, there is a lack of research on simultaneous identification of elastic properties and interlaminar properties in previous research. The reasons are as follows: (1) The key to identifying multiple properties through an experiment is that the experiment is sensitive to multiple properties. Therefore, how to design a single experiment sensitive to both elastic properties and interlaminar properties is a problem. (2) There are many elastic properties and interlaminate properties of laminates; identifying multiple properties at the same time can introduce the ill-posed problem due to the nature of the inverse problem. This means that the proper selection of the number of parameters to be identified at the same time is an important problem.
In this paper, a method of simultaneous identification of elastic properties and interlaminar properties is proposed, by using the characteristics of the test of the laminate with embedded delamination under compressive load. Whether the embedded delamination of the laminate extends under the compressive load is regarded as two different stages of parameter identification. Based on the results of sensitivity analysis (SA), the key elastic properties can be identified before the extension of the delamination under compressive load, and the key interlaminar properties can be identified after the extension of the delamination. In this way, two kinds of properties can be identified in one experiment, but the elastic and interlaminar properties are not identified at the same time, so as to alleviate the risk of an ill-posed problem of multiple properties identification.

Experimental Study
In this section, the response characteristics of the laminate with embedded delamination is illustrated by experiments, which can be regarded as the basis of the methodology proposed in Section 3.

Sample Preparation
The sample material of the laminate with delamination damage is carbon fiber reinforced epoxy composite. The laminate consists of 20 layers, each thickness is 0.104mm, and the layering sequence is [45/0//−45/0/45/90/−45/0/45/0/0/45/0/−45/90/45/0/−45/0/45] (//represents the position of the embedded delamination). The dimension properties of the laminate are shown in Figure 1. The delamination is a square with sides of 10 mm in the middle of the laminate. The embedded delamination was produced by placing polytetrafluoroethylene (PTFE) film (10 mm × 10 mm, 0.01 mm thickness) in the preset position before the material was cured and the material was cured by hot pressing method. The ultrasonic scanning results of the delamination areas of three samples are shown in Figure 2. Considering the clamping requirements, 50 mm aluminum reinforcing sheets were pasted on both sides of the sample, and strain gauges were pasted along the loading direction at the designed sites, as shown in Figure 1. Strain gauges 1 and 7 were used to observe the strain characteristics at the end of the sample. Strain gauges 2-6 and 8-12 were used to observe the strain characteristics of the surface of the embedded delamination region and its adjacent areas under compressive loading.

Experiment Process
The experiment was carried out on a universal strength testing machine Zwick/Z100 (ZwickRoell, UIm, Germany). The experiment partly referred to ASTM D7137 [21] such as clamping strategy and loading rate, and the experimental device is shown in Figure 3. Both sides of the sample were clamped by a wedge device, and then the wedge device was assembled into a rectangular fixture with a wedge groove to facilitate the testing machine to apply the compressive load. The complete experimental process is shown in Figure 4. The whole test adopted displacement control loading, the loading rate was 0.1 mm/min, and the preload was 50 N. The strain information at different points of the sample was obtained by strain collection device, and the sampling frequency was 100 Hz.

Experiment Result
The load-displacement curves of three experiments are shown in Figure 5. Taking sample 2 as an example, the strain-displacement curves at different measuring points (i.e., MPs in Figure 6) and the load-displacement curve are simultaneously given in Figure 6 for comparative analysis. It can be found that:

•
Although the slope of the load-displacement curves of each sample was slightly deviate in the early stage due to the difference of fixture installation and the small variations in the production procedure of the samples, the peak load of all samples was close to 4000 N.

•
The whole loading process can be divided into two stages: in stage 1, the load increased continuously, there was no local bucking, and the delamination did not extend; in stage 2, the load decreased gradually, local bucking began to occur, and the delamination began to extend in the direction perpendicular to the load. The information of local buckling can be proved by the strain characteristics of the MP4 pasted on the surface of the layered area as shown in Figure 1 (strain gage 04). When the load is near its peak, the strain of the MP4 changes rapidly from negative to positive, which represents a sudden change in the deformation characteristics of a region, from compression to tension, that is, outward buckling. Meanwhile, the slopes of strain-displacement curves of other MPs also changed greatly, representing the further intensification of the bending characteristics of the structure.

Framework
Considering the structural state of the laminate with embedded delamination under compressive load can be divided into two stages. In stage 1, there is no delamination expansion; the carrying capacity of the laminate is mainly reflected in the elastic properties of the material. In stage 2, due to the delamination expansion caused by local buckling, the carrying capacity is reflected in the elastic properties and interlaminar properties of the material. Therefore, based on such characteristics, a segmented method for identifying mechanical properties of the laminate is proposed. The elastic properties of the material are identified in stage 1, the calibrated properties are passed to stage 2, and the interlaminar properties of the material are identified in stage 2. The framework of parameter identification method is shown in Figure 7. The steps are as follows: 1.
Based on the prior information of the material properties, the initial distributions of m elastic properties p 1 of the material are given (elastic properties including Young's modulus, shear modulus, Poisson's ratio, etc.). Then, n samples are randomly sampled according to the distribution characteristics of each property and randomly grouped as G 1 m×n . Then, G 1 m×n is put into the finite element model (FEM) and the load response F 1 and the strain response R 1 of the measured points in Figure 1 are obtained.

2.
The load sensitivity varying with the strain of each property is calculated on the basis of method of SA according to the sample (G 1 m×n ) and the outputs (F 1 and R 1 ). Thus, the main sensitive properties p * 1 are determined.

3.
Based on the dynamic Bayesian network (DBN), p * 1 is identified using the time varying load-strain data of stage 1 in the experiment.

4.
The identified properties are regarded as the real mechanical properties of the material.
Step 1 is repeated among the l interlaminar properties p 2 (interlaminar properties including fracture toughness, interlaminar strength, etc.), forming the group G 2 l×n and outputs (load response F 2 and the strain response R 2 ). 5.
Determine the main sensitivity properties p * 2 of interlaminar properties according to step 2. Then, p * 2 is identified using the time varying load-strain data of stage 2 in the experiment based on the DBN.

Sensitivity Analysis Method
In this paper, the global sensitivity analysis (GSA) is adopted because it can consider the influence of the distribution interval of inputs on the outputs. There are many GSA methods [22][23][24][25][26][27]; the random balance design Fourier amplitude sensitivity test (RBD-FAST) [25] has been adopted in this paper. The FAST method is briefly introduced as follows: Let Y = f (x) be a model output with m random inputs (x = [x 1 , x 2 , . . . , x m ]), and assume that all the inputs can form a unit hypercube. Each input is represented by a given search function (G i (·)): where s is a scalar variable varying over the range −π < s < π, ω i is a specific frequency for each input variable, and G i (·) is the common search function as follows [27]: The outputs (Y) can be expanded into a Fourier series: For any positive integer (ω i ), the period T is 2π. The corresponding Fourier coefficients are defined as follows: The spectrum of the Fourier series is defined as follows: The first order of the variance-based global sensitivity of x i can be defined as: In the FAST method, the ω i needs to satisfy the linearly independent, which makes the computational cost of high-dimensional problems unbearable. To reduce the computational cost, the random balance design method is combined with the FAST [25]. In RBD-FAST, the ω i takes the same value for each input variable. The method distinguishes the characteristics of the different inputs by controlling the ordering of the inputs without calculating a large number of samples, which can reduce the computing cost effectively. The process can be summarized as follows: 1.
Determine s= [s 1 , s 2 , . . . , s N ] and sort s randomly. The disordered sample of x is then obtained by random permutation of s, and the model outputs (Y) are calculated by a disordered sample of x.

2.
The outputs (Y) are rearranged according to the original sequence of the sample of x i . The Fourier transform is then performed according to Equation (5). 3.
The first-order sensitivity index corresponding to x i can be obtained through Equations (6)- (8). Return to step 2 to determine the sensitivity of other inputs.

Dynamic Bayesian Network
There are many classical methods for parameter identification, such as the iterative regularization method [28], the Tikhonov regularization method [29], and Bayesian method [30]. When the problem of parameter identification has the characteristics of large number parameters to be identified, model complex and time-dependent data, dynamic Bayesian networks (DBNs) are an excellent choice [31]. Therefore, a DBN is adopted in this paper for parameter identification.
A DBN is an extension of a Bayesian network (BN) in time domain. A BN is a directed acyclic graph model used for uncertainty inference. In a BN, random variables are represented vertices and their dependence are represented by directed edges. The quantitative characterization of dependence can be expressed by a probability density function or a definite function [32]. Compared to a BN, a DBN has additional lines connecting the same variable between two adjacent time points. This allows a DBN not only to integrate a BN content, but also to accumulate previous knowledge. There are many inference algorithms that can be used for a DBN, including multiple Kalman filters [33][34][35] and particle filters (PFs) [36]. Considering the no-linear characteristics of the model studied in this paper, PFs is selected as the inference algorithm of the DBN. The sequential importance resampling (SIR) algorithm, which is one of the PF algorithms, is employed in this paper.
The DBN depicted in Figure 8 is used to briefly introduce the SIR algorithm. As shown in Figure 7, the state variable X evolves over time according to the state function: where ε is noise terms in the state function. The measurement variable Z is calculated by the measurement function: where η is noise terms in the measurement function. The SIR algorithm is as follows [37]: 1.
The initial particles x i 0 N i=1 are generated according to the prior probability density function (PDF).

2.
Loop the following steps at t = 1, 2, . . . , : (1) Sampling from the proposal PDF: generating particles and calculating the corresponding weights ω i t according to Resampling and estimating: Particles set is resampled to

Finite Element Simulations
The finite element model of the laminate has been established based on the commercial software ABAQUS (Dassault Systemes SIMULIA, Paris, France). The modeling procedure are shown in Figure 9. The 8-nodes 3D solid composite elements are used to model the laminate (C3D8RC3), which can provide accurate interlaminar stress and transverse shear effect [38]. To predict the delamination evolution, the interface between two sub-laminates is modeled by 8-nodes 3D cohesive elements (COH3D8). The volume thickness of the cohesive element is zero, and the bilinear traction-separation law has been employed herein (refer to [39]). In Figure 9, the left side of the laminate is completely fixed (U x,y,z = 0), and the compressive displacement is applied to the right side (U y,z = 0). The laminate is then divided into two zones artificially to save the computing cost. In zone 1, where sub-laminate local buckling occurs and delamination may propagate, the cohesive element has been positioned in the area as shown in Figure 9. Surface-to-surface contact element has been placed in the delamination zone to avoid overlaps between elements. The minimum mesh size is 0.5 mm. In zone 2, only global buckling occurs and the minimum mesh size is 0.75 mm. The quadratic stress failure criterion (QUADs) has been used as the delamination initiation criterion to calculate the delamination evolution under mixed stress state [38], which is expressed as: where T, S are the maximum tractions in the normal, and shear directions, respectively. T n , T s , T t are the traction components in the normal, first, and second shear directions. When the damage occurs, the damage evolution continues in the cohesive zone. An energy-based evolution model was adopted in this paper. Considering the mixed-mode behavior, the power-law criterion was used [38], which is expressed as: where G IC , G IIC , G IIIC are the critical energy release rates (fracture toughness) for mode I, II, III, respectively. G I , G II , G III are the corresponding release rate values in analysis and α is a coefficient. The nonlinear solution of the problems presented here is performed using standard ABAQUS procedures. The grid independence verification has been verified by 1.4 times number of grids, and the results of load-displacement curves are consistent. Note that the initial size and position of the embedded delamination are fixed in the simulation of this paper. The mechanical properties discussed in Section 3.1.1, including elastic properties and interlaminate properties, are the variable inputs of the model, which affect the response characteristics of the structure under displacement load. Other parameters remain unchanged during the simulation.

Verification Example
The verification example is based on the experiment in Section 2, and the process of the verification example is as follows: 1.
The initial distribution intervals of 12 material properties of the laminate are shown in Table 1. A total of 1000 samples are randomly selected for each property, and the samples of each property are randomly combined to form a 1000 sample set G 1 12×1000 . Then, bring G 1 12×1000 into the finite element model in Section 3.2 to obtain the strain response of the concerned measuring points under different loads.

2.
The sensitivities of properties of the laminate at MP9 (just for illustration, other measuring points can also be used) under different strains in stage 1 are analyzed to obtain the key property p * 1 . Then, p * 1 is identified based on the DBN along the loading process as shown in Figure 6. The observation points in DBN are shown in Table 2 and the observation error is set as 1% of the observation data (Observation data is obtained from MP9 in stage 1).

3.
Similar to step 1, resample 1000 samples of interface properties to form the set G 2 6×1000 . The sensitivities of properties of the laminate under different strains in stage 2 are analyzed to obtain the key property p * 2 . Then, p * 2 is brought into the finite element model to calculate the strain response of MP9. At this time, the input elastic properties of the model are the properties just identified. p * 2 is identified along the loading process in DBN as shown in Figure 7. The observation points are shown in Table 2 and the observation error is set as 1% of the observation value (Observation data is obtained from MP9 in stage 2).  There are some notes as follows: (1) All the mechanical properties are used in the sampling in step 1 for SA in the example, the purpose of which is to further prove that the interface parameters have no effect in the stage 1.

Results and Discussion
First, a sample (E 1 = 165.56 GPa, G 12 = 4.17 GPa, G IC = 240.94 N/m, T = 10.09 MPa, S = 53.85 MPa) in step 1 of Section 3.3 was randomly selected to show the results of its finite element calculation, as shown in Figure 10. The characteristics of the load-displacement curve and the strain-displacement curve at MP9 are completely consistent with the corresponding curves of the experiment in Figure 6. In the simulation, it can be clearly seen that near the inflection point of the load-displacement curve, the delamination expansion caused by local buckling begins to occur, which further verifies the conclusion given by the experiment. As the load gradually increases, the delamination gradually extends along the direction perpendicular to the load, which is consistent with the conclusions of [38,39].
Then, the property sensitivity in stage 1 is shown in Figure 11. It can be seen that the sensitivity characteristics of properties do not change significantly at different MPs. Therefore, the calculation example only uses the data of MP9 to illustrate the method. In stage 1, the dominant properties are E 1 and G 12 , and E 1 is absolutely dominant in the whole mode. This is mainly because the laminate does not occur delamination extension in stage 1, the deformation of the structure along the loading direction is more obvious under compression, and the direction of strain extracted is the same as the loading direction. Therefore, the modulus along the loading direction E 1 and the in-plane shear modulus related to the loading direction G 12 have important influence on the load.  According to the results of SA, only properties E 1 and G 12 need to be identified in stage 1. Then, based on the DBN, data of five observation nodes have been used to identify E 1 and G 12 . In order to ensure the robustness of the identification results, a total of five identifications were carried out and their mean values were compared with those obtained by standard tests as shown in Table 3. It can be seen that the coefficient of variation of the two properties is less than 0.73%, and the maximum error with the properties obtained from the standard tests is no more than 1.57%, which proves the robustness and accuracy of the identification method. The property sensitivity in stage 2 is shown in Figure 12. At this time, only the sensitivity of interlaminar properties is considered, and the elastic properties obtained by stage 1 are taken as the true values of the model. In stage 2, the dominant properties are G IC , T, and S. With the increase of the strain (i.e., the gradual extension of the delamination), the sensitivity of G IC gradually increases while G IIC does not. This shows that the delamination characteristics of the laminate under compressive load are mainly related to the mode of G IC . To our surprise, the sensitivity of penalty stiffness K and Power-law parameter α can be ignored. The identification results of the three interlaminar properties are shown in Tables 3 and 4. The coefficient of variation of the three properties is less than 1.11%, and the maximum error with the properties obtained from the standard tests is no more than 5.44%. The posterior distribution of the five identified properties is shown in Figure 13. It can be found that the posterior distribution characteristics of the properties of five times identification are similar, which further proves the sampling 1000 particles can meet the robustness requirements of the method.   In order to show that the segmented parameter identification method proposed in this paper can improve the accuracy of parameter identification, the identification results of parameters were compared with the results of Levenberg-Marquardt (L-M) method, which is a classical method in iterative regularization [9]. When using the L-M method, five main parameters are identified at the same time. The comparison results are shown in Table 4. It can be found that the accuracy of parameter identification using the proposed method is higher than that using the L-M method, which identifies five parameters at the same time. The maximum error of the identification results of the two method is 5.44% to 8.61%.
Compare the load-strain curve of the first identification result (i.e., E 1 = 142.85 GPa) with the experiment as shown in Figure 14. In Figure 14, Experiment means the experimental curve; Initial bounds represent the upper and lower bounds of load-strain curves for the initial distribution properties; Mean (Stage 1) and Mean (Stage 2) represent the mean value of the identification results of stage 1 and stage 2, respectively; UB (Stage 1) and UB (Stage 2) represent the uncertainty bounds of the identification results in stage 1 and stage 2, which is calculated by mean ± 2 × standard deviation and mean ± 4 × standard deviation, respectively. It can be seen that the load-strain curves obtained by property identification are in good agreement with the experimental data. The uncertainty interval of load-strain curve formed by identification can effectively envelop the experiment data. After determining the key elastic properties (E 1 ,G 12 ) in stage 1, the whole uncertainty interval is significantly reduced. This means the elastic properties play a crucial role in the load-strain curve, and further explains the significance of obtaining the elastic properties of the material before identifying the interlaminar properties.
Furthermore, the mean value of the identification results is brought into the finite element model to obtain the load-strain curves at other measuring points, and the curves are compared with corresponding experimental data as shown in Figure 15. In Figure 15, Cal(MP3) and Cal(MP8) represent the results of FEM at MP3 and MP8, respectively; Real(MP3) and Real(MP8) represent the experiment data at MP3 and MP8, respectively. It can be found that the identified properties can effectively characterize the strain-load curves of other measuring points, but there is a certain deviation when the strain reaches around 1250µε (inflection point of the strain-load curve). This is mainly due to the following two reason: (1) the discreteness of material properties at different positions; (2) the finite element model is too ideal to fully reflect all the real state when the structural state sudden changes.

Conclusions
In this paper, a multiple mechanical properties identification method is developed based on the features that the response of defective laminate under compressive load has two stages. In stage 1, there is no extension of the delamination. The results of SA show that only two elastic properties E 1 , G 12 need to be identified, and E 1 has a dominant position in the whole mode. In stage 2, the delamination begins to extend, and the results of SA show that the key interlaminar properties to be identified are fracture toughness G IC , interlaminar strength T, S. The sensitivity of G IC increases with the increase of delamination expansion. The penalty stiffness K and power-law parameter α have little effect in the whole stage 2.
In order to show the robustness and accuracy of the property identification method, five attempts at property identification were carried out and the results were compared with the results of the standard tests. The results show that the maximum coefficient of variation of the five identified properties was less than 1.11%, and the maximum error of the mean of identification results compared with the standard tests was less than 5.44%. Furthermore, in order to prove that the segmented parameter identification method proposed in this paper can improve the accuracy of parameter identification, the identification results of parameters were compared with the results of L-M method, which identifies five parameters at the same time. The results show that the accuracy of parameter identification using the proposed method was better (maximum identification error 5.44% vs 8.61%).
After determining the key elastic properties E 1 , G 12 in stage 1, the whole uncertainty interval of strain-load curve tracking is significantly reduced in stage 2. This means the elastic properties play a crucial role in both two stages. Identifying the interlaminar properties based on the identified elastic properties can improve the identification accuracy.
In future work, specific regions could be selected for property identification according to sensitivity characteristics based on the measurement technology of digital image correlation, so as to improve the sensitivity of multiple properties and the identification accuracy at the same time.