Micromechanical Modeling of Fiber Reinforced Composites Based on Elastoplasticity and its Application for 3D Braided Glass/Kevlar Composites

Ji Hoon Kim,1 Hansun Ryou,1 Myoung-Gyu Lee,2 Kwansoo Chung,3 Jae Ryoun Youn,3 Tae Jin Kang3 1 Department of Materials Science and Engineering, Seoul National University, Gwanak-gu, Seoul 151-744, Korea 2

Department of Materials Science and Engineering, The Ohio State University, Columbus, Ohio 43210

3

Department of Materials Science and Engineering, Intelligent Textile System Research Center, Seoul National University, Seoul, 151-744, Korea

Micromechanical modeling to calculate the mechanical properties of fiber reinforced composites is proposed. To describe the mechanical behavior of the yarn and the matrix, which are the main constituents of fiber reinforced composites, the elastoplastic constitutive law was adopted. In particular, anisotropic elastoplasticity based on Hill’s orthotropic yield function and anisotropic kinematic hardening was utilized for the yarn, while the isotropic elastoplastic constitutive law was applied for the matrix. The effective properties of the unit cell in fiber reinforced composites were then calculated based on the finite element method. For verification, the method was successfully applied for 3D braided glass/Kevlar fiber reinforced composites in both linear elastic and nonlinear inelastic ranges. POLYM. COMPOS., 28:722–732, 2007. ª 2007 Society of Plastics Engineers

INTRODUCTION Fiber reinforced composites have been widely used in aerospace, automobile, and defense industries because of their high stiffness to weight ratios, good impact resistance, high damage tolerance, and good durability. Fiber reinforced composites with textile preforms have good transverse mechanical properties for a single ply, good impact, and crack resistance in the out-of-plane direction as well as dimensional stability. In particular, 3D braided composites have demonstrated significant improvement in

Correspondence to: Kwansoo Chung; e-mail: [email protected] Contract grant sponsor: KOSEF; contract grant number: R11-2005-065. DOI 10.1002/pc.20357 Published online in Wiley InterScience (www.interscience.wiley.com). C 2007 Society of Plastics Engineers V

POLYMER COMPOSITES—-2007

damage tolerance [1]. Also, braids can be easily made into various shapes, including tubes [2, 3]. In addition, the hybridization of yarns makes braided composites have the wide ranges of mechanical properties [4, 5]. Their complex microstructures, however, make the prediction of mechanical properties difficult, especially in the nonlinear inelastic range beyond the linear elastic range. Efforts to understand the mechanical behavior of fiber reinforced composites follow two extremes: continuum and micromechanical approaches. In the continuum approach, the composite is considered as a homogeneous material having a uniform (average) property, ignoring the property difference between each constituent [6]. The continuum approach is common for structural analysis involving numerous unit-cells. On the other hand, the micromechanical approach considers the composite as the combination of various materials and derives overall properties utilizing homogenization procedures based on the individual properties of constituents. The micromechanical approach is useful to optimize the design parameters of composites associated with textile structures, such as the fiber volume fraction, to meet the target properties of composites. Many have been working on the micromechanical analysis of fiber reinforced composites. In the response average method (or the volume average method), stresses and strains are obtained by volume averaging the stiffness or compliance matrices. The response average method has been used to predict the mechanical properties of woven [7], 2D biaxially braided [8, 9], and triaxially braided composites [10–15]. The micromechanical analysis model based on the classical laminate plate theory (CLPT) has also been developed to predict the stiffness and strength properties of textile composites, including the mosaic, fiber undulation,

FIG. 1. The schematic hierarchical structure of fiber reinforced composites.

and fiber bridging models [16–18]. Carey et al. [19–21] applied the CLPT to predict the longitudinal elastic modulus of the 2D braided fiber composite. Raju and Wang [22] improved the CLPT by removing the assumptions used to derive the fiber undulation model and the bridging model. Huang and coworkers [23, 24] applied the bridging model to simulate the mechanical properties of composites reinforced with woven and braided fabrics. Redman and Douglas [25] used a unique combination of the basic rule-of-mixture and CLPT for the prediction of tensile elastic properties of braided composites. Falzon et al. [26] extended the point-wise lamination approach to analyze the compaction problem for plain woven composites. The micromechanical analysis model is known to give a good estimate of in-plane stiffness but is less accurate in predicting shear stiffness, Poisson’s ratio, and out-of-plane properties. The Mori-Tanaka micromechanical method has been used to predict the effective properties of composites because the method can give explicit formulae for the stiffness matrix. Gommers et al. [27] applied the method to calculate the elastic properties of textile composite materials. The finite element analysis has also been used to predict the mechanical properties of fiber reinforced composites, in which the nonlinear mechanical properties of fiber reinforced composites were accounted for by introducing nonlinear constitutive models for constituent materials such as nonlinear moduli [28, 29] and failure criteria [30]. Since the finite element analysis of 3D braided composites is difficult because of their complex microstructures [31], Bigaud and Hamelin [32] applied a multiscale energetic approach to woven fabric and 3D braided composites by introducing a three-dimensional aggregate of subcells to represent a fabric unit cell. Bogdanovich and Pastore [33] developed a procedure for the structural analysis of three-dimensional elasticity problems by incorporating local displacements of subelements. Goyal et al. [34] applied three-dimensional micromechanical finite element models to investigate the effect of various paramDOI 10.1002/pc

eters such as the braid angle, the waviness ratio, material properties, and the cross-sectional shape on the effective engineering properties of 2
2 braids. Ivanov and Tabiei [35] developed a computational material model for the plain-woven fabric composite for use in finite element analysis. Queka et al. [36] investigated the tensile and compressive properties of the carbon 2D triaxially braided composite based on the micromechanical finite element method. Tan et al. [37] developed analytical and finite element models based on sinusoidal beams to investigate the elastic and failure properties of plain woven composites. In this paper, a micromechanical procedure based on the finite element analysis is developed utilizing elastoplasticity, which is capable of handling the anisotropic and nonlinear behavior of composites. In particular, anisotropic elastoplasticity based on Hill’s orthotropic yield function and anisotropic kinematic hardening was utilized for the yarn (impregnated in the matrix), while the isotropic elastoplastic constitutive law was applied for the matrix. The effective properties of the unit cell in fiber reinforced composites were then calculated based on the finite element method. For verification, the method was applied for 3D braided glass/Kevlar fiber reinforced composites and the predicted properties of composites were compared with experimental results for simple tension tests performed along various sample directions. Note that the application of the proposed approach is limited to a stretching mode without unloading.

MECHANICAL PROPERTIES OF CONSTITUENTS Fiber reinforced composites are composed of yarns, which are impregnated in the matrix, and the matrix, which fills the space between yarns. The hierarchical structure of composites is shown in Fig. 1. Here, the mechanical properties of the yarn and the matrix are modeled using the elastoplastic constitutive law. POLYMER COMPOSITES—-2007 723

Yarn Properties

In elastoplasticity, the strain increment is additively decoupled as

Yarns (or one of unidirectional fiber reinforced composites) have high anisotropic material properties. Yarns usually have much higher longitudinal stiffness than transverse stiffness because of the alignment of fibers. Also, the stress–strain behavior of yarns is generally nonlinear because of microcracks formed during stretching. Therefore, anisotropic elastic constitutive laws, which have been used to model the behavior of fiber reinforced composites, are not sufficient to model the nonlinear behavior of yarns. In this paper, elastoplastic constitutive equations are used to properly describe the anisotropic and nonlinear behavior of yarns.

0

1

2

1n223 2 6 E2 D 6 n21 ð1þn23 Þ 6 E2 D 2 6 6 n21 ð1þn23 Þ 6 E2 D 2 6

ds11 B ds22 C C B B ds33 C C B B ds23 C ¼ C B @ ds31 A 6 6 4 ds12

n21 ð1þn23 Þ E22 D 1n21 n12 E1 E2 D n23 þn21 n12 E1 E2 D

0 0 0

(2)

where C is the elastic stiffness matrix. If yarns are assumed transversely isotropic, Eq. 2 can be written with five independent elastic constants in the symmetric matrix form as

3

0

0

0

0

0

0

0

0

0 0 0

E2 2ð1þn23 Þ

0 G12 0

0 0 G12

F ¼ Fðs22  s33 Þ2 þ Gðs33  s11 Þ2 þ Hðs11  s22 Þ2

0 0

1 0 7 de11 7B 7 B de22 C C 7B 7 B de33 C C 7B 7 B 2de23 C C 7@ 7 2de31 A 5 2de12

(3)

0

F ¼ Fðs22  s33 Þ2 þ Gðs33  s11 Þ2 þ Gðs11  s22 Þ2 þ 2ð2F þ GÞs223 þ 2Ms231 þ 2Ms212  s2 ¼ 0

ð5Þ

The three material parameters, F, G, and M, can be Y Y determined from two tensile yield stresses sY 22 ð¼ s33 Þ, s11 Y and one pure shear yield stress sY ð¼ s Þ, respectively. 12 31 Beyond the initial yield point, yarns show different nonlinear stress–strain behavior in different loading directions. To account for the directional difference in hardening, the anisotropic kinematic hardening law based on the modified Chaboche back stress evolution rule is introduced [6, 41]: da ¼ C1 

ðs  aÞ de  G2  ade aiso

(6)

where s is the Cauchy stress, a is the back stress measuring the position of the yield surface, siso is the size of the yield surface, and de is the effective plastic strain increment, which is determined from the consistency condition and the associated flow rule:

ð4Þ

where sij is the component of the Cauchy stress and s is the effective yield stress representing the yield surface size, while F, G, H, L, M, and N are material constants characterizing anisotropic yielding behavior. For transversely isotropic materials, G ¼ H, M ¼ N, and L ¼ 2F þ G so that Eq. 4 becomes 724 POLYMER COMPOSITES—-2007

dr ¼ C  dee

0

0 0 0

(1)

where de, dee, and dep are total, elastic and plastic strain increments, respectively. The Cauchy stress increment ds is given by the generalized Hooke’s law as

n21 ð1þn23 Þ E22 D n23 þn21 n12 E1 E2 D 1n21 n12 E1 E2 D

where D ¼ ð1 þ n23 Þð1  n23  2n12 n21 Þ=E22 E1 , while E, G, and n are Young’s modulus, the shear modulus, and Poisson’s ratio, respectively (with E2 n12 ¼ E1 n21 ). The subscripts 1, 2, and 3 denote the three orthogonal axes of the fiber coordinate system defined in Fig. 1. The values of elastic stiffness constants in Eq. 3 are calculated using micromechanical models. For yarns made of the isotropic fiber and the isotropic matrix, the equations developed by Hashin and Rosen [38] are used for calculation. For yarns made of the transversely isotropic fiber and the isotropic matrix, the equations suggested by Chamis [39] are used. Examples of the isotropic fiber, the transversely isotropic fiber, and the isotropic matrix are the glass fiber, the Kevlar fiber, and the epoxy, respectively. These equations used for micromechanical models are summarized in Appendix A. Hill’s orthotropic yield function is used to account for anisotropic yielding or the linearity limit of yarns [40], which is given as

þ 2Ls223 þ 2Ms231 þ 2Ns212  s2 ¼ 0

de ¼ dee þ dep

dep ¼ de

qsiso qr

(7)

In Eq. 6, C1 and C2 are the fourth-order tensors containing parameters to be determined from the hardening curves. Note that when Eq. 5 is used along with Eq. 6 for the kinematic hardening law, s becomes siso by replacing DOI 10.1002/pc

s with (sa). Note that Eq. 6 allows only six proportional loading directions along which the hardening behavior can be prescribed. If off-diagonal terms are ignored for simplification, Eq. 6 becomes, in the component form, daij ¼ ðgij nij  hij aij Þde ðno summation on ijÞ

(8)

where n ¼ ðs  aÞ=siso , while aij is the components of a, gij, and hij are the nonvanishing components of G1 and G2, respectively. When all gij’s and hij’s share the same values, respectively, the traditional isotropic kinematic hardening law is recovered. It is further assumed here that g22 ¼ g33 ¼ g23 , g13 ¼ g12 , h22 ¼ h33 ¼ h23 , and h13 ¼ h12 to be consistent with the condition of the initial transverse isotropy and the plastic work equivalence principle (see Appendix B for derivation). Since Eq. 8 is the first order ordinary differential equation with constant coefficients for the back stress, the following relationship is obtained for monotonously proportional loading: Y sij ¼ sY ij þ aij ¼ sij þ

FIG. 2. The 3D circular braiding machine. [Color figure can be viewed in the online issue, which is available at www.interscience.wiley.com.]

p gij nij ð1  ehij eij =mij Þ hij

(9)

where m ¼ qsiso =qðs  aÞ and sY ij is the initial yield stress. Therefore, the constants gij and hij are obtained from the curve fitting of measured hardening curves.

FIG. 3. Carrier movements during the 4-step braiding process. DOI 10.1002/pc

POLYMER COMPOSITES—-2007 725

FIG. 4. 3D braided glass/Kevlar fiber reinforced composite samples: (a) the glass side and (b) the Kevlar side. [Color figure can be viewed in the online issue, which is available at www.interscience.wiley.com.]

Matrix Properties The matrix material such as the epoxy shows isotropic but nonlinear stress–strain behavior. Therefore, an isotropic elastoplastic constitutive law with the isotropic elastic stiffness and the Mises yield function was used to describe the matrix behavior. As for hardening, the isotropic hardening rule by Voce [42] was used: Be s ¼ sY Þ m þ Að1  e

(10)

where sY m is the yield stress of the matrix, e is the accumulative effective plastic strain, while A and B are material constants determining the isotropic hardening curve. On the basis of the elastoplastic theory, the constitutive laws of the yarns and the matrix were implemented into the general purpose finite element program ABAQUS/ Standard using the material user subroutine UMAT [43].

obtain the mechanical responses of composites, boundary conditions are imposed on the unit-cell descretized by three-dimensional solid elements. Local stresses and strains satisfying the equilibrium condition and the geometric compatibility condition are calculated by solving the boundary value problem. Since the unit cell consists of various constituents, resulting solution is inhomogeneous and the overall stress and strain are obtained by averaging local stresses and strains, i.e., Z Z 1 1 Sij ¼ sij dV and Dij ¼ eij dV (11) V V V V where V is the volume of a unit-cell, while Sij and Dij are the overall stress and strain components, respectively. The average elastic constants and nonlinear stress–strain behavior can be obtained using the overall stress and strain with appropriate boundary conditions.

HOMOGENIZATION On the basis of the behavior of constituent materials, the overall average properties of composites are calculated utilizing the following procedure. First, a unit-cell representing the composite structure is identified. Then, the average mechanical properties of composites are obtained by applying physical assumptions to the unit-cell. There are various homogenization methods depending on physical assumptions applied to the unit-cell. In this paper, a homogenization method based on the finite element analysis was used. For comparison, other homogenization methods such as the response average (in stiffness and compliance) method [44] and the method of cells [45] were also considered for the elastic behavior. Finite Element Method (FEM) Finite element method (FEM) can be effectively utilized to predict the properties of composites [46]. To 726 POLYMER COMPOSITES—-2007

Response Average (in stiffness and compliance) Method If all local strains in a unit-cell are assumed identical (the iso-strain condition), the overall elastic stiffness matrix C is obtained by averaging local stiffness matrices Cl.

TABLE 1. Basic properties of the fibers and the matrix [48, 52–54]. Parameter E1 (GPa) n12 n23 E2 (GPa) G12 (GPa) Density (g/cm3) a

S2 glass fiber

Kevlar 29 fiber

Epoxy resin

86.9 0.20 – – – 2.45

62.0 0.36 0.36a 2.4 1.8 1.44

2.1 0.37 – – – 1.12

Assumed.

DOI 10.1002/pc

FIG. 5. The measured and calculated longitudinal stress–strain curves of the glass and Kevlar yarns in the (a) longitudinal and (b) transverse directions.

1 C¼ V

Z Cl dV

(12)

V

Also, if all local stresses in a unit-cell are assumed identical (the iso-stress condition), the effective elastic compliance matrix S is obtained by averaging local compliance matrices Sl. Z 1 Sl dV (13) S¼ V V In calculating average values using Eqs. 12 and 13, the matrices of local stiffness or compliance with respect to a local (or fiber) coordinate system should be transformed into the global coordinate system. In this method, the equilibrium and compatibility conditions are ignored and responses by the constitutive law are averaged.

TABLE 2. Material parameters of the glass and Kevlar yarns.

E1 (GPa) E2 (GPa) n12 n23 G12 (GPa) Y s11 (MPa) Y s11 (MPa) F G M g11 (MPa) gij (MPa, otherwise) h11 hij (otherwise)

DOI 10.1002/pc

Glass yarn

Kevlar yarn

74.2 24.0 0.22 0.30 7.6 75.7 32.3 81.2 0.5 1.5 536,538 0 563.7 0

53.0 2.4 0.36 0.36 1.6 56.6 8.1 48.3 0.5 1.5 194,986 0 220 0

Method of Cells (MC) The method of cells (MC) was first developed by Aboudi [45] and extended to the general 3D case by Paley and Aboudi [47]. Later, the method was formulated based on local stress and strain tensors by Bigaud and Hamelin [48]. In the MC, a unit-cell is discretized into the subcells of rectangular parallelepiped shapes, in which each subcell has its own local strain and local stress (as in the finite element method). Local strains and stresses are obtained by solving the boundary value problem, considering the displacement continuity (or geometric compatibility) condition and traction continuity (or equilibrium) condition but only in an average sense (therefore, not as strictly as in the finite element method). Therefore, the method of cells is computationally more efficient when compared with the finite element method but its results may not be as accurate as the finite element method [49]. VERIFICATION Material Preparation The preform of the 3D braided glass/Kevlar fiber reinforced composite was fabricated by the 4-step braiding process using the 3D circular braiding machine with 2,014 carriers and 104 pistons shown in Fig. 2. Glass and TABLE 3. Material parameters for the epoxy resin. Paramters

Value

smY (MPa) A (MPa) B

25.0 24.5 402.4

POLYMER COMPOSITES—-2007 727

FIG. 6. The measured and calculated stress–strain curves of the epoxy.

Kevlar fibers were used as reinforcing fibers. The movements of carriers during the 4-step braiding process are shown in Fig. 3. The white circles denote the carriers with glass yarns. Two carriers are marked ‘‘a’’ and ‘‘b’’ to show the carrier movements clearly. The black circles are the carriers with Kevlar yarns, which move only along the circumferential direction, not along the radial direction. Ultimately, Kevlar fibers align along the axial direction, while glass fibers are angled about 278 with the axial direction. Also, one side is rich with Kevlar fibers (therefore, called the Kevlar side, while the other side is called the glass side, for convenience). At Step 1, the carriers move along the radial direction except for the carriers with Kevlar fibers. At Step 2, the carriers move along the circumferential direction except for the inner- and outermost carriers. In Steps 3 and 4, the carriers move in the same manner with Steps 1 and 2 but to opposite directions. The cylindrical preform with a radius of 30 mm prepared by the circular braiding machine was cut along the axial direction and flattened to be a rectangular one. Then, the preform was placed in the mold and the epoxy resin was transferred. The mold was heated at 70–1308C for 8 h for curing. After the curing process, the final

FIG. 7. Yarn paths and unit-cell of the braided composite.

728 POLYMER COMPOSITES—-2007

FIG. 8. Subcells of the braided composite: subcell (a) A, (b) A0 , (c) B and (d) B0 .

FIG. 9. Finite element model of the unit-cell. [Color figure can be viewed in the online issue, which is available at www.interscience. wiley.com.]

DOI 10.1002/pc

ters, the perfect plastic behavior was assumed after the initial yielding (so that aij ¼ 0 without hardening). The material parameters for the glass and Kevlar yarns are listed in Table 2. The material parameters for the epoxy resin were determined from the simple tension test, and results are summarized in Table 3. Using the material parameters of the elastoplastic constitutive law, the stress–strain curves were recalculated as shown in Figs. 5 and 6. The results show that not only linear elastic but also nonlinear behavior of the glass, Kevlar yarns, and the epoxy matrix were successfully described using the elastoplastic constitutive law.

FIG. 10. Dimensions of the test specimen.

sample with 3.3 mm thickness was obtained as shown in Fig. 4. The fiber volume fractions were measured for the final sample. The glass fiber volume fraction was determined from the net weight of the glass fibers after burning off the Kevlar fibers and the epoxy. The Kevlar fiber volume fraction was estimated based on the fiber feed rate of the braiding process. The volume fractions were 43 and 7% for the glass and Kevlar fibers (not yarns), respectively.

Material Characterization On the basis of the elastic constants of the glass, Kevlar fibers, and the epoxy resin as listed in Table 1, the elastic constants for yarns were calculated using equations listed in Appendix A. Note that the fiber volume fraction of 85% in yarns was assumed for both the glass and Kevlar yarns. The material parameters for the yield criterion and the hardening law of the yarns were determined from experiments. Since preparing unidirectional composite samples with the high (85%) fiber volume fraction was difficult, the stress–strain curves were obtained by extrapolating those of the low (20%) fiber volume fraction composites using the rule-of-mixture. In the longitudinal direction, the initial yield point was identified as the limit of linearity as shown in Fig. 5. In the transverse direction, the initial yield point was assumed to be the failure point. The initial shear yield stress was obtained assuming the pffiffiffi Y Mises material property (sY ¼ s = 3 ). On the basis of 12 11 these initial yield stresses, the three yield criterion parameters were calculated using Eq. 5. As for the anisotropic hardening parameters, g11 and h11 were calculated by fitting the longitudinal hardening curve using Eq. 9. For the transverse hardening parame-

Unit-Cell Structure The unit-cell structure of the composite was defined by analyzing the braiding process. First, points which yarns must pass during the braiding process were identified, and they were connected to form a yarn path considering other yarn paths. After the unit-cell was identified, the uni-cell was further divided into seven subcells of four kinds as shown in Fig. 7. The dimension of the unit-cell was 20
10
3:3 ðmm3 Þ. Figure 8 shows schematic glass yarn paths in each subcell. Note that subcell A0 can be obtained by translating subcell A along the x (axial-) direction, while subcell B0 can be obtained by transforming subcell B using the xy-plane as the mirror plane and translating it along the y direction. This idealized yarn paths may not be realistic. Deviations of actual structures from the ideal ones for woven and braided composites have been investigated elsewhere [50]. For the more accurate description of yarn paths, the unit cell structure can be obtained by analyzing the photomicrographs of braided samples [51]. After yarn paths are determined, the unit-cell was discretized with 8,000 3D brick elements as shown in Fig. 9. Elements are classified into the glass yarn type, the Kevlar yarn type, and the epoxy matrix type. Each element has its own orientation along which material properties are defined.

RESULTS AND DISCUSSION Simple tension tests were carried out for the braided composite along the axial (x-) and transverse (y-) directions by the standard procedure ASTM D3039-76 using the 10-ton test machine Instron 8516 system. The specimen geometry is shown in Fig. 10. The axial elastic mod-

TABLE 4. Measured and calculated elastic constants of the 3D braided composite. Method Response average in stiffness Response average in compliance MC FEM EXP

DOI 10.1002/pc

Ex (GPa)

Ey (GPa)

Ez (GPa)

nxy

nxz

nyz

Gxy (GPa)

Gzx (GPa)

Gyz (GPa)

5.1 32.4 16.9 24.3 24.4

4.0 12.4 9.6 9.5 9.3

4.1 13.6 13.1 11.4 –

0.38 0.54 0.58 0.42 –

0.35 0.14 0.15 0.17 –

0.35 0.27 0.23 0.20 –

2.8 8.9 3.1 4.4 –

2.7 5.6 1.8 3.6 –

2.9 9.2 2.0 3.4 –

POLYMER COMPOSITES—-2007 729

FIG. 11. Measured and calculated stress–strain curves of the 3D braided composite in the (a) axial and (b) transverse directions.

ulus Ex and the transverse elastic modulus Ey were obtained from the initial slopes of the axial and transverse stress–strain curves, respectively. Using the homogenization method based on the finite element analysis, the orthotropic elastic constants were calculated for the 3D glass/Kevlar braided fiber reinforced composite. For comparison, the response average method in stiffness and compliance and the method of cells were also used to calculate the elastic constants. The calculated and measured elastic constants are summarized in Table 4. The FEM and MC predictions of elastic constants were between the values predicted by the response average method in stiffness and that in compliance, which correspond to the upper and lower bounds, respectively. The axial modulus, Ex, predicted by the FEM showed good agreement with the experimental result, while the MC underestimated the axial modulus. As for the transverse modulus, Ey, both the FEM and MC results showed good agreement with the measured value. In Fig. 11, the stress–strain curves predicted by FEM are compared with measured curves in the axial and transverse directions, which also show reasonably good agreement.

stitutive law was applied based on the Mises yield function. The unit-cell was constructed and the finite element method was used to calculate the homogenized properties of the unit-cell. For verification, the 3D braided glass/ Kevlar fiber reinforced composite was fabricated and the method was applied following the proposed procedure. The results showed that not only the elastic constants but also anisotropic nonlinear stress–strain behavior gave good agreement with experiment results, confirming that the method can be effectively utilized in designing (or optimizing) the microstructures of fiber reinforced composites so as to meet the target properties of final products. APPENDIX A: EQUATIONS FOR THE ELASTIC CONSTANTS OF YARNS Equations for the elastic constants of yarns by Hashin and Rosen [38] are: E 1 ¼ E f V f þ Em V m þ

E2 ¼ E3 ¼

SUMMARY A homogenization method based on the finite element analysis incorporating elastoplastic modeling was developed to effectively predict the mechanical behavior of fiber reinforced composites. In this paper, the individual mechanical properties of the yarn and the matrix were described by applying elastoplastic constitutive laws. As for the yarn properties, the elastic parameters were calculated by applying micromechanical models, while the anisotropy and nonlinearity of the stress–strain hardening relation were described using Hill’s orthotropic yield criterion and the anisotropic kinematic hardening rule, respectively. For the matrix, the isotropic elastoplastic con730 POLYMER COMPOSITES—-2007

4V f V m ðnf  nm Þ Vm kf

þ kVmf þ G1m

4kt Gt   4k v2 kt þ Gt 1 þ Et 1112

2

(A1)

(A2)

V m Gm þ ð1 þ V f ÞGf (A3) ð1 þ V f ÞGm þ V m Gf   V f V m ðnf  nm Þ k1m  k1f ¼ nf V f þ nm V m þ Vm Vf 1 kf þ km þ Gm

G12 ¼ G13 ¼ Gm

n12 ¼ n13

(A4)

n23 ¼

E2 1 2Gt

(A5)

DOI 10.1002/pc

where

kt Gt ¼

kf ¼ Ef =2ð1  nf  2n2f Þ

(A6)

km ¼ Em =2ð1  nm  2n2m Þ

(A7)

km kf þ ðV f kf þ V m km ÞGm ¼ V m kf þ V f km þ G m

3 2 ða þ bm V f Þð1 þ rV f Þ  3V f V m b2m Gm 3 2 ða  V f Þð1 þ rV f Þ  3V f V m b2m

a ¼ ðg þ bm Þ=ðg  1Þ bm ¼

1 3  4nm

(A8)

(A9) (A10)

r ¼ ðbm  gbf Þ=ð1 þ gbf Þ

(A13)

g ¼ Gf =Gm

(A14)

where V is the volume fraction. Note that the subscripts f and m indicate the fiber and the matrix, respectively, while the constants with the subscript t and the superscript * are dummy constants used for Eqs. A2 and A5. Equations for the elastic constants of yarns by Chamis [39] are as follows:

G12 ¼ G13 ¼

p g22 n22 ð1  eh22 e22 =m22 Þ h22

(A15)

(B2) (B3)

where m ¼ qs=qðs  aÞ and n ¼ ðs  aÞ=siso . For the pure shear test on the 2–3 plane, s23 ¼ sY 23 þ

(A11) (A12)

E2 ¼ E3 ¼

s22 ¼ sY 22 þ

dep22 ¼ dem22

1 bf ¼ 3  4nf

E1 ¼ E1f V f þ Em V m

Here, a is the effective back stress whose value is obtained from the effective stress by replacing the Cauchy stress with the back stress. For the simple tension test in the 2 direction,

p g23 n23 ð1  eh23 e23 =m23 Þ h23

dep23 ¼ dem23

(B4) (B5)

From Eqs. B3 and B5 dep23 ¼

m23 m22 dep22

(B6)

Applying Eq. B6 for the plastic work equivalence principle shown in Eq. B1 leads to s23 ¼

m22 s22 m23

(B7)

Also, from the definition of n, sY 23 ¼

n23 Y s ; n22 22

s23 ¼

n23 s22 n22

(B8)

Em 1  V f ð1  Em =E2f Þ

(A16)

for proportional loading. Comparing Eqs. B2 and B4, while considering Eqs. B7 and B8, leads to

Gm 1  V f ð1  Gm =G12f Þ

(A17)

h23 ¼ h22

n12 ¼ n13 ¼ n12f V f þ nm V m n23 ¼

E2 1 2G23

g23 ¼

(A18)

n22 m22 g22 ¼ g22 n23 m23

(B9) (B10)

(A19) REFERENCES

where G23

Gm ¼ 1  V f ð1  Gm =G23f Þ

(A20)

APPENDIX B: RELATIONSHIP BETWEEN TRANSVERSE ISOTROPIC HARDENING PARAMETERS BASED ON THE PLASTIC WORK EQUIVALENCE PRINCIPLE For proportional loading, the plastic work equivalence becomes dW ¼ siso de þ ade ¼ sij depij

DOI 10.1002/pc

(B1)

1. J.-H. Byun and T.-W. Chou, J. Strain Anal. Eng. Des., 24, 253 (1989). 2. A. Nakai, A. Fujita, A. Yokoyama, and H. Hamada, Composite Structures, 32, 501 (1995). 3. J.S. Tsai, S.J. Li, and L.J. Lee, J. Compos. Mater., 32, 829 (1998). 4. R.E. Powell and H.-Y. Yeh, J. Reinforc. Plast. Compos., 14, 164 (1995). 5. T.D. Kostar, T.-W. Chou, and P. Popper, J. Mater. Sci., 35, 2175 (2000). 6. M.G. Lee, D. Kim, K. Chung, J.R. Youn, and T.J. Kang, Polym. Polym. Compos., 12, 225 (2004). 7. K.L. Branch, K.N. Shivakumar, and V.S. Avva, Threedimensional tow inclination model for calculating elastic constants of three-dimensional triaxial braided composites, AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dy-

POLYMER COMPOSITES—-2007 731

8. 9. 10. 11.

12.

13. 14. 15. 16. 17. 18. 19. 20. 21. 22. 23. 24.

namics, and Materials Conference and Exhibit, 37th, Salt Lake City, UT, USA (Apr. 15–17, 1996), 1788 (1996). A. Aggarwal, S. Ramakrishna, and V.K. Ganesh, J. Compos. Mater., 35, 665 (2001). R.A. Naik, P.G. Ifju, and J.E. Masters, J. Compos. Mater., 28, 656 (1994). N.K. Naik and V.K. Ganesh, Composites, 26, 281 (1995). R.A. Naik, Analysis of 2D triaxial and 3D multi-interlock braided textile composites, AIAA ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference and Exhibit, 37th, Salt Lake City, UT, USA (Apr. 15–17, 1996), Technical Papers. Pt. 3 (1966). R.A. Naik, Multiaxial stiffness and strength analysis of woven and braided composites, AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference and Exhibit, 38th, and AIAA/ASME/AHS Adaptive Structures Forum, Kissimmee, FL, USA (Apr. 7–10, 1997), 1148 (1997). M.S. Dadkhah, J.G. Flintoff, T. Kniveton, and B.N. Cox, Composites, 26, 561 (1995). J.-H. Byun, Compos. Sci. Technol., 60, 705 (2000). Y. Yan and S.V. Hoa, J. Compos. Mater., 36, 963 (2002). T. Ishikawa and T.W. Chou, J. Compos. Mater., 16, 2 (1982). T. Ishikawa and T.W. Chou, J. Mater. Sci., 17, 3211 (1982). T. Ishikawa and T.W. Chou, AIAA J., 21, 1714 (1983). J. Carey, M. Munro, and A. Fahim, J. Reinforc. Plast. Compos., 22, 813 (2003). J. Carey, A. Fahim, and M. Munro, J. Reinforc. Plast. Compos., 23, 1845 (2004). J. Carey, M. Munro, and A. Fahim, Polym. Compos., 26, 152 (2005). I.S. Raju and J.T. Wang, J. Compos. Technol. Res., 16, 15 (1994). Z.M. Huang, Compos. Sci. Technol., 60, 479 (2000). Z.M. Huang, K. Fujihara, and S. Ramakrishna, J. Compos. Technol. Res., 25, 35 (2003).

25. C.J. Redman and C.D. Douglas, ‘‘Advanced materials: Performance through technology insertion,’’ in Proceedings of the 38th International SAMPE Symposium and Exhibition, Anaheim, CA, May 10–13, 719–727 (1993). 26. P.J. Falzon, I. Herszberg, and V.M. Karbhari, J. Compos. Mater., 30, 1210 (1996). 27. B. Gommers, I. Verpoest, and P.V. Houtte, Acta Mater., 46, 2223 (1998). 28. A. Dasgupta, R.K. Agarwal, and S.M. Bhandarkar, Compos. Sci. Technol., 56, 209 (1996). 29. A. Tabiei and I. Ivanov, Int. J. Non Lin. Mech., 39, 175 (2004).

732 POLYMER COMPOSITES—-2007

30. V. Carvelli and C. Poggi, Compos. A, 32, 1425 (2001). 31. T. Zeng, L.-Z. Wu, and L.-C. Guo, Compos. Struct., 64, 399 (2004). 32. D. Bigaud and P. Hamelin, Compos. Struct., 38, 361 (1997). 33. A.E. Bogdanovich and C.M. Pastore, Compos. Sci. Technol., 56, 291 (1996). 34. D. Goyal, X. Tang, J.D. Whitcomb, and A.D. Kelkar, Mech. Adv. Mater. Struct., 12, 113 (2005). 35. I. Ivanov and A. Tabiei, Compos. Struct., 54, 489 (2001). 36. S.C. Quek, A.M. Waas, K.W. Shahwan, and V. Agaram, International Journal of Mechanical Sciences, 45, 1077 (2003). 37. P. Tan, L. Tong, and G.P. Steven. Composite Structures, 47, 797 (1999). 38. Z. Hashin and B.W. Rosen, J. Appl. Mech., 31, 223 (1964). 39. C.C. Chamis, NASA Technical Memorandum 83320, Presented at the 38th Annual Conference of the Society of Plastics Industry, Houston, TX, February 7–11 (1983). 40. R. Hill, Proc. R. Soc. London Ser. A, 193, 281 (1948). 41. J.L. Chaboche, Int. J. Plast., 2, 149 (1986). 42. E. Voce, J. Inst. Metals, 74, 537 (1948). 43. Hibbitt, Karlsson, and Sorensen Inc., ABAQUS User’s Manual, Version 6.5, Hibbitt, Karlsson, and Sorensen Inc., Providence, RI (2003). 44. A.E. Bogdanovich and C.M. Pastore, Mechanics of Textile and Laminated Composites. Chapman and Hall, London, UK (1996). 45. J. Aboudi, Int. J. Eng. Sci., 20, 605 (1982). 46. C.T. Sun and R.S. Vaidya, Compos. Sci. Technol., 56, 171 (1996). 47. M. Paley and J. Aboudi, Mech. Mater., 14, 127 (1992). 48. D. Bigaud and P. Hamelin, ‘‘Modeling of 2D and 3D textile-reinforced composites by means of imbricate-type elements approaches,’’ in Computational Techniques for Materials, Civil-Comp Press, Edinburgh, UK, 73 (2000). 49. T.O. Williams and J. Aboudi, Acta Mech., 138, 151 (1999). 50. C.M. Pastore, Compos. Manuf., 4, 217 (1993). 51. J.E. Masters, R.L. Foye, C.M. Pastore, and Y.A. Gowayed, J. Compos. Technol. Res., 15, 112 (1993). 52. W. Zhao, L. Minnetyan, and C.C. Chamis, ‘‘Design of smart composite structures for durability,’’ in Proceedings of the SAMPE 2003 Symposium and Exhibition, Long Beach, CA, May 11–15 (2003). 53. J. Singletary, H. Davis, W. Knoff, and M.K. Ramasubramanian, ‘‘Transverse compression of PPTA fibers,’’ in New Frontiers in Fiber Science (2001). Abstract. The Fiber Society Spring Meeting 2001, Raleigh, NC, May 23–25 (2001). 54. H.H. Yang, Kevlar Aramid Fiber, Wiley, New York (1993).

DOI 10.1002/pc