Damage and failure of a laminated composite plate

A composite plate with a hole in the center is subjected to in-plane compression. The plate is made of 24 plies of T300/976 graphite-epoxy in a [(−45/+45)${}_6$]${}_s$ layup. Each ply has a thickness of 0.1429 mm (0.005625 in); thus, the total plate thickness is 3.429 mm (0.135 in). The plate has a length of 101.6 mm (4.0 in) and a width of 25.4 mm (1.0 in), and the diameter of the hole is 6.35 mm (0.25 in). The plate is loaded in compression in the length direction. The thickness of the plate is sufficient that out-of-plane displacements of the plate can be ignored. The compressive load is measured, as well as the length change between two points, originally a distance of 25.4 mm (1.0 in) apart, above and below the hole. The plate geometry is shown in Figure 1.

The material behavior of each ply is described in detail by Chang and Lessard. The initial elastic ply properties are longitudinal modulus $E_x$=156512 MPa (22700 ksi), transverse modulus $E_y$=12962 MPa (1880 ksi), shear modulus $G_{xy}$=6964 MPa (1010 ksi), and Poisson's ratio ${\nu }_{xy}$=0.23. The material accumulates damage in shear, leading to a nonlinear stress-strain relation of the form

where $G_{xy}$ is the (initial) ply shear modulus and the nonlinearity is characterized by the factor $\alpha$=2.44×10⁻⁸ MPa⁻³ (0.8×1⁻⁵ ksi⁻³).

Failure modes in laminated composites are strongly dependent on geometry, loading direction, and ply orientation. Typically, one distinguishes in-plane failure modes and transverse failure modes (associated with interlaminar shear or peel stress). Since this composite is loaded in-plane, only in-plane failure modes need to be considered, which can be done for each ply individually. For a unidirectional ply as used here, five failure modes can be considered: matrix tensile cracking, matrix compression, fiber breakage, fiber matrix shearing, and fiber buckling. All the mechanisms, with the exception of fiber breakage, can cause compression failure in laminated composites.

The failure strength in laminates also depends on the ply layup. The effective failure strength of the layup is at a maximum if neighboring plies are orthogonal to each other. The effective strength decreases as the angle between plies decreases and is at a minimum if plies have the same direction. (This is called a ply cluster.) Chang and Lessard have obtained some empirical formulas for the effective transverse tensile strength; however, in this model we ignore such effects. Instead, we use the following strength properties for the T300/976 laminate: transverse tensile strength $Y_t$=102.4 MPa (14.86 ksi), ply shear strength $S_c$=106.9 MPa (15.5 ksi), matrix compressive strength $Y_c$=253.0 MPa (36.7 ksi), and fiber buckling strength $X_c$=2707.6 MPa (392.7 ksi).

The strength parameters can be combined into failure criteria for multiaxial loading. Four different failure modes are considered in the model analyzed here.

• Matrix tensile cracking can result from a combination of transverse tensile stress, ${\sigma }_y$, and shear stress, ${\sigma }_{xy}$. The failure index, $e_m$, can be defined in terms of these stresses and the strength parameters, $Y_t$ and $S_c$. When the index exceeds 1.0, failure is assumed to occur. Without nonlinear material behavior, the failure index has the simple form,

  $$e_m^2={\left(\frac{{\sigma }_y}{Y_t}\right)}^2+{\left(\frac{{\sigma }_{xy}}{S_c}\right)}^2.$$

  With nonlinear shear behavior taken into consideration, the failure index takes the more complex form,

  $$e_m^2={\left(\frac{{\sigma }_y}{Y_t}\right)}^2+\frac{2{\sigma }_{xy}^2/G_{xy}+3\alpha {\sigma }_{xy}^4}{2S_c^2/G_{xy}+3\alpha S_c^4}.$$

• Matrix compressive failure results from a combination of transverse compressive stress and shear stress. The failure criterion has the same form as that for matrix tensile cracking:

  $$e_m^2={\left(\frac{{\sigma }_y}{Y_c}\right)}^2+\frac{2{\sigma }_{xy}^2/G_{xy}+3\alpha {\sigma }_{xy}^4}{2S_c^2/G_{xy}+3\alpha S_c^4}.$$

  The same failure index is used since the previous two failure mechanisms cannot occur simultaneously at the same point. After the failure index exceeds 1.0, both the transverse stiffness and Poisson's ratio of the ply drop to zero.

• Fiber-matrix shearing failure results from a combination of fiber compression and matrix shearing. The failure criterion has essentially the same form as the other two criteria:

  $$e_{fs}^2={\left(\frac{{\sigma }_x}{X_c}\right)}^2+\frac{2{\sigma }_{xy}^2/G_{xy}+3\alpha {\sigma }_{xy}^4}{2S_c^2/G_{xy}+3\alpha S_c^4}.$$

  This mechanism can occur simultaneously with the other two criteria; hence, a separate failure index is used. Shear stresses are no longer supported after the failure index exceeds 1.0, but direct stresses in the fiber and transverse directions continue to be supported.

• Fiber buckling failure occurs when the maximum compressive stress in the fiber direction ($-{\sigma }_x$) exceeds the fiber buckling strength, $X_c$, independent of the other stress components:

  $$e_b=-\frac{{\sigma }_x}{X_c}.$$

  It is obvious that, unless the shear stress vanishes exactly, fiber-matrix shearing failure occurs prior to fiber buckling. However, fiber buckling may follow subsequent to fiber shearing because only the shear stiffness degrades after fiber-matrix shearing failure. Fiber buckling in a layer is a catastrophic mode of failure. Hence, after this failure index exceeds 1.0, it is assumed that the material at this point can no longer support any loads.

In this example the primary loading mode is shear. Therefore, failure of the plate occurs well before the fiber stresses can develop to a level where fiber buckling takes place, and this failure mode need not be taken into consideration.

Chang and Lessard assume that after failure occurs, the stresses in the failed directions drop to zero immediately, which corresponds to brittle failure with no energy absorption. This kind of failure model usually leads to immediate, unstable failure of the composite. This assumption is not very realistic: in reality, the stress-carrying capacity degrades gradually with increasing strain after failure occurs. Hence, the behavior of the composite after onset of failure is not likely to be captured well by this model. Moreover, the instantaneous loss of stress-carrying capacity also makes the postfailure analysis results strongly dependent on the refinement of the finite element mesh and the finite element type used.

To simulate the shear nonlinearity and the failure modes (matrix failure in tension or compression and fiber-matrix shear failure), the elastic properties are made linearly dependent on three field variables. The first field variable represents the matrix failure index, the second represents the fiber-matrix shear failure index, and the third represents the shear nonlinearity (damage) prior to failure. The dependence of the elastic material properties on the field variables is shown in Table 1.

To account for the nonlinearity, the nonlinear stress-strain relation must be expressed in a different form: the stress at the end of the increment must be given as a linear function of the strain. The most obvious way to do this is to linearize the nonlinear term, leading to the relation

where i represents the increment number. This relation can be written in inverted form as

thus providing an algorithm to define the effective shear modulus.

However, this algorithm is not very suitable because it is unstable at higher strain levels, which is readily demonstrated by stability analysis. Consider an increment where the strain does not change; i.e., ${\gamma }_{xy}^{\left(i+1\right)}={\gamma }_{xy}^{\left(i\right)}={\gamma }_{xy}.$ Let the stress at increment i have a small perturbation from $\sigma _{xy}^e{}^{\left(i\right)}$, the exact solution at that increment: ${\sigma }_{xy}^{\left(i\right)}=\sigma _{xy}^e{}^{\left(i\right)}+\mathrm{\Delta }{\sigma }_{xy}^{\left(i\right)}$. Similarly, at increment i+1, ${\sigma }_{xy}^{\left(i+1\right)}=\sigma _{xy}^e{}^{\left(i+1\right)}+\mathrm{\Delta }{\sigma }_{xy}^{\left(i+1\right)}$. For the algorithm to be stable, $\mathrm{\Delta }{\sigma }_{xy}^{\left(i+1\right)}$ should not be larger than $\mathrm{\Delta }{\sigma }_{xy}^{\left(i\right)}$. The perturbation in increment i+1 is calculated by substituting ${\sigma }_{xy}^{\left(i\right)}$ in the effective shear modulus equation and linearizing it about $\sigma _{xy}^e{}^{\left(i\right)}$:

where ${\sigma }_{xy}=\sigma _{xy}^e{}^{\left(i\right)}.$ The perturbation in increment i+1 is larger than the perturbation in increment i if

which, after elimination of ${\gamma }_{xy}$, reduces to the expression

Hence, instability occurs when the “nonlinear” part of the shear strain is larger than the “linear” part of the shear strain.

To obtain a more stable algorithm, we write the nonlinear stress-strain law in the form

where $\beta$ is an as yet unknown coefficient. In linearized form this leads to the update algorithm

or, in inverted form,

Following the same procedure as that for the original update algorithm, it is readily derived that a small perturbation, $\mathrm{\Delta }{\sigma }_{xy}^{\left(i\right)}$, in increment i reduces to zero in increment i+1 if $\beta =2\alpha$. Hence, the optimal algorithm appears to be

Finally, this relation is written in terms of the damage parameter d:

where

This relation is implemented in user subroutines USDFLD and VUSDFLD, and the value of the damage parameter is assigned directly to the third field variable used for definition of the elastic properties.

The failure indices are calculated with the expressions discussed earlier, based on the stresses at the start of the increment:

The values of the failure indices are not assigned directly to the field variables: instead, they are stored as solution-dependent state variables. Only if the value of a failure index exceeds 1.0 is the corresponding user-defined field variable set equal to 1.0. After the failure index has exceeded 1.0, the associated user-defined field variable continues to have the value 1.0 even though the stresses may reduce significantly, which ensures that the material does not “heal” after it has become damaged.

The plate consists of 24 plies of T300/976 graphite-epoxy in a [(−45/+45)${}_6$]${}_s$ layup. Instead of modeling each ply individually, we combine all plies in the −45° direction and all plies in the +45° direction. Consequently, only two layers need to be modeled separately:

1. A layer in the −45° direction with a thickness of 1.715 mm (0.0675 in).

2. A layer in the +45° direction with a thickness of 1.715 mm (0.0675 in).

The corresponding finite element model consists of two layers of CPS4 plane stress elements, with thicknesses and properties as previously discussed. The quarter-symmetry finite element model is shown in Figure 1.

The implementation of nonlinear material behavior with user-defined field variables is explicit: the nonlinearity is based on the state at the start of the increment. Hence, in Abaqus/Standard analyses the user must ensure that the time increments are sufficiently small, which is particularly important because the automatic time increment control in Abaqus/Standard is ineffective with the explicit nonlinearity implemented in USDFLD. If automatic time incrementation is used, the maximum time increment can be controlled from within subroutine USDFLD with the variable PNEWDT. This capability is useful if there are other nonlinearities that require automatic time incrementation. In this example the only significant nonlinearity is the result of the material behavior. Hence, fixed time incrementation can be used effectively. In Abaqus/Explicit analyses the stable time increment is usually sufficiently small to ensure good accuracy.