Jtam-A4.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 51, 4, pp. 883-890, Warsaw 2013 NUMERICAL HOMOGENIZATION OF FIBER-REINFORCED COMPOSITES WITH COMPLEX MICROSTRUCTURAL FEATURES Marek Romanowicz Bialystok University of Technology, Department of Mechanical Engineering, Białystok, Poland e-mail address: m.romanowicz@pb.edu.pl Amicromechanicalmodel for solution of the problemof transverse stiffness of unidirectional fiber reinforced composites taking into account the presence of the interphase and a random distribution of fibers over the transverse cross-section has been developed. The influence of interphase properties on the elasticmodulus in the transverse direction of glass/epoxy com- posites has been investigated. It has been found that numerical computations of transverse Young’s modulus obtained from the unit cell models for different fiber volume fractions are in a good agreement with experimental data. Key words: fiber-reinforced composites, micro-mechanics, finite element analysis 1. Introduction The third phase between original constituents of a composite, known as an interphase, forms as a result of irreversible chemical reaction in the processing of thermosetting-matrix composites. It has been observed by using the atomic force microscope nanoindentation that the interphase has spatially varying rather than uniform properties (Gao and Mader, 2002). It is difficult to investigate experimentally the effect of such thin interphases between fibers and a matrix on composite properties. An alternative, or at least a complement to experimentation, is making use of numericalmicromechanicalmodelswhich take into account the presence of the interphase. Application of a representative unit cell approach in the modeling of the interphase is based on the assumption that it is possible to isolate a small volume element of a composite that is statistically representative of thematerial as a whole. Such an approach for themodeling of the elastic behavior of fiber reinforced composites with interphases was proposed by Wacker et al. (1998), Young and Pitchumani (2004) andWang et al. (2006). In the present study, both the inhomogeneous properties of the interphase and the random fiber distribution over the transverse cross-section are incorporated into a onemicromechanical model. Application of such a unit cell in the modeling of the mechanical properties of fiber reinforced composites is novel. It can be expected that the consideration of these two aspects simultaneously allows us to identify the composite properties accurately. 2. Description of the micromechanical model 2.1. Unit cell No manufacturing process can guarantee perfect uniformity of fiber placement in the com- posite volume. Therefore, the analysis of composites with a randommicrostructure is desirable. In this paper, the unit cell models of a randomly distributedfiber composite are generated using Wongsto and Li’s algorithm (Wongsto and Li, 2005). The generation of a microstructure in this algorithm is based on disturbing an initially hexagonal periodic array of fibers packed in 884 M. Romanowicz a two-dimensional plane. The main input data are the number of fibers N, the fiber radius rf and the fiber volume fraction Vf. A random angle α between 0 ◦ and 360◦ is first generated to determine the direction for the fiber to shift. Along this direction, amaximumpossible distance of the shift d is found. Next, the second random number k is generated between 0 and 1, and the actual distance of shift is calculated as kd. This process is performed for every fiber to complete one iteration. The stirring of fibers terminates after the required number of iterations is completed and when the fiber volume fraction Vf in the central region of the unit cell is the same as in the outskirts region. Wongsto and Li’s algorithm was implemented in Free Pascal programming language. Figure 1 shows a typical unit cell with Vf =0.52 containing N =105 fibers. Using an ANSYS finite element code, plain strain models made of isoparametric trian- gular finite elements with 6 nodes (PLANE183) were generated. It was assumed that perfect bonding between the matrix and the fibers exists. To ensure accurate displacement and stress field representation within each unit cell, sufficiently dense meshes consisting of approximately 2 ·106 elements with a quadratic displacement field approximation were used. Themeshes com- prised 400 elements equally spaced around the circumference of each fiber. It was found that such a discretization was fine enough to give accurate results from the numerical viewpoint and any increased mesh density provided no further accuracy in the solution obtained. Fig. 1. Typical unit cell model with N =105 fibers distributed at random for the fiber volume fraction Vf =0.52 2.2. Interphase model Following the approach by Anifantis (2000), it has been assumed that elastic properties of the interphase are non-uniform across its thickness and vary exponentially. The empirical relationships between the radial distance from the fiber boundary and elastic moduli of the interphase Ei, νi, are Ei(r)=Em [ 1+ ( p Ef Em −1 ) 1− r ri e 1− r ri 1− rf ri e 1− rf ri ] νi(r)= νm [ 1+ ( q νf νm −1 ) 1− r ri e 1− r ri 1− rf ri e 1− rf ri ] (2.1) where r denotes the radial component of the polar coordinate system whose origin is placed in the centre of a fiber, ri and rf are radii of the interphase and fiber, respectively. The factor q can be approximated by a function of the factor p, which, in turn, determines the bonding Numerical homogenization of fiber-reinforced composites ... 885 conditions on the fiber surface and is defined by the ratio between the interphase modulus at r= rf and the fiber modulus q= 1− νf νm ( Ef Em −1 ) νf νm 1−p p +1 (2.2) where p= Ei(r= rf) Ef (2.3) where 0 < p < 1. Note, that a smaller p leads to a lower stiffness of the interphase. In nu- merical implementation, the radially graded interphase is discretized into 10 concentric layers characterized by homogeneous elastic properties as shown in Fig. 2. Fig. 2. Close-up view of the finite element mesh in the region of individual consistuents 2.3. Model parameters Numerical calculations have been performed for 3 different thicknesses of the interphase t = 300, 500 and 700nm and 4 different fractions of the fiber volume Vf = 0.333, 0.436, 0.531, 0.631. In this paper, 5 different non-uniform distributions of the elastic moduli within the interphase characterized by factors p=0, 0.15, 0.30, 0.50, 1.00 have been assumed for each value of t and Vf. Note, that when p= 0, a two-phase model is considered and no interphase occurs. And, in turn, when p=1, a three-phasemodel without jumps in both elastic properties on the fiber surface is considered. Finite element computations have been carried out for a glass/epoxy composite reported by Wacker et al. (1998). Thematerial parameters for the glass fibers sizedbyusing thepolyurethanebinder (PU)andanhydride-curedepoxy resin (LY556/917) are summarized in Table 1. Table 1.Properties of the constituent materials Consitituent Properties glass fiber (E-glass) Ef =77GPa, νf =0.2, rf =6µm epoxy resin (556/917) Em =3.11GPa, νm =0.34 886 M. Romanowicz 2.4. Homogenization According to Hill’s principle (Hill, 1963), the volume average of the strain energy in a unit cell should be equal to the energy produced by the average strain and stress in an equivalent material U =Ueq (2.4) where the strain energy of the equivalent material is given by Ueq = 1 2 ∫ V CijklΦijΦkl dV = 1 2 CijklΦijΦklV (2.5) and the strain energy stored in the unit cell is as follows U = 1 2 ∫ V σijεij dV (2.6) where σij, εij are components of the microstress and microstrain, respectively, Φij are compo- nents of themacrostrain and Cijkl are components of the stiffness tensor. The elastic constants ν23 and E2 are determined by applying normal components of the macrostrain to the unit cell, i.e. first Φ22 (state I), and next simultaneously Φ22 and Φ33 (state II). By using Eq. (2.4), an appropriate system of linear algebraic equations for the unknown stiffnesses can be formed. Under plane strain conditions, the solution for the stiffnesses is C2222 = 2UI V (ΦI22) 2 C2233 = UII(ΦI22) 2 −UII[(ΦII22) 2 − (ΦII33) 2] V (ΦI22) 2ΦII22Φ II 33 (2.7) The relationships between the stiffnesses and the elastic constants are as follows ν23 = C2233 C2222+C2233 E2 =C2222 (1+ν23)(1−2ν23) 1−ν23 (2.8) 2.5. Boundary conditions In the numerical homogenization, the periodic boundary conditionsmust be imposed on the unit cell to reflect the repeatability of the microstructure and to ensure the compatibility of the displacement fields. On the assumption of periodicity, each displacement field ui may be decomposed in a linear part associated with themacrostrain Φij and a periodic one u p i ui(x1,x2)=Φijxj +u p i(x1,x2) (2.9) These relations are implemented at each periodic pair of nodes to link the displacements of the top and the bottomboundaries and the displacements of the right and left boundaries of the unit cell. Because of a huge number of nodes at the opposite boundary surfaces, a APDL program (Ansys, 2007) has been used to generate all required constraint conditions (2.9) automatically. 3. Validation of random fiber distribution An efficientmethod of characterizing the regularity or irregularity in the distributions generated by this code is based on the computing of the pair distribution function g(r) (Pyrz, 1994). This statistical function is defined as the probability of finding the centre of a fiber inside an annulus Numerical homogenization of fiber-reinforced composites ... 887 of internal radius r and thickness dr with the centre at a randomly selected fiber. It can be expressed by the following equation g(r)= 1 2πr dK(r) dr (3.1) where K(r) is the second order intensity function. It gives the average number of fiber centroids expected to lie within the radius r about an arbitrary fiber centroid. For observations within the finite window of area A, the second order intensity function is defined as (Pyrz, 1994) K(r)= A N2 N ∑ k=1 Ik(r) Rp (3.2) where N is the number of fibers within A, Ik(r) is the number of fiber centroids within a circle of the radius r not counting the central fiber, and Rp is the ratio of the circumference of a circle of the radius r inside the viewingarea to the entire circumference. For regular fiberdistributions, g(r) takes very sharp peaks and strongly oscillates with r. For random fiber distributions, on the other hand, g(r) approaches a constant value of 1 as r increases. 4. Results and discussion In Figs. 3 and 4, comparative analyzes of g(r) in dependence on the non-dimensional radial distance are presented for two fiber volume fractions about 40% and 60%.The results presented in these figures are obtained by the averaging over 10 samples containing N =39, 105 and 538 fibers. The comparison reveals that the code used in the present paper is, in general, able to produce random fiber distributions. Due to edge effects, the condition that g(r) = 1 as r increases cannot be fully satisfied for small N. Although, the fiber arrangements for small Vf and large N aremore homogeneous and statistically valid, a tendency toward the random fiber distribution is visible for all tested samples. Therefore, taking the above into account, the unit cells involving the number N =105 have been chosen for further study. Fig. 3. Variation of the pair distribution function g(r) with normalized radius r/rf for the fiber volume fraction about 40%; (a) random fiber arrangement; (b) initial hexagonal fiber arrangement 888 M. Romanowicz Fig. 4. Variation of the pair distribution function g(r) with normalized radius r/rf for the fiber volume fraction about 60%; (a) random fiber arrangement; (b) initial hexagonal fiber arrangement Thenumerical calculations of transverseYoung’smodulus in dependence on the fiber volume fraction for the glass/epoxy composite with three values of the thickness of the interphase are presented in Figs. 5-7. In order to validate the numerical results, transverse Young’s moduli obtained through the unit cell models have been compared to the experimental data reported byWacker et al. (1998). It can be seen from these figures that the interphase parameters t and p have the significant effect on transverse Young’s modulus. Note, that unit cell models without the interphase (p = 0) always produce unreliable results that lie under experimental points. For the case of models of too thin interphase (Fig. 5), numerical points have been also found to localize under experimental points even for an extremely stiff interphase (p= 1). Models with the interphase of amoderate thickness (Fig. 6) guarantee thatmean values of the experimental moduli liewithin the range of p from0 to 1. For the case ofmodels of a thick interphase (Fig. 7), not only themean values but also standard deviations of themeans are within this scope. Thus, the interphase with t= 700nm is the most reliable. For the thickness t= 700nm, the models with the interphase parameters p=0.15 and 0.30 are in a very good agreement with themean values of experimental data. Fig. 5. Comparison between the model predictions for t=300nm and the experimental data (Wacker et al., 1998) Numerical homogenization of fiber-reinforced composites ... 889 Fig. 6. Comparison between the model predictions for t=500nm and the experimental data (Wacker et al., 1998) Fig. 7. Comparison between the model predictions for t=700nm and the experimental data (Wacker et al., 1998) 5. Concluding remarks By using the unit cell approach and the finite element method, a micromechanical model of fiber-reinforced composite for predicting transverse Young’s modulus has been proposed. In particular, the effect of complex features such as the presence of the third phase between the reinforcement and the matrix and its inhomogeneous as well as the random distribution of the fibers over the transverse cross-section have been incorporated into the model. Also, a method of identifying the interphase properties has been presented and validated by comparison with the experimental data. Acknowledgments The financial support of the National Science Centre of Poland under contract DEC- 2011/03/D/ST8/04817 is thankfully acknowledged. References 1. AnifantisN.K., 2000,Micromechanical stress analysis of closely packedfibrous composites,Com- posites Science and Technology, 60, 1241-1248 2. ANSYS 11, 2007,ANSYS Users’s Manual, Cannonsburg, PA 890 M. Romanowicz 3. Gao S.L.,Mader E., 2002, Characterization of interphase nanoscale property variations in glass fibre reinforced polypropylene and epoxy resin composites,Composites Part A, 33, 559-576 4. Hill R., 1963, Elastic properties of reinforced solids: some theoretical principles, Journal of the Mechanics and Physics of Solids, 11, 357-372 5. Pyrz R., 1994, Correlation of microstructure variability and local stress field in two-phase mate- rials,Materials Science and Engineering A, 177, 253-259 6. Wacker G., Bledzki A.K., Chate A., 1998, Effect of interphase on the transverse Young’s modulus of glass/epoxy composites,Composites Part A, 29, 619-626 7. Wang J., Crouch S.L.,Mogilevskaya S.G., 2006,Numerical modeling of the elastic behavior offiber-reinforcedcompositeswith inhomogeneous interphases,Composites Science andTechnology, 66, 1-18 8. Wongsto A., Li S., 2005, Micromechanical FE analysis of UD fibre-reinforced composites with fibres distributed at random over the transverse cross-section,Composites Part A, 36, 1246-1266 9. YangF.,PitchumaniR., 2004,Effects of interphase formationon themodulus and stress concen- tration factor of fiber-reinforced thermosetting-matrix composites, Composites Science and Tech- nology, 64, 1437-1452 Manuscript received November 28, 2012; accepted for print February 14, 2013