Mechanical Homogenization of Transversely Isotropic CNT/GNP Reinforced Biocomposite for Wind Turbine Blades: Numerical and Analytical Study

Abstract

One of the biggest problems facing the use of carbon nanotubes in reinforced composites is agglomeration within the matrix phase. This phenomenon—caused by Van der Waals forces—leads to dispersion problems and weakens the properties of the composites. This research presents a multi-stage homogenization approach used to investigate the influence of the aspect ratio, volume fraction, and agglomeration of the nanofillers on the effective mechanical properties of a polymer biocomposite containing randomly dispersed carbon nanotubes and graphene nanoplatelets. The first stage consisted in evaluating the properties of the reinforced polymers by the CNT/GNP. The second step consisted in combining the reinforced polymers with different natural and synthetic unidirectionally oriented fibers. It was found that agglomeration has a huge influence on the mechanical properties of the composite. The novelty of this work consisted of the consideration of the parameters influencing the elastic properties using different micromechanics approaches and numerical techniques.

1. Introduction

Several factors are considered while selecting materials for wind energy turbines, mainly stiffness and density. High rigidity prevents blades from bending, and low density allows for lighter blades, which mean higher turbine efficiency [1,2,3,4]. Traditionally, turbines were made from metal, although these could not resist the permanent stress cycles. Fiber-reinforced composites (FRC) presented a decent replacement for metallic blades, with glass and carbon fibers being the most commonly used fibers in the industry because of their superior properties; however, their unrecyclable nature and price (about $10.0 per lb. for carbon fibers) [5,6] are the biggest inconveniences that need to be overcome. Biocomposites, which consist of natural fibers such as flax, sisal or alfa immersed in a polymer, offer the best solution to the waste and price issue as they can behave in a mechanically similar way to synthetic fibers [7]. As they are a combination of two materials or more, biocomposites are considered to be heterogeneous materials, therefore their effective properties need to be evaluated for a more comprehensive understanding of these materials.

Recently, nanofillers have made considerable advances in the fiber-reinforced composites industry because of their ability to significantly improve mechanical, electrical, and thermal properties [8,9,10]. Since their discovery by Iijima [11], CNT’s exceptional structure and physical properties, which include their remarkable conductivity, have attracted enormous interest in scientific research [7,12]. For practical applications, it is necessary to establish structure–property relationships for a reliable design of these materials, which has inspired several research works to evaluate physical properties of CNT-reinforced composites, such as the study Thostenson et al. [13] which investigated the influence of an MWCNT-reinforced polystyrene composite and found that the mechanical properties were significantly improved. Bonnet et al. [14] studied thermal properties of single-walled carbon nanotube (SWNT)/polymethyl methacrylate composite thick films. The composite’s thermal conductivity was improved by 55% and its electrical conductivity was increased by orders of magnitude with the addition of 7% SWNT to the polymer matrix. Fidelus et al. [15] examined the thermo-mechanical properties of epoxy-based nanocomposites based on low weight fractions (from 0.01 to 0.5 wt.%) of randomly oriented single- and multi-walled carbon nanotubes and observed important improvements in the elastic modulus and tensile strength. On the other hand, graphene nanoplatelets (GNP), which are composed of several 2D layers of graphene, can provide better mechanical properties than CNTs [16,17]. In Rafiee et al.’s investigation [18], the mechanical characteristics of epoxy nanocomposites containing single-walled carbon nanotubes, multi-walled carbon nanotubes, and graphene platelets were compared. Young’s modulus, ultimate tensile strength, fracture toughness, fracture energy, and the material’s resistance to fatigue crack propagation were the mechanical parameters that were measured. The findings showed that carbon nanotube additives drastically underperform graphene platelets in terms of performance. In Gao’s [17] study, a micromechanical model was proposed to predict the thermomechanical properties of graphene reinforced metal matrix. Results of the study showed that the yield strength of the nanocomposites can be dramatically enhanced by adding a small volume percent of GNP, for example, an increment of 92.0% by 0.7 v% GNP/Al2024, 128.0% by 1.5 v% GNP/Al and 78% by 2.5 v% GNP/Cu, respectively.

However, the main challenge facing the use of CNTs and GNPs is their poor dispersion, which leads to the formation of agglomerated regions due to the high length-to-diameter ratio and the small elastic modulus in the radial direction of the nanofillers that generate van der Waals attractive forces the main cause of this phenomenon. Several research investigated the effect of the agglomeration of nanofillers. Narh et al. [19] studied the effect of agglomerated and de-agglomerated MWCNT on the thermal and mechanical characteristics of polyethylene oxide and found that nanocomposites containing deagglomerated nanotubes appear to show marginal improvement in toughness over those containing agglomerated MWCNTs. Alian et al. [20] proposed a multi-scale model based on molecular dynamics and micromechanical models to determine the effect of waviness and agglomeration of CNTs on the bulk elastic properties of epoxy composites and found that both limited their reinforcement effect. Daghigh et al. [21] investigated the nonlocal bending and buckling behaviors of agglomerated CNTRC nanoplates numerically and concluded that the elastic properties of the nanoplates are extremely affected by any type of agglomeration of CNTs, and also that ignoring clustering in the numerical model can lead to significant errors. Shi et al. [22] examined the effect of the agglomeration and found that clustering significantly reduces the stiffening effects of CNTs by 80% in the most severe cases. In Ji’s study [23], GNP agglomerations were investigated theoretically and it was shown that any type of agglomeration damages the stiffening effect. Maghsoudlou et al. [9] studied the elastic modulus of single-walled carbon nanotubes/epoxy composites computationally and experimentally and found that Young’s modulus increased while increasing the SWCNT weight fraction until 0.3 wt% and dropped sharply. This was explained by the agglomeration effect on the nanofillers. García-Macias et al. [8] presented a comparison of mean-field homogenization approaches for randomly and aligned oriented CNT-reinforced polymers and found that agglomeration caused detrimental effects for the mechanical properties, and more severely for the aligned CNTs. From this short review, the importance of, and the need for, considering the agglomeration effect in the modelling of CNT/GNP-reinforced composites can be concluded. For this purpose, we proposed a two-stage homogenization procedure (Figure 1) considering the agglomeration of the nanofillers. First, we evaluated the mechanical properties of the reinforced polymer by CNTs and GNPs using the Mori–Tanaka analytical approach and the numerical study was conducted using Digimat MF and FE. The obtained results were used in a second homogenization using Chamis’s, Hashin–Rosen and Halpin–Tsai approaches for isotropic alfa and glass fibers. For transversely isotropic carbon fibers, we used Hahn’s and Halpin–Tsai approaches. Results were compared with the numerical results calculated using Digimat.

2. First Homogenization

2.1. Mathematical Model

2.1.1. CNT Inclusion

CNTs can be employed as a nanofiller to enhance the mechanical properties of a polymeric matrix. For this purpose, the matrix is supposed as isotropic and characterized by the Young’s modulus $E_m$ and the Poisson’s ratio ${\nu }_m$. The CNT inclusion is considered as transversely straight with an infinite aspect ratio. In this work, the Eshelby–Mori–Tanaka homogenization scheme (Figure 1) was used to determine the elastic properties of the CNT-reinforced matrix [1,24,25]. For a CNT embedded in the composite, it is related to a local coordinate system (O, x₁, x₂, x₃), where the x₃-axis lies along its longitudinal direction and the (x₁, x₂) plane coincides with its transversely isotropic plane. The carbon nanotube CNT is considered as a rolled graphene sheet and is only strong in its length direction as seen in Figure 2a. Comparably, a single-layered graphene sheet GNP is strong along all its in-plane directions, as seen in Figure 2b, and offers a large surface/volume ratio [26]. When the nanofillers are perfectly randomly distributed in the matrix, the overall elastic property of the reinforced composite is isotropic as an average effect of orientation.

It should be noted that the superscript (*) is employed to identify the property related to the reinforced matrix. The bulk $K_m^{*}$ and shear modulus $G_m^{*}$ of the reinforced composite can be given by Equations (1) and (2).

$$K_m^{*}=K_m+\frac{V_r \left({\delta }_r-3K_m{\alpha }_r\right)}{3 \left(V_m+V_r{\alpha }_r\right)}$$

(1)

$$G_m^{*}=G_m+\frac{V_r \left({\eta }_r-2G_m{\beta }_r\right)}{2 \left(V_m+V_r{\beta }_r\right)}$$

(2)

where

$${\alpha }_r=\frac{3 \left(K_m+G_m\right)+k_r-l_r }{3 \left(G_m+k_r\right)}$$

(3)

$${\beta }_r=\frac{1}{5}\left[\frac{4G_m+2k_r+l_r }{3\left(G_m+k_r\right)}+\frac{4G_m }{G_m+p_r}+\frac{2G_m\left(3K_m+G_m\right)+2G_m\left(3K_m+7G_m\right)}{G_m\left(3K_m+G_m\right)+m_r\left(3K_m+7G_m\right)}\right]$$

(4)

$${\delta }_r=\frac{1}{3}\left[n_r+2l_r+\frac{\left(2k_r+l_r\right)\left(3K_m+2G_m-l_r\right)}{G_m+k_r}\right]$$

(5)

$$\begin{gathered} {\eta }_r=\frac{1}{5}\left[\frac{2}{3}\left(n_r-l_r\right)+\frac{8G_m{p}_r}{G_m+p_r}+\frac{2\left(k_r-l_r\right)\left(2G_m+l_r\right)}{3\left(G_m+k_r\right)} \\ +\frac{8m_r{G}_m\left(3K_m+4G_m\right)}{3K_m\left(m_r+G_m\right)+G_m\left(7m_r+G_m\right)}\right] \end{gathered}$$

(6)

It should be pointed out that $K_m$ and $G_m$ specify the bulk and the shear moduli of isotropic matrix, which can be evaluated from Equations (7) and (8).

$$K_m=\frac{E_m}{3\left(1-2{\nu }_m\right)}$$

(7)

$$G_m=\frac{E_m}{2\left(1+{\nu }_m\right)}$$

(8)

The reinforced matrix is evidently isotropic. Therefore, the Young’s modulus $E_m^{*}$ and the Poisson’s ratio ${\nu }_m^{*}$ are required to define its properties. They can be computed in Equations (9) and (10).

$$E_m^{*}=\frac{9K_m^{*}{G}_m^{*}}{3K_m^{*}+G_m^{*}}$$

(9)

$${\nu }_m^{*}=\frac{9K_m^{*}-2G_m^{*}}{6K_m^{*}+2G_m^{*}}$$

(10)

Finally, the density ${\rho }_m^{*}$ of the reinforced matrix can be defined by using the rule of mixture as

$${\rho }_m^{*}=\left({\rho }_r-{\rho }_m\right)V_r+V_m$$

(11)

2.1.2. GNP Inclusion

For the case of the nanoplatelet-reinforced composites, the expressions for effective bulk modulus and shear modulus are in Equations (12) and (13).

$$K_m^{*}=K_m+\frac{V_r \left({\delta }_r-3K_m{\alpha }_r\right) }{3\left(c_m+V_r{\alpha }_r\right)}$$

(12)

$$G_m^{*}=G_m+\frac{V_r \left({\eta }_r-2G_m{\beta }_r\right) }{2\left(c_m+V_r{\beta }_r\right)}$$

(13)

where

$${\alpha }_r=\frac{3k_m+2n_r-2l_r }{3n_r}$$

(14)

$${\beta }_r=\frac{4G_m+7n_r+2l_r }{15n_r}+\frac{2G_m }{5p_r}$$

(15)

$${\delta }_r=\frac{3k_m\left(n_r+2l_r\right)+4\left(k_r{n}_r-l_r{}^2\right) }{3n_r}$$

(16)

$${\eta }_r=\frac{2}{15}\left(k_r+6m_r+6G_m-\frac{l_r{}^2+2G_m{l}_r }{n_r}\right)$$

(17)

Then, the Young’s modulus $E_m^{*}$ and the Poisson’s ratio ${\nu }_m^{*}$ for the GNP-reinforced matrix can be defined from Equations (9) and (10).

2.2. Numerical Model

2.2.1. RVE Geometry

In this research, an RVE method was used to predict the overall properties of the reinforced polymers. The different microstructures were created using Digimat FE.; Digimat is a commercial code from e-Xstream Engineering that includes a number of micromechanics tools for both homogenization and dehomogenization. In our study, Digimat MF (semi-analytical approach) and Digimat FE (numerical approach) were used.

For the smooth running of the calculation, some assumptions are made:

-

The CNT inclusions are considered to have an ellipsoidal (prolate) shape. GNP inclusions are considered as disc-shaped spheroidal inclusions with aligned and random orientations.

-

The phases are considered as perfectly bounded.

The dimension of the representative volume element (RVE) is considered as (1, 1, 1) units, as seen in Figure 3, and for each one, a uniaxial loading strain test was performed to evaluate the elastic properties using the Mori–Tanaka homogenization scheme (Figure 1).

2.2.2. Materials Properties

The elastic parameters of CNT and GNP are taken from [27] and written in terms of Hill’s elastic moduli in

Table 1. Mechanical properties of nanofillers in form of Hill’s elastic moduli.

	${\mathit{k}}_{\mathit{r} }\left(\mathrm{GPa}\right)$	${\mathit{l}}_{\mathit{r}} \left(\mathrm{GPa}\right)$	${\mathit{m}}_{\mathit{r}} \left(\mathrm{GPa}\right)$	${\mathit{n}}_{\mathit{r}} \left(\mathrm{GPa}\right)$	${\mathit{p}}_{\mathit{r}} \left(\mathrm{GPa}\right)$
CNT	536	184	132	2143	791
GNP	850	6.8	369	102,000	102,000

.

Analogously, the engineering constants $E_{11}^r, E_{22}^r, E_{33}^r,{\nu }_{12}^r,{\nu }_{13}^r, {\nu }_{23}^r, G_{12}^r,G_{13}^r,G_{23}^r$ are employed to define the elastic behavior of the equivalent continuum model from Equations (18)–(23). Then, the five transversely isotropic constants are calculated in

Table 2. Mechanical properties of CNT and GNP in terms of engineering constants.

	${\mathbf{E}}_{11}^{\mathbf{r}}$ $\left(\mathbf{GPa}\right)$	${\mathbf{E}}_{22}^{\mathbf{r}}={\mathbf{E}}_{33}^{\mathbf{r}}$ $\left(\mathbf{GPa}\right)$	${\mathsf{\nu }}_{23}^{\mathbf{r}}$	${\mathsf{\nu }}_{12}^{\mathbf{r}}={\mathsf{\nu }}_{13}^{\mathbf{r}}$	${\mathbf{G}}_{23}^{\mathbf{r}}$ $\left(\mathbf{GPa}\right)$	${\mathbf{G}}_{12}^{\mathbf{r}}={\mathbf{G}}_{13}^{\mathbf{r}}$ $\left(\mathbf{GPa}\right)$
CNT	421.14	2079.8	0.17164	0.59522	132	791
GNP	1029.204	102,000	0.004	0.4	369	102,000

.

$$E_{11}^r=\frac{4m_r \left(k_r{n}_r-l_r{}^2\right)}{k_r{n}_r-l_r{}^2+m_r{n}_r}$$

(18)

$$E_{33}^r=E_{22}^r=n_r-\frac{l_r{}^2}{k_r}$$

(19)

$${\nu }_{13}^r={\nu }_{12}^r=\frac{l_r}{2k_r}$$

(20)

$${\nu }_{23}^r=\frac{n_r \left(k_r-m_r\right)-l_r{}^2}{n_r \left(k_r+m_r\right)-l_r{}^2}$$

(21)

$$G_{13}^r=G_{12}^r=p_r$$

(22)

$$G_{23}^r=m_r$$

(23)

The remaining engineering constants must be found through the following relations in Equation (23).

$$\frac{{\nu }_{ij}^r}{E_i^r}=\frac{{\nu }_{ji}^r}{E_j^r}, G_{ij}^r=G_{ji}^r \mathrm{for} i, j=1, 2, 3$$

(24)

The mechanical properties of polypropylene PP and unsaturated polyester UP polymers are described in

Table 3. Mechanical properties of polymers.

Materials	PP	UP
Model	Elastic	Elastic
Symmetry	Isotropic	Isotropic
Density ρ (g/cm3)	0.9	1.3
Young’s modulus E (GPa)	1.4	3.8
Poisson’s ratio v	0.45	0.42

[1].

The mechanical properties used in Digimat MF are depicted in

Table 4. Mechanical properties of CNT and GNP as inputs in Digimat MF.

Elastic Properties	CNT	GNP
Axial Young’s modulus (GPa)	2079.8	102,000
In plane Young’s modulus (GPa)	421.14	1029.204
In plane Poisson’s ratio	0.59522	0.4
Transverse Poisson’s ratio	0.17164	0.004
Transverse shear modulus (GPa)	791	102,000
In plane shear modulus (GPa)	132	369
Density (g/cm3)	1.2	2.2

, by using 10³ and 10⁻⁴ as aspect ratios for CNT and GNP, respectively.

2.3. First Homogenization Results

The elastic moduli result of ($E_m^{*}$, ${\nu }_m^{*}$) with Digimat MF and the Mori-Tanaka-Eshelby model were validated for CNT embedded in a PP polymer matrix in Figure 4a,b, to know the viability of the numerical model. Then, by just using the numerical results of Digimat MF, we compared the mechanical properties of CNT embedded in PP and UP, with GNP embedded in PP and UP as seen in Figure 5a,b.

As seen in Figure 4a, the two approaches show an acceptable agreement between the Mori–Tanaka scheme and Digimat MF. The computational tool offers accurate results with an error in the range of 1% for small volume fraction (≤10%), and as the volume fraction goes up, the error for the effective Young’s modulus reaches 36%, and for the Poisson’s ratio, 7%.

It is also apparent from the figures that increasing the volume fraction of the inclusions increases the effective Young’s modulus and decreases the Poisson’s ratio. The best properties were obtained from the GNP-UP mixture, which reached 260 GPa for the elastic modulus, and the lowest properties were obtained from the CNT-PP. The combination of any matrix with GNP offered better elastic properties because of their excellent mechanical properties, contrary to previous works where we usually find that the CNT inclusions offers better enhancement. These results can be explained by the non-validity of the isotropy assumptions of CNTs that can lead to an overestimation of the effective properties of the reinforced matrix [8].

3. Second Homogenization

3.1. Mathematical Model

3.1.1. Isotropic Fiber

The elastic properties of alfa and glass fibers are considered as isotropic, as described in

Table 5. Mechanical properties of isotropic fibers.

Materials	Alfa	E-Glass
Model	Elastic	Elastic
Symmetry	Isotropic	Isotropic
Density ρ (g/cm3)	1.52	2.54
Young’s modulus E (GPa)	19.4	73
Poisson’s ratio v	0.34	0.23

.

a. 

Chamis approach

The method of Chamis [28] represents the most used homogenization technique, with five independent constants used to describe the mechanical behavior of the composite. The equations defining the elastic properties of a unidirectional lamina made of anisotropic fibers in an isotropic matrix are formulated by Equations (25)–(30). Alfa and e-glass fibers are considered as isotropic; the axial and in-plane Young’s and shear modulus are the same $E_{11}^f=E_{22}^f$) and ($G_{12}^f=G_{23}^f$).

The relations shown below define the elastic properties of a unidirectional lamina made of anisotropic fibers in an isotropic matrix.

$$E_{11}=E_1^f{V}_f+E_m^{*}{V}_m^{*}$$

(25)

$$E_{22}=E_{33}=\frac{E_m^{*}}{1-V_f \left(1-E_m^{*}/E_2^f\right)}$$

(26)

$$G_{23}=\frac{G_m^{*}}{1-V_f \left(1-G_m^{*}/G_{23}^f\right)}$$

(27)

$$G_{12}=G_{13}=\frac{G_m^{*}}{1-V_f \left(1-G_m^{*}/G_{12}^f\right)}$$

(28)

$${\nu }_{23}=\frac{E_{22}}{2G_{23}}-1$$

(29)

$${\nu }_{12}={\nu }_{13}={\nu }_{12}^f{V}_f+{\nu }_m^{*}{V}_m^{*}$$

(30)

It should be noted that the properties along the third direction coincide with the ones in the second direction due to the assumption of a transversely isotropic medium. Therefore, five independent constants are required to describe the mechanical behavior of the composite. Due to its simplicity, it is one of the most exploited homogenization techniques. It is easy to find many recent works in the literature where the term “rule of mixture” is introduced to define the current approach.

b. 

Hashin–Rosen model

The homogenization technique based on a variational method illustrated in the work by Hashin and Rosen [29] represents a simple approach to compute the mechanical properties of a unidirectional fiber-reinforced composite in which both the fibers and the matrix have isotropic features. As explained in the previous sections, this assumption is valid if the fibers are made of e-glass. On the contrary, this approach cannot be employed for carbon or Kevlar fibers since these kinds of reinforcing phases are transversely isotropic. The engineering constants for the current method assume the following aspect.

$$E_{11}=E_1^f{V}_f+E_m^{*}{V}_m^{*}+\frac{4V_f{V}_m^{*} {\left({\nu }_f-{\nu }_m^{*}\right)}^2 }{V_m^{*}/{\overline{k}}_f+V_f/{\overline{k}}_m^{*}+1/G_m^{*}}$$

(31)

$$E_{22}=E_{33}=\frac{4 {\overline{k}}_t{G}_t}{{\overline{k}}_t+G_t \left(1+4{\overline{k}}_t{\left({\nu }_{12}\right)}^2/E_{11}\right)}$$

(32)

$${\nu }_{23}=\frac{E_{22}}{2G_t}-1$$

(33)

$${\nu }_{12}={\nu }_{13}={\nu }_f{V}_f+{\nu }_m^{*}{V}_m^{*}+\frac{V_f{V}_m^{*}\left({\nu }_f-{\nu }_m^{*}\right)\left(1/{\overline{k}}_m^{*}+1/{\overline{k}}_f\right)}{V_m^{*}/{\overline{k}}_f+V_f/{\overline{k}}_m^{*}+1/G_m^{*}}$$

(34)

$$G_{23}=\frac{E_{22}}{2\left(1+{\nu }_{23}\right)}$$

(35)

$$G_{12}=G_{13}=G_m^{*}\frac{V_m^{*}{G}_m^{*}+\left(1+V_f\right)G_f}{\left(1-V_f\right)G_m^{*}+V_m^{*}{G}_f}$$

(36)

where the following quantities are required

$${\overline{k}}_f=\frac{E_f}{2(1-{\nu }_f-{\left({\nu }_f\right)}^2}$$

(37)

$${\overline{k}}_m^{*}=\frac{E_m^{*}}{2(1-{\nu }_m^{*}-{\left( {\nu }_m^{*}\right)}^2}$$

(38)

$${\overline{k}}_t=\frac{{\overline{k}}_m^{*}{\overline{k}}_f+\left(V_f{\overline{k}}_f+V_m^{*}{\overline{k}}_m^{*}\right)G_m^{*}}{V_m^{*}{\overline{k}}_f+V_f{\overline{k}}_m^{*}+G_m^{*}}$$

(39)

$$G_t=G_m^{*}\frac{\left(\alpha +{\beta }_m^{*}{V}_f\right) \left(1+\xi {\left(V_f\right)}^3\right)-3V_f{\left(V_m^{*}{\beta }_m^{*}\right)}^2}{\left(\alpha -V_f\right)\left(1+\xi {\left(V_f\right)}^3\right)-3V_f{\left(V_m^{*}{\beta }_m^{*}\right)}^2}$$

(40)

$$\alpha =\frac{G_t/G_m^{*}+{\beta }_m^{*}}{G_t/G_m^{*}-1}$$

(41)

$${\beta }_m^{*}=\frac{1}{3-4{\nu }_m^{*}}$$

(42)

$${\beta }_f=\frac{1}{3-{\nu }_f}$$

(43)

$$\xi =\frac{{\beta }_m^{*}-{\beta }_f{G}_f/G_m^{*}}{1+{\beta }_f{G}_f/G_m^{*}}$$

(44)

For the sake of completeness, it should be noted that ${\overline{k}}_f$ and ${\overline{k}}_m^{*}$ denote the plane strain bulk moduli for the fibers and the hybrid matrix, respectively.

c. 

Halpin–Tsai model

The Halpin–Tsai model equations have been used for a long time to predict the properties of composites reinforced with short fibers. These equations were originally developed for long-fiber composites [30]. For composites with aligned short fibers, the Halpin–Tsai model describes the modulus of elasticity as

$$\frac{P}{P_m}=\frac{1+\zeta \eta V_f}{1-\eta V_f}$$

(45)

$$\eta =\frac{P_f/P_m-1}{P_f/P_m+\zeta }$$

(46)

where P represents the property of the composite while $P_f$ and $P_m$ are the corresponding fibers and matrix properties, which can be Young’s modulus E, Shear modulus G or Poisson’s ratio ν. Where $V_f$ denotes the fiber volume fraction and ζ is geometry parameter that depends on the geometry of the filler as follows:

• ζ = 2 L/d for calculation of the longitudinal modulus.
• ζ = 2 for calculation of the transversal modulus.

3.1.2. Transversely Isotropic Fiber

a. 

Hahn model

The homogenization technique proposed by Hahn in his work is used to obtain the mechanical properties of a fibrous composite made of fibers with a circular cross-section randomly distributed in a plane normal to the direction of the oriented fibers. As a consequence, the composite turns out to be macroscopically transversely isotropic and five independent elastic constants are needed to characterize it. The following quantities below are required

$${\Delta }_1=\frac{1+G_m^{*}/G_{12}^f}{2}$$

(47)

$${\Delta }_2=\frac{3-4{\nu }_m^{*}+G_m^{*}/G_{23}^f}{4\left(1-{\nu }_m^{*}\right)}$$

(48)

$${\Delta }_K=\frac{1+G_m^{*}/K_f}{2\left(1-{\nu }_m^{*}\right)}$$

(49)

where $K_f$ is the bulk modulus of the reinforcing fibers defined as follows

$$K_f=\frac{E_2^f}{3 \left(1-2{\nu }_{23}^f\right)}$$

(50)

The plane strain bulk modulus of the composite$K_T$ can be now evaluated

$$K_T=\frac{V_f+{\mathsf{\Delta }}_K{V}_m^{*}}{V_f/K_f+{\mathsf{\Delta }}_K{V}_m^{*}/K_m^{*}}$$

(51)

where $K_m^{*}$ is the bulk modulus of the hybrid matrix that assumes the following aspect

$$K_m^{*}=\frac{E_m^{*}}{3 \left(1-2{\nu }_m^{*}\right)}$$

(52)

First, the following elastic parameters are obtained

$$E_1=E_1^f{V}_f+E_m^{*}{V}_m^{*}$$

(53)

$${\nu }_{12}={\nu }_{13}={\nu }_{12}^f{V}_f+{\nu }_m^{*}{V}_m^{*}$$

(54)

$$G_{12}=G_{13}=\frac{V_f+{\mathsf{\Delta }}_1{V}_m^{*}}{V_f/G_{12}^f+{\mathsf{\Delta }}_1{V}_m^{*}/G_m^{*}}$$

(55)

$$G_{23}=\frac{V_f+{\mathsf{\Delta }}_2{V}_m^{*}}{V_f/G_{23}^f+{\mathsf{\Delta }}_2{V}_m^{*}/G_m^{*}}$$

(56)

Once the five independent constants $K_T$, $E_1$, ${\nu }_{12}={\nu }_{13}$, $G_{12}=G_{13}$, $G_{23}$ are evaluated, the remaining quantities can be computed as follows

$$E_2=E_3=\frac{4E_1{K}_T{G}_{23}}{E_1{K}_T+G_{23} \left(E_1+4K_T{\left({\nu }_{12}\right)}^2\right)}$$

(57)

$${\nu }_{23}={\nu }_{12}^f{V}_f+{\nu }_m^{*}{V}_m^{*}\left(\frac{1+{\nu }_m^{*}-{\nu }_{12}{E}_m^{*}/E_1}{1+{\left({\nu }_m^{*}\right)}^2+{\nu }_m^{*}{\nu }_{12}{E}_m^{*}/E_1}\right)$$

(58)

It should be noted that relations (57) and (58) are equal to the corresponding formulas of the previous approach. The methodology proposed by Hahn belongs to the mechanics of material methods as well.

b. 

Halpin-Tsai model

The homogenization method, which takes its name from the works by Halpin and by Tsai [30], evaluates the mechanical properties of the three-phase multiscale composite by using the Hill’s elastic moduli and a semi-empirical approach. As far as the reinforcing fibers are concerned, the Hills’s elastic moduli takes the following form:

$$k_f=\frac{E_2^f}{2 \left(1-{\nu }_{23}^f-2{\nu }_{21}^f{\nu }_{12}^f\right)}$$

(59)

$$l_f=\frac{{\nu }_{12}^f{E}_2^f}{2 \left(1-{\nu }_{23}^f-2{\nu }_{21}^f{\nu }_{12}^f\right)}=2{\nu }_{12}^f{k}_f$$

(60)

$$m_f=\frac{E_2^f}{2 \left(1+{\nu }_{23}^f\right)}=\frac{1-{\nu }_{23}^f-2{\nu }_{21}^f{\nu }_{12}^f}{1+{\nu }_{23}^f}{k}_f$$

(61)

$$n_f=\frac{E_1^f\left(1-{\nu }_{23}^f\right)}{\left(1-{\nu }_{23}^f-2{\nu }_{21}^f{\nu }_{12}^f\right)}=2 \left(1-{\nu }_{23}^f\right)\frac{E_1^f}{E_2^f}{k}_f$$

(62)

$$p_r=G_{12}^f$$

(63)

where $k_f$, $l_f$, $m_f$, $n_f$, $p_f$ denote the properties of the fibers in terms of Hill’s moduli. On the other hand, if the isotropic hybrid matrix is considered, one obtains

In which $k_m^{*}$, $l_m^{*}$, $m_m^{*}$, $n_m^{*}$, $p_m^{*}$are the Hill’s moduli of the hybrid matrix. The overall mechanical properties of the composite in terms of the Hill’s moduli can be now evaluated

$$k=\frac{k_m^{*}\left(k_f+m_m^{*}\right)V_m^{*}+k_f\left(k_m^{*}+m_m^{*}\right)V_f}{\left(k_f+m_m^{*}\right)V_m^{*}+\left(k_m^{*}+m_m^{*}\right)V_f}$$

(64)

$$l=V_f{l}_f+V_m^{*}{l}_m^{*}+\frac{l_f-l_m^{*}}{k_f-k_m^{*}} \left(k-V_f{k}_f-V_m^{*}{k}_m^{*}\right)$$

(65)

$$m=m_m^{*}\frac{2V_f{m}_f\left(k_m^{*}+m_m^{*}\right)+2V_m^{*}{m}_f{m}_m^{*}+V_m^{*}{k}_m^{*}\left(m_f+m_m^{*}\right)}{2V_f{m}_m^{*}\left(k_m^{*}+m_m^{*}\right)+2V_m^{*}{m}_f{m}_m^{*}+V_m^{*}{k}_m^{*}\left(m_f+m_m^{*}\right)}$$

(66)

$$n=V_f{n}_f+V_m^{*}{n}_m^{*}+{\left(\frac{l_f-l_m^{*}}{k_f-k_m^{*}}\right)}^2 \left(k-V_f{k}_f-V_m^{*}{k}_m^{*}\right)$$

(67)

$$p=\frac{\left(p_f+p_m^{*}\right)p_m^{*}{V}_m^{*}+2p_f{p}_m^{*}{V}_f}{\left(p_f+p_m^{*}\right)V_m^{*}+2p_m^{*}{V}_f}$$

(68)

Finally, the engineering constants of the composite are found

$$E_1=n-\frac{l^2}{k}$$

(69)

$$E_2=E_3=\frac{4m \left(kn-l^2\right)}{kn-l^2+mn}$$

(70)

$${\nu }_{12}={\nu }_{13}=\frac{l}{2k}$$

(71)

$${\nu }_{23}=\frac{n \left(k-m\right)-l^2}{n \left(k+m\right)-l^2}$$

(72)

$$G_{12}=G_{13}=p$$

(73)

$$G_{23}=m$$

(74)

It should be noticed that expressions (69)–(74) coincide with the definitions shown in Equations (18)–(23) for the CNT-equivalent continuum model.

3.2. Numerical Model

The second homogenization of reinforced matrix and fiber was performed by using Digimat FE (see Figure 6). The RVE size was automatically defined, and by applying a uniaxial loading, the transversely isotropic properties were computed. The diameter of the alfa fiber is 95 µm. At this point, this reinforced matrix can be combined with the reinforcing fibers to obtain the desired three-phase multiscale composite.

3.3. Second Homogenization Results

To evaluate the validity of the approaches used in this study, we tried to compare results from the analytical methods (Chamis, Hashin–Rosen, Halpin–Tsai for isotropic alfa and e-glass fibers and Hahn and Halpin–Tsai for transversely isotropic carbon fibers) and the numerical methods using Digimat MF and FE. In Figure 7 and Figure 8, a comparison between the approaches was illustrated. The use of different micromechanics approaches provides in some cases different results, so it is important to discuss these variations. It can be observed for all the elastic parameters that all the approaches show an excellent agreement for small volume fractions. However, the more the amount of the nanofillers is increased, the more differences are observed. In particular, the axial Young’s modulus is underestimated when the Chamis approach is employed. Contrarywise, when evaluating the transverse Young’s modulus, it is observed that this property is underestimated by the numerical methods and the various analytical approaches give similar and higher values, thus two curve groups can be noticed, the first group curve is composed of analytical approaches and is above the second group by 30%.

The longitudinal Young’s modulus increases as a function of CNT and GNP volume fraction-reinforced alfa and glass fibers with PP and UP polymers. The glass-reinforced matrix has a longitudinal Young’s modulus greater than the alfa-reinforced matrix, especially for UP polymer (Figure 9a). Alfa fiber-reinforced polymer has maximum values of longitudinal Young’s modulus in the order of 70 GPa, 68 GPa, 64 GPa and 62 GPa for CNT-UP, GNP-UP, CNT-PP and GNP-PP, respectively. Concerning the glass fiber-reinforced polymer, the longitudinal Young’s modulus has maximum values in the order of 140 GPa, 138 GPa, 112 GPa and 110 GPa for CNT-UP, GNP-UP, CNT-PP and GNP-PP, respectively. Enhanced elastic properties were observed for the CNT-UP effective matrix and the worst were observed for GNP-PP. In addition, the transversal Young’s modulus was lower than the longitudinal Young’s modulus (Figure 9b). The alfa-reinforced composite has a higher transverse Poisson’s ratio ν12 than the glass-reinforced composite (Figure 9c). The in-plane Poisson’s ratio ν23 of the glass-reinforced composite decreased for all volume fractions, while that of the alfa-reinforced composite decreased until 30% and then remained constant for the UP matrix (Figure 9d). It was concluded that CNT inclusion and glass fibers have a greater effect than GNP inclusion and alfa fibers on the mechanical performance of the biocomposite. When the volume fraction of CNT and GNP inclusions is large, their mechanical properties dominate over those of the alfa and glass fibers. Due to their high elastic properties, UP polymer, CNT inclusion and glass fiber can improve the composite properties more than PP polymer, GNP inclusion and alfa fibers. GNPs have a higher Young’s modulus than CNTs, but because of their plate shape, the composite does not have increased mechanical stiffness. Focusing on the alfa fiber, for the polymer, E11 of alfa CNT-UP has a higher value than alfa CNT-PP by 23%. For inclusion, alfa CNT-UP has an E11 value greater than alfa GNP-UP by 17.01%. In the same way, the longitudinal Young’s modulus of glass CNT-UP is 45% greater than alfa CNT-UP. The same observations were found for the other elastic properties of transversal Young’s modulus, out-of-plane Poisson’s ratio, and in-plane Poisson’s ratio. It was observed that in the three phases of the biocomposite, CNT and GNP inclusions dominated the mechanical properties at high volume fractions.

4. Effect of Aspect Ratio on the Mechanical Properties

To investigate the geometry effect of the inclusions, we studied the effect of the aspect ratio on the elastic properties of the reinforced matrix. These results were calculated using the Digimat MF tool, considering the aspect ratio of CNTs to vary between 50 and 1000, and of GNP to vary between 1 and 10⁻⁴.

For CNT inclusions (Figure 10a), increasing the aspect ratio means increasing the effective stiffness of the reinforced matrix and decreasing Poisson’s ratio. Increasing the aspect ratio by 100 from 100 to 200 doubles the effective Young’s modulus. However, the more we increase the aspect ratio, the slower the increase rate, until reaching a Young’s modulus of 180 GPa, where increasing the aspect ratio has no effect. For GNP inclusions (Figure 10b), the effective Young’s modulus increases with increasing the volume fraction and the aspect ratio; however, for the Poisson’s ratio, increasing the aspect ratio decreases ${\nu }_m^{*}$, and we can also notice that for low aspect ratios, Poisson’s ratio varies linearly.

5. Effect of Agglomeration on the Mechanical Properties

5.1. Agglomeration of CNTs

First, the nanosized reinforcing phase is apt to agglomerate, causing a spatially non-uniform distribution of the fillers in composites. To theoretically examine the agglomeration effect, we assumed that some graphene sheets are concentrated in some spherical regions in the matrix, while the rest keep a desired uniform dispersion, as illustrated in Figure 11. Then, the composite was divided into two parts with different loadings of reinforcement which can be considered as two different phases in calculation. We gave the names agglomeration phase and effective matrix phase to the regions inside and outside the spheres, respectively, and further recommend that in both phases the graphene sheets are randomly oriented. In our model, a two-parameter description was employed to characterize such a special distribution [23], where the agglomeration parameters are written as

$$\mu =\frac{W^{in}}{W} ; \eta =\frac{W_r^{in}}{W_r}$$

(75)

If $W_r$ denotes the total volume of CNTs inserted in the matrix, it can be defined as the sum of two contributes

$$W_r=W_r^{in}+W_r^m$$

(76)

In which $W_r^{in}$ and $W_r^m$ represent the volume of nanofillers within the inclusions and scattered in the matrix, respectively. Thus, the total volume W of a representative element is given by the following relation

$$W=W_r+W_m$$

(77)

where $W_m$ stands for the matrix volume.

The volume fraction of CNTs $V_r$ and of the matrix $V_m$ can be expressed as follows

$$V_r=\frac{W_r}{W} \mathrm{and} V_m=\frac{W_m}{W}$$

(78)

In which the relation $V_r$+$V_m$ = 1 allows to relate these two quantities. The reinforcing phase $V_r$ can be seen as the sum of the volume fraction of CNT within the inclusions $V_r^{in}$ and the volume fraction of the nanoparticles scattered in the matrix $V_r^m$ as shown below

$$V_r=V_r^{in}+V_r^m$$

(79)

Then the volume fraction of nanofillers in the agglomeration phase $V_r^{in}$ and that in the effective matrix phase $V_r^m$ can be formulated as

$$V_r^{in}=\frac{W_r^{in}}{W^{in}}=\frac{\eta V_r}{\mu }$$

(80)

$$V_r^m=\frac{W_r-W_r^{in}}{W-W^{in}}=\frac{\left(1-\eta \right)V_r}{\left(1-\mu \right)}$$

(81)

A two-step procedure was applied to reckon the above model. First, the overall elastic properties of the agglomeration and the effective matrix phases are calculated from Equations (80) and (85) by replacing $V_r$ with $V_r^{in}$ and $V_r^{out}$given in Equation (84), respectively. The bulk modulus of the spherical inclusions $K_{in}^{*}$, and the effective matrix $K_{out}^{*}$ are computed as

$$K_{in}^{*}=K_m+\frac{\eta V_r \left({\delta }_r-3K_m{\alpha }_r\right)}{3\left[\mu +\eta V_r\left(1-{\alpha }_r\right)\right]}$$

(82)

$$K_{out}^{*}=K_m+\frac{\left(1-\eta \right)V_r \left({\delta }_r-3K_m{\alpha }_r\right)}{3\left[1-\mu +V_r\left(1-\eta \right)\left({\alpha }_r-1\right)\right]}$$

(83)

Whereas the corresponding shear modulus of the spherical inclusions $G_{in}^{*}$, and the effective matrix $G_{out}^{*}$ are computed as

$$G_{in}^{*}=G_m+\frac{\eta V_r \left({\eta }_r-2G_m{\beta }_r\right)}{2\left[\mu +\eta V_r\left( {\beta }_r-1\right)\right]}$$

(84)

$$G_{out}^{*}=G_m+\frac{\left(1-\eta \right)V_r \left({\eta }_r-2G_m{\beta }_r\right)}{2\left[1-\mu +V_r\left(1-\eta \right)\left({\beta }_r-1\right)\right]}$$

(85)

It should be pointed out that $K_m$ and $G_m$ specify the bulk and the shear moduli of the sole isotropic matrix, can be formulated from Equation (90). For the elastic Hill’s moduli, CNT inclusions are computed from Equations (3)–(6).

$$G_m=\frac{E_m}{2 \left(1+{\nu }_m\right)} K_m=\frac{E_m}{3 \left(1-2{\nu }_m\right)}$$

(86)

Secondly, the agglomeration phase is considered as spherical inclusions embedded in the effective matrix. The effective bulk moduli $K_m^{*}$ and shear moduli $G_m^{*}$ of this hybrid matrix, enriched by CNTs or GNP both included in the spherical inclusions and scattered in the matrix can be computed as

$$K_m^{*}=K_{out}^{*}\left(1+\frac{\mu \left(\frac{K_{in}^{*}}{K_{out}^{*}}-1\right)}{1+\left(1-\mu \right)\left(\frac{K_{in}^{*}}{K_{out}^{*}}-1\right)\left(\frac{1+{\nu }_{out}^{*}}{3-3{\nu }_{out}^{*}}\right)}\right)$$

(87)

$$G_m^{*}=G_{out}^{*}\left(1+\frac{\mu \left(\frac{G_{in}^{*}}{G_{out}^{*}}-1\right)}{1+\left(1-\mu \right)\left(\frac{G_{in}^{*}}{G_{out}^{*}}-1\right)\left(\frac{8-10{\nu }_{out}^{*}}{15-15{\nu }_{out}^{*}}\right)}\right)$$

(88)

In which ${\nu }_{out}^{*}$ stands for the Poisson’s ratio defined as follows

$${\nu }_{out}^{*}=\frac{3K_{out}^{*}-2G_{out}^{*}}{6K_{out}^{*}+2G_{out}^{*}}$$

(89)

The hybrid matrix is evidently isotropic. Therefore, the Young’s modulus $E_m^{*}$, the Poisson’s ratio ${\nu }_m^{*}$ and the density ${\rho }_m^{*}$ are required to define its properties. They can be computed from Equations (9)–(10), respectively

For $\eta \ge \mu$, with $\mu , \eta \in \left[0,1\right].$ The first parameter μ defines the width of the spherical inclusion $W^{in }$ with respect to the total volume W. On the other hand, the second agglomeration parameter η quantifies the volume of CNTs inside the inclusions $W_r^{in}$ with respect to total volume of CNTs $W_r$. In general, three limit cases of the agglomeration ratio can be defined depending on the values of $\mu , \eta$.

For $\mu <\eta$ and $\mu <1$ th, the agglomeration is partial, and the nanofibers are both included in the inclusions and scattered in the matrix (Figure 11b). It should be noted that the heterogeneity of CNTs is amplified for $\mu <1$, by increasing the value of the parameter η. In contrast, if $\mu =\eta =1$ is set, the volume of CNTs is completely concentrated in the inclusions $W_r=W_r^{in}$, which coincide with the whole reference domain ($W_{in}=W$). This circumstance is depicted in Figure 11a and represents a null level of agglomeration. Finally, the complete agglomeration of CNTs is defined by $\mu <\eta$ and $\eta =1$ (Figure 11c). In other words, all the CNTs are located within the spherical inclusions.

5.1.1. Method Comparison

Before studying the agglomeration effect using Shi’s model [22], we investigated the validity of the model by comparing results from the analytical model (

Table 6. Mechanical properties used in the Matlab code and Digimat MF for CNT and GNP.

	CNT	GNP
Axial Young’s modulus (GPa)	2079.8	102,000
In plane Young’s modulus (GPa)	421.14	1029.204
In plane Poisson’s ratio	0.59522	0.4
Transverse Poisson’s ratio	0.17164	0.004
Transverse shear modulus (GPa)	791	102,000
In plane shear modulus (GPa)	132	369
Density (g/cm3)	1.2	2.2

) with semi-numerical results from Digimat MF.

It was supposed that the inclusions had a prolate ellipsoid with an aspect ratio of 10⁴, and the agglomeration to have a spherical shape. The volume fraction of CNTs was assumed to be equal to 1%.

As seen in Figure 12, there is an acceptable agreement between the two approaches: until $\mu =0.4$, the results from Digimat MF tend to be linear. The errors between the two approaches can be explained by several explanations: errors can come from the fact that inclusions in the analytical model are assumed to have an infinite aspect ratio, in addition the low value of $\eta$ used in this study plays a role in the validity of this comparison as when we increased the latter parameter, agreement was no longer found.

5.1.2. Results and Discussions

The influences of the parameters $\mu$ and $\eta$ on the effective modulus were investigated individually by fixing another for only a small volume fraction, as in reality the use of nanofillers is limited to low levels. First, we considered the most severe case of agglomeration where all the inclusions were concentrated in the same place which we can express mathematically by $\eta =1$. We plotted the effective Young’s modulus and Poisson’s ratio under different volume fractions in function of $\mu$ (Figure 13).

The effective elastic modulus has a maximum value when the CNTs are uniformly distributed in the composite. As the agglomeration parameter decreases, the effective stiffness also decreases rapidly. It is also noted that for $\mu$ ≤ 0.6, the addition of CNTs has no effect on the elastic properties of the reinforced polymer. Both agglomeration parameters are required to describe this phenomenon. The first η indicates the amount of inclusions located in the agglomeration and the other indicates the size of the agglomerations.

The effective properties are shown in Figure 14 for $\mu$ = 0.2, where it is seen that with the increase of the relative amount of CNTs in the concentrated regions, the effective Young’s modulus increases slightly to a maximum value of 20 GPa for a volume fraction of 4% and $\eta$ = 0.15 before dropping sharply to reach its minimum value when η reaches the unity. For $\mu$ = 0.5 (Figure 15), we can observe that the curves keep the similar tendency but with greater magnitudes, so it can be concluded that increasing the volume of the concentrated regions leads to an increase in the effective properties.

5.2. Agglomeration of GNPs

For the agglomeration of GNPs, a similar mathematical model was used for modelling the different agglomeration phases with a slight difference in the mathematical expression of the bulk and shear modulus where [23]

$$K_m^{*}=K_{out}^{*}\left(1+\frac{\mu \left(\frac{K_{in}^{*}}{K_{out}^{*}}-1\right)}{1+\left(1-\mu \right)\left(\frac{K_{in}^{*}}{K_{out}^{*}}-1\right)\left(\frac{3K_{out}^{*}}{3K_{out}^{*}+4G_{out}^{*}}\right)}\right)$$

(90)

$$G_m^{*}=G_{out}^{*}\left(1+\frac{\mu \left(\frac{G_{in}^{*}}{G_{out}^{*}}-1\right)}{1+\left(1-\mu \right)\left(\frac{G_{in}^{*}}{G_{out}^{*}}-1\right)\left(\frac{6 \left(K_{out}^{*}+2G_{out}^{*}\right)}{5 \left(3K_{out}^{*}+4G_{out}^{*}\right)}\right)}\right)$$

(91)

Using the same technique for CNT agglomeration, we investigated the influence of the two agglomeration parameters by fixing one and varying the other. Figure 13 shows the variation of effective Young’s modulus and Poisson’s ratio in function of µ under $\eta$ = 1 which indicated that all the nanofillers are agglomerated in some spherical areas. The variation of the same properties with $\eta$ under $\mu$ = 0.2 and $\mu$ = 0.5 are presented in Figure 14 and Figure 15.

Results and Discussions

In the extreme case where $\eta$ = 1 (Figure 16), we can observe that the graphene sheets have the same agglomeration behavior (the curves have similar shapes with a difference in the magnitudes) as the carbon nanotubes, where the elastic property of the reinforced matrix drops as the agglomeration parameter decreases. It is concluded from Figure 17 and Figure 18 that increasing $\mu$ indicates that the amount of the fillers located in the agglomeration areas reduce the effective Young’s modulus and Poisson’s ratio. However, increasing the second agglomeration parameter $\eta$ enhances the elastic properties of the matrix, which can be explained by the fact that $\eta$ signifies the size of the agglomeration sub-regions; consequently, the bigger the concentrated areas are, the more uniform the distribution.

6. Conclusions

In this paper, we studied the elastic behavior of a biocomposite based on Alfa fibers and CNT/GNP-reinforced polymer compared with synthetic fibers used in the wind energy field. As the elastic behavior is characterized by the Young’s modulus and Poisson’s ratio, we evaluated those two properties using analytical and numerical approaches, and also investigated the effect of the volume fraction, aspect ratio, and agglomeration of the inclusions.

The key conclusions of this work can be summarized as follows:

Computational and numerical approaches show a good agreement for small volume fraction up to 10%;

Due its shape, the CNT-reinforced composite has better overall properties in the axial direction and the GNP-reinforced composite has better properties in the transversal direction. The effect of the aspect ratio was also studied, and by increasing the aspect ratio from 1 t 0.0001, the elastic modulus increases for the GNP and reaches a maximum of 260 GPa. By increasing the aspect ratio from 50 to 1000, the elastic modulus increases for the CNT and reach a maximum of 180 GPa;

In the second homogenization, the axial and transverse Young’s modulus increases with increases in the volume fraction of CNT and GNP. The glass reinforcement shows the best mechanical behavior with a maximum axial Young’s modulus of 140 GPa; surprisingly, mixing alfa fibers with GNP inclusions gave us better elastic properties than that of the carbon fibers, with a maximum of axial Young’s modulus of 90 GPa.

Finally, as the agglomeration parameter decreases, the effective stiffness also decreases rapidly. It was also noted that for μ ≤ 0.6, the addition of CNTs has no effect on the elastic properties of the reinforced polymer. Results indicated that with the increase of the relative amount of CNTs in the concentrated regions, the effective Young’s modulus increases slightly to a maximum value of 20 GPa for a volume fraction of 4% and η = 0.15 before dropping sharply to reach its minimum value when η reaches the unity.