Computation of Effective Elastic Properties Using a Three-Dimensional Semi-Analytical Approach for Transversely Isotropic Nanocomposites

: A three-dimensional semi-analytical ﬁnite element method (SAFEM-3D) is implemented in this work to calculate the effective properties of periodic elastic-reinforced nanocomposites. Different inclusions are also considered, such as discs, ellipsoidals, spheres, carbon nanotubes (CNT) and carbon nanowires (CNW). The nanocomposites are assumed to have isotropic or transversely isotropic inclusions embedded in an isotropic matrix. The SAFEM-3D approach is developed by combining the two-scale asymptotic homogenization method (AHM) and the ﬁnite element method (FEM). Statements regarding the homogenized local problems on the periodic cell and analytical expressions of the effective elastic coefﬁcients are provided. Homogenized local problems are transformed into boundary problems over one-eighth of the cell. The FEM is implemented based on the principle of the minimum potential energy. The three-dimensional region (periodic cell) is divided into a ﬁnite number of 10-node tetrahedral elements. In addition, the effect of the inclusion’s geometrical shape, volume fraction and length on the effective elastic properties of the composite with aligned or random distributions is studied. Numerical computations are developed and comparisons with other theoretical results are reported. A comparison with experimental values for CNW nanocomposites is also provided, and good agreement is obtained.


Introduction
Single-walled carbon nanotube (SWNT) composites embedded in a polymer matrix are called nanocomposite materials. They show superior mechanical, thermal and processing properties in comparison with other materials that do not have reinforcement. Nanocomposite materials are suitable for aerospace, automotive, aeronautic and sporting goods applications [1][2][3][4]. Current global expectations regarding the creation of high-performance, lightweight and low cost products requires the design of specific properties [5,6]. In this sense, the effective property estimation of multi-phase composites, based on a technological approach supported by modeling and numerical simulation, allows a better understanding and analysis of new materials with higher standards, as well as reducing costs and production times [7][8][9][10]. The effective properties of nanocomposites are influenced by factors such as their structure, orientation, dispersion, fiber diameter and matrix stiffness [11][12][13][14]. Different techniques have been developed to study the elastic behavior of multi-phase fiber-reinforced composites (FRCs) [15,16].
Homogenization techniques have always been a useful tool to describe the structureproperty relationship of composite materials [17][18][19][20][21]. For example, Bravo-Castillero et al. [22] determined the effective properties of elastic FRCs with square and hexagonal fiber distributions by means of the asymptotic homogenization method (AHM) for bone mechanics applications. Rodríguez-Ramos et al. [23] calculated the effective properties of elastic FRCs with transverse isotropic constituents and parallelogram-like periodic unit cells using the AHM. Zhao et al. [24] predicted the effective elastic properties of polymer nanocomposites using a novel numerical implementation of AHM. Baek et al. [25] proposed a two-step homogenization approach to analyze the mechanical behavior of polymeric nanocomposites based on multiscale homogenization and the clustering density concept. In the literature, different homogenization implementations and computer modeling techniques to predict the effective properties of nanocomposites can be found; see for instance [26][27][28]. They represent a rewarding advantage through validations by comparing different models.
In addition, Otero et al. [29] proposed a semi-analytical method to find the effective coefficients of elastic FRCs with aligned inclusions and transversely isotropic constituents. An extension of this 2D semi-analytical method for coupled field problems was reported by Espinosa-Almeyda et al. [30]. Gabbert et al. [31] developed a semi-analytical approach based on FEM linked with the homogenization method to solve the inverse problem to determine microscale distributions. Nasution et al. [32] presented a new asymptotic expansion homogenization model to study the thermodynamic properties of 3D composites. In their work, the mathematical treatment takes into account the benefits of 2D homogenization and periodicity in the in-plane direction. The mechanical behaviors of three-dimensional short fiber composites, including interface debonding, were studied using an extended finite-element method (XFEM) [33]. Fisher et al. [34,35] studied the effect of carbon nanotube (CNT) inclusions on the effective elastic properties of composites by numerical simulation techniques, using finite elements coupled with the Mori-Tanaka model. A computational approach based on FEM was implemented to predict the elastic behavior of a hybrid nanocomposite reinforced by CNTs and nano-clay platelets in an epoxy system [36]. A review of FEM applications in the analysis of composite materials was reported in [37]. More recently, many research works have reported on the field of nanocomposites reinforced by CNT inclusions; for example, Tarfaoui et al. [38] analyzed the effect of CNT additives on the mechanical behavior of three-phase polymer-matrix composites by experimental techniques. Choi et al. [39] proposed a multi-scale method to obtain the mechanical properties of interfacial regions of SWNT/epoxy nanocomposites. E. Moumen et al. [40] estimated the elastic properties of laminated composites reinforced by CNT using a numerical homogenization approach, and comparisons with experimental data were reported. A review of the most important and novel theoretical and experimental methods applied to study the mechanical behavior of 1D, 2D and 3D CNT-reinforced nanocomposites was presented in [41].
In this work, the effective coefficients of periodic multi-phase elastic FRCs are calculated by the three-dimensional semi-analytical finite element method (SAFEM-3D). Herein, the elastic reinforced nanocomposites are assumed to have isotropic or transversely isotropic CNT inclusions within an isotropic matrix. The SAFEM-3D approach is based on a combination of the AHM and the finite element method (FEM). Using the AHM, the homogenized local problems are stated and the analytical expressions of the effective elastic properties are obtained. The FEM allows local problems to be solved via the minimum potential energy principle. Local problems over the unit cell are transformed into boundary problems over one-eighth of the unit cell. Numerical values are reported for different nanocomposite materials with an aligned or random distribution of the inclusions. The effective coefficients of periodic elastic FRCs with different inclusions in the form of discs, ellipsoidals, spheres, carbon nanotubes (CNT) and carbon nanowires (CNW) are calculated. The effect of the inclusion geometry shape, volume fraction and inclusion length on the effective elastic properties is analyzed. Comparisons with some theoretical methods and experimental data reported in the literature show good agreements.

Statement of the Problem for Heterogeneous Media
A heterogeneous medium occupying a volume Ω ∈ 3 with a boundary ∂Ω = ∂ 1 Ω ∪ ∂ 2 Ω is considered as a periodic two-phase elastic FRC (inclusion/matrix) at the Cartesian coordinate system x = (x 1 , x 2 , x 3 ) (see Figure 1). In particular, an elastic FRC with transversely isotropic CNT inclusions embedded into an isotropic matrix is studied. The CNT's inclusions have a periodic distribution along the axes Ox 1 , Ox 2 and Ox 3 without overlapping within the homogeneous matrix. The periodic unit cell S is taken as a cubic matrix with a centered CNT inclusion in the form of discs, spheres, ellipsoids and hollow cylinders at the Cartesian coordinate system y = (y 1 , y 2 , y 3 ). The regions occupied by the matrix S 1 with volume V 1 and the CNT inclusions S 2 with volume V 2 satisfy that S = S 1 ∪ S 2 with S 1 ∩ S 2 = ∅, and the unit cell volume is V = V 1 + V 2 . Perfect contact conditions between the matrix and the inclusions are assumed; that is, mechanical displacements and stresses are continuous across the interface, denoted by Γ (see, for instance, [23]). The values associated with the matrix and the inclusion are indicated by superscripts in parentheses (2) and (1), respectively. The static governing equations for a heterogeneous multi-phase elastic composite associated with the theory of linear elasticity [42,43] are defined by together with the boundary conditions where i, j = 1, 2, 3 and σ ij and u i represent the components of the second-order stress tensor and the displacement vector, respectively. In addition, n j is the normal unit vector to ∂ 2 Ω. Furthermore, the constitutive equations are defined by Hooke's law: where C ijkl (k, l = 1, 2, 3) is the elasticity modulus and ε kl represents the strain components, such as Besides, the perfect contact conditions between the matrix and inclusion phases are assumed: u

Asymptotic Homogenization Method: Homogeneous Problem, Local Problems and Effective Coefficients
The homogenized local problems over S, the homogeneous static problem and effective coefficients linked to the periodic elastic FRC are derived from Equations (1)-(5) using the AHM [17,22,23] by means of the asymptotic expansion of the displacement u i in powers of the parameter , resulting in the following: where u (0) i is a function only of the slow variable x, u (k) i (x, y) (k = 1, 2, · · · ) are S-periodic functions with respect to the fast variable y and the superscripts denote the order of terms in the expansions. Here, the two-scales, denoted by x (global or slow) and y (local or fast), characterize the global behavior of the composite and the heterogeneities at a local level over S, respectively. The scales are related by y = x/ , where is a small enough parameter which represents the ratio between the dimensions of the periodic cell S and the characteristic dimensions of the nanocomposite.
The mathematical statement of the pq L local problems over the unit cell S are stated as subject to the perfect interface conditions, as shown in Equation (5), which in this case are ijpq (y) + pq τ (2) ij (y)n (2) j , on Γ where the components of the Y-periodic local stress tensor are pq τ ij (y) = C ijkl (y) ∂ pq N k (y) ∂y l . (10) and pq N k (k, p, q = 1, 2, 3) are the components of the S-periodic local pseudo-displacements vector, providing the solution of the pq L local problems. Additionally, the following condition, over S, is needed to guarantee a unique solution of the local problems: pq N k = 0. Here, · = (1/|S|) (·)dy is the average operator over the periodic cell S and |S| is the volume of the cell. The homogeneous static problem equivalent to Equation (1) is stated as follows: where the constitutive equations of the homogeneous problem are where u (0) p (x) is the solution of the homogenized problem in Equation (11) and C * ijpq are the effective coefficients, which are defined as follows: The boundary conditions expressed in Equation (2) are transformed into The mathematical procedures to obtain the pq L local problems (Equations (7)-(10)) and the homogeneous static problem (Equations (11)-(12)) can be seen in Appendix A. More details explaining he AHM procedure and the rigorous theoretical mathematical background can be found in the classical works [17,20,44,45].
Six local problems pq L (p, q = 1, 2, 3) in terms of the local variable y are defined in the unit cell S. Each problem is decoupled into an independent set of equations according to the plane and out-of-plane deformation state assumed in linear elasticity. Table 1 shows the correspondence between the effective elastic coefficients and the local problems pq L. Table 1. Local problems and associated effective coefficients.
The local problems can be presented as follows: Problems λ3 L with λ, β = 1, 2 and λ = β : Problems ββ L with β = 1, 2, 3 : Problem 12 L : In the Equations (15)-(17), the elastic coefficients C ijkl (y) = C ijkl (y 1 , y 2 , y 3 ) are assumed to be even functions with respect to the variables y 1 , y 2 and y 3 . Then, the following conditions are satisfied: where d h = {0, l h } (h = 1, 2, 3), l h is the semi-length of the cell in the y h direction and the function δ hi is defined as follows: The local problems pq L over S can be transformed into boundary problems over oneeighth of the periodic unit cell S by applying the conditions in Equations (18) and (19) [46]. Thus, the local functions are related to the new variable pq M k as follows: where p, q, l = 1, 2, 3. From Equations (7), (10) and (11), the following relations can be obtained: and the effective elastic coefficients shown by Equation (13) can be expressed as The local problems pq L can be restated as equivalent boundary value local problems pq L over one-eighth of the cell S as follows: where pq σ pq σ where ββ = 11, 22, 33 and λ, Υ = 1, 2. The boundary conditions associated with the boundary value local problems over one-eighth of S are summarized in Appendix B.
The effective coefficients are calculated with Equation (25) and they are listed according to the relevant problem as follows: Problems ββ L (ββ = 11, 22, 33): Note that the local problems subject to periodic conditions on the cell (S) (Equations (7)-(10)) are transformed into boundary problems ((Equations (26)-(29)) on one-eighth of S, considering the symmetries of the cell through Equation (21). The boundary conditions are given in Appendix B. The FEM implementation via the minimum potential energy principle for the boundary problems is explained in the following section.

Solution of Local Problems through the Finite Element Method
For many cases of FRCs with complex geometrical structures, the determination of analytical solutions for the local problems and the exact expressions of the effective coefficients through the analytic AHM are difficult or even impossible to solve. Only for specific inclusion geometries, the exact analytical solutions have been reported; see for instance [23]. In this sense, an alternative approach is the finite element method. In the present work, the SAFEM-3D approach is developed to find the solution for the local problems and the effective elastic properties. The SAFEM-3D is implemented by combining the AHM and the FEM; i.e., the local problems obtained by the AHM are solved using the FEM through the minimum potential energy principle, in a similar manner to the methodology developed by Otero et al. [29,46] applied to a bi-dimensional elastic twophase FRC.
The total corresponding potential energy Π for an elastic solid medium Ω-see for instance [47]-is defined as where Π involves the deformation energy per unit volume of the body and the potential energies associated with the body forces f , tensile forces T and point loads P i respectively.
In the present work, the forces f and T and loads P i are assumed to be zero. The relationships given by Equations (26)- (29) are related as follows: where pq σ = pq σ 11 pq σ 22 pq σ 33 pq σ 23 pq σ 13 pq σ 21 pq ε = pq ε 11 pq ε 22 pq ε 33 pq ε 23 pq ε 13 pq ε 21 and The region of study (one-eighth of the unit cell S) is divided into a finite number of tetrahedral elements. Each element is characterized by 10 nodes located at the vertices, and the remaining ones are at the midpoints of the edges. Additionally, each node is associated with a pseudo-displacement vector connected with the pq L local problem to be solved. The pseudo-displacements pq M k are defined as a linear combination of the shape functions Ψ i whose coefficients are the nodal values of the unknown pseudo-displacement where pq w k(s) represents the k-th component of the pseudo-displacement at the s-th node (s = 1, · · · , 10) for the local problem pq L.
The shape functions Ψ i are defined through natural coordinates [48], as follows: where ξ 1 , ξ 2 , ξ 3 and ξ 4 = 1 − ξ 1 − ξ 2 − ξ 3 represent the element's natural coordinates. The element node vector pq w is defined as Taking into account Equations (41) and (42), the relationship between deformations and pseudo-displacements as shown in Equation (39) can be written in matrix form as functions of the natural coordinates and the element's shape functions, as follows: where B is the element strain-displacement matrix: ij represent the inverse matrix coefficients of the Jacobian transformation J between the Cartesian system (y 1 , y 2 , y 3 ) and the natural coordinate system (ξ 1 , ξ 2 , ξ 3 ), and the matrix H 9×30 is defined in Appendix C.
The stress values are determined as functions of the natural coordinates and the element's shape functions by means of the relations shown in Equations (37)-(42), such that Replacing Equations (44) and (46) into Equation (36), we obtain the deformation energy Π e associated with an element e in the following form: where K e is the stiffness matrix of the element: Therefore, taking into account the connectivity of the elements, the total energy is equal to where K and pq W represent the global rigid matrix and the global pseudo-displacement vector on S, respectively. The minimum potential energy is found by means of an algebraic system of equations obtained through ∂Π and by applying the corresponding boundary conditions of each pq L local problem; see Equations (A15)-(A20).
The effective properties of the element can be computed as a function of the natural coordinates and the element's shape functions, replacing the derivates of the local functions pq M k into Equations (30)- (35), thus: for each local problem. Therefore, we can summarize: Problems ββ L : The effective elastic coefficient is The steps of the numerical implementation of the SAFEM-3D are shown in the flowchart of Figure A1, Appendix D. Mesh generation is performed using ANSYS (Solid187 element), and the solution to the local problems and effective coefficients is determined by an algorithm implemented in Matlab. The suitable number of elements generated by ANSYS to ensure the accuracy of SAFEM-3D is given by the following steps.
Step 1: The initiation of the generation of a mesh of the periodic cell by ANSYS (with few elements); then, the mesh is imported to a MATLAB program to solve local problems and calculate the effective coefficients.
Step 2: A new mesh with a greater number of elements is generated and the effective coefficients are recalculated.
Step 3: The difference between the values of the effective coefficients from steps 1 and 2 is calculated. Then, the adequate number of elements required is reached when the difference is less than 0.01%; otherwise, repeat step 2 and 3 for as long as necessary.

Results and Discussion
The effective coefficients in Equations (51)-(60) are calculated for two and three-phase periodic FRCs with different orientations of inclusions (aligned or random). The inclusions are defined by their aspect ratios (α), where discs (α < 1), spheres (α = 1) and ellipsoids (α > 1) are assumed. For the numerical calculations, the mechanical properties of the constituents are defined in Table 2 and taken from different works, in which the effective engineering and the Hill's constants in the y 1 -direction are related by where m = C * 2323 = (C * 3333 − C * 2233 )/2 and p = C * 1313 = C * 1212 are the transverse and longitudinal shear moduli, k = (C * 3333 + C * 2233 )/2 is the transverse bulk modulus, l = C * 1122 = C * 1133 is the associated cross modulus and n = C * 1111 is related to uniaxial straining in the y 1direction [49].

Aligned NTC-Reinforced Composites
In this section, the effect properties of two and three-phase FRCs with aligned inclusions (discs, spheres, ellipsoids and hollow and massive cylinders) are studied. The discretization for one-eighth of the periodic cells using ANSYS (Solid187 element) is shown; see Figure 2. In Figure 3, the effective longitudinal Young's modulus (E * 11 ), transverse Young's modulus (E * 22 ), transverse shear modulus (µ * 23 ), Poisson's ratio (ν * 12 ), plain strain bulk modulus (κ * 23 ) and axial shear modulus (µ * 12 ) are reported for a two-phase RFC (graphite/epoxy) [11] as a function of the volume fraction. Furthermore, comparisons between the SAFEM-3D and the Mori-Tanaka method (MT) [11] are shown for three different types of inclusions-discs, spheres and ellipsoids-with aspect ratios of α = 0.5, α = 1 and α = 2, respectively. As can be seen from these comparisons, the effective properties E * 11 , E * 22 , µ * 23 , µ * 12 and κ * 23 are in good agreement for volume fractions less than 0.3; however, a notable difference is shown when the volume fractions are greater than 0.3. On the one hand, there are small differences in the whole range of volume fraction for ν * 12 . In addition, a remarkable difference between both models is observed for increasing values of the aspect ratio. Additionally, an increasingly monotonic behavior of the effective properties can be observed as the reinforcement volume fraction increases, except for ν * 12 .  The discretization of hollow and massive cylinders is shown in Figure 4. Figure 5 illustrates the behavior of the normalized longitudinal Young's modulus (E * 11 /E 11 ) for a two-phase FRC (SWNT/LaRC-SI) as a function of the inclusion's volume fraction. In addi-tion, a comparison between the SAFEM-3D and three different homogenization schemes based on the Mori-Tanaka method (MT) [14] are reported. Here, the SWNT isotropic reinforcement [14] embedded in an isotropic matrix has an aspect ratio of α = 12.3. The SWNT reinforcement is designed with a hollow cylindrical geometry (see Figure 4a) simulating the carbon nanotube. It can be observed that the SAFEM-3D results are close to those reported by the MT models. Furthermore, there are two types of concavity in the effective behaviors according to different approaches.
11 ) as a function of the inclusion's volume fraction for a two-phase fiber-reinforced composite (FRC) (LaRC-SI/SWNT) with hollow and massive nanotubes in an aligned configuration and aspect ratio of α = 12.3. Table 3 shows the effective elastic moduli of a two-phase FRC (Inclusion [46]/matrix [46]) for three different values of the inclusion's volume fraction: Here, the constituent properties of the composite are considered as transversely isotropic materials, and cylindrical inclusions (see Figure 4b) are assumed. In addition, a comparison between the analytical AHM and the semi-analytical approaches of SAFEM-2D [46] and SAFEM-3D are reported. An analysis of the relative error is also shown. As can be seen, good agreement between the AHM and the SAFEM-2D is observed. However, a major difference is observed for the SAFEM-3D approach. This is because the SAFEM-3D calculations are developed considering that the fiber's length in the y 3 axial direction is equal to 0.49 within one-eighth of the cubic cell S with an edge length of 0.5. Furthermore, it is observed that, as the inclusion volume fraction increases, the relative errors increase, and this is more noticeable for C * 3333 . Here, the SAFEM-3D model is an approximation to the long-fiber case. On the other hand, as expected, an increase of the inclusion volume fraction leads to an increment in the effective properties, because the elastic properties of the inclusions are higher than those of the matrix. Table 3. Effective elastic coefficients (GPa) and relative errors:  Table 4 reports the effective elastic moduli of a two-phase FRC computed by a semianalytical method (SAM) [50], the self-consistent method (SCM) [50] and the SAFEM-3D as a function of the aspect ratios (α) and the volume fractions (v f ). In addition, an analysis of the relative error between the SAFEM-3D and the other two models (SAM and SCM) is also shown. Here, a two-phase FRC with transversely isotropic ellipsoidal inclusions (Al 6061) embedded in a matrix [50] is assumed; see Table 2. A good agreement among the three models is observed, with a relative error below 1% in every case. This result is to be expected because small inclusion volume fractions are assumed. On the other hand, it is observed that all the effective elastic properties increase as the volume fraction of the ellipsoidal inclusions decreases. This fact is caused by a much stiffer matrix and more matrix material within the system. Table 4. Effective elastic coefficients (GPa) for a two-phase FRC reinforced by ellipsoidal inclusions with four different aspect ratios (α) and volume fractions (v f ). Here, the relative errors:  Table 5 illustrates the effective elastic coefficients of a three-phase FRC for different volume fractions of the inclusions and computed relative errors (Error 1 and Error 2), as in Table 4. Here, the composite material is reinforced by ellipsoidal (α = 0.666) and spherical (α = 1) inclusions. For the computations, the constituents materials of the SiC matrix, spherical inclusions (Al 6061) and ellipsoidal inclusions (matrix [50]) are given in Table 2. In this analysis, the configuration of these inclusions within the composite is performed to validate the SAFEM-3D model within the range aspect ratios and volume fractions used in this work. Additionally, it is worth mentioning that the distribution and shape of the fibers provide a broader view on the scope of the SAFEM-3D model. As can be seen in Table 5, a good agreement is observed between the models. Moreover, the effective elastic properties decrease as the volume fraction of both inclusions increase. This fact is due to a stiffer SiC matrix being considered; therefore, there is less matrix material in the system. Here, the relative errors are in general below 2.2%. Table 5. Effective elastic coefficients (GPa) for a three-phase FRC reinforced by ellipsoidal and spherical inclusions for different aspect ratios (α) and volume fractions (v f ).  Figure 6 illustrates the longitudinal Young's modulus (E * 11 ) as a function of the nanotube length, for a two-phase composite reinforced by aligned isotropic SWNT nanotubes [12] in a LaRC-SI polymeric matrix. Here, the SWNT nanotubes are simulated by ellipsoidal inclusions within a matrix. In addition, a comparison between the SAFEM-3D and the continuous-equivalent method [12] is reported. The numerical calculations are performed considering a SWNT nanotube volume fraction equal to 0.01. As can be seen, very good agreement is observed. Figure 6. Effective longitudinal Young's modulus for a two-phase FRC with an aligned isotropic SWNT nanotube [12] with a 1% volume fraction.

Random NTC Reinforced Composites
In this section, the effect of the distribution of inclusions on the effective longitudinal Young's modulus (E * 11 ) property is studied. For this purpose, the SAFEM-3D is implemented over different two-phase FRCs, in which the periodic unit cell has 72 similar inclusions embedded in the matrix. The effect of non-aligned inclusions is determined by their rotations with respect to the axial axis y 1 . Different fiber orientations (θ = 0 • (aligned inclusions), 10 • , 15 • , 20 • , 25 • , 30 • , and 35 • (random inclusions)) are considered. Non-alignment is implemented by dividing one-eighth of the cell into nine prisms containing a single inclusion, which are assumed to have different orientations. No contact between inclusions is allowed. Then, the inclusions are randomly oriented around the axial axis, but all of them have a fixed angle with respect to that axis. In this way, the effect of non-aligned inclusions can be described. For transversely isotropic inclusions, the elastic tensor should be rotated according to the inclusion orientation. In Figure 7, the discretization for one-eighth of the periodic cells using ANSYS (solid187 element) is shown for different inclusion shapes (discs, spheres and ellipsoids) with 1% of the inclusion volume fraction.
In Figure 8, the effective longitudinal Young's modulus (E * 11 ) as a function of the inclusion aspect ratio α is shown for a two-phase FRC with aligned (θ = 0 • ) and random (θ = 35 • ) transversely isotropic SWNT [13] inclusions in an isotropic LaRC-SI matrix. Here, the numerical calculation is determined considering a 1% SWNT volume fraction for different aspect ratio values 0 ≤ α ≤ 14. The inclusion shapes for α = 0.6, 1 and 10 can be seen in Figure 7. Additionally, a comparison between the SAFEM-3D and the deformation energy method based on the Mori-Tanaka model [12,13] is reported. A good agreement is observed between both models. From Figure 8, it is observed that the aligned nanotubes array shows a maximum value of the longitudinal Young's modulus in the whole range of aspect ratio values.  In Figure 9, the effective longitudinal Young's modulus (E * 11 ) as a function of the inclusion's aspect ratio α for different fiber orientations (θ = 0 • (aligned), 10 • , 15 • , 20 • , 25 • , 30 • and 35 • (random)) is shown. Here, the periodic FRC is defined by transversely isotropic SWNT [13] inclusions (see Figure 7) in an isotropic LaRC-SI matrix, where a 1% SWNT volume fraction is considered. Besides, a comparison between the SAFEM-3D and the deformation energy method based on the Mori-Tanaka model [12,13] is given for aligned and random fiber distributions. The comparisons show good agreement. From Figure 9, it is observed that the effective longitudinal Young's modulus (E * 11 ) has a strong dependence on the fiber orientation. As θ increases, the value of the longitudinal Young's modulus decreases. Maximum values of E * 11 are obtained for the aligned distribution of the inclusions. The effect of carbon nanowire (CNW) inclusions on the effective longitudinal Young's modulus (E * 11 ) is now presented. The CNWs are simulated by ellipsoidal inclusions with different orientations, as seen in Figure 10. In Figure 11, the behavior of E * 11 as a function of the CNW weight content (CNW wt content)(%) for different CNW orientations θ = 0 • (Aligned), 10 • , 20 • , 30 • and 35 • (Random) is reported. In addition, a comparison between the SAFEM-3D results and experimental data [10] is shown. For the numerical computation, the material constituents of the matrix (polystyrene) and the CNW inclusions are assumed to be isotropic materials, such as E = 3.349 GPa and ν = 0.34 for polystyrene [10] and E = 169 GPa and ν = 0.185 for CNW [51]. As can be observed, the theoretical results using the SAFEM-3D show that E * 11 decreases as θ increases; i.e., E * 11 for aligned CNW inclusions (θ = 0) is greater than that for random CNW inclusions (θ = 35 • ) in the y 1 axial direction. It is remarkable that the SAFEM-3D values for random inclusions are within the expected error range when the CNW inclusions are between 1% and 5% CNW wt content. From a physical point of view, as the CNT wt content increases, E * 11 should grow, as can be seen in the SAFEM-3D results ( Figure 11). The difference between the SAFEM-3D values and the experimental values when the CNW inclusions have 5% CNW wt content is also reported in [10].

Conclusions
In this work, the estimation of the effective elastic properties of periodic reinforced nanocomposites is developed. The nanocomposites are defined with isotropic or transversely isotropic inclusions embedded in an isotropic matrix. Different inclusions in the form of discs, ellipsoidals, spheres, carbon nanotubes and carbon nanowires are also considered. Two configurations of inclusions are taken into account: aligned and random. A three-dimensional semi-analytical approach (SAFEM-3D) is implemented, using a combination of the asymptotic homogenization method and the finite element method. The asymptotic homogenization method provides the local problems and the effective coefficients, whereas the finite element method gives the solution to the local problem by using the minimum potential energy.
The results have been compared with theoretical models validating the efficiency and accuracy of the SAFEM-3D for the analysis of the composites. Based on the numerical results, we have concluded that fiber-reinforced composites with aligned nanotube inclusions show a higher effective longitudinal Young's modulus than random nanotube arrays. Random carbon nanotube arrays show that the nanotube length has a significant effect on the composite elastic properties in the axial direction. Furthermore, according to the values of the computed relative errors , it is concluded that the SAFEM-3D values have a relative error less than 2.2% when compared with the other reported theoretical models; see Tables 4 and 5. These results show the efficiency and accuracy of the SAFEM-3D model. For the case of long fibers, a slight difference is observed between SAFEM-3D and AHM (or SAFEM-2D) values when the inclusion volume fraction increases; see Table 3. This error increment is not attributed to the accuracy of SAFEM-3D; it is rather a result of the fact that the long fibers are not considered in the model.
A comparison with experimental data for CNW nanocomposites is reported and a good agreement is obtained. An increase in the CNW weight content enhances the magnitude of E * 11 . This effect is described by the SAFEM-3D results for different orientations of inclusions. This is also observed in the experiment when the CNW weight content is 1% and 3%; however, when the CNW weight content is equal to 5%, a decrease in the property is observed. From the physical point of view, we consider that this last result is not consistent with the expected behavior; thus, we assume that it is a consequence of the composite's manufacture.
An important advantage of the SAFEM-3D approach is the possibility to design any kind of geometry reinforcement for any type of material, including isotropic and transversely isotropic nanocomposites in three-dimensional configurations for a huge range of random angle values.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.
Substituting Equation (A11) into Equation (A10), the pq L local problems statement in Equation (7) is obtained. Next, replacing Equation (A11) into Equations (A4) and (A6): Applying the average operator and the periodicity property of pq N(y) in (A13), the homogeneous deformation tensor can be written as From the average of Equation (A5) and taking into account the periodicity of σ (1) ij , the homogeneous static problem in Equation (11) is obtained.

Appendix C. Matrix H
The matrix H 9×30 defined in Equation (45) related to the strain-displacement matrix B is presented through the form