A finite element based orientation averaging method for predicting elastic properties of short fiber reinforced composites

Short fiber reinforced composites have a variety of micro-structural parameters that affect their macromechanical performance. A modeling methodology, capable of accommodating a broad range of these parameters, is desirable. This paper describes a micro-mechanical model which is developed using Finite Element Analysis and Orientation Averaging. The model is applicable to short fiber reinforced composites with a wide variety of micro-structural parameters such as arbitrary fiber volume fractions, fiber aspect ratios and fiber orientation distributions. In addition to the Voigt and Reuss assumptions, an interaction model is developed based on the self-consistent assumption. Comparisons with experimental results, and direct numerical simulations of Representative Volume Elements show the capability of the model for fair predictions.


Introduction
Short Fiber Reinforced Composites (SFRCs) are becoming more and more interesting for different industries due to their high strength-todensity and stiffness-to-density ratios in comparison to unreinforced polymers. Besides, the manufacturing process of these materials are quick and low cost [1,2]. In order to obtain efficient and optimized designs, the ability to predict the behavior of SFRCs in a quantitative manner is crucial. To do so, and more specifically, to capture the effect of a wide variety of micro-structural parameters of SFRCs affecting their macro-mechanical behavior, it is necessary to use micro-mechanics based models. As a result, a large number of studies have been devoted to this topic (see e.g. Ref. [3][4][5][6][7][8]).
A micro-mechanical modeling approach which has been frequently used for SFRCs is computational homogenization (see e.g. Ref. [9][10][11][12]). In this approach, a numerical Representative Volume Element (RVE) of an SFRC is analyzed by numerical methods (most often by the Finite Element Method) and the homogenized material response is obtained by volume averaging. This modeling approach has a very strong predictive capability. Nevertheless, it is not always feasible to use this method for SFRCs mainly due to: computationally expensive simulations, and more importantly, difficult RVE generation [13,14]. Mirkhalaf et al. [15] found out that generating realistic RVEs for SFRCs with high fiber volume fractions and particularly, high fiber aspect ratios is very challenging.
An alternative approach to computational homogenization for modeling SFRCs is two-step homogenization techniques (see e.g. Ref. [6,16,17]). In this modeling approach, the first step typically uses a mean-field model such as Mori-Tanaka [16,17], and in the second step, the fiber orientation distribution is taken into account using an interaction rule such as the Voigt assumption. Within a two-step homogenization approach, another possibility is to use Finite Element Analysis in the first step [6], and performing Orientation Averaging (OA) in the second step. Modniks and Andersons [6] used this method to estimate anisotropic elastic properties of short flax fiber reinforced Polypropylene. A Finite Element model was developed to predict the elastic properties of a Unit Cell (herein considered as a single fiber embedded in the matrix material), still respecting the respective volume fractions of fibers and matrix. The elastic properties of a randomly distributed SFRC was then obtained using the OA approach assuming Voigt interaction.
In this study, we present a two-step homogenization approach for SFRCs, using Finite Element Analysis and Orientation Averaging. FE calculations of a UC are used to obtain UC homogenized properties. Using FEA not only gives accurate UC properties, but also provides the opportunity of including further phenomena such as matrix inelasticity and fiber-matrix debonding. For the orientation averaging phase (second step), in addition to the Voigt assumption (upper bound), a Reuss interaction is developed to obtain the lower bound of stiffness, as well. Simulations in this study show that they are strong assumptions, and not always resulting in accurate predictions. Thus, the self-consistent assumption was used to develop a novel interaction model as an inbetween approach. Comparisons to computational homogenization analyses on realistic RVEs (performed in this study) show that for some cases, the self-consistent model provides a more accurate prediction of the composite stiffness properties than the Voigt and Reuss models. The presented method in this paper is applicable to almost any SFRC, with any desired fiber aspect ratio, fiber volume fraction and fiber orientation distribution. This is an important advantage over computational homogenization, since RVE generation is, for some cases, very challenging.
The remainder of this paper is structured as follows. Section 2 describes the modeling approach and the Voigt and Reuss interactions. Section 3 illustrates the implementation of the model and FE calculations required for elastic predictions. Some initial results, and comparison to experiments (adopted from literature), are presented in Section 4. Section 5 describes the interaction model developed based on the selfconsistent assumption, and its implementation. Section 6 gives final results and comparisons to experiments, and RVE simulations. Finally, Section 7 summarizes the contribution of this study and gives some concluding remarks.

Orientation averaging
In this section, the micro-mechanical model, developed based on an orientation averaging method, is illustrated. First, the homogenized properties of a UC (including a single fiber) are obtained. Different approaches could be used for that purpose (see Ref. [18] for an overview). In this study, FE simulations are used due to two main reasons, first: to have a more accurate description of a UC, and second: to be able to incorporate other phenomena such as inelasticity and fiber-matrix debonding in future studies. Once the UC homogenized properties are known, the mechanical response of the studied SFRC is calculated based on the orientation distribution of the fibres. Fig. 1 shows schematically different phases of the method.
Two configurations are considered: one at the composite level (global), and one at the UC level (local). Fig. 2 shows a schematic representation of a Unit Cell and the two aforementioned configurations.
If R is the rotation from the composite configuration to the UC configuration, we have: where, e i and e L i indicate the global and local orthonormal base vectors, respectively. A convenient and practical way of parameterizing the orientation of a fibre in a 3D space is using two angles [19], as it is shown in Fig. 3.
In Fig. 3, p is a unit vector representing the fibre orientation. Then, the rotation tensor (represented in matrix format) is obtained as In order to obtain the rotation tensor, first, the global configuration is rotated by an angle of θ around axis-2 and then, it is rotated by an angle of φ around axis-3 (considering fiber extension in 3-direction). The UC stress is obtained by where, σ L U is the UC stress at the local configuration. We note that Equation (3) represents the change from local to global coordinates. However, in order to formally use the notation of tensors (as opposed to matrices), we identify each coordinate system as a configuration. The symbol ⊗ represents a non-standard open product. 1 By weighted integration over the unit sphere, the volume averaged composite stress becomes:    3. Angles to describe a fibre orientation in a 3D configuration. 1 The index notation for the non-standard open product (⊗) is given where σ C is the composite stress, and ψ{p} is the probability distribution function of the orientation [19]. We define the integration over the unit sphere Ω as ∫ and note that distribution function (ψ{p}) has to be normalized, i.e., ∫ In the same fashion, we may introduce the transformation of the UC strain: where ε L U is the UC strain in the local configuration. The composite strain is then obtained as the average of UC strains: In the elastic regime (where both fibre and matrix are described by the Hooke's law), the UC stress at the local configuration is given by In order to determine ε L U and σ L U for each UC (pertaining to each direction p), we need to introduce a modeling assumption for the local global interaction. Below, we introduce the two limiting assumptions of Voigt and Reuss.

Voigt assumption
The Voigt interaction can be referred to as the uniform strain assumption. In other words, using this interaction implies full kinematic compatibility between the interacting components (assuming they are in parallel). Using Voigt interaction could be interpreted as a global isostrain situation. Thus, we assume the uniform strain ε C throughout the composite, and the local strain of each UC is obtained as Using Equation (10) together with the constitutive relation (9) in relation (4) results in the composite stress: Where the composite stiffness using the Voigt assumption is identified as where the superscript V stands for the Voigt assumption.

Reuss assumption
The Reuss interaction is in fact a uniform stress assumption. This assumption implies that all UCs have the same global stress state (UCs connected in series), where σ C is constant, and Inserting Equation (9) with the transformed composite stress in Equation (13) into Equation (8) gives the flexibility relation: where the composite stiffness using the Reuss assumption is identified as where the superscript R stands for the Reuss assumption.

Implementation and unit cell FE calculations
In this section, the implementation of the OA approach, and FE calculations for the elastic predictions are described.

Voigt interaction
For an actual case of finite number of fibres, the integral (in relation (12)) converts to a summation over the number of fibres or Unit Cells: where N f refers to the number of fibres. It will be explained in Section 3.3 how to obtain the UC stiffness (C L U ) by FE simulations.

Reuss interaction
In case of Reuss interaction assumption (Equation (15)), the composite stiffness is obtained as

FE calculations for elastic predictions
To compute the composite stiffness (Equation (16) and (17)), it is needed to calculate the stiffness of a UC (C L U ). FE analyses are performed for that purpose. To obtain the UC dimensions, the same distance is considered from the fiber to all sides of the UC. This is schematically shown in Fig. 4.
Thus, knowing the fiber dimensions and the fiber volume fraction, the UC dimensions are obtained. The UC with a unidirectional fibre is assumed to be transversely isotropic, for which the matrix representation of the compliance tensor is given by (note that direction 3 is the fiber direction and plane 12 is the isotropic plane):

Fig. 4.
Obtaining the UC dimensions by assuming the same distance from the fiber to all sides of a UC.
The Voigt notation has been used where the stress and strain components are ordered as σ T = [σ 11 σ 22 σ 33 σ 12 σ 13 σ 23 ] ε T = [ε 11 ε 22 ε 33 2ε 12 2ε 13 2ε 23 ]. (19) Due to the symmetry of the compliance matrix, and considering the assumption of UC transverse isotropy, there are 5 independent elastic components to calculate the UC stiffness (C L U ). To obtain these elastic properties from FEA, it is needed to conduct three simulations under three different loading conditions. Application of a uniaxial stress loading in the fiber direction (see Fig. 5(a)), the Young's modulus in the fiber direction (E 33 ) and the Poisson's ratios ν 31 , ν 32 (which are equal due to the UC symmetry in 1 and 2 directions) are obtained. Uniaxial stress on transverse direction (see Fig. 5(b)), the transverse Young's modulus (E 11 ) and the Poisson's ratio ν 12 are obtained. Finally, for the shear modulus in planes parallel to the fiber direction (G 13 or G 23 ) a shear loading parallel to the fiber direction is needed (see Fig. 5(c)).
Abaqus was used to perform FE simulations on UCs. In order to enforce Periodic Boundary Conditions in Abaqus simulations, a plugin developed by Omairey et al. [20] is used.
Remark 1: By packing UCs in a periodic structure, the transverse directions (all possible directions in a plane perpendicular to the fiber direction) do not have the exact same properties. This is because the distances between the fibers (in a plane perpendicular to the fiber) in different directions are different. Hence, the UC structure used in this study (see Fig. 2(a)) is not a perfectly transversely isotropic material. Nevertheless, since deviations from a perfect transverse isotropy are negligible, it is a reasonable assumption.
Remark 2: To check the validity of the suggested UC ( Fig. 2(a)), unidirectional RVEs were generated for a polyamide/glass composite, and the homogenized properties of the RVEs and the UC were compared. Spatial discretization of these RVEs are shown in Fig. 6. Homogenized elastic properties of these RVEs and the UC are given in Table 1. It is seen that there are good agreements between RVEs and UC homogenized properties, and thus, it is reasonable to use the suggested UC structure.
Remark 3: Mesh sensitivity analyses were performed to make sure about the convergence of the behavior of UCs. Fig. 7 shows the transverse cross section of three FE discretizations of a UC for a polypropylene/flax SFRC with 13% of fiber volume fraction (this SFRC is modeled in Section 6). The mesh size in the longitudinal direction is similar to the cross section mesh size. The elastic properties of the UC with different FE meshes are shown in Table 2.
Remark 4. The number of fibres (or equivalently the number of UCs) is important since it should be a representative number. For each of the simulations performed in this study, the number of fibres (N f ) has been increased until a convergence in the macroscopic behavior is obtained. Table 3 gives the values of the Young's modulus obtained for a magnesium/carbon SFRC with 10% of Fiber volume fraction (this composite is modeled in Section 4) considering different number of UCs in the simulations, and assuming Voigt interaction.

Initial results
In this section, different SFRCs are modeled and the results are compared against experimental results and computational homogenization performed on Representative Volume Elements (RVEs). Information about these composites are summarized in Table 4.

A polyamide/glass composite
Elastic stiffness of an SFRC made from Polyamide 6.6 (PA 6.6) matrix reinforced with V f = 10% of short glass fiber is obtained and compared to experiments taken form [1]. The fiber length and diameter are l f = 240μm and d f = 10μm which result in an aspect ratio of λ = 24. Fibers have a planar distribution with a preferred direction which is caused by the injection molding fabrication process. Fig. 8(a) shows the preferentially planar orientation distribution of fibers in the composite. Kammoun et al. [1] cut samples from the injection molded plate with different angles with respect to the Injection Flow Direction (IFD). This is schematically shown in Fig. 8(b). Fig. 9 shows a comparison between stiffness prediction of the model and experimental results.

A magnesium/carbon composite
A composite made from AZ91D magnesium alloy matrix and T300 short carbon fibers is analyzed in this section. Fiber volume fraction is equal to 10%. Fibers are randomly distributed, and they have a length of 105 μm and a diameter of 7 μm which lead to an aspect ratio of 15.
Experimental results are taken from Ref. [21] and compared to the predictions obtained by the OA method (A random distribution of fibers is considered). Also, computational homogenization is performed on a realistic micro-structural sample and homogenized elastic properties are obtained. Digimat-FE was used for RVE generation and the pertinent FE analysis. Fig. 10 shows an RVE (with 3D randomly distributed fibers) and its spatial discretization.
It was very challenging to generate an RVE which fulfill both of intended fiber orientation distribution and fiber volume fraction. As a result, several attempts were required to reach that goal and create such an RVE. Hence, results obtained for only one realization is presented here.
In order to have a representative micro-structural sample, it is needed to obtain a representative size for the sample. Mirkhalaf et al. [22] proposed an approach for determining the required RVE size for heterogeneous materials at finite strains. For short fiber composites, different values are suggested by different authors (see e.g. Refs. [ [23][24][25]). We chose the RVE size as L RVE = 210 μm which is two time the fiber length. A comparison of the obtained results in the simulations and experiments is presented in Table 5.

Polyamide/glass composites
In this section, polyamide matrix SFRCs reinforced with randomly distributed short glass fibers are considered, with two fiber volume fractions, namely 7% and 10%. Using the chosen software, generating and simulating RVEs with higher fiber volume fractions were found to be very difficult and thus, a maximum of 10% was considered. Fig. 11 depicts an RVE of this composite with 10% of fiber volume fraction, distributed randomly. The length of the RVE was considered L RVE = 240 μm which is the same as the fiber length. A few efforts were given to generate bigger samples but they failed at either meshing or solution phases.. The obtained elastic properties are shown in Fig. 12. It is seen that RVE results for this composite are between the Voigt and Reuss bounds. In the next section, an intermediate interaction model (between the upper and lower bounds) will be presented based on the self-consistent assumption.

Self-consistent interaction
Theoretically, the Voigt and Reuss interactions will result in upper and lower bounds of the stiffness, respectively. Therefore, the model was extended to also include the self-consistent interaction assumption, which is an intermediate approach between the upper and lower bounds. We seek the effective stiffness (C SC C ) such that: under the assumption that each of the unit cells are embedded in an equivalent homogeneous medium of stiffness C SC C . In other words, using this assumption means that the state of each UC in the aggregate is equivalent to the state of the UC embedded in a matrix with properties equivalent to the whole aggregate. Obviously, the properties of this equivalent homogeneous medium is not known a priori and are meant to be obtained.
To obtain the stiffness of a SFRC assuming the self-consistent interaction, we start with a matrix-inclusion problem [26], as schematically shown in Fig. 13. Two different domains are distinguished denoted by Ω m (matrix domain) and Ω r (reinforcement domain). In the present case, the reinforcement domain Ω r corresponds to the homogenized UC ( Fig. 1  (c)), and the matrix domain Ω m corresponds to the homogenized composite ( Fig. 1(d)). The stresses in the matrix and reinforcement domains are given by σ m = C m : ε, (21) σ r = C r : ε r , (22) where, C m and C r are assumed constant. Furthermore, for elliptic inclusions, ε r (and thus σ r ) will be constant over the elliptic reinforcement [26]. The strain in the reinforcement is related to the average strain as follows: Fig. 6. Four unidirectional sample RVEs for a polyamide/glass composite.  Fig. 7. Transverse cross section of three FE mesh of a UC for a polypropylene/flax SFRC.

Table 2
Homogenized elastic properties of a polypropylene/flax UC with different spatial discretizations (see Fig. 7).   Table 4 Characteristics of short fiber reinforced composites analyzed in Section 4.
where, A is the fourth order strain concentration tensor. Assuming perfectly bonded single ellipsoidal inclusion in an infinite elastic medium, the Eshelby solution [26] is used. The strain concentration tensor is the obtained as (see the mathematical manipulations in Appendix A): where I is the fourth order identity tensor and E is the fourth order Eshelby tensor. Assuming the self-consistent interaction implies that the matrix domain has the properties of the equivalent homogeneous medium: Using the homogenized stiffness (C) in the strain concentration tensor (Equation (24)) results in: and by using Equation (23), we obtain: : ε. (27) In this study, the inclusions are the Unit Cells. Thus, Equation (27) can be re-written as where, the superscript "SC" stands for the Self-Consistent assumption. Re-writing Equation (26) and using Equation (28) results in: The composite stress (in global configuration) is given by Using Equations (29) and (30), the composite stiffness, using the selfconsistent assumption, is identified as In Equation (31), the UC stiffness is given by   9. Comparison between model predictions and experimental results taken from Ref. [1]. Note that the all the results in Fig. 9 are normalized results. For all three sets of results (experimental, Voigt and Reuss) the obtained stiffness at different angles are compared to their corresponding stiffness at angle θ = 0 ∘ . It basically shows that the ratio of Voigt results at different cutting angles to the Voigt reference angle (θ = 0 ∘ ) has a better agreement to the corresponding experimental ratios. The comparison is presented in this way because absolute stiffness values were not given in the experiments.
It should be noted that Equation (31) does not have an explicit solution, and it has to be solved through an iterative procedure. We chose a fixed point iterative procedure which is explained in Appendix B. Also, since the UC used in this study, is not isotropic, the standard Eshelby tensor for isotropic inclusion can not be used. Therefore, for completeness, it is in Appendix C described how to calculate the Eshelby tensor for an anisotropic medium.

Implementation
Assuming self-consistent interaction between UCs, the composite stiffness is obtained by (see Equation (31)) In Equation (33), each UC stiffness C U,i is obtained as

Results including self-consistent interaction
An aspect ratio of 1 is considered for calculation of the Eshelby tensor for the simulations of this study (in the fourth order tensor D, calculated at the global configuration). It should be, however, noted that the actual aspect ratio was considered in the UC Finite Element simulations.

Polyamide/glass composites
The polyamide/glass composites, modeled in 4.3, is analyzed with the self-consistent model, as well. The obtained properties are shown in Fig. 14. It is seen the self-consistent model predictions are closer to the RVE results for this composite. The Orientation Averaging simulations take 49 s for 1000 UCs, and 107 s for 10,000 UCs, for all three interactions (for volume fraction of V f = 7%, and using a personal laptop). These are simulation times once the UC stiffness is known.
To make comparisons to RVE computational homogenization simulations for SFRCs with higher fiber volume fractions, glass fibers with a lower aspect ratio (λ = 5) was considered as well. The same fiber diameter (d f = 10 μm) is considered, but with a length of l f = 50 μm. With this aspect ratio, a fiber volume fraction of V f = 20% was achieved in the RVE analyses. Fig. 15 shows an RVE of this composite with its spatial discretization.

Polypropylene/flax composites
The last set of SFRCs modeled in this study are bio-composites made from Polypropylene matrix and short flax fibers. The details of these composite are given in Table 7.
Elastic properties of these composites are obtained using the OA method and compared against experimental results taken from Ref. [6]. It should be mentioned that we assumed isotropic fibers in this study. Due to the high aspect ratio of fibers and random distributions, the desired fiber volume fractions were not achieved for comparative RVE analyses. Fig. 16 shows a comparison between the model predictions and experiments. It is seen that with the Voigt assumption, good predictions of the elastic modulus is obtained for fiber volume fraction V f = 13%,21%. For the highest volume fraction (V f = 29%), the model (with the Voigt assumption) slightly overestimates the experiments. This might be a consequence of assuming a single length for all fiber reinforcements. It is shown by Andersons et al. [27] there is a fiber length distribution in extruded flax composites which is not accounted for in the simulations.    Table 6 gives obtained elastic modulus using the OA method and RVE analysis.

Conclusions
In this paper, a micro-mechanics based approach, combining Finite Element Analysis and Orientation Averaging, was developed to predict elastic properties of Short Fiber Reinforced Composites. A self-consistent interaction was developed as an intermediate approach between the upper and lower bounds. The authors believe the method is favorable with the following motivation: • Comparisons between the method results, experiments and computational homogenization of realistic RVEs show the capability of the method for adequate predictions; • The method is applicable to almost any short fiber composite with an arbitrary fiber volume fraction, fiber aspect ratio and fiber orientation distribution (among other properties); • The method is computationally efficient; • It is possible to extend the method and include other micro-structural phenomena such as inelasticity and matrix-fiber debonding.
As a result of the aforementioned advantages, this modeling approach can be used for real-life engineering components. According to the results obtained in this study, more comprehensive numerical examples and comparisons to experiments are needed to investigate the appropriateness of each interaction model for different composites. It should, however, be emphasized that theoretically, each interaction which results in a closer agreement with computational homogenization simulations on realistic RVEs, is a more accurate interaction.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix B
In this Appendix, the implicit solution to obtain the composite stiffness using the self-consistent interaction (Equation (31)) is explained. In an actual case of finite number of fibers, Equation (31) gets a summation form: In an iterative scheme, the composite stiffness at iteration (m +1) is given by For the initial guess, the composite stiffness obtained using Voigt interaction is used. To make sure that the model predictions using the self-consistent interaction is independent of the initial guess, the Reuss stiffness was also tried as the initial guess, and the same results were obtained. However, we stress that the Voigt stiffness is obtained at a (slightly) lower computational cost. The iterative procedure continues until the following convergence criterion is satisfied for all the components of the stiffness tensor: where Tol is a pre-defined value which is considered to be 0.001 in this study.

Appendix C
To obtain the composite stiffness using relation (31), it is needed to have the Eshelby tensor (E ijkl ). Since the UCs are not isotropic, the standard Eshelby tensor for isotropic inclusions can not be used. Instead, the Eshelby tensor for an anisotropic medium with stiffness tensor C ijkl , is given by Ref. [28]: where, the components of the fourth order tensor D are given by