Micromechanical properties of unidirectional composites filled with single and clustered shaped fibers

1 Introduction

Apart from classic circular fibers, polygonal fibers such as steel fibers [1] and glass fibers [2] are also of interest in the context of fiber-reinforced composites. It is known that the noncircular cross sections of fibers are likely to generate high stress concentration in composites in response to external loading and thus to affect stress transfer from matrix to fibers. Besides, the cluster distribution of fillers in the matrix is usually unavoidable because of manufacturing conditions for composites [3], [4]. It is therefore essential to understand the effects of the cross-sectional shape and clustering of fiber on the overall mechanical properties of composites to improve their macroscopic and microscopic material behavior.

In the past decade, research efforts have been made in the subject of inclusion shape. For instance, the mean-field theory was developed to determine the overall stress-strain relation of two-phase composites caused by different inclusion shapes [5]. Wu established a self-consistent model to calculate the elastic moduli of a two-phase material for investigating the influence of inclusion shape [6]. Qin and Swain [7] developed a dilute micromechanics model – a circular hollow cylinder filled with liquid or gas phase, which is surrounded by two circular cylindrical shells, a thin shell and a matrix phase. Tsukrov and Novak [8] studied the effect of irregularly shaped inclusions to the effective elastic moduli of two-dimensional composites by the conformal mapping procedure. However, it is very difficult to obtain closed-form solutions for elastic media containing noncircular cross-sectional inclusions. Alternatively, approximation solutions by numerical computational methods such as the finite element method (FEM) [9], [10], [11], [12] and the boundary element method (BEM) [13], [14], [15] can be obtained. Some investigators have studied the influence of inclusion shape on the effective Young’s moduli of fiber/particle reinforced composites by computational micromechanics. For example, Brockenbrough et al. [16] used the FEM tool to study the effects of fiber distribution and fiber’s cross-sectional geometry with circular, hexagonal, and square shapes in metal-matrix composites. Huang and Talreja performed finite element analysis in a representative volume cell including void microstructure to investigate the effect of void geometry on the effective elastic constants of unidirectional fiber reinforced composites [17]. The void can be viewed as a special inclusion with zero Young’s moduli. Banks-Sills et al. [18] investigated the influence of four ceramic particles with spherical, cylindrical, cubic, and rectangular parallelepiped shapes on elastic properties of metal-matrix composites by using FEM. Beyerlein et al. [19] analyzed the influence of bone-shaped fiber in short-fiber composites by using FEM. Herráez et al. [20] developed computational micromechanics with FEM for analyzing the effect of fiber shape on the effective mechanical properties of composites. In their study, circular and noncircular fibers with lobular, polygonal (three and four edges), and elliptical shapes were considered. Besides, different to the FEM based on the domain discretization mentioned previously, the BEM based on boundary-only discretization has also been developed to demonstrate the influence of inclusion shape. Dong established the submodel boundary integral formulations for the inclusion and the matrix material phase to determine the effective elastic properties of doubly periodic composites filled with circular hole, square inclusion, and hexagonal inclusion [21]. Chen and Liu [22] applied the BEM to model multiple cells of composites filled with coated circular fibers. Pan et al. [23] applied the BEM for the analysis of three-dimensional composite laminates with holes. Moreover, other numerical methods such as radial basis collocation methods with radial integration technique [24] for strains approximation in four cross-sectional inclusions [25] and the Voronoi cell FEM [26] were developed for composite analysis. However, these studies seldom involved the systematic and comprehensive analysis of inclusion/fiber shapes. Besides, in the context of filler clustering analysis, most of the work was restricted to composites filled with spherical particles. For example, a three-dimensional clustering microstructure model with multiple simple particles in composites was established to study the effects of the microstructure of the composite, such as particle size, shape, and distribution [27]. Later, the effect of sphere particle clustering in polymer composites was investigated by Segurado et al. [3], [28]. More recently, the thermal and elastic behavior of composite with multiple circular fibers or holes was studied by special hybrid finite element models to achieve the purpose of mesh reduction [29], [30], [31], [32], [33]. Besides, the particle-based discrete element method was developed for microbond and random arrangement of fiber-reinforced composites [34], [35].

In the present paper, the comprehensive contribution of variously shaped fibers (inclusions) and their clusters to the effective elastic moduli of two-dimensional unidirectional composites is studied in the representative unit cell of periodic composites, followed by micromechanical FEM simulation coupling the homogenization technology [21], [36], [37]. The cross-sectional shapes of polygonal inclusions are triangular, square, pentagonal, hexagonal, octagonal, and circular. Here, the circular boundary can be viewed as a special polygon with infinite number of edges. To determine the overall elastic constants of composites, the representative unit cell subject to periodic conditions is used, which are given to keep the edges of the cell straight after the deformation, and then is solved by the FEM to calculate the stress distribution along the cell edge. With the obtained stress results, the five elastic constants in the orthotropic elastic constitutive equation can be evaluated by solving three groups of linear equations of system in two unknowns. Then systematic and quantitative relationships are established between the inclusion volume concentration and the elastic moduli for various shapes and cluster pattern of inclusions.

The paper is organized as follows: In Section 2, formulations of mechanical properties of periodic unidirectional composites under consideration are described, and then in Section 3, solution procedure on the effective elastic properties of the composite are presented by micromechanical finite element analysis. In Section 4, numerical results are given and discussed to demonstrate the effect of inclusion shape and clustering. Finally, some conclusions are drawn in Section 5.

2 Mechanical properties of composites

For unidirectional composites with multiple polygonal inclusions and in periodic square arrangement, as displayed in Figure 1, we assume that both inclusions and matrix are isotropic elastic, and the corresponding elastic parameters are denoted by E^(I), v^(I), and E^(M), v^(M). In addition, perfect bonding at the interface between the inclusions and the matrix is assumed here.

Figure 1:

Sketch of (A) composite with doubly periodic inclusions, (B) representative unit cell, and (C) equivalent medium.

For such periodic structure, a representative unit cell can be isolated from the composite, and the unit cell has the same elastic constants and inclusion volume fraction as the composite under consideration [38]. Figure 1B shows a unit cell with edge length 2l and 2h along the coordinate directions x₁ and x₂. On the cell outer boundary, the suitable periodic boundary conditions meeting the requirement of straight-line constraint after deformation are applied for different loading conditions, including uniform tension and shear deformation [21], [39], as displayed in Figure 2:

Figure 2:

Periodic boundary conditions along the cell boundary.

(A) Biaxial tension and (B) shear.

1. Biaxial uniform tension

   (1)u1(l, x2)=−u1(−l, x2)=u¯l1t2(l, x2)=t2(−l, x2)=0u2(x1, h)=−u2(x1, −h)=u¯h2t1(x1, h)=t1(x1, −h)=0

2. In-plane shear

   (2)u2(l, x2)=−u2(−l, x2)=u¯l2t1(l, x2)=t1(−l, x2)=0u1(x1, h)=−u1(x1, −h)=u¯h1t2(x1, h)=t2(x1, −h)=0

where u̅_l1, u̅_l2, u̅_h1, and u̅_h2 are unknown constant displacement constraints to be determined together with other unknown variables discussed in the following paragraph. These unknown constants are meant to keep the opposite edges straight after the deformation to represent the constraint of the neighboring cells to the one under study.

Because the composite containing periodic array of inclusions can be modeled as a homogeneous orthotropic medium, shown in Figure 1C, with certain effective moduli that describes the average material properties of the composite, the corresponding orthotropic constitutive relationship associated with the average stress σ¯ij and the strain ε¯ij over the unit cell can be written as [40], [41]

where E̅_i (i=1, 2) is the effective orthotropic elastic modulus constant along the coordinate axis x_i. ν¯ij is the Poisson’s ratio corresponding to a contraction in the direction j when an extension is applied in the direction i. G̅_ij is the shear modulus in the direction j on the plane whose normal is in the direction i. According to the material symmetry, one has

We know that a uniform stress and strain state will exist under uniform loading in a homogeneous material. To determine the elastic constants E̅₁ and ν¯12 in the effective homogeneous orthotropic medium, a uniform tension along the x₁ direction is considered. The average stress state in this case can be described by

and the average strain state is given by the conventional definition of strain [38],

where σ₀ is a known stress value.

Substituting the stress and strain states in Equations (5) and (6) into Equation (3) yields

Similarly, to determine the effective elastic properties E̅₂ and ν¯21, the uniform tension along the x₂ direction can be applied and the corresponding average stress and strain states are, respectively, written by

and

Substituting them into the orthotropic constitutive Equation (3) yields

Finally, to determine the effective shear modulus G̅₁₂, one needs to apply the uniform shear loading in the x₁−x₂ plane, that is,

then the corresponding average strain state can be represented by

The substitution of Equations (11) and (12) into Equation (3) gives

3 Micromechanical finite element analysis

From the procedure previously mentioned, it can be seen that the specified boundary displacements are involved to evaluate the effective elastic properties of the composite for the three loading cases. In this section, the finite element micromechanical analysis is performed to determine the specified values of boundary displacement constraints according to the applied stress component.

First, according to the superposition principle of linear elasticity, the micromechanical model established by applying the boundary conditions [Equation (1)] on the cell boundary shown in Figure 2A can be decomposed into two submodels: submodel 1 and submodel 2 (see Figure 3). Let’s assume that the stress components related to submodel 1 under the unit displacement condition u̅_h2=1 and submodel 2 under the unit displacement condition u̅_l1=1 are represented as σ11(1),σ22(1),τ12(1) and σ11(2),σ22(2),τ12(2), respectively, where the upper index denotes the number of the submodels.

Figure 3:

Submodels 1 and 2 corresponding to the biaxial tension deformation.

(A) Submodel 1 and (B) submodel 2.

If the uniaxial uniform tension along the x₁ direction (σ¯11=σ0,   σ¯22=0,   σ¯12=0) is specified, the equilibrium relationship of the cell can be written as [21], [39]

in which the quantities in the square brackets represent the average values of the stress components on the top and right surfaces, which can be evaluated by simple trapezoidal integral.

Solving the linear system of equations [Equation (14)], we can obtain displacements u̅_l1 and u̅_h2 corresponding to the given value σ₀. Subsequently, we can obtain the effective elastic constants E̅₁ and ν¯12 by using Equation (7).

In a similar manner, the equilibrium relationship of the cell subject to the uniform tension along the x₂ direction (σ¯11=0,   σ¯22=σ0,   σ¯12=0) can be written as [21], [39]

from which we can evaluate one set of unknown displacements u̅_l1 and u̅_h2, and further the effective properties E̅₂ and ν¯21 can be calculated using Equation (10).

Second, the micromechanical model established by applying the in-plane shear boundary conditions [Equation (2)] on the outer boundary of the cell shown in Figure 2B can be decomposed into two submodels, submodel 3 and submodel 4 (see Figure 4), provided that the superposition principle of linear elasticity is used. Also we assume that τ12(3) and τ12(4) are, respectively, shear stress components under the unit displacement condition u¯h1=1 in submodel 3 and u¯l2=1 in submodel 4.

Figure 4:

Submodels 3 and 4 corresponding to the shear deformation.

(A) Submodel 3 and (B) submodel 4.

Then the equilibrium relationship of the cell subject to the uniform shear (σ¯11=0,   σ¯22=0,   σ¯12=τ0) can be written as [21], [39]

from which we can obtain the third set of periodicity displacement constraints u̅_h1 and u̅_l2. Furthermore, the effective shear constant G̅₁₂ can be calculated by Equation (13).

From the previously mentioned procedure for predicting the elastic constants of the composite, it is found that the four submodels under specified unit displacement constraints need to be solved for stress solutions along the top and right surfaces of the cell. In this study, the finite element simulation is performed for solving the four submodels under unit displacement constraint because of its advantages of flexibility in inhomogeneity, complex geometry, and boundary conditions.

4 Results and discussion

During the micromechanical computation, each inclusion is assumed to be isotropic and has same elastic modulus E^I=10 and Poisson’s ratio ν^I=0.3. Besides, the elastic modulus and the Poisson’s ratio of the matrix medium are taken as E^M=1 and ν^M=0.3, respectively. Moreover, for the sake of simplification, the square packing is set by l=h=1. In addition, both the matrix and the inclusions are modeled by the 8-node quadratic quadrilateral isoparametric plane stress elements in FEM simulation. In the practical computation, the mesh size is set to be small enough so that the effect of stress concentration of the polygonal cross-sectional inclusion and the heterogeneous nature of microstructure can be captured to acquire accurate and stable results. In the study, the global element size can be set to be as small as 0.01 times the side length of inclusion, and the maximum numbers of nodes are 42,749, 31,008, 36,603, 33,902, 25,308, and 42,783 for triangular, square, pentagonal, hexagonal, octagonal, and circular inclusions, respectively.

4.2 Single inclusion

After the verification of the present micromechanical finite element model, six types of regular inclusions with triangular, square, pentagonal, hexagonal, octagonal, and circular cross sections are studied in this subsection. The side length of the inclusion can be evaluated by means of the specified value of the inclusion volume concentration. Generally, the geometrical center of the inclusion is set to coincide with that of the square cell domain. However, for the triangular and pentagonal shapes, the midpoint of height is set to coincide with that of the square domain to achieve large fiber volume concentration. Moreover, it is seen that because of the limitation of the inclusion shape, the maximum values of inclusion volume concentration are different for inclusions with different shapes. For example, the maximum theoretical volume concentration of the triangular inclusion is 43.3%, whereas for the circular inclusion, the maximum theoretical value of the volume concentration reaches 78.5%.

Because of mesh dependence of the FEM, the convergence test is first performed by considering the triangular inclusion under the condition of 10% fiber volume concentration. Figure 5 displays the variation of von-Mises stress at the midpoint of the upper edge of the cell as the number of nodes changes. It is clear that the result becomes relatively stable when the mesh density increases, and a smaller relative error (0.02%) appears when the number of nodes increases from 6537 to 10,028. Here, the result with 10,028 nodes is referred for the evaluation of relative error.

Figure 5:

Convergence illustration for triangular fiber with 10% fiber volume concentration.

Figure 6 displays the von Mises stress distribution in the cell for each inclusion for the case of submodel 1, which represents the uniform tension along the axis 2-direction. The inclusion volume concentration is taken to be 30%. It is noted that the angular inclusions including triangular, square, pentagonal, hexagonal, and octagonal inclusions are under a much higher stress than the circular inclusion, indicating more load transfer to the angular inclusions. The stress in the circular inclusion is quite uniform, whereas that in the angular inclusion is not. Moreover, it is found from Figure 6 that the stress concentration becomes weak as the number of edges of the polygonal inclusion increases. The effective elastic properties of the homogeneous orthotropic medium for various inclusion volume concentrations are shown in Figures 7–10, in which results for different shaped inclusions are given for comparison. It is observed that the two effective elastic moduli, E₁ and E₂, and the effective shear modulus, G₁₂, have an evident increase, as the inclusion volume concentration increases for each inclusion, while the Poisson’s ratio μ₂₁ decreases. More importantly, it can be seen from Figures 7 to 10 that there is no significantly difference of effective elastic constants E₁, E₂, and G₁₂ for different shape of inclusions, except for the triangular inclusion. Although, the difference of Poisson’s ratio caused by geometrical shape variation is significant. Also, it is found that the triangular cross section brings the largest E₁ and G₁₂ when the fiber volume concentration is 40%. The main reason is that the triangular cross section has the largest sharpness, and its top vertex is too close to the upper edge of the cell, so the edge effect becomes strong and leads to higher elastic modulus E₁. In addition, it is observed that the circular inclusion can produce nearly isotropic mechanical properties, among all polygonal inclusions considered in the study.

Figure 6:

The von Mises stress distribution of six types of inclusion for submodel 1 with 30% inclusion volume concentration.

Figure 7:

Effective elastic modulus E₁ of the composite for different shaped inclusions.

Figure 8:

Effective elastic modulus E₂ of the composite for different shaped inclusions.

Figure 9:

Effective Poisson’s ratio µ₂₁ of the composite for different shaped inclusions.

Figure 10:

Effective shear modulus G₁₂ of the composite for different shaped inclusions.

4.3 Clustered inclusions

In this subsection, the clusters of hexagonal and circular inclusion are investigated for two different inclusion volume concentrations, 10% and 40%. The cluster pattern is displayed in Figure 11. In the figure, the clustering degree parameter λ is introduced to depict the distance between inclusions. For the case of hexagonal inclusion, Table 2 indicates the changes of overall mechanical properties of the composite for different λ and v_f, respectively. For the small filler volume concentration, i.e. v_f=10%, it is found that there is no significant change of mechanical properties when varying normalized clustering degree parameters λ/a. Moreover, compared with the results of single hexagonal inclusion, the difference is almost negligible. When the filler volume concentration increases to 40%, a large value for clustered inclusions, the clustering effect of hexagonal inclusion becomes evident for two different values of clustering degree ratio (λ/a). For the Poisson’s ratio μ₁₂, there is 22% deviation from that for the single hexagonal inclusion. Similar conclusions are drawn for clustered circular inclusions (see Table 3). Moreover, it is clear from Tables 2 and 3 that the circular inclusion cluster has more averaged performance than the hexagonal inclusion cluster.

Figure 11:

Schematic illustrations of the clustered models for the hexagonal and circular inclusions packed in square arrays.

Table 2:

Results of mechanical properties for the clustered hexagonal inclusions.

λ/a	E1	E2	μ12	μ21	G12
vf=10%
1.50	1.1449	1.1428	0.3027	0.3021	0.4420
3.00	1.1425	1.1425	0.3018	0.3018	0.4398
Single inclusion	1.1502	1.1502	0.2970	0.2970	0.4356
vf=40%
1.40	1.9765	1.9457	0.2864	0.2820	0.6991
1.52	1.9520	1.9105	0.3048	0.2986	0.7193
Single inclusion	1.9462	1.9439	0.2498	0.2496	0.6384

Table 3:

Results of mechanical properties for the clustered circular inclusions.

λ/a	E1	E2	μ12	μ21	G12
vf=10%
1.50	1.1463	1.1463	0.30186	0.30186	0.44103
3.00	1.1397	1.1397	0.30263	0.30263	0.43908
Single inclusion	1.1467	1.1466	0.29803	0.29799	0.43477
vf=40%
1.40	2.0621	2.0621	0.2553	0.2553	0.6832
1.52	1.9009	1.9009	0.3246	0.3246	0.7412
Single inclusion	1.9222	1.922	0.2536	0.2536	0.6324

5 Conclusions

The effective in-plane mechanical behavior of unidirectional composite filled with inclusions has been simulated by means of computational micromechanics. In the present computational micromechanical modeling, the comprehensive contributions of variously shaped and clustered polygonal inclusions are taken into consideration using finite element analysis coupling with the homogenization technology based on mean-field theory, which is introduced to meet the straight-line boundary constraints during the deformation of the unit cell of composite. Initially, the cross-sectional shapes of inclusions are taken to be triangular, square, pentagonal, hexagonal, octagonal, and circular. It is found that for all shaped inclusions, the inclusion volume concentration plays a key role to the stiffness parameters of composite. The elastic moduli and the shear modulus increase with the increase of inclusion volume concentration (v_f). However, the Poisson’s ratio μ₁₂ decreases with the increase of v_f. Moreover, circular inclusions show the best averaged performance than angular polygonal inclusions. It is also observed that there is just slight derivation of elastic properties for them, except for the Poisson’s ratio, despite the dramatically changed local stress distribution caused by different shaped inclusions. In addition, the effects of clustered hexagonal and circular inclusions on the effective elastic properties of composite materials are demonstrated. It is observed that among all effective elastic properties, the Poisson’s ratio always produces larger derivation to that in the single inclusion case for specified two inclusion volume concentrations. Besides, because the method for calculating the effective constants is based on linear superposition, it is currently not suitable for hyperelastic soft composites [42].