Abstract
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 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.
1. 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] 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–9] and significantly improved by Cox [10]. Shear-lag model relies on the basic principle that 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 stress-perturbation 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 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 Barbero [21–23], 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–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] 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.
2. 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 laminate element composed of two 0°-degree outer plies of thickness and a 90°-degree inner ply of thickness . Total laminate thickness is thus . In addition, is loading direction, is crack spacing, and ply interface is located along coordinate .
(a)
(b)
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 90-degree plies.where , , and are integration constants. These profiles are then introduced within the framework of classic equations for elasticity: first kinematics (strain-displacement), then compliance, and finally equilibrium within the unit cell shown in Figure 1(b) [12]. This produces the boundary value problem (BVP) with governing equation in terms of average longitudinal stress in the 90-degree layer ()where , and the laminate stiffness constants for shear and Young’s modulus , subject to a global load .
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 bywhere (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 () is found out.where the perturbation function is evaluated by the following formulation, provided that in most epoxy/carbon materials, :where , , and . The material constants and may be evaluated as follows:For the equations governing case , the reader is referred to Talreja [11]. As in the previous method, a stress field is recovered and then applied in a degradation rule that has an effect on material’s elasticity. where 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.
2.3. 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. Then, even though precision is similar to Hashin, McCartney model includes explicit expressions for displacement field. Finally, in contrast to shear-lag, the McCartney model satisfies exactly a more comprehensive set of boundary conditions. The basic unit cell of the laminate between two cracks, subjected to uniaxial stress, is shown in Figure 2(a).
(a)
(b)
Then, from basic elasticity and including thermal strains, the field equations for generalized plane strain conditions are expressed by 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 and are independent of and are of the formwhere 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 obtainedwhere is the -displacement at the interface of inner and outer pliesAnd perturbation function , similar to Hashin’s, is defined for most composites asHere , , 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 , , and 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.
2.4. 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] 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.
(a)
(b)
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 , 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 . Meanwhile, the restriction parameter accounts for the constraining effect on ply cracks caused by adjacent plies in the laminate. where 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 laminates where the values of , , , and 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] laminate, the total damage tensor is given byOnce the model is ready, (24a) and (24b) is applied again, but now for the current crack separation parameter present in the laminate. The final results are the degraded elasticity engineering constants of the material. Subsequently, the resulting stress tensor is recalculatedAs 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], it is first necessary to perform simulations on cross-ply laminates. As in all other stiffness degradation models, SDM stiffness reduction for 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 ( was also prepared.Boundary conditions for each crack spacing case remained the same, with applied displacement at right face of = 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 is the volume of the whole unit cell, , where and are the normal axial stress and volume of element ; 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 .
3. 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 = 0.203 mm, stacking parameter , and aspect ratio of cracking = 2.5 and = 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 shear-lag assumption forces the model not to follow the zero-traction 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 by Hashin method is illustrated in Figure 7(b), from where it is observed that becomes progressively higher as it approaches the interlaminar interface.
(a)
(b)
Afterwards, degradation of longitudinal Young’s modulus against crack density , (27) for glass/epoxy laminates (), is plotted for all four stiffness reduction models (Figure 8). For shear-lag model used by Berthelot [12], Ogin model stiffness reduction curve is applied [13]. Meanwhile, for the Hashin variational method [14], Reifsneider formulation was employed instead [5]. McCartney [18] and SDM [25] incorporate their own stiffness reduction methods. As a first verification step, convergence for all four curves is verified against (29) proposed by Reifsneider [5]; there is very good agreement between all model results.Regarding experimental results, a comparison was drawn between all four models and Smith and Wood experimentation [32]. Again, from Figure 8 it is seen that only shear-lag and SDM follow good agreement at initial cracking stages. Smith and Wood do not provide experimental results for advanced ( cracks/mm), nonetheless convergence has been proved with (29).
4. Conclusions
A comparison between four micromechanical models has been presented. From degradation rules, it is noticed that all four methods covered in this paper predicted convergence value quite accurately. From Hashin’s variational method it is important to remark that it not only predicts correctly the boundary conditions at the surfaces of the crack but also readily gives out the complete stress field. Although limited by model assumptions, there are novel methods capable of analyzing arbitrary stacking [28] based on Hashin’s approach. In contrast, shear-lag degradation model is more accurate when compared against experimental results [32] but it takes additional differentiation and integration steps in order to obtain the final stress field. However, the assumptions under which it is governed limit its application to cross-ply laminates.
McCartney’s approach shows further improvement over the last two methods because it takes the generalized plane strain assumption, which soundly resonates with laminate mechanics. This model is also unique because it takes into account thermal strains and stresses. Furthermore, this model is applied to multiaxial loading conditions. By last, McCartney provides explicit expressions for computing displacements throughout the laminate cell. However, McCartney model is computationally expensive because of higher order differentiation involved in the computation of stiffness reduction.
For the case of SDM, the model not only converges but also replicates experimental results with relative accuracy. Furthermore, SDM is designed for numerical experimentation, implying a more economical alternative in comparison to traditional physical experimentation. Damage mechanics models present the advantage of being flexible because they may be applied to a wider range of laminate stacking, going beyond the classical cross-ply case and thus proving more useful in practical cases that involve laminates of arbitrary stacking [23, 25], such as the widely used [0/±45/90] laminates. However, these modelling approaches may be computationally expensive and heavily dependent on mesh size for the case of virtual simulations.
Nomenclature
| : | Interlaminar region displacement function (only McCartney model) |
| : | Integration constants (variational method) |
| : | Integration constants for parabolic displacement profile in 0-degree plies (shear-lag) |
| : | Integration constants for parabolic displacement profile in 90-degree plies (shear-lag) |
| : | McCartney perturbation function |
| : | Variational method material constants |
| : | Damaged stiffness tensor |
| : | Undamaged stiffness tensor |
| : | SDM damage tensor |
| : | Ply longitudinal Young’s modulus |
| : | Ply transverse Young’s modulus |
| : | Intact laminate longitudinal Young’s modulus |
| : | Intact laminate transverse Young’s modulus |
| : | Cracked laminate longitudinal Young’s modulus |
| : | McCartney material constants |
| : | Ply in-plane shear modulus |
| : | Ply out-of-plane shear modulus |
| : | McCartney integration constants |
| : | Volume of unit cell |
| : | Crack aspect ratio |
| : | SDM material constants for ply |
| : | Crack density |
| : | Half-laminate thickness |
| : | Stress ratio in cross-ply laminates, obtained from classical lamination theory (CLT) |
| : | Half crack spacing |
| : | Normal orientation vectors |
| : | Variational method material constants |
| : | McCartney integration constants |
| : | Crack spacing, that is, |
| : | Ply thickness |
| : | Displacements |
| : | Applied displacement |
| : | Normal average crack opening displacement |
| : | Positions |
| : | Temperature increase |
| : | Variational method integration constant |
| : | McCartney integration constant |
| : | McCartney material degradation constants |
| : | Laminate stacking parameter |
| : | Variational method integration constants |
| : | Longitudinal and transverse thermal expansion coefficients |
| : | Longitudinal and transverse strains in undamaged laminates |
| : | Longitudinal and transverse strains in cracked laminates |
| : | Normal and shear strain |
| : | Shear-lag integration constant |
| : | Ply orientation respect loading axis |
| : | SDM interlaminar restriction parameter |
| : | Shear-lag integration constant |
| : | In-plane ply Poisson’s ratio |
| : | Out-of-plane ply Poisson’s ratio |
| : | In-plane intact laminate Poisson’s ration |
| : | Normalized -position coordinate |
| : | Stress in intact 0- and 90-degree plies |
| : | Longitudinal load applied on laminate |
| : | Normal stress |
| : | Shear stress |
| : | Average longitudinal stress |
| : | Interlaminar shear stress function |
| : | Variation method perturbation function. |
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.