Improving failure sub-models in an orthotropic plasticity-based material model

Abstract

Theoretical details of two failure criteria implemented in an orthotropic plasticity model are presented. Improvements to the well-known Puck Failure criterion and a recently developed Generalized Tabulated Failure criterion are used to illustrate how to link a failure sub-model to existing deformation and damage sub-models in the context of explicit finite element analysis. These models are implemented in LS-DYNA, a commercial transient dynamic finite element code. Two validation tests are used to evaluate the failure sub-model implementation and improvements - a stacked-ply test carried out at room temperature under quasi-static tensile and compressive loadings, and a high-speed, projectile impact test where there is significant damage and material failure of the impacted panel. Results indicate that developed procedures and improvements provide the analyst with a reasonable and systematic approach to building predictive impact simulation models.

Keywords

Orthotropic composite plasticity failure modeling impact loads

Introduction

Amongst several key topics in building an end-to-end predictive modeling framework for composites, three focus areas are: predictions of elastic and plastic deformation, damage initiation and progression, and the failure of structural systems subjected to quasi-static and impact loads. In addition to the complexity of accurately simulating the deformation and damage response of composites, accurately predicting composite failure, particularly in the context of a structural analysis, is a huge challenge. For example, composites fail in compression in a variety of ways – fiber crushing, splitting, elastic microbuckling, matrix failure and plastic microbuckling, buckle delamination, shear-band formation, etc. There are other modes of failure involving fracture of the material that is brought about by fiber breakage, fiber micro-buckling, fiber pullout, matrix cracking, delamination, debonding, or any combination of these mechanisms. In fact, after three Worldwide Failure Exercises that started in the early 1990s, the consensus is that there exist several damage and failure theories none of whom are capable of “…predicting all the strength and deformation of the lamina and laminates within 10% of the measured data”.¹

The architecture of even the simplest composite leads to a complex state of stress in structural systems. While unidirectionally loaded test specimens yield important stress-strain data in the Principal Material Directions (PMDs), the damage and failure states obtained from the test results are for the most part, for simple states of stress. The interaction between different states of stress are captured by higher level theories, such as Tsai-Wu, for predicting when the material would yield and how the yield surface would evolve.² However, the damage and failure theories lack the resolution to accurately predict damage and especially failure, at the continuum scale where the composite architectural details are missing. While there is widespread consensus that these architectural details need to be a part of the structural analysis process, there is no consensus as to how to define them and then how to use them at the appropriate scale primarily because of the enormous computational expense in carrying out multi-scale and multi-level finite element (FE) computations.

An orthotropic elasto-plastic damage material model suitable for impact analysis of composite materials has been developed through a joint research project funded by the Federal Aviation Administration and the National Aeronautics and Space Administration.^3–6 The material model is implemented into LS-DYNA⁷ as MAT_ COMPOSITE_TABULATED_PLASTICITY_ DAMAGE, also known as MAT_213. MAT_213 is a constitutive model with a modular architecture comprised of three sub-models - deformation, damage and failure. The deformation sub-model predicts the linear and the non-linear response of a composite and is driven by a set of tabulated stress-strain data set. The damage sub-model predicts the reduction in the stiffness of the material. The failure sub-model predicts when there is no more load carrying capacity in the FE, carries out the computations for a gradual decrease in the stiffness of the element, and tags the element for erosion from the FE model.

In an earlier work,³ the well-known Puck Failure Criterion (PFC) was discussed, verified and validated. PFC was implemented as a part of the failure sub-model in MAT_213. While obtaining all the required material parameters proved to be difficult, the results from the validation exercises were encouraging. Some numerical calibration was necessary to obtain realistic results. In this paper, some of the shortcomings of PFC are discussed and modifications to the theory and implementation are introduced especially with respect to post-peak/post-failure simulation so that a stable post-peak behavior can be achieved. In addition, a (new) second failure model, the Generalized Tabulated Failure Criterion (GTFC)⁸ as implemented in MAT_213, is discussed. This model is built with user input that define the two modes of failure – in-plane and out-of-plane modes. In the context of this paper, the term failure onset implies the instant at which there is no further increase in the load carrying capacity exhibited by the finite element, and the term erosion refers to the deletion of the finite element from further model calculations. The first part of the paper discusses the relevant parts of the deformation and damage sub-models. It is shown that the damage sub-model can be used not only for handling pre-peak damage state, but also can be used to affect post-peak behavior. The second part of the paper shows the improvements made to both PFC and GTFC. The discussions center around the required input, onset of failure detection, and post-peak behavior including (finite) element erosion. Numerical examples form the next section of the paper. A quasi-static stacked ply test and a high-velocity impact test of a flat panel are used as validation tests for the developed and improved constitutive models. The paper concludes with discussions on how some of the modeling challenges were overcome and what can be done to further improve these models.

Deformation and damage Sub-Models

A minimum of twelve stress-strain input curves are required to drive the deformation sub-model. The stress-strain input curves include (a) tension and compression curves in the 1, 2 and 3 PMDs, (b) shear curves in the 1–2, 2–3 and 1–3 Principal Material Planes, and (c) three 45° off-axis tension/compression curves in the Principal Material Planes. When available, additional stress-strain curves at different strain-rate and temperature combinations can also be specified. A general three-dimensional stiffness matrix is used where the components are obtained from the user input stress-strain curves for the current strain-rate and temperature in the finite element.^4–6,9 The onset of plasticity is predicted using a yield function which is a modified form of the Tsai-Wu failure criterion. The yield function coefficients are computed based on the current flow stress. The direction of the plastic flow is dictated by a quadratic plastic potential function.

The damage sub-model predicts the reduction in the stiffness of the material assuming strain equivalence between the true and the effective stress spaces.^10,11 The true stress space corresponds to the damaged state whereas the effective space corresponds to the undamaged state. Strain equivalence assumption allows the damage sub-model to predict the reduction in the stiffness without affecting the deformation sub-model. The true stress $(σ_{i j}^{})$ and the effective stress $(σ_{i j}^{eff})$ tensor components are expressed using equations (1) and (4). The normal (tensile and compressive) components of the stress are related as

σ_{i i}^{eff} = {\begin{matrix} \frac{σ_{i i}^{}}{1 - c_{i}^{d}} & i f σ_{i i}^{} \geq 0 \\ \frac{σ_{i i}^{}}{1 - c_{i + 3}^{d}} & i f σ_{i i}^{} < 0 \end{matrix}

(1)

with

\begin{array}{l} c_{i}^{d} = 1 - {(1 - d_{11_{T}}^{i i_{T}}) (1 - d_{11_{C}}^{i i_{T}}) (1 - d_{22_{T}}^{i i_{T}}) (1 - d_{22_{C}}^{i i_{T}}) \\ (1 - d_{33_{T}}^{i i_{T}}) (1 - d_{33_{C}}^{i i_{T}}) (1 - d_{12}^{i i_{T}}) (1 - d_{23}^{i i_{T}}) (1 - d_{13}^{i i_{T}})} \end{array}

(2)

\begin{array}{l} c_{i + 3}^{d} = 1 - {(1 - d_{11_{T}}^{i i_{C}}) (1 - d_{11_{C}}^{i i_{C}}) (1 - d_{22_{T}}^{i i_{C}}) (1 - d_{22_{C}}^{i i_{C}}) \\ (1 - d_{33_{T}}^{i i_{C}}) (1 - d_{33_{C}}^{i i_{C}}) (1 - d_{12}^{i i_{C}}) (1 - d_{23}^{i i_{C}}) (1 - d_{13}^{i i_{C}})} \end{array}

(3)

where

i

= 1, 2 and 3. The shear components are given as

σ_{i j}^{eff} = \frac{σ_{i j}^{}}{1 - c_{k}^{d}}

(4)

with

\begin{array}{l} c_{k}^{d} = 1 - {(1 - d_{11_{T}}^{i j}) (1 - d_{11_{C}}^{i j}) (1 - d_{22_{T}}^{i j}) (1 - d_{22_{C}}^{i j}) \\ (1 - d_{33_{T}}^{i j}) (1 - d_{33_{C}}^{i j}) (1 - d_{12}^{i j}) (1 - d_{23}^{i j}) (1 - d_{13}^{i j})} \end{array}

(5)

where,

i j

= 12, 23 and 13, and

k

= 7, 8 and 9, respectively. The parameter,

c_{i}^{d}

, is referred to as the effective damage parameter and is a measure of damage in each component in the principal material direction/plane. As will be shown later, this parameter can also be used in the post-peak regime in the failure sub-model since it relates the true and the effective stress tensor components.

d_{i j}^{k l}

is referred to as the damage parameter and accounts for damage in the

k l

direction due to loading in the

i j

direction. If the

i j

direction and the

k l

direction are the same, then

d_{i j}^{k l}

is referred to as the uncoupled damage parameter; else, it is referred to as coupled damage parameter. For example,

d_{22_{T}}^{22_{T}}

is an uncoupled damage parameter obtained by loading the specimen in 2-direction tension and interrogating the reduction in stiffness in the 2-direction tension. Similarly,

d_{22_{C}}^{22_{T}}

is a coupled damage parameter obtained by loading the specimen in 2-direction compression and interrogating the reduction in the stiffness in the 2-direction tension.

The plasticity algorithm works in the effective stress-space in order to uncouple the deformation and the damage sub-models.¹¹ Therefore, all the stress-strain input curves which are in the true stress space are converted into the effective stress space using only the associated uncoupled damage parameters. The plastic multiplier is computed during the plasticity algorithm after which the stress is updated using a radial return scheme in the effective stress space. The effective stress tensor is ultimately converted back to the true stress space. This approach will be of importance when discussing the failure sub-models.

Failure modeling

In this section, two improved versions of failure theories are discussed - Puck Failure Criterion (PFC) and Generalized Tabulated Failure Criterion (GTFC). The reason for incorporating multiple failure theories in the same material model is that, as will be discussed later, each one of these criteria has its strengths and weaknesses. GTFC provides a relatively new concept for handling failure and is driven primarily by a user-defined failure envelope.

Puck failure criterion

Puck failure theory is a physically based failure model for unidirectional composites.^12–14 It assumes that there is transverse isotropy in the material at the lamina level. The failure theory considers two modes of failure - Fiber Fracture (FF) and Inter Fiber Fracture (IFF). The implementation of PFC in MAT_213 has been discussed in an earlier publication.³ The failure modes are used to predict the failure onset in the finite element. Furthermore, a stress degradation model (SDM) which gradually degrades the material behavior after failure onset, and ultimately predicts when element erosion needs to take place, is used. SDM helps avoid a sudden drop in the stress after failure onset thus providing not only numerical stability but also serves as a mesh regularization scheme.

The SDM discussed earlier³ has been modified. In the previous implementation, at the element Gauss point under consideration, once failure onset is detected, the deformation and the damage calculations were skipped. The value of true stress corresponding to the failure onset was stored and was scaled down to compute the degraded stress. Sometimes this approach resulted in excessively elongated elements, smaller critical time step, and thereby longer simulation run time. In the improved implementation presented here, the computations related to deformation and damage are carried out even after failure onset. The effective stress is held constant, but the true stress is gradually reduced. This approach provides not only post-peak numerical stability but the ability to control the residual strength in the finite element. In the context of this paper, residual strength, $c^{R S}$ , is expressed as a fraction of the peak stress that the material has when the element is eroded, i.e. $σ^{f} = c^{R S} σ^{\max}$ , $0 \leq c^{R S} \leq 1.0$ . This gradual degradation of stress is explained in Figure 1(a) with the corresponding effective damage parameter $(c_{}^{d})$ shown in Figure 1(b). Let $σ_{}^{0}$ and $ε_{}^{0}$ represent the stress and strain values at failure onset, respectively. If the FF criterion is satisfied before the IFF criterion, only the stress in the 1-direction is degraded. It is assumed that the material has load carrying capacity in the 2-direction, 3-direction, 1–2 plane, 2–3 plane and 1–3 plane even after FF. The stress components other than the 1-direction stress are degraded only after the IFF criterion is satisfied. Conversely, if IFF is satisfied before FF, it is assumed that the composite has load carrying capacity in the 1-direction. Therefore, the stress is not degraded in the 1-direction till the FF criterion is satisfied. It should be noted that $σ_{}^{0}$ is in the true stress space, whereas the corresponding effective stress counterpart is denoted as $σ_{}^{eff, 0}$ . At the onset of failure, the initial value of the effective damage parameter, $c_{0}^{d}$ , is computed from the damage parameter $d_{k l}^{i j}$ . As can be seen in Figure 1(a), the effective stress is kept constant at $σ_{}^{eff, 0}$ since a non-negative slope of the effective stress is required by the plasticity algorithm. The true stress is brought down to the residual strength, $σ_{}^{f}$ , using the deformation and the damage sub-models with the corresponding strain value of $ε_{}^{1}$ . $σ_{}^{f}$ is computed based on the value of Post Peak Residual Damage ( $PPRD$ ) and the corresponding effective damage parameter stays constant at $PPRD$ as shown in Figure 1(b). $PPRD$ is a user input parameter for PFC. The residual strength is maintained constant till the element is eroded at the erosion strain value of $ε_{}^{f}$ that is computed as

ε_{}^{f} = ε_{}^{0} + \frac{2 Γ}{σ_{}^{0} L}

(6)

where,

Γ

is the fracture energy in the mode under consideration, and

L

is the characteristic length of the finite element (taken as the cube root of the volume of the element).

Figure 1.

(a) Stress vs strain response using PFC (note $σ^{\max} = σ^{0}$ ) (b) Effective damage parameter vs strain using PFC.

Equation (6) is based on the assumption that the area under the stress-strain curve between $ε_{}^{0}$ and $ε_{}^{f}$ , ignoring the residual strength, corresponds to fracture energy. The use of the characteristic length provides the mesh regularization feature in the failure sub-model. The effective damage parameter after failure onset can be expressed as

c_{}^{d} = {\begin{matrix} (1 - c_{0}^{d}) (\frac{ε_{} - ε_{}^{f}}{ε_{}^{f} - ε_{}^{0}}) + 1 & ε_{}^{0} \leq ε_{} \leq ε_{}^{1} \\ PPRD & ε > ε_{}^{1} \end{matrix}

(7)

For $ε < ε_{}^{0}$ , $c_{}^{d}$ is computed using equations (2), (3) and (5).

Twenty parameters are used in the PFC to detect failure onset, implement the SDM algorithm, and carry out element erosion. Most, not all, can be readily found from experiments. The tensile, compressive and shear strength parameters, ${R_{∥}^{t}, R_{⊥}^{A t}, R_{| |}^{c}, R_{⊥}^{c}, R_{⊥ ∥}^{A}}$ can be obtained from experiments.¹⁵ The procedure to determine the inclination parameters, ${p_{⊥ ∥}^{c}, p_{⊥ ∥}^{t}, p_{⊥ ⊥}^{t}, p_{⊥ ⊥}^{c}}$ , is discussed in Puck et al.¹⁶ The magnification factor, $m_{σ_{f}}$ , is used in the FF criterion and accounts for stress magnification in the fiber due to inhomogeneity in the matrix stress distribution in the vicinity of the fiber.¹³ In addition, in the FF criterion, the fiber Poisson’s ratio, $ν_{∥ ⊥ f}$ , and modulus of elasticity, $E_{∥ f}$ are used.¹⁷ Next, the fiber direction fracture energy, $Γ_{f}$ , can be determined through experiments.¹⁸ This value has been estimated in this paper since no experimental data is available for the composite. Mode I and Mode II interlaminar fracture energies, ${Γ_{I}^{}, Γ_{I I}^{}}$ , are determined through experiments.¹⁹ In the absence of experimental data, the tension, compression and shear-related post-peak residual damage parameters, ${PPR D_{T ∥}^{}, PPR D_{C ∥}^{}, PPR D_{T ⊥}^{}, PPR D_{C ⊥}^{}, PPR D_{S}^{}}$ , are estimated values.

The algorithm used in PFC with the new SDM implementation is as follows.

Input: The strain tensor for the current time step, stress tensor from the previous time step (true stress space), and failure-related parameters for a specific finite element. The current failure state is also required as input and is one of the following – (a) the Gauss point has not satisfied any of the failure criterion, or (b) the Gauss point has satisfied a failure criterion and is undergoing post-failure stress degradation.

Output: Updated effective damage parameters and erosion status of the element.

Step 1.1(a): If the current point is not in FF failure state, update the stresses using the effective damage parameters $c_{1}^{d}$ and $c_{4}^{d}$ (equations (2) and (3)).

Step 1.1(b): If the current point is in FF state, update the stresses using the following effective damage parameters.

If $σ_{11}^{} \geq 0$

\begin{array}{l} c_{1, t}^{d} = c_{1, t - 1}^{d} + Δ c_{1, t}^{d} where \\ Δ c_{1, t}^{d} = {\begin{matrix} (1 - c_{1, t - 1}^{d}) (\frac{| ε_{11, t}^{} | - | ε_{11, t - 1}^{} |}{| ε_{11}^{f} | - | ε_{11, t - 1}^{} |}) & ε_{11}^{0} \leq ε_{11} \leq ε_{11}^{1} \\ 0 & c_{1, t}^{d} > PPR D_{T ∥}^{} \end{matrix} \end{array}

(8)

else

\begin{array}{l} c_{4, t}^{d} = c_{4, t - 1}^{d} + Δ c_{4, t}^{d} where \\ Δ c_{4, t}^{d} = {\begin{matrix} (1 - c_{4, t - 1}^{d}) (\frac{| ε_{11, t}^{} | - | ε_{11, t - 1}^{} |}{| ε_{11}^{f} | - | ε_{11, t - 1}^{} |}) & ε_{11}^{0} \leq ε_{11} \leq ε_{11}^{1} \\ 0 & c_{4, t}^{d} > PPR D_{C ∥}^{} \end{matrix} \end{array}

(9)

where

ε_{i j}^{}

are the strain tensor components. The subscript

t

is used to denote the current time step.

Step 1.2(a): If the current point is not in IFF failure state, update the stresses using the effective damage parameters, $c_{2}^{d}$ , $c_{5}^{d}$ , $c_{3}^{d}$ , $c_{6}^{d}$ , $c_{7}^{d}$ , $c_{8}^{d}$ and $c_{9}^{d}$ (equations (2), (3) and (5)).

Step 1.2(b): If the current point is in IFF state, update the stresses using the effective damage parameters as follows.

For $i$ -direction normal components (i = 2, 3):

If $σ_{i i}^{} \geq 0$

c_{i, t}^{d} = c_{i, t - 1}^{d} + Δ c_{i, t}^{d} where Δ c_{i, t}^{d} = {\begin{matrix} (1 - c_{i, t - 1}^{d}) (\frac{| ε_{i i, t}^{} | - | ε_{i i, t - 1}^{} |}{| ε_{i i}^{f} | - | ε_{i i, t - 1}^{} |}) & ε_{i i}^{0} \leq ε_{i i} \leq ε_{i i}^{1} \\ 0 & c_{i, t}^{d} > PPR D_{T ⊥}^{} \end{matrix}

(10)

else

\begin{array}{l} c_{i + 3, t}^{d} = c_{i + 3, t - 1}^{d} + Δ c_{i + 3, t}^{d} where \\ Δ c_{i + 3, t}^{d} = {\begin{matrix} (1 - c_{i + 3, t - 1}^{d}) (\frac{| ε_{i i, t}^{} | - | ε_{i i, t - 1}^{} |}{| ε_{i i}^{f} | - | ε_{i i, t - 1}^{} |}) & ε_{i i}^{0} \leq ε_{i i} \leq ε_{i i}^{1} \\ 0 & c_{i + 3, t}^{d} > PPR D_{C ⊥}^{} \end{matrix} \end{array}

(11)

For $i j$ shear component (ij = 12, 23, 13, k = 7, 8, 9):

c_{k, t}^{d} = c_{k, t - 1}^{d} + Δ c_{k, t}^{d} where Δ c_{k, t}^{d} = {\begin{matrix} (1 - c_{k, t - 1}^{d}) (\frac{| ε_{i j, t}^{} | - | ε_{i j, t - 1}^{} |}{| ε_{i j}^{f} | - | ε_{i j, t - 1}^{} |}) & ε_{i j}^{0} \leq ε_{i j} \leq ε_{i j}^{1} \\ 0 & c_{k, t}^{d} > PPR D_{S}^{} \end{matrix}

(12)

Step 2.1(a): If the current point is not in FF state, check the FF criterion. FF onset is assumed to take place if the stress exposure, $f_{EFF}^{}$ , is equal to or greater than one¹³

f_{EFF}^{} = \frac{| σ_{∥}^{e q} |}{R}

(13)

σ_{∥}^{e q} = σ_{11}^{} - (ν_{∥ ⊥} - ν_{∥ ⊥ f} \frac{E_{∥}}{E_{∥ f}} m_{σ_{f}}) (σ_{22} + σ_{33})

(14)

where

{\begin{matrix} R = R_{| |}^{t}, i f σ_{∥}^{e q} \geq 0 \\ R = R_{| |}^{c}, i f σ_{∥}^{e q} < 0 \end{matrix}

(15)

When $f_{EFF}^{} \geq 1$ , compute the failure strain as

ε_{11}^{f} = ε_{11}^{0} + \frac{2 Γ_{f}^{}}{σ_{11}^{0} L}

(16)

Step 2.1(b): If the current point is in FF state, erode the element if one of the following criteria are satisfied:

(i) $| ε_{11}^{} | \geq | ε_{11}^{f} |$

(ii) $(σ_{11}^{0}) (σ_{11}^{}) < 0$

Step 2.2(a): If the current point is not in IFF state, check the IFF criterion. IFF onset is assumed to take place if the stress exposure, $f_{EIFF}^{}$ , computed for $θ = θ_{f p},$ is equal to or greater than one using equation (17).¹³ $θ_{f p}$ is the orientation of the fracture plane with respect to the 3-direction. $f_{EIFF}^{}$ is computed by varying $θ$ from $- 90 °$ to $90 °$ . The angle $θ$ which yields the maximum value of $f_{EIFF}^{}$ is denoted $θ_{f p}$ .

f_{EIFF}^{} = {\begin{matrix} \sqrt{{[(\frac{1}{R_{⊥}^{A t}} - \frac{p_{⊥ ψ}^{t}}{R_{⊥ ψ}^{A}}) σ_{n}^{} (θ)]}_{}^{2} + {(\frac{τ_{n t}^{} (θ)}{R_{⊥ ⊥}^{A}})}_{}^{2} + {(\frac{τ_{n 1}^{} (θ)}{R_{∥ ⊥}^{A}})}_{}^{2}} + \frac{p_{⊥ ψ}^{t}}{R_{⊥ ψ}^{A}} σ_{n}^{} (θ) & σ_{n}^{} \geq 0 \\ \sqrt{{(\frac{τ_{n t}^{} (θ)}{R_{⊥ ⊥}^{A}})}_{}^{2} + {(\frac{τ_{n 1}^{} (θ)}{R_{∥ ⊥}^{A}})}_{}^{2} + {(\frac{p_{⊥ ψ}^{c}}{R_{⊥ ψ}^{A}} σ_{n}^{} (θ))}_{}^{2}} + \frac{p_{⊥ ψ}^{c}}{R_{⊥ ψ}^{A}} σ_{n}^{} (θ) & σ_{n}^{} < 0 \end{matrix}

(17)

where

σ_{n}^{} = σ_{22}^{} \cos_{}^{2} θ + σ_{33}^{} \sin_{}^{2} θ + 2 σ_{23}^{} \sin θ \cos θ

(18)

τ_{n t}^{} = - σ_{22}^{} \sin θ \cos θ + σ_{33}^{} \sin θ \cos θ + σ_{23}^{} (\cos_{}^{2} θ - \sin_{}^{2} θ)

(19)

τ_{n 1}^{} = σ_{13}^{} \sin θ + σ_{12}^{} \cos θ

(20)

R_{⊥ ⊥}^{A} = \frac{R_{⊥}^{c}}{2 (1 + p_{⊥ ⊥}^{c})}

(21)

(\frac{p_{⊥ ψ}^{t, c}}{R_{⊥ ψ}^{A}}) = (\frac{p_{⊥ ⊥}^{t, c}}{R_{⊥ ⊥}^{A}}) \cos^{2} (ψ) + (\frac{p_{⊥ ∥}^{t, c}}{R_{∥ ⊥}^{A}}) \sin^{2} (ψ)

(22)

\cos^{2} (ψ) = \frac{τ_{n t}^{2}}{τ_{n t}^{2} + τ_{n l}^{2}}

(23)

When $f_{EIFF}^{} = 1$ , compute the following failure strains:

ε_{i j}^{f} = ε_{i j}^{0} + \frac{2 Γ_{I}^{}}{σ_{i j}^{0} L}, for i j = 22, 33

(24)

ε_{i j}^{f} = ε_{i j}^{0} + \frac{2 Γ_{I I}^{}}{σ_{i j}^{0} L}, for i j = 12, 23, 13

(25)

Step 2.2(b): If the current point is in the IFF state, erode the element if at least one of the following criteria are satisfied:

(i) $| ε_{i j, t}^{} | \geq | ε_{i j}^{f} |, i j = 22, 33, 12, 23, 13$

(ii) ${\begin{cases} i f σ_{22}^{0} \geq σ_{33}^{0} \\ {\begin{cases} i f (σ_{22}^{0} \geq 0 and σ_{22}^{0} \geq F_{S}^{} R_{⊥}^{A t}) \\ o r i f (σ_{22}^{0} < 0 and σ_{22}^{0} \geq F_{S}^{} R_{⊥}^{c}) \\ e rode element i f σ_{22}^{0} \times σ_{22}^{} < 0 \end{cases} \\ e lse \\ {\begin{cases} i f (σ_{33}^{0} \geq 0 and σ_{33}^{0} \geq F_{S}^{} R_{⊥}^{A t}) \\ o r i f (σ_{33}^{0} < 0 and σ_{33}^{0} \geq F_{S}^{} R_{⊥}^{c}) \\ e rode element i f σ_{33}^{0} \times σ_{33}^{} < 0 \end{cases} \end{cases}$

(iii) ${\begin{cases} i f (σ_{12}^{0} \geq σ_{23}^{0}) and (σ_{12}^{0} \geq σ_{13}^{0}) and (σ_{12}^{0} \geq F_{S}^{} R_{∥ ⊥}^{A}) \\ {erode element i f σ_{12}^{0} \times σ_{12}^{} < 0 \\ e lse i f (σ_{23}^{0} \geq σ_{12}^{0}) and (σ_{23}^{0} \geq σ_{13}^{0}) and (σ_{23}^{0} \geq F_{S}^{} R_{⊥ ⊥}^{A}) \\ {erode element i f σ_{23}^{0} \times σ_{23}^{} < 0 \\ e lse i f (σ_{13}^{0} \geq σ_{12}^{0}) and (σ_{13}^{0} \geq σ_{23}^{0}) and (σ_{13}^{0} \geq F_{S}^{} R_{∥ ⊥}^{A}) \\ {erode element i f σ_{13}^{0} \times σ_{13}^{} < 0 \end{cases}$

Sometimes after an element is eroded, the neighboring elements undergoing the stress degradation process may experience sudden changes in the strain field due to stress redistribution. This process may lead to undesirable stress oscillations (stress reversals, tension to compression or vice versa) taking place in a very short duration of time. Criteria (ii) and (iii) in Step 2.2(b) erode the elements in which stress reversal with respect to the stress at the failure onset takes place. Since the PFC assumes transverse isotropy, care needs to be taken regarding which stress component ( $σ_{22}^{}$ or $σ_{33}^{}$ ) has to be chosen to keep track of stress reversal. This process is taken care of by the criterion in (ii). Similar strategy is used for the shear components. $F_{S}^{}$ is used here to ensure that the elements in which the stress components are oscillating are eroded. Numerical experimentation has shown that a value of $F_{S}^{} = 0.8$ has worked reasonably well and is currently used for all the verification and validation tests presented in this paper.

Although Puck failure theory has been considered one of the better models in World-Wide Failure Exercise I (WWFE-I) and World-Wide Failure Exercise II (WWFE-II),^20,21 there are limitations that should be noted. First, the theory is valid for transversely isotropic material at the lamina level. Second, the Master Fracture Body is not a closed envelope.¹³ Third, the failure criteria are not functions of rate and temperature. This factor is important for unidirectional polymeric composites that exhibit rate and temperature sensitivities and form the impacted structural system.

Generalized tabulated failure criterion

A failure model where an arbitrarily shaped failure surface in the stress space can be used to predict the failure of a composite is highly desirable. One would expect the failure surface of a composite to be complex such that a mathematical expression cannot be defined and used as is done with traditional failure criteria.^22,23 In the earlier version of GTFC,⁸ the in-plane, the out-of-plane and the combined in-and-out-of-plane failure states were defined in terms of stresses. The current version of GTFC assumes that failure states are strain rather than stress-based so that post-peak behavior can be handled seamlessly. Failure surfaces need to be defined for each of these failure states in the equivalent failure strain and failure angle space.²⁴

Figure 2 shows the in-plane failure surface where $ε_{I P_{}^{FAIL}}^{e q}$ is a function of the failure angle, $θ_{I P}^{}$ , and can take any shape depending on the composite properties. The equivalent strain, $ε_{I P}^{e q}$ and failure angle, $θ_{I P}^{}$ are computed for each time step and element Gauss point as

ε_{I P}^{e q} = \sqrt{ε_{11}^{2} + ε_{22}^{2} + 2 ε_{12}^{2}}

(26)

θ_{I P}^{} = \cos_{}^{- 1} (\frac{σ_{22}}{\sqrt{σ_{22}^{2} + σ_{12}^{2}}})

(27)

Figure 2.

General form of in-plane failure surface.

Note that the failure angle is expressed in terms of $σ_{22}^{}$ and $σ_{12}^{}$ as in unidirectional composites, the stress in the 1-direction, $σ_{11}^{}$ , is usually of a magnitude higher than $σ_{22}^{}$ and $σ_{12}^{}$ . Using $σ_{11}^{}$ in the computation of the failure angle⁸ instead of using $σ_{12}^{}$ may make the points on the failure surface agglomerate as illustrated in Figure 3 where most of the data points are concentrated at $θ_{I P}^{}$ values of $-^{180} \circ$ , $^{0} \circ$ and $^{180} \circ$ . This result is not desirable since the interpolation between unequally spaced points may yield inaccurate equivalent failure strain values.

Figure 3.

In-plane failure surface where failure angle is computed using $σ_{11}^{}$ .

Let the in-plane failure state $d_{1}^{}$ be computed as

d_{1}^{} = \frac{ε_{I P}^{e q}}{ε_{I P_{}^{FAIL}}^{e q}}

(28)

An out-of-plane failure surface can be constructed similar to Figure 2 with the failure surface expressed in terms of $ε_{O O P_{}^{FAIL}}^{e q}$ as a function of $θ_{OOP}^{}$ for a given $σ_{33}^{}$ . $ε_{OOP}^{e q}$ and $θ_{OOP}^{}$ are computed for each time step and element Gauss point as

ε_{OOP}^{e q} = \sqrt{ε_{33}^{2} + 2 ε_{13}^{2} + 2 ε_{23}^{2}}

(29)

θ_{OOP}^{} = \cos_{}^{- 1} (\frac{σ_{13}}{\sqrt{σ_{13}^{2} + σ_{23}^{2}}})

(30)

Similarly, the out-of-plane failure state $d_{2}^{}$ is computed as

d_{2}^{} = \frac{ε_{OOP}^{e q}}{ε_{O O P_{}^{FAIL}}^{e q}}

(31)

An element is eroded if $d \geq 1$ where

d = \{\begin{cases} \max (_{d}^{1},_{d}^{2}) & if n = 0 \\ \sqrt[n]{_{(_{d}^{1})}^{n} +_{(_{d}^{2})}^{n}} & if n > 0 \end{cases}

(32)

where

n

is a user defined interaction parameter that can be used to couple the in-plane state of stress/strain to the out-of-plane state of stress/strain. GTFC requires two sets of tabulated input - the equivalent failure strain and the corresponding failure angle for in-plane

(ε_{I P_{}^{FAIL}}^{e q}, θ_{I P}^{})

and the out-of-plane

(ε_{O O P_{}^{FAIL}}^{e q}, θ_{OOP}^{})

failure modes. In this paper, these values are based on data obtained from uniaxial testing,¹⁵ though it is desirable to generate a richer set of data either via biaxial/triaxial laboratory tests or virtual tests. A minimum of 21 parameters (9 stress-strain curves that are augmented with post-peak data, 9 uncoupled damage parameter curves, set of in-plane

(ε_{I P_{}^{FAIL}}^{e q}, θ_{I P}^{})

values and set of out-of-plane

(ε_{O O P_{}^{FAIL}}^{e q}, θ_{OOP}^{})

values, and

n

) are used in GTFC to compute the post-peak response including element erosion.

The algorithm used in GTFC is as follows.

Input: The strain tensor for the current time step, stress tensor from the previous time step which are in the true stress space, and failure-related parameters.

Output: Element erosion status.

Step 1: For a given time step, update the element Gauss point stress tensor using the deformation and the damage sub-models.

Step 2.1: Compute $ε_{I P}^{e q}$ and $θ_{I P}^{}$ using equations (26) and (27). Using the interpolation function, obtain $ε_{I P_{}^{FAIL}}^{e q}$ corresponding to the computed $θ_{I P}^{}$ .

Step 2.2: Compute $d_{1}^{}$ using equation (28).

Step 3.1: Compute $ε_{OOP}^{e q}$ and $θ_{OOP}^{}$ using equations (29) and (30). Using the interpolation function, obtain $ε_{O O P_{}^{FAIL}}^{e q}$ corresponding to the computed $θ_{OOP}^{}$ .

Step 3.2: Compute $d_{2}^{}$ using equation (31).

Step 4: Compute $d$ using equation (32). Erode the element if $d \geq 1$ .

Unlike PFC, GTFC does not require failure onset to be detected to compute the post-peak response. Each input stress-strain curve has a pre-peak and post-peak regime that are used to capture the peak stress via the deformation and damage sub-models, and handle the post-peak behavior via the damage sub-model similar to what is done with the stress degradation model in PFC. Figure 4(a) shows a typical stress-strain input curve for a given component. The portion of the stress data between $ε_{}^{0}$ and $ε_{}^{1}$ can be artificially generated if experimental data is not available. Note that $c_{i}^{R S}$ corresponding to $i = 1, . ., 6$ represent the ratio of the peak stress to the residual strength for the normal tensile (1, 2, 3)/compressive (4, 5, 6) components in the PMDs, and $i = 7, 8, 9$ for the three principal shear components. For the deformation sub-model to carry out the plasticity-related computations, the input stress in the effective stress space cannot have a negative slope. To ensure that this condition is satisfied, a damage parameter curve is used as input (Figure 4(b)). This damage parameter data is computed for the respective component by keeping the effective stress constant and equal to the peak effective stress. The uncoupled damage parameter is the same as the effective damage parameter as no coupled-damage parameters are used in the post-peak calculations.

Figure 4.

(a) Stress-strain input augmented with post-peak data (note $σ^{\max} = σ^{0}$ ) (b) Effective damage parameter input for handling post-peak behavior.

Like PFC, the current implementation of GTFC is not rate and temperature dependent. Obtaining experimentally a rich set of failure surface data is neither easy nor feasible with the current state-of-the-art. The combination of experimental testing along the principal material directions and virtual testing for combined states of stress appears to be a viable option. To alleviate this problem, in this paper, the in-plane and out-of-plane responses are computed separately and then combined via the interaction term, $n$ .

It should be noted that computation of the failure strain (strain at which element is eroded) is different for PFC and GTFC. In the case of the PFC, the failure strains are computed in the material model as soon as the element reaches the failure onset (Figure 1), and are computed separately for the FF mode (equation (16)) and the IFF mode (equations (24) and (25)). In the case of GTFC, the failure strains are specified as user input in the form of equivalent failure strains in the two modes: in-plane, $ε_{I P_{}^{FAIL}}^{e q}$ , and out-of-plane, $ε_{O O P_{}^{FAIL}}^{e q}$ .

Numerical results

General Details: The verification and validation for both GTFC and PFC implementations were done using the T800/F3900²⁵ unidirectional composite. T800/F3900 is a carbon fiber composite with epoxy matrix and a fiber volume fraction of roughly 60%. Experiments show that the composite exhibits brittle failure with little or no post-peak strength.^15,26 One-element verification tests are used to explain the implementation nuances of both failure criteria. Two validation tests are used to determine the efficiency and accuracy of the failure theories – a stacked-ply test carried out at Arizona State University (ASU) with loading in the quasi-static regime, and a high-velocity impact test carried out at National Aeronautics and Space Administration – Glenn Research Center (NASA-GRC). Three different mesh sizes have been used for the stacked-ply validation tests – coarse, medium, and fine with increasing mesh densities to study mesh dependencies and convergence properties. For the sake of brevity, details of the convergence analyses used in generating the optimal mesh (one that balances accuracy and computational cost) for the impact validation test are not presented here but can be found in Hoffarth et al.⁹ In the impact validation tests, cohesive zone elements (CZE) have been used for modeling the delamination behavior and the details can be found in an earlier publication.¹⁹

PFC Parameters: The same input parameters used in our previous research³ are used here: $R_{∥}^{t} = 366100 psi (2524 MPa),$ $R_{⊥}^{A t} = 6500 psi (44.8 MPa),$ $R_{| |}^{c} = 106000 psi (731 MPa),$ $R_{⊥}^{c} = 25500 psi (176 MPa),$ $R_{⊥ ∥}^{A} = 18600 psi (128 MPa),$ $p_{⊥ ∥}^{c} = 0.30,$ $p_{⊥ ∥}^{t} = 0.35$ , $p_{⊥ ⊥}^{t} = 0.25,$ $p_{⊥ ⊥}^{c} = 0.30,$ $ν_{⊥ ∥} = 0.0168,$ $ν_{⊥ ∥ f} = 0.2,$ $m_{σ_{f}} = 1.1,$ $Γ_{f} = 400 l b / i n (7 \times 10^{4} N / m)$ , $Γ_{I}^{} = 2.15 l b / i n (376 N / m)$ , $Γ_{I I}^{} = 10.4 l b / i n (1821 N / m)$ , $E_{∥} = 2.37 \times 10^{7} p s i (158 GPa)$ in tension, $E_{∥} = 1.7 \times 10^{7} p s i (117 GPa)$ in compression, and $E_{∥ f} = 1.8 \times 10^{7} p s i (124 GPa)$ . The only new input parameters are the $PPRD$ parameters whose values are estimated. First, it is assumed that there is no residual strength in tension and shear based on the observations made from the uniaxial stress experiments conducted on T800/F3900 composite specimens that were loaded until failure.¹⁵ Second, $PPR D_{C ∥}^{}$ is estimated assuming the compressive strength in the 1-direction after the failure onset is approximately equal to the compressive strength in the 2-direction and there is a 50% residual strength in the transverse compressive direction. Composite panels still have some load carrying capacity after failure onset since the loading is in compression. In summary, unless otherwise stated, the following values are used in the FE models: $PPR D_{T ∥}^{} = 0.99$ , $PPR D_{C ∥}^{} = 0.255$ , $PPR D_{T ⊥}^{} = 0.99,$ $PPR D_{C ⊥}^{} = 0.5,$ and $PPR D_{S}^{} = 0.99$ .

GTFC Parameters: For all the simulations using GTFC, it has been assumed that the compression residual strength parameters are equal to the respective peaks in the stress-strain input curves, $c_{i}^{R S} = 1.0, i = 4, 5, 6$ . The residual strength used for the tension and the shear components are discussed for each test case separately. The equivalent failure strain is a constant function of the failure angle for both the in-plane and the out-of-plane modes. In the absence of experimental data, it is assumed that $ε_{I P_{}^{FAIL}}^{e q} = ε_{O O P_{}^{FAIL}}^{e q}$ . The interaction term, $n$ is taken as 2 thereby coupling the in-plane and the out-of-plane states of stress that are expected in an impact event.

Numerical Calibration of Failure-related Parameters: As stated earlier, there are numerous failure-related parameters that can be altered to achieve better predictions. Amongst all failure-related parameters, the following parameters were used to improve the finite element predictions via a few trial-and-error runs: PFC: $R_{| |}^{c}$ , $Γ_{f}$ , $Γ_{I}^{}$ , $Γ_{I I}^{},$ $PPR D_{T ∥}^{},$ $PPR D_{C ∥}^{}$ , $PPR D_{T ⊥}^{}$ , $PPR D_{C ⊥}^{}$ , and $PPR D_{S}^{}$ . GTFC: $c_{i}^{R S}, i = 1, . ., 9,$ $ε_{I P_{}^{FAIL}}^{e q}, ε_{O O P_{}^{FAIL}}^{e q}$ .

Nomenclature: The following nomenclature is used in the numerical example graph legends: (a) MAT213-GTFC and MAT213-PFC imply that the simulation was run with GTFC and PFC, respectively, (b) the post-fixes -F, -M and -C imply that the simulation was run with fine, medium and coarse mesh model, respectively, (c) the post-fix -Modified Strength implies that one or more failure-related values have been calibrated, (d) Experiment implies experimental data from a single test replicate (e.g., impact test), (e) for all MAT213-PFC simulations, Model represents the processed data from the experimentally obtained stress-strain curves,¹⁵ and (f) for all the MAT213-GTFC simulations, Model implies that calibrated post-peak stress-strain data is added to the experimentally obtained pre-peak data.

PFC single element verification test

Background: One single-element example is discussed to illustrate the new SDM implementation. The example considered is 2-direction tension model. The test is executed with two different sets of input parameters that influence the post-peak behavior (Table 1). $Γ_{I I}^{}$ shown in the Set 2 column is a fictitious value and is used to show the significance of the parameter.

Table 1.

Input parameters for verification.

Component	Set 1	Set 2
$PPR D_{T ⊥}^{}$	0.5	0.7
$Γ_{I I}^{}$	10.4 lb/in (1821 N/m)	40.4 lb/in (7074 N/m)

Finite Element Modeling: The model involving the 8-noded hexahedral element is shown in Figure 5. The boundary conditions are applied in order to obtain a uniaxial state of stress in the 2-direction. The pin support represents restraints in all the three translational displacements. The crossed red color arrows represent a restraint in the translational displacement along the direction of the arrow. The black color arrows represent applied displacements.

Figure 5.

Schematic of 2-direction tension test model.

Results: The Model curve as well as the stress-strain response of the FE models are shown in Figure 6 where the effective damage parameter is shown using the secondary axis.

Figure 6.

PFC input and results: 2-direction tension test stress and effective damage parameter versus strain response.

Discussions: The pre-peak stress-strain responses remain the same for Set 1 and Set 2. The difference is observed in the post-peak regime. The larger the value of $Γ_{I I}^{}$ , the larger the area is under the curve in the post-peak regime. It can also be observed that, the smaller the value of $PPR D_{T ⊥}^{}$ , the larger the residual strength. These observations help verify the implementation of PFC in the failure sub-model.

GTFC single element verification test

Background: This test is the same single element test used with the PFC, but now executed with GTFC.

Finite Element Modeling: As can be seen in Figure 7(a), there is an additional post-peak regime in the Model curve. As explained earlier, the stress-strain data in the pre-peak regime are obtained from experiments, whereas the stress-strain data in the post-peak regime are synthetic data that is used to control the post-peak stress degradation behavior. The figure also shows the effective damage parameter curve which increases from 0.0 to a value of 0.9 (see secondary axis), and the corresponding effective stress which remains constant. In this test case, the residual strength in all the tension and the shear components is 10% of the respective peaks in the stress-strain input curve. The values of $ε_{I P, FAIL}^{e q} = 0.01$ and $ε_{OOP, FAIL}^{e q} = 0.01$ are calibrated values that are obtained with the objective to ensure that the element erodes after stress degradation.

Figure 7.

GTFC input and results: 2-direction tension test (a) Stress and effective damage parameter versus strain response (b) Equivalent strain versus strain (c) d, d₁ and d₂ versus strain.

Results: The drop in the stress value at a strain of ∼0.008 is because of element erosion. Figure 7(b) shows $ε_{I P}^{e q}$ and $ε_{OOP}^{e q}$ increasing with respect to the 2-direction strain with the same trend as $d_{1}^{}$ and $d_{2}^{}$ shown in Figure 7(c), respectively.

Discussions: The simulation stress-strain response matches the Model curve in both the pre-peak and the post-peak regimes. As is evident from Figure 7(c), the element is eroded when $d = 1$ . These observations help verify the implementation of GTFC in the failure sub-model.

Stacked-Ply validation test with PFC and GTFC

Background: The first set of validation tests was done using stacked-ply specimens consisting of 8 plies with a lay-up of [0/90/45/–45]_s as shown in Figure 8. Two different tests were carried out – (a) specimen subjected to tension in the y-direction, and (b) specimen subjected to compression in the y-direction. The grey colored regions represent the portion of the specimen which were tabbed using glass fiber tabs and gripped in the test fixture.^3,27

Figure 8.

Schematic diagram (a) tension test (b) compression test.³

Finite Element Modeling: The tabbed regions in the specimen were excluded from the finite element model to create a computationally efficient model. Figure 9 shows the FE models with the fine mesh used for the stacked-ply validation tension test (SPVTT) and stacked-ply compression validation test (SPVCT). 8-node hexahedral elements were used. There are eight elements through the thickness corresponding to the eight plies in the experimental specimen. Table 2 shows FE model details for the fine, medium and the coarse meshes. The SPVTT boundary conditions are shown in Figure 10 using the highlighted nodes in the fine mesh model. All of the nodes on the support face were restrained from displacement in the y-direction as shown in Figure 10(a). The center line nodes through the z-direction were restrained from displacement in the x-direction as shown in Figure 10(b). This process was done to avoid strain localization on the support face. Since only the gage section is modeled, the displacement at the end of the gage section from the experiment was taken from the Digital Image Correlation (DIC) analysis²⁷ and applied to the FE model along the positive y-direction on the loading face nodes as shown in Figure 10(c). The nodes on the loading face were also restrained in the x-direction. To avoid rigid body motion in the z-direction (out-of-plane), the center line nodes through the x-direction on either end are restrained in the z-direction direction as shown in Figure 10(d). The boundary conditions are similar for both SPVTT and SPVCT models except for the applied loading. Also, in the case of SPVCT, all the nodes at the support face were restrained in displacements in both the x and the y-directions. LS-DYNA’s⁷ reduced integration option has been used with appropriate hourglass control.

Figure 9.

FE model used for (a) stacked-ply validation tension test (SPVTT) (b) stacked-ply validation compression test (SPVCT).

Table 2.

Number of elements in the FE models.

Test	Number of elements (largest aspect ratio)
Test	Fine	Medium	Coarse
SPVTT	38,400 (3.2)	9600 (6.4)	2400 (12.8)
SPVCT	25,600 (1.6)	6400 (3.2)	1600 (6.4)

Figure 10.

SPVTT boundary conditions (a) y-displacement restraint (b) x-displacement restraint (c) x-displacement restraint and applied displacement along y-direction (d) z-displacement restraint.

Results: The load versus time graph is used to compare FE and experimental results. For the FE model, the nodal reactions at the support face along the loading direction were added to obtain the load value. The experimental values were obtained from the test frame load cell. The SPVTT load vs time graphs with PFC and GTFC are shown in Figure 11(a) and (b), respectively. The energy vs time graph with the fine mesh is shown in Figure 12.

Figure 11.

SPVTT load versus time plot with (a) PFC and (b) GTFC.

Figure 12.

SPVTT energy versus time plot (fine mesh model) with (a) PFC and (b) GTFC.

The SPVCT results are shown in Figure 13(a) using PFC with the fine mesh before failure-related parameters were calibrated as well as after calibration. The GTFC results are shown in Figure 13(b). No failure-related parameter calibration was performed. The corresponding energy plots are shown in Figure 14(a) and (b), respectively.

Figure 13.

SPVCT load versus time plot with (a) PFC and (b) GTFC.

Figure 14.

SPVCT energy versus time (fine mesh) with (a) PFC and (b) GTFC.

Discussion SPVTT Models: The peak load predicted using PFC agrees well with the experimental results. The overall FE response is a little stiffer. The distribution of the longitudinal stress (y-stress) is uniform in the model. IFF first initiates in the ±45° plies, then in the 90° plies, after which FF initiates in the 0° plies. This process is immediately followed by a uniform erosion of the elements throughout the entire model. The differences when there is complete loss of load carrying capacity (zero load value) amongst the three mesh sizes is probably due to the different aspect ratios used in these models. The aspect ratios of the element play an important role in element erosion (equation (6)). Overall, the prediction is an improvement compared to results using the earlier implementation of PFC.³ In the simulations using GTFC, the equivalent failure strain was assumed to be 0.8. The peak load predictions are lower than the experimental values - the coarse and the medium model are about 15% less and the fine model is about 20% less. It should be noted that after the true stress reaches the peak in each component of input stress-strain curves, the true stress starts to degrade. This finding is in contrast to PFC where the true stress degrades only after the failure criterion is satisfied. The difference in the prediction amongst the different mesh sizes (about 5%) can be attributed to the fact that GTFC is not mesh independent.

Discussion SPVCT Models: In the model using the fine mesh with PFC, the peak load prediction was about 20% higher than the experimental result. Subsequently, a few failure-related parameters were calibrated. The 1-direction compressive strength, $R_{∥}^{c}$ was reduced slightly by about 2% to 102899 psi (709 MPa), a value that is within the experimental coefficient of variation of 9.7%.¹⁵ In addition, $Γ_{f}$ was reduced to 20 lb/in (3502 N/m). The value of $PPR D_{C ∥}^{}$ was increased to 0.75 since the previously assumed value (0.255) made the FE model too stiff. After calibration, the predictions agreed with the experimentally obtained value. IFF first initiates in the ±45° plies and 90° plies, followed by FF in the 0° plies. Like the SPVTT case, differences in the point of time at which the model completely loses its of load carrying capacity compared to the experimental values can be observed. This finding can be attributed to the difference in the aspect ratios of the elements in the three models. Similarly, the GTFC failure-related parameters were calibrated. The tension and shear residual strength components were taken as 10% of the respective peak stress, and the equivalent failure strains were reduced from 0.8 to $ε_{I P_{}^{FAIL}}^{e q} = ε_{O O P_{}^{FAIL}}^{e q} = 0.016 (F), 0.011 (M), 0.008 (C)$ making the peak load prediction more accurate but mesh dependent.

In all the stacked-ply simulations, the load is predominantly carried by the 0° plies. As can be seen from the energy plots, as required, the kinetic energy and the hourglass energy are negligibly small for both PFC and GTFC. Some improvements can be made to the models. First, only the gage section of the specimens is modeled. Modeling the whole specimen with the fiberglass tabs will lead to a much larger model but may improve the predictions made for the compression test since boundary effects such as stress concentrations can be avoided. This change will also require the modeling of the interface between the tabs and the composite. In the ±45° plies in the FE model, the edges of the elements are aligned along the loading direction rather than the along the fibers. Modeling the edges of the elements in ±45° plies along the fiber direction with the use of tie-break may improve the results as failure of elements is in this case, dependent on element orientation.

Impact validation test

Background: A set of ballistic impact tests was conducted at NASA-Glenn Research Center.²⁸ These tests were designed to produce validation test data for the developed MAT_213 constitutive model. In each of the tests, the target is a composite panel impacted with a 50 gm (0.122 lbm) aluminum projectile at different velocities ranging from 119 ft/s (36.3 m/s) to 530 ft/s (161.5 m/s). The composite (flat) panel is made of 16 plies of T800/F3900 composite with a lay-up of [(0/90/45/–45)₂]_S. The dimensions of each panel are 12 in. x 12 in. x 0.122 in. (30.48 cm x 30.48 cm x 0.3099 cm). DIC data gathered from high-speed cameras were used with the ARAMIS software system^28,29 to obtain both projectile information as well as the displacement of the front and back sides of the impacted panel.

In this paper, the LVG1075 test with the projectile velocity as 385 ft/s (117.3 m/s) is chosen for validating the improved failure sub-models. Figure 15(a) shows that the panel is clamped and bolted, and Figure 15(b) shows a close-up view of the rear side of the panel with visible damage. The projectile rebounded after the impact with an average velocity of 46.4 ft/s (14.15 m/s).

Figure 15.

(a) Back view of the composite panel after the test clamped (b) Zoomed in picture of the cracks formed on the back side of the panel.³ Closer examination shows a through crack.

Finite Element Modeling: The FE model is shown in Figure 16. Table 3 shows the details of the FE model. 8-node hexahedral elements were used. The panel is modeled using 16 elements through the thickness with a one-element layer for each ply. In between the plies, cohesive zone elements are used. There are 15 CZE layers modeled using LS-DYNA’s MAT_COHESIVE_GENERAL,³⁰ also known as MAT_186. MAT_186 can be used for modeling cohesive zone elements with an arbitrary traction-separation law in the form of tabulated data. Interaction between fracture modes I and II is accounted for via a damage formulation. The aluminum projectile was modeled using LS-DYNA’s MAT_PIECEWISE_LINEAR_PLASTICITY,³⁰ also known as MAT_024. MAT_024 is an elastic-plastic isotropic material model where the yield stress can be scaled as a function of strain rate. The nodes at the location of the bolts are restrained in-plane, and the nodes at the clamps are restrained in the out-of-plane direction.

Figure 16.

LVG1075 simulation showing impacted panel at the final time step with the projectile hidden from view (a) using PFC (b) using GTFC.

Table 3.

FE model details.

Number of elements			Number of nodes
MAT_213 (plies)	MAT_186 (between plies)	MAT_024 (projectile)	Number of nodes
373184	349860	17040	775238

Numerical calibration was used with the PFC model. The post-peak residual damage parameters used were: $PPR D_{T ∥}^{} = 0.64,$ $PPR D_{C ∥}^{} = 0.0,$ $PPR D_{T ⊥}^{} = 0.64,$ $PPR D_{C ⊥}^{} = 0.0,$ and $PPR D_{S}^{} = 0.64$ . The fracture energies were also calibrated, and the values used were: $Γ_{I}^{} = 400.0 l b / i n (7 \times 10^{4} N / m)$ and $Γ_{I I}^{} = 400.0 l b / i n (7 \times 10^{4} N / m)$ . Numerical calibration was also used with the GTFC model. It is assumed that the in-plane residual strengths in the 1-direction tension, 2-direction tension and 1-2 shear components are 30%. The through thickness residual strengths (3-direction tension, 2-3 shear and the 1-3 shear components) are assumed to be 34%. The equivalent failure strain has been assumed to be 0.8.

Results: Figure 16(a) shows the last frame from the LVG1075 simulation run with PFC and Figure 16(b) with GTFC. The simulations were terminated when the projectile rebound velocity became constant. Node identified as N402688 in Figure 16(b) is used to track and compare the out-of-plane displacement (Figure 17(a)). Figure 17(b) shows the projectile velocity as a function of time. It should be noted that (a) negative velocity indicates projectile rebound, and (b) the average rebound velocity of the projectile obtained from the experiment is shown as an average value.

Figure 17.

(a) Out-of-plane displacement at node 4,02,688 (b) Projectile velocity plotted against time. Experimental rebound velocity estimated and averaged over a period of 1.5 ms.

Discussions: As can be seen from Figure 16(a), there is widespread surface damage including elements eroded around the clamped region in the FE model contrary to the experimental results (Figure 12). The element erosion is due to the combination of both FF and IFF. The out-of-plane displacement is captured very well for the first 1 ms that includes both the first positive and the first negative peaks. It should be recognized that node 4,02,688 in the region with considerable surface damage starting around 1 ms and hence both the FE computed response and the experimental values are reliable only for that time duration. The calibration of the six failure-related parameters yields a flakier panel with a different damage and crack patterns and a higher rebound velocity. Perhaps, regression analysis involving these parameters may yield improved results, but the primary objective of the validation test is to show the difficulties in modeling impact tests without resorting to an involved parameter calibration process.

The GTFC results show qualitatively that the crack pattern obtained from the FE simulation reasonably matches the experiment. In addition, the crack is a through crack similar to the experiment. The maximum out-of-plane displacement prediction agrees well with the experimentally obtained value though there is difference in the wavelength of the displacement profile. The rebound velocity it overestimated similar to the PFC results. It should be noted that the calibration of the input residual strength was done with multiple objectives – matching crack pattern, out-of-plane displacement and rebound velocity, and it is very likely that decreasing the residual strengths in the tension and the shear components can lead to a less stiff panel response.

Concluding remarks

A three-dimensional general orthotropic material model involving elastic and inelastic deformations, damage and failure has been presented. The focus is on handling the challenges of a continuum level model in predicting when failure would take place in two loading scenarios - quasi-static loading, and high velocity impact events. Improvements were made to the post-peak behavior of the Puck Failure Criterion. These changes resulted in much better prediction for the staked-ply tension and compression validation tests. Minor calibration was needed to improve the results for the stacked-ply compressive test. Linking the post-peak behavior to fracture energy makes the strength predictions less dependent on the finite element mesh as evidenced from the convergence analyses. The impact test simulation required calibration of the post-peak residual damage parameters to obtain a reasonable prediction. The out-of-plane displacement is predicted very well for the first 1 ms following the impact. The rebound velocity is within 10% of the experimental result after normalization with respect to impact velocity. However, the damage and crack patterns are different. The material behavior appears to be flaky in comparison to the tested composite panel.

The Generalized Tabulated Failure Criterion is a new failure criterion implemented in MAT_213. It permits an arbitrary-shaped failure surface to be defined and used. In its current form, the failure surface is decomposed into in-plane and out-of-plane failure modes thereby making it suitable for modeling unidirectional composites. The stacked-ply tension test results show the peak load to be underpredicted, but the compressive test results are very close to the experimental results. Minimal calibration was done with the residual strength and failure strain values. The strength of GTFC is demonstrated in the impact test simulation. The rebound velocity is within 10% of the experimental result after being normalized with the impact velocity. The crack pattern qualitatively matches the experimental results. The out-of-plane displacement is close to the first positive peak, but appears to diverge from the experimental results after that.

It is well understood that damage and failure of composite systems is driven by the composite architecture and the constituent materials, and that while multi-scale modeling may make the numerical model more predictive, the enormous computational expense makes this approach untenable. What is desirable from an analyst’s viewpoint, is a reasonable and systematic approach that requires (a) input data that is largely experimental (physical and virtual) driven, (b) a small number of damage and failure-related parameters, and (c) an even smaller number of parameters that may require some tuning/calibration to yield consistent and somewhat conservative results. Further improvements that can be made to make the predictions more accurate include (a) using rate-dependent stress-strain data in the deformation sub-model, (b) making damage and failure rate-dependent, and (c) defining a richer failure surface for use in GTFC. Analysis of the impact FE simulations show that the computed strain-rate values are very high (in the order of 10³). There is ample experimental evidence showing that for unidirectional carbon/epoxy composites the strength in the in-plane shear and the transverse components increase with increase in strain-rate.³¹ These issues are being addressed in ongoing research work.

Footnotes

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: Authors Shyamsunder, Khaled, and Rajan gratefully acknowledge the support of the Federal Aviation Administration through Grant #12-G-001 titled “Composite Material Model for Impact Analysis” and #17-G-005 titled “Enhancing the Capabilities of MAT213 for Impact Analysis”, William Emmerling and Dan Cordasco, Technical Monitors.

ORCID iD

Subramaniam D Rajan

References

Kaddour

Hinton

MJ.

Failure criteria for composites. Amsterdam: Elsevier Ltd, 2018.

Sun

Chen

JL.

A simple flow rule for characterizing nonlinear behavior of fiber composites. J Compos Mater 1989; 23: 1009–1020.

Shyamsunder

Khaled

Rajan

, et al. Implementing deformation, damage and failure in an orthotropic plastic material model. J Compos Mater 2020; 54: 463–484.

Goldberg

Carney

DuBois

, et al. Development of an orthotropic elasto-plastic generalized composite material model suitable for impact problems. ASCE J Aerosp Eng 2016; 29: 04015083. DOI: 10.1061/(ASCE)AS.1943-5525.0000580.

Hoffarth

Rajan

Goldberg

, et al. Implementation and validation of a three-dimensional plasticity-based deformation model for orthotropic composites. Composites A 2016; 91: 336–350.

Hoffarth

Khaled

Shyamsunder

, et al. Verification and validation of a three-dimensional orthotropic plasticity constitutive model using a unidirectional composite. Fibers 2017; 5: 12–13.

Ansys-LST. LS-DYNA R11, www.lstc.com/products/ls-dyna (accessed 10 December 2019).

Goldberg

Carney

DuBois

, et al. Implementation of a tabulated failure model into a generalized composite material model. J Compos Mater 2018; 52: 3445–3460.

Hoffarth

Khaled

Shyamsunder

, et al. Development of a tabulated material model for composite material failure, MAT213 part 1: theory, implementation, verification & validation. Technical Report DOT/FAA/TC-19/50, P1, 2020.

10.

Goldberg

Carney

DuBois

, et al. Analysis and characterization of damage using a generalized composite material model suitable for impact problems. ASCE J Aerosp Eng 2018; 31: 04018025. DOI: 10.1061/(ASCE)AS.1943-5525.0000854.

11.

Khaled

Shyamsunder

Hoffarth

, et al. Damage characterization of composites to support an orthotropic plasticity model. J Compos Mater 2019; 53: 941–967.

12.

Puck

Schurmann

Failure analysis of FRP laminates by means of physically based phenomenological models. Compos Sci Technol 2002; 62: 1633–1662.

13.

Deuschle

Kroplin

BH.

Finite element implementation of puck's failure theory for fiber-reinforced composites under three-dimensional stress. J Compos Mater 2012; 46: 2485–2513.

14.

Deuschle

Puck

Application of the puck failure theory for fiber-reinforced composites under three-dimensional stress: comparison with experimental results. J Compos Mater 2013; 47: 827–846.

15.

Khaled

Shyamsunder

Hoffarth

, et al. Experimental characterization of composites to support an orthotropic plasticity material model. J Compos Mater 2018; 52: 1847–1872.

16.

Puck

Kopp

Knops

Guidelines for the determination of the parameters in puck’s action plane strength criterion. Compos Sci Technol 2002; 62: 371–378.

17.

Kodagali

Tessema

Kidane

Progressive failure analysis of a composite lamina using puck failure criteria. In: Thirty-Second technical conference, American Society for Composites 2017, West Lafayette, IN, 23–25 October 2017.

18.

Pinho

Robinson

Iannucci

Fracture toughness of the tensile and compressive fiber failure modes in laminated composites. Compos Sci Technol 2006; 66: 2069–2079.

19.

Khaled

Shyamsunder

Holt

, et al. Enhancing the predictive capabilities of a composite plasticity model using cohesive zone modeling. Composites A 2019; 121: 1–17.

20.

Hinton

Kaddour

and Soden

PD.

The world-wide failure exercise: its origin, concept and content. In: MJ Hinton, AS Kaddour and PD Soden (eds) failure criteria in fibre reinforced polymer composites: the world-wide failure exercise, Edited by MJ Hinton, AS Kaddour and PD Soden, Amsterdam: Elsevier, 2004.

21.

Kaddour

Hinton

MJ.

Benchmarking of triaxial failure criteria for composite laminates: comparison between models of part (A) of ‘WWFE-II. J Compos Mater 2012; 46: 2593–2632.

22.

Tsai

A general theory of strength for anisotropic materials. J Compos Mater 1971; 5: 58–80.

23.

Hashin

Failure criteria for unidirectional fiber composites. J Appl Mech 1980; 47: 329–334.

24.

Shyamsunder

Khaled

Rajan

SD.

MAT_213 v1.3.5 Status. In: FAA research project monthly teleconference, Tempe, AZ, 23 October 2019.

25.

Toray Carbon Fibers America, www.toraycma.com/ (2018, accessed 28 February 2018).

26.

Khaled

Shyamsunder

Schmidt

, et al. Development of a tabulated material model for composite material failure, MAT213 Part 2: experimental tests to characterize the behavior and properties of T800-F3900 Toray composite. Technical Report DOT/FAA/TC-19/50, P2, 2020.

27.

Holt

NT.

Experimental and simulation validation tests for MAT213, MS Thesis, School of Sustainable Engineering and the Built Environment, Arizona State University, 2018.

28.

Melis

Pereira

Goldberg

, et al. Dynamic impact testing and model development in support of NASA’s advanced composites program. In: Structures, structural dynamics and materials conference, AIAA 2018. West Lafayette, IN, 8–12 January 2018.

29.

GOM. ARAMIS Professional, www.gom.com/3d-software/gom-system-software/aramis-professional.html (accessed 18 December 2020).

30.

Ansys-LST. LS-DYNA keyword user’s manual volume II – material models. Livermore, California, USA: Author, 2020.

31.

Daniel

Werner

Fenner

JS.

Strain-rate-dependent failure criteria for composites. Compos Sci Technol 2011; 71: 357–364.