Composites: Part A 33 (2002) 1095–1105 www.elsevier.com/locate/compositesa
Non-orthogonal constitutive equation for woven fabric reinforced thermoplastic composites Woong Ryeol Yua, Farhang Pourboghrata,*, Kwansoo Chungb, Michael Zampalonia, Tae Jin Kangb a
Department of Mechanical Engineering, Michigan State University, East Lansing, MI 48824, USA b School of Materials Science and Engineering, Seoul National University, Seoul, South Korea Received 7 November 2001; revised 14 May 2002; accepted 30 May 2002
Abstract One of the ultimate objectives of this study was to investigate the feasibility of shaping preconsolidated woven FRT (fabric reinforced thermoplastics) using stamp thermo-hydroforming, a new forming method for composite manufacturing. A new constitutive model has been developed based on a homogenization method by considering the microstructures of composites including both the mechanical and structural properties of fabric reinforcement. In particular, the current model aims to account for the effect of the fiber strength difference and orientation on anisotropy and also to simulate shear deformation without significant length change, common in FRT composite forming. For validation purposes, the model was implemented in an explicit dynamic finite element code and tested for in-plane simple shear, pure shear, uniaxial tension, and draping behavior of woven composites. q 2002 Elsevier Science Ltd. All rights reserved. Keywords: A. Fabrics/textiles; B. Anisotropy; C. Computational modelling; E. Prepeg; Sheet hydroforming
1. Introduction Woven FRT composites are now gaining popularity in manufacturing automotive and aerospace parts due to their superior reinforcing properties, ease of handling and well established textile technologies [1]. Fabric reinforcements fabricated into composites with the thermoplastic resin not only improve mechanical and physical properties but also enable rapid, automated production of composite structures. Numerous manufacturing processes have been proposed over the years to shape thermoplastic composites, ranging from injection molding to diaphragm forming and filament winding. Due to their high success with metals, various attempts have also been made to apply sheet-stamping techniques to composites [2,3]. A difficulty to overcome while applying sheet-stamping for the formation of fabric reinforced composites is the limited draping capability of fabrics, which are prone to wrinkling. Therefore, there exists a current need to develop appropriate constitutive equations for the proper analysis and design of manufacturing processes of FRT composites with to reduce wrinkling. The research on constitutive modeling for woven FRT * Corresponding author. Tel.: þ 1-517-353-0819; fax: þ1-517-353-1750. E-mail address:
[email protected] (F. Pourboghrat).
can be classified into two general categories; continuum mechanics and non-continuum mechanics based models (Fig. 1). The non-continuum approach is used for the modeling of textile preforms without resin. This enables draping simulations before resin injection during resin transfer molding, with particle dynamics representing the intersection point between warp and weft threads [4]. Kinematic analysis [5,6] has been performed without consideration of the constitutive equation for the simulation of the deformation of a single ply of an ideal fabric. In this procedure, the draping simulation was performed by generating a surface, computing fabric nodal point location and displaying the draped surface. In these papers, the authors showed that the entire fabric draping pattern depends solely on the choice of the initial path. Therefore, kinematic analysis was shown to provide a quick evaluation of the wrinkling possibility due to fabric locking during the process design phase for manufacturing composite parts using dry woven fabrics. For woven FRT, the continuum-based model is preferred due to its complex deformation behavior during manufacturing. In general, for composite modeling the material can be treated as either homogenous or non-homogeneous. Based on non-homogeneous concepts, truss and shell
1359-835X/02/$ - see front matter q 2002 Elsevier Science Ltd. All rights reserved. PII: S 1 3 5 9 - 8 3 5 X ( 0 2 ) 0 0 0 5 3 - 2
1096
W.R. Yu et al. / Composites: Part A 33 (2002) 1095–1105
Fig. 1. Hierarchical review of composite modeling research.
(membrane) elements are used to model textile structure and resin matrix, respectively [7,8]. In this meso-cell approach, the constitutive equation considers the microstructural interaction between the matrix and the reinforcement using trusses and shell (or membrane) elements. Therefore, the meso-cell approach can account for the evolution of the microstructures such as fiber angles. However, since this approach takes a large amount of computation and preprocessing time in analyzing the forming processes, the homogenous continuum approach is advantageous over the meso-cell approach. Many models assuming homogeneous material properties have been developed utilizing the well-defined mathematical theory (e.g. orthotropic constitutive equation) thereby reducing the computational cost [9,10]. These studies assume the preservation of the initial orthotropy of woven FRT during all forming processes, which is a common practice in sheet metal forming analysis. This assumption is not proper for the analysis of woven FRT, since the change in fiber angles is significant [11]. Therefore, it needs proper internal variables that account for major microstructural evolution. Several efforts have been made to develop numerical models capable of capturing the fiber angle evolution. This resulted in non-orthogonal continuum-based constitutive equations for unidirectional composites that maintained fiber inextensibility and incompressibility assumptions [12] based on an ‘Idealized Fiber Reinforced Materials’ theory
[13,14]. This approach has been applied to the diaphragm forming analysis of unidirectional thermoplastic composites [15,16] and recently this approach has been extended to thermoplastic composites with woven fabric reinforcements [17]. There was no example offered to show that this model was suitable for the finite element analysis of FRT forming. A great deal of research has been devoted to using homogenization methods to account for the evolution of microstructural parameters, such as fiber angle, to predict the macroscopic deformed shape of composites. Homogenization methods consist of modeling a unit cell in 3D and then solving a set of governing equations to obtain the homogenized material properties of the composite based on the geometrical description of the unit cell [18]. In spite of the computational disadvantages, including long computational time, the homogenization method has been used to model a multitude of domains such as the microlevel of the fiber and the meso-level of the fabric [19]. Recently, a new method has been developed that reduces computational costs and accounts for microstructural effects on the forming behavior of thermoplastic composites. The method utilizes numerical experimental techniques to extract the non-linear elastic material properties based on the homogenization method, all under a continuum mechanics framework [20]. Most models assuming homogenous material properties have used either the orthotropic constitutive equation with an orthogonal material axis or the computationally expensive homogenization method for calculating the homogenized
W.R. Yu et al. / Composites: Part A 33 (2002) 1095–1105
material properties of composites. The aim of this study is to accurately represent the deformation behavior of a woven FRT using a new constitutive equation that is developed based on non-orthogonal material axis under a continuum framework. The current study uses a unit cell structure of woven FRT and a type of homogenization method to extract macromaterial properties by deriving an analytic form for the material properties relating stress and strain. By predicting the preferred (fiber) orientation evolution, the current model can account for the effect of orientation on anisotropy as well as the effect of the fiber strength difference.
2. Constitutive equation 2.1. Modeling of the structural unit of fabric To develop a convenient tool to better understand the physical source of anisotropy in the constitutive equations of composites, a structural net comprised of several warp and weft fibers is considered as shown in Fig. 2. The net has many unit cells separated by the intersection points of the warp and weft fibers, at which warp and weft yarns can rotate freely. Here, fibers are attached to the resin matrix and the net shares a common uniform strain. By applying mechanical properties to the fibers and the matrix, the total forces and stresses are calculated based on the information from the structural parameters (thickness, structural unit cell dimension and fiber angle) and fiber properties. The material stiffness (relating strains and stresses) is then calculated using the kinematics and the force equilibrium of the related structural net in concern. 2.2. Derivation of stress and strain relationship Consider tensile strains imposed along the a fiber direction (the x-direction in the materially embedded coordinate system) and also along the vertical to the a fiber direction (the y-direction). For the prescribed strain under the small deformation assumption, the force acting on the structure can be derived considering the kinematics and fiber properties and then, the stresses are related to the
1097
strain. The small deformation assumption for fibers is not severe in the sense that fiber strain is very small during forming and the resin is soft at the forming temperature. To calculate forces for a prescribed strain, the Lagrangian strain measure and its simplified form are used as follows [21] ds2 2 dS2 ¼ dXT ðFT F 2 IÞdX ¼ dXð2EÞdX
ð1Þ
ðds þ dSÞðds 2 dSÞ ¼ dXð2EÞdX ðds þ dSÞ ø 2dS dXðEÞdX ø ðds 2 dSÞdS
ð2Þ
where ds and dS are the length of a material line element ðdXÞ after and before deformation, respectively, while F and E are the deformation gradient and strain tensor, respectively. Using Eqs. (1) and (2), strains on fibers in specific direction are calculated as follows. 1¼
ds 2 dS dX·E·dX ¼ dS ðdSÞ2
ð3Þ
Consider the uniaxial strain 1x imposed in the a fiber direction with the condition that 1y ¼ 1xy ¼ 0: Note that the formulation given below is based on a local coordinate system consisting of x-axis along a fiber direction and y-axis orthogonal to the x-axis. Through Eq. (3) and notations in Fig. 2, the strains in fiber direction are " # " # " # a 1x 0 l~al a g ; dX ¼ ; dX ¼ ð4Þ E¼ 0 0 b 0 1a ¼ 1g ¼
dXa ·E·dXa ldXa l2 dXg ·E·dXg ldXg l2
¼ 1x ¼
a2 1x ða2 þ b2 Þ
ð5Þ ð6Þ
Then, the forces on the a and g fibers are obtained using linear relationship between fiber strain and stress as follows F a ¼ Aa E a 1x F g ¼ Ag E g
a2 1x a2 þ b2
ð7Þ ð8Þ
where A a and E a are the cross-sectional area and elastic modulus of a yarn in the a direction, respectively, while A g
Fig. 2. A structural net (left) and a unit cell (right) of woven fabric reinforced composites in materially embedded coordinate system.
1098
W.R. Yu et al. / Composites: Part A 33 (2002) 1095–1105
and E g are those of a yarn in the g direction. Note that ða; bÞ is the material element ðdXÞ of the g fiber based on materially embedded coordinate system (Fig. 2), while a~ and b~ are vectors showing the magnitude and direction of the a and g fibers in the unit cell consisting of a pair of warp and weft. To calculate the total forces in the structural net, the total number of fibers or yarns contributing to each direction is needed and these are denoted by nax ; nay ; ngx and ngy : Here, the superscript in the notation indicates specific fibers while the subscript refers to the direction that is normal to a boundary plane. For example, the ngx means the number of g fibers facing the boundary plane normal to x-axis (note that, nay ¼ 0). With these notations, total forces in the structural net due to the prescribed strain become t
a 1x a pffiffiffiffiffiffiffiffiffiffi a2 þ b2 a2 þ b2
a2 1x b pffiffiffiffiffiffiffiffiffiffi 2 2 2 a þb a þ b2
t
Fy ¼ ngy F g sin F ¼ ngy Ag Eg
t
Fyx ¼ ngy F g cos F ¼ ngy Ag Eg
a2 1x a pffiffiffiffiffiffiffiffiffiffi 2 2 2 a þb a þ b2
ð10Þ ð11Þ ð12Þ
ð13Þ
Fxy ¼ ngx F g sin F ¼ ngx Ag Eg
b2 1y b pffiffiffiffiffiffiffiffiffiffi a2 þ b2 a2 þ b2
ð14Þ
Fy ¼ ngy F g sin F ¼ ngy Ag Eg
b2 1y b pffiffiffiffiffiffiffiffiffiffi þ b2 a2 þ b2
ð15Þ
b2 1y a pffiffiffiffiffiffiffiffiffiffi 2 2 2 a þb a þ b2
ð16Þ
a2
Fyx ¼ ngy F g cos F ¼ ngy Ag Eg
Using Eqs. (1) – (3), fiber strain due to prescribed shear strain ð1xy Þ is calculated when 1x ¼ 1y ¼ 0: For a fibers ! ! ! a 0 0 1xy ¼ a 0 ¼0 ð17Þ a 0 1xy 0 0 0 therefore, no strain develops. For g fibers whose direction is ða; bÞ ! ! ! a b1xy 0 1xy ¼ 2ab1xy ð18Þ ¼ a b a b 1xy 0 a1xy b
s
Fxy ¼ ngx Ag Eg
s
Fy ¼ ngy Ag Eg
2ab1xy b pffiffiffiffiffiffiffiffiffiffi a2 þ b2 a2 þ b2
2ab1xy b pffiffiffiffiffiffiffiffiffiffi ; 2 2 2 a þb a þ b2
Fyx ¼ ngy Ag Eg
ð19Þ
2ab1xy a pffiffiffiffiffiffiffiffiffiffi 2 2 2 a þb a þ b2
Now, consider the number of fibers contributing to the total force assuming that there are n joints in the diagonal direction of the two fiber families as shown in Fig. 2. For the purpose of explanation, the way of counting the number of fibers facing a specific plane was described in Fig. 3 in the case of a net containing six joints. Based on the scheme, by generalizing the six joints into n joints, nax and nay could be expressed as follows ðn 2 1Þ‘ sin f ‘ sin f þ 1 ¼ ðn 2 1Þ þ 1 ¼ n ð22Þ b ‘ sin f
while nay ¼ 0: For g fibers ðn 2 1Þ‘ sin f ngx ¼ MaxInt þ 1 ; h ð23Þ
ðn 2 1Þ‘ cos f g þ1 ny ¼ MaxInt a
where MaxInt{A} indicates the largest integer among those smaller than the value A and a ¼ l~al: After substituting Eqs. (22) and (23) into Eqs. (9) – (16), (20) and (21), the total force for the 1x, 1y and 1xy are obtained by adding up Eqs. (9) –(16), (20) and (21). By taking the limit value for n ! 1 of the normalized total force, the stress becomes
sxx ¼ lim
n!1
Fx ; cn‘ sin f
sxy ¼ lim
n!1
Fxy ; cn‘ sin f ð24Þ
therefore, the strain on the g fiber becomes 2ab1xy ¼ dSðds 2 dSÞ
2ab1xy a pffiffiffiffiffiffiffiffiffiffi ; 2 2 2 a þb a þ b2 ð20Þ
nax ¼
b2 1y a t pffiffiffiffiffiffiffiffiffiffi Fx ¼ ngx F g cos F ¼ ngx Ag Eg 2 2 2 a þb a þ b2
t
Fx ¼ ngx Ag Eg
ð9Þ
Total forces for strain 1y with 1x ¼ 1xy ¼ 0 are also similarly obtained as
t
s
s 2
a2 1 b t ffi Fxy ¼ ngx F g sin F ¼ ngx Ag Eg 2 x 2 pffiffiffiffiffiffiffiffiffi 2 a þb a þ b2
t
The total forces contributed by the shear strain become
ð21Þ
Fx ¼ nax F a þ ngx F g cos F ¼ nax Aa Ea 1x þ ngx Ag Eg
2ab1xy ds 2 dS ¼ 2 dS a þ b2
syy ¼ lim
n!1
Fy ; cn‘ cos f
syx ¼ lim
n!1
Fyx cn‘ cos f
W.R. Yu et al. / Composites: Part A 33 (2002) 1095–1105
1099
Fig. 3. Schematic diagram for calculating nax ; nay ; ngx and ngy :
which results in 2 3 sxx 6 7 6 syy 7 4 5 sxy 2
2! a a E~ a 6 þG 6 bc c h 6 6 2! 6 6 b a ¼6 G 6 c a 6 6 ! 6 6 b a2 4 G c h
2
1xx
the jamming and frictional effect has originated from the fact that the fictitious fibers, orthogonal to the a fibers, can prevent the a and g fibers from being unrealistically too close (Fig. 4). The equation can be derived by following the procedure mentioned, Eq. (27), and incorporated into Eq. (24) to construct Eq. (25). 2! a b G c h 2! E~ b b b þ G a c c a 2! b b G c h
3 a ab 7 G c 7 h 7 7 7 b ab 7 7 G a c 7 7 7 7 b ab 7 5 G c h
3
6 7 7 £6 4 1yy 5 21xy
ð25Þ
f
Fy ¼ Eb Ab 1y ¼ E~ b 1y ;
syy ¼ lim
n!1
Fy ¼ nby E~ b 1y ;
nE~ b 1yy E~ b ¼ 1 nca ca yy
Therefore, large value of E~ b brings large shear resistance in the woven structure and thus results in small shear deformation. The values of E~ b may be experimentally determined as a non-linear function of the fiber angle from experiments such as the picture frame shear test [22] but here a specific value is selected to fit the experimental results well.
where E~ a ¼ Ea Aa ; G ¼ Ag Eg =ða2 þ b2 Þ3=2 and c is the thickness of composite. Note that the constitutive equation is symmetric since 2! 2! a b b a ¼G ; G a c c h ð26Þ 2! a ab b a G ¼G c c h h In Eq. (25), E~ b is added to account for the jamming and frictional effect of the fabric structure. This idea for treating
ð27Þ
Fig. 4. Physical meaning of arbitrary constant, Eb :
1100
W.R. Yu et al. / Composites: Part A 33 (2002) 1095–1105
2.3. Modeling of resin material behavior Existence of the resin could be ignored in modeling composites because the resin stiffness is much less compared to the reinforcement stiffness. In the thermoforming of thermoplastic composites, the role of the resin is more important since the heat generated during forming changes the phase of the resin into a molten state, thus enabling the rheological flow of the resin and the reinforcement. Therefore, it might be necessary to model the resin material as a viscous material. Since the current study focuses on the development of the constitutive equation and its validation, the resin property is assumed as an isotropic elastic material under the plane stress conditions. Using this assumption, the constitutive equation for the resin matrix becomes
s11 ¼ s22
EM ðð1 2 nÞ111 þ n122 þ n133 Þ ð1 þ nÞð1 2 2nÞ
with n ¼ 0:5: Therefore, the constitutive equation for the resin matrix under the plane stress condition is
s22 ¼
2EM 1 ; 3 22
s12 ¼
ð32Þ
The two direction vectors determine important parameters such as the fiber angle and direction of the fiber as shown in Fig. 2. Eq. (25) is updated prior to the next stress calculation with the new values of a and b. Note also that Eq. (25) is expressed in terms of the materially embedded coordinate system aligned with the a fiber direction.
3.1. In-plane shear deformation behavior
ð29Þ
EM 2EM ðð1 2 2nÞ111 Þ ¼ 1 ð1 þ nÞð1 2 2nÞ 3 11
~ old b~ new ¼ F·b
For validation purposes, several simulations were performed and the results were compared with either experiments or other published results.
where EM and n are the elastic modulus and Poisson’s ratio of the resin matrix, respectively. From Eq. (28) and the incompressibility condition, the thickness strain is
s11 ¼
ð31Þ
3. Results and discussion
EM ðn111 þ n122 þ ð1 2 nÞ133 Þ ¼ 0 ð1 þ nÞð1 2 2nÞ
133 ¼ 2ð111 þ 122 Þ
a~ new ¼ F·~aold
ð28Þ
EM ðn111 þ ð1 2 nÞ122 þ n133 Þ ¼ ð1 þ nÞð1 2 2nÞ
s33 ¼
These vectors are updated using the deformation gradient tensor from the incremental deformation as shown in Eqs. (31) and (32). Note that the updating of these vectors using the rotation tensor would preserve orthogonality.
ð30Þ
EM E ð2112 Þ ¼ M g12 ; 3 3
s33 ¼ 0 This equation is combined with Eq. (25) to obtain the constitutive equation for FRT composites.
For the first example, a square sheet of woven FRT composite was subjected to a shear boundary force on its right-hand side edge while the left-hand side edge was constrained not to move in the x-direction but freely move in the y-direction. To prevent the rigid body motion of the composite, the lower left corner was completely fixed as shown in Fig. 5. A total of 10 £ 10 rectangular elements (CPS4R element in ABAQUS) were used with the material properties summarized in Table 1 and a constant value 10 ðE~ b ¼ 10 NÞ was used in the current constitutive model. To compare the predicted shape by using current nonorthogonal constitutive equation with the orthotropic constitutive equation, two results were shown in Fig. 6(a) and (b). For the orthotropic constitutive model, a small shear modulus was inputted with large Young’s moduli in both x and y material directions. The large difference between
2.4. Implementation of the constitutive equation into FEM code To predict the overall deformed shape of a FRT during forming, the developed constitutive equation was implemented into the ABAQUS explicit code using the user material subroutine (VUMAT). The constitutive equation described so far is suitable for implementation in any finite element analysis package in that the equation is given as an explicit matrix form. This feature is crucial for incorporating contact and strain rate effects into forming analysis. ~ expressing the magnitude Two vector quantities, a~ and b; and direction of two families of fibers in a unit cell, define the main internal variables for the constitutive equation.
Fig. 5. Geometry and boundary condition for simple shear simulation.
W.R. Yu et al. / Composites: Part A 33 (2002) 1095–1105
1101
Table 1 Properties of the woven fabric prepreg used in this study [6] Material properties Young’s modulus of fibers (S-2 glass)
69 GPa
Structural parameters of fabric Yarn thickness Yarn width Yarn spacing Ratio 0 –908
0.41 mm 5.33 mm 6.40 mm 53:47
Parameters in Eq. (21) E~ a
1.5 £ 105 N
~g
1.3 £ 105 N
E EM a~ b~ Aa b
A
0.001 MPa (0.0064, 0.0) m (0.0, 0.0064) m 2.19 mm2 1.94 mm2
tensile and shear modulus, which represents two strong threads in warp and weft, and weak shear resistance, are often used in predicting the deformation of woven composites using the orthotropic material assumption. The deformed shape obtained using current non-orthogonal constitutive equation shows no significant length change in warp and weft directions, i.e. shear dominant deformation shape on cross points between warp and weft threads. Fig. 6(a) shows that the nature of the woven FRT composite is properly simulated, i.e. in simple shear test the woven FRT deforms into a parallelogram without significant length change. In the simulation using orthotropic constitutive equation, the same deformation shape could be expected due to the material properties used, i.e. very small shear modulus compared to tensile modulus, however, the predicted shape in Fig. 6(b) shows significant length change in the left edge. This length change comes from the assumption that the material axes remain orthogonal during deformation process. In finite element analysis, orthotropic
Fig. 6. Simple shear simulation (a) with non-orthogonal constitutive equation, (b) orthotropic constitutive equation (E11 ¼ 19 GPa; E22 ¼ 19 GPa; n ¼ 0:1; G ¼ 0:01 MPa).
Fig. 7. Update scheme of material axis under orthotropic constitutive assumption in finite element analysis.
material axes are updated by the rotation tensor, which results in wrong identification of material axes during deformation process (Fig. 7). This issue is also discussed in the pure shear example below. To investigate the pure shear deformation behavior of the new constitutive equation, 10 £ 10 finite element (CPS4R element in ABAQUS) meshes were used for a square piece of the FRT composite. The shear force was loaded on the four sides of the FRT composite. The lower left corner node was constrained in all directions to prevent rigid body rotation (Fig. 8). Dong et al. [11] used this problem to compare their material law, which is comprised of an averaging of two orthotropic constitutive equations, to model the fabric reinforced composite. Fig. 9 shows the deformed shape calculated using the orthotropic material property provided as a material option by the ABAQUS Explicit code. As shown in Fig. 9(a), due to the elongation of the boundary, it is apparent that this material law could not simulate pure shear deformation. Fig. 9(b) shows the
Fig. 8. Geometry and boundary condition for pure shear simulation.
1102
W.R. Yu et al. / Composites: Part A 33 (2002) 1095–1105
Fig. 9. Pure shear simulation using the orthotropic constitutive equation with the unidirectional lamina properties (E11 ¼ 19 GPa; E22 ¼ 79 MPa; n ¼ 0:1; G ¼ 0:01 MPa): (a) deformed shape and (b) material direction.
Fig. 11. Pure shear deformation with the new constitutive equation: (a) strain in the 1-fiber direction and (b) strain in the 2-fiber direction after deformation.
orthogonal material directions on the mesh, which was updated using the rotation tensor, thus remaining orthogonal for the entire deformation, though material axes were changed by pure shear deformation. Fig. 10 shows the deformed shape obtained under pure shear loading using the current non-orthogonal constitutive equation. With the new constitutive equation, no significant length change is observed in the specimen boundary (Fig. 10(a)). This was possible because fiber directions were properly updated during the shear deformation as shown in Fig. 10(b). In Fig. 11(a) and (b), strains in two fiber directions, always non-orthogonal, are plotted with near zero values. The success with the pure shear simulation demonstrates that the current constitutive equation can predict the deformation behavior of fabric reinforced thermoplastics at high temperatures. At high temperatures, the main deformation mode is the shearing behavior of two fiber families at intersection points.
symmetry condition, only one quarter of the specimen was modeled using 900 shell elements (S3R element in ABAQUS) with the material properties given in Table 1. The right edge of the preform was constrained by assigning zero displacement to all the nodes in the vertical direction. The friction effect at the intersection points within the preform was determined by fitting the fiber angle for the center point (W in Fig. 13) to the experimental measurement. Thus two E~ b values were determined numerically by fitting the fiber angle for the center point (W in Fig. 13) to the experimental measurement. Before fiber jamming angle [8], E~ b ¼ 15 was used while it increased to E~ b ¼ 2 £ 104 after fiber jamming angle. Fig. 12 shows that the comparison between the predicted and measured center point fiber angles at different percentage elongations are in good agreement. Fig. 13 shows the simulated and the actual deformed shapes of the preform at 25% elongation. The simulated shape shows three distinct deformation zones with severe local change in fiber angles. Fig. 14 shows more clearly the contour of fiber angle distribution for both the experiment and the analysis. Sidhu et al. [8] identified various modes of deformation taking place in uniaxial tensile test in three different zones.
3.2. Uniaxial extension of woven fabric preform Uniaxial extensions were numerically performed on a plain weave specimen without resin, i.e. EM ¼ 0; to quantitatively compare the results of the current constitutive model with experiments [8]. The preform size was 203.2 mm £ 406.4 mm (8 in. £ 16 in.) and fiber orientation of the preform was ^ 458 based on elongation direction. Due to the
Fig. 10. Pure shear deformation with the new constitutive equation: (a) deformed shape and (b) two fiber direction.
Fig. 12. Comparison of center fiber angle between experiment [8] and analysis.
W.R. Yu et al. / Composites: Part A 33 (2002) 1095–1105
1103
Fig. 13. Deformed shape of tensile specimen at 25% elongation. (a) Experiment [8]: whole domain, (b) analysis: quarter domain.
In Zone I, there is no significant deformation and thus the fiber angle stays the same as in the undeformed specimen. Materials in Zone II exhibit both shearing and sliding between yarns and some fiber change. In Zone III, the main deformation is by shearing only. An almost exact deformation behavior of the tensile specimen can be observed from the simulated results as shown in Fig. 14. 3.3. Draping behavior of woven fabric preform This simulation was done to investigate the drapability and forming behavior of woven plain fabric preform by using the non-orthotropic constitutive equation. The simulated results were compared with experiments with respect to the boundary profile and fiber angle distribution of the draped preform. A finite element model for the draping simulation was constructed as shown in Fig. 15(b) by using the features of the original apparatus in Fig. 15(a) used in Ref. [23]. The radius of the upper mold in this model was 100 mm while the lower one was 97 mm. As the upper mold moved downwards, the fabric was draped over the lower mold, which consisted of a hemispherical dome surrounded by a flat base. The distance from the lower mold’s center to the edge of the flat base is 260 mm. A square 320 £ 320 mm2 of woven plain fabric was used for the experimental study of draping behavior [23]. For the simulation, 1250 shell elements (S3R in ABAQUS) were used to model a quarter section of the woven preform due to the symmetry. The material properties for the fabric blank were calculated based on the structure and material properties of the woven preform used in the experiment [23]. The friction coefficient between the molds and the fabric blank were chosen from literature [10], which experimentally determined those values according to mold material type such as aluminum and Perspex. The material and process variables used in the analysis are listed in Table 2.
Fig. 14. Comparison between fiber angle distribution from (a) experiment [8] and (b) analysis (angle unit: radian).
Fig. 16 shows the edge profile of the deformed woven preform. The simulated profile is very close to the experimental one, although there is a little disagreement in the weft direction. This is due to the assumption that there is no material and structural differences between warp and weft directions in the current simulation. However, the real fabric may have had different structural and material properties in the warp and weft directions, resulting in the minor disagreement in the weft direction. The local fiber orientation is shown in Fig. 17. The fiber orientation was contoured experimentally by measuring the orientation of the deformed gridlines inscribed on the fabric
1104
W.R. Yu et al. / Composites: Part A 33 (2002) 1095–1105
Fig. 16. Comparison between boundary profiles of draped preform from experiment [23] and analysis.
Fig. 15. (a) Experimental apparatus for the draping of woven fabric. This picture was copied and pasted in this paper from Ref. [23], page 1410. (b) Finite element model for draping simulation of woven preform in current work.
prior to draping [23]. When considering the fiber angle contour for the simulated draped fabric, a dark portion of 908 fiber angles is observed around the apex of the lower mold, which corresponds to the zero-shear zone. At the equator, where the hemispherical dome changes to the flat Table 2 Material and structural properties of woven fabric used in draping simulation Sheet material
E-glass plain weave
a Fiber strength ðE~ a Þ g Fiber strength ðE~ g Þ Hardening parameter ðE~ b Þ Resin modulus (E M)
90 kN
0.0
Structural parameters Thickness Yarn space
0.58 mm 6.25 mm
Friction coefficient Upper mold-sheet Lower mold-sheet
0.5 0.2
90 kN 10 N
Fig. 17. Comparison between fiber angle distribution from (a) experiment [21] and (b) analysis.
W.R. Yu et al. / Composites: Part A 33 (2002) 1095–1105
surface, the fabric deformation is maximum and a fiber angle of 35.28 was predicted, compared to an experimentally measured fiber angle of 348 [23]. From these comparisons, it could be concluded that the current constitutive equation can accurately simulate the draping behavior of woven FRT composites over a rigid mold.
[4]
[5]
[6]
4. Conclusions [7]
Based on the homogenization method, a constitutive equation suitable for the thermoforming analysis of fabric reinforced thermoplastic composites has been developed, considering the microstructure of fabric reinforcement. In this new constitutive equation, the microstructural information such as the fiber angle was introduced using internal variables. For validation purposes, in-plane simple shear, pure shear, uniaxial extension, and 3D draping problems were investigated. The solutions were obtained using the ABAQUS explicit finite element code by incorporating the current constitutive model into the code using the user material subroutine (VUMAT). The results demonstrated that the current model properly accounts for the effect of the differences in fiber strength and orientation on anisotropic behavior. As a result, shear dominant behaviors with no length change were successfully simulated with the new constitutive equation. In a future paper, the authors will describe a non-orthogonal material model which includes the effect of fiber undulations in the weft and warp directions.
[8]
[9]
[10]
[11]
[12]
[13] [14] [15]
[16]
Acknowledgements This work was supported by the postdoctoral fellowship program of Korea Science and Engineering Foundation (KOSEF) and by the State of Michigan Research Excellence Fund through Composite Materials and Structures Center (CMSC) at Michigan State University. This work was also partially supported by the Ministry of Science and Technology in Korea through the National Research laboratory for which the authors feel grateful.
References
[17] [18]
[19]
[20]
[21] [22]
[1] Long A. Process modelling for textile composite. EUROPAM 2000; 2000. [2] O’Bradaigh CM. Sheet forming of composite materials. In: Advani SG, editor. Flow and rheology in polymer composites manufacturing. Amsterdam: Elsevier; 1994. p. 517–69. [3] Krebs J, Friedrich K, Bhattacharyya D. A direct comparison of
[23]
1105
matched-die versus diaphragm forming. Composites Part A 1998;29: 183–8. Breen DE, House DH, Wozny MJ. A particle-based model for simulating the draping behavior of woven cloth. Textile Res J 1994; 64(11):663 –85. VanWest BP, Pipes RB, Keefe M, Advani SG. The draping and consolidation of commingled fabrics. Compos Manufact 1991;2(1): 9–22. VanWest BP, Pipes RB, Keefe M. A simulation of the draping of bidirectional fabrics over arbitrary surfaces. J Text Inst 1991;81(4): 448–60. Billoet JL, Cherouat A. New numerical model of composite fabric behavior: simulation of manufacturing of thin composite part. Adv Compos Lett 2000;9(3):167 –79. Sidhu RMJS, Averill RC, Riaz M, Pourboghrat F. Finite element analysis of textile composite preform stamping. Compos Struct 2001; 52:483–97. Luca P, Lefebure P, Pickett AK. Numerical and experimental investigation of some press forming parameters of two fiber reinforced thermoplastics: APC2-AS4 and PEI-CETEX. Composites, Part A 1998;29:101– 10. Dong L, Kekakou C, Bader MG. Solid-mechanics finite element simulations of the draping of fabrics: a sensitivity analysis. Composites Part A 2000;31:639–52. Dong L, Lekakou C, Bader MG. Processing of composites: simulations of the draping of fabrics with updated material behavior law. J Compos Mater 2001;35(2):138–63. O’Bradaigh CM, Pipes RB. A finite element formulation for highly anisotropic incompressible elastic solids. Int J Numer Meth Engng 1992;33:1573 –96. Rogers TG. Rheological characterization of anisotropic materials. Composites 1989;20(1):21–7. Spencer AJM. Deformations of fiber reinforced materials. Oxford: Clarendon Press; 1972. McGuinness GB, O’Bradaigh CM. Effect of preform shape on buckling of quasi-isotropic thermoplastic composite laminates during sheet forming. Compos Manufact 1995;6(3–4):269–80. O’Bradaigh CM, McGuinness GB, Pipes RB. Numerical analysis of stresses and deformations in composite materials sheet forming: central indentation of a circular sheet. Compos Manufact 1993;4(2): 67–83. Spencer AJM. Theory of fabric-reinforced viscous fluids. Composites, Part A 2000;31:1311–21. Hsiao SW, Kikuchi N. Numerical analysis and optimal design of composite thermoforming process. Comput Meth Appl Mech Engng 1999;177:1 –34. Takano N, Uetsuji Y, Kahiwagi Y, Zako M. Hierarchical modeling of textile composite materials and structures by the homogenization method. Modeling Simul Mater Sci Engng 1999;7:207 –31. Peng XQ, Cao J. Numerical determination of mechanical elastic constants of textile composites. Proceedings of International Symposium on Affordable Composites Manufacturing, 2000. p. 677 –86. Lai WM, Rubin D, Krempl E. Introduction to continuum mechanics. New York: Pergamon Press; 1982. McGuinness GB, O’Bradaigh CM. Development of rheological models for forming flows and picture-frame shear testing of fabric reinforced thermoplastic sheets. J Non-Newtonian Fluid Mech 1997; 73:1–28. Mohammed U, Lekakou C, Bader MG. Experimental studies and analysis of the draping of woven fabrics. Composites, Part A 2000;31: 1409– 20.
