Composites: Part B 33 (2002) 45±56

www.elsevier.com/locate/compositesb

A dual homogenization and ®nite element approach for material

characterization of textile composites
Xiongqi Peng, Jian Cao*
Department of Mechanical Engineering, Northwestern University, 2145 Sheridan Road, Evanston, IL 60208, USA
Received 11 February 2001; accepted 25 June 2001

Abstract
A novel procedure for predicting the effective nonlinear elastic moduli of textile composites through a combined approach of the
homogenization method and the ®nite element formulation is presented. The homogenization method is ®rst applied to investigate the
meso-microscopic material behavior of a single ®ber yarn based on the properties of the constituent phases. The obtained results are
compared to existing analytical and experimental results to validate the homogenization method. Very good agreements have been obtained.
A unit cell is then built to enclose the characteristic periodic pattern in the textile composites. Various numerical tests such as uni-axial and
bi-axial extension and trellising tests are performed by 3D ®nite element analysis on the unit cell. Characteristic behaviors of force versus
displacement are obtained. Meanwhile, trial mechanical elastic constants are imposed on a four-node shell element with the same outer size
as the unit cell to match the force±displacement curves. The effective nonlinear mechanical stiffness tensor is thus obtained numerically as
functions of elemental strains. The procedure is exempli®ed on a plain weave glass composite and is validated by comparing to experimental
data. Using the proposed approach, the nonlinear behavior of textile composites can be anticipated accurately and ef®ciently. q 2002
Elsevier Science Ltd. All rights reserved.
Keywords: A. Fiber; A. Hybrid; C. Finite element analysis; Homogenization method

1. Introduction Among the manufacturing processes for composite

materials, stamping is a potential one to solve the presently
Textile composite materials have recently received existing bottle-neck in the mass production of composites.
considerable attention, due to their structural advantages Stamping is a high-speed, highly automated manufacturing
of high speci®c-strength and high speci®c-stiffness as well process that has been used extensively in the high-volume
as their improved resistance to impact. Compared with automotive industry, primarily for sheet metal forming.
unidirectional composites, the interlacing of ®ber bundles Recently, researchers and engineers are trying to extend it
in textile composites prevents the growth of damage and to composite materials in order to achieve a low-cost mass
hence provides an increase in impact toughness. Besides production of textile composites [1±6]. To fully understand
their advantageous mechanical properties, textile compo- the mechanical behavior of textile composites during
sites are easy to handle and have excellent formability and processing and for the purpose of optimal design, it is essen-
hence are widely employed in aircraft, boat and defense tial to obtain the effective material properties, such as elastic
industries. However, they remain at the periphery of high constants, thermal conductivity and coef®cient of viscosity
volume industries such as construction, automotive, and of fabric composites from the known material properties of
consumer goods because of long process cycle time. The the constituent phases.
lack of an effective and economic manufacturing process for In recent years, many efforts have been given to the
textile composites seriously limits the usage of this kind of estimation of effective material properties of composite
materials. Hence, extensive simulation and experimental materials. Boisse et al. [2±4] developed a constitutive
research on the development of a new fabrication technol- model for a single thread in tension and then built a
ogy for cost-effective, mass production of composite membrane model based on potential energy minimization
materials is needed. to simulate glass fabric shaping. Chen and Cheng [7]
de®ned an exponential function to represent the orientation
* Corresponding author. Tel.: 11-847-467-1032; fax: 11-847-491-3915. distribution of ®bers in a mis-oriented short-®ber composite.
E-mail address: [email protected] (J. Cao). They then predicted the effective elastic moduli of the
1359-8368/02/$ - see front matter q 2002 Elsevier Science Ltd. All rights reserved.
PII: S 1359-836 8(01)00052-X
46 X. Peng, J. Cao / Composites: Part B 33 (2002) 45±56

composite by combining the mean stress concept of [14] simulated the deep drawing process for thermoplastic
Mori and Tanaka and the eigenstrain idea of Eshelby. composite laminates by including ¯ow, heat transfer and
By applying the micro-mechanical techniques on a sub- residual stress analysis. They obtained the governing equa-
volume of a unit cell of textile composites, Weissenbach tions and material properties of the composites at the form-
et al. [8] provided basic data to re¯ect the variations in ing temperature by the homogenization method. Yi et al.
the elastic stiffness within the composite. The data was [15] obtained the effective viscoelastic moduli of a compo-
interpolated and then a stiffness function was obtained site with a periodic microstructure in the time and frequency
by curve-®tting. They ®nally obtained an averaged stiff- domain. A summarization of the hierarchical homogeniza-
ness by mapping the stiffness function over the whole tion models for composite material characterization was
part. Gowayed and Yi [9] constructed a model to predict given by Takano et al. [16].
the inelastic behavior, damage progression and failure of Though the homogenization method is effective in pre-
textile reinforced composites by integrating with a 3D dicting material properties, the huge computational cost
hybrid ®nite element approach. They applied an expanded limits its application in simulating the forming of complex
Hahn and Tsai model to characterize the inelastic behavior structures as textile composites. By integrating the advant-
of the textile composites and took the maximum stress ages of the conventional ®nite element formulation and the
criterion to predict the initial failure. Wang and Sun homogenization method, this paper presents a novel proce-
et al. [10,11] investigated the deformation behavior of dure for predicting the effective nonlinear elastic moduli of
PEEK unidirectional thermoplastic composite at ele- textile composites. These nonlinear effective moduli can be
vated temperature during forming. They regarded the incorporated in a user material subroutine associated with
composite ply as an orthotropic and homogeneous con- shell elements in commercial FEM codes so that an ef®cient
tinuum. An anisotropic elastic/viscoplastic constitutive and accurate FEM simulation of composite sheet forming is
model was employed to describe the inelastic and rate- feasible. First, based on the properties of the constituent
dependent material response. phases, the homogenization method is employed to predict
Due to the immense variety of available composite the effective elastic constants of the ®ber yarn of the textile
materials, possible fabric construction geometry, and the composites, which is regarded as a unidirectional compo-
change of composite material properties with processing site. A unit cell is then built to enclose the characteristic
temperatures, it is very time-consuming to obtain material periodic pattern in the textile composites. Using the unit
characterizations of various composites by experimental cell, various numerical tests can be performed. By correlat-
approach. Analytical methods, on the other hand, cannot ing the force versus displacement curves of the unit cell and
deal with complex fabric architecture. A potential approach a four-node shell element with the same outer size as the
to obtain the effective material properties of composites unit cell, the effective nonlinear mechanical stiffness tensor
ef®ciently is the homogenization method. can be obtained numerically as functions of strain tensor.
As a systematic mathematical approach, the homogeniza- The procedure is exempli®ed on a plain weave glass com-
tion method has several attractive features. It is capable of posite and is validated by comparing with experimental
dealing with complex microstructures in a periodic order, results. The macroscopic response in terms of the effective
which may be too computationally intensive to be discre- material moduli of textile composites will be the inputs to
tized into a conventional FEM model. Furthermore, the our future stamping simulation. The entire approach is
conventional numerical tools such as FEM can be easily shown in Fig. 1.
incorporated with this method with small modi®cations.
Most importantly, the homogenization method guarantees
convergence, which implies that the solution of a periodic
microstructure will not be signi®cantly affected by the
speci®c boundary conditions imposed on the global system.
The homogenization method has been successfully applied
in the estimation of the effective material properties of
composites and the simulation of the manufacturing process
of composites. Dasgupta and Agarwa [12] predicted the
orthotropic thermal conductivity of plain-weave ®ber re-
inforced composites by the homogenization method. They
built a unit cell to enclose the characteristic periodic pattern
in the composite. A newly developed 3D series±parallel
thermal resistance network was then used to solve the
steady-state heat transfer problem for this unit cell. Liew
et al. [13] formulated a periodic model to estimate the effect
of porosity on the moduli of woven fabric composites and
validated their model by experiments. Hsiao and Kikuchi Fig. 1. Multi-scale material characterization approach.
 X. Peng, J. Cao / Composites: Part B 33 (2002) 45±56 47

and
! !
1 2vi 2vj 1 2vi 2vj
eijx v† ˆ 1 ; eijy v† ˆ 1
2 2xj 2xi 2 2yj 2yi
8†
The elastic coef®cients Cijkl are periodic functions of x and
depend on 1 :
1
Cijkl ˆ Cijkl x=1† ; Cijkl y† 9†
Thus the stress can be expressed as:
1 21
s ij1 ˆ s x; y† 1 s ij0 x; y† 1 1s ij1 x; y† 1 ¼ 10†
1 ij
where
1 n
s ijn x; y† ˆ Cijkl ekl x; y†; n ˆ 21; 0; 1; ¼ 11†

Fig. 2. Global structure and meso-microstructure. Therefore, the elastic problem with a periodic microstruc-
ture as shown in Fig. 2 is presented as:
2. Basic equations of the homogenization method for s ij;1 j 1 fi ˆ 0 in V 12†
elastic cases
s ij1 nj ˆ Ti on G 1 13†
We assume that an elastic body as shown in Fig. 2 is an
assembly of a periodic microscopic unit cell. The global u1i ˆ 0 on G 2 14†
coordinate x is related to the local coordinate y as:
Substituting Eq. (11) in (12), and equaling the powers of 1;
y ˆ x=1 1† we obtain:

where 1 is the ratio between the dimension of a unit cell and 2s ij21
ˆ0 15†
a structure body, which is a very small positive number. 2yj
The important and essential postulation in the homogeni-
zation method is that the displacement ®eld ui can be 2s ij0
expressed as an asymptotic expansion: ˆ0 16†
2yj
ui ˆ u1i x† ˆ u0i x; y† 1 1u1i x; y† 1 12 u2i x; y† 1 ¼;
2† 2s ij0 2s ij1
y ˆ x=1 1 1 fi ˆ 0 17†
2xj 2yj
Noting that, The variational form for Eq. (15) is:
1 !
2f x; y† 2f x; y† 1 2f x; y† Z 2s 21 Z 0
ˆ 1 3† ij 0 1 2uk
xi 2xi 1 2yi dui dV ˆ Cijkl du0 dV ˆ 0 18†
V 1 2yj V1 2yl ;j i
where f is a general function. We have, for the strain tensor For a Y-periodic function F y†; we de®ne a mean operator
eij : as follows:
!  
1 2ui 2uj 1 Z x 1 Z Z
eij ˆ 1 ˆ e21 x; y† lim1 F d V ˆ F y†dY dV 19†
2 2xj 2xi 1 ij 1 !0 V1 1 uYu V Y
where uYu stands for the volume of the unit cell.
1e0ij x; y† 1 1e1ij x; y† 1 ¼ (4)
Since the homogenization method consists of ®nding the
limit of solutions to Eqs. (15)±(17) as 1 tends to zero, we
where
have for Eq. (18):
!
e21 0
ij x; y† ˆ eijy u † 5† Z 0
1 2uk
lim Cijkl du0 dV
1 !0 1 V 1 yl ;j i
e0ij x; y† ˆ eijx u0 † 1 eijy u1 † 6† !
1 Z Z 2u0
ˆ C ijkl k du0i dY dV ˆ 0 20†
e1ij x; y† ˆ eijx u1 † 1 eijy u2 † 7† uYu V Y yl ;j
48 X. Peng, J. Cao / Composites: Part B 33 (2002) 45±56

Using the divergence theorem on Eq. (20) yields: Substituting Eq. (28) in (27) yields
!
1 Z Z 2u0k Z Z 2u1k 2du1i
C ijkl du0 dY dV C ijkl dY dV
uYu V Y yl ;j i V Y 2yl 2yj
Z 2u 0 Z 2C pkl 2du1i
1 Z I 2u0 k
ˆ C ijkl k nj du0i ds dV ˆ 0 21† 1 Cijpq dY dV
uYu V 2yl V 2xl Y 2yq 2yj
s
ˆ0 29†
Thus
Applying the divergence theorem leads to
2u k0
ˆ0 22† Z I
2yj 2du1i
C ijpq u1p nq ds dV
V 2yj
i.e. u0i is a function of x only, and Eq. (2) can be rewritten as: s

ui ˆ u0i x† 1 1u1i x; y† 1 12 u2i x; y† 1 ¼; y ˆ x=1 23† Z I 2u0k 2du1i

1 C ijpq C pkl nq ds dV
We can regard u0i as the macroscopic displacement, while V 2xl 2yj
s
u1i ; u2i ; ¼ are the microscopic displacements. The physical
meaning of Eq. (23) thus is that the real displacement ui ˆ0 30†
is rapidly oscillating around u0i due to the inhomogenity
from the microscopic point of view. u1i ; u2i ; ¼ are the Hence we get
perturbing displacements according to the microstructure. 2u0k
Hence, the scale factor 1 is multiplied to the microscopic u i1 ˆ 2C ikl x; y† 31†
2xl
displacements.
From Eqs. (6), (8) and (11), we have Substituting Eq. (31) in (24) yields:
! !
2u0k 2u1k 0 2C mkl 2u0k
0
s ij ˆ Cijkl 1
1 24† s ij ˆ Cijkl 2 Cijmn 32†
2xl 2yl 2yn 2xl

Substituting Eq. (24) into the variational form of Eq. (16) Then integrating on Y leads to the stress±strain relations for
yields: the homogenized elastic body:
!
Z 2u0k 2u1k 2u0k
1
Cijkl 1 du1 dV ˆ 0 25† s~ ij0 ˆ Cijkl
H
33†
V1 2xl 2yl ;j i 2xl
where
As 1 ! 0 1 ; we have,
! 1 Z 0
Z 2u0k 2u1 s~ ij0 ˆ s ij x; y†dY 34†
lim1 1
Cijkl 1 k du1i dV uYu Y
1!0 V 1 2xl 2yl ;j " #
! 1 Z 2C mkl
1 Z Z
H
2u0k 2u 1 Cijkl ˆ Cijkl y† 2 Cijmn y† dY 35†
ˆ C ijkl 1 k du1i dY dV ˆ 0 26† uYu Y 2yn
uYu V Y 2xl 2yl ;j
H
Cijkl denotes the homogenized elastic constants.
Integrating by parts, and noting that ˆ 0 at the bound- du i1 Eq. (35) can be obtained numerically using the ®nite
ary of Y, and u0i is a function of x only, we obtain: element discretization with speci®c boundary conditions
Z Z 2u1k 2du1i on a unit cell. According to the standard conventions in
C ijkl dY dV the ®nite element analysis, the microscopic Eq. (28) can
V Y 2yl 2yj
be written as follows:
! Z 
Z 2u 0 Z 2du1i Z
1 k
Cijkl dY dV BT CB dY C ij ˆ BT Cij dY 36†
V 2xl Y 2yj Y Y

where C is the stress±strain matrix of the constituents of the

ˆ0 27†
composite, B is the discrete displacement±strain matrix
Introducing C which satis®es depending on the elemental shape functions. C ij is a vector
of column of ij ij ˆ 11; 22; 33; 23; 31; 12† of the stress±
Z 2C pkl 2du1i Z 2du1i
C ijpq dY ˆ Cijkl dY 28† strain matrix C. C ij is the characteristic displacement
Y 2yq 2yj Y 2yj vector associated with the ij mode. A conventional ®nite
 X. Peng, J. Cao / Composites: Part B 33 (2002) 45±56 49

element solver such as ABAQU/Standard can be used to Table 1

calculate Eq. (36) by imposing appropriate multi-point Material properties of E-glass/PP composite
constrains to satisfy the periodic boundary conditions of Property Unit E-Glass Epoxy Polypropylene
the unit cell.
The homogenized elastic constants de®ned by Eq. (35), Tensile modulus GPa 73.1 3.45 1±1.4
Poisson's ratio ± 0.22 0.35 0.3
therefore, can be expressed as follows:
Axial shear modulus GPa 30.19 1.83
1 Z Density kg/m 3 2540 900
‰CH Š ˆ C I 2 C†dY 37†
uYu Y
where The in-plane shear modulus G12 ; in the mechanics approach,
is given as:
C ˆ C 11 ; C 22 ; C 33 ; C 32 ; C 31 ; C 12 † 38†
1 V V
ˆ f 1 m 42†
G12 Gf Gm
3. Elastic constants of unidirectional ®ber reinforced The Halpin±Tsai equations for the longitudinal elastic
composites modulus and major Poisson's ratio are the same as those
in the mechanics approach. While the transverse elastic
A detailed ®nite element model (FEM) was developed to modulus E2 given by the H±T equation is:
model a unit cell of the unidirectional E-glass/epoxy
E2 1 1 jhV f
composite (Fig. 3). The ®ber bundle is shown in gray with ˆ 43†
its cross-section represented by a circle. The matrix is Em 1 2 h Vf
shown in darker gray in Fig. 3. The elastic constants for where j is the reinforcing factor depending on ®ber geo-
the glass ®ber and epoxy matrix are given in Table 1. metry, packing geometry and loading conditions. For a
The homogenization method is then employed to deter- circular ®ber in square array packing geometry, j has the
mine the elastic constants of the present composite. The value of 2. h is given as:
obtained elastic constants under various ®ber volume frac-
E f =Em † 2 1
tion Vf are compared with those from the mechanics of hˆ 44†
materials approach [17] and Halpin±Tsai's equation [18] Ef =Em † 1 j
as well as experimental data [19,20]. The H±T equations for the in-plane shear modulus G 12 is
The longitudinal elastic modulus E1 given by the
G12 1 1 jhV f
mechanics approach in Ref. [17] is: ˆ 45†
Gm 1 2 hVf
E1 ˆ Ef Vf 1 E m Vm 39†
where
where Vf is the ®ber volume fraction, Vm is the matrix
Gf =Gm † 2 1
volume fraction. hˆ 46†
The transverse elastic modulus E2 is given as: Gf =Gm † 1 j

1 V V Fig. 4 shows the effective longitudinal elastic modulus

ˆ f 1 m 40†
E2 Ef Em
The major Poisson's ratio n12 is
n12 ˆ nf Vf 1 nm Vm 41†

Fig. 3. FEM model of the unit cell for ®ber yarn. Fig. 4. Axial elastic modulus of the composite.
50 X. Peng, J. Cao / Composites: Part B 33 (2002) 45±56

Fig. 5. Transverse elastic modulus of the composite. Fig. 7. Axial Poisson's ratio.

of the composite E 1 with different ®ber volume fractions. attributed to the imperfect adhesion between ®bers and
The longitudinal elastic modulus E1 obtained from the matrices, the possible separation of ®bers in one bundle
mechanics approach is usually in a good agreement with under compression and shear, and some possible voids in
experimental data. As shown in Fig. 4, E1 values calculated the composite.
from the homogenization method are very close with those Fig. 6 shows the in-plane shear modulus of the composite
from the mechanics approach. as a function of ®ber volume fraction. Again, the mechanics
In contrast, as shown in Fig. 5, the experimental data approach underestimates the in-plane shear modulus. Both
for the transverse elastic modulus [19] are quite higher the homogenization method and H±T equation are in good
compared to those obtained from the mechanics approach, agreements with the experimental data [20].
and are slightly lower than those predicted by the Fig. 7 shows the axial Poisson's ratio obtained from the
homogenization method. The H±T equation, on the other homogenization method and the mechanics approach. As
hand, provides a good approximation due to the fact that shown in Fig. 7, the homogenization method provides a
the H±T equation is a semi-analytical empirical func- lower axial Poisson's ratio. This is partially due to the
tion. The mechanics approach provides a low bound for higher transverse elastic modulus estimated by this method.
the transverse elastic modulus. The difference between the No experimental data of axial Poison's ratio is available for
experimental data and the homogenization method can be comparison.

4. Effective mechanical stiffness tensor for textile

composites

The homogenization method presented above demon-

strates its accuracy and reliability in obtaining the equiva-
lent elastic moduli of a unidirectional composite. In theory,
the method can be used in the consequent simulation and
optimization of fabrication processes of textile composites,
which are made of the unidirectional ®ber bundles.
However, the computation cost will be too high for real-
life engineering. Here, we propose an equivalent homo-
genization method, which would be applicable to any
complicated ®ber architecture. In the following, we will
demonstrate the approach of determining the equivalent
nonlinear material properties of a plain weave E-glass/PP
composite.
The objective is to develop a means that can characterize
Fig. 6. In-plane shear modulus. the global material behaviors of textile composites using a
 X. Peng, J. Cao / Composites: Part B 33 (2002) 45±56 51

cost-ef®cient shell element as in the sheet metal forming presented by McBride and Chen [22] is used in this paper
simulation. However, the behavior of textile composites to model a unit cell for the plain weave composite. The unit
during the shaping process is very different from that of a cell, as shown in Fig. 8, is de®ned by four sinusoidal curves
sheet metal. During the sheet-metal forming process, the in terms of yarn width w, yarn spacing s, and fabric thick-
blank is usually subjected to large extension strains. The ness h to represent the periodic pattern in the composite. The
fabric yarns, on the contrary, undergo small extension cross-sections of ®ber yarns are approximated as circular
along the yarn directions while experiencing large angular arcs. The characteristic values of w and s are averaged
deformation between weft and warp yarns (a phenomenon from 30 measurements as w ˆ 3:72 mm; s ˆ 5:14 mm:
sometimes called trellising) [2,3]. Consequently, the model The thickness h is obtained by solving the following implicit
used for simulating the fabric behavior will account for equation:
large deformation and consider geometric nonlinearity. m
The material constitutive equation is given as the following Vunit-cell ˆ 4AL ˆ unit-cell 49†
rcomposite
with the assumption of an orthotropic material at the plane
stress condition: where Vunit-cell is the actual occupied volume of the unit cell,
2 3 A is the cross-section area of the ®ber yarn, L is the undu-
8 9 E1 n12 E2 8 9
s1 > 6 1 2 n12 n21 1 2 n12 n21 0 7> e1 > lated length of any yarn in the unit cell, munit-cell is the
>
> > 6 7 > >
< = 6 7< = averaged unit cell mass, rcomposite is the density of the
s2 ˆ 6 6 n E
12 2 E 2
7
7 e 2 composite. A, L and rcomposite are given as:
> > 0 7> >
>
: ; 6
> 4 1 2 n12 n21 1 2 n12 n21 5>: >
; !2  
s
12 g 12 h2 1 w2 2wh w w 2 2 h2
21
0 0 G12 Aˆ sin 2 50†
4h h2 1 w 2 2 4h
47†
where h2 p2
Lˆs1 51†
E1 E 16s
ˆ 2 48†
n12 n21
rcomposite ˆ vf rfiber 1 vm rmatrix 52†
To account for large deformation, the logarithmic strain is
employed for the shell elements. where vf is the volume fraction of the ®ber and vm is the
The nonlinear elastic constants in Eq. (47) will be volume fraction of the matrix. For our plain weave compo-
obtained by matching the force±displacement curves of site, the calculated thickness h ˆ 0:39 mm:
numerical tests, such as uni-axial, bi-axial extensions and The material properties of the constituent phases of the
trellising [21], between a detailed unit cell model in Fig. 8 ®ber yarn (E-glass ®ber) and resin (Polypropylene), are
and a four-node shell element with the same outer size of the listed in Table 1. The volume fraction of the E-glass is
unit cell. 70%. The ®ber yarns of the plain weave composites are
regarded as unidirectional composites. The elastic moduli
4.1. Unit cell for plain weave composites are assumed to be linear and orthotropic, and are obtained
from the above homogenization method. The predicted
The geometric description for plain weave composites elastic constants for the ®ber yarns are:
El ˆ 51:92 GPa; Et ˆ 21:97 GPa; nlt ˆ 0:2489

ntt ˆ 0:2143; Glt ˆ 8:856 GPa; Gtt ˆ 6:250 GPa

where subscript l represents the longitudinal direction and t
denotes the transverse direction of the ®ber yarn. Notice that
the yarn crushing due to the compression between ®ber
yarns is not considered here.
The unit cell is discretized by 3D 8-node continuum
elements. Each ®ber yarn is modeled by 304 elements.
The above effective elastic constants obtained from the
homogenization method presented in Section 2 are imposed
on the ®ber yarns. The pin-jointed net idealization is
assumed along the four corner lines of the unit cell. Contact
conditions are prescribed between the possible interlacing
surfaces of the ®ber yarns under loading. Unless speci®ed,
the friction coef®cient in the surface contacting is assumed
Fig. 8. Unit cell for plain weave composites. to be 0.05 as from Ref. [21]. Boundary conditions will be
52 X. Peng, J. Cao / Composites: Part B 33 (2002) 45±56

Fig. 9. Equivalent shell element in extension test.

modi®ed to re¯ect the periodic boundary conditions of the

macrostructure under different loading conditions.

4.2. Numerical tests on plain weave unit cell

Numerical tests as uni-axial, bi-axial extensions and trel- Fig. 11. Shrinkage in uni-axial extension.
lising [21] are conducted on the unit cell for the determina-
tion of the equivalent elastic constants for the shell element. the ratio between the warp direction major strain e1 and the
The nonlinear mechanical stiffness tensor is iteratively weft direction major strain e2 at the end of each increment as
modi®ed and imposed on the shell element in each incre- the value of Poisson's ratio at the strain level of e1 : The thus
ment of the FEM analysis to match the force±displacement obtained Poisson's ratio n12 is plotted in Fig. 13 as circles.
curves of the unit cell obtained from numerical tests, respec- Finally the equation for Poisson's ratio n12 is obtained via
tively, as shown in Figs. 9 and 10. The shell element will be curve ®tting as:
prescribed by the same extension displacement in each (
0:02 e1 , 0
increment as the unit-cell. Besides, for the uni-axial exten- n12 ˆ 2
sion case, a same amount of resulting shrinkage in the trans- 346:634e1 2 25:9684e1 1 0:54416 0 # e1 , 0:04
verse direction as that in the unit-cell will also be imposed 53†
on the shell element.
The unit cell is ®rst subjected to uni-biaxial and equal bi- The tensile elastic moduli E1 and E2 are determined by the
axial extension tests for the determination of Poisson's ratio combination of the uni-axial and equal bi-axial extension
and the tensile elastic moduli. In the uni-axial extension test, tests. The obtained curve of extension displacements versus
the resulting shrinkages along the weft direction are shown loading forces in equal bi-axial extension test is shown as
in Fig. 11 with circles. The reaction forces due to the exten- circles in Fig. 14. Since in extension tests, the shear strain is
sions along the ®ber warp direction are shown in Fig. 12 zero and there is no coupling between extensions and shear-
with circles. ing, an arbitrary small positive constant is assigned to the
The equivalent Poisson's ratio n12 for the shell element is
obtained solely from the uni-axial extension test. We take

Fig. 10. Equivalent shell element in trellising. Fig. 12. Reaction force in uni-axial extension tests.
 X. Peng, J. Cao / Composites: Part B 33 (2002) 45±56 53

Fig. 15. Tensile elastic modulus in extension tests.

Fig. 13. Equivalent Poisson's ratio for shell element.
symmetry in the equal bi-axial extension test. E1 for the
shear modulus in the mechanical stiffness tensor for the equal bi-axial test is shown in Fig. 15 with a solid line,
shell element. Furthermore, in the uni-axial extension test, which can be described using E1 ˆ g e1 †:
the strain in the weft (shrinkage) direction e2 is negative. From Fig. 15, it is observed that the tensile elastic
Experimental studies [2,3] in textile composites showed that modulus will increase with the increment of difference
the yarn buckles immediately under a compression load. between the warp direction major strain e1 and the weft
Hence, the fabric compressive stiffness is taken to be negli- direction major strain e2 : Hence we choose the equation
gible. We assume that the tensile elastic modulus E2 in the g e1 † in the equal bi-axial extension test as a foundation
uni-axial extension test is very small and is taken to be 0.1. and add a correction to include the effect of other strain
The tensile elastic modulus E1 is then numerically deter- distributions. The ®nal form for the tensile elastic modulus
mined by matching the force±displacement curve of unit obtained is:
cell as in Fig. 12 with circles to those of the corresponding 8
single shell element. The resulting E1 is shown in Fig. 15 >
> 10g212 1 0:1 ei , 0
<
with a dotted line, which can be represented by E1 ˆ f e1 †: 3 2
Ei MPa† ˆ A3 ei 1 A2 ei 1 A1 ei 1 A0 1 0 # ei , 0:04
Secondly, the tensile elastic modulus E1 for the equal bi- >
>
:
axial extension test is numerically determined by the similar ue1 2 e2 u exp B1 ei 1 B0 †
approach as in the uni-axial extension test, except that E1 54†
and E2 are equal under the same strain level due to the
where i ˆ 1; 2 and

A3 ˆ 7:60033 £ 107 ; A2 ˆ 29:81521 £ 106 ;

A1 ˆ 4:35522 £ 105 ; A0 ˆ 4:64523 £ 103 ;

B1 ˆ 287:34587; B0 ˆ 13:02737

The equivalent shear modulus G12 is determined from

the trellising test. In the trellising test, the edges of the
shell element will be clamped, but can rotate freely along
the corner points. This is to consist with the pin-jointed net
idealization of the unit-cell. The deformed shape of the unit-
cell is shown in Fig. 16. The force±displacement curve of
the unit-cell in the trellising is shown with circles in Fig. 17.
In the trellising, the major strains are negligible comparing
with the shear strain and there is no coupling between exten-
sion and shearing, we hence assign small positive constants
to the tensile elastic moduli, E1 and E2 ; and Poisson's ratio
n12 : Again, the shear modulus G12 is numerically deter-
Fig. 14. Reaction force in equal bi-axial extension. mined by matching the force±displacement curve of the
54 X. Peng, J. Cao / Composites: Part B 33 (2002) 45±56

Fig. 18. Equivalent shear modulus G12 for shell element.

Fig. 16. Deformed unit cell in trellising.

unit cell as in Fig. 17 with circles to those of the correspond- 5. Validation of material model
ing single shell element. The obtained G12 is shown in Fig.
18 with circles and can be represented as: The previous obtained material model is validated on two
( scales. First, the user material subroutine is integrated into
1:7444g12 1 4:96 £ 1022 g12 # 0:207
G12 MPa† ˆ the single shell element used in the material characteriza-
3 2
C3 g12 1 C2 g12 1 C1 g12 1 C0 g12 . 0:207 tion. The uni-axial and equal bi-axial extension tests and the
55† trellising tests are repeated by using this new material
model. The obtained load±displacement curves of the
where
shell element under different numerical tests are compared
C3 ˆ 11:3103; C2 ˆ 22:3030; with those corresponding to the unit-cell in Figs. 12, 14
and 17, respectively, with solid lines. As shown in those
C1 ˆ 28:9429; C0 ˆ 1:2059 ®gures, the single shell element with the new material
model can provide very good predictions of reaction forces
Finally, the above representations of the Poisson's ratio n12 in the unit-cell under extension and trellising tests. The
(Eq. (53)), the tensile elastic moduli E1 and E2 (Eq. (54)), shrinkage of the shell element in the uni-axial extension
and the shear modulus G12 (Eq. (55)) as well as the consti- test is shown in Fig. 11 with a solid line. A good agree-
tutive Eq. (47) are implemented in a user material sub- ment is obtained by comparing with the shrinkage in the
routine and integrated to the ABAQUS/Standard input ®le unit-cell.
for further simulation.

Fig. 17. Comparison of stretch force in trellising. Fig. 19. Comparison of Fx in unequal bi-axial extension.
 X. Peng, J. Cao / Composites: Part B 33 (2002) 45±56 55

Fig. 22. Comparison with experiment results.

Fig. 20. Comparison of Fy in unequal bi-axial extension. versus displacement curves along the warp and directions
are drawn with solid lines in Figs. 19 and 20, respectively.
As shown in Figs. 19 and 20, good agreements are obtained
To further validate the material model, an unequal bi- by comparing with the unit cell results.
axial extension test is conducted on the unit cell. The exten- Secondly, multi-element FEM model is used to simulate
sions along the ®ber warp direction ewarp and the weft the experimental trellising tests conducted on different areas
direction eweft are chosen to be ewarp ˆ 2eweft : The obtained of composite patches [22]. Here, we have three cases: 0.012,
force versus displacement curve along the ®ber warp 0.025 and 0.044 m 2. The composite patch is clamped on
direction is shown in Fig. 19 with circles. The other four sides in the trellis shear frame with the warp and
force±displacement curve along the weft direction is weft tows parallel to the side of the frame, as shown in
shown with circles in Fig. 20. Again, the shell element Fig. 21. Hinges at each corner allow the sides of the
with the above user material subroutine is used to simulate frame to rotate freely. Fig. 22 shows the load normalized
this unequal bi-axial extension. The corresponding force by the composite patch area with respect to the shear angle
calculated from the cross-corner displacement at the room
temperature. FEM simulation is implemented on a compo-
site patch with a size of 24 cm £ 24 cm: The composite
patch is discretized by 10 £ 10 identical shell elements.
The boundary conditions are prescribed to re¯ect the beha-
vior of the frame, as stated before. Again, the user material
subroutine is imposed on each shell element. The resulting
normalized load±shear angle curve is shown in Fig. 22. As
can be seen from Fig. 22, the FEM simulation results are
very close to the experimental data. This validates the effec-
tiveness of our procedure for predicting the mechanical
stiffness tensor for plain weave composites.

6. Summary

A novel procedure is presented in this paper for predict-

ing the effective nonlinear elastic moduli of textile compo-
sites by using the ®nite element formulation and the
homogenization method. The procedure is exempli®ed on
a plain weave composite. The general procedure can be
summarized as:

Step 1. build a unit cell for the ®ber yarn.

Fig. 21. Trellis shear frame in experiment. Step 2. apply the homogenization method illustrated in
56 X. Peng, J. Cao / Composites: Part B 33 (2002) 45±56

Section 2 to obtain the elastic constants of the ®ber yarn. investigation of some press forming parameters of two ®bre re-
Step 3. build a unit cell for the textile composite. inforced thermoplastics: APC2-AS4 and PEI-CETEX. Composite
Part A 1998;29:101±10.
Step 4. characterize the material behavior of force versus [7] Chen CH, Cheng CH. Effective elastic moduli of misoriented short-
displacement by numerical trellising and extension tests ®ber composites. Int J Solids Struct 1996;33(17):2519±39.
on the unit cell via FEM analysis. [8] Weissenbach G, Limmer L, Brown D. Representation of local stiff-
Step 5. obtain the equivalent elastic constants for shell ness variation in textile composites. Polym Polym Compos 1997;5(2):
elements by correlating the force versus displacement 95±101.
[9] Gowayed Y, Yi L. Mechanical behavior of textile composite materials
curves of the unit cell and a shell element with the
using a hybrid ®nite element approach. Polym Compos 1997;18(3):
same outer size as the unit cell. 313±9.
[10] Wang C, Sun CT, Gates TS. Elastic/viscoplastic behavior of ®ber-
Comparison with experimental results validates the effec- reinforced thermolplastic composites. J Reinforc Plast Compos 1996;
tiveness of this procedure. The procedure can be easily 15:360±77.
extended to other textile composites by building the corre- [11] Wang C, Sun CT. Experimental characterization of constitutive
models for PEEK thermoplastic composite at heating stage during
sponding unit cell. Future work will include the implemen-
forming. J Compos Mater 1997;31:1480±506.
tation of the material model in a general stamping [12] Dasgupta A, Agarwal RK. Orthotropic thermal conductivity of plain-
simulation and process optimization. weave fabric composite using a homogenization technique. J Compos
Mater 1992;26(18):2736±58.
[13] Liaw PK, Yu N, Hsu DK, Miriyala N, Saini V, Snead LL, McHargue
Acknowledgements CJ, Lowden RA. Moduli determination of continuous ®ber ceramic
composites (CFCCs). J Nucl Mater 1995;219:93±100.
The authors would like to thank the NSF Division of [14] Hsiao SW, Kikuchi N. Numerical analysis of deep drawing process
for thermoplastic composite laminates. J Engng Mater Technol
Design, Manufacture, and Industrial Innovation (DMI-
1997;119:314±8.
9900185) and Ford for their support of this research, and [15] Yi YM, Park SH, Youn SK. Asymptotic homogenization of viscoe-
Prof. Julie Chen for helpful discussion. lastic composites with periodic microstructures. Int J Solids Struct
1998;35(17):2039±55.
[16] Takano N, Uetsuji Y, Kashiwagi Y, Zako M. Hierarchical modeling
References of textile Composite materials and structures by the homogenization
method. Model Simul Mater Sci Engng 1999;7:207±31.
[1] Bradaigh CM, McGuinness GB, Pipes RB. Numerical analysis of [17] Kaw AK. Mechanics of composite materials. New York: CRC Press,
stress and deformations in composite materials sheet forming: central 1997.
indentation of a circular sheet. Compos Manufact 1993;4(2):67±83. [18] Halpin JC, Tsai SW. Effects of environmental factors on composite
[2] Boisse P, Cherouat A, Gelin JC, Sabhi H. Experimental study and materials. Air Force Technical Report AFML-TR-67-243, June 1969.
®nite element simulation of a glass ®ber fabric shaping process. [19] Tsai SW. Structural behavior of composite materials. NASA CR-71
Polym Compos 1995;16(1):83±95. 1964.
[3] Gelin JC, Cherouat A, Boisse P, Sabhi H. Manufacture of thin compo- [20] Noyes JV, Jones BH. Analytical design procedures for the strength
site structures by the RTM process: numerical simulation of the and elastic properties of multilayer ®ber composites. Proceedings of
shaping operation. Compos Sci Technol 1996;56:711±8. AIAA/ASME Ninth Structure Dynamics and Material Conference,
[4] Boisse P, Borr M, Buet K, Cherouat A. Finite element simulations of 1968, Paper 68±336.
textile composite forming including the biaxial fabric behavior. [21] Prodromou AG, Chen J. On the relationship between shear angle and
Compos Part B: Engng 1997;28(4):453±64. wrinkling of textile composite performs. Composites Part A
[5] Hou M. Stamping forming of continuous glass ®ber reinforced poly- 1997;28:491±503.
propylene. Composite Part A 1997;28:695±702. [22] McBride TM, Chen J. Unit-cell geometry in plain-weave fabrics
[6] De Luca P, Lefebure P, Pickett AK. Numerical and experimental during shear deformation. Compos Sci Technol 1997;57:345±51.