# A Numerical Model to Predict the Anisotropy of Polymer Composites Reinforced with High-Aspect-Ratio Short Aramid Fibers

Authors: Jianhong Gao; Xiaoxiang Yang; Lihong Huang
Journal: Advances in Polymer Technology
 (2019)

---

## Abstract

Some fiber types have a high aspect ratio and it is very difficult to predict their composites using traditional finite element (FE) modeling. In this study, an FE model was developed to predict the anisotropy of composites reinforced by short aramid fibers. Three fiber distribution types were studied as follows: perfectly aligned, normally distributed, and randomly distributed fibers. The elastic constants were obtained, and, for different alignment angles and parameters in the fiber orientation distribution function, their numerical results were compared to those of the Mori–Tanaka model. Good agreement was obtained; thus, the employed FE model is an excellent and simple method to predict the isotropy and anisotropy of a composite with high-aspect-ratio fibers. Therefore, the FE model was employed to predict the orientation distribution of a composite fiber with a nonlinear matrix. The predicted and experimental results agree well.


---

## Body

## 1. Introduction

Short-fiber-reinforced polymer composites have been widely applied in industrial fields. In addition to the low cost and easy processing, versatility is an important characteristic that is determined by the performance of the composite’s constituents, fiber volume fraction, short-fiber aspect ratio, fiber orientation, and interface between the matrix and fiber [1]. A composite with aligned fibers can be recognized as an anisotropic material; however, a composite with random fibers is isotropic and has completely different mechanical properties. Generally, the fiber orientation varies because of different flow patterns inside the mold, processing conditions, and the rheological properties of the material [2]. Therefore, it is difficult for the fibers in a composite to be completely or randomly aligned. The mechanical properties of a real composite with short fibers are always intermediate between those of completely and randomly aligned fibers.Many analyses and numerical methods have been used to predict the elastic properties of short-fiber-reinforced composites. Among the analysis methods, the Mori-Tanaka model has been proven to produce the most accurate predictions [3]. To assess the effect of fiber orientation on the elastic properties, classical models, such as the Eshelby equivalent inclusion, Halpin–Tsai equation, self-consistent, shear lag, and bounding models and their extensions, have been applied to various materials. Mortazavian et al. [4] used experimental methods along with an analytical approach to study the tensile strength and elastic modulus of short-fiber-reinforced polymer composites. They proved that a laminate analogy and modified Tsai–Hill criteria provide satisfactory predictions. Tian et al. [5] investigated the effect of fiber orientation on the effective elastic properties of short-fiber-reinforced metal matrix composites using two-step mean-field homogenization procedures including the D-I/Reuss, M-T/Reuss, D-I/Reuss, M-T/Voigt, and D-M/V-R models. Jiang et al. [1] discussed the effects of fiber orientation distribution on the overall elastic properties of composites using Mori–Tanaka’s method. Ogierman et al. [6] studied the influence of process-induced orientation of the reinforcement on mechanical properties using Mori–Tanaka’s micromechanical model and a numerical method. Laspalas et al. [2] investigated the mechanical behavior of short-fiber-reinforced plastic composites by applying the Tandon–Weng micromechanical model considering the local fiber orientation distribution. Chen et al. [7] determined the effective moduli of composites containing misoriented fibers based on the Eshelby–Mori–Tanaka theory. Huang [8] presented an analytical approach to evaluate the orientation effects on the elastic properties of a composite containing randomly oriented fibers. Fu et al. [9] modeled fiber orientation distributions using a laminate analogy approach.It is convenient to use analytical approaches; however, compared to a numerical method, the microstructure field of deformation in composites cannot be directly obtained. Therefore, numerical methods which predict the anisotropy have been developed for different materials. Tian et al. [5] employed finite element (FE) methods with solid elements for short fibers to predict the effects of fiber orientation on the effective elastic properties of metal matrix composites. Ogierman et al. [6] investigated the effect of fiber orientation using a macroscale FE model with different boundary conditions. Laspalas et al. [2] studied the elasticity and failure of short-fiber-reinforced composites using FE simulation including an orientation average. Jansson et al. [10] employed the FE method to predict the anisotropic and nonlinear behavior of plastics reinforced by glass fibers. However, in FE modeling with multiple inclusions, “jamming” of fibers [11] easily occurs, in which the generated new fibers intersect the existing fibers. Extensive efforts have been devoted to developing methods for generating the representative volume elements (RVEs) of composites reinforced by inclusions of different shapes. Yi et al. [12] proposed modified random sequential adsorption (RSA) algorithms to generate an RVE of a random chopped fiber-reinforced composite material with straight and curved fibers. Harper et al. [11] presented a method to generate RVEs with random discontinuous carbon fibers in the plane and no limitation to the volume fraction of the fiber. Li et al. [13, 14] presented modified RSA algorithms for generating RVEs with complex microstructures such as fibers and spatially randomly distributed pores based on the microstructure information of carbon/carbon (C/C) composites.However, most numerical models have only focused on composites with a low-aspect-ratio fiber. The aspect ratio refers to the ratio of the fiber length to the diameter. For a high aspect ratio, particularly when the fiber radius is at a micron scale and the fiber length is at a millimeter scale, such as an aramid fiber, it is very difficult to establish a numerical FE model using a traditional method particularly when the volume fraction of the fibers is high. Only a few papers have reported on this topic [15, 16]. In our previous study [17], a numerical model was proposed to investigate the mechanical properties of an aramid-fiber-reinforced rubber composite with a high aspect ratio. It was proved that the FE model could simulate large deformation mechanical behavior well. However, the applicability of this proposed FE model in describing the anisotropy of a composite remains unknown.In this study, the FE model was further developed to predict the anisotropy of short-fiber-reinforced composites. This paper is structured as follows: In Section2, the analytical methods for predicting the elastic modulus of the composites with perfectly aligned, completely random, and probability orientation distribution fibers are introduced. In Section 3, the FE model and the algorithms of the different fiber orientation distributions are presented. In Section 4, the boundary conditions for the different mechanical behaviors are introduced and a means of obtaining the effective elastic constants is presented. In Section 5, the numerical results of the elastic constants with a different orientation distribution are compared to those of the Mori–Tanaka model and a high coincidence is shown. Finally, based on the aforementioned study, the FE model was applied to composites with a nonlinear matrix to predict the fiber orientation distribution. The numerical results were compared to the experimental results and found to agree well. This study verifies that the employed FE model is excellent to predict the elastic properties of a short-fiber-reinforced composite.


## 2. Analytical Models

### 2.1. Mori–Tanaka Model

In the Eshelby theory of equivalent inclusion, by replacing the elastic moduli of the inclusion with those of the matrix, an inhomogeneity problem can be equivalent to a homogeneity problem. Meanwhile, an eigenstrain is introduced into the inclusion domain. Based on the Eshelby equivalent inclusion theory [18], the Mori–Tanaka model [19] considers the mutual effect between the inclusions in the composite, thereby improving the accuracy of the prediction results. It has been proven that the Mori–Tanaka model is a reliable and simple method for inclusion problems and has been widely used to predict effective elastic moduli.It is supposed that the fiber axes coincide with the principal direction of the composite matrix. The composite is subjected to a uniform stressσ0 on the far field boundary, which is regarded as the average composite stress. For a homogeneous material, the constitutive relation is σ0=Lε0, where L is the elastic constant tensor of the composite. The average strain of the composite is as follows [19]:(1)ε0=1-C1ε0+C1ε1=I+C1AL0-1σ0where ε(0) and ε(1) are the average strain of the matrix and inclusions, respectively. L0 is the elastic constant tensor of the matrix.C1 is the volume fraction of the inclusions.A is expressed as follows:(2)A=L0+L1-L0C1I+1-C1S-1L0-L1where I is the unit tensor of four orders andS is the Eshelby tensor of four orders which is related to the elastic property of the matrix and the shape of the inclusions. The analytical expressions ofS for an ellipsoid inclusion problem can be found in many previous studies [20, 21]. Thus, the effective modulus of the composite can be obtained as follows:(3)L=L0I+C1A-1Hence, the effective elastic modulus L is an explicit function that is easy to obtain.


### 2.2. Analytical Model with Orientation Effects

In the previous section, the deduction was based on the assumption that the fibers are oriented along the principal matrix axes. However, in practice, because of processing conditions, it is difficult to achieve a perfect alignment and there is always a fiber orientation distribution in a composite. Therefore, to determine the real effective modulus, the effect of the fiber orientation distribution should be considered. In this section, we discuss the effective moduli for three cases: a perfectly aligned fiber whose principal axes present an angle referring to the principal matrix axes, a probability fiber orientation distribution, and completely random fibers.Generally, the orientation vector of fiberp can be expressed as p=cosθ,sinθcosφ,sinθsinφ, where θ and φ are the Euler angles as shown in Figure 1.x1,x2, andx3 are the three directions of the global coordinates, while x1′, x2′, and x3′ are the three directions of the local coordinates. In this case, the general framework in the previous section remains valid. The only change was to determine the tensor of Eq. (3) in the local coordinates referring to the global structural geometry. We denoted ui and uj′ as the unit vectors of the global coordinates (0,x1,x2, andx3) and local coordinates (0, x1′, x2′, and x3′), respectively. A relationship between them can be constructed as follows [1]:(4)ui=αjiuj′where the transformation matrix αji is as follows:(5)αji=cosφsinθcosφsinθsinφ-sinθcosθcosφcosθsinφ0-sinφcosφTherefore, for the case of perfectly aligned fibers, we obtained a relationship between the elastic tensor Lijkl of the global coordinates and Li′j′k′l′ of the local coordinates as follows:(6)Lijkl=αi′iαj′jαk′kαl′lLi′j′k′l′where αi′i,αj′j,αk′k,andαl′l are the elements of the transformation matrix αji. Lijkl and Li′j′k′l′ are the elements of the corresponding elastic tensor.Figure 1
Orientation vectorp of a spatially oriented short fiber.In the case of a misaligned fiber, the orientation average of the elastic tensor of the composite is provided by the following integration [1]:(7)L=∫θ=0π∫φ=02πLθ,φρθ,φsinθdθdφwhere ρ(θ,φ) is the probability density function of the fiber orientation. For a planar problem, Eq. (7) can be simplified combined with a discrete algorithm [5, 22] of the orientation average as follows:(8)L=∫θ=0πLθρθsin⁡θdθ=∑i=0NθLθi+Lθi+12ωθiwhere ∆θ=π/Nθ is the incremental angle, Nθ is the number of the angle incrementation, θi=i∆θ, andi is the integer within (0, Nθ). ω(θi) is the probability of the fiber orientation within the range of (θi,θi+1). L(θi) is the elastic constant when θ=θi. Eq. (8) is a simplified algorithm that avoids the difficulty of an integral, particularly for the complex probability density function ρ(θ). The precision of the result of Eq. (8) depends on Nθ because an increase in Nθ will lead to high calculation precision. Previous studies show that if Nθ is appropriately selected, high-precision calculation results for elastic constants can be achieved.In the case of randomly distributed fibers,ω(θi) in each angle range (θi,θi+1) is recognized to be equal. It can be calculated by ωθi=1/Nθ. Substituting it into Eq. (8), the elastic tensor is obtained.


## 2.1. Mori–Tanaka Model

In the Eshelby theory of equivalent inclusion, by replacing the elastic moduli of the inclusion with those of the matrix, an inhomogeneity problem can be equivalent to a homogeneity problem. Meanwhile, an eigenstrain is introduced into the inclusion domain. Based on the Eshelby equivalent inclusion theory [18], the Mori–Tanaka model [19] considers the mutual effect between the inclusions in the composite, thereby improving the accuracy of the prediction results. It has been proven that the Mori–Tanaka model is a reliable and simple method for inclusion problems and has been widely used to predict effective elastic moduli.It is supposed that the fiber axes coincide with the principal direction of the composite matrix. The composite is subjected to a uniform stressσ0 on the far field boundary, which is regarded as the average composite stress. For a homogeneous material, the constitutive relation is σ0=Lε0, where L is the elastic constant tensor of the composite. The average strain of the composite is as follows [19]:(1)ε0=1-C1ε0+C1ε1=I+C1AL0-1σ0where ε(0) and ε(1) are the average strain of the matrix and inclusions, respectively. L0 is the elastic constant tensor of the matrix.C1 is the volume fraction of the inclusions.A is expressed as follows:(2)A=L0+L1-L0C1I+1-C1S-1L0-L1where I is the unit tensor of four orders andS is the Eshelby tensor of four orders which is related to the elastic property of the matrix and the shape of the inclusions. The analytical expressions ofS for an ellipsoid inclusion problem can be found in many previous studies [20, 21]. Thus, the effective modulus of the composite can be obtained as follows:(3)L=L0I+C1A-1Hence, the effective elastic modulus L is an explicit function that is easy to obtain.


## 2.2. Analytical Model with Orientation Effects

In the previous section, the deduction was based on the assumption that the fibers are oriented along the principal matrix axes. However, in practice, because of processing conditions, it is difficult to achieve a perfect alignment and there is always a fiber orientation distribution in a composite. Therefore, to determine the real effective modulus, the effect of the fiber orientation distribution should be considered. In this section, we discuss the effective moduli for three cases: a perfectly aligned fiber whose principal axes present an angle referring to the principal matrix axes, a probability fiber orientation distribution, and completely random fibers.Generally, the orientation vector of fiberp can be expressed as p=cosθ,sinθcosφ,sinθsinφ, where θ and φ are the Euler angles as shown in Figure 1.x1,x2, andx3 are the three directions of the global coordinates, while x1′, x2′, and x3′ are the three directions of the local coordinates. In this case, the general framework in the previous section remains valid. The only change was to determine the tensor of Eq. (3) in the local coordinates referring to the global structural geometry. We denoted ui and uj′ as the unit vectors of the global coordinates (0,x1,x2, andx3) and local coordinates (0, x1′, x2′, and x3′), respectively. A relationship between them can be constructed as follows [1]:(4)ui=αjiuj′where the transformation matrix αji is as follows:(5)αji=cosφsinθcosφsinθsinφ-sinθcosθcosφcosθsinφ0-sinφcosφTherefore, for the case of perfectly aligned fibers, we obtained a relationship between the elastic tensor Lijkl of the global coordinates and Li′j′k′l′ of the local coordinates as follows:(6)Lijkl=αi′iαj′jαk′kαl′lLi′j′k′l′where αi′i,αj′j,αk′k,andαl′l are the elements of the transformation matrix αji. Lijkl and Li′j′k′l′ are the elements of the corresponding elastic tensor.Figure 1
Orientation vectorp of a spatially oriented short fiber.In the case of a misaligned fiber, the orientation average of the elastic tensor of the composite is provided by the following integration [1]:(7)L=∫θ=0π∫φ=02πLθ,φρθ,φsinθdθdφwhere ρ(θ,φ) is the probability density function of the fiber orientation. For a planar problem, Eq. (7) can be simplified combined with a discrete algorithm [5, 22] of the orientation average as follows:(8)L=∫θ=0πLθρθsin⁡θdθ=∑i=0NθLθi+Lθi+12ωθiwhere ∆θ=π/Nθ is the incremental angle, Nθ is the number of the angle incrementation, θi=i∆θ, andi is the integer within (0, Nθ). ω(θi) is the probability of the fiber orientation within the range of (θi,θi+1). L(θi) is the elastic constant when θ=θi. Eq. (8) is a simplified algorithm that avoids the difficulty of an integral, particularly for the complex probability density function ρ(θ). The precision of the result of Eq. (8) depends on Nθ because an increase in Nθ will lead to high calculation precision. Previous studies show that if Nθ is appropriately selected, high-precision calculation results for elastic constants can be achieved.In the case of randomly distributed fibers,ω(θi) in each angle range (θi,θi+1) is recognized to be equal. It can be calculated by ωθi=1/Nθ. Substituting it into Eq. (8), the elastic tensor is obtained.


## 3. Numerical Model Development

### 3.1. Method

In the traditional method of finite element analysis (FEA), the solid element is always used to simulate the mechanical behavior of the fiber. However, using this method, it is very difficult to establish a composite model with a high fiber volume fraction [12]. Therefore, in [17], an FE model was introduced for a planar problem of the elastic properties of a short-fiber-reinforced rubber composite. However, thus far, we did not know the applicability of this model for the analysis of the anisotropic properties of a composite.The FE model was established as a square whose size isD×D [17], as shown in Figure 2. Short fibers were generated and the linear hybrid beam element (BH21H) was employed in the ABAQUS FE software. For the matrix, plane stress continuum elements (CPS8R) were used. Using the embedded technique, which has been verified to yield quite similar results with those for the traditional method [23], the mesh of the matrix and the fibers are independent of each other, as shown in Figure 3, and are more regular, leading to more precise numerical results. Another advantage is that a model with high volume fraction and high-aspect-ratio fibers is realized.Figure 2
FE model of a short-fiber-reinforced composite, square shape with sizeD×D.Figure 3
FE models for three fiber distribution types:D=20 mm, lf=1 mm.
(a)
Perfectly aligned fibers (θ=50°) (b)
Normally distributed fibers (c)
Completely random distribution


### 3.2. Algorithm for Fiber Generation

In this study, planar problems with three fiber distribution types, i.e., perfectly aligned, probability function distributed, and completely random fibers, were discussed. The detailed flow chart can be referred to in the literature [17]. To avoid fiber intersection, the RSA algorithm was employed using Python language in the ABAQUS software. Three random numbers,x,y, andθ, which represent the moving distance along thex direction, moving distance along they direction, and the rotation angle of the fiber, respectively, were used. To ensure higher accuracy of the fiber volume fraction, each fiber was generated with an equal length in the model, including the fibers around the model edge. For this reason, the ranges ofx andy were set as follows: x∈(-D/2+lf/2(1-sin2θ)-1/2,D/2-lf/2(1-sin2θ)-1/2), y∈(-D/2+lf/2(1-cos2θ)-1/2,D/2-lf/2(1-cos2θ)-1/2), where lf is the fiber length. For different fiber orientation distributions, the FE models are shown in Figure 3 and the ranges ofθare interpreted as follows:(1)
Perfectly aligned fibers:θ is the orientation angle of the fibers which is invariable for all the fibers in the model.(2)
Completely random fibers:θ isa random number within (-π/2, π/2).(3)
Misaligned fibers with a probability function distribution: the incremental angle∆θ=π/Nθ is set where Nθ is the number of the angle incrementation. First, nfiiscalculated according to a given probability function of the fiber orientation ρ(θi). nfi and ρ(θi) are the fiber amount and orientation probability, respectively, distributed within an angle range (∆θ×i, ∆θ×(i+1)). Subsequently, fibers are generated with random angles within the range (∆θ×i, ∆θ×(i+1)) and the corresponding generated fiber amount is nfi.i is an integer between 0 and Nθ. Therefore, this algorithm is a discretized method whose precision depends on ∆θ. Decreasing the ∆θ enhances the calculation precision.


### 3.3. Probability Density Function of the Fiber Distribution

For most real composites reinforced with aramid fibers of a fixed aspect ratio, there are always different fiber orientation distributions that can be described by a symmetric function, such as a normal function, and a nonsymmetric function, e.g., a Gamma, Weibull, or exponential function. To study the effect of fiber orientation distribution on the elastic properties of the composite, a normal distribution function was considered in this study. For the planar problem, only one variable, the Euler angleθ, is included such that the function can be expressed as follows [1]:(9)ρθ=1σ2πexp⁡-12θ-μσ2where μ (-∞<μ<∞) is the mean and σ2 (σ > 0) is the variance of the distribution. Figure 4 shows the fiber distribution when in Eq. (9) θ=0; μ=0; and σ2=30, 100, and 500 (degree measure was used for σ2). The diagram shows that as σ2 increases, the fiber distribution is increasingly disordered. There are two extremes: when σ2=0, the fibers perfectly align, and when σ2 is sufficiently large, the fiber distribution is completely random. Indeed, when σ2 increases to a large value, the range of θ, (-π/2, π/2) is unable to accommodate all the fibers because of the normal distribution characteristics. However, fortunately, the effect of the orientation angle of the fiber on the mechanical properties of the composite is periodic, which indicates that all the fibers of a θi±nπ orientation angle exert the same effect. According to this, ρ(θi) can be flexibly calculated for large σ2.Figure 4
Effect of varianceσ2 of the normal probability density function on fiber distribution: θ=0, μ=0.


## 3.1. Method

In the traditional method of finite element analysis (FEA), the solid element is always used to simulate the mechanical behavior of the fiber. However, using this method, it is very difficult to establish a composite model with a high fiber volume fraction [12]. Therefore, in [17], an FE model was introduced for a planar problem of the elastic properties of a short-fiber-reinforced rubber composite. However, thus far, we did not know the applicability of this model for the analysis of the anisotropic properties of a composite.The FE model was established as a square whose size isD×D [17], as shown in Figure 2. Short fibers were generated and the linear hybrid beam element (BH21H) was employed in the ABAQUS FE software. For the matrix, plane stress continuum elements (CPS8R) were used. Using the embedded technique, which has been verified to yield quite similar results with those for the traditional method [23], the mesh of the matrix and the fibers are independent of each other, as shown in Figure 3, and are more regular, leading to more precise numerical results. Another advantage is that a model with high volume fraction and high-aspect-ratio fibers is realized.Figure 2
FE model of a short-fiber-reinforced composite, square shape with sizeD×D.Figure 3
FE models for three fiber distribution types:D=20 mm, lf=1 mm.
(a)
Perfectly aligned fibers (θ=50°) (b)
Normally distributed fibers (c)
Completely random distribution


## 3.2. Algorithm for Fiber Generation

In this study, planar problems with three fiber distribution types, i.e., perfectly aligned, probability function distributed, and completely random fibers, were discussed. The detailed flow chart can be referred to in the literature [17]. To avoid fiber intersection, the RSA algorithm was employed using Python language in the ABAQUS software. Three random numbers,x,y, andθ, which represent the moving distance along thex direction, moving distance along they direction, and the rotation angle of the fiber, respectively, were used. To ensure higher accuracy of the fiber volume fraction, each fiber was generated with an equal length in the model, including the fibers around the model edge. For this reason, the ranges ofx andy were set as follows: x∈(-D/2+lf/2(1-sin2θ)-1/2,D/2-lf/2(1-sin2θ)-1/2), y∈(-D/2+lf/2(1-cos2θ)-1/2,D/2-lf/2(1-cos2θ)-1/2), where lf is the fiber length. For different fiber orientation distributions, the FE models are shown in Figure 3 and the ranges ofθare interpreted as follows:(1)
Perfectly aligned fibers:θ is the orientation angle of the fibers which is invariable for all the fibers in the model.(2)
Completely random fibers:θ isa random number within (-π/2, π/2).(3)
Misaligned fibers with a probability function distribution: the incremental angle∆θ=π/Nθ is set where Nθ is the number of the angle incrementation. First, nfiiscalculated according to a given probability function of the fiber orientation ρ(θi). nfi and ρ(θi) are the fiber amount and orientation probability, respectively, distributed within an angle range (∆θ×i, ∆θ×(i+1)). Subsequently, fibers are generated with random angles within the range (∆θ×i, ∆θ×(i+1)) and the corresponding generated fiber amount is nfi.i is an integer between 0 and Nθ. Therefore, this algorithm is a discretized method whose precision depends on ∆θ. Decreasing the ∆θ enhances the calculation precision.


## 3.3. Probability Density Function of the Fiber Distribution

For most real composites reinforced with aramid fibers of a fixed aspect ratio, there are always different fiber orientation distributions that can be described by a symmetric function, such as a normal function, and a nonsymmetric function, e.g., a Gamma, Weibull, or exponential function. To study the effect of fiber orientation distribution on the elastic properties of the composite, a normal distribution function was considered in this study. For the planar problem, only one variable, the Euler angleθ, is included such that the function can be expressed as follows [1]:(9)ρθ=1σ2πexp⁡-12θ-μσ2where μ (-∞<μ<∞) is the mean and σ2 (σ > 0) is the variance of the distribution. Figure 4 shows the fiber distribution when in Eq. (9) θ=0; μ=0; and σ2=30, 100, and 500 (degree measure was used for σ2). The diagram shows that as σ2 increases, the fiber distribution is increasingly disordered. There are two extremes: when σ2=0, the fibers perfectly align, and when σ2 is sufficiently large, the fiber distribution is completely random. Indeed, when σ2 increases to a large value, the range of θ, (-π/2, π/2) is unable to accommodate all the fibers because of the normal distribution characteristics. However, fortunately, the effect of the orientation angle of the fiber on the mechanical properties of the composite is periodic, which indicates that all the fibers of a θi±nπ orientation angle exert the same effect. According to this, ρ(θi) can be flexibly calculated for large σ2.Figure 4
Effect of varianceσ2 of the normal probability density function on fiber distribution: θ=0, μ=0.


## 4. Numerical Implementation

### 4.1. Model Parameters

In this study, aramid fiber, which is widely applied in industrial products, was employed for the investigation. Typically, aramid fiber has a high aspect ratio because its diameter is very small at approximately 0.012 mm. For chopped fibers, the length is changed to meet product requirements. Because an aramid fiber is short and rigid, it is straight at the beginning. Therefore, a straight line can be employed to simulate its initial elastic property in the FE model. According to the product manual, a fiber of modulusEf= 70 GPa was selected for this study and the Poisson ratio was 0.3. For the polymer composite, a matrix modulus of 1 GPa was employed. The mechanical properties of the composites containing the aforementioned components were similar to those of a composite with a resin matrix, which is widely used in the industry. The volume fraction of a fiber was defined by vf=πrf2lfnf/D×D×t.tis the model thickness, set to be 0.1, and nf is the total fiber amount in the model. According to [17], with increasing D/lf, the numerical results will converge. D/lf=30 was used in this study.


### 4.2. Boundary Conditions

The FE model in this study can be regarded as at a mesoscopic scale. If the model with a minimal size has the same effective mechanical properties as the bulk material, it can be recognized as the RVE. Generally, for a microscopic model, a periodic boundary condition will be employed [24]. In the literature [11], Harper et al. adopted periodic boundary conditions and provided effective material properties within an inner RVE region. In [17], considering the Hill condition [25], four boundary conditions were imposed, i.e., uniform traction, uniform displacement, mixed boundary condition, and periodic boundary condition, for the model uniaxial tension. It has been proven that, with increasing model size, results with different boundary condition tend to be the same.Therefore, we believed the uniform displacement boundary condition would provide the most reasonable and correct numerical results for the elastic property of uniaxial tension. The boundary conditions for the different elastic behaviors are described in Table1. The variableu is the displacement load of the node on the corresponding edge; the superscript represents the load direction. The subscripti orj represents any node on the corresponding edge and thus the expression uiy=ujy implies that the displacement of each node is equal. For simple shear deformation, ui,x=L/2x=ui,x=-L/2x implies that the displacement along thex direction of the node on the x=L/2 edge is equal to that of the corresponding node on the x=-L/2 edge. Among all the variables listed in Table 1, only δx and δy are provided, which are fixed values for each elastic behavior. For different elastic behaviors, the deformations are shown in Figure 5. Their elastic constants can be calculated by extracting the results using Python as follows:Table 1
Boundary conditions for the different elastic constants.
Elastic behavior Elastic constant(s) Boundary conditions Uniaxial tension along thex direction E x, vx x=D/2, ux=δx; x=- D/2, ux=-δx; y= D/2, uiy=ujy; y=- D/2, uiy=ujy Uniaxial tension along they direction E y, vy y=D/2, uy=δy; y=- D/2, uy=-δy; x= D/2, uix=ujx; x=- D/2, uix=ujx Simple shear E x y y=D/2, ux=δx, uy=0; y=- D/2, ux=uy=0; u i , x = L / 2 x = u i , x = - L / 2 x; ui,x=L/2y=ui,x=-L/2yFigure 5
Different elastic deformations of short-fiber-reinforced composite.
(a)
Uniaxial tensile deformation along the x direction (b)
Uniaxial tensile deformation along the y direction (c)
Shear deformation(  1) Elastic Moduli E x and v x (10) σ x = R F D × t = ∑ i = 1 n f x = D / 2 i D × t = - ∑ i = 1 n f x = - D / 2 i D × t (11) ε x = U D = ∑ i = 1 n u x = D / 2 i - ∑ i = 1 n u x = - D / 2 i n D (12) ε y = U D = ∑ i = 1 n u y = D / 2 i - ∑ i = 1 n u y = - D / 2 i n Dwhere σx is the stress in thex direction of the model imposed by uniaxial tension. εx and εy are the strain in thex andy directions, respectively.RF is the reaction force of the model edge which is loaded on andn is the node amount of one edge, which is the same for each edge of the square model. fx=D/2i and ux=D/2i are the reaction force and displacement of the node on the x=D/2 edge, respectively. The elastic modulus and Poisson ratio in thex direction can be expressed, respectively, as follows: Ex=σx/εx and vx=εy/εx.(  2) Elastic Moduli E y and v y. Replacingx withy in Eqs. (10), (11), and (12), Ey and vy are obtained.(  3) Shear Modulus E x y (13) E x y = τ x y γ x y = ∑ i = 1 n f y = D / 2 i arctan ⁡ δ x / Dwhere τxy and γxy are the model shear stress and strain. Other variables have been previously defined.


## 4.1. Model Parameters

In this study, aramid fiber, which is widely applied in industrial products, was employed for the investigation. Typically, aramid fiber has a high aspect ratio because its diameter is very small at approximately 0.012 mm. For chopped fibers, the length is changed to meet product requirements. Because an aramid fiber is short and rigid, it is straight at the beginning. Therefore, a straight line can be employed to simulate its initial elastic property in the FE model. According to the product manual, a fiber of modulusEf= 70 GPa was selected for this study and the Poisson ratio was 0.3. For the polymer composite, a matrix modulus of 1 GPa was employed. The mechanical properties of the composites containing the aforementioned components were similar to those of a composite with a resin matrix, which is widely used in the industry. The volume fraction of a fiber was defined by vf=πrf2lfnf/D×D×t.tis the model thickness, set to be 0.1, and nf is the total fiber amount in the model. According to [17], with increasing D/lf, the numerical results will converge. D/lf=30 was used in this study.


## 4.2. Boundary Conditions

The FE model in this study can be regarded as at a mesoscopic scale. If the model with a minimal size has the same effective mechanical properties as the bulk material, it can be recognized as the RVE. Generally, for a microscopic model, a periodic boundary condition will be employed [24]. In the literature [11], Harper et al. adopted periodic boundary conditions and provided effective material properties within an inner RVE region. In [17], considering the Hill condition [25], four boundary conditions were imposed, i.e., uniform traction, uniform displacement, mixed boundary condition, and periodic boundary condition, for the model uniaxial tension. It has been proven that, with increasing model size, results with different boundary condition tend to be the same.Therefore, we believed the uniform displacement boundary condition would provide the most reasonable and correct numerical results for the elastic property of uniaxial tension. The boundary conditions for the different elastic behaviors are described in Table1. The variableu is the displacement load of the node on the corresponding edge; the superscript represents the load direction. The subscripti orj represents any node on the corresponding edge and thus the expression uiy=ujy implies that the displacement of each node is equal. For simple shear deformation, ui,x=L/2x=ui,x=-L/2x implies that the displacement along thex direction of the node on the x=L/2 edge is equal to that of the corresponding node on the x=-L/2 edge. Among all the variables listed in Table 1, only δx and δy are provided, which are fixed values for each elastic behavior. For different elastic behaviors, the deformations are shown in Figure 5. Their elastic constants can be calculated by extracting the results using Python as follows:Table 1
Boundary conditions for the different elastic constants.
Elastic behavior Elastic constant(s) Boundary conditions Uniaxial tension along thex direction E x, vx x=D/2, ux=δx; x=- D/2, ux=-δx; y= D/2, uiy=ujy; y=- D/2, uiy=ujy Uniaxial tension along they direction E y, vy y=D/2, uy=δy; y=- D/2, uy=-δy; x= D/2, uix=ujx; x=- D/2, uix=ujx Simple shear E x y y=D/2, ux=δx, uy=0; y=- D/2, ux=uy=0; u i , x = L / 2 x = u i , x = - L / 2 x; ui,x=L/2y=ui,x=-L/2yFigure 5
Different elastic deformations of short-fiber-reinforced composite.
(a)
Uniaxial tensile deformation along the x direction (b)
Uniaxial tensile deformation along the y direction (c)
Shear deformation(  1) Elastic Moduli E x and v x (10) σ x = R F D × t = ∑ i = 1 n f x = D / 2 i D × t = - ∑ i = 1 n f x = - D / 2 i D × t (11) ε x = U D = ∑ i = 1 n u x = D / 2 i - ∑ i = 1 n u x = - D / 2 i n D (12) ε y = U D = ∑ i = 1 n u y = D / 2 i - ∑ i = 1 n u y = - D / 2 i n Dwhere σx is the stress in thex direction of the model imposed by uniaxial tension. εx and εy are the strain in thex andy directions, respectively.RF is the reaction force of the model edge which is loaded on andn is the node amount of one edge, which is the same for each edge of the square model. fx=D/2i and ux=D/2i are the reaction force and displacement of the node on the x=D/2 edge, respectively. The elastic modulus and Poisson ratio in thex direction can be expressed, respectively, as follows: Ex=σx/εx and vx=εy/εx.(  2) Elastic Moduli E y and v y. Replacingx withy in Eqs. (10), (11), and (12), Ey and vy are obtained.(  3) Shear Modulus E x y (13) E x y = τ x y γ x y = ∑ i = 1 n f y = D / 2 i arctan ⁡ δ x / Dwhere τxy and γxy are the model shear stress and strain. Other variables have been previously defined.


## 5. Verification and Discussion

### 5.1. Perfectly Aligned Fibers

To investigate the applicability of the proposed model for predicting the anisotropic properties of composites with perfectly aligned fibers, models with a different fiber orientation angleθ relative to the x axis were established. The boundary conditions shown in Table 1 were imposed for the different elastic deformations and the results were compared to those of the Mori–Tanaka model. Figure 6 shows this comparison when lf=1 mm and vf=0.01. It can be seen that all the elastic constants in the numerical results, Ex, Ey, Exy, vx, and vy, agree well with the results of the Mori–Tanaka model.Figure 6
Comparison between the numerical and Mori–Tanaka model results for different elastic constants of a composite withlf=1 mm and vf=0.01.
(a)
E x, Ey, and Exy (b)
v x and vyFor the elastic modulusEx, assuming that the uniaxial tension along thex direction is imposed, the minimum occurs at approximately θ=60° but not at 90°. To understand this, the matrix and fiber stresses were investigated. By averaging the stress of the matrix and fiber of the RVE, the contributions from the resin and fiber were both considered as follows [11]:(14)σxave=σxfiber-ave+σxmatrix-ave=1V∫vσxfiberdV+∫vσxmatrixdV=1D2t∑a∑iSFbiAbSVOLbicos⁡θ+∑b∑ismciIVOLciwhere σxave, σxfiber-ave, and σxmatrix-ave are the average stress of the model, beam, and matrix, respectively, andV is the RVE volume. The subscriptsa andb are the number of beams and matrix elements, respectively, and the subscripti is their integration point number.SF is the section force of the beam which is multiplied by cos⁡θ; its component in the loading directionx can then be obtained. sm is the stress component of each integrated point of the matrix material in the loading direction. Ab is the cross-sectional area of the beam element. SVOL is the integrated section volume of each beam element.IVOLis the parameter of the integration point volume of the matrix element.For the composite withlf=1 mm and vf=0.01, σxave, σxfiber-ave, and σxmatrix-ave were obtained as shown in Figure 7. It can be seen that when θ<60°, as θ increases, the fiber stresses decrease, which weakens the fiber reinforcement effect. Until θ increases to approximately 60°, the stress in the fibers is compressive. Meanwhile, for the average matrix stress, it decreases when θ<50° and then increases from nearly 60°. Thus, the total stress minimum occurs at approximately θ=60°, which explains the minimum EX of θ=60°.Figure 7
Average stress of the total model, matrix, and fibers of a composite withlf=1 mm and vf=0.01 when the uniaxial tension of thex direction is imposed.


### 5.2. Different Fiber Orientation Distributions

To express the anisotropy of the composites, a parameterαij was introduced that refers to the deviation degree from isotropy. The isotropic factor can be defined as follows [26]:(15)αij=L-ijT-L-ij2L-9-i-j9-i-j,withi,j=1,2,and3,where L-ij is the corresponding element in the stiffness matrix and L-ijT is expressed by L-ijT=L-ii+L-jj/2. For the planar problem, only the factor α12 was used to describe the deviation degree from isotropy. The closer α12 is to 1.0, the more isotropic the composite is.Models with differentvariancesofthedistributionσ2 (σ > 0) were established when θ in Eq. (9) was 0°. σ2=0 implies that the fibers are perfectly aligned. If σ2 increases to a sufficiently large value, the fibers are considered to be randomly distributed. For one σ2, in the same manner, the elastic constants, Ex, Ey, Exy, vx, andvy, are calculated via simulation and an analytical model. Then, α12 is obtained. The relationship between α12 and the five elastic constants is plotted as shown in Figure 8. The comparisons between the numerical and Mori–Tanaka model results are shown to agree well. As shown in Figure 8, the minimum longitudinal modulus Ex occurs when α12 is near 1 indicating a randomly distributed fiber. Meanwhile, the transverse modulus, Ey, and shear modulus, Exy, both show a maximum value. In addition, with decreasing α12, Ey seldom changes at first but the effect of α12 on Ex is greater than that on Ey.Figure 8
Relationship between elastic constants and isotropic factorα12.
(a)
E x, Ey, and Exy (b)
v x and vyIn addition, the relationship between the isotropic factorα12 and the variancesσ2 is shown in Figure 9. It can be seen that when σ2 increases to approximately 2000, α12 is approximately 1.0 and tends to remain constant, which indicates a randomly distributed fiber. Figure 9 can serve as a reference to determine the degree of the anisotropy of composites and is applicable to all short-fiber-reinforced composites of normal orientation distribution. In addition, other anisotropic models with different fiber microstructural parameters were calculated and there was good agreement as shown in Figures 6 and 7. Therefore, we believe that the employed FE model predicts the anisotropy of a composite very well.Figure 9
Relationship between the isotropic factorα12 and variancesσ2.


## 5.1. Perfectly Aligned Fibers

To investigate the applicability of the proposed model for predicting the anisotropic properties of composites with perfectly aligned fibers, models with a different fiber orientation angleθ relative to the x axis were established. The boundary conditions shown in Table 1 were imposed for the different elastic deformations and the results were compared to those of the Mori–Tanaka model. Figure 6 shows this comparison when lf=1 mm and vf=0.01. It can be seen that all the elastic constants in the numerical results, Ex, Ey, Exy, vx, and vy, agree well with the results of the Mori–Tanaka model.Figure 6
Comparison between the numerical and Mori–Tanaka model results for different elastic constants of a composite withlf=1 mm and vf=0.01.
(a)
E x, Ey, and Exy (b)
v x and vyFor the elastic modulusEx, assuming that the uniaxial tension along thex direction is imposed, the minimum occurs at approximately θ=60° but not at 90°. To understand this, the matrix and fiber stresses were investigated. By averaging the stress of the matrix and fiber of the RVE, the contributions from the resin and fiber were both considered as follows [11]:(14)σxave=σxfiber-ave+σxmatrix-ave=1V∫vσxfiberdV+∫vσxmatrixdV=1D2t∑a∑iSFbiAbSVOLbicos⁡θ+∑b∑ismciIVOLciwhere σxave, σxfiber-ave, and σxmatrix-ave are the average stress of the model, beam, and matrix, respectively, andV is the RVE volume. The subscriptsa andb are the number of beams and matrix elements, respectively, and the subscripti is their integration point number.SF is the section force of the beam which is multiplied by cos⁡θ; its component in the loading directionx can then be obtained. sm is the stress component of each integrated point of the matrix material in the loading direction. Ab is the cross-sectional area of the beam element. SVOL is the integrated section volume of each beam element.IVOLis the parameter of the integration point volume of the matrix element.For the composite withlf=1 mm and vf=0.01, σxave, σxfiber-ave, and σxmatrix-ave were obtained as shown in Figure 7. It can be seen that when θ<60°, as θ increases, the fiber stresses decrease, which weakens the fiber reinforcement effect. Until θ increases to approximately 60°, the stress in the fibers is compressive. Meanwhile, for the average matrix stress, it decreases when θ<50° and then increases from nearly 60°. Thus, the total stress minimum occurs at approximately θ=60°, which explains the minimum EX of θ=60°.Figure 7
Average stress of the total model, matrix, and fibers of a composite withlf=1 mm and vf=0.01 when the uniaxial tension of thex direction is imposed.


## 5.2. Different Fiber Orientation Distributions

To express the anisotropy of the composites, a parameterαij was introduced that refers to the deviation degree from isotropy. The isotropic factor can be defined as follows [26]:(15)αij=L-ijT-L-ij2L-9-i-j9-i-j,withi,j=1,2,and3,where L-ij is the corresponding element in the stiffness matrix and L-ijT is expressed by L-ijT=L-ii+L-jj/2. For the planar problem, only the factor α12 was used to describe the deviation degree from isotropy. The closer α12 is to 1.0, the more isotropic the composite is.Models with differentvariancesofthedistributionσ2 (σ > 0) were established when θ in Eq. (9) was 0°. σ2=0 implies that the fibers are perfectly aligned. If σ2 increases to a sufficiently large value, the fibers are considered to be randomly distributed. For one σ2, in the same manner, the elastic constants, Ex, Ey, Exy, vx, andvy, are calculated via simulation and an analytical model. Then, α12 is obtained. The relationship between α12 and the five elastic constants is plotted as shown in Figure 8. The comparisons between the numerical and Mori–Tanaka model results are shown to agree well. As shown in Figure 8, the minimum longitudinal modulus Ex occurs when α12 is near 1 indicating a randomly distributed fiber. Meanwhile, the transverse modulus, Ey, and shear modulus, Exy, both show a maximum value. In addition, with decreasing α12, Ey seldom changes at first but the effect of α12 on Ex is greater than that on Ey.Figure 8
Relationship between elastic constants and isotropic factorα12.
(a)
E x, Ey, and Exy (b)
v x and vyIn addition, the relationship between the isotropic factorα12 and the variancesσ2 is shown in Figure 9. It can be seen that when σ2 increases to approximately 2000, α12 is approximately 1.0 and tends to remain constant, which indicates a randomly distributed fiber. Figure 9 can serve as a reference to determine the degree of the anisotropy of composites and is applicable to all short-fiber-reinforced composites of normal orientation distribution. In addition, other anisotropic models with different fiber microstructural parameters were calculated and there was good agreement as shown in Figures 6 and 7. Therefore, we believe that the employed FE model predicts the anisotropy of a composite very well.Figure 9
Relationship between the isotropic factorα12 and variancesσ2.


## 6. Application to Nonlinear Material

### 6.1. Effect of Orientation Distribution

The hyperelastic property of rubber composites reinforced by aramid fibers has been studied [17]. By averaging the experimental stress–strain curves along the x-direction (along the flowing direction) and y-direction (vertical to the flowing direction) specimens, the averaged curve can be recognized as the mechanical response of a composite with randomly distributed fibers. By comparing to the experimental data, the numerical model of isotropy with random fibers predicts the nonlinear elastic response of stress–strain well. For real specimens, because the short fibers in the composite are always oriented along the compound flowing direction, anisotropy always exists; however, only by means of macromechanics, it is difficult to determine the orientation degree. As described in this section, the anisotropic FE model was extended to nonlinear material.The material parameters of the two components, rubber and short aramid fibers, are the same as [17]. The Ogden N3 hyperelastic constitutive model was employed to describe the rubber matrix elastic behavior. Based on the tensile experiment results, the material parameters were obtained as follows: μ1=-8.116, α1=1.975, μ2=5.725, α2=2.388, μ3=5.467, and α3=-4.495. For the aramid fibers, the elastic modulus is 132 GPa, the fiber diameter is 0.012 mm, and the studied fiber length is 1.5 mm. The difference in this study was the composite fiber distribution, which was established with an anisotropic morphology. Still, a perfect interface and straight morphology for the aramid fibers at the initial deformation were supposed. For fiber length lf=1.5 mm and fiber volume fraction vf=0.0049, using the anisotropic FE model, the stress–strain response of the composites with different perfectly aligned angles was obtained as shown in Figure 10. The curve shapes with blue scatters shown in Figures 10(a) and 10(b) are similar to those of Ex and Ey of the linear material shown in Figure 6, respectively. The minimum longitudinal modulus shown in Figure 11 also occurs when θ=60°. The stress–strain responses of the composites with different variancesσ2 from the numerical results are shown in Figure 11, which corresponds to Figure 8 for the linear material.Figure 10
Response of stress–strain relationship of composites with an aligned angleθ of fibers, lf=1.5 mm, and vf=0.0049.
(a)
Stress–strain response, when the composites were loaded along the x direction (b)
Stress–strain response, when the composites were loaded along the y directionFigure 11
Response of the stress–strain of composites with differentvariancesσ2 from the numerical results: lf=1.5 mm, vf=0.0049.
(a)
Stress–strain response when the composites were loaded along the x direction (b)
Stress-strain response when the composites were loaded along the y direction


### 6.2. Determination of the Orientation Distribution

For the hyperelastic material, the constitutive behavior is always expressed in the form of strain energyW. The hyperelastic constitutive model of an isotropic material has been comprehensively studied [27, 28]. However, for an anisotropic material, many terms remain to be explored [29–31]. To describe the orientation degree of the fibers as the isotropic factor α12 of the linear material, a variable, β=Wx/Wy was introduced, where Wx and Wy are the strain energy of the elastic response when the tensile load is the along x and y directions, respectively. According to the numerical results, for fiber length lf=1.5 mm and fiber volume fraction vf=0.0049, the relationships between β and Wx and Wy are shown in Figure 12. It can be seen that when β is near 1.0, Ex and Ey coincide.Figure 12
Relationship betweenWx (Wy) and the isotropic factor β for a rubber composite reinforced with aramid fibers: lf=1.5 mm, vf=0.0049.From the effect ofσ2 on β as shown in Figure 13, according to the experimental result β=1.09, the corresponding σ2 of the fiber orientation distribution was determined. σ2, which was approximately 1400, is the most reasonable value. Similarly, σ2 for lf=1.5 mm and vf=0.0146 was determined to be approximately 3000. The experimental results [12] and the corresponding stress–strain curves of the determined σ2 are compared in Figure 14. There are good agreements, indicating that the proposed method predicting the fiber orientation distribution is applicable.Figure 13
Relationship between the isotropic factorβ and variancesσ2 for a rubber composite reinforced with aramid fibers: lf=1.5 mm, vf=0.0049.Figure 14
Comparison between the experimental and numerical results of the model with the determinedσ2.
(a)
l f=1.5 mm, vf=0.0049 (b)
l f=1.5 mm, vf=0.0146


## 6.1. Effect of Orientation Distribution

The hyperelastic property of rubber composites reinforced by aramid fibers has been studied [17]. By averaging the experimental stress–strain curves along the x-direction (along the flowing direction) and y-direction (vertical to the flowing direction) specimens, the averaged curve can be recognized as the mechanical response of a composite with randomly distributed fibers. By comparing to the experimental data, the numerical model of isotropy with random fibers predicts the nonlinear elastic response of stress–strain well. For real specimens, because the short fibers in the composite are always oriented along the compound flowing direction, anisotropy always exists; however, only by means of macromechanics, it is difficult to determine the orientation degree. As described in this section, the anisotropic FE model was extended to nonlinear material.The material parameters of the two components, rubber and short aramid fibers, are the same as [17]. The Ogden N3 hyperelastic constitutive model was employed to describe the rubber matrix elastic behavior. Based on the tensile experiment results, the material parameters were obtained as follows: μ1=-8.116, α1=1.975, μ2=5.725, α2=2.388, μ3=5.467, and α3=-4.495. For the aramid fibers, the elastic modulus is 132 GPa, the fiber diameter is 0.012 mm, and the studied fiber length is 1.5 mm. The difference in this study was the composite fiber distribution, which was established with an anisotropic morphology. Still, a perfect interface and straight morphology for the aramid fibers at the initial deformation were supposed. For fiber length lf=1.5 mm and fiber volume fraction vf=0.0049, using the anisotropic FE model, the stress–strain response of the composites with different perfectly aligned angles was obtained as shown in Figure 10. The curve shapes with blue scatters shown in Figures 10(a) and 10(b) are similar to those of Ex and Ey of the linear material shown in Figure 6, respectively. The minimum longitudinal modulus shown in Figure 11 also occurs when θ=60°. The stress–strain responses of the composites with different variancesσ2 from the numerical results are shown in Figure 11, which corresponds to Figure 8 for the linear material.Figure 10
Response of stress–strain relationship of composites with an aligned angleθ of fibers, lf=1.5 mm, and vf=0.0049.
(a)
Stress–strain response, when the composites were loaded along the x direction (b)
Stress–strain response, when the composites were loaded along the y directionFigure 11
Response of the stress–strain of composites with differentvariancesσ2 from the numerical results: lf=1.5 mm, vf=0.0049.
(a)
Stress–strain response when the composites were loaded along the x direction (b)
Stress-strain response when the composites were loaded along the y direction


## 6.2. Determination of the Orientation Distribution

For the hyperelastic material, the constitutive behavior is always expressed in the form of strain energyW. The hyperelastic constitutive model of an isotropic material has been comprehensively studied [27, 28]. However, for an anisotropic material, many terms remain to be explored [29–31]. To describe the orientation degree of the fibers as the isotropic factor α12 of the linear material, a variable, β=Wx/Wy was introduced, where Wx and Wy are the strain energy of the elastic response when the tensile load is the along x and y directions, respectively. According to the numerical results, for fiber length lf=1.5 mm and fiber volume fraction vf=0.0049, the relationships between β and Wx and Wy are shown in Figure 12. It can be seen that when β is near 1.0, Ex and Ey coincide.Figure 12
Relationship betweenWx (Wy) and the isotropic factor β for a rubber composite reinforced with aramid fibers: lf=1.5 mm, vf=0.0049.From the effect ofσ2 on β as shown in Figure 13, according to the experimental result β=1.09, the corresponding σ2 of the fiber orientation distribution was determined. σ2, which was approximately 1400, is the most reasonable value. Similarly, σ2 for lf=1.5 mm and vf=0.0146 was determined to be approximately 3000. The experimental results [12] and the corresponding stress–strain curves of the determined σ2 are compared in Figure 14. There are good agreements, indicating that the proposed method predicting the fiber orientation distribution is applicable.Figure 13
Relationship between the isotropic factorβ and variancesσ2 for a rubber composite reinforced with aramid fibers: lf=1.5 mm, vf=0.0049.Figure 14
Comparison between the experimental and numerical results of the model with the determinedσ2.
(a)
l f=1.5 mm, vf=0.0049 (b)
l f=1.5 mm, vf=0.0146


## 7. Conclusion

In this study, an FE model with an embedded technique was employed to predict the anisotropy of polymer composites reinforced with short aramid fibers. Using this FE method, one can obtain the RVE with high aspect ratio fibers in the composites. To obtain the elastic properties of the real material, the RVEs of three types of fiber orientation distributions were established using Python language: perfectly aligned, normally distributed, and randomly distributed fibers. The five elastic properties, tensile elastic moduliEx and Ey, shear elastic modulus Exy, and Poisson ratios vx and vy, were obtained by different elastic deformations with corresponding boundary conditions. For different fiber orientation distributions, the numerical results were compared to those of the Mori–Tanaka model and found to agree well.Based on the aforementioned conclusion, the FE model was applied to predict the fiber orientation distribution. As a nonlinear material, an isotropic factor in the form of strain energy was proposed to present the deviation degree from isotropy. According to the relationships between the isotropic factor and variance, the most reasonable variances used to describe the fiber orientation distribution in the real composite were determined. Comparison between the experimental results and numerical results of the determined variance indicated good agreements. The FE model used in this study is a simple and convenient method. The numerical prediction for a composite with a high-aspect-ratio fiber has its advantages. In addition, it can be applied to simulate the large deformation of a hyperelastic material that is difficult to realize using the traditional method because of the easy occurrence of numerical nonconvergence. The proposed FE model aids in predicting the anisotropy of short-fiber-reinforced composites and has a shorter experimental period and lower cost.



---

*Source: 5484675-2019-07-02.xml*