Verification and Validation of a Three-Dimensional Orthotropic Plasticity Constitutive Model Using a Unidirectional Composite

Abstract

A three-dimensional constitutive model has been developed for modeling orthotropic composites subject to impact loads. It has three distinct components—a deformation model involving elastic and plastic deformations; a damage model; and a failure model. The model is driven by tabular data that is generated either using laboratory tests or via virtual testing. A unidirectional composite—T800/F3900, commonly used in the aerospace industry, is used in the verification and validation tests. While the failure model is under development, these tests indicate that the implementation of the deformation and damage models in a commercial finite element program, LS-DYNA, is efficient, robust and accurate.

1. Introduction

Modeling and simulation of composites during an impact event is a huge challenge given the wide variety of composites and the associated complexity in characterization of the behavior of the constituent materials as well as the interaction between them. While composites have been in use for decades in a variety of industries such as civil structures, automotive and aerospace, building a predictive model is still daunting. Kaddour et al. [1] list some of the challenges facing the industry. They include (a) shorter life cycles, (b) automated manufacture, (c) production of high volumes, (d) integrating 3D structures into 3D architectures, (e) development of alternative materials, and (f) meeting climate change targets. In the United States, several governmental agencies (including NASA and the FAA) have recognized the importance of building a framework for composites system by forming a public–private consortium. A press release [2] states that “NASA formed the consortium in support of the Advanced Composites Project, which is part of the Advanced Air Vehicles Program in the agency’s Aeronautics Research Mission Directorate. The project’s goal is to reduce product development and certification timelines by 30 percent for composites infused into aeronautics applications.” A major reason for these challenges is the lack of mature material models that should be able to predict, with some degree of certainty, the deformation, damage and failure of composite systems.

Several approaches can be used to develop such numerical models. Vaziri et al. [3] categorize the continuum mechanics approaches into four groups—(i) nonlinear elasticity theories; (ii) damage theories coupled with elasticity; (iii) classical incremental plasticity theory; and (iv) endochronic plasticity theory. Examples of these approaches and models can be found in [4,5,6,7,8,9]. From the viewpoint of engineering design, models should be driven with varying levels of input depending on the complexity of the composite so that the material model can be tailored for the required response. Currently, however, many material models are driven entirely by a few sets of material properties capturing orthotropic deformation and damage behavior. These include Young’s moduli, Poisson’s ratios, shear moduli, failure strain, and reduction factors for the strength of constituent parts (fibers) based on the failure of other constituent parts (matrix). Others require additional parameters to characterize failure using classical failure theories at the macroscopic level, e.g., Chang–Chang or Tsai–Wu, or a mixture of macroscopic level and local (fiber/matrix) driven failure, e.g., Hashin. This type of approach leads either to models with only a few parameters, which provide a crude approximation at best to the actual stress–strain curve, or to models with many parameters which require a large number of complex tests to characterize [12]. Often, the lack of flexibility or sophistication in these models require that designers adjust the values of some of the model parameters using data from full-scale structural tests.

In the present approach, the composite impact model is divided into a deformation model, a damage model and a failure model (Figure 1). Such a partitioning allows for the elastic and plastic deformations to be captured by the deformation model and the reduction in stiffness to be captured by the damage model. In other words, the deformation model simulates the nonlinear material response of the composite (due to either deformation or damage mechanisms), and the damage model simulates the nonlinear unloading due to stiffness reduction. The macroscopic tabulated failure model captures failure and does not rely on characterization of mechanisms leading to failure. The response of a system can then be monitored as the finite element calculations are processed through the deformation and damage models so that the failure model can be used to carry out failure predictions. The deformation model generalizes the Tsai–Wu failure criteria. A non-associative flow rule is used and the yield surface is based on strain-hardening orthotropic yield function. A strain equivalent formulation is utilized in the damage model that permits the plasticity and damage calculations to be uncoupled and captures the nonlinear unloading and local softening of the stress–strain response. A diagonal damage tensor is defined to account for the directional-dependent variation of damage. However, in composites, it has been found that loading in one direction can lead to damage in multiple coordinate directions [19]. To account for this phenomenon, the terms in the damage matrix are semi-coupled as explained later. The overall framework is driven by experimentally obtained tabulated temperature- and rate-dependent stress–strain data as well as data that characterizes the damage matrix and failure [13].

This paper discusses the details of the orthotropic elasto-plastic constitutive model. The current version supports solid finite elements with plans to extend the support to shell elements. The theoretical and implementation details are followed by numerical studies using a custom version of LS-DYNA where the constitutive model is implemented as MAT213. A commonly used unidirectional laminate, T800/F3900 graphite/epoxy composite, is used in the numerical studies.

2. Theoretical Details

The developed constitutive model is designed to model composites under impact loading. The input to the model consists of 12 tabulated stress–strain curves generated from actual experimental or virtual tests [10]—1, 2 and 3-direction tension curves, 1, 2 and 3-direction compression curves, 1-2, 2-3 and 3-1 plane shear curves, and 45° off-axis tension/compression curves in the 1-2, 2-3 and 3-1 planes. A brief summary of the deformation model is presented here and further details can be found in our earlier publications [11,12,13]. The yield function is a modified form of the commonly used Tsai–Wu composite failure model and is defined as

$$f(\sigma )=a+\left(\begin{array}{cccc}{F}_1& F_2& F_3& \begin{array}{ccc}0& 0& 0\end{array}\end{array}\right)\left[\begin{array}{c}{\sigma }_{11}\\ {\sigma }_{22}\\ {\sigma }_{33}\\ \begin{array}{c}{\sigma }_{12}\\ {\sigma }_{23}\\ {\sigma }_{31}\end{array}\end{array}\right]+{\left[\begin{array}{c}{\sigma }_{11}\\ {\sigma }_{22}\\ {\sigma }_{33}\\ \begin{array}{c}{\sigma }_{12}\\ {\sigma }_{23}\\ {\sigma }_{31}\end{array}\end{array}\right]}^T\left[\begin{array}{cccc}\begin{array}{c}{F}_{11}\\ F_{12}\\ F_{13}\\ \begin{array}{c}0\\ 0\\ 0\end{array}\end{array}& \begin{array}{c}{F}_{12}\\ F_{22}\\ F_{23}\\ \begin{array}{c}0\\ 0\\ 0\end{array}\end{array}& \begin{array}{c}{F}_{13}\\ F_{23}\\ F_{33}\\ \begin{array}{c}0\\ 0\\ 0\end{array}\end{array}& \begin{array}{c}\begin{array}{ccc}0& 0& 0\end{array}\\ \begin{array}{ccc}0& 0& 0\end{array}\\ \begin{array}{ccc}0& 0& 0\end{array}\\ \begin{array}{c}\begin{array}{ccc}{F}_{44}& 0& 0\end{array}\\ \begin{array}{ccc}0& F_{55}& 0\end{array}\\ \begin{array}{ccc}0& 0& F_{66}\end{array}\end{array}\end{array}\end{array}\right]\left[\begin{array}{c}{\sigma }_{11}\\ {\sigma }_{22}\\ {\sigma }_{33}\\ \begin{array}{c}{\sigma }_{12}\\ {\sigma }_{23}\\ {\sigma }_{31}\end{array}\end{array}\right],$$

(1)

where $a=-1$. The yield function coefficients, $F_{ij}$, depend on the current yield stress values and are calculated as

$$\begin{array}{ccccc}\begin{array}{l}{F}_1=\frac{1}{{\sigma }_{11}^T}-\frac{1}{{\sigma }_{11}^C}\\ F_2=\frac{1}{{\sigma }_{22}^T}-\frac{1}{{\sigma }_{22}^C}\\ F_3=\frac{1}{{\sigma }_{33}^T}-\frac{1}{{\sigma }_{33}^C}\end{array}& & \begin{array}{l}{F}_{11}=\frac{1}{{\sigma }_{11}^T{\sigma }_{11}^C}\\ F_{22}=\frac{1}{{\sigma }_{22}^T{\sigma }_{22}^C}\\ F_{33}=\frac{1}{{\sigma }_{33}^T{\sigma }_{33}^C}\end{array}& & \begin{array}{l}{F}_{44}=\frac{1}{{\sigma }_{12}^2}\\ F_{55}=\frac{1}{{\sigma }_{23}^2}\\ F_{66}=\frac{1}{{\sigma }_{31}^2}\end{array}\end{array}$$

(2)

$$F_{12}=\frac{2}{{\left({\sigma }_{12}^{45}\right)}^2}-\frac{F_1+F_2}{{\sigma }_{12}^{45}}-\frac{1}{2}\left(F_{11}+F_{22}+F_{44}\right)$$

(3)

$$F_{23}=\frac{2}{{\left({\sigma }_{23}^{45}\right)}^2}-\frac{F_2+F_3}{{\sigma }_{23}^{45}}-\frac{1}{2}\left(F_{22}+F_{33}+F_{55}\right)$$

(4)

$$F_{13}=\frac{2}{{\left({\sigma }_{31}^{45}\right)}^2}-\frac{F_1+F_3}{{\sigma }_{31}^{45}}-\frac{1}{2}\left(F_{11}+F_{33}+F_{66}\right)$$

(5)

where the superscripts T, C and 45 denote data obtained from tension, compression and 45-degree off-axis tests, respectively. A non-associative flow rule is used to compute the evolution of the components of plastic strain and the plastic potential function is defined as

$$h=\sqrt{H_{11}{\sigma }_{11}^2+H_{22}{\sigma }_{22}^2+H_{33}{\sigma }_{33}^2+2H_{12}{\sigma }_{11}{\sigma }_{22}+2H_{23}{\sigma }_{22}{\sigma }_{33}+2H_{31}{\sigma }_{33}{\sigma }_{11}+H_{44}{\sigma }_{12}^2+H_{55}{\sigma }_{23}^2+H_{66}{\sigma }_{31}^2}$$

(6)

where the $H_{ij}$ terms are a set of constant coefficients which are functions of the plastic Poisson ratios, with the coefficients defined as input parameters in the model.

2.1. Non-Associated Plastic Flow

Classical plasticity theory using associated plasticity assumes that the increment of plastic strain is normal to the yield surface, where the plastic potential function is defined as

$${\dot{\mathsf{\epsilon }}}^p=\dot{\lambda }\frac{\partial f}{\partial \mathsf{\sigma }}$$

(7)

based on Drucker’s stability postulate which works well for metals [14] (where ${\dot{\mathsf{\epsilon }}}^p$ is the plastic strain rate, $\dot{\lambda }$ is the plastic multiplier, $f$ is the yield function and $\sigma$ is the stress tensor in Voigt notation). However, associated plasticity is not ideal for most composites with various degrees of plastic anisotropy. For example, some composites exhibit linear elastic behavior in one or more directions, e.g., a unidirectional composite when loaded in the fiber direction. This would indicate that there is no plastic flow, or plastic strain accumulation with respect to the fiber direction. However, the yield function in Equation (1) cannot accommodate this. Thus, non-associated plasticity is required in creating a generalized composite material model as the flow law coefficients for the non-associated plastic potential function can be determined independently from the yield function terms. Mathematically, this situation (that associated plasticity cannot be used) can be described as follows. Consider an isotropic Tsai-Wu flow rule:

$$f={\sigma }_{vm}^2-A_{1}p-A_2{p}^2$$

(8)

where ${\sigma }_{vm}$ is the von Mises stress and $p$ is the pressure. The condition for the uniaxial tension case can be derived as follows ($s_{xx}$ is the deviatoric stress):

$${\dot{\epsilon }}_{xx}^p\approx \frac{\partial f}{\partial {\sigma }_{xx}}=2{\sigma }_{vm}\frac{3s_{xx}}{2{\sigma }_{vm}}-A_1\left(-\frac{1}{3}\right)-2p{A}_2\left(-\frac{1}{3}\right)=3s_{xx}+\frac{A_1}{3}+\frac{2p{A}_2}{3}$$

(9)

$${\dot{\epsilon }}_{vol}^p\approx \frac{\partial f}{\partial {\sigma }_{xx}}+\frac{\partial f}{\partial {\sigma }_{yy}}+\frac{\partial f}{\partial {\sigma }_{zz}}=A_1+2p{A}_2$$

(10)

$${\dot{\epsilon }}_{vol}^p=\left(1-2{\nu }_p\right){\dot{\epsilon }}_{xx}^p\mathrm{and}{s}_{xx}={\sigma }_{xx}+p=\frac{2}{3}{\sigma }_{xx}=-2p$$

(11)

$$A_1+2p{A}_2=\frac{\left(1-2{\nu }_p\right)}{3}\left[9s_{xx}+A_1+2p{A}_2\right]$$

(12)

Assuming that the linear term is zero, the plastic Poisson’s ratio can be derived as

$$A_1=0\Rightarrow 2A_2=\frac{\left(1-2{\nu }_p\right)}{3}\left[-18+2A_2\right]\Rightarrow {\nu }_p=\frac{9+2A_2}{18-2A_2}$$

(13)

Equation (13) shows that the plastic Poisson’s ratio can assume any value which is non-physical. Hence a non-associative flow rule is appropriate for the problem under consideration in this paper and is given by

$${\dot{\mathsf{\epsilon }}}^p=\dot{\lambda }\frac{\partial h}{\partial \sigma }$$

(14)

where $h$ is the plastic potential function.

2.2. Flow Rule and Determination of Flow Rule Coefficients

The deformation model uses a flow surface that is defined to determine the plastic flow direction (flow rule) as shown in Equation (14). In the plastic potential function (Equation (6)), the $H_{ij}$ terms are the flow rule coefficients, assumed to be constant, and ${\sigma }_{jj}$ are the current stress values. To accurately represent the plastic flow of the material, the flow rule coefficients are calculated directly using the experimental data. Six of the coefficients are determined from unidirectional testing of the composite in the three principal material directions. In each case, the plastic potential function can be simplified resulting in set of functions from which the plastic Poisson’s ratios, ${\nu }_{ij}^p$ can be computed as

$$\begin{array}{ccc}\begin{array}{l}{\sigma }_{11}\ne 0\\ {\nu }_{12}^p=-\frac{{\dot{\epsilon }}_{22}^p}{{\dot{\epsilon }}_{11}^p}=-\frac{H_{12}}{H_{11}}\\ {\nu }_{13}^p=-\frac{{\dot{\epsilon }}_{33}^p}{{\dot{\epsilon }}_{11}^p}=-\frac{H_{13}}{H_{11}}\end{array}& \begin{array}{l}{\sigma }_{22}\ne 0\\ {\nu }_{21}^p=-\frac{{\dot{\epsilon }}_{11}^p}{{\dot{\epsilon }}_{22}^p}=-\frac{H_{12}}{H_{22}}\\ {\nu }_{23}^p=-\frac{{\dot{\epsilon }}_{33}^p}{{\dot{\epsilon }}_{22}^p}=-\frac{H_{23}}{H_{22}}\end{array}& \begin{array}{l}{\sigma }_{33}\ne 0\\ {\nu }_{32}^p=-\frac{{\dot{\epsilon }}_{22}^p}{{\dot{\epsilon }}_{33}^p}=-\frac{H_{23}}{H_{33}}\\ {\nu }_{31}^p=-\frac{{\dot{\epsilon }}_{11}^p}{{\dot{\epsilon }}_{33}^p}=-\frac{H_{13}}{H_{33}}\end{array}\end{array}$$

(15)

where ${\dot{\mathsf{\epsilon }}}_{ij}^p$ are the plastic strain rate components. It should be noted that since the plastic strains and their ratios are constantly changing, the average values are used in computing the coefficients. An examination of the ratios (the ratios compare the degree of plasticity in each material direction) in Equation (15) shows that the flow rule coefficients are not linearly independent. This requires that one of the flow rule coefficients be assumed. Usually the direction with the most plasticity is chosen and the flow rule coefficient value is assumed, e.g., $H_{ij}=1.0$. The stress–strain curve in that direction is tagged as the “master” curve. As an example, consider $H_{11}=1.0$. This results in the following simplified equations.

$$\begin{array}{l}{H}_{12}=-{\nu }_{12}^p{H}_{11}=-{\nu }_{12}^p\\ H_{13}=-{\nu }_{13}^p{H}_{11}=-{\nu }_{13}^p\\ H_{22}=\frac{{\nu }_{12}^p}{{\nu }_{21}^p}{H}_{11}=\frac{{\nu }_{12}^p}{{\nu }_{21}^p}\\ H_{23}=\frac{-{\nu }_{23}^p{\nu }_{12}^p}{{\nu }_{21}^p}{H}_{11}=\frac{-{\nu }_{23}^p{\nu }_{12}^p}{{\nu }_{21}^p}\\ H_{33}=\frac{{\nu }_{13}^p}{{\nu }_{31}^p}{H}_{11}=\frac{{\nu }_{13}^p}{{\nu }_{31}^p}\end{array}$$

(16)

Similarly, assuming the 2-direction as the master curve, $H_{22}=1.0$, and

$$\begin{array}{l}{H}_{11}=\frac{{\nu }_{21}^p}{{\nu }_{12}^p}{H}_{22}=\frac{{\nu }_{21}^p}{{\nu }_{12}^p}\\ H_{12}=-{\nu }_{21}^p{H}_{22}=-{\nu }_{21}^p\\ H_{13}=\frac{-{\nu }_{13}^p{\nu }_{21}^p}{{\nu }_{12}^p}{H}_{22}=\frac{-{\nu }_{13}^p{\nu }_{21}^p}{{\nu }_{12}^p}\\ H_{23}=-{\nu }_{23}^p{H}_{22}=-{\nu }_{23}^p\\ H_{33}=\frac{{\nu }_{13}^p{\nu }_{21}^p}{{\nu }_{31}^p{\nu }_{12}^p}{H}_{22}=\frac{{\nu }_{13}^p{\nu }_{21}^p}{{\nu }_{31}^p{\nu }_{12}^p}\end{array}$$

(17)

The first six flow rule coefficients are calculated with the experimentally obtained plastic Poisson’s ratios depending on the identified master curve. The final three terms –$H_{44}$, $H_{55}$ and $H_{66}$ must be calculated differently. These terms are calculated by fitting the shear curves with the master curve in the effective stress versus effective plastic strain space. This can be achieved by minimizing the difference between the master curve and shear curves, thus determining the optimal H value. The optimal flow rule coefficient value, $H_{ll}^{*}$, is obtained by computing the difference between the two curves and the problem is posed as a least-squares problem:

$$f({\mathrm{H}}_{ll})={{\displaystyle \sum _{k=1}^n\left[{\left({\widehat{\sigma }}_{ii}\right)}_k-{\left({\widehat{\sigma }}_{ij}\right)}_k\right]}}^2$$

(18)

such that

$$H_{ll}^{\min }\le H_{ll}^{*}\le H_{ll}^{\max }$$

(19)

where n is the number of data points in the master curve, ${\left({\widehat{\sigma }}_{ii}\right)}_k$ is the kth effective stress value from the master curve and ${\left({\widehat{\sigma }}_{ij}\right)}_k$ is the effective stress value from the shear curve. An alternative process for calculating these final three coefficients is available as fitting the off-axis curves to the master curve instead of the shear curves [15].

2.3. Modeling Damage

Damage in composite material models is usually used to track the progression of micro-cracking. The continuum damage model developed by Matzenmiller et al. [16] uses the initiation and accumulation of damage to model the nonlinear response of a composite by implementing damage parameters to reduce the effective stiffness of the material in a particular coordinate direction. The full damage tensor can be used to relate the true stress to the effective stress as

$$\left(\begin{array}{c}{\sigma }_{11}\\ {\sigma }_{22}\\ {\sigma }_{33}\\ {\sigma }_{12}\\ {\sigma }_{23}\\ {\sigma }_{13}\end{array}\right)=\left[\begin{array}{cccccc}{M}_{11}& M_{12}& M_{13}& M_{14}& M_{15}& M_{16}\\ M_{21}& M_{22}& M_{23}& M_{24}& M_{25}& M_{26}\\ M_{31}& M_{32}& M_{33}& M_{34}& M_{35}& M_{36}\\ M_{41}& M_{42}& M_{43}& M_{44}& M_{45}& M_{46}\\ M_{51}& M_{52}& M_{53}& M_{54}& M_{55}& M_{56}\\ M_{61}& M_{62}& M_{63}& M_{64}& M_{65}& M_{66}\end{array}\right]\left(\begin{array}{c}{\sigma }_{11}^{eff}\\ {\sigma }_{22}^{eff}\\ {\sigma }_{33}^{eff}\\ {\sigma }_{12}^{eff}\\ {\sigma }_{23}^{eff}\\ {\sigma }_{13}^{eff}\end{array}\right)$$

(20)

However, this full matrix is rarely employed. Typically, damage models use a diagonal damage tensor and assume that the components are a function of a single plastic strain, e.g., $M_{11}=M_{11}\left({\epsilon }_p^{11}\right)$. However, in the current work, a semi-coupled, directional-dependent relationship, e.g., $M_{11}=M_{11}\left({\epsilon }_p^{11},{\epsilon }_p^{22},{\epsilon }_p^{12}\right)$ indicating that the components are functions of plastic strains in multiple directions, is used to construct the modified relationship as

$$\left(\begin{array}{c}{\sigma }_{11}\\ {\sigma }_{22}\\ {\sigma }_{33}\\ {\sigma }_{12}\\ {\sigma }_{23}\\ {\sigma }_{13}\end{array}\right)=\left[\begin{array}{cccccc}{M}_{11}& 0& 0& 0& 0& 0\\ 0& M_{22}& 0& 0& 0& 0\\ 0& 0& M_{33}& 0& 0& 0\\ 0& 0& 0& M_{44}& 0& 0\\ 0& 0& 0& 0& M_{55}& 0\\ 0& 0& 0& 0& 0& M_{66}\end{array}\right]\left(\begin{array}{c}{\sigma }_{11}^{eff}\\ {\sigma }_{22}^{eff}\\ {\sigma }_{33}^{eff}\\ {\sigma }_{12}^{eff}\\ {\sigma }_{23}^{eff}\\ {\sigma }_{13}^{eff}\end{array}\right)$$

(21)

The damage parameter $d_{ij}^{kl}$ is defined as a function of the plastic strain and represents damage in kl due to loading along ij. By assuming that normal damage is due to all normal and shear terms, and shear damage is also due to all normal and shear terms, the damage model terms in Equation (21) can be expressed explicitly as

$$\begin{array}{l}{M}_{11}=\left(1-d_{11}^{11}\left({\epsilon }_e^p\right)\right)\left(1-d_{22}^{11}\left({\epsilon }_e^p\right)\right)\left(1-d_{33}^{11}\left({\epsilon }_e^p\right)\right)\left(1-d_{12}^{11}\left({\epsilon }_e^p\right)\right)\left(1-d_{23}^{11}\left({\epsilon }_e^p\right)\right)\left(1-d_{13}^{11}\left({\epsilon }_e^p\right)\right)\\ M_{22}=\left(1-d_{11}^{22}\left({\epsilon }_e^p\right)\right)\left(1-d_{22}^{22}\left({\epsilon }_e^p\right)\right)\left(1-d_{33}^{22}\left({\epsilon }_e^p\right)\right)\left(1-d_{12}^{22}\left({\epsilon }_e^p\right)\right)\left(1-d_{23}^{22}\left({\epsilon }_e^p\right)\right)\left(1-d_{13}^{22}\left({\epsilon }_e^p\right)\right)\\ M_{33}=\left(1-d_{11}^{33}\left({\epsilon }_e^p\right)\right)\left(1-d_{22}^{33}\left({\epsilon }_e^p\right)\right)\left(1-d_{33}^{33}\left({\epsilon }_e^p\right)\right)\left(1-d_{12}^{33}\left({\epsilon }_e^p\right)\right)\left(1-d_{23}^{33}\left({\epsilon }_e^p\right)\right)\left(1-d_{13}^{33}\left({\epsilon }_e^p\right)\right)\\ M_{44}=\left(1-d_{11}^{12}\left({\epsilon }_e^p\right)\right)\left(1-d_{22}^{12}\left({\epsilon }_e^p\right)\right)\left(1-d_{33}^{12}\left({\epsilon }_e^p\right)\right)\left(1-d_{12}^{12}\left({\epsilon }_e^p\right)\right)\left(1-d_{23}^{12}\left({\epsilon }_e^p\right)\right)\left(1-d_{13}^{12}\left({\epsilon }_e^p\right)\right)\\ M_{55}=\left(1-d_{11}^{23}\left({\epsilon }_e^p\right)\right)\left(1-d_{22}^{23}\left({\epsilon }_e^p\right)\right)\left(1-d_{33}^{23}\left({\epsilon }_e^p\right)\right)\left(1-d_{12}^{23}\left({\epsilon }_e^p\right)\right)\left(1-d_{23}^{23}\left({\epsilon }_e^p\right)\right)\left(1-d_{13}^{23}\left({\epsilon }_e^p\right)\right)\\ M_{66}=\left(1-d_{11}^{13}\left({\epsilon }_e^p\right)\right)\left(1-d_{22}^{13}\left({\epsilon }_e^p\right)\right)\left(1-d_{33}^{13}\left({\epsilon }_e^p\right)\right)\left(1-d_{12}^{13}\left({\epsilon }_e^p\right)\right)\left(1-d_{23}^{13}\left({\epsilon }_e^p\right)\right)\left(1-d_{13}^{13}\left({\epsilon }_e^p\right)\right)\end{array}$$

(22)

Since even the simplified form shown in Equations (21)–(22) is rather complex, one must carry out damage related testing to see how many of these terms are really necessary to characterize the damage of a specific composite. Determining the damage parameters requires the measurement of either the damaged modulus or plastic strain (by unloading the material to a state of zero stress) at several total strain values. One approach to using the experimental data is to compute the damaged modulus as shown in Equation (23) for the normal stress–strain relationship in the 11 direction.

$$\begin{array}{l}{d}_{11}^{11}=1-\frac{E_{11}^{d11}}{E_{11}}\\ M_{11}^{11}=1-d_{11}^{11}\\ {\sigma }_{11}^{eff}=\frac{{\sigma }_{11}}{M_{11}^{11}}\\ {\epsilon }_{11}^p={\epsilon }_{11}-\frac{{\sigma }_{11}^{eff}}{E_{11}}\end{array}$$

(23)

where $E_{11}^{d11}$ is the damaged modulus of elasticity in the 1-direction, and $E_{11}$ is the undamaged modulus of elasticity. This approach makes it possible to create a simple damage wrapper around the deformation model and details can be found in [13,19].

The theory as described in this section, is implemented in a special version of LS-DYNA as material model 213 (MAT213) and is used in the numerical studies.

3. Experimental and Numerical Modeling Studies

An impact test (LVG905) was conducted at National Aeronautics and Space Administration’s Glenn Research Center to validate the developed model [12]. The test involves a 30.48 cm × 30.48 cm × 0.3099 cm (12 in × 12 in × 0.122 in) thick composite flat panel subjected to a low velocity impact. The panel is made of T800S/3900-2B[P2352W-19] BMS8-276 Rev-H-Unitape fiber/resin unidirectional composite [17]. Toray describes T800S as an intermediate modulus, high tensile strength graphite fiber. The epoxy resin system is labeled F3900 where a toughened epoxy is combined with small elastomeric particles to form a compliant interface or interleaf between fiber plies to resist impact damage and delamination [17]. Magnified views of the composite are shown in Figure 2 and Figure 3.

The schematics of the impact test are shown in Figure 4a,b. Figure 4c shows the 50.8 gm (0.112 lbm) hollow Al-2024 projectile with radiused front face. Figure 4d shows the drawing of the panel and the support plate. The unidirectional panel was fabricated with 16 plies with the fibers in the panel oriented vertically corresponding to the Y-direction in Figure 5a. Figure 4e,f shows the panels before and after impact (no visible damage and cracks were observed). Each panel was tested before the impact test using ultrasound techniques (pulse echo scans with the sample immersed in water, a 5 MHz focused transducer was used using 1 mm scan and step increment) to determine if there were any initial damage or flaws. The tests revealed no flaws or damage. Digital Image Correlation (DIC) data was gathered from high-speed cameras and analyzed to yield full-field displacements and strains. It should be noted that the impact velocity was chosen to be low enough to cause deformation with little or no damage, and certainly no localized failure. After the tests, the panels were carefully investigated and the examinations showed no damage and failure.

In the current version of MAT213 implementation, (a) the deformation model requires 12 stress–strain curves for each set of temperatures and strain rates: 1, 2 and 3-direction tension curves, 1, 2 and 3-direction compression curves, 1-2, 2-3 and 3-1 plane shear curves, and 45° off-axis tension/compression curves in the 1-2, 2-3 and 3-1 planes [10,11,12]; (b) the damage model requires $d_{ij}^{kl}$ versus total strain curves [18,19]; and (c) a simple failure model (principal strain theory) is currently implemented with the intent to flag models where one or more strain components has exceeded limiting values, rather than to actually predict failure. For the impact test under consideration, quasi-static, room temperature data was used.

The 50.8 gm (0.112 lbm) projectile (Figure 4c) was fired at the panel at a velocity of 710.184 cm/s (23.3 ft/s). The low velocity impact was designed to cause elastic and inelastic deformations with little or no damage. Examination of the panel before (Figure 4e) and after impact (Figure 4f) showed no visible damage or cracks. Full-field displacements and strain data on the back face of the composite panel were obtained using digital image correlation (DIC) analysis of high-speed images.

The special version of LS-DYNA containing the material model (MAT213) was used for an explicit finite element (FE) analysis. To compare MAT123 with an existing composite material model in LS-DYNA, MAT22 was used (see

Table 1. Composite Material Properties used in MAT22 FE Model.

Model Parameter	Value
Mass density (g/mm3)	1.55 (10−3)
Ea (GPa)	1.51 (102)
Eb (GPa)	7.89
Ec (GPa)	7.89
νba	0.01385
νbc	0.01385
νca	0.3792
Gab (GPa)	4.00
Gbc (GPa)	2.24
Gca (GPa)	4.00
Shear Strength (MPa)	94.87
Longitudinal Tensile Strength (GPa)	2.84
Transverse Tensile Strength (GPa)	0.060
Transverse Compressive Strength (GPa)	0.168
Alpha	0.0
Normal Tensile Strength (MPa)	60.12
Transverse Shear Strength (MPa)	1.03 (102)
Transverse Shear Strength (MPa)	94.87

). Figure 5 shows various views and details of the finite element (FE) model. The boundary conditions were applied to the plate in a way that mimicked the manner the plate was supported in the test frame. The bolted assembly shown in Figure 4d was modeled by fixing the X, Y, Z translational displacements of the nodes in the gripping region of the panel. A detailed FE model convergence analysis was carried out to find the most efficient and accurate model—the optimal mesh. This optimal mesh contained 288,000 8-noded hexahedral elements with one-integration point and hourglass control type 6. A typical element is 0.127 cm × 0.127 cm × 0.0635 cm (0.05 in × 0.05 in × 0.025 in) with five elements through the thickness of the panel (the thickness of each element is roughly equal to the thickness of three plies). The aluminum impactor was modeled with 27,200 8-noded hexahedral elements. The aluminum impactor was modeled using a piecewise linear plasticity model (MAT24) with the material properties given in [12]. No damping parameters were used in the FE models.

The runs were made on a 4-node cluster—4 Dell Precision Workstations (2 are Dell T3600 with Intel E5-1607 Xeon 3 GHz with 20 GB RAM, 2 are Dell T3610 with Intel E5-1607 Xeon v2 3 GHz with 20 GB RAM) (Dell Corporation, Round Rock, TX, USA). All CPUs are quad-core. The four nodes are connected to each other using a Dell Gigabit switch. The operating system is Red Hat 6.7 and the LS-DYNA version uses HP-MPI library.

Two comparison metrics were used—the maximum out-of-plane (Z-direction) displacement and the contour of the out-of-plane displacements, both on the back face of the plate. The contour plots are shown in Figure 6 for t = 0.0007 s. The MAT213 results (Figure 6b) are very similar to that of impact test (Figure 6a), with a rounded shape that is slightly elongated in the fiber direction especially below the center of impact (COI).

Figure 7 shows the out-of-plane displacement versus time results for the test and the simulations. The results are discussed in the next section.

4. Discussion

MAT213 results show that both the positive peak and the negative peak displacement values are close to the test values for the first cycle (Figure 7). Thereafter, the peak displacements are overestimated. In addition, there is a difference in the period. Nevertheless, the maximum positive and negative displacements are within a few percent of the experimental values. MAT22, on the other hand, over-predicts the maximum positive displacement (by about 35%) and under-predicts the maximum negative displacement (by about 80%). The displacement contours for MAT213 (Figure 6b) are very similar to that of the experiment (Figure 6a)—an ellipsoidal shape with the major axis in the fiber direction. MAT22 (Figure 6c) shows a more pronounced elongated distribution of displacements in the fiber direction. It should be noted that MAT22 has been developed for use with composites that essentially exhibit linear behavior with sudden, brittle failure.

Finally, it should be noted that there are a few differences between the impact test and finite element models. First, in the impact test, the projectile impact was not a direct hit, i.e., the roll, pitch and yaw angles were not all zero. However, zero roll, pitch and yaw angles were assumed in the finite element models. While the desired center of impact was the center of the panel, due to the low velocity, the center of impact was below the center of the panel (Figure 6a); second, it is likely that the period differences in the test and FE models are partly due to no damping parameters being used in the FE models.

Figure 8 shows the energy checks are used to gage the quality of the finite element solution. The total kinetic energy and the total internal energy appear to vary reasonable well with time. The hourglass energy values are relatively small (less than 10% of the total energy).

The total wall clock time for the MAT213 run was about 13 h and 38 min and for the MAT22 run 2 h and 3 min.

5. Conclusions

A new material model and computational framework have been developed for modeling orthotropic composites subject to impact loads. The algorithm has been implemented in a special version of LS-DYNA as MAT213 material model. The deformation model is driven by a set of 12 tabular data representing stress–strain curves for one temperature and strain rate combination. A semi-coupled, directional-dependent damage tensor is used to capture the effects of damage. The failure model is under development. The behavior of a unidirectional composite, T800/F3900, is used to verify and validate the developed material model. Preliminary results from modeling an impact problem indicate that the implementation is robust and accurate.

Several efforts are underway to improve the material model. First, a failure model is being developed so that element failure can be predicted and if necessary, the element can be eroded from the FE model. This will allow for modeling of a wider array of actual impact tests where deformation, damage and failure take place. Second, the code implementation has not been optimized for speed. It is expected that significant increase in speed would be possible with code optimization (vectorization, loop unrolling etc.). Third, the code has not been validated against impact tests with sufficiently high velocity where rate, and possibly temperature-dependent material behavior, is an important part of the FE behavioral model. While MAT213 has the capabilities to model these scenarios, experimental data is not currently available and the research team is exploring ways to conduct material characterization and impact tests to yield the data.