Numerical generation of a random chopped fiber composite RVE and its elastic properties

Composites Science and Technology 68 (2008) 2792–2798

Contents lists available at ScienceDirect

Composites Science and Technology journal homepage: www.elsevier.com/locate/compscitech

Numerical generation of a random chopped fiber composite RVE and its elastic properties Yi Pan, Lucian Iorga, Assimina A. Pelegri * Mechanical and Aerospace Engineering, Rutgers, State University of New Jersey, 98 Brett Road, Piscataway, NJ 08854, USA

a r t i c l e

i n f o

Article history: Received 6 February 2008 Received in revised form 12 May 2008 Accepted 3 June 2008 Available online 12 June 2008 Keywords: A. Short-fiber composites B. Elastic properties C. Finite element analysis (FEA) C. Modeling C. Multiscale modeling

a b s t r a c t The elastic properties of random chopped fiber-reinforced composites (RAFCs) are of paramount importance for their sound application in lightweight structures. Mass production of random chopped fiberreinforced composite (RaFC) at a fraction of the cost of composite laminates establishes RaFCs as alternate candidate materials for manufacturing lightweight components in the automotive industry. Nevertheless, understanding and modeling of their mechanical and fracture properties are still fields of active research, yet to be exhausted. In this paper, methods to generate an RVE for random fiber or particle reinforced composites numerically are reviewed. A modified random sequential absorption algorithm is proposed to generate a representative volume element (RVE) of a random chopped fiber-reinforced composite (RaFC) material. It is assumed that the RVE represents the composite material within the framework of elasticity. The RVE thus created is analyzed to obtain the mechanical properties of the composite material by using finite element analysis (FEA). RVE generation uses both straight and curved fibers so as to achieve high fiber volume fractions (VFs) that are extremely difficult to obtain by using straight fibers alone. This work extends the capability of RVE generation of RaFCs to higher volume fractions, here 35.1% is illustrated, which are in the range of values employed in industrial applications. Ó 2008 Elsevier Ltd. All rights reserved.

1. Introduction Significant efforts have been devoted toward the use of lightweight structures to increase energy efficiencies in various industrial and commercial sectors [1–4]. Fiber-reinforced composites have found numerous applications in aerospace industry for their high specific strength and specific stiffness [5]. However, the cost of traditional composite materials is also considerable. Random chopped fiber-reinforced composites (RaFCs) have emerged as promising alternative materials for lightweight structures due to their low cost and mass production capabilities. Their potential application in, for example, automotive industry has been documented [1,2,4]. In order to expand their use, accurate material characterization is required. The main difficulty in fully exploring the capabilities of the RaFCs lies in the apparent impediment to effectively model their geometry at the micro-level for high fiber volume ratios (35–40%). This difficulty becomes even more obvious at high aspect ratio (AR) fibers. Among several notable methods to predict the constitutive behavior of random-fiber composites the analytical homogenization method based on the Eshelby’s strain concentration tensor [6], where the fiber is treated as a second phase inclusion, enjoys * Corresponding author. Tel.: +1 732 445 0691; fax: +1 732 445 3124. E-mail address: [email protected] (A.A. Pelegri). 0266-3538/$ - see front matter Ó 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.compscitech.2008.06.007

considerable support. Significant success has been enjoyed by techniques based on the Mori-Tanaka’s mean field method [7]. The aforementioned methods are employed to determine the material properties for unidirectional composites with fiber volume fraction (VF) and aspect ratio similar to those of the random composite analyzed. Orientation averaging, as developed by Advani and Tucker [8], has to be performed in order to account for the random orientations. A different set of methodologies for the characterization of RaFC response is based on the classical laminated plate theory (CLPT), as illustrated by ‘‘Laminated Random Strand Method” proposed by Ionita and Weistman [9] for the analysis of high volume fraction carbon fiber/urethane composites. Due to their computational efficiency, such methods can readily be employed in moving-window analyses to statistically assess the analysis sample size effects on the simulated response of the material. However, this method does not provide any information on the micro-mechanical stress–strain state. A third approach is enabled by the availability of powerful computational hardware and advanced finite element packages. It is possible to simulate the material behavior directly using direct 3D finite element analysis (FEA) [11–19] on a statistically representative geometric entity, referred to as representative volume element (RVE). Of the aforementioned material characterization methods, namely, FEM, field based methods and quasi-analytical methods, the FEM based schemes are the most promising. The key issue of

Y. Pan et al. / Composites Science and Technology 68 (2008) 2792–2798

the FEM approach is identification and generation of an RVE1 through which the effective homogeneous material properties are derived. The first definition of RVE is that an RVE is a statistical representation of material [20]. In this context, an RVE must be sufficient large such that it contains a large number of inclusions in a heterogeneous material, and the effective properties derived from the RVE represent the true material properties in macroscopic scale. As pointed out by Drugan and Willis [21], another pragmatic definition of RVE is that ‘‘the smallest material volume element of the composite for which the usual spatially constant ‘overall modulus’ macroscopic constitutive representation is a sufficiently accurate model to represent mean constitutive response”. Adopting the second definition of RVE, Drugan and Willis [21] showed that this definition yield small quantitative estimates of RVE size while in contrast to qualitatively large RVE size implied by the first definition. The determination of RVE size remains an open question. In FEM analyses toward the determination of the effective material properties, both of the above mentioned definitions come into effect. On one hand, the first definition guarantees the RVE is large enough to represent the material macroscopically. On the other hand, a large RVE is computationally expensive and the RVE size should be minimized to increase the computational efficiency. In practical FEM simulation, which can be seen in works by Iorga et al. [10], Gusev et al. [13], and Trias et al. [22] and Kari et al. [15], a series of various sizes of RVEs are generated and the material properties obtained from which are compared to determine a relatively small size for the RVE within certain accuracy. The remaining question is how the RVE micro-level geometry on which FEM analyses are performed is identified and generated numerically. It depends on the microstructure of the material. While identification of an RVE (or UC) is immediate in woven fiber [23,24], and laminate [25] composites since they have a repeatable architecture. For composite materials with random inclusion phase arrangements the recognition and numerical construction of an RVE is not straightforward [13–16,18,26]. Knowledge of the shape, orientation distribution and/or location of the reinforcing fiber/ particle is required a priori. Geometric shapes corresponding to the reinforced phase particles are inserted in a volume box in such a way that they follow the predefined location and/or orientation distribution. In this paper, the composite is a two phase material and the reinforced phase is referred to fiber. There are three families of methods for numerically generating RaFC (random-fiber composite) RVEs, namely, the random sequential adsorption (RSA) schemes [15,16,26–29], the Monte Carlo procedures [13,14,18] and image reconstruction technique [17]. An RSA-type scheme sequentially adds fibers to a region by randomly generating its location and orientation angle. The new fiber is not allowed to intersect fibers previously accepted [15,16,26,30]. Böhm et al. [16] employed a modified RSA approach to generate RVEs of metal matrix composites (MMC) with 15% volume fraction (VF) reinforced by random short fibers with aspect ratio (AR) equal to five (AR = 5) and by random spherical particles. Fibers or particles were not allowed to overlap and the minimum distance between neighboring fibers or particles was set at 0.0075 times the side length of the RVE. Geometrical periodicity was also applied, that is, parts of fiber that exceeded the surfaces of the cube are cut and shifted so that they enter the cube from the opposite surfaces. Tu et al. [26] also used a similar RSA approach to generate RVEs for studying the thermo-conductivity and elastic modulus of composites reinforced by fibers (AR = 7) and spherical particles. Pan et al. [30] applied a modified RSA algorithm for a random-fiber 1

Note that, there is no consensus on the terminology used for the RaFCs in the literature. This statistical depiction is referenced as a representative volume element (RVE), Kari et al. [15], Gusev [13,14], or unit cell (UC), Böhm et al. [16], and Duschlbauer et al. [18].

2793

composite with composite with fiber AR = 10 and 13.5% fiber volume fraction. We note that the volume fraction achievable through RSA is much smaller than that predicted by existing analytical and phenomenological models. This phenomenon is called ‘‘Jamming” [29]. To overcome the jamming problem and to achieve higher volume fraction, Kari et al. [15] used different sizes of fibers by depositing them inside the RVE in a descending manner. They first deposited the largest aspect ratio fibers and after reaching the jamming limit, then deposited the next largest possible aspect ratio fibers in the RVE. However, such an approach is not applicable for composites with fibers of fixed aspect ratio. Monte Carlo methods are two-step procedures that start from a configuration with arbitrary fiber locations and orientations within a large box and rearrange the fibers’ locations and orientation, without accepting intersections, until a predetermined desired orientation state is reached. Finally, it decreases the size of the box toward the designated fiber volume fraction without altering the orientation of the fibers. Intersection of fibers is not allowed at either step during the process [13,14,18]. The MC procedure has been extensively used by Gusev and his coworkers [13,14,31] and Duschlbauer et al. [18] for studying the material properties of short fiber and particle reinforced composites. A ‘‘jamming” limit similar to that occurring in RS algorithms was also observed in this method. The fiber volume fraction of the RVE generated by MC was up to 21% for fiber aspect ratio 10 [18]. It is noted that for larger fiber aspect ratio, the maximum achievable volume fraction is even smaller. A recent advancement in RVE generation is employing X-ray tomography and 3D image analysis and reconstruction methods as demonstrated for wood-based fiberboards [17]. In this method, microtomographic images of the material’s microstructure are obtained with a synchrotron radiation tomograph. 3D fiber network are reconstructed and generated based on the X-ray absorption radiographs of a sample using specialized software. However, the use of such an approach for generating numerical representations of carbon or glass random composites has yet to be demonstrated. The main drawback of the RSA or Monte Carlo procedures is the small fiber volume fraction that can be achieved. Consequently, they cannot be applied directly for the study of composites with straight fibers of fixed aspect ratio. This calls for a solution to overcome the jamming limit and in the present work we introduce a method to account for fiber kinks that occur in high fiber volume fraction composites by making use of curved fiber geometries in the RVE generation algorithm. This allows for a significant increase in the achievable fiber VF, close to values of 35–40%, representative of typical production composites [32]. Such figures have so far proved elusive using the geometry generation approaches currently employed in the literature. In the second part of the study we employ a homogenization scheme to estimate the macro-scale material properties of the simulated material. The RVE is generated by using an RSA-inspired scheme, which accepts new fibers in the volume according to the fiber geometry (aspect ratio 20) and orientation. The scheme is modified to account for curved fibers in order to achieve high fiber volume fractions (35.1% in the present study). For similar values of the aspect ratio, existing analytical [33,34], experimental [35] and numerical [28] models predict a smaller maximum achievable volume fraction (18.5–30%) for RVEs employing straight cylindrical fibers. The numerical results are compared and validated using Halpin–Tsai’s analytical method [36].

2. Generation of 3D RVE using modified RSA The generation of 3D representative volume element (RVE) of a random chopped fiber composite requires the development of a

2794

Y. Pan et al. / Composites Science and Technology 68 (2008) 2792–2798

comprehensive algorithm in order to capture the ‘‘random” and ‘‘curved” characteristics of fibers existing in the real material. The protocol to generate the RVE is summarized in Section 2.1. Penetration of any two fibers is not allowed physically and that is reflected in the protocol. The mathematical and algorithmic implementation of 2D and 3D fiber penetration tests is discussed in Sections 2.2 and 2.3, respectively. 2.1. RVE generation protocol In the study of particle packing and particle reinforced composites the random sequential adsorption (RSA) technique has been used for spheres [16,26,28,37], spherocylinders [16,28], ellipsoids and rods [15,16,26,28]. For short fiber-reinforced composites with small aspect ratio (<10) and low volume fraction (<25%), the fibers, or fiber bundles, can be modeled as straight cylinders. The relation between the fiber aspect ratio and the maximum achievable fiber volume fraction has been studied by Evans and Gibson [33], Parkhouse and Kelly [35], Toll [34] and Williams and Phillipse [28], among others. Their results indicate that increase of the fiber aspect ratio results in decrease of the maximum achievable fiber volume fraction. For a fiber with aspect ratio 20, these models predict maximum fiber volume fraction of 20% by Evans and Gibson [33], 30% by Parkhouse and Kelly [35], 18.5% by Toll [34], and 27% by Williams and Phillipse [28]. However, the above predicted volume fractions are relatively small compared to the values observed in industrial applications (35–40%) [32,38]. To increase the fiber volume fraction, the bends of the fiber bundle has to be accounted for. Based on our knowledge about the manufacturing process – during which E-glass fiber bundles (also referred to as fibers in this text) are chopped and sprayed into a mold – as well as on optical observations of material specimens, we assume that the fiber is uniformly distributed in the x–y plane and that the fibers have an elliptical cross-section. The in-plane angle of a fiber is in the range of (0, 2p) while the location of its midpoint is randomly generated within a 3D box with uniform probability. The portion of a fiber that is outside the box will be trimmed. The in-plane (x–y) and through thickness (along z-direction) views of the E-glass/epoxy composite specimen (3 mm in thickness) are shown in Fig. 1a and b. The crossing fibers in the highlighted region of Fig. 1a illustrate fibers that bend at the crossing region with other bundles because they cannot physically penetrate each other. After the crossing, the two fibers can continue to be lying within the same layer again (as shown in the highlighted region.) Fig. 1a also indicates the irregular orientation of the fibers. The elliptical shape of the fibers crosssection is demonstrated in Fig. 1b. It is further demonstrated that the majority of the fibers are lying in a similar manner with their minor axis point-

ing in the z-direction. Therefore, we assume the fibers can be either curved or straight, as well as, that a fiber can only bend into and out of the layer such that the project ion of the fiber on to x–y plane remains straight. Based on the above observations and assumptions, the spine of the fiber is approximated by connected straight line segments while its elliptical cross-section is approximated by a dodecagon inscribed to an ellipse with major axis 2a and minor axis 2b, see Fig. 2. This way, a fiber is generated by sweeping the cross-sectional profile along the spine. The protocol to generate the RVE is described next: (1) A fiber-rich sub-layer of 2b thickness (bundle’s minor axis) accommodates straight fibers and portions of the curved fibers (other portions are with in a neighboring sub-layer). A thin matrix-rich sub-layer with the thickness of b/10 is introduced to separate fiber-rich sub-layers. (2) A fiber-mat layer has two fiber-rich sub-layers and three matrix-rich spacing sub-layers, as shown in Fig. 3. The curved fibers can bend away from one layer into the opposite sub-layer and bend back to the previous sub-layer when necessary (for example, to avoid intersecting other existing fibers). In order to have a fiber population sufficiently large, three fiber-mat layers are stacked up to obtain the final RVE, as shown in Fig. 3, ensuring that the two largest eigen values of the second order orientation tensor [8] for all fibers in the three layers are as close to 0.5 as possible. We note

Fig. 2. Dodecagonal approximation of fiber bundle’s elliptical cross-section with major axis 2a and minor axis 2b.

Fig. 1. (a) In-plane and (b) through thickness views of the random chopped E-glass fiber-reinforced composite. The highlighted region in (a) shows the fiber curve to avoid intersecting (highlighted region). In (b) the elliptical cross-section of fibers lying with their minor axis pointing along the z-direction.

2795

Y. Pan et al. / Composites Science and Technology 68 (2008) 2792–2798

Fiber-rich layer

a

Matrix-rich layer

S1

b/10

IA

Top fiber-rich sub-layer

2b Fiber-mat-1

Bottom fiber-rich sub-layer

IB

S2 z Fiber A

Fiber B

x

Fiber-mat-2

b

nB

P4 Fiber-mat-2

(ξB=1)

Δφ

PA(ξA) r_a

IA

Fig. 3. Layer arrangement of the RVE. The RVE consists of three fiber-mat layers while each fiber-mat layers consists of two fiber-rich sub-layers and three matrixrich sub-layers. The thickness of the fiber-rich and matrix-rich sub-layers are equal to 2b, and b/10, respectively, where b is the semi-minor axis of the ellipse.

IB

P1

P2 (ξA=1)

(ξA=0)

nA

y

P3

(3)

(4)

(5)

(6)

that a second order fiber orientation tensor with eigenvalues (0.5, 0.5, and 0) corresponds to an ideally transversely isotropic RVE [14]. The fiber is generated randomly in these two fiber-rich sublayers for each fiber-mat layer. Initially, the spine of a fiber is represented by a line segment with center point coordinates (x, y, and z) and orientation (u, h), where u e (0, 2p) being the in-plane angle and h being the out-of-plane angle, which is equal to zero in this paper. The z coordinate can be either 0 or t, with t being the sum of the thickness of a fiber-rich sub-layer and a matrix-rich sub-layer. The fiber can be straight or curved. Straight fibers reside entirely in either the top or bottom fiber-rich sub-layers. If a newly generated fiber potentially intersects existing fibers, the new one has to bend away from the existing fiber to the opposite sub-layer to avoid intersection. Thus, the fiber becomes a curved one. The curved fiber ‘‘interweaves” the straight fibers within the fiber-mat layer. The generation of a curved fiber implies, first, the creation of a straight line-segment (referred to as the generic line of the fiber) with random orientation and mid-point location. All fibers previously included in the configuration are projected on to the x–y plane and the minimum distance between the current fiber and all the previous ones is computed. The projection reduces the 3D geometric problem into a 2D one. Then, all the intersections on the generic lines are marked. Finally, the current generic line is ‘‘curved” by moving the intersecting points to the opposite sub-layer (i.e., by changing the z coordinate of the intersecting point). Two additional side points are added to each side of the shifted point, one on the top and one on the bottom sub-layer, to simulate the local bending of the fiber. This process generates the spine of the fiber as illustrated in Fig. 4a. Each 3D fiber is obtained by sweeping the cross-sectional profile a long the spine. Each straight fiber is approximated by a prism with two dodecagonal bases, while a curved fiber is approximated by a series of connected irregular polyhedra. In order to avoid interpenetration of fibers, a 3D intersection test is performed upon generation of a curved fiber. The distance between each polyhedral segments of the newly generated fiber and those of the existing ones is calculated and

x

PB(ξB)

(ξB=0)

Fig. 4. Schematic illustration of side (a) and top (b) views of intersecting fibers. (a) One fiber over-crossing another (side view). The generic lines of two fibers, A and B, originally intersect at point IB on the bottom layer. This point of fiber A was moved to the top sub-layer (noted as IA) and two side points S1 and S2 were also assigned to fiber A introduce a bend to avoid intersecting fiber B. (b) Relative position of two potentially intersecting fibers A and B at points IA and IB on the projection x–y plane (top view).

checked so as not to exceed the minimum separation value imposed. The mathematical implementation of the 2D and 3D fiber penetration tests is described in Section 2.2. The complete procedure for generating the RVE is outlined in the flowchart in Fig. 5. 2.2. Mathematical implementation In the following, the 2D fiber intersection will be discussed first, and second, the 3D intersection analysis algorithm will be presented. For the 2D analysis, we will refer to Fig. 4b, which illustrates the relative position of two fibers’ projection, A and B, on the x–y plane. The length of the fibers has been normalized to one. Given one end point (~ P E ) and the unit direction (~ n) of a fiber’s generic line segment, points on the line segment can be parameterized as

~ PðxÞ ¼ ~ P E þ nL~ n;

0 6 n 6 1;

ð1Þ

where L is the length of the fiber and the scalar n is the normalized distance measured along the segment from one end point. Thus, the distance between two points on two different generic line segments is given by

dðn1 ; n2 Þ ¼ k~ P 1 ðn1 Þ  ~ P 2 ðn2 Þk;

0 6 n1 6 1 and 0 6 n2 6 1:

ð2Þ

To determine the minimum distance between two fibers the constrained on linear two-variable function d(n1, n2) described in Eq. (2) needs to be minimized, i.e., solve

dmin 2D ¼ minðdðn1 ; n2 ÞÞ subject to the constraints 0 6 n1 6 1 and 0 6 n2 6 1.

ð3Þ

2796

Y. Pan et al. / Composites Science and Technology 68 (2008) 2792–2798

Error! Start X, φ, Layer

Separation? (2D computation) Yes

Adjust the fiber to a Curved one

No

Save data and update variables

No

Separation? (3D computation) Yes

Vf < Vf_obj Yes Done! Fig. 5. Flow chart of the RVE generation algorithm.

For the 3D analysis, we note first that a straight fiber is represented as a convex prism with two dodecagon end-faces, which a curve done consists of several convex irregular polyhedra. A convex polyhedron with s faces can be defined algebraically as the set of solutions to a system of linear in equalities:

A~ x 6 b;

ð4Þ

where A is a real s  3 matrix, b is a real s  1 vector and ~ x is a vector corresponding to the point coordinates in 3D space [39]. The Linear Matrix Inequality (LMI) in Eq. (4) can be computed from the vertices of the polyhedron. The distance between any two points from polyhedron A and polyhedron B is of the form:

dð~ x1 ; ~ x2 Þ ¼ k~ x1  ~ x2 k;

A1~ x1 6 b1 and A2~ x2 6 b2 :

ð5Þ

The minimum distance between two convex polyhedra can be determined by solving the convex optimization problem:

dmin 3D ¼ minðdð~ x1 ; ~ x2 ÞÞ

ð6Þ

subject to the constraints A1~ x1 6 b1 and A2~ x2 6 b2 . The numerical solution of the two convex optimization problems described in Eqs. (3) and (6) can be obtained by applying the existing algorithm described in next section. 2.3. Algorithm implementation The RVE generation algorithm is implemented in Matlab [40]. During the 2D geometry calculation, the constrained nonlinear optimization function fmincon from the Optimization Toolbox is called to solve the problem described in Eq. (3). In the 3D case, the CVX modeling system developed by Grant et al. [41] for Matlab is employed to formulate and solve the convex optimization problem described by Eq. (6). To choose the RVE size, we rely on results from previous studies for straight fiber RVEs with varying sizes [10] which suggest that L=‘ ¼ 2 yields satisfactory results, where L and ‘ are the length of RVE and length of fibers, respectively. Similar results were also reported by Kari et al. [15] for transversely

Fig. 6. RVE of the random chopped fiber-reinforced composite with curved fiber bundles in ABAQUS, fiber volume fraction of which is 35.1%: (a) RVE with curved fibers, (b) fiber orientation distribution, and (c) mesh of the RVE with 4-node tetrahedron elements.

randomly distributed short-fiber composites. In this work, we select L=‘ ¼ 2 as the RVE size. The RVE has dimensions 80b  80b  12.7b, where b is the semi-minor axis of the cross-sectional ellipse, and contains 174 straight and curved fibers with aspect ratio, ‘=2b ¼ 20. The cross section of a fiber is a dodecagon approximation of an ellipse with a/b = 2, where a is the major semi-axis. Following the protocol, a RVE with fiber volume fraction of 35.1% is generated as shown in Fig. 6a. The in-plane fiber orientation distribution is shown in Fig. 6b. 3. Finite element analysis results Finite element analysis (FEA) is employed to determine the composite material’s linear elastic response. The FE model is built

2797

Y. Pan et al. / Composites Science and Technology 68 (2008) 2792–2798

and solved using the commercial FEA package ABAQUS [42], in which the Matlab geometry data is imported via the Python scripting interface. It is assumed that the fiber and matrix are both linear elastic and isotropic and that they are perfectly bonded at their interface. Due to the in-plane random distribution of the fibers, the composite behaves, at the macroscopic scale, similar to a homogeneous, transversely isotropic material. The condition for equivalence is formulated as the equality of the strain energy in the two media subject to the same boundary conditions. The determination of the equivalent homogeneous material properties implies a homogenization procedure (for detail see [38]) that requires the application of six independent loading conditions on the analyzed RVE. Each loading case consists of specifying displacement fields that render null all but one of the six independent components of the strain tensor. Resolving the stress field in the heterogeneous material through a static equilibrium analysis allows for the calculation of the average stress field com i . Since the average strain ej is imposed, the stiffness tenponents, r  ij , can be directly sor components for the equivalent material, C calculated from the generalized Hooke’s law:

r i ¼ C ijej ; i; j ¼ 1; 2; . . . ; 6:

14:47

6 5:06 6 6 6 3:84  C¼6 6 0 6 6 4 0 0

5:06

3

3:84

0

0

0

14:22 3:87

0

0

3:87

9:49

0

0

0

0

5:09

0

0 0

0 0

0 0

0 7 7 7 0 7 7: 0 7 7 7 0 5

ð8Þ

3:06 0 3:10

We note that the computed values of the off-diagonal entries in the  are two or three orders 4th–6th rows and columns of the matrix C of magnitude smaller than the other entries and, thus, can be considered null for all practical purposes. To investigate the in-plane isotropy, we employ the parameter defined in Ref. [43]:

2C 44 aXY ¼  XY  ; C 11  C 12

FEM micro-mechanical analysis

Halpin–Tsai equations based empirical predictions [31]

E

(GPa)

11.83

12.56 G

(GPa)

5.09

4.48

ð7Þ

For the composite considered in this study, the fiber and matrix Young modulus and Poisson ratios are Ef = 70 GPa and mf = 0.2, and Em = 3 GPa and mm = 0.35, respectively (note that the values provided for E and m are generic values representative of the materials used). The meshing of the model is performed using ABAQUS’s built-in free meshing algorithm. The RVE geometry does not lend itself, due to its complexity, to the use of hexahedra elements. Thus the meshing is performed with 4-node tetrahedron elements (C3D4 in ABAQUS), with a total element count of 492,647, corresponding to 253,965 DOFs. The meshed RVE is shown in Fig. 6c. Following the application of the homogenization procedure, the resulting material elastic constants of the homogeneous material  ij are (in GPa) C

2

Table 1 Comparison of random-fiber composite in-plane elastic constants

ð9Þ

Fig. 7. Curved fiber formed by connected polyhedra.

obtained using the traditional equations for short random-fiber composites [36], based on the Halpin–Tsai estimations of the longitudinal and transverse moduli of the equivalent aligned short-fiber composite. The results are listed in Table 1. Note that this result is based on empirical relations and the model does not account for the curved fibers, thus justifying the higher value for Young’s modulus [9]. The higher value of the in-plane shear modulus for the FEA result may be due to the interweaving regions formed by straight and curved fibers, which are stiffer when subjected to in-plane shear loading. As mentioned before, the 3D volume of a fiber is generated by sweeping the cross-sectional profile along a fiber spine. In the case of curved fibers, the spine is formed by a series of connected line segments, with some of them oriented in out-of-plane directions. Since the cross-sectional profile always remains normal to the path during sweeping [42], the presence of the out-of-plane spine segments leads to the creation of sharp edges in the region of the two added side-points, as can be seen in Figs. 4(a) and 7. However, this is not true in the real sample, where a fiber smoothly bends away from another at the crossing region. The artificial sharp turns of the fiber generated from the sweeping about the connected line might introduce local stress concentrations which may influence the accuracy of the numerical solution to a loading simulation. Thus, it may be desirable to avoid such sharp ‘‘turns” in the fiber geometry. The influence of the local stress concentration on the overall RVE response is a subject of on going research, as is the interrogation of the results for similar ‘‘smooth” geometries.

where

C 11 þ C 22 C XY 11 ¼ 2

4. Conclusions

ð10Þ

which takes the value aXY = 1 for a transversely isotropic material whose symmetry axis coincides with the z axis. Substituting the en from Eq. (8) into Eqs. (9) and (10), we have aXY ¼ 91:2%, tries in C confirming that the behavior of the generated RVE is close to that of a transversely isotropic material. Results validation is performed by comparing Young’s and Shear moduli for the equivalent homogeneous material, computed  ij matrix obtained from previous section, against those using the C

By intermixing curved and straight fibers, the fiber volume fraction of the RVE can be increased drastically, comparable to values achieved in industrial-grade materials. This micro-geometry generation technique allows for very detailed simulations of the response of high volume fraction RaFCs, which can capture and quantify effects at the micro-scale level. The representativeness of the geometry generated using the described approach was investigated via a FE analysis based homogenization approach that allows the characterization of the macro-

2798

Y. Pan et al. / Composites Science and Technology 68 (2008) 2792–2798

scale composite response. The finite element analysis results indicate that the generated RVE satisfactorily approximates the elastic transversely isotropic behavior of mostly in-plane random chopped fiber-reinforced composites. These first results demonstrate that micro-scale geometries generated with the presented algorithm can be employed in in-depth studies such as the characterization of RaFC response beyond the elastic range or their damage initiation and propagation behavior. Acknowledgements This work was funded by NSF through the CMS-0409282 Grant and partially supported by the Department of Energy Cooperative agreement No. DE-FC05-950R22363. Such support does not constitute an endorsement by the Department of Energy of the views expressed herein. The authors would like to gratefully acknowledge the support of NSF Program Manager Dr. Ken Chong, as well as the insights and assistance of Drs. Libby Berger and Stephen Harris of GM. References [1] Jacob GC, Starbuck JM, Fellers JF, Simunovic S, Boeman RG. Fracture toughness in tandom-chopped fiber-reinforced composites and their strain rate dependence. J Appl Polym Sci 2006;100:695–701. [2] Jacob GC, Starbuck JM, Fellers JF, Simunovic S, Beoman RG. Crashworthiness of various rancom chopped carbon fiber reinforced epoxy composite materials and their strain rate dependence. J Appl Polym Sci 2006;101:1477–86. [3] Corum JM, Battiste RL, Ruggles MB, Ren W. Durability-based design criteria for a chopped-glass-fiber automotive structural composite. Compos Sci Technol 2001;61:1083–95. [4] Dahl JS, Smith GL, Houston DQ. The influence of fiber tow size on the performance of chopped carbon fiber reinforced composites. In: materials and processing technologies for revolutionary applications, 37th ISTC conference proceedings. Seattle, Washington: 2005. [5] Daniel IM, Ishai O. Engineering mechanics of composite materials. Oxford University Press; 1994. [6] Eshelby JD. The determining of the elastic field of an ellipsoidal inclusion and related problem. Proc Roy Soc Lond Ser A 1957;241:376–96. [7] Mori T, Tanaka K. Average stress in matrix and average elastic energy of materials with misfitting inclusions. Acta Metall Mater 1973;21:271–4. [8] Advani SG, Tucker CL. The use of tensors to describe and predict fiber orientation in short fiber composites. J Rheol 1987;31:751–84. [9] Ionita A, Weitsman WJ. On the mechanical response of randomly reinforced chopped-fibers composites: data and model. Compos Sci Technol 2006;66:2566–79. [10] Iorga L, Pelegri AA, Pan Y. Numerical characterization of material elastic properties for random fiber composites, in press. [11] Garnich MR, Karami G. Finite element micromechanics for stiffness and strength of wavy fiber composites. J Compos Mater 2004;38(4):273–92. [12] Meraghni F, Desrumaux F, Benzeggagh MF. Implementation of a constitutive micromechanical model for damage analysis in glass mat reinforced composite structures. Compos Sci Technol 2002;62:2087–97. [13] Gusev AA. Representative volume element size for elastic composites: a numerical study. J Mech Phys Sol 1997;45:1449–59. [14] Gusev AA, Heggli M, Lusti HR, Hine PJ. Orientation averaging for stiffness and thermal expansion of short fiber composites. Adv Eng Mat 2002;4:931–3. [15] Kari S, Berger H, Gabbert U. Numerical evaluation of effective material properties of randomly distributed short cylindrical fibre composites. Comput Mat Sci 2007;39:198–204. [16] Böhm HJ, Eckschlager A, Han W. Multi-inclusion unit cell models for metal matrix composites with randomly oriented discontinuous reinforcements. Comput Mat Sci 2002;25:42–53.

[17] Faessel M, Delisée C, Bos F, Castéra P. 3D modelling of random cellulosic fibrous networks based on X-ray tomography and image analysis. Compos Sci Technol 2005;65:1931–40. [18] Duschlbauer D, Böhm HJ, Pettermann HE. Computational simulation of composites reinforced by planar random fibers: homogenization and localization by unit cell and mean field approaches. J Compos Mater 2006;40:2217–34. [19] Eason TG, Ochoa OO. Material behavior of structural reaction injection molded composites under thermomechanical loading. J Compos Mater 2001;34:432. [20] Hill R. Elastic properties of reinforce solid: some theoretical principles. J Mech Phys Sol 1963;11:357–72. [21] Drugan WJ, Willis JR. A mircomechanics-based nonlocal constitutive equation and estimates of representative volume element size for elastic composites. J Mech Phys Sol 1996;44(4):497–524. [22] Trias D, Costa J, Turon A, Hurtado JF. Determination of the critical size of a statistical representative volume element (SRVE) for carbon reinforced polymers. Acta Mater 2006;54:3471–84. [23] Ivanov I, Tabiei A. Three-dimensional computational micro-mechanical model for woven fabric composite. Compos Struct 2001;54:489–96. [24] Quek SC, Waas A, Shahwan KW, Agaram V. Compressive response and failure of braided textile composite: part 2 – computations. Int J Non-Linear Mech 2004;39:649–63. [25] Caizzo AA, Constanzo F. On the constitutive relations of materials with evolving microstructure due to microcracking. Int J Sol Struct 2000;37:3375–98. [26] Tu ST, Cai WZ, Yin Y, Ling X. Numerical simulation of saturation behavior of physical properties in composites with randomly distributed second-phase. J Compos Mater 2005;39(7):617–31. [27] Torquato S. Random heterogeneous materials: microstructure and macroscopic properties. In: Antman SS, Marsden JE, Sirovich L, Wiggins S, editors. Interdisciplinary applied mathematics, vol. 16. NewYork: Springer; 2002. [28] Williams SR, Philipse AP. Random packings of spheres and spherocylinders simulated by mechanical contraction. Phys Rev E 2003;67. 051301-1-9. [29] Widom B. Random sequential addition of hard sphere to a volume. J Chem Phys 1966;44(10):3888–94. [30] Pan Y, Iorga L, Pelegri AA. Analysis of 3D random chopped fiber reinforced composites using FEM and random sequential adsorption. Comput Mat Sci 2008. doi:10.1016/j.commatsci.2007.12.016. [31] Gusev AA, Hine PJ, Ward IM. Fiber packing and elastic properties of a transversely random unidirectional glass/epoxy composite. Compos Sci Technol 2000;60:535–41. [32] Berger L. General motors. Personal communication; 2007. [33] Evans KE, Gibson AG. Prediction of the maximum packing fraction achievable in randomly orientated short-fibre composites. Compos Sci Technol 1986;25:149–62. [34] Toll S. Packing mechanics of fiber reinforcements. Polym Eng Sci 1998;38(8):1337–50. [35] Parkhouse JG, Kelly A. The random packing of fibers in three dimensions. Proc Math Phys Sci 1995;451(1943):737–46. [36] Agarwal BD, Broutman LJ. Analysis and performance of fiber composites. New York: John Wiley & Sons; 1990. [37] Coelho D, Thovert JF, Adler PM. Geometrical and transport properties of random packings of spheres and asperical particles. Phys Rev E 1997;55(2):1959–78. [38] Pan Y, Iorga L, Pelegri AA. Analysis of 3D random chopped fiber reinforced composites using FEM and random sequential adsorption. Comput Mater Sci 2007. [39] Weisstein EW. Convex polyhedron. From mathworld – a Wolfram web resource; 2002, . [40] The Mathworks, Matlab, The Language of Technical Computing, Version 7.5; 2007. [41] Grant M, Boyd S, Ye Y. CVX users’ guide, Version. 1.1; 2007. [42] ABAQUS/Standard, ABAQUS analysis user’s manual, Version 6.6. RI: ABAQUS Inc.; 2006. [43] Kanit T, N’Guyen F, Forest S, Jeulin D, Reed M, Singleton S. Apparent and effective physical properties of heterogeneous materials: representativity of samples of two materials from food industry. Comput Methods Appl Mech Eng 2006;195:3960–82.