Micromechanics of Cracked Laminates under Uniaxial Load: A Comparison between Approaches

This paper compares stiffness degradation models of cross-ply glass fiber/epoxy laminates based on four of the most commonly used approaches to micromechanical modelling: shear-lag, variational, McCartney, and synergistic damage mechanics (SDM). All of these include the process of defining [0/90]s laminate unit cell, from which governing differential equations and corresponding boundary conditions are stated. Afterwards, these boundary value problems (BVP) are solved in order to obtain a stress function which couples the initial and perturbation stresses, the latter being in function of crack density, thus related to material stiffness reduction. When compared against experimental results, shear-lag model presented accurate results however, additional differentiation and integration steps were required in order to obtain the final stress field. Hashin’s variational method predicts correctly the boundary conditions at crack surfaces and gives out the complete stress field. McCartney’s approach shows further improvement over the previous two models, taking into account thermal strains and stresses. Finally, SDM, which is designed for numerical experimentation, implying a more economical alternative in comparison to traditional physical experimentation, also presented very good agreement with experimental results and can be extended to arbitrary laminate stackings, going beyond the classical cross-ply.


Introduction
Composite laminates have found increased use in mechanical applications over the last years, especially in areas related to aerospace, civil, and mechanical systems.Efforts have been particularly intense on damage assessment of large composite structures.There are many models such as the ones developed by Shokrieh [1] that focus on the degradation of mechanical properties, relying on experimental curve fits which do not necessarily have a direct physical interpretation.This can represent a large obstacle for understanding the mechanisms of onset and propagation of damage over a component or structure.
Of these damage modes, matrix cracking has received the most attention, followed by delamination.After the material reaches certain level of crack density, there is progressively less space for transverse crack formation.This phenomenon is known as crack saturation [2] and the corresponding density is known as Characteristic Damage State (CDS), but before this happens, delamination may occur.Delamination [3] is the mechanism of separation between plies caused by an interfacial weakening partially due to previously formed intralaminar cracks.Delamination proves to be more catastrophic than matrix cracking, because it could render the structure useless in a short amount of time driving it to failure [4].
Matrix cracking in 90-degree plies was first studied by Reifsneider [5], who experimentally investigated this phenomenon on [0/±45/90] s laminates subjected to uniaxial quasistatic and cyclic tests.Although practical, physical experiments can become relatively expensive when less assumption is considered, especially when compared with ever-increasing computer capacity [6].
The first micromechanical model was shear-lag, introduced by Garrett and Bailey [7][8][9] and significantly improved by Cox [10].Shear-lag model relies on the basic principle that The crack is supposed to span all the laminate width, adapted from [12,14].
relaxed stresses arising from the formation of intralaminar cracks transfer to ply interface [11] in form of shear stress, increasing in turn the occurrence of delamination.In the present review, Berthelot complete parabolic shear-lag model is used [12], which means a parabolic variation of axial displacement is introduced for both 0-and 90-degree layers (see (1a)-(1b)).Then, corresponding stress functions are obtained, which ultimately translate into stress accumulation and degradation mechanical properties, particularly stiffness [13].
Another family of micromechanical approaches is proposed by Hashin [14] who proposes a variational stressperturbation function, departing from crack-free stress state.In order to achieve this, Hashin applied the minimum complementary energy principle and optimized it.Subsequently, the total stress field (crack-free plus perturbation) was recovered for cross-ply (0/90 m ) s laminates.Similar to the first approach, the stresses formed cracks, which in turn translated in material degradation, chiefly stiffness reduction.Hashin's variational analysis has become a seminal work in the field of micromechanical damage analysis on composite materials [14].For example, it has been used and adapted for predicting fatigue [15], random crack growth [16], and delamination [17].Thereafter, McCartney [18] introduced his model based on Hashin's variational approach and classical elasticity, presenting three novelties: inclusion of the thermal expansion strains, explicit expressions for displacements, and a further expansion to three-dimensional stress evaluation that will not be reviewed in this work.
Continuum damage mechanics (CDM) is another family that has become more commonplace in literature over the last years due to increased computational capacity.These models are based on updating stiffness matrices in function of crack opening displacement (COD) or density.Their scheme may be iterative, such as the so-called self-consistent methods [19,20] and the proposed by , or departing from an initial set of real-life [24] or virtual finite element simulations such as the Synergistic Damage Mechanics (SDM) model presented by Singh and Talreja [25][26][27], which expanded the scope of micromechanical modelling to laminates other than cross-ply.
Driven by aerospace and wind energy applications in which other laminate configurations such as the [0/±45/90] s are preferred, recent models such as the ones proposed by Hajikazemi [28] and Adumitroaie [29] have also expanded to laminates other than cross-ply.However, as these models are usually extensions of previous models (Hashin's variational for the case of Hajikazemi's model and SDM for Adumitroaie's model), the main goal of this work is drawing a comparison between these four approaches (shear-lag, Hashin's variational, McCartney, and SDM) and highlighting the advantages and setbacks for each one.This comparison will be helpful for selecting the most appropriate matrix cracking approach for future works related to damage propagation in fiber-reinforced polymers for different applications.The nomenclature used throughout the article is summarized in the Nomenclature.

Models
2.1.Shear-Lag Model.As mentioned above, the shear-lag models were the first approaches on intralaminar cracking.The basic governing principle consists in stress redistribution along the loading direction causing an increase in interfacial shear-load due to crack onset.The basic unit cell is shown in Figure 1, representing a symmetric (0/90) s laminate element composed of two 0 ∘ -degree outer plies of thickness  0 and a 90 ∘ -degree inner ply of thickness  90 .Total laminate thickness is thus ℎ = 2( 0 +  90 ).In addition,  is loading direction, 2 is crack spacing, and ply interface is located along coordinate  =  90 .
Shear-lag model is considerably simple compared to the other approaches presented below because of the following assumptions: (a) plane stress and (b) longitudinal stresses   are assumed to be uniform across the thickness of 0 ∘ and 90 ∘ plies and (c) the tensile load applied to 0 ∘ plies is considered to be transmitted into the 90 ∘ ply by shear of the transverse ply [12].
However, these assumptions may compromise the accuracy of the model.For this reason, Berthelot [12] introduced parabolic profiles for -direction displacements in 0-and 90degree plies.
This produces the boundary value problem (BVP) with governing equation in terms of average longitudinal stress in the 90-degree layer ( 90  ) where This governing equation is subject to the boundary conditions (BCs) shown in Table 1.
Thus, the complete solution for -direction stress field in 0-degree and 90-degree plies is given out by where (4e) expresses the aspect ratio .This parameter is of great importance as it reflects the relationship of crack spacing to ply thickness, having a profound impact on the degradation model.Meanwhile, (4g) shows the stacking parameter .Shear-lag method calculates in a straightforward fashion the interfacial shear-load that may be applied later to delamination models.The complete model for interlaminar shear stress is also obtained as follows: The shear stress field is obtained introducing (4c)-(4d) into (2) and solving for   .Subsequently the -direction axial stress field is obtained introducing this solution into equilibrium equation ( 6) and solving for   .This procedure is repeated for the 0-degree ply.
Finally, these stresses have an effect on the material, which is reflected in degradation rules such as the stiffness reduction formula obtained by Ogin [13]: 2.2.Variational Method.In 1985, Hashin proposed another approach to solve the crack problem in a cross-ply laminate, already shown in Figure 1.The variational method is an energetic approach that consists in posing an internal energy formulation for a perturbation stress, which is considered a variation from the intact material stress, caused by crack onset.Because this is an equilibrium situation, total energy change is considered near zero, with this method being therefore soundly based on virtual work statics.
The assumptions taken for this approach are as follows: (i) The normal stress in external load direction is constant over ply thickness.
(ii) Shear stresses develop only within a boundary layer of unknown thickness in between plies.
(iii) Cracks remain sufficiently far apart so that their mutual interaction can be neglected.
Said this, the stresses in the cracked laminate are where ,  subscripts stand for ,  direction of the stresses.The index  represents the ply, whether 0-or 90-degree.
Finally, the subscript 0 represents the stress of the intact matrix and , the perturbation stress generated by the crack.After applying the equilibrium equations obtained from basic elasticity over the region illustrated in Figure 1 and the boundary conditions listed in Table 1, the perturbation stress field in terms of intact 90-degree ply stress ( 90 ) is found out.
where the perturbation function () is evaluated by the following formulation, provided that in most epoxy/carbon materials, 4 >  2 : where  1 =  ) .

(11c)
For the equations governing case 4 <  2 , the reader is referred to Talreja [11].As in the previous method, a stress  1. Adapted from [18].field is recovered and then applied in a degradation rule that has an effect on material's elasticity.
where  1 =  90 /  and is obtained from traditional laminate analysis.The other functions are evaluated as follows: Finally, it is important to take into account the fact that this formulation is only applicable to equally spaced cracks.If the separation is arbitrary, a probabilistic distribution such as those used by Vinogradov [16] must be used.

McCartney Method.
McCartney proposed another analytical method in 1992 [18].Similar to variational approach, the method models the same phenomenon by optimizing an energetic expression.As in the previous two models, it only considers cross-ply laminates.Nevertheless, there are three noteworthy novelties.First of all, components of thermal expansion are taken into account for stain and stress modelling.
It is important to keep in mind that (13) actually represents two sets of equations, one for outer 0-degree plies and another one for inner 90-degree plies.Engineering elastic constants that feature an apostrophe are the modified moduli, which are shown in more details in [18].The original moduli appear with no apostrophe.Finally, the development for the uniform transverse strain in the cracked laminate  *  is also shown in the referred paper.Assuming perfect interlaminar bonding, McCartney model is governed by the boundary conditions shown in Table 1 (cf.Figure 2(b)), where   is the longitudinal strain experienced by the cracked composite and is also explained in further details in [18].Now, assume that longitudinal stress components  0  and  90  are independent of  and are of the form where  and  * are the corresponding longitudinal and transverse strains for the undamaged laminate.We now proceed to introduce (14) in field (1a) and (1b) with boundary conditions (Table 1) in order to obtain the following stress field: It is interesting to note the similarities between (10a), (10b), (10c), and (10d) from Hashin and ( 14) and ( 15) from McCartney.Now, the corresponding displacement field (16a), (16b), (16c), and (16d) is obtained where () is the -displacement at the interface of inner and outer plies And perturbation function (), similar to Hashin's, is defined for most composites as Figure 3: (a) Representative unit cell for a laminate under SDM micromechanics analysis, adapted from [25].(b) Special case for cross-ply analysis, adapted from [24].
Here , , and  are constants dependent on geometry and material and are defined thoroughly in [18].Then, the function () is computed from the following integral: Once the perturbation function is obtained it is used to calculate stiffness reduction in the cracked laminate.
where ] 0  ,  0  , and  0  are the longitudinal Poisson ratio and longitudinal and transverse Young's moduli of the undamaged laminate, respectively.Finally, in the same work, McCartney introduces a novel 3D version of the formerly explained model, which is not included in this analysis for the sake of comparison, and because the plane strain assumption under which the 2D model operates is enough for adequately modelling the mechanics of cracked composite plies.

Synergistic Damage Mechanics (SDM).
Up to now, micromechanical models shown in Sections 2.1 to 2.3 refer to cracked composites limited to the cross-ply case, which has very limited applications.For most aerospace and wind energy applications, a [0/±45/90] s configuration is more common.For this reason, Singh and Talreja developed the Synergistic Damage Mechanics (SDM) model [25], which is explained as follows.
The basic cell under which the SDM model operates is shown in Figure 3(a).The picture clearly shows that the model is flexible for stacking more arbitrary ply orientations than those presented previously.For the sake of comparison against previous shown models, the basic cell used for simulations shown in Section 3 is also cross-ply [24], Figure 3(b); nevertheless, in this section the general case is described in order to discuss all the innovations the SDM model poses.
Of these innovations, the most important over the previous micromechanical models is that SDM can initiate and multiply cracks in any of the layers of the laminate, even if the probability of this occurring at 0 ∘ is relatively low.Each of these crack orientations is defined in SDM as damage modes, denoted by .As it will be described, each of these modes will have their particular stress fields, stiffness changes, and damage tensors as each one of them contributes individually to the total damage state and stiffness reduction in the laminate.
The model starts with a damage state tensor D  ij , with a given spacing (  ), thickness (  ), restriction parameter (  ), and orientation     for the damage mode .Finally, the total laminate thickness is represented by .
The normal orientation vector is expressed by   = (sin , cos , 0).Meanwhile, the restriction parameter   accounts for the constraining effect on ply cracks caused by adjacent plies in the laminate.
where  22 is the transformed applied strain component normal to the crack surface and ũ  is the normal average crack opening displacement (COD) for a given crack spacing.These averaged CODs are computed previously from a series of FE simulations using software such as ANSYS [30].The micromechanical unit cell (Figure 3) is modelled by replicating the material stacking and elastic constants, defined by the size of crack spacing, and three boundary conditions [25] shown in Table 1.
Another important quantity resulting from the model is the stiffness tensor and its reduction caused by the presence of the microcracks.For the case being [25], the damage stiffness tensors for each mode    are expressed by This stiffness reduction scheme needs the determination of material constants    .In order to obtain these constants, stiffness reduction curves are taken from previously described ANSYS simulations over [0/90] s laminates where the values of  11 ,  22 ,  12 , and V 12 are plotted against the value of crack spacing   for each ply.Finally, a curve fitting procedure is conducted in order to satisfy the following set of equations: where the 0-superscripted engineering constants correspond to those of the virgin laminate.For a [0/±/90] s laminate, the total damage tensor is given by Subsequently, the resulting stress tensor is recalculated As stated previously, even if Synergistic Damage Mechanics (SDM) model has the capability of analyzing stresses and stiffness degradation in laminates such as [0/±45/90] s , it is first necessary to perform simulations on cross-ply laminates.
As in all other stiffness degradation models, SDM stiffness reduction for [0/90] s glass/epoxy laminates is plotted against crack density, which is also in terms of the reciprocal of crack spacing as per (27).This means that for each case of crack density (  = 16, 8, 4, 2, 1, 0.75, 0.5, 0.2, and 0.1 mm) FE models were prepared and analyzed using ANSYS by applying loading, geometry, and boundary conditions as shown in Figure 3(b).An intact case simulated with highly spaced cracking (  = 32) was also prepared.
Boundary conditions for each crack spacing case remained the same, with applied displacement at right face of () =2 = 5 m.Left face is fixed and there is a symmetry condition in the bottom face, along -axis.Cracks were supposed to have grown throughout the width of the cell.Subsequently, each model was meshed using 10,000-25,000 3D tetrahedral elements.Then, stress field was recovered and averaged along the volume.At last, longitudinal Young's modulus is then obtained, according to [31], by applying the following equation: where  RVE is the volume of the whole unit cell, ∫     ∼ ∑  ( ,   ), where  , and   are the normal axial stress and volume of element ;  0 is the applied displacement in the right face, as well.Finally, these curves can later be used for obtaining through fitting, materials constants    , which can later be used in (24a) and (24b) in order to calculate stiffness degradation in stacking cases such as [0/45/90] s .

Results and Discussion
From the formulation explained in Section 2, stress distribution is drawn for longitudinal stress   in 90-degree layer for an -glass fiber/epoxy laminate with material constants (Table 2) obtained from Highsmith and Reifsneider [5], with 90-ply thickness  90 = 0.203 mm, stacking parameter  = 1, and aspect ratio of cracking  = 2.5 and  1 = 0.464.Cracking aspect ratio and 90-ply thickness will be kept constant through results analysis; for other values their variation will be presented as needed.
It is seen from Figure 4 that minimum normal stress is zero and is found where the cracks are located; meanwhile the maximum longitudinal stress is found midway between consecutive cracks.Furthermore, it is observed from Figure 5 that crack onset has increased interlaminar shear.
Nevertheless, shear-lag method is at disadvantage when compared to variational micromechanics because the shearlag assumption forces the model not to follow the zerotraction boundary condition, in stark comparison to Hashin formulation.This brings considerable error for the shear-lag method as it approaches the crack surface (Figure 6).
Taking the shear stress variation for the interface located between 0-and 90-degree plies from Hashin variational method (Figure 7(a)) [14], it is observed from Figure 7 that shear stress redistributes to somewhere in the middle between crack surface and the -axis of symmetry because of the restriction imposed on crack surface.Again, this contrasts with Berthelot shear-lag model, which puts the maximum location of shear redistribution just above the cracking surface (Figure 5).Finally, a complete stress distribution field for  90  by Hashin method is illustrated in Figure 7(b), from where it is observed that  90  becomes progressively higher as it approaches the interlaminar interface.

Figure 2 :
Figure 2: (a) Unit cell used to derive the equations of the McCartney model.(b) Coordinates for boundary conditions shown in Table1.Adapted from[18].

Table 1 :
Comparison of boundary conditions of the four models analyzed in the present article.