An approach for nonlinear modeling of polymer matrix composites

Abstract

In the present work, a mechanistic modeling approach is pursued for material characterization and for modeling inelastic deformation of polymer matrix composite components. The model attempts to capture the dominant micromechanical deformation mechanisms in laminated composites caused by matrix inelasticity at elevated temperatures. Given material characteristics of the constituent materials, the model can be used in predicting stress, time and temperature-dependent response of a composite under a broad range of thermal and mechanical load conditions. This article describes the modeling approach and examples of its use in a finite element analysis framework. Examples include analyses of simple test specimen coupons, stress concentration at holes and a structural element configuration of a polymer matrix composite. In each case, the model predictions are compared with the experimental measurements.

Keywords

Polymer matrix composite viscoplasticity micromechanics deformation modeling

Introduction

Polymeric matrix composites (PMCs) are enabling materials in modern high-performance turbine engines. In several component applications, temperature and environment conditions are such that the materials are prone to localized nonlinear deformations and damage, such as at stress concentrations. The nonlinear behavior of these composites under these conditions is primarily due to the rate sensitive behavior of the polymeric matrix used. In addition, damage in the form of matrix microcracks, fiber fracture and fiber–matrix debonding can also cause nonlinear deformation in PMCs. But most of the design and life prediction methods available to designers are based on the assumption of linear elastic material behavior. Use of these available methods gives insufficiently accurate predictions of component failure and lifetimes. This can in turn result in over-design, of weight critical parts.

A number of models that describe the nonlinear behavior of PMCs exist in the literature. For example, Ha and Springer¹ have proposed a phenomenological-based viscoplastic constitutive model for unidirectional plies. Similarly a one-parameter flow rule has been used by Sun and Chen² and by Sun and Yoon³ to model orthotropic plasticity in polymeric matrix-based unidirectional plies and laminates. Several researchers^4–8 have used a similar approach to model plasticity and creep in PMCs. In all of the above approaches, the macroscopic behavior of unidirectional plies is modeled and lamination theory is used to extend the model to laminates. Determination of the model parameters typically requires the characterization of the inelastic behavior of several angle-ply laminates. A limitation of such models is that they cannot be used to predict the nonlinear structural response of laminates that contain one or more plies that are different from the ones used to characterize the model.

The unit-cell approach has also been used to model the nonlinear behavior of PMCs. In this approach, the fiber is assumed to be elastic and all the nonlinearity is attributed to the matrix. Sun and Chen⁹ have modeled the matrix as an elastic–plastic material while Tsai and Chen¹⁰ have used a rate sensitive model to describe matrix behavior.

In this article, we propose a micromechanics-based continuum approach to model nonlinear behavior of PMCs. This method is analogous to that used to model inelastic behavior in metal matrix composites.¹¹ While the method includes explicit consideration of inelasticity and damage in composites, in the present work we assume that there is no damage, such as fiber–matrix debonding and matrix cracking, in the composite. We further assume the fiber to be linear elastic and the matrix to exhibit inelastic behavior. Depending on the choice of the nonlinear models, both plasticity and/or creep of the matrix can be modeled. Once the in situ matrix is characterized the model can be used to predict the nonlinear response of the composite.

In the following work, all the observed nonlinearity in the composite behavior is attributed to the matrix. This assumption is consistent with that of the other researchers whose works have been reported here. Damage in the form of matrix microcracks, fiber fracture and fiber–matrix debonding are not considered in the present work. The model described here can nevertheless be easily extended to include consideration of damage by following the approach described in reference.¹¹

Theory

With 1, 2 and 3 as the orthogonal principal material directions of a composite, the second invariant of the deviatoric stress tensor can be defined as

\begin{matrix} J 2 A = \sum_{i = 1}^{3} \frac{1}{3} (N ii σ ii + R ii) 2 + (N 12 σ 12 + R 12) 2 \\ + (N 23 σ 23 + R 23) 2 + (N 31 σ 31 + R 31) 2 \\ - \frac{1}{3} [(N 11 σ 11 + R 11) (N 22 σ 22 + R 22) \\ + (N 22 σ 22 + R 22) (N 33 σ 33 + R 33) \\ + (N 33 σ 33 + R 33) (N 11 σ 11 + R 11)] \end{matrix}

(1)

where

N ij

are dimensionless constants with specific physical meanings. Their values are chosen such that the product

N ii σ ii

represents the component

(σ ij m)

of average stress in the matrix material imposed by external sources. For a given composite system, the values of

N ij

can be established either by employing micromechanics considerations, or by conducting a limited number of laboratory tests. The parameters

R ij

are the components of the average internal (residual) stress in the matrix. These are induced by the consolidation process used in fabricating the composite.

σ ij

are components of the average stress imposed on the composite by external sources, such as mechanical or thermal straining.

An effective stress is defined as

σ e = \sqrt{3 J 2 A}

(2)

and the deviatoric stress components are

\begin{matrix} S ii A = N ii σ ii + R ii - \frac{1}{3} \sum_{i = 1}^{3} (N ii σ ii + R ii) and \\ S ij A = 2 (N ij σ ij + R ij) for i \neq j . \end{matrix}

(3)

Then, the inelastic strain rate in the composite can be expressed in the form of the Prandtl–Reuss relations:

\overset{\cdot}{ɛ} ij I = λ ij S ij A

(4)

where

λ ij

are scalar parameters. To determine

λ ij

, it is assumed that the fiber material is linearly elastic and that the inelastic strain rate in the matrix material (

\overset{\cdot}{ɛ} ij m, I

) is given by

\overset{\cdot}{ɛ} ij m, I = λ S ij m

(5)

Here $S ij m$ is the matrix deviatoric stress tensor. For isotropic solids, $λ$ is a scalar parameter, which depending on the constitutive theory, may be a function of the second invariant $(J 2)$ of the deviatoric stress, temperature and other variables.

With these assumptions, for the case of plane stress (i, j = 1, 2)

λ ij = λ N ij (1 - V f)

(6)

where

V f

is the fiber volume fraction in the composite and

λ

is associated with the viscoplastic theory used in modeling the viscoplastic behavior of the matrix material.

The total strain rate $(\overset{\cdot}{ɛ} ij)$ due to an external stress increment is composed of elastic strain rate $(\overset{\cdot}{ɛ} ij e)$ and inelastic strain rate $(\overset{\cdot}{ɛ} ij I)$ caused by matrix material inelasticity. The elastic strain rate, $\overset{\cdot}{ɛ} ij e$ can be found by using the generalized Hooke’s law for orthotropic solids, and time rates of the elastic constants of the composite assumed to be a continuum.

The parameters $λ ij$ corresponding to each inelastic strain rate component are found by using micromechanics considerations. For example, consider the case when all stress components except $σ 11$ and $R 11$ are zero. Then, according to equations (1) and (2), the composite yields when

N 11 σ 11 + R 11 = σ 0

(7)

The left-hand side of the above equation is the total average stress in the matrix material in direction 1, and the right-hand side is the matrix material’s yield strength in simple tension

(σ 0)

. In direction 1, the load on the composite is shared between the fiber and matrix (Figure 1), and, under elastic loading, the iso-strain condition yields the strain in direction 1 as

ɛ 11 = σ 11 / E 11 = σ 11 m / E m

(8a)

It follows that

N 11 = σ 11 m / σ 11 = E m / E 11

(8b)

Here

E m

is Young’s modulus of the matrix material and

E 11

is the elastic modulus of the composite in the fiber direction given by

E 11 = V f E f + (1 - V f) E m

(9)

In the above equation,

E f

represents Young’s modulus of the fiber. Thus, an estimate of

N 11

is given by equation (8). Alternatively, one can establish

N 11

by using tensile test data on the composite in the fiber direction.

Figure 1.

Two-dimensional stresses in unidirectional two-phase composite.

By assuming, under purely transverse and under purely shear loading, that the average matrix stress is the same as the applied stress (Figure 1), we can get the following expressions for the other dimensionless parameters in equation (1):

N 22 = N 33 = N 12 = N 13 = N 23 ≅ 1.0

(10)

The residual stress,

R ij

, can be found by detailed micromechanics analysis, or estimated by using simpler solutions available in the literature. For example, using an elastic cylinder within a cylinder model of a unit composite cell,¹² estimates for the residual stress components in the matrix are

R 11 = \frac{E m E f V f}{E 11} (\bar{α} f - \bar{α} m) (T - T p),

(11a)

R 22 = \frac{E 11 (1 - V f)}{V f 2} ξ R 11 and

(11b)

R 12 = 0.0

(11c)

where

ξ = \frac{g + ν m E f - h ν f E m}{{\begin{matrix} hE f [g (1 + ν m) + E m (ν m - ν f)] - [h ν f E m - ν m E f - g] \\ [(2 ν m - 1) E f - 2 h ν f E m - g] \end{matrix}}}

g = E f - hE m and h = (V f - 1) / V f

In the above expressions, T is the temperature and

T p

is the consolidation temperature of the composite where it is considered to be free of internal stress. The parameters

\bar{α} f

and

\bar{α} m

are temperature-averaged values of the coefficient of thermal expansion of the fiber and the matrix materials, respectively, and

ν f

and

ν m

are Poisson’s ratio of the fiber and the matrix, respectively.

The parameter $λ$ in equations (5) and (6) depends on the choice of a particular theory for representing the constitutive behavior of the matrix material. Here we use the following power-law equation to describe matrix inelastic behavior:

λ = A (σ e m) n - 1 (ɛ e m, I) p,

(12)

where A(T), n(T) and p(T) are temperature (T)-dependent material specific parameters. Thus according to this model the inelastic strain rate in the matrix depends both on the effective stress,

σ e m

, and on the effective accumulated inelastic strain,

ɛ e m, I

, in the matrix.

The nonlinear problem represented by the above equations requires numerical solution. In the present work, the solutions were obtained using the finite element method.

Model verification

The experimental data selected for the present study are on two PMCs, the AS4/3501-5 thermoset system and the AS4/PEEK thermoplastic system, both of which exhibit nonlinear behavior at room temperature (RT). The test data for these material systems at room temperature were obtained from references ² and³, respectively.

The first step in applying the model is to determine the material parameters relevant to the model. The elastic properties of the constituents and the fiber volume fraction are usually available in the literature. Typical values of the elastic properties of the AS4 graphite fiber and that of the two matrix materials that were analyzed in this study and are listed in Table 1. Since the fiber volume fraction is dependent on the specification of each composite and since it was not reported in references² and ³ it was estimated from the 0° unidirectional modulus using the rule-of-mixtures equation (9).

Table 1.

Material properties of PMCs at room temperature.

AS4/3501–5	AS4/PEEK
E_f (GPa)	234	234
v_f	0.3	0.3
E_m (GPa)	3.9	4.9
v_m	0.3	0.35
V_f	0.58	0.53
A (Mpaⁿs⁻¹)	7.02 × 10⁻⁸	8.06 × 10⁻⁸
n	2.52	2.55
P	−0.12	−0.119
E₁₁ (GPa)	138.0	127.6

PMCs: polymer matrix composites.

Ideally, the nonlinear model constants of the matrix should be determined from test data obtained from the in situ matrix that has gone through the same processing cycle as the composite. In the absence of such data they should be determined from tension test data obtained from 90° unidirectional laminates in which the matrix experiences the full applied load. Since such data were not available in the present case, data obtained from alternate matrix dominated load cases was used in the parameter estimation. The model constants in equation (12) for 3501-5 and PEEK were determined from tension data obtained from the 45° angle-ply and the (±45)_S laminate, respectively. The constituent properties and the estimated values of the material specific parameters for these material systems are given in Table 1. The fits of the data using these parameters in the above model for the two PMCs are shown in Figures 2 and 3.

Figure 2.

Modeling nonlinear response of AS4/3501-5 PMC (RT).

Figure 3.

Modeling nonlinear response of AS4/PEEK PMC (RT). PMC: polymer matrix composite; RT: room temperature.

The PMC nonlinear model described above does include consideration of residual stresses. Thermoset composites are usually processed at relatively lower temperatures than thermoplastic composites and consequently the residual stresses are expected to be significantly lower. In the present study, we assumed very small residual stress (0.3 MPa) in the matrix in the AS4/3501-5 composite. A stress-free temperature $(T p)$ of 251.7℃ has been reported in reference³ for the AS4/PEEK composite. The associated residual stress in the matrix was determined to 37.9 MPa.

The predicted tension response of several off-axis AS4/3501-5 laminates along with corresponding test data is shown in Figure 4. The predictions, obtained using the $J 2 A$ theory and the material constants in Table 1, capture the observed nonlinear behavior in the stress–strain data reasonably well. One consequence of plastic deformation in the composite is that the apparent Poisson’s ratio of the material increases. This can be determined from the finite element analysis (FEA) solutions. As shown in Figure 5, the model predicts the change in the plastic Poisson’s ratio, the ratio of the plastic strains, with the angle of the unidirectional ply. Note that in all cases data from a single off-axis laminate were used to predict the response of the other test cases.

Figure 4.

Prediction of nonlinear response of angle-ply AS4/3501-5 PMC (RT).

Figure 5.

Prediction of plastic Poisson’s ratio in angle-ply AS4/3501-5 PMC (RT).

In the off-axis composite, the stress in the ply coordinate system is affected by the angle of the ply. Before the onset of matrix yielding and damage in the composite the variation of stress components with fiber angle is shown in Figure 6. In the figure all the stresses have been normalized with the applied stress. As expected it is seen that the longitudinal stress, $σ 11$ , is maximum in the 0^o unidirectional composite and minimum in the 90° composite; the transverse stress, $σ 22$ , shows the exact opposite behavior. The shear stress, $σ 12$ , is dominant at the in-between angles and peaks in the 45° composite.

Figure 6.

Predicted ply stresses in angle-ply AS4/3501-5 PMC (RT).

In the present model, expression (8b) for $N 11$ is used to determine the matrix stress $(σ 11 m)$ from the fiber direction composite stress $(σ 11)$ . In addition, the transverse stress and the shear stress components are considered to be fully seen by the matrix. The matrix stress components thus considered by the model and corresponding to the ply stresses in Figure 6 are shown in Figure 7. The $J 2$ equivalent stress in the matrix $(σ e m)$ given by

σ e m = \sqrt{(σ 11 m) 2 + (σ 22 m) 2 + 3 (σ 12 m) 2 - σ 11 m σ 22 m}

(13)

is also plotted in the figure. In the three-dimensional FEA of the off-axis specimens all the out-of-plane stresses were relatively insignificant in magnitude.

Figure 7.

Predicted matrix stresses in angle-ply AS4/3501-5 PMC (RT).

Figure 7 shows that for very small angles, since the fiber carries most of the stress, the equivalent stress in the matrix is relatively small. With increase in the ply angle the matrix equivalent stress increases till it peaks at about 50° to 60° after which the magnitude of the equivalent stress decreases gradually; the relative change in the matrix equivalent stress after 45° is not very significant. Thus according to these results the model predicts that the inelastic strain rate in the matrix, determined by equations (5) and (12), would peak between 50° and 60° and that after about 45° it would not change significantly. This prediction is consistent with the data in Figure 4.

By comparing Figures 6 and 7 it can be seen that the transverse $(σ 22 m)$ and shear $(σ 12 m)$ stress in the matrix are equal in value to the corresponding stress in the composite as is assumed in the model. The longitudinal stress in the composite $(σ 11)$ is shared between the fiber and matrix and so the matrix stress in the direction of the fibers $(σ 11 m)$ is not equal to $σ 11$ . Note that when the load is applied transverse to the fibers (90° off-axis case), the composite stress $σ 11$ is zero but the corresponding stress in the matrix is non-zero. The value of the latter depends both on the Poisson effect caused by the stress applied in the 90° direction and on the matrix residual stress in the 0° direction.

The above example shows that the micromechanics-based model described above is able to predict composite response under multiaxial stresses. However, in all the above test cases the stress state in the specimen is uniform. We therefore next consider the modeling of PMC structural configurations that contain stress gradients.

Tension data from tapered 30° and 45° off-axis specimens and 0° and ±45 rectangular specimens with a central hole have been reported in reference.³ In this work the specimens, tested at room temperature, were instrumented with strain gages and the strains recorded. The geometry of these AS4/PEEK specimens along with the locations of the strain gages is described in the reference. Tension data from standard ±45 specimens, also reported in reference,³ were first used to determine the nonlinear material constants listed in Table 1.

The planar view of the three-dimensional FEA mesh of the tapered and the hole specimens are shown in Figure 8(a) and (b), respectively. Along the out-of-plane direction each ply of the composite was modeled using one layer of elements. The FEA model of the tapered and hole specimens that were analyzed each contained 4299 nodes and 800 twenty-node brick elements.

Figure 8.

FEA mesh of PMC: (a) tapered specimen and (b) hole specimen.

The predicted nonlinear stress–strain behavior of the tapered 30° and 45° specimens is compared with the corresponding test data in Figure 9. In this specimen the stress state varies, both across the width (due to the ply orientation), and also along its length (due to the geometry). The strains predicted at the gage location are seen to agree with the test data.

Figure 9.

Strains in angle-ply tapered AS4/PEEK specimen (RT).

Stress analysis of the hole specimens shows that the largest composite stress is in the loading direction and situated at the top-edge of the hole in the net-section plane. The presence of the central hole in these specimens results in the stress concentration and stress gradients seen close to the hole. Comparison of the hoop strain predicted at the strain gage locations close to the hole in the 0° and ±45 specimens with the corresponding test data is shown in Figures 10 and 11, respectively. The strain gages were located at different locations in the two specimens: in the ±45 specimen, the gage was located where the load–direction stress is maximum, while in the 0° specimen, two gages were located at an angle of about 56° to the load axis. In both specimens ,the model is able to predict the strains at the gage locations reasonably well.

Figure 10.

Hoop strain in unidirectional AS4/PEEK hole specimen (RT).

Figure 11.

Hoop strain in (±45)_S AS4/PEEK hole specimen (RT).

The $J 2$ equivalent stress in the composite and in the matrix at the edge of the hole in the unidirectional and ±45 specimens as determined by the model, and corresponding to a value of applied stress that does not cause inelastic deformation, is shown in Figure 12. It is seen that in both specimens, 0° and ±45, the composite stress is maximum at the top of the hole (ϕ = 90°). In the 0° specimen, the largest component of the matrix stress is the shear stress. Along the edge of the hole the shear stress peaks at about ϕ = 77° and this is also the location where the matrix equivalent stress is maximum. In the ±45 specimen all the matrix stress components peak at ϕ = 90° and therefore this is also the location where the matrix equivalent stress is the largest. Thus according to equations (5) and (12), these are also the locations where the inelastic strain in the matrix will be the greatest in the two types of specimens.

Figure 12.

Variation of matrix stress along hole edge.

Figure 13 shows the variation of the maximum stress-concentration defined as

K t = \frac{σ max}{σ 0}

(14)

in the unidirectional and ±45 hole specimens as they are loaded up. Here

σ max

and

σ 0

are the maximum stress in the direction of the load in the composite and the average remote stress applied on the specimen, respectively. As mentioned earlier, the maximum composite stress is in the direction of the load at the top-edge of the hole. Figure 12 shows that in the unidirectional specimen the matrix stress at this location is relatively small and so locally not much inelastic deformation is expected. Thus the local stress concentration is not expected to decrease significantly as the remote load is increased; in fact Figure 13 shows a slight increase in the stress concentration in the 0° specimen due to the Poisson effect. In the ±45 specimen, the equivalent stress in the matrix is also the highest at the top-edge of the hole (Figure 12) and therefore the local material is expected to see the largest amount of inelastic strain. As the material undergoes inelastic deformation, some of the load will be transferred to the neighboring region. Thus, Figure 13 shows that in the ±45 specimen the stress-concentration decreases as the remote load is increased.

Figure 13.

Change in hole stress-concentration in composite.

Conclusion

The present study has introduced a micromechanics-based approach to model the inelastic response of PMCs. In this approach, once the matrix material is characterized the model can be used to predict structural response of composites with other laminate configurations and geometry. The results of the present study demonstrate the application of the model to predict structural response of PMCs in the presence of stress gradients.

Footnotes

Funding

This research received no specific grant from any funding agency in the public, commercial or not-for-profit sectors.

Conflict of interest

None declared.

References

Springer

. Time dependent behavior of laminated composite at elevated temperature. J Compos Mater 1989; 23: 1159–1197.

Sun

Chen

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

Sun

Yoon

. Elastic-plastic analysis of AS4/PEEK composite laminate using a one-parameter plasticity model. J Compos Mater 1992; 26: 293–308.

Gates

Sun

. Elastic/viscoplastic constitutive model for fiber, reinforced thermoplastic composites. AIAA J 1991; 29(3): 457–463.

Weeks

Sun

. Modeling non-linear rate-dependent behavior in fiber-reinforced composites. Compos Sci Technol 1998; 58: 603–611.

Spathis

Kontou

. Non-linear viscoplastic behavior of fiber-reinforced polymer composites. Compos Sci Technol 2004; 64: 2333–2340.

Kawai

Masuko

. Creep behavior of unidirectional and angle-ply T800H/3631 laminates at high temperature and simulations using a phenomenological viscoplasticity model. Compos Sci Technol 2004; 64: 2373–2384.

Yokozeki

Ogihara

Yoshida

. Simple constitutive model for nonlinear response of fiber-reinforced composites with loading-directional dependence. Compos Sci Technol 2007; 67: 111–118.

Sun

Chen

. A micromechanical model for plastic behavior of fibrous composites. Compos Sci Technol 1991; 40(2): 115–129.

10.

Tsai

J-L

Chen

K-H

. Characterizing nonlinear rate-dependent behaviors of graphite/epoxy composites using a micromechanical approach. J Compos Mater 2007; 41(10): 1253–1273.

11.

Ahmad

Newaz

Nicholas

Analysis of characterization methods for inelastic composite material deformation under multiaxial stresses. In: Kalluri

Bonacuse

(eds). ASTM STP 1387-EB: multiaxial fatigue and deformation: testing and prediction 2000, pp. 41–53.

12.

Gatewood

. Thermal stresses, New York, NY: McGraw-Hill, 1957.