A nonlinear strain rate and pressure-dependent micro-mechanical composite material model for impact problems

Abstract

A micro-mechanical composite material model is developed to simulate the behavior of unidirectional composites under impact loading conditions in the nonlinear finite element solver (LS-DYNA^®). The nonlinear strain rate and pressure dependency in the composite material model is accounted by the resin, which uses previously developed state variable viscoplastic equations. These equations have been originally developed for metals; however, these are modified to account for the significant contributions of hydrostatic stresses typically observed in polymers. The material model also uses a continuum damage mechanics (CDM) based failure model to incorporate the progressive post-failure behavior. A set of Weibull distribution functions are used to quantify this behavior, and a methodology of assigning physical significance to the choice of damage/softening parameters used in these functions is presented. The impact response of composite laminate plates has been simulated and compared to the experiments. It has been observed that the predicted results compare favorably to the experiments.

Introduction

In the past few decades, laminated composites have become the most favored material system for a wide range of industrial applications. Factors such as potential for structural tailoring, superior strength and stiffness for a given weight along with higher toughness and mechanical damping, and so on, have justified their increased use. Due to this high usage, the behavior of laminated composites under impact remains a concern, mainly as impacts are unavoidable during manufacturing, day-to-day usage, and maintenance.

Predicting an accurate behavior of the laminated composite, especially for impacts which induce a significant amount of internal damage, is a significant challenge reported by Collombert et al. ¹ and reviewed by Abrate. ^2,3 With the advancements made in computing technology, finite element (FE) simulations have become a substantial tool to study impact problems, especially with the use of an accurate material model as detailed information on the spatial and temporal distribution of damage can be made available during impact. ⁴ The popular transient nonlinear FE code LS-DYNA has emerged as a powerful tool to carry out such analysis, as it offers a fast and simple user interface to implement own complex material models.

In numerous experiments over the years, many researchers have observed that the modulus and strength of the composites increases with strain rate. ^{5 –7} This dependency of a composite response on strain rate is important in impact simulations. ¹⁷  Furthermore, many experimental works ^{9 –12} have shown that composite sensitivity to strain rate is usually driven by the resin behavior. Hence, to accurately capture the response of a composite lamina, there is a need to accurately represent the nonlinear behavior of the polymer constituent.

Polymeric materials have been typically observed to exhibit nonlinear behavior beyond 1–2% strain. Traditionally, this behavior was typically represented with the use of viscoelastic material models. ^13,14 However, the use of state variable-based viscoplastic material models, typically used for metals, has gained significant popularity ^15,16 in the last two decades. Even though the deformation response in polymers is due to the nonlinear response of long chain molecules as opposed to the propagation of dislocations for metals, a similar unified inelastic strain variable approach can still be utilized with the state variables representing a different physical meaning. Pioneering work in developing numerical material models for polymers has been carried out by Goldberg over the years. ^8,17

In comparison to metals, polymers in general exhibit a strong pressure-dependent behavior, ^{18 –20} which results in different yielding under uniaxial tension, compression, shear, and biaxial loadings. Furthermore, the assumption of volume constancy that typically holds in metals during plastification does not hold for thermoplastic polymers. This effect is particularly noticed in the tensile range (uniaxial and biaxial tensile stress states) and cannot be neglected. The reorientation of molecule chains and the presence of damages lead to anisotropic behavior of the polymers. The aforementioned characteristics make modeling polymers a challenge much less as a part of a micro-mechanics-based composite material model. A material model that encapsulates these behaviors of polymers has been initially proposed by Goldberg et al. ²¹ and further enhanced by Zheng and Binienda. ²²

Capturing the failure modes accurately is another important aspect which enhances the complexity involved in the numerical modeling of composites. As observed from the work of Mayes and Hansen, ²³ a wide range of failure modes are to be captured for the accurate representation of a composite material model. Typically, many researchers ^4,24,25 have accomplished this using a macro-mechanical approach (where the composite material is modeled as an anisotropic, homogenous material), due to the simplicity and computational efficiency it offers. However, a more realistic and accurate way of doing this is using micro-failure criterion in a micro-mechanical approach (where a representative volume cell (RVC) composed of the fiber and matrix is used to represent the composite).

Typically, damage growth in a polymeric composite manifests itself in the form of strain softening of the material. Matzenmiller et al. ²⁶ developed a continuum damage mechanics (CDM) based model which utilizes the Weibull damage functions to progressively damage the material, thereby introducing strain softening into the constitutive response. This model has been observed to improve the prediction of impact damage on composites significantly. ^{4,27 –30} However, these works use ad hoc softening parameter values in the Weibull functions to predict damage depending on the composite and the structure that is being analyzed, thereby lacking a proper way to characterize the composite material models.

The objective of the present work is to develop and enhance on the three-dimensional material model built by Tabiei and Aminjikarai ³⁰ for predicting the impact response of unidirectional composites (UDC). As explained in the following sections, this is achieved by

implementing strain rate sensitive and pressure-dependent constitutive relations (modified Bodner–Partom viscoplastic state variable model ^22,31 ) for the resin in a micro-mechanical-based RVC which utilizes isostrain boundary conditions and

implementing a CDM-based progressive post-failure model and develop a methodology that assigns physical significance to the choice of damage/softening parameters.

Micro-mechanics of the UDC

The RVC used to develop the micro-mechanical relations is shown in Figure 1. This RVC is the same as the one originally proposed by Pecknold and Rahman ³² and further used in various micro-mechanical models in the literature ^29,30,33 and in the work of Medikonda et al. ³⁴ However, for completeness, a brief description of the RVC is provided here. The fibers are assumed to be of square cross section for computational efficiency; this is especially beneficial since a model implemented in an explicit FE code uses very small time steps to satisfy the stability criterion. The unit cell is divided into three subcells: one fiber subcell, denoted by f , and two matrix subcells, denoted by M_A and M_B , respectively. The three subcells are grouped into two parts: material part A consists of the fiber subcell f and the matrix subcell M_A , and material part B consists of the remaining matrix M_B . The dimensions of the unit cell are 1 × 1 unit square. The dimensions of the fiber and matrix subcells are denoted by W_f and W_m, respectively, as shown in Figure 1 and defined as shown below

Figure 1.

A representative volume cell of unidirectional fiber-reinforced polymer composite.

$W_f=\sqrt{V_f};\hspace{0.17em}{W}_m=1-W_f$

1

where V_f is the fiber volume fraction. As explained in detail in “RVC stress calculations” section below, effective stresses in the RVC are determined from the subcell values in two phases: first, stresses in fiber f and matrix M_A are combined to obtain effective stresses in part A which are then combined with stresses in matrix M_B to obtain the effective RVC stresses.

Viscoplastic constitutive relations for matrix material

As discussed in the Introduction, the strain rate-dependent behavior of a polymer matrix composite is mainly attributed to the viscoplastic nature of the resin component. Hence, strain rate dependency is incorporated in the current model using the viscoplastic relationship developed by Goldberg et al. ²¹ and further enhanced by Zheng and Binenda. ²² It has to be noted that the model proposed by Goldberg et al. ²¹ is a modified form of the Bodner–Partom viscoplastic state variable model, originally developed for metals. For completeness, the modified viscoplastic relations are discussed briefly in this section. Further details about the relations can be found in the study by Goldberg et al. ²¹

The total strain rate is assumed to be the sum of elastic and inelastic strain rates. The elastic strain rate is equal to the ratio of stress rate to Young’s modulus of the material, while the inelastic strain rate is a function of the deviatoric stress, S_ij, the second invariant of the deviatoric stress tensor J₂, and an isotropic state variable Z, representing the resistance to the molecular flow

${\dot{\epsilon }}_{ij}^I=2D_0\text{exp}\left[-\frac{1}{2}{\left(\frac{Z}{{\sigma }_e}\right)}^{2n}\right]\left(\frac{S_{ij}}{2\sqrt{J_2}}+\alpha {\delta }_{ij}\right)$

2

where D₀ and n are the material constants representing the maximum inelastic strain rate and the rate dependence of the deformation response, respectively. The effective stress σ_e is defined as

${\sigma }_e=\sqrt{3J_2}+\sqrt{3}\alpha {\sigma }_{kk}$

3

where α is the state variable controlling the level of hydrostatic stress effects and the two state variables (Z, α) are considered to be proportional to the effective deviatoric inelastic strain rate ${\dot{e}}_e^I$ .

$\dot{Z}=q(Z_1-Z){\dot{e}}_e^I$

4

$\dot{\alpha }=q({\alpha }_1-\alpha ){\dot{e}}_e^I$

5

The effective deviatoric inelastic strain rate, ${\dot{e}}_e^I$ , and the mean inelastic strain rate, ${\dot{\epsilon }}_m^I$ , are given by the following relations

$\begin{array}{l}{\dot{e}}_e^I=\sqrt{\frac{2}{3}{\dot{e}}_{ij}^I{\dot{e}}_{ij}^I}\text{; }\\ {\dot{e}}_e^I={\dot{\epsilon }}_{ij}^I-{\dot{\epsilon }}_m^I\&{\dot{\epsilon }}_m^I=({\dot{\epsilon }}_{11}^I+{\dot{\epsilon }}_{22}^I+{\dot{\epsilon }}_{33}^I)/3\end{array}$

6

It has to be noted that Z₁ and α₁ representing the maximum values of Z and α, respectively, are input as material constants along with their initial values Z₀ and α₀. In addition, q is a material constant representing the “hardening” rate.

One procedure for determining the resin material constants can be found in the original work of Goldberg et al. ²¹ ; alternatively, the material constants can also be determined by fitting the constitutive model to the experimental results using an error minimization scheme such as the least squares method.

It is well known that the elastic modulus significantly increases with increasing strain rate in polymers and subsequently in polymer matrix composites. Zheng and Binenda, ²² drawing inspiration from the work of Yen, ²⁷ have proposed the use of dynamic elastic moduli that account for the strain rate effect

$\begin{array}{l}E=E_0\left(1+C\text{ln}\frac{\dot{\epsilon }}{{\dot{\epsilon }}_0}\right)\text{, where}\\ \dot{\epsilon }=\sqrt{\frac{2}{3}[{({\dot{\epsilon }}_{11}-{\dot{\epsilon }}_m)}^2+{({\dot{\epsilon }}_{22}-{\dot{\epsilon }}_m)}^2+{({\dot{\epsilon }}_{33}-{\dot{\epsilon }}_m)}^2+2{\dot{\epsilon }}_{12}^2+2{\dot{\epsilon }}_{23}^2+2{\dot{\epsilon }}_{31}^2]}\\ {\dot{\epsilon }}_m=({\dot{\epsilon }}_{11}+{\dot{\epsilon }}_{22}+{\dot{\epsilon }}_{33}\text{)/3}\end{array}$

7

Since there is no closed form solution for the first-order differential equations shown in equations (2) through (7), which formulate one differential equation per component, a numerical solution is obtained at each explicit FE time step.

The incremental forms of equations (2) and (4) to (6) are as shown below

$\text{d}{\epsilon }_{ij}^I=\left(2D_0\text{exp}\left[-\frac{1}{2}{\left(\frac{Z}{{\sigma }_e}\right)}^{2n}\right]\left(\frac{S_{ij}}{2\sqrt{J_2}}+\alpha {\delta }_{ij}\right)\right)\text{d}t$

8

$\begin{array}{l}\text{d}{e}_e^I=\sqrt{\frac{2}{3}\text{d}{e}_{ij}^I\text{d}{e}_{ij}^I},\\ \text{where d}{e}_e^I=\text{d}{\epsilon }_{ij}^I-\text{d}{\epsilon }_m^I\text{ and d}{\epsilon }_m^I=(\text{d}{\epsilon }_{11}^I+\text{d}{\epsilon }_{22}^I+\text{d}{\epsilon }_{33}^I)/3\end{array}$

9

$\text{d}Z=q(Z_1-Z)\text{d}{e}_e^I$

10

$\text{d}\alpha =q({\alpha }_1-\alpha )\text{d}{e}_e^I$

11

where the prefix “d” represents the respective incremental values and the inelastic strain increment is multiplied by the time step “dt.” It has to be noted that in an explicit FE code, at time step n + 1, the following values are available and are all used as input parameters for the matrix material stress response calculation: stress of the matrix material from the previous time step n, ${\sigma }_{ij}^{m(n)}$ , strain at the current time step, ${\epsilon }_{ij}^{(n+1)}$ , strain rate at the current time step, ${\dot{\epsilon }}_{ij}^{(n+1)}$ , inelastic strain from the previous time step, ${\epsilon }_{ij}^{I(n)}$ , and lastly the state variables from the previous time step, $Z^{m(n)}$ and ${\alpha }^{m(n)}$ . Since the current material model is being developed to study impact problems, there is a high degree of possibility that the model might encounter high strain increments. Hence, in order to ensure stability, a four-step Runge–Kutta time integration has been performed as shown in Figure 2 for the resin material.

Figure 2.

Illustration of a four-step Runge–Kutta method implemented for the resin material.

Constitutive relations of fibers

The fibers are considered to be linearly elastic materials which are initially transversely isotropic but become orthotropic with damage evolution. Damages to the fibers are assumed to be a result of only the direct stresses and the shear stresses do not contribute here. This causes the damages to be oriented in the material directions of the fibers and act independently. The constitutive relations of the fibers can be written in matrix form as

${\left\{\sigma \right\}}^f=\left[C_f\right]{\left\{\epsilon \right\}}^f$

12

where $\left[C_f\right]$ is the stiffness matrix which can be partitioned into direct and shear stress stiffness matrices as follows

$\left[C_f\right]={\left[S_f\right]}^{-1}=\left[\begin{array}{cc}{\left[S_{fd}\right]}^{-1}& {\left[0\right]}_{3\times 3}\\ {\left[0\right]}_{3\times 3}& {\left[S_{fs}\right]}^{-1}\end{array}\right]$

13

The direct stress compliance matrix, whose inverse is the direct stress stiffness matrix, should be symmetric and the following relationship should be obeyed

$\frac{{\nu }_{ij}}{E_i}=\frac{{\nu }_{ji}}{E_j}\hspace{0.17em},\hspace{0.17em}i,j=1,2,3\text{ and }i\ne j\hspace{0.17em}\left(\text{no\hspace{0.17em}summation}\right)$

14

The direct and shear stress compliance matrices in terms of the properties of the fibers are

$\left[S_{fd}\right]=\left[\begin{array}{ccc}\frac{1}{(1-d_1)E_1}& -\sqrt{\frac{{\nu }_{12}}{(1-d_1)E_1}\frac{{\nu }_{21}}{(1-d_2)E_2}}& -\sqrt{\frac{{\nu }_{12}}{(1-d_1)E_1}\frac{{\nu }_{21}}{(1-d_3)E_2}}\\ & \frac{1}{(1-d_2)E_2}& -\frac{{\nu }_{23}}{\sqrt{(1-d_2)E_2(1-d_3)E_2}}\\ \text{Symm}.& & \frac{1}{(1-d_3)E_2}\end{array}\right]$

15

$\left[S_{fs}\right]=\left[\begin{array}{ccc}\frac{1}{G_{12}}& 0& 0\\ & \frac{1}{G_{23}}& 0\\ \text{Symm}.& & \frac{1}{G_{o12}}\end{array}\right]$

16

where E₁ and E₂ are the longitudinal and transverse moduli of the fibers, respectively; ν_ij, $i,j=1,2,3$ , and $i\ne j$ are its Poisson’s ratios; $G_{o12}$ and $G_{12}$ are its initial and strain rate-dependent in-plane shear moduli, respectively; $G_{23}$ is its transverse shear moduli; and d_i, $i=1,2,3$ , are damage parameters which are discussed in the following section on progressive failure modeling. In-plane shear is considered as a parameter which is strain rate dependent. Hence, the shear modulus $G_{12}$ is given by the following relation

$G_{12}=a_G{s}_s+G_{o12}$

17

$s_s=\frac{1}{t}{\displaystyle \underset{0}{\overset{t}{\int }}\text{log}\left(\frac{\left|{\dot{\epsilon }}_{12}\right|}{{\dot{\epsilon }}_o}\right)\text{dt}}$

18

where a_G is a material constant which expresses the strain rate sensitivity of $G_{12}$ , t is the time elapsed, $G_{o12}$ is the initial in-plane shear modulus of the fibers, ${\dot{\epsilon }}_{12}$ is the in-plane shear strain rate, and ${\dot{\epsilon }}_o$ is a basic strain rate with which the current strain rates are compared and is accepted as the strain rate of the static loading for a given working strain rate range. A time integration of the strain rates is deemed necessary as they are not constant in impact simulations and also because the constitutive response for the fibers is based on their secant stiffness and not the tangential stiffness.

Damage evolution in constituents

Preserving the thermodynamic constraints in damage mechanics, damage is considered an irreversible process and has been implemented as such in the current model. In their work, Matzenmiller et al. ²⁶ proposed an approach which used a Weibull distribution function to describe the damage evolution as a function of strain over the entire load range. Depending on the failure mode considered, a slightly modified version of this model is used here. The evolution function for damage describing fiber breakage at time step n + 1 is expressed as follows

$d_1^{(n+1)}=\text{max}\left\{1-\text{exp}\left[-\frac{1}{m_{1}e}{\left(\frac{{\sigma }_{1t|c}^{\text{fud}}}{{\sigma }_{1t|c}}\right)}^{m_1}\right]\hspace{0.17em},\hspace{0.17em}{d}_1^{(n)}\hspace{0.17em}\right\}$

19

where $t|c$ denotes tension or compression. When ${\epsilon }_{11}>0$ , the parameters for tension are utilized; otherwise the parameters provided for compression are used. When the damage d ₁ reaches 0.01 in tension, the FE is considered to have totally failed and is deleted from the simulation. ${\sigma }_{1t|c}^{\text{fud}}$ is the undamaged stress in the fibers, which is calculated using the same equations given in equations (12) to (15) except that the damages in those equations are considered to be zero, that is, d₁ = d₂ = d₃ = 0. Also, ${\sigma }_{1t|c}$ is not simply the strength of the fibers but is a reduced value given by the following relation

${\sigma }_{1t|c}=\frac{X_{t|c}}{b_{t|c}}$

20

where $X_{t|c}$ is the tensile/compressive strength of the pure fibers and $b_{t|c}$ is a reduction factor. This is done as the fibers generally display reduced strengths in UDCs as compared to stand-alone tests done on themselves, especially for uniaxial longitudinal tension.

In the transverse directions, the damage evolution functions (equations (21) and (22)) are considered to independent of the strain rate. The properties of the fibers in both transverse directions are the same; therefore, the evolution functions as well as their parameters are the same and only the history of the loading is different.

$d_2^{(n+1)}=\text{max}\left\{1-\text{exp}\left[-\frac{1}{m_{1}e}{\left(\frac{{\sigma }_{2t|c}^{\text{fud}}}{{\sigma }_{2t|c}}\right)}^{m_1}\right]\hspace{0.17em},\hspace{0.17em}{d}_2^{(n)}\hspace{0.17em}\right\}$

21

$d_3^{(n+1)}=\text{max}\left\{1-\text{exp}\left[-\frac{1}{m_{1}e}{\left(\frac{{\sigma }_{3t|c}^{\text{fud}}}{{\sigma }_{2t|c}}\right)}^{m_1}\right]\hspace{0.17em},\hspace{0.17em}{d}_3^{(n)}\hspace{0.17em}\right\}$

22

The damages in transverse directions of the fibers are constrained to not exceed 0.10.

Damages are imposed on the matrix material, but they affect only the shear stresses of the resin, which is considered to be the main contributor to the shear stresses of the RVC. A single Weibull distribution function is accepted again as an evolution function of the damages but it involves the ultimate strain for the damage development rather than the ultimate stress. The damage evolution function for in-plane shear is as follows

$d_4^{(n+1)}=\text{max}\left\{1-\text{exp}\left[-{\left(\frac{\left|{\epsilon }_{12}\right|}{{\epsilon }_{4m}}\right)}^{m_s}\right]\hspace{0.17em},\hspace{0.17em}{d}_4^{(n)}\hspace{0.17em}\right\}$

23

The damage evolution functions for the matrix material in the out of plane directions for a time step n + 1 are calculated as follows

$d_5^{(n+1)}=\text{max}\left\{1-\text{exp}\left[-{\left(\frac{\left|{\epsilon }_{23}\right|}{{\epsilon }_{5m}}\right)}^{m_s}\right]\hspace{0.17em},\hspace{0.17em}{d}_5^{(n)}\hspace{0.17em}\right\}$

24

$d_6^{(n+1)}=\text{max}\left\{1-\text{exp}\left[-{\left(\frac{\left|{\epsilon }_{31}\right|}{{\epsilon }_{4m}}\right)}^{m_s}\right]\hspace{0.17em},\hspace{0.17em}{d}_6^{(n)}\hspace{0.17em}\right\}$

25

The damages of the matrix material are constrained to not exceed 0.20. These damages are applied on the matrix material shear stresses when the stress response of the subcells is calculated. It has to be noted that damages for the matrix material are being applied on the stresses unlike what was done for the fibers, where the elastic moduli were reduced. The reason for this is considered to be two-fold; one is because the matrix material model is assumed to be isotropic while the damages in the material are not and secondly, the viscoplastic material model calculates the actual stresses of the material and not the effective stresses of the damaged material. ^30,35 Lastly, the primary reason for constraining the damages is to account for the numerical instabilities that arise when stress in an element goes to zero.

Delamination

Delamination is a failure mode often considered to be a consequence of the quadratic interaction between the out-of-plane stresses of a lamina and is assumed to be mainly a lamina failure. The loading criterion for this failure mode has the following form

$S^2\left\{{\left(\frac{E_3〈{\epsilon }_{33}〉}{S_{3t}}\right)}^2+{\left(\frac{G_{23}{\gamma }_{23}}{S_{230}+S_{SR}}\right)}^2+{\left(\frac{G_{31}{\gamma }_{31}}{S_{310}+S_{SR}}\right)}^2\right\}-r^2=0$

26

Here, E₃ is the normal tensile modulus of the lamina, $G_{23}$ and $G_{31}$ are the transverse shear moduli of the lamina, $S_{3t}$ is the through-the-thickness tensile strength of the lamina, $S_{230}$ and $S_{310}$ are the transverse shear strengths of the lamina for tensile ${\epsilon }_{33}$ , r is the damage threshold, and S is a scale factor introduced to provide better correlation of delamination area with experiments which can be determined by fitting analytical prediction to experimental data for the delamination area. It has to be noted that the compressive normal traction does not normally contribute to the delamination onset and this explains the use of the Macaulay brackets, $〈•〉=\left(|•|+•\right)/2$ , in equation (26). However, under compressive through-the-thickness strain, that is, ${\epsilon }_{33}<0$ , the damaged/delaminated surface is considered to be “closed” and the out-of-plane shear strengths are assumed to increase due to the presence of friction. This is accounted for in the current model using the Coulomb–Mohr theory

$S_{SR}=E_3\text{tan}\phi 〈-{\epsilon }_z〉$

27

where φ is the Coulomb’s friction angle. The normal tensile and the out-of-plane shear moduli of the lamina are computed at the first time step and stored as a material property.

$d_{\text{lam}}^{(n+1)}$ is the damage variable associated with this failure mode and its evolution is given by the relation

$d_{\text{lam}}^{(n+1)}=\text{max}\left\{1-\text{exp}\left[\frac{1}{m_d}\left(1-r^{m_d}\right)\right],d_{lam}^n\right\}$

28

where r is the damage threshold as given in equation (26) and m_d is the damage exponent for delamination. Delamination damage is constrained at 0.10 to avoid numerical difficulties, and when this maximum value is reached in an element, it is considered to be fully delaminated.

When delamination failure given by equation (26) occurs in an element, depending on the opening or closing of the damage surface, the damages variables ( $d_z^{(n+1)}$ , $d_{yz}^{(n+1)}$ , and $d_{zx}^{(n+1)}$ ) specified in equations (29) and (30) are applied on the effective RVC stresses.

$\text{For tensile\hspace{0.17em}mode},\hspace{0.17em}{\epsilon }_{33}>0:\hspace{0.17em}{d}_z^{(n+1)},d_{yz}^{(n+1)},d_{zx}^{(n+1)}=d_{\text{lam}}^{(n+1)}$

29

$\text{For compressive mode},{\epsilon }_{33}<0:\hspace{0.17em}{d}_{yz}^{(n+1)},d_{zx}^{(n+1)}=d_{\text{lam}}^{(n+1)}$

30

For tensile mode, all the through-the-thickness stress components ${\sigma }_{33}$ , ${\sigma }_{23}$ , and ${\sigma }_{31}$ of RVC are reduced, while for compressive mode, since the damage surface is considered to be closed, only ${\sigma }_{23}$ and ${\sigma }_{31}$ are reduced.

RVC stress calculations

The effective stresses in the RVC are determined from the subcell values in two phases: first, stresses in the fiber f and matrix M_A are used to determine the effective stresses of part A ; then these stresses along with the stresses in matrix M_B are used to determine the effective stresses in the RVC.

In the homogenization process, isostrain boundary conditions are assumed for all the three subcells of the RVC, which implies that the rule of mixture is applied for the stress calculations. In a physical sense, this means that the fiber and matrix materials act in parallel in all directions under loading, which isn’t quite realistic. ³⁴ However, this assumption is made in order to simplify the micro-mechanical relations.

The direct stresses of part A are calculated from the direct stresses of the fiber subcell f and the matrix subcell M_A using the following relations

${\sigma }_{11}^A=W_f{\sigma }_{11}^f+(1-W_f){\sigma }_{11}^{m_A}$

31

${\sigma }_{22}^A=W_f{\sigma }_{22}^f+(1-W_f){\sigma }_{22}^{m_A}$

32

${\sigma }_{33}^A=W_f{\sigma }_{33}^f+(1-W_f){\sigma }_{33}^{m_A}$

33

The assumption of isostrain boundary conditions in all directions mainly poses problems in accurately estimating the shear behavior. Now it is well known that the behavior of UDCs under shear is dominated by the behavior of the matrix material. Hence, in order to accurately represent this, ad hoc volume fraction coefficients are implemented for shear and a rule of mixture involving them is applied. Then, the shear stress in part A is determined by applying the damages to the matrix material introduced in the previous section.

${\sigma }_{12}^A=V_{s4}{\sigma }_{12}^f+(1-V_{s4})(1-d_4){\sigma }_{12}^{m_A}$

34

${\sigma }_{23}^A=V_{s5}{\sigma }_{23}^f+(1-V_{s5})(1-d_5){\sigma }_{23}^{m_A}$

35

${\sigma }_{31}^A=V_{s4}{\sigma }_{31}^f+(1-V_{s4})(1-d_6){\sigma }_{31}^{m_A}$

36

The different shear volume fraction coefficients, V_{s 4} and V_{s 5}, for the in-plane and transverse shear have values quite lower than the volume fraction of the fibers. Since the matrix material is modeled as viscoplastic and the fibers are modeled as elastic, after the saturation of the plasticity in the matrix material (shear directions), the fibers tend to absorb the load and assume the role of strain hardening.

Finally, the effective stresses in the RVC are obtained by applying the rule of mixtures again which yields the following relations including the softening of through-the-thickness components due to delamination failure as follows

${\left\{\begin{array}{c}{\sigma }_{11}\\ {\sigma }_{22}\\ {\sigma }_{33}\\ {\sigma }_{12}\\ {\sigma }_{23}\\ {\sigma }_{31}\end{array}\right\}}_{\text{RVC}}=\left[\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ & 1& 0& 0& 0& 0\\ & & d_z^{(n+1)}& 0& 0& 0\\ & & & 1& 0& 0\\ & \text{symm}& & & d_{yz}^{(n+1)}& 0\\ & & & & & d_{zx}^{(n+1)}\end{array}\right]\left(W_f{\left\{\begin{array}{c}{\sigma }_{11}\\ {\sigma }_{22}\\ {\sigma }_{33}\\ {\sigma }_{12}\\ {\sigma }_{23}\\ {\sigma }_{31}\end{array}\right\}}_A+W_m{\left\{\begin{array}{c}{\sigma }_{11}\\ {\sigma }_{22}\\ {\sigma }_{33}\\ {\sigma }_{12}\\ {\sigma }_{23}\\ {\sigma }_{31}\end{array}\right\}}_{M_B}\right)$

37

The total stresses of the fiber, the internal state variables along with the total stresses of the matrix material, and the progressive damage variables are stored as history variables at each time step of the explicit time integration process for the next time step calculations.

Damage parameter estimation

The damage parameter that has been used in the Weibull functions above (equations (19), (21), (22), etc.) has been observed to be very problem dependent and difficult to characterize. A small value of the damage exponent (e.g. m₁ in equation (19)) makes the material behave in a very ductile manner and the behavior becomes increasingly brittle as this value increases. Hence, it is difficult to obtain the softening response of most quasi-brittle materials. The softening response heavily depends on the setup and test machines, which can lead to very different results. The choice of damage parameters for each mode has been debated by Williams and Vaziri ⁴ and Tabiei and Aminjikarai. ³⁰ A procedure for the calculation of softening parameter used in the damage evolution function has been presented in the current study.

In order to avoid undesirable localization, Bazant and Oh ³⁶ proposed the crack band theory. According to this theory, regardless of the choice of mesh, the overall energy dissipation due to damage must remain constant. This concept was further used by Pinho et al. ³⁷ and subsequently Cousigné et al., ³⁸ who have discussed a smeared formulation involving damage variables, where a length parameter has been introduced into the damage constitutive law. Their work primarily relies on the calculation of failure strain based on the fracture toughness, which is subsequently used to calculate the damage in the composite. This calculation of a failure strain is of little consequence to the damage function discussed in the current work. Hence, a numerical algorithm which iteratively solves for the softening parameter has been presented.

The total energy absorbed by one element after failure can be expressed as

$U=V\times U_{\text{vol}}$

38

where $U_{\text{vol}}$ is the total energy absorbed per unit volume; in simpler terms, it is area captured under the stress–strain curve before and after the onset of failure. The volume of the element is denoted by V. Based on Figure 3, the total energy can also be represented in terms of the fracture plane area A and the energy per unit area or the fracture toughness Γ.

Figure 3.

Fracture plane for a load. ³⁸

$U=A\times \Gamma$

39

Equations (38) and (39) can be equated and the energy dissipated per unit volume absorbed by an element is given by the following relation

$U_{\text{vol}}=\Gamma \times \frac{A}{V}=\frac{\Gamma }{l}$

40

where l is the characteristic element length. For square elements, it is expressed in terms of the area at the integration point by the following relation ³⁸

$l=1.12\times \sqrt{A_{ip}}$

41

Since the fracture toughness Γ is a constant known value for a particular material, we can solve for the softening parameter m_f numerically as $U_{\text{vol}}$ is known for a fixed mesh size. The procedure of accomplishing this is presented in the flowchart shown in Figure 4 and has been implemented as a user subroutine in the material which is called at the start of any simulation.

Figure 4.

Flowchart demonstrating the calculation of the damage parameters.

It has to be noted that once the softening parameter m_f has been calculated for a certain $U_{\text{vol}}$ and mesh size, it becomes the default value. For a new mesh/element size, $U_{\text{vol}}$ is scaled appropriately depending on the element’s characteristic length and the softening parameter is recalculated using the subroutine discussed above

$U_{\text{vol}}^{\text{new}}=U_{\text{vol}}^{\text{def}}\times \frac{l^{\text{def}}}{l^{\text{new}}}$

42

The validity of the subroutine has been first verified by implementing it in MATLAB^® for a T300/913 composite laminate and different element sizes ranging from 0.25 mm to 4 mm. The necessary initial inputs are $E=132$ GPa, ${\sigma }^{\text{onset}}=2005$ MPa, and $\Gamma =91.6$ MPa mm. The default element size for this problem has been considered to 1 mm; hence, Γ will be equal to $U_{\text{vol}}$ for this element size and $U_{\text{vol}}$ will be appropriately scaled according to equation (42) for the remaining element sizes.

Figure 5 illustrates that the algorithm calculates different damage parameters for different element sizes. It has to be noted that the damage exponent increases with increasing element length. The corresponding stress versus strain response has been illustrated in Figure 6.

Figure 5.

Calculated damage versus strain.

Lastly, validity of the whole procedure can be verified by back calculating the fracture toughness for each element size. This can be done by measuring the area under each curve in Figure 6 and multiplying it with the corresponding element length. Figure 7 demonstrates that a constant value of fracture toughness that has been input into the subroutine is reproduced for every element length considered, that is, $\Gamma =91.6$ MPa mm.

Figure 6.

Stress versus strain response for the different damages computed.

Figure 7.

Corresponding fracture toughness back calculated.

The damage exponent is calculated individually for each element inside the mesh and this can be quite taxing for large-scale problems. As it can be observed in the flowchart presented in Figure 4, the algorithm draws inspiration from the bisection method and thus determines the damage parameter rather quickly. The speed at which the damage parameter is determined is further illustrated in Figure 8, for an input/target fracture toughness of 21.7 MPa mm and an initial guess of 2 for the damage exponent. The solution converges in just nine iterations. Of course, it has to be noted that the speed of the solution is slightly influenced by the initial guess provided to the subroutine.

Figure 8.

Fracture toughness versus damage exponent for an element with characteristic length = 1.

Results and discussion

Impact on a 10-ply E-glass/epoxy laminate

In the current section, impact loading on E-glass/epoxy laminate plates has been simulated using the composite material model discussed so far. The experiments have been carried out in a drop-weight setup by Collombet et al. ¹ The plates are circular with a diameter of 200 mm and a thickness of 1.8 mm. They are put together by stacking 10 unidirectional plies in one of the three configurations $\left[0_n/{90}_m/0_n\right]$ , where $2n+m=10$ with $n=2,3,4$ . The cylindrical projectile impacting the laminate plates with an incident energy of 27 J has a hemispherical end of mass 2.3 kg, density 7.8 kg m⁻³, length 600 mm, diameter 25 mm, and the following properties: $E=210$ GPa and $\nu =0.3$ . The plates are clamped on their periphery over a 20-mm wide ring and a uniform pressure is applied on the ring by means of springs. Additional details of the setup can be found by Collombet et al. ¹

The FE models used in the simulations are shown in Figure 9. As seen is these figures, only a quarter of the experimental setup is modeled due to symmetry, and this helps reduce the computational effort. In addition, simulations are carried out for a coarse mesh (Figure 9(a)) and a fine mesh (Figure 9(b)). As seen in Figure 9(c), each lamina is modeled using a single layer of solid elements. The support on which the laminate rests is completely constrained at the bottom and is modeled as a rigid material. The ring is constrained from moving in the vertical Z-direction to account for the pressure applied on it, and this replicates the uniform pressure applied on the springs. ¹ Lastly, the impactor is modeled as a rigid material as well and only the contact surface of the hemispherical end is modeled using shell elements. Eroding surface-to-surface contact is defined between the impactor and the plates, while automatic surface-to-surface contact is defined between the support, the ring, and the plates. The constituent properties and parameters used in the simulations are given in Table 1. It has to noted that some of the required material parameters which are not available in the study by Collombet et al. ¹ have either been obtained from the literature ^17,22,30,35 or valid estimates have been used. For example, the parameters required for the viscoplastic constitutive relations of the epoxy resin are not available directly in the literature and neither are their uniaxial and shear responses at different strain rates from which they can be determined. Hence, the parameters of the epoxy used are based on the values provided for the E-862 resin ²² which shares similar elastic properties with the epoxy resin.

Figure 9.

FE model used for FE simulation of a 27 J impact on $\left[0_n/{90}_m/0_n\right]$ , n = 2, 3, 4, E-glass/epoxy laminates—(a) full model with coarse mesh, (b) full model with fine mesh, and (c) close-up view of the laminate.

Table 1.

Properties and parameters of E-glass/epoxy materials used in verification examples.

Vf	${\dot{\epsilon }}_o$ s−1	E1 GPa	${\nu }_{12}$	Xt MPa	Xc MPa	E2 GPa	${\nu }_{23}$
0.43	1 × 10–3	74	0.25	3500	1600	4.93	0.25
${\sigma }_{2t}$ MPa	${\sigma }_{2c}$ MPa	$G_{o12}$ GPa	$G_{23}$ GPa	$a_{s4}$	${\epsilon }_{4m}$	${\epsilon }_{5m}$	$V_{s4}$
80	255	30	1.97	0.011	0.15	0.15	0.014
$V_{s5}$	aG GPa	bt	bc	Em GPa	νm	$D_0(s^{-1})$	n
0.020	0.900	1.00	1.00	3.00	0.4	1 × 10 6	0.80
Z0 MPa	Z1 MPa	q	α 0	α 1	C	$S_{3t}$ MPa	$S_{230}$ MPa
420.0	820.0	120.0	0.15	0.06	0.0	100	86
$S_{310}$ MPa	φ°
64	20

Rate sensitivity of the composite material is predominantly accounted by the resin in the proposed material model. In the parametric studies carried out by Zheng, ³⁹ it has been observed that Z₁ and α₁ have a significant effect on the strain rate. Increasing value of Z₁ and decreasing value of α₁(in tension) tend to increase the maximum/yield stress. Parameters Z₀ and α₀ on the other hand have been observed to control the starting point of the nonlinearity.

Ad hoc volume fractions control the effect that the resin behavior has on the overall composite which in turn effects its rate-dependent behavior in shear. However, this is subjective and has to be determined experimentally depending on the stress–strain responses of the composite and its constituents. In general, small changes in the ad hoc volume fraction parameters will not have a great influence on the answer.

It is important to note that the criterion for element erosion (fiber breakage), discussed in “Damage evolution in constituents” section, has been deactivated for the cases studied. Mainly because of the experimental reference paper by Collombet et al., ¹ there was no reference to any fiber breakage.

Figures 10 –12 show the contact force versus time history for all three-laminates simulated along with the experimental results. It must be noted that the contact force obtained from the simulations has been multiplied by a factor of four since only a quarter of the model has been subjected to impact. Based on the observed results, it can be said that the effect of the mesh sizes considered has been minimum.

Figure 10.

Contact force–time history predicted by the current model for a 27 J impact on a $\left[0_2/{90}_6/0_2\right]$ E-glass/epoxy laminate.

Figure 11.

Contact force–time history predicted by the current model for a 27 J impact on a $\left[0_3/{90}_4/0_3\right]$ E-glass/epoxy laminate.

Figure 12.

Contact force–time history predicted by the current model for a 27 J impact on a $\left[0_4/{90}_2/0_4\right]$ E-glass/epoxy laminate.

The peak contact forces obtained for $\left[0_2/{90}_6/0_2\right]$ and $\left[0_4/{90}_2/0_4\right]$ laminates are found to be the closest to their respective experimental results, whereas the peak force predicted for the $\left[0_3/{90}_4/0_3\right]$ laminate is found to be slightly higher than the experimental value. However, it has to be noted that this error in peak force (which is approximately 5%) is still smaller than the error reported in Tabiei and Aminjikarai ³⁰ (approximately 7%). It must be noted that all six numerical simulations predict a slightly bigger mismatch during unloading possibly due to the assumed boundary conditions, contact algorithm, or the approximate material properties. Overall, the predicted contact force history is found to agree reasonably well with the experimental result in each case.

Effect of hydrostatic stresses on impact problems

As discussed in the “Introduction” section, hydrostatic stresses have been observed to have a significant contribution to the response of a polymer. In the current model, hydrostatic stress effects are accounted through the inclusion of a state variable “α”, and this can be turned off by setting ${\alpha }_0={\alpha }_1=0.0$ . A single element test carried out using the material properties given in Table 1 illustrates the importance of including hydrostatic stress effects (Figure 13), as observed there is a difference of nearly 33% in the maximum stresses observed.

Figure 13.

Comparing the stress response of the epoxy resin with and without the inclusion of the hydrostatic stresses (“α”) for a strain rate of 1.0/s.

In this section, the importance of having hydrostatic stresses in analyzing the impact response has been further investigated. In order to do this, the fine mesh models discussed in the previous section have been simulated with and without the inclusion of “α” and the results have been presented in Figures 14 to 16.

Figure 14.

Contact force–time history predicted by the current model for a 27 J impact on a $\left[0_2/{90}_6/0_2\right]$ E-glass/epoxy laminate with and without “α.”

Figure 15.

Contact force–time history predicted by the current model for a 27 J impact on a $\left[0_3/{90}_4/0_3\right]$ E-glass/epoxy laminate with and without “α.”

Figure 16.

Contact force–time history predicted by the current model for a 27 J impact on a $\left[0_4/{90}_2/0_4\right]$ E-glass/epoxy laminate with and without “α.”

Based on these results, it can be said that the peak contact force predicted by the simulations with the inclusion of hydrostatic stress effects is closer to the experimental values. The exact amount of error for each case is reported in Table 2 and the maximum error has been reported in the simulations of $\left[0_3/{90}_4/0_3\right]$ E-glass/epoxy laminate without “α.”

Table 2.

Comparison of error in the peak contact force reported for each laminate configuration.

Laminate configuration	Error with experiments
	With α (%)	Without α (%)
$\left[0_2/{90}_6/0_2\right]$	2.672	3.178
$\left[0_3/{90}_4/0_3\right]$	5.268	11.213
$\left[0_4/{90}_2/0_4\right]$	0.168	1.992

Impact event on a 22-ply E-glass/epoxy laminate

In order to better understand the hydrostatic stress effects, impact events have also been carried out on a much thicker E-glass/epoxy plate. The FE model used in the simulations is shown in Figure 17. Once again, in order to reduce the computational effort, only a quarter of the laminate is modeled due to symmetry and each lamina is modeled using a single layer of solid elements. The model consists of 22 unidirectional plies with a stacking sequence of $\left[{45}_{22}\right]$ . The laminate model is square in shape with a side length of 198 mm and a thickness of 13.2 mm, whereas the impactor is modeled as a rigid part which has a length of 40 mm and a radius of 6.35 mm. In addition, the laminate is constrained on the top and bottom with rigid supports with a radius of 12.7 mm cut out at the corner where the impact event occurs.

Figure 17.

Full model and close-up view of the FE model used for simulation of impact on a thick plate.

Unlike the low-velocity boundary condition used in the earlier case, the impactor in the current case is being displaced by 10 mm (in the downward direction) to observe the response of the plate at a higher state of deformation. Simulations have been carried out with and without the inclusion of “α” to compare the effect of hydrostatic stress in predicting the contact force and the results have been presented in Figure 18. It has to be noted that due to the lack of experimental data, only the numerical results have been compared in this figure.

Figure 18.

Contact force–time history predicted by the thick plate model on a $\left[{45}_{22}\right]$ E-glass/epoxy laminate with and without “α.”

Results presented in Figure 18 show that at higher states of deformation, the contribution resulting from the inclusion of “α” shows a significant effect on the contact force. For the current case, the simulation with “α” predicts a peak force of approximately 9981 N, whereas the case without “α” predicts a value of approximately 8518 N.

Conclusions

A micro-mechanics based material model is implemented into an explicit dynamic FE code LS-DYNA and is used for simulating the behavior of composite structures under impact loading. The model assumes an RVC and the effective stresses of the unit cell are calculated based on the isostrain boundary conditions between the fiber and resin subcells. The nonlinear behavior of the resin is accounted with the help of a viscoplastic state variable model proposed by Goldberg et al. ²¹ and Zheng and Binenda. ²² This model represents both the strain rate dependency and the pressure dependency usually encountered in polymeric materials.

Progressive post-failure behavior is modeled using a CDM-based damage model. Different Weibull damage functions are used for various failure modes such as fiber failure, matrix cracking, and delamination, and the response of the composite is softened depending on the failure mode. The problem of determining the damage/softening exponents in the damage functions has been addressed using a fast iterative algorithm from within the constitutive model itself. The procedure discussed in this algorithm assigns physical significance to these exponents and even accounts for the strain localization typically observed in FE models.

All of the above discussed characteristics of the micro-mechanical model make it unique and very suitable for impact simulations, validation of which is presented in the current study for impact on E-glass/epoxy laminates with different configurations. Contact force–time histories predicted by the model are found to agree very well with experimental results. In addition, the effect of hydrostatic stresses on the peak impact force is studied in detail, and the inclusion of this characteristic of the polymer has been observed to aid in predicting more accurate results, thereby establishing its significance.