Experimental Study of the Probabilistic Fatigue Residual Strength of a Carbon Fiber-Reinforced Polymer Matrix Composite

Degradation of the mechanical properties of fiber-reinforced polymer matrix composites (PMCs) subjected to cyclic loading is crucial to the long-term load-carrying capability of PMC structures in practice. This paper reports the experimental study of fatigue residual tensile strength and its probabilistic distribution in a carbon fiber-reinforced PMC laminate made of unidirectional (UD) carbon-fiber/epoxy prepregs (Hexcel T2G190/F263) with the ply layup [0/±45/90]S after certain cycles of cyclic loading. The residual tensile strengths of the PMC laminates after cyclic loading of 1 (quasistatic), 2000, and 10,000 cycles were determined. Statistical analysis of the experimental data shows that the fatigue residual tensile strength of the PMC laminate follows a two-parameter Weibull distribution model with the credibility ≥ 95%. With increasing fatigue cycles, the mean value of the fatigue residual strength of the PMC specimens decreased while its deviation increased. A free-edge stress model is further adopted to explain the fatigue failure initiation of the composite laminate. The present experimental study is valuable for understanding the fatigue durability of PMC laminates as well as reliable design and performance prediction of composite structures.


Introduction
Fiber-reinforced polymer matrix composites (PMCs) made of high modulus microfibers (e.g., carbon, glass, boron, natural plant fibers, etc.) in compliant polymeric resins represent a large class of manmade lightweight structural materials that have been finding rapidly expanding applications in aerospace, aeronautical and ground vehicles, sports utilities, gas and liquid pipeline transportation, offshore and marine construction, and civil infrastructure, among others [1][2][3][4][5][6][7][8][9][10][11]. Such broad applications are resulting mainly from the unique properties of PMCs including their high specific stiffness and strength, excellent fatigue tolerance and immunity to corrosion, amicable manufacturability, etc., [1][2][3]5]. As microstructures, PMCs are highly heterogeneous materials, and when subjected to external loading, high stress inhomogeneities and stress singularities exist in PMCs, which result in their progressive damage [12]. Typically, the evolution of structural damage in PMCs includes microcrack nucleation, growth and propagation in resins, delamination between neighboring laminas, fiber debonding, fiber break, and, finally, catastrophic failure of the PMC specimens or structures [12][13][14]. Obviously, the evolution of this microscale to macroscale damage in PMCs is thermodynamically irreversible, which gradually degrades the load-carrying capacity of the PMCs until catastrophic failure happens. With the ever expanding applications of PMCs in transportation and other engineering sectors in 2 of 22 recent years, the long-term mechanical and other physical behaviors of PMCs subjected to cyclic loading have become one of the important research topics in composite design, manufacturing, testing, and applications [15][16][17]. To date, substantial experimental and theoretical investigations have been conducted to understand the fatigue behaviors of composite materials subjected to cyclic loadings.
To mention a few, Diao et al. [18][19][20] studied the fatigue response of carbon-fiber/PEEK composite laminates and formulated a statistical physics-based fatigue model to interpret their observations. Wu [21] conducted the thermomechanical fatigue tests of advanced composite laminates made of unidirectional (UD) carbon-fiber/epoxy prepregs and detected two opposite thermal effects in the mechanical behavior of such composites, i.e., thermally triggered fatigue failure and resin toughening, which can be related to the rate effect of PMC fracture [22]. Dzenis et al. [23][24][25] performed the fatigue tests on a UD carbon-fiber/epoxy composite laminate and studied the cycle-based damage accumulation in the composite laminates by using an acoustic emission (AE) technique. In the specific case of adhesively bonded composite joints, they further identified the fatigue fracture modes by comparing the collected AE signals with those obtained in pure mode fatigue fracture tests. Rudov-Clark and Mouritz [26] conducted the tension-tension fatigue tests of a three-dimensional (3D) orthogonal woven composite and found that the fatigue response of such a 3D woven composite was more complicated than that of composite laminates due to its 3D woven fiber architecture and related inhomogeneous stress and strain states. In order to understand the stiffness degradation of UD fiber-reinforced composite laminates subjected to cyclic loads, Cista et al. [27] formulated a progressive damage model based on the concept of fiber fragmentation, in which fatigue damage was simplified as cumulative fiber breakage. Moreover, significant experimental studies have also been performed to explore the fatigue fracture of fiber-reinforced composite laminates with precracks (delamination) and notches [28][29][30][31][32][33][34], which were focused mainly on the pure or mixed mode fatigue crack growth rate with respect to the mean value and amplitude of the stress intensity factor (SIF) and frequency of the applied cyclic loadings. So far, standard test methods for pure and mixed mode fatigue tests of fiber-reinforced composite laminates have been formulated [35,36]. Detailed reviews on recent experimental research of fatigue delamination of fiber-reinforced composite laminates and related fatigue damage modeling are available in the literature [37][38][39][40].
In addition, it is critically important to understand and predict the reliable lifetime and progressive evolution of the mechanical properties of composite materials subjected to cyclic loading. In the past three decades, substantial experimental and theoretical studies have been devoted to determining and modeling the evolution of the residual strength of composite laminates in the literature [41][42][43][44][45][46][47][48][49][50][51][52][53][54][55].
To mention a few, in an earlier study of the fatigue residual strength of composite materials, Yang [41] proposed a phenomenological fatigue residual strength model with a linear strength degradation with respect to the fatigue cycle, which was largely validated by the detailed fatigue tests of a quasi-isotropic graphite/epoxy composite laminate with a layup of [0/45/90/-45 2 /90/45/0] 2 with fitted model parameters. Hashin [42] utilized cumulative damage theory to formulate a residual life and residual strength model based on the master SN curves, which can deal with multistage loading cases. Diao et al. [44] formulated a stochastic model of the residual strength and fatigue life of composite laminates to take into account the governing damage modes in the laminate such as matrix cracking and delamination. This model was validated by their fatigue test data of a UD glass-fiber/epoxy laminate. Paepegem and Degrieck [46] utilized the concept of damage growth rate to formulate a kinetic fatigue damage model to predict the residual strength of composite laminates. Philippidis and Passipoularidis [48] performed the comparative studies and review of several existing residual strength models and their experimental validation in the literature. They concluded that the use of complicated phenomenological models required large experimental data sets for implementation and did not necessarily pay back in terms of accurate predictions and that simple models requiring limited experimental efforts could be preferable. Mejri et al. [51] conducted the bending fatigue tests of natural fiber-reinforced composites subjected to cyclic loading of constant strain level. They defined the residual fatigue strength as the maximum flexural stress of the sample, which progressively degraded with cycles as the result of cumulative specimen can only be used to generate one sampling data point. So far, the literature database of fatigue residual strength of composite laminates after fatigue loading of certain cycles is still limited. Therefore, in this work we conducted the experimental study to determine the fatigue residual tensile strength of an angle-ply carbon-fiber/epoxy composite laminate made of UD carbon-fiber/epoxy prepregs (Hexcel T2G190/F263) with the ply layup [0/±45/90] S after fatigue loading of certain cycles. The quasistatic tensile strength and the residual tensile strengths of the PMC laminate after cyclic loading of 2000 and 10,000 cycles were determined, respectively. Two-parameter Weibull model was used for statistical analysis of the experimental data to examine the variation of the model parameters with fatigue cycles. A free-edge stress model was further adopted to explain the fatigue failure initiation of the composite specimens. Finally, conclusions of the present study were drawn.

Composite Laminate Processing and Specimen Preparation for Quasistatic Tension and Fatigue Tests
The thermosetting polymer composite laminate used for this study was made of Hexcel T2G190/F263 supplied by Hexcel Co., Stamford, CT, USA. The UD carbon-fiber/epoxy prepregs were reinforced with UD T300 carbon fibers. The in-plane mechanical properties of the resulting UD carbon-fiber/epoxy laminate are listed in Table 1 [14]. In this experimental investigation, four laminated panels with a quasi-isotropic layup of [0/±45/90] S were assembled following a hand lay-up procedure. Each laminated panel was cured within a vacuum-assisted two-chamber press-clave (as illustrated in Figure 1), which was assembled into a 15-ton hydraulic hot-press installed with an electrical heating unit (the Carve Inc., Wabash, IN, USA) under controlled temperature, pressure and vacuum environment. The manufacturer recommended curing cycle was applied as (2) 120 min at 350 ± 5 • F (177 ± 3 • C); (3) −5 • F/min to 140 ± 5 • F (± 3 • C); (4) 80 ± 5 psi (0.55 ± 0.03 MPa). Table 1. In-plane mechanical properties of the UD carbon-fiber/epoxy composite laminate made of Hexcel T2G190/F263 prepregs [14].

Mechanical Properties in Material Coordinate System Values
Young's modulus in the fiber direction, E 1 (GPa) 132.7 Young's modulus transverse to fiber direction, E 2 (GPa) 8     Schematic assembly of the two-chamber press-clave for PMC curing [14].
The curing temperature profile of the Hexcel T2G190/F263 prepregs is shown in Figure 2. The laminated panels after curing were tabbed using composite strips, which were cut from a commercial woven glass-fabric/epoxy composite laminate, and a Miller Stephenson two-part adhesive 907. The tabbing can prevent premature failure of the composite near the specimen grips at the two ends and ensure that the final failure happened at the middle segment of the specimens. Rectangular composite specimens were machined from the tabbed composite panels at room temperature using a high-speed diamond saw installed with a water-cooling system, as shown in Figure 3. The specimen dimensions were 135 mm in length, 15 mm in width, and 1.35 mm in thickness. All the free-edges of the composite specimens were carefully polished using fine sandpaper to avoid premature edge delamination.

Quasistatic Tension and Tension−Tension Fatigue Tests of the PMC Specimens
Quasistatic tension and tension-tension fatigue tests were performed on an MTS servo-hydraulic testing machine installed with an Instron data acquisition and reduction system. Quasistatic uniaxial tension tests based on ASME standard D3039/D3039M-17 [113] were utilized to determine the static ultimate tensile strength and its probabilistic distribution of the composite laminate used in this study. During the quasistatic tension tests, displacement-control method with the loading rate of 1 mm/min was specified and 27 specimens were tested.
For the tension-tension fatigue tests of the composite specimens, the ASTM standard D 3479/D3479M-19 was maintained [114]. The fatigue peak loading stress was specified as 289. 85 MPa, which is close to 0.85 of the mean static tensile strength of the composite laminate that was determined from the quasistatic tension tests. Such fatigue peak loading stress does not evoke the first ply failure of the composite laminate according to the classic laminate theory (CLT, See Section 3). During the tension-tension fatigue tests, triangular loading wave was specified, the loading ratio R (=σ min /σ max ) was set as 0.1, and the fatigue frequency was selected as 5 Hz. Two sets of 12 composite specimens with the dimensions and processing procedure the same as those for quasistatic tension tests were used for fatigue tests of fixed cycles of 2000 and 10,000, respectively. After each fatigue test, the quasistatic tension test was further performed to determine the residual ultimate tensile strength of the surviving specimen using a displacement-control method with the loading rate of 1 mm/min. Figure 4 shows the typical quasistatic tensile stress-strain diagrams of the composite specimens without fatigue loading, after 2000-cycle fatigue loading, and after 10,000-cycle loading. It was found that the effective tensile modulus of the composite specimens slightly decreased after certain cycles of fatigue loading due to the cumulative damage. Stress drops before the final failure could be detected in each stress-strain diagram due to the progressive damage with increasing loading level. Experimental data of the quasistatic tension tests show that the tensile strength of the composite specimens has a wide range of distribution for either quasistatic loading or fatigue loading for 2000 or 10,000 cycles. As a matter of fact, for an angle-ply fiber-reinforced PMC laminate, a variety of random factors, such as fiber strength, resin strength and toughness, fiber alignment, fiber/matrix bonding strength, ply orientation, voids, etc., influence the effective mechanical properties of the laminate. Thus, with sufficient testing specimens, it is more suitable to adopt probabilistic models for data reduction of the quasistatic tensile strength and fatigue residual tensile strength of the composite laminate, which would provide more specific information to understand the statistical behavior of the material properties for reliable structural design and analysis. Hereafter, we introduce the probabilistic models of the quasistatic and fatigue residual tensile strengths of the composite laminate for experimental data reduction. In the view of reliability theory, the survival probabilities of the quasistatic tensile failure and fatigue residual tensile failure of the composite laminate are obtained through the median rank formula [51,106,115]:

Data Reduction and Probabilistic Residual Strength Model Fitting
where i is the i-th specimen for a sample size of n specimens in an increasing tensile or fatigue residual tensile strength sequence. It is assumed that both the quasistatic tensile strength and the fatigue residual tensile strength of the tested composite laminate obey three two-parameter Weibull distributions with different parameters, each of which corresponds to the survival probability at tensile stress σ in a quasistatic tensile failure (either before or after fatigue loading of certain cycles) as and the corresponding mean value of the strength is In the above, m and σ0 are the shape and scale parameters of a two-parameter Weibull distribution, respectively, and Γ( ) is the Γ-function defined as The two unknown parameters m and σ0 for each two-parameter Weibull model of either the quasistatic tensile strength or the fatigue residual strength after fatigue loading of certain cycles can Hereafter, we introduce the probabilistic models of the quasistatic and fatigue residual tensile strengths of the composite laminate for experimental data reduction. In the view of reliability theory, the survival probabilities of the quasistatic tensile failure and fatigue residual tensile failure of the composite laminate are obtained through the median rank formula [51,106,115]: where i is the i-th specimen for a sample size of n specimens in an increasing tensile or fatigue residual tensile strength sequence. It is assumed that both the quasistatic tensile strength and the fatigue residual tensile strength of the tested composite laminate obey three two-parameter Weibull distributions with different parameters, each of which corresponds to the survival probability at tensile stress σ in a quasistatic tensile failure (either before or after fatigue loading of certain cycles) as and the corresponding mean value of the strength is 8 of 22 In the above, m and σ 0 are the shape and scale parameters of a two-parameter Weibull distribution, respectively, and Γ( ) is the Γ-function defined as The two unknown parameters m and σ 0 for each two-parameter Weibull model of either the quasistatic tensile strength or the fatigue residual strength after fatigue loading of certain cycles can be determined by means of the maximum-likelihood estimation method, as shown in Table 2. In each case, the two-parameter Weibull distribution of the tensile strength has an acceptable credibility R ≥ 0.95, which indicates that both the quasistatic tensile strength and the fatigue residual strength of the composite laminate obey the two-parameter Weibull distribution very well. It can be also found from Table 2 that the values of both the shape parameter m and scale parameter σ 0 of the two-parameter Weibull model decrease with increasing cycle of fatigue loading. Table 2.
Model parameters of the two-parameter Weibull distribution models based on maximum-likelihood estimation. Survival probabilities of the quasistatic tensile strength and the fatigue residual strengths of the composite laminate after fatigue loading of 2000 and 10,000 cycles are plotted in symbols with significant scattering in Figures 5 and 6. Correspondingly, the survival probabilities predicted by the three Weibull distribution models are plotted with curves in Figure 5 and lines in Figure 6 as well. The failure probability densities of the quasistatic ultimate tensile failure and the fatigue residual ultimate tensile failure are shown in Figure 7. It can be found from Figure 7 that given the fatigue loading parameters (e.g., the peak loading stress, stress ratio, and loading wave shape), the mean value of the fatigue residual tensile strength decreases with increasing cycles of fatigue loading; in contrast the distribution of the fatigue residual tensile strength becomes more and more scattered with increasing cycles of fatigue loading. The statistical results based on the above two-parameter Weibull distribution models are reasonable to correlate to the experimental observations. In addition, the probabilistic models of the quasistatic tensile strength and fatigue residual tensile strength of the composite laminates proposed in this study could be used for the reliable design of composite structures and predication of the reliable lifetime of composite structures subjected to external loading.

No. of Samples
cycles of fatigue loading. The statistical results based on the above two-parameter Weibull distribution models are reasonable to correlate to the experimental observations. In addition, the probabilistic models of the quasistatic tensile strength and fatigue residual tensile strength of the composite laminates proposed in this study could be used for the reliable design of composite structures and predication of the reliable lifetime of composite structures subjected to external loading.

Model-Based Tensile Strength and Free-Edge Stress Predictions
Based on the quasistatic mechanical properties of the UD composite laminates made of Hexcel T2G190/F263 prepregs as tabulated in Table 1, the theoretical quasistatic tensile stress-strain diagram of the present [0/±45/90] S composite laminate under uniaxial tension is predicted according to the classic laminate theory (CLT) [1] as shown in Figure 8, from which the predicted first ply failure occurs at the 90 • plies, corresponding to the effective tensile stress of 311.11 MPa, provided that the ply failure of the composite laminate obeys the conservative maximum stress failure criterion as adopted in the CLT-based modeling. The theoretical predictions of the effective modulus and quasistatic tensile strength of the composite laminate are reasonably close to those determined in the present experiments, as shown in Figures 4 and 5.  In addition, both the quasistatic tension and fatigue tests showed that edge-delamination first initiated at the 45°/90° interfaces of the composite laminate, i.e., the maximum normal stress σzz or out-of-plane shear τxz first appeared at the 45°/90° interfaces. To validate such maximum stresses at the 45°/90° interfaces of the laminate, the high-efficiency stress-function variational method initially formulated by Yin [79,80] and further refined and extended by Wu [14,45] is adopted for the freeedge stress analysis of this composite laminate (See Appendix A for detailed derivations). The mechanical properties of the UD Hexcel composite laminate provided in Table 1 are used in this freeedge stress analysis, and additional mechanical properties to be used are derived or reasonably assumed as G13 = 4.76 × 10 3 MPa, G23 = 3.024 × 10 3 MPa, ν23 = 0.46. The composite laminate with the layup of [0/±45/90]S is assumed to have a straight free-edge and a constant uniaxial strain εxx = 0.01 is applied along the global x-axis, as shown in Figure A1 in Appendix A. The predicted free-edge stresses of the laminate are shown in Figure 9, in which Figure 9a-d show the stress variations along the 45°/−45° interface within L-distance (L: Laminate thickness) from the free edge (left) and Figure  9e,f show the variations of the normal stress σzz and out-of-plane shear stress τxz across the laminate thickness at the free edge. It can be found from Figure 9e,f that the highly localized free-edge normal stress σzz and out-of-plane shear stress τxz reach the maximum values at the 45°/90° and −45°/45° interfaces, respectively, where possible edge-delamination initiation could appear when the laminate subjected to uniaxial tension along the x-axis as observed in the experiments. Therefore, the above stress-function variational method is a powerful tool that can be used conveniently for prediction of the edge-delamination initiation in angle-ply composite laminates subjected to quasistatic tension or tension−tension fatigue loading. In addition, this semianalytic method can be used for laminate layup sequence design and optimization for effective suppression of the high free-edge stresses in composite laminates. In addition, both the quasistatic tension and fatigue tests showed that edge-delamination first initiated at the 45 • /90 • interfaces of the composite laminate, i.e., the maximum normal stress σ zz or out-of-plane shear τ xz first appeared at the 45 • /90 • interfaces. To validate such maximum stresses at the 45 • /90 • interfaces of the laminate, the high-efficiency stress-function variational method initially formulated by Yin [79,80] and further refined and extended by Wu [14,45] is adopted for the free-edge stress analysis of this composite laminate (See Appendix A for detailed derivations). The mechanical properties of the UD Hexcel composite laminate provided in Table 1 are used in this free-edge stress analysis, and additional mechanical properties to be used are derived or reasonably assumed as G 13 = 4.76 × 10 3 MPa, G 23 = 3.024 × 10 3 MPa, ν 23 = 0.46. The composite laminate with the layup of [0/±45/90] S is assumed to have a straight free-edge and a constant uniaxial strain ε xx = 0.01 is applied along the global x-axis, as shown in Figure A1 in Appendix A. The predicted free-edge stresses of the laminate are shown in Figure 9

Concluding Remarks
In this experimental study, quasistatic tension and tension−tension fatigue tests have been conducted for the determination of the quasistatic tensile strength and fatigue residual tensile strength and related probabilistic distributions of a composite laminate with a ply layup of [0/±45/90]S. Several conclusions can be drawn from the presented test results, data reduction and model-based simulations as follows.
(1) The quasistatic tensile strength and fatigue residual tensile strength of the tested Hexcel composite laminate follow two-parameter Weibull distributions very well with the shape and scale parameters decreasing with increasing fatigue cycles.

Concluding Remarks
In this experimental study, quasistatic tension and tension-tension fatigue tests have been conducted for the determination of the quasistatic tensile strength and fatigue residual tensile strength and related probabilistic distributions of a composite laminate with a ply layup of [0/±45/90] S . Several conclusions can be drawn from the presented test results, data reduction and model-based simulations as follows.
(1) The quasistatic tensile strength and fatigue residual tensile strength of the tested Hexcel composite laminate follow two-parameter Weibull distributions very well with the shape and scale parameters decreasing with increasing fatigue cycles. Furthermore, the present experimental study and related data reduction as well as strength and free-edge stress analysis of composite laminates indicate that the mechanical properties of PMC laminates have large variations due to various random factors involved in composite microstructures and processing. It is necessary to introduce the concept of probabilistic strength distribution in order to more accurately describe the mechanical behavior of PMC laminates and to predict their reliable performance and lifetime when integrated in composite structures.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Stress-Function Variational Method for Free-Edge Stress Analysis of Composite Laminates
Wu [14,45] refined and extended Yin's stress-function variational method [79,80] for high-efficiency, accurate free-edge stress analysis of general angle-ply composite laminates made of an arbitrary number of UD composite plies. Consider a composite laminate with traction-free edges, as shown in Figure A1. Assume the laminate to be made of laminae with uniform thickness. Each ply (lamina) may be oriented at an arbitrary angle in the x-y plane. In the laminate coordinate system, the constitutive equation for each individual ply can be expressed in terms of the generalized Hooke's law in a contracted notation [1,116,117] as where S ij (i,j = 1,2,3 . . . ,6) are elements of the transformed compliance matrix, α i (i = 1,2,3) are the transformed thermal expansion coefficients, ∆T is the temperature change from a reference, stress-free condition, and subscript indices (1, 2, 3) stand for (x, y, z). Furthermore, the present experimental study and related data reduction as well as strength and free-edge stress analysis of composite laminates indicate that the mechanical properties of PMC laminates have large variations due to various random factors involved in composite microstructures and processing. It is necessary to introduce the concept of probabilistic strength distribution in order to more accurately describe the mechanical behavior of PMC laminates and to predict their reliable performance and lifetime when integrated in composite structures.
Author Contributions: X-F.W. conceived and performed the experimental and numerical studies. X-F.W. and O.Z. analyzed the data and wrote the paper.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Stress-Function Variational Method for Free-Edge Stress Analysis of Composite Laminates [14,60]
Wu [14,45] refined and extended Yin's stress-function variational method [79,80] for highefficiency, accurate free-edge stress analysis of general angle-ply composite laminates made of an arbitrary number of UD composite plies. Consider a composite laminate with traction-free edges, as shown in Figure A1. Assume the laminate to be made of laminae with uniform thickness. Each ply (lamina) may be oriented at an arbitrary angle in the x-y plane. In the laminate coordinate system, the constitutive equation for each individual ply can be expressed in terms of the generalized Hooke's law in a contracted notation [1,116,117] where Sij (i,j = 1,2,3…,6) are elements of the transformed compliance matrix, αi (I = 1,2,3) are the transformed thermal expansion coefficients, ∆T is the temperature change from a reference, stressfree condition, and subscript indices (1, 2, 3) stand for (x, y, z). Figure A1. Geometry and coordinate system for a laminate free-edge problem. Figure A1. Geometry and coordinate system for a laminate free-edge problem.
In Equation (A1), the following contracted notation is adopted and S ij is the inverse transformed stiffness matrix in the (x, y, z)-coordinate system as Here, introduce the transform of a ply stiffness matrix, which is used for stiffness transform between the ply-material coordinate system and the global coordinate system of the laminate. Consider a rotation of the original (x, y, z)-coordinate system along the z-axis to the current coordinate system (x / , y / , z / ), as shown in Figure A2. In Equation (A1), the following contracted notation is adopted 1  1 1  1  1 1   2  2 2  2  2 2   3  3 3  3  3 3   4  2 3  4  2 3   5  3 2  5  3 1   6  1 2  6  1 2   ,  , , and Sij is the inverse transformed stiffness matrix in the (x, y, z)-coordinate system as ij ij Here, introduce the transform of a ply stiffness matrix, which is used for stiffness transform between the ply-material coordinate system and the global coordinate system of the laminate. Consider a rotation of the original (x, y, z)-coordinate system along the z-axis to the current coordinate system (x / , y / , z / ), as shown in Figure A2.
. Figure A2. Rotation of the coordinate system along z-direction.
The nonzero stiffness elements C / ij in the new coordinate system can be expressed in terms of C 11 in the original coordinate system as [116] C / 11 = m 4 C 11 + 2m 2 n 2 (C 12 + 2C 66 ) + 4mn(m 2 C 16 + n 2 C 26 ) + n 4 C 22 , Traction-free boundary conditions of the laminate described in Figure A1 can be expressed as The laminate is assumed to be in the state of plain strain in x-direction, which is characterized by specifying a constant strain value of ε 1 . Thus, the corresponding axial stress σ 1 can be extracted from relation (A1) as where the tensor summation rule is assumed. The constant strain value of ε 1 can be reflected in the generalized Hooke's law (A1) as Since the laminate is in the plane strain state, two Lekhnitskii's stress potentials F(y, z) and Ψ(y, z) [118] are introduced as In the i-th layer, a nondimensional thickness coordinate η is defined as η = (z − z i−1 )/(z i − z i−1 ). In the classic plate theory, the in-plane stresses vary linearly across the plate thickness, and this feature is approximately valid in each thin layer of a composite laminate. Thus, the stress function F in each ply can be approached by a polynomial function of degree three in η-direction (i.e., the thickness direction), while the stress function Ψ in each ply is approximated by another polynomial function of degree two in η-direction (i.e., the laminate width direction). In the i-th ply, the stress functions are assumed as t 2 F i (y) and tΨ i (y) on the interface z = z i and t 2 F i−1 (y) and tΨ i−1 (y) on z = z i−1 , where t is the ply thickness. Thus, from (A9), the cubic polynomial approximation of the in-plane stress function in the i-th ply must have the following expressions where G i (y) and G i−1 (y) are the derivatives of F i (y) and F i−1 (y), and F 0 (y) = G 0 (y) = F n+1 (y) = G n+1 (y) = 0 are related to the traction-free conditions on the lower and upper surfaces of the composite laminate. The quadric polynomial approximation of the out-of-plane function in the i-th ply can be expressed as where H i (y) is the derivative of G i (y) on the upper side of the same interface, which is generally discontinuous across the interface, and Ψ 0 (y) = Ψ n+1 (y) = 0 relates to the traction-free conditions on the lower and upper surfaces of the composite laminate. The above two Lekhnitskii's stress potentials F(y, z) and Ψ(y, z) are obtained by modifying those introduced by Yin [94,95] (introducing thickness t and then ensuring F, G, ψ and H with the same units). These modifications lead to a well-conditioned generalized eigenvalue problem, which can be solved efficiently even for a quite large matrix rank.
The governing equations are achieved by taking the stationary value of the complementary strain energy of the laminate per unit length as By substituting (A10) and (A11) into (A12) and then integrating by parts and ply by ply, a system of variational equations is yielded as where and the global matrices W, V, and U are constant real-valued symmetric square matrices and b is a column relating the initial specified external and thermal strains, which are assembled from the element matrices W e , V e , U e , and b e in each ply as where matrices (A15), (A16), and (A18) relate to vectors: [F (4) i , F i−1 , G i , G i−1 , G i , G i−1 , Ψ i , Ψ i−1 , H i ] T , and [F (2) i , F i−1 , G i , G i−1 , Ψ i , Ψ i−1 , H It should be mentioned that the global coefficient matrices W, V, and U are real-valued symmetric square matrices, and coefficient matrix W of the differential operator d 4 /dx 4 has the nontrivial elements only in a 2n × 2n square submatrix in the upper left corner. Equation (A13) is satisfied for arbitrary variations δH 0 , δF 1 , δG 1 , δΨ 1 , . . . , and δH n if and only if {Y} satisfies the following Euler-Lagrange equations (i.e., the weak solution) with the homogenous boundary conditions at the free edges 1, 2, . . . , n) at y= 0 and y= 2b. (A20) To solve Equation (A19), one can first solve the corresponding homogenous equation by determining its eigenvalues and related eigenvectors. Thus, by assuming {Y} = {X 0 }exp(λy/t), one can obtain the corresponding generalized eigenvalue problem as The particular solution to Equation (A19) can be determined as where {X 10 } = {F 1 , F 2 , . . . , F n , G 1 , G 2 , . . . , G n } T , and {X 20 } = {ψ 1 , ψ 2 , . . . , ψ n , H 0 , H 1 , . . . , H n } T . By introducing an auxiliary vector {X 1 } = λ 2 {X 10 } T , eigenvalue problem (A21) can be cast into a generalized eigenvalue problem as and I is a diagonal unit matrix, and 0 is a zero matrix. Generalized eigenvalue problem (A24) is well conditioned, and its eigenvalues and corresponding eigenvectors can be obtained using a standard mathematic program. The general eigenvalues and eigenvectors are complex-valued since matrices A and B are unsymmetrical. Numerical results show that a standard algorithm provided by MatLab is robust for regular laminates with the ply number n up to 70, which corresponds to a huge generalized unsymmetrical eigenvalue problem of the order 421. In the case of Equation (A21), assume the obtained eigenvalues are λ i (i = 1, 2 . . . 6n + 1) and the corresponding eigenvectors are Φ i (i = 1, 2, . . . , 6n + 1). Thus, the general solution to Equation (A19) can be obtained from a linear combination of the eigenvectors as [α (1) i exp(λ i y/t) + α (2) i exp(−λ i y/t)]Φ i + Y 0 , (y ≥ 0) (A26) where (12n + 2) unknowns α (1) i , α i (i= 1, 2, . . . , 6n+1) may be determined by (12n + 2) traction-free conditions as described in (A20). Because of the rapidly decaying exponential functions, in the calculations only a semi-infinite laminate with one free-edge needs to be considered. Thus, the general solution (A26) can be further simplified as where λ i (i = 1, 2, . . . , 6n + 1) are the eigenvalues with negative real parts, and the coefficients α i (i = 1, 2, . . . , 6n + 1) are determined by (6n + 1) traction-free conditions at y = 0 as described in (A20). Finally, the stress components can be extracted by substituting (A27) and (A14) into (A9).