- Email: [email protected]
FEM analysis of the elastic behavior of composites and nanocomposites with arbitrarily oriented reinforcements
Journal Pre-proofs FEM analysis of the elastic behavior of composites and nanocomposites with arbitrarily oriented reinforcements A. Greco PII: DOI: Reference:
S0263-8223(19)34105-4 https://doi.org/10.1016/j.compstruct.2020.112095 COST 112095
To appear in:
Composite Structures
Received Date: Revised Date: Accepted Date:
28 October 2019 8 January 2020 20 February 2020
Please cite this article as: Greco, A., FEM analysis of the elastic behavior of composites and nanocomposites with arbitrarily oriented reinforcements, Composite Structures (2020), doi: https://doi.org/10.1016/j.compstruct. 2020.112095
This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
© 2020 Published by Elsevier Ltd.
FEM analysis of the elastic behavior of composites and nanocomposites with arbitrarily oriented reinforcements
A.Greco Department of Engineering for Innovation, University of Salento, Via per Monteroni, 73100, Lecce, Italy. *Corresponding author Email: [email protected]
Author statement Antonio Greco: Conceptualization, Methodology, Software, Data curation, Writing- Original draft preparation, Investigation
Abstract This work is aimed at studying the effect of different morphological features on the stiffness of discontinuous fiber reinforced composites and nanocomposites by means of FEM analysis. The morphological features include not only volume fraction and aspect ratio of reinforcement, but also its orientation. The FEM approach, based on continuum mechanics, required the definition of a proper representative volume element, in which reinforcement is randomly dispersed in the simulation domain. The developed model has been used in order to calculate the stiffness matrix of unidirectional discontinuous glass fiber reinforced composite (DGFRC) as well as carbon nanotubes reinforced nanocomposites (CNTRNC). Comparison with closed form analytical models showed that the HT model is able to predict the evolution of the elastic constants over a wide range of volume fractions and aspect ratio. The TW model, in contrast, shows a poor agreement with the simulation results. On the other hand, the developed FEM model also allowed accounting for reinforcement orientation. The calculated elastic properties are in excellent good agreement with the corresponding values calculated according to the classical lamination theory, indicating the accuracy of the developed FEM model for estimation of stiffness of composites and nanocomposites, also accounting for reinforcement orientation. Keywords: FEM, discontinuous fibers, orientation
Introduction Carbon nanotubes (CNTs) are regarded as one of the most promising reinforcement materials for the next generation of high-performance structural and multifunctional composites [1]. Some CNTs are stronger than steel, lighter than aluminum, and more conductive than copper [2]. Theoretical and experimental studies have shown that CNTs exhibit extremely high tensile modulus (1 TPa) and strength (150 GPa). Due to this unique combination of physical and mechanical properties, CNTs have emerged as excellent candidates for use in next generation nanocomposites. However, in order to attain the best properties for the bulk composite, several issues related to the alignment, dispersion, aspect ratio, orientation, and load transfer need to be optimized. Two approaches have generally been adopted in treating this class of materials, namely, experimental and theoretical. Although there exists plenty of notable experimental efforts, these are severely hindered by the technical difficulties encountered during the manipulation, fabrication, and processing of CNTs and their composites. As such, there exists a huge variability amongst experimental results. This has motivated the use of theoretical and computational approaches to predict the potential mechanical properties of nano-reinforced composite systems [3]. In order to simulate the mechanical behavior of nanocomposites, several peculiar effects should be taken into account, which are not available in common microscale particles reinforced composites. Among them, one should mention the peculiar relevance of interface/interphases, of nanoreinforcement distribution and clustering and of aspect ratio on nanocomposite properties [4]. Therefore, due to the wide ranges of scales associated with CNT, a multi-scale modeling is required to predict mechanical properties of investigated materials. All involved length scales, starting from atomic scale and lasting to bulk material level, are scanned in this approach. Explanation of the relationships between physical conditions governing on atomic and molecular scale and the behavior of materials at the macro level depends on the development of models, in which all the nano, micro, meso and macro scales are considered [5]. In this regard, several nano-mechanical models are developed to calculate the mechanical properties of CNT reinforced composites. Computational approaches like molecular dynamics (MD) and continuum mechanics play significant role in the development of the CNT/polymer composites. The atomistic modeling techniques like classical MD simulations and density functional theory have focused on small length scale calculation. Molecular modeling enables the simulation of single CNTs using different approaches with varying degrees of simplification. This allows to calculate the properties of the equivalent cylinders [6]; however, MD cannot be directly transferred to multi-scale models.
On the other hand, continuum mechanics based methods can be applied to parametrically evaluate the effective elastic properties of CNT-composites. These methods can be used to correlate nanoscale to macro-scale by incorporating the necessary information from the atomic scale such as properties of chemical bonds and the interphase between the different materials. The majority of the reported continuum-mechanics methods use finite element analysis (FEA) models, which are based on Representative Volume Element (RVE) concept [7]. Modeling of nanocomposites is performed by choice of proper RVE, formed by a single or few cylinders surrounded by the matrix, and numerical estimation of the nanocomposite mechanical properties [8, 9, 10, 11]; the phenomena of load transfer from the host matrix to CNT can be accomplished through an interface region between the two [12]. Therefore, simulation of interface region plays a critical role in understanding the behavior of CNT based nanocomposites under different loading conditions [13]. However, multi scale modeling have also been developed [14], which directly correlates the fully atomistic region with the surrounding region modeled using continuum mechanics approach and FEM. However, this approach is severely limited by the choice of a very small RVE, which is usually represented as single CNT surrounded by the matrix. In a third, higher scale step, an homogenization procedure is developed in order to consider the nanocomposite as homogeneous material, requiring the mathematical modeling of the overall composite formed by CNTs surrounded by the matrix and the numerical prediction of the mechanical properties of the overall composite [15]. A more detailed review of the numerical approaches used for simulation of CNT nanocomposites, is presented in [16]. Traditional continuum techniques, which neglect the details of atomic structure of the materials, were also proposed. Although the traditional continuum models are incapable of accurately describing the interactions at the CNT-polymer interface, they are able to provide valuable information concerning the influence of the microscale morphology on the properties of CNT-reinforced composites, such as curvature, aspect ratio and volume fraction of the dispersed CNTs [17]. Both micromechanics and finite element methods are widely used in modeling the effective bulk properties of nanocomposites at the microscale level [18]. However, the developed continuum techniques are suitable to study perfectly aligned (parallel to the direction of loading) or randomly oriented nanoparticles [19]. All the developed models recognize that fiber volume fraction and aspect ratio are the parameters that mainly influence the elastic properties. The effect of partial fiber orientation, which can be related to processing method, is still not well understood. On the other hand, in recent works, Greco et al. [20, 21, 22,] introduced a FEM model, based on a random array of lamellae, coupled with a non-interpenetrating condition, which is able to study at a
meso-scale the diffusion behavior of uniformly dispersed, randomly positioned nanocomposites based on montmorillonite (MMT) lamellae, also accounting for the orientation of nanoparticles. The model has also been adapted to study multi-scale diffusion in partially intercalated MMT nanocomposites, also accounting for intra-stack diffusion [23, 24, 25]. Therefore, this work is aimed at studying by continuum mechanics, the evolution of the mechanical properties of a generally discontinuous fiber reinforced composite. This analysis, therefore, includes both discontinuous glass fiber reinforced composite (DGFRC) as well as carbon nanotubes reinforced nanocomposites (CNTRNC). The main morphological parameters are accounted in this analysis, also including reinforcement orientation. The obtained results have been compared with closed form analytical models. Mechanical modeling In order to model the mechanical properties of DGFRC and CNTRNC, the finite element (FE) solutions were obtained by using the Comsol Multiphysics software (version 3.5, Comsol AB, Sweden). 3D geometries were adopted, modeling the composite as an array of high aspect ratio particles, uniformly and randomly dispersed in the integration domain. A parallelepiped domain was used, and reinforcing fibers were represented as cylinders of diameter D and length T. The position of fibers or CNTs inside the domain was randomized by means of a properly developed subroutine of Matlab 6.5. The program, based on a Monte Carlo stochastic approach, produces random numbers and provides as an output a file that contains the center of reinforcing particles, including non interpenetrating conditions between them. The composite/nanocomposite is characterized by three parameters:
Aspect ratio of the reinforcing particles AR
Volume fraction of the reinforcing phase, F
Angle of orientation, defined as the angle formed by the cylinder axis with the x axis; o x=0, corresponding to unidirectional composite; all the reinforcing particles are aligned along the x axis, which is therefore regarded as the principal axis of the composite; o 𝜃𝑥 ≠ 0, corresponding to oriented composite. In order to account for random orientation in the y-z plane, the azimuthal angle 𝜃𝑦𝑧 is randomized, which leads to a quasi-isotropic behavior in the y-z plane; o x=54°, corresponding to completely random orientation, accomplished by random orientation of the azimuthal angle 𝜃𝑦𝑧 , leads to a fully isotropic behavior.
A sketch of the simulation domain is reported in Figure 1 a). The domain is composed of “external” boundaries, which are made by the external contour of the domain, and “internal” boundaries, which
are made by the interface between the reinforcing particles and the polymer. The Matlab subroutine also includes a non intersecting condition between each cylinder and the boundaries of the external domain.
z y
x Figure 1. sketch of the FEM domain The problem is solved by means of the 3D structural mechanics module, solid stress strain. For each simulation, a force per unit area P was applied on two opposing faces of the external boundary. Therefore, by applying P in the x, y, and z direction, the corresponding values of the modulus, E xx, Eyy and Ezz can be obtained according to: eq. 1)
E𝑖𝑖 =
𝑃𝐿0𝑖 ∆𝐿𝑖
Where L0i is the initial length along the ii direction, and Li the corresponding elongation. When reinforcement is oriented along the x axis, x=0, the axis x coincides with the principal direction, and is therefore labeled as 1, whereas directions y and z are labeled as 2 and 3. For each loading condition along the ii axis, the boundaries parallel to the direction of the applied load are considered to be stress-free: eq. 2)
𝜎𝑖𝑖 = ±𝑃 indicates loading applied in the ii direction, perpendicular to the i axis
eq. 3)
𝜎𝑗𝑖 = 0 indicates shear stress free on the faces perpendicular to the i axis
eq. 4)
𝜎𝑗𝑗 = 0 indicates normal stress free on the faces perpendicular to the j axis;
The assumption of stress free boundaries parallel to the direction of loading has been verified and will be discussed in the next section. The stress free boundary condition eq. 3) and eq. 4) allow to calculate, for each loading condition, the corresponding Poisson modulus: eq. 5)
𝜀
𝜈𝑗𝑖 = − 𝜀𝑗𝑗 = 𝑖𝑖
∆𝐿𝑗 𝐿0𝑖 𝐿0𝑗 ∆𝐿𝑖
In addition, simulations were also run in order to evaluate the shear modulus of the composite. In this case, the boundary conditions are as follows: eq. 6)
𝜎𝑖𝑗 = ±𝑃 indicates loading applied in the j direction on a plane perpendicular to the i
axis; in order to attain equilibrium to rotation, 𝜎𝑗𝑖 = ±𝑃 is also applied. eq. 7)
𝜎𝑖𝑖 = 𝜎𝑗𝑗 = 0 indicates normal stress free on the face parallel to the direction of the
applied load eq. 8)
𝑢𝑘𝑘 = 0 indicates no displacement of the domain in the third, i.e. k direction; this
condition is necessary in order to make sure that displacement of the domain only occurs in the i and the j direction; Accordingly, the shear modulus is obtained as: eq. 9)
G𝑖𝑗 =
𝑃𝐿0𝑖 ∆𝐿𝑗
The developed model also allows to calculate the I and II order coefficients of mutual influence, as defined by Jones in [26]. For modeling DGFRC, the modulus of fibers was taken equal to that of glass fibers, with isotropic properties EF=72 GPa and F=0.22. The modulus of the matrix was taken as EM=1 GPa and F=0.33; it is worth observing that, according to all the proposed models, the modulus of the matrix is expected to have no influence on the relative modulus of the nanocomposite. The geometry representation of CNTs neglects the internal cavity modeling them as cylinders; the values of the Hill’s elastic moduli [27] were used for calculating the equivalent engineering constants reported in Table I. Table I: elastic constants for CNT
E11CNT=(GPa)
2100
E22CNT=E33CNT
G12CNT=G13CNT
G23CNT
(GPa)
(GPa)
(GPa)
420
790
130
13CNT=12CNT 23CNT
0.1716
0.595
The developed FEM model was initially tested by simulating a continuous unidirectional, transversely isotropic composite; the theoretical stiffness matrix was obtained using the materials elastic constant and widely recognized classical micromechanics theory [26]. Each property of the composite, PC, was therefore obtained from the corresponding property of the CNT, PCNT, and the property of the matrix, Pm, by means of the well known Halpin-Tsai equations: eq. 10)
𝑃𝐶 = 𝑃𝑚
𝑃𝐶𝑁𝑇 (1+𝜉𝑣𝐶𝑁𝑇 )−𝜉𝑣𝐶𝑁𝑇 +𝜉 𝑃𝑀 𝑃𝐶𝑁𝑇 (1−𝑣𝐶𝑁𝑇 )+𝑣𝐶𝑁𝑇 +𝜉 𝑃𝑀
𝜈31 = 𝜈21 = 𝜈13 𝜈23 = 𝜈32 = 𝜈13 𝐺23 =
𝐸33 𝐸11
1 − 𝜈31 1 − 𝜈13
𝐸22 1 − 𝜈23
In eq. 10), according to [26], is a function of aspect ratio of the reinforcement, AR; =2AR for calculation of E11, and 12=13; =2 for E22=E33; =1 for G12= G13. For a volume fraction of fibers F=0.2 and fibers aligned in the x direction, the stiffness matrix determined by FEM analysis was in excellent agreement to the stiffness matrix built with the material elastic constant, as reported in Table II (coefficient of regression R2=0.99). In order to further confirm the accuracy of the developed FEM model, simulations were also run on continuous fibers aligned with an orientation angle measured with respect to the x axis, x=10° (anisotropic behavior of the composite). The theoretical stiffness matrix was obtained by the elastic constant calculated according to eq. 10) and further rotation of the stiffness matrix according to the classical lamination theory. Again, a comparison between the two stiffness matrix, reported in Table III, shows an excellent agreement between the results from simulation and theoretical predictions (coefficient of regression was R2=0.985) which again indicates the suitableness of the developed approach for the prediction of the elastic constant of composites. In particular, the very good agreement found between the off-diagonal terms (Q12, Q13,Q23) from FEM and analytical modeling clearly indicates that the boundary conditions of eq. 2)-eq. 4) and eq. 6)-eq. 8) provide an affordable estimate of the elastic constants of multi-phase composites. Table II: stiffness matrix for continuous unidirectional composite. Black: Halpi-Tsai model. Red: simulation results
𝑄11 = 2.3
𝑄12 = 0.83
𝑄13 = 0.83
𝑄14 = 0
𝑄15 = 0
𝑄16 = 0
𝑄11 = 2.2
𝑄12 = 0.83
𝑄13 = 0.81
𝑄14 < 2 ∗ 10−6
𝑄15 < 2 ∗ 10−6
𝑄16 < 2 ∗ 10−6
𝑄21 = 0.83
𝑄22 = 1.8
𝑄23 = 0.83
𝑄24 = 0
𝑄25 = 0
𝑄26 = 0
𝑄21 = 0.83
𝑄22 = 1.7
𝑄23 = 0.81
𝑄31 = 0.83
𝑄32 = 0.83
𝑄33 = 1.8
𝑄31 = 0.81
𝑄32 = 0.81
𝑄33 = 1.7
𝑄41 = 0
𝑄42 = 0
𝑄43 = 0
𝑄41 < 2 ∗
10−6
𝑄51 = 0 𝑄51 < 2 ∗
10−6
𝑄61 = 0 𝑄61 < 2 ∗
10−6
𝑄42 < 2 ∗
10−6
𝑄52 = 0 𝑄52 < 2 ∗
10−6
𝑄62 = 0 𝑄62 < 2 ∗
10−6
𝑄43 < 2
∗ 10−6
𝑄53 = 0 𝑄53 < 2
∗ 10−6
𝑄63 = 0 𝑄63 < 2
∗ 10−6
𝑄24 < 2
∗ 10−6
𝑄25 < 2
𝑄34 = 0 𝑄34 < 2 ∗ 10
𝑄35 = 0 −6
𝑄35 < 2 ∗ 10
𝑄44 = 0.46 𝑄44 = 0.46
𝑄45 < 2
∗ 10−6
𝑄64 = 0 𝑄64 < 2
∗ 10−6
𝑄26 < 2 ∗ 10−6 𝑄36 = 0
−6
𝑄45 = 0
𝑄54 = 0 𝑄54 < 2
∗ 10−6
∗ 10−6
𝑄36 < 2 ∗ 10−6 𝑄46 = 0 𝑄46 < 2 ∗ 10−6
𝑄55 = 0.46
𝑄56 = 0
𝑄55 = 0.46
𝑄56 < 2 ∗ 10−6
𝑄65 = 0
𝑄66 = 0.49
𝑄65 < 2
∗ 10−6
𝑄66 = 0.48
Table III: stiffness matrix for continuous oriented composite. Black: Halpi-Tsai model and stiffness matrix rotaion according to lamination theory. Red: simulation results 𝑄̅11 = 2.3
𝑄̅12 = 0.83
𝑄̅13 = 0.83
𝑄̅14 = 0
𝑄̅15 = 0
𝑄̅16 = 0
𝑄̅11 = 2.2
𝑄̅12 = 0.84
𝑄̅13 = 0.84
𝑄̅14 < 3 ∗ 10−6
𝑄̅15 < 3 ∗ 10−6
𝑄̅16 < 3 ∗ 10−6
𝑄̅21 = 0.83
𝑄̅22 = 1.8
𝑄̅23 = 0.85
𝑄̅24 = 0
𝑄̅25 = 0
𝑄̅26 = 0.0083
𝑄̅21 = 0.84
𝑄̅22 = 1.7
𝑄̅23 = 0.84
𝑄̅24 < 3 ∗ 10−6
𝑄̅25 < 3 ∗ 10−6
𝑄̅26 = 0.0051
𝑄̅31 = 0.83
𝑄̅32 = 0.85
𝑄̅33 = 1.8
𝑄̅34 = 0
𝑄̅35 = 0
𝑄̅36 = −0.09
𝑄̅31 = 0.82
𝑄̅32 = 0.83
𝑄̅33 = 1.7
𝑄̅41 = 0
𝑄̅42 = 0
𝑄̅43 = 0
𝑄̅41 < 3 ∗
10−6
𝑄̅42 < 3 ∗
10−6
𝑄̅43 < 3 ∗
10−6
𝑄̅34 < 3 ∗
10−6
𝑄̅35 < 3 ∗
10−6
𝑄̅44 = 0.48
𝑄̅45 = −6 ∗
10−3
𝑄̅44 = 0.46
𝑄̅45 = −3 ∗
10−3
𝑄̅36 = −0.07 𝑄̅46 = 0 𝑄̅46 < 3 ∗ 10−6
𝑄̅51 = 0
𝑄̅52 = 0
𝑄̅53 = 0
𝑄̅54 = −6 ∗ 10−3
𝑄̅55 = 0.46
𝑄̅56 = 0
𝑄̅51 < 3 ∗ 10−6
𝑄̅52 < 3 ∗ 10−6
𝑄̅53 < 3 ∗ 10−6
𝑄̅54 = −3 ∗ 10−3
𝑄̅55 = 0.45
𝑄̅56 < 3 ∗ 10−6
𝑄̅61 = 0
𝑄̅62 = 0.0083
𝑄̅63 = −0.09
𝑄̅64 = 0
𝑄̅65 = 0
𝑄̅66 = 0.49
𝑄̅61 < 3 ∗ 10−6
𝑄̅62 = 0.0048
𝑄̅63 = −0.06
𝑄̅64 < 3 ∗ 10−6
𝑄̅65 < 3 ∗ 10−6
𝑄̅66 = 0.48
Preliminary tests Initially, in order to properly develop the FEM model, some optimization trials were performed. As demonstrated by Chen and Papathanasiou [28], the simulation data are affected by a high scattering due to the stochastic nature of the reinforcing particle distribution. In order to reduce this scattering, the number of random distributions (Nrd) and the number of reinforcing paticles (Nf) were optimized by keeping a constant vf= 0.01 and the aspect ratio equal to 20 for perfectly aligned reinforcing particles. At first, the attention was focused on the minimum number of particles to insert in the domain in order to obtain a stable solution. The approach is similar to those adopted for studying diffusion in nanocomposites; basically, since stress free boundary conditions are applied on the boundaries parallel to the direction of the applied stress, the number of particles must be adequately high in order to neglect the border effects [20]. Results reported in Figure 2 clearly show that, with increasing
number of reinforcing particles in the domain, an increasing trend of the modulus is observed the modulus increases. However, the FEM results tend to an asymptotic value; keeping in mind the necessity to reduce the computation times, the number of particles was chosen, in each simulation, equal to 25. 1.60
E11 (GPa)
1.56 1.52 1.48 1.44 1.40
10
15 20 25 30 number of reinforcing particles
35
Figure 2: effect of the number of reinforcing particles on simulation results Afterwards, the number of random distribution was optimized; as highlighted in Figure 3, a low number of Nrd involves a very high standard deviation of the simulated longitudinal modulus; as the number of random distributions increases, the standard deviation decreases, and also the average value tends to level off, reaching a stable value at Nrd=10. 1.60 1.55
E11 (GPa)
1.50 1.45 1.40 1.35 1.30
0
5 10 number of random distributions
15
Figure 3: effect of the number of random distributions on simulation results Therefore, for each simulation, results reported are obtained as the average of 10 simulations, each one run with 25 reinforcing particles.
Finally, we focused on the effect of the domain shape; in order to explain its effect on the FEM results, two different bidimensional simulation were performed on the domains reported in Figure 4: in the first case, Figure 4 a), the length of the domain is equal to the length of reinforcing particle LF, increased by an arbitrarily low amount of matrix along the x axis, LV (necessary in order to avoid continuous fiber condition). The additional matrix layer is also highlighted in Figure 4 a) in light red. In the second domain, the length of the domain is increased in order to allow for positioning of two reinforcing particles along the x axis. The interparticle distance along the x axis, for symmetry reasons, is equal to 2LV. Also, in Figure 4 b), the interparticle region is highlighted in light green. A plot of the deformation x along the red line of Figure 4 a) and blue and green lines of Figure 4 b) is reported in Figure 5, together with the average value of the deformation and the stress along the x axis. If one single fiber would represent a satisfactory RVE, then, for symmetry reasons, the three plots would be expected to be perfectly overlapped. Instead, in the interparticle region, a higher deformation is observed; in facts, due to stress continuity at the fiber/matrix interface, a higher stress is also found in the same region, compared to the region LV. However, the higher stress found in the interparticle region contributes to reduce the stress, and therefore the deformation, far away from the fibers; for example, a plot along the purple lines of Figure 4 a) and b), reported in Figure 5 b), clearly shows the reduced stresses and deformation, particularly in the area which is highlighted in light yellow in Figure 4 a). As a consequence of this, the modulus increases with increasing number of fibers along the x axis, Nx. Simulations run on domains correspondent to those reported in Figure 4 a) and b), with increasing length and increasing number of fibers, are reported in Figure 6 a), clearly showing that at least 5 fibers are required in the x direction in order to obtain a stable solution. No error bar is reported in Figure 6 a), since the simulations were obtained on a regular array of fibers arranged according to the scheme of Figure 4; therefore, one single simulation was performed for each value of Nx.
Ldom=LF+2LV
LV
a)
2Ldom=2LF+4LV
LV
2LV
b) x Figure 4: 2D simulation domain with one or two fibers along the principal direction
1 fiber along LV
-7
1.8x10
-8
-7
1.6x10
5.0x10
x,avg=76.5 Pa
1.4x10
2 fibers along 2LV
-7
1.2x10
-8
4.0x10
-8
x
x,avg=10.0*10
-7
1.0x10
x,avg=97 Pa
-8
3.5x10
-8
-8
3.0x10
-8
2.5x10
8.0x10
-8
6.0x10
-8
-8
4.0x10
1 fiber 2 fibers
-8
4.5x10
2 fibers along LV
-7
x
-8
x,avg=8.2*10
0.0
0.2
0.4
0.6
0.8
2.0x10
1.0
a)
0.0
0.2
0.4
0.6
0.8
1.0
x/Ldom
x/LV
b)
Figure 5: deformation along the x axis calculated from 2D simulations; a) along the axis of reinforcing fiber; b) far away from the axis of reinforcing fiber In analyzing the effect of the shape of the domain in 3D simulation, it must be initially emphasized that in any case, for oriented nanoparticles, the minimum length of the domain the x direction, Ltotx, must satisfy: eq. 11)
𝐿𝑡𝑜𝑡𝑥 > 𝐴𝑅 ∗ 𝐷 ∗ 𝑐𝑜𝑠𝜃𝑥
Where D is the diameter of the reinforcing particle. On the other hand, the total volume of the domain is given as: eq. 12)
𝜋𝐷 3
𝑁
𝑉𝑡𝑜𝑡 = 𝜙𝑓 ∗ 𝐴𝑅 ∗
4
𝑓
Where F is the volume fraction of reinforcing phase. Under the constraint of eq. 11), the domain will also satisfy eq. 13)
𝑁 ∗𝜋𝐷 2
𝑉
𝐿𝑡𝑜𝑡𝑦 = 𝐿𝑡𝑜𝑡𝑧 < √𝐿 𝑡𝑜𝑡 = √4∗𝜙𝑓
𝑓 ∗𝑐𝑜𝑠𝜃
𝑡𝑜𝑡𝑥
Therefore, in general, the simulation domain will be a squared base parallelepiped with its own minimum aspect ratio, given as: eq. 14)
𝐿
𝐴𝑅𝑚𝑖𝑛,𝑑𝑜𝑚𝑎𝑖𝑛 = 𝐿𝑡𝑜𝑡𝑥 = 2
𝐴𝑅∗(𝜙𝐹 )1/2 ∗(𝑐𝑜𝑠𝜃)3/2
𝑡𝑜𝑡𝑦
√𝑁𝑓 ∗𝜋
On the other hand, if, according to the previous discussion, Nx particles need to be aligned in the x direction, then eq. 11) is modified as: eq. 15)
𝐿𝑡𝑜𝑡𝑥 = 𝑁𝑥 𝐴𝑅 ∗ 𝐷 ∗ 𝑐𝑜𝑠𝜃𝑥
Which substituted into eq. 12) and eq. 13) yields an actual aspect ratio of the domain: 1
eq. 16)
𝐴𝑅𝑑𝑜𝑚𝑎𝑖𝑛 = 2
𝐴𝑅∗(𝜙𝐹 )2 ∗(𝑁𝑥 )3/2 ∗(𝑐𝑜𝑠𝜃)3/2 √𝑁𝑓 ∗𝜋
The results of Figure 6 b) are a plot of the modulus as a function of the ratio between the actual aspect ratio of the domain and the minimum aspect ratio of the domain with randomly dispersed discontinuous reinforcing cylinders in 3D simulations. As clearly observed, the modulus increases with increasing aspect ratio of the domain. Even in this case, the modulus reaches a plateau value for sufficiently high values of 𝐴𝑅𝑑𝑜𝑚𝑎𝑖𝑛 , which is roughly 10 times the minimum aspect ratio. It is worth observing that, according to eq. 16),
𝐴𝑅𝑑𝑜𝑚𝑎𝑖𝑛 𝐴𝑅𝑚𝑖𝑛,𝑑𝑜𝑚𝑎𝑖𝑛
= (𝑁𝑥 )3/2, and therefore a ratio of 10 roughly
corresponds to 𝑁𝑥 = 4.7, which is a value very close to that observed in Figure 6 a) for bi-dimensional simulation. Based on the results of Figure 6 -Figure 6, in each domain the minimum length of the domain along the x axis was taken according to eq. 15). Finally, it is worth observing that, for high values of the aspect ratio of the reinforcing particles and volume fractions, the aspect ratio of the domain can reach very high values (about 200). In order to test the stability of the developed approach towards such an elongated domain, the simulation on unidirectional continuous fiber reinforced composite, was also performed with an aspect ratio of the domain equal to 200. The stability of the solution is confirmed by the fact that the coefficient of regression between the simulated stiffness matrix and the stiffness matrix estimated by Halpin-Tsai, eq. 10), is R2=0.9999.
The results obtained indicate that with an adequately high number of nanoparticles, boundary effects do not significantly influence the FEM results; therefore, the choice of a “free” boundary condition on the faces parallel to the direction of the applied load is justified. This is a distinct advantage of the proposed approach. In most of the reported works, in facts, the boundary condition on the faces parallel to the direction of the applied load is considered to be symmetric, which, however, does not allow for its contraction, yielding a null Poisson modulus. Instead the proposed approach allows for estimating the whole stiffness matrix of the composite. 4.04
7.2
4.02
6.8 6.4 6.0
3.98
E11 (GPa)
E11 (GPa)
4.00
3.96 3.94
5.2 4.8
3.92 3.90
5.6
4.4 1
2
3
Nx
4
5
6
4.0
2
4
6 8 10 ARdomain/ARmin,domain
12
14
Figure 6: effect of domain shape on longitudical modulus; a) for 2D simualtions; b) for 3 D simulations Simulation in DGFRC In Figure 7, the longitudinal modulus along the principal direction 1 (E11) is reported for unidirectional discontinuous glass fibers. Also, model prediction according to the Halpin-Tsai theory, obtained with 𝜉 = 2𝐴𝑅 [29], clearly shows an excellent good agreement with the FEM simulation data, as highlighted by the very high coefficient of regression between FEM and HT, R2=0.92. On the other hand, as reported in Figure 7 a), the Tandon-Weng (TW) model was also applied. This approach, which represents a simplification of the Mori-Tanaka method, significantly overestimates the longitudinal tensile modulus. For the calculation of the TW model, the Eshleby’s tensor coefficients were calculated according to the procedure for prolate ellipsoid reported in [30]. From these, the elastic coefficients of the composite were obtained according to the formula reported in [31]. The principal Poisson ratio was estimated by the non-iterative approach proposed in [32]. However, the TW model shows a very good agreement with the FEM simulation results at low values of F, but significantly overpredicts the FEM simulation results at high values of F. The other elastic constants, E22, G12 and 12 from FEM modeling are reported in Figure 7 b), c) and d), together with Halpin-Tsai and Tandon-Weng models prediction. The HT values were initially obtained with =2 for E22, =1 for G12 and =2AR for 12.
The transversal modulus E22 from TW is in quite good agreement with FEM simulation data (R2=0.85), whereas the HT model, not shown in Figure 7 b), significantly underpredicts the E22 from FEM (R2=0.27). On the other hand, when using as a fitting parameter, the best fitting of the FEM simulation data was found for =4.4, (R2=0.94) as shown in Figure 7 b). As reported in Figure 7 c), the HT model (=1) significantly underpredicts the G12 from FEM simulation (R2=0.49). The TW model, not shown in Figure 7 c), provides lower estimate than the HT model (R2=-0.44). In this case, the best fitting with the HT model was found using =2.5 (R2=0.99). Finally, as reported Figure 7 d), the 12 from FEM is much higher than the corresponding HT (=2AR) and TW models. In both cases the poor agreement between analytical models and simulation is highlighted by the very low R2=-75. However, in this case, using as a fitting parameter did not allow to improve the accuracy of model prediction. FEM simulation Halpin-Tsai =2AR Tandon-Weng
4
1.18 1.14
F=0.019
2
F=0.0064 F=0.0032
1.12
F=0.019
1.10 1.08 1.06 1.04
F=0.0064
1.02
F=0.0032
1.00 20
40
60
80
100 AR
120
140
20
160
a)
40
60
80
100 AR
120
140
160
b) FEM simulation Halpin-Tsai =1 Halpin-Tsai =2.5
0.43
0.330
F=0.032
F=0.0032 F=0.0064
0.329 F=0.019
0.40
12
G12 (GPa)
0.41
0.39 F=0.0064
0.38 0.37
FEM simulation Halpin-Tsai =2AR Tandon-Weng
0.331
0.42
c)
F=0.032
1.16
F=0.032
E22 (GPa)
E11 (GPa)
3
1
FEM simulation Halpin-Tsai =4.4 Tandon-Weng
1.20
40
60
80
100 AR
120
140
F=0.019
0.327
F=0.0032
20
0.328
F=0.032
0.326
160
20
40
60
80
100 AR
120
140
160
d) Figure 7: elastic constants for unidirectional DGFRC; a) longitudinal modulus; b) transverse modulus; c) longitudinal shear modulus; d) principal Poisson ratio
Afterward, the effect of the orientation angle was taken into account. Therefore, the main orientation angle (measured with respect to the x axis, and therefore corresponding to a rotation around the y or the z axis) was changed between 0 and 54 °. However, the simple rotation along the z axis, would
lead to a composite with Eyy>Ezz, as clearly observed in Figure 8 a); though this is a condition which could be encountered in polymer composites or nanocomposites, in many cases the existence of a single direction of flow during processing (extrusion, injection molding) causes orientation of reinforcing particles along one single direction (the x axis in our example), but quasi-isotropic properties in the other two directions (y and z in the present case). In order to attain this condition, the azimuthal angle yz (measured with respect to the y or the z axis, and therefore corresponding to a rotation around the x axis, as highlighted in Figure 8), was randomly varied between 0 and 180°. The resulting domain reported in Figure 8 b) guarantees quasi-isotropic behavior in the y-z plane. x x z
a)
yz
y
b)
Figure 8: sketch of the rotation angles in oriented composites The results obtained for the different elastic constants are reported in Figure 9 a)-d). In this case, no comparison with closed form analytical model is possible. In facts, both Halpi-Tsai and TandonWeng are not suitable to accurately predict the evolution of the elastic constants in unidirectional composites, which is the basis for calculating the elastic constants in oriented composites; however, for each F and AR, the elastic constants of oriented composites were calculated from the elastic constant of unidirectional composites, obtained from FEM analysis, followed by rotating the compliance matrix by an angle x around the y axis and further rotating by an angle yz=45° around the x axis, which corresponds to a random orientation also in the y-z plane [33]. The azimuthal angle does not affect the longitudinal modulus Exx, but it affects the other elastic constants reported in Figure 9. In all cases, a very good agreement is observed between the FEM simulated constants and constants obtained from properties calculated in the principal direction followed by matrix rotation.
In particular, for Exx, a monotonic behavior with increasing orientation angle is found; for x=54 °, the value of Exx is lower than the corresponding value of E22; this is due to the fact that for all the simulated GFRP, the value of G12 is lower than
𝐸11 [26]. 𝐸 2( 11 +𝜈12 ) 𝐸22
However, the FEM simulation is able
to capture this behavior of the GFRP. For the same reason, in the range between 0 ° and 45 °, the value of Eyy decreases with increasing orientation angle. Even in this case, the FEM simulation is able to highlight the decrease of Eyy. The FEM model is also able to predict the evolution of Gxy, and of xy with their corresponding maximum. FEM simulation prediction from elastic constants along the principal axes and formula for rotation
1.8
AR=20
F=0.0064
AR=80
1.4
F=0.0064
AR=20 =0.0032 F
1.2
F=0.0019
1.08
Eyy (GPa)
Exx (GPa)
1.10
F=0.0019
1.6
FEM simulation prediction from elastic constants along the principal axes and formula for rotation
AR=40
F=0.0032
1.06
AR=20
F=0.0064
AR=80
F=0.0064
AR=20
1.04
F=0.00032
AR=40
F=0.0032
AR=20
1.02
AR=20 1.0 0
10
20 30 40 orientation angle (°)
50
a)
0
10
20 30 40 orientation angle (°)
50
60
b) FEM simulation prediction from elastic constants along the principal axes and formula for rotation
0.350
F=0.019
AR=20
0.435 0.420
FEM simulation prediction from elastic constants along the principal axes and formula for rotation
0.345 F=0.0064
F=0.0064
AR=80
AR=20 F=0.0032
0.405
AR=40
0.390
F=0.0032
0
10
20 30 40 orientation angle (°)
50
c)
F=0.0064
AR=20 0.335 F=0.0032 =0.0032 F
0.330 0.325
60
AR=80
AR=20
AR=40
AR=20
0.375
F=0.0064
F=0.019
0.340 xy
0.450
Gxy (GPa)
1.00
60
0
10
AR=20
20 30 40 orientation angle (°)
50
60
d)
Figure 9: elastic constants for GFRP at different orientation angles. a) longitudinal modulus; b) transverse modulus; c) longitudinal shear modulus; d) principal Poisson ratio. black F=0.0032, AR=20; red F=0.0032, AR=40; blue F=0.0064, AR=20; magenta F=0.0064, AR=80; olive F=0.019, AR=20 Simulation in CNT nanocomposites The longitudinal modulus for fully oriented nanoparticles, together with Halpin-Tsai prediction obtained with 𝜉 = 2𝐴𝑅, is reported in Figure 10 a). As in the case of glass fibers, the HT theory is in reasonable good agreement with the FEM simulation data. For CNTRNC, no comparison is provided
with Tandon-Weng, since closed form solutions for this approach are only available for isotropic reinforcement [34]. The results for the other elastic constant, reported in Figure 10 b), c) and d), substantially confirm the behavior already discussed for glass fibers: the HT model underpredicts both E22 (=2) and G12 (=1) from FEM; 12 from FEM is also much higher than the corresponding HT (=2AR). For E22, when considering as a fitting parameter, the best fit is obtained with =5.4; for G12, the best fit is obtained with =2.2. FEM simulation Halpin-Tsai =2AR
FEM simulation Halpin-Tsai =2 Halpin-Tsai =5.5
1.25 F=0.032
1.20
F=0.032
1.15
F=0.019
F=0.019
4
F=0.0064
3
E22 (GPa)
E11 (GPa)
10 9 8 7 6 5
F=0.0032
2
F=0.0064
1.10
F=0.0032
1.05 1
20
40
60
80
100 AR
120
140
a)
40
60
80
100 AR
120
140
0.330
F=0.032
160
0.405 0.400
F=0.0032
0.329
FEM simulation Halpin-Tsai =1 Halpin-Tsai =2.2
0.410
F=0.0064
FEM simulation Halpin-Tsai =2AR
0.328 F=0.019 12
G12 (GPa)
20
b) 0.415
0.395
0.327
F=0.019
0.326
0.390 0.385
F=0.0064
0.380
F=0.0032
20
c)
1.00
160
40
60
80
100 AR
120
140
0.325 0.324
160
F=0.032
20
40
60
80
100 AR
120
140
160
d)
Figure 10: elastic constants for unidirectional CNTRNC; a) transverse modulus; b) longitudinal shear modulus; c) principal Poisson ratio Also in the case of CNT nanocomposites, the effect of the orientation angle was taken into account. Following the same procedure used for GFRP, the results obtained for the different elastic constants are reported in Figure 11 a)-d). For comparison purposes, the elastic constant were also calculated from the elastic constant of unidirectional composite, followed by rotating the compliance matrix by an angle x around the y axis and further rotating by an angle yz=45 ° around the x axis. As in the case of GFRP, a very good agreement is observed between the FEM simulated constants and constants obtained from properties calculated in the principal direction followed by matrix rotation.
FEM simulation prediction from elastic constants along the principal axes and formula for rotation
3.6 3.2
F=0.019
1.14
AR=40
1.12
F=0.0064
F=0.019
1.10
AR=80
AR=40
2.8
F=0.0064
AR=80
2.4
F=0.0064
AR=40
2.0
Eyy (GPa)
Exx (GPa)
FEM simulation prediction from elastic constants along the principal axes and formula for rotation
1.16
F=0.0032
AR=40
1.6
F=0.0032
F=0.0064
1.06
AR=40
1.04
AR=20
1.2
1.08 F=0.0032
F=0.0032
AR=40
AR=20
1.02 0
10
20 30 40 orientation angle (°)
50
0
60
a)
10
20 30 40 orientation angle (°)
50
60
b)
0.52
FEM simulation prediction from elastic constants along the principal axes and formula for rotation F=0.019
0.50
AR=40
Gxy (GPa)
0.48 0.46
0.365
F=0.0064
0.360
AR=80
0.355
F=0.0064
0.350
AR=40
0.44
0.370
xy
0.54
F=0.0032
0.42
F=0.0032
0.38
AR=20 0
c)
10
20 30 40 orientation angle (°)
50
F=0.019
AR=40 F=0.0064
F=0.0064
AR=40
AR=80
0.340
AR=40
0.40
0.345
FEM simulation prediction from elastic constants along the principal axes and formula for rotation
0.335
F=0.0032
0.330
AR=40
0.325
60
0
10
20 30 40 orientation angle(°)
F=0.0032
AR=20 50
60
d)
Figure 11: elastic constants for CNTRNC at different orientation angles: a) longitudinal modulus; b) transverse modulus; c) longitudinal shear modulus; d) principal Poisson ratio. black F=0.005, AR=20; red F=0.005, AR=40; blue F=0.01, AR=40; magenta F=0.01, AR=80; olive F=0.03, AR=40 Conclusions In this work, a geometrical FEM model previously derived to study diffusion in polymer nanocomposites has been properly modified in order to study the effect of different morphological features on the stiffness evolution in short fiber reinforced composites and nanocomposites. Compared to similar works reported in literature, the developed model can account for different morphological features, including processing induced fiber orientation. The accuracy of the developed approach has been tested by comparing the simulation results with closed form solutions for continuous fiber reinforced composites, also accounting for orientation of the fibers. After, the FEM model has been used in order to calculated the stiffness matrix of discontinuous glass fiber reinforced composite (DGFRC) as well as carbon nanotubes reinforced nanocomposites (CNTRNC).
For DGFRC, comparison with closed form analytical models (Halpin-Tsai and Tandon-Weng models) showed that HT is able to predict the evolution of elastic constants, provided that the geometric factor is sued as a fitting parameter. On the other hand, the TW model does not provide an adequate estimate for the elastic constants of the DGFRC. For CNTRNC, comparison with Tandon-Weng model was not performed, since closed form solutions were only found for isotropic reinforcement; however, comparison with the Halpin-Tsai model substantially confirmed the results found for glass fibers. Finally, the developed FEM model was used for calculating the stiffness at different orientation of the reinforcement. In this case, results were not compared with classical micromechanical approaches, since these latter were proven to be unsuitable to predict the evolution of longitudinal modulus even for unidirectional reinforcement. Instead, for each value of volume fraction and aspect ratio of the reinforcement, the longitudinal, transversal, shear and Poisson modulus were calculated from the corresponding values obtained for uni-directional composites from FEM, followed by stiffness matrix rotation according to classical lamination theory. The elastic properties from FEM analysis are in excellent good agreement with the corresponding values calculated according to the classical lamination theory, indicating the accuracy of the developed FEM model for estimation of stiffness of composites and nanocomposites, also accounting for reinforcement orientation.
References 1 Endo M, Hayashi T, Kim YA, Terrones M, Dresselhaus MS. Applications of carbon nanotubes in the twenty-first century. Philos Trans R Soc London 2004; A 362: 2223–2238 2 Moniruzzaman M, Winey KI. Polymer nanocomposites containing carbon nanotubes. Macromolecules 2006; 39: 5194–5205. 3 Wernik JM, Meguid SA. Multiscale micromechanical modeling of the constitutive response of carbon nanotube-reinforced structural adhesives. Int J Solids Struct 2014; 51: 2575–2589 4 Wang HW, Zhou HW, Peng RD, Mishnaevsky L Jr. Nanoreinforced polymer composites: 3D FEM modeling with effective interface concept. Compos Sci Technol 2011; 71: 980–988 5 Rafiee R, Firouzbakht V. Multi-scale modeling of carbon nanotube reinforced polymers using irregular tessellation technique. Mech Mater 2014; 78: 74–84 6 Schiebold M, Schmidt H, Mehner J. A Finite Element Approach for Modeling and Simulation of CNT/Polymer Composites. Phys Status Solidi A 2019; 1800952 7 Kumar P, Srinivas J. Numerical evaluation of effective elastic properties of CNT-reinforced polymers for interphase effects. Comput Mater Sci 2014; 88: 139–144
8 Tserpes KI, Chanteli A. Parametric numerical evaluation of the effective elastic properties of carbon nanotube-reinforced polymers. Compos Struct 2013; 99: 366–374 9 Garcia-Macias E, Guzman CF, Saavedra Flores EI. Castro-Triguero R. Multiscale modeling of the elastic moduli of CNT-reinforced polymers and fitting of efficiency parameters for the use of the extended rule-of-mixtures. Compos Part B-Eng 2019; 159: 114–131 10 Wang JF, Liew KM. On the study of elastic properties of CNT-reinforced composites based on element-free MLS method with nanoscale cylindrical representative volume element. Compos Struct 2015; 124: 1–9 11 Joshi P, Upadhyay SH. Evaluation of elastic properties of multi walled carbon nanotube reinforced composite. Comput Mater Sci 2014; 81: 332–338 12 Shahzad Zuberi MJ, Esat V. Investigating the mechanical properties of single walled carbon nanotube reinforced epoxy composite through finite element modelling. Compos Part B-Eng 2015; 71: 1–9 13 Pantano A, Mantione P. A numerical-analytical model for the characterization of composites reinforced by carbon nanotubes. Appl Phys A 2010; 99: 895–902 14 Gupta AK, Harsh SP. Analysis of mechanical properties of carbon nanotube reinforced polymer composites using multi-scale finite element modeling approach. Compos Part B-Eng 2016; 95: 172178 15 Nie J, Jia Y, Qu P, Shi Q. Carbon Nanotube/Carbon Fiber Multiscale Composite: Influence of Interfacial Strength on Mechanical Properties. J Inorg Organomet Polym 2011; 21: 937–940 16 Yengejeh SI, Kazemi SA, Öchsner A. Carbon nanotubes as reinforcement in composites: A review of the analytical, numerical and experimental approaches. Comput Mater Sci 2017; 136: 85–101 17 Lusti HR, Gusev AA. Finite element predictions for the thermoelastic properties of nanotube reinforced polymers. Model Simul Mater Sci Eng 2004; 12: S107-S119 18 Chen X, Alian AR, Meguid SA. Modeling of CNT-reinforced nanocomposite with complex morphologies using modified embedded finite element technique. Compos Struct 2019; 227: 111329 19 Babu KP, Mohite PM, Upadhyay CS. Development of an RVE and its stiffness predictions based on mathematical homogenization theory for short fibre composites. Int J Solids Struct 2018; 130– 131: 80–104 20 Greco A, Maffezzoli A. Two-dimensional and three-dimensional simulation of diffusion in nanocomposite with arbitrarily oriented lamellae. J Membr Sci 2013; 442: 238-244 21 Greco A. Numerical simulation and mathematical modeling of 2D multi-scale diffusion in lamellar nanocomposite. Comput Mater Sci 2014; 90: 203–209
22 Greco A. Simulation and modeling of diffusion in oriented lamellar nanocomposites. Comput Mater Sci 2014; 83: 164–170 23 Greco A, Maffezzoli A. Finite Element Modeling of Multiscale Diffusion in Intercalated Nanocomposites. J Nanomater 2015; 482698 24 Greco A, Maffezzoli A. Finite element simulation and analytical modeling of 3D multi scale diffusion in nanocomposites with permeable stacks. Model Simul Mater Sci Eng 2016; 24 (1): 015003 25 Greco A, Esposito Corcione C, Maffezzoli A. Effect of multi scale diffusion on the permeability behaviour of intercalated nanocomposites. J Membr Sci 2016; 505: 92–99 26 Jones RM, Mechanics of Composite Materials, 2nd ed, 1999, Taylor & Francis, Philadelphia, US 27 Tornabene F, Bacciocchi M. Fantuzzi N, Reddy JN. Multiscale Approach for Three-Phase CNT/Polymer/Fiber Laminated Nanocomposite Structures. Polym Compos 2019; 40: E102-E126. 28 Chen X, Papathanasiou TD. Micro-scale modeling of axial flow through unidirectional disordered fiber arrays, Compos Sci Technol 2007; 67: 1286-1293 29 Halpin JC, Kardos JL. The Halpin-Tsai equations: a review. Polym Eng Sci 1976; 16(5): 344-352 30 David EC, Zimmerman RW, Compressibility and shear compliance of spheroidal pores: Exact derivation via the Eshelby tensor, and asymptotic expressions in limiting cases. Int J Solids Struct 2011; 48: 680–686 31 Zare Y, Development of simplified Tandon-Weng solutions of Mori-Tanaka theory for Young’s modulus of polymer nanocomposites considering the interphase. J Appl Polym Sci 2016: 43816 32 Tucker III CL, Liang E. Stiffness predictions for unidirectional short-fiber composites: Review and evaluation, Compos Sci Technol 1999; 59 (5): 655-671 33 Mortazavian S, Fatemi A. Effects of fiber orientation and anisotropy on tensile strength and elastic modulus of short fiber reinforced polymer composites. Compos Part B-Eng 2015; 72: 116–129 34 Tandon GP, Weng GJ. The effect of aspect ratio of inclusions on the elastic properties of unidirectionally aligned composites, Polym Compos 1984; 5(4): 327-333
S0263-8223(19)34105-4 https://doi.org/10.1016/j.compstruct.2020.112095 COST 112095
To appear in:
Composite Structures
Received Date: Revised Date: Accepted Date:
28 October 2019 8 January 2020 20 February 2020
Please cite this article as: Greco, A., FEM analysis of the elastic behavior of composites and nanocomposites with arbitrarily oriented reinforcements, Composite Structures (2020), doi: https://doi.org/10.1016/j.compstruct. 2020.112095
This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
© 2020 Published by Elsevier Ltd.
FEM analysis of the elastic behavior of composites and nanocomposites with arbitrarily oriented reinforcements
A.Greco Department of Engineering for Innovation, University of Salento, Via per Monteroni, 73100, Lecce, Italy. *Corresponding author Email: [email protected]
Author statement Antonio Greco: Conceptualization, Methodology, Software, Data curation, Writing- Original draft preparation, Investigation
Abstract This work is aimed at studying the effect of different morphological features on the stiffness of discontinuous fiber reinforced composites and nanocomposites by means of FEM analysis. The morphological features include not only volume fraction and aspect ratio of reinforcement, but also its orientation. The FEM approach, based on continuum mechanics, required the definition of a proper representative volume element, in which reinforcement is randomly dispersed in the simulation domain. The developed model has been used in order to calculate the stiffness matrix of unidirectional discontinuous glass fiber reinforced composite (DGFRC) as well as carbon nanotubes reinforced nanocomposites (CNTRNC). Comparison with closed form analytical models showed that the HT model is able to predict the evolution of the elastic constants over a wide range of volume fractions and aspect ratio. The TW model, in contrast, shows a poor agreement with the simulation results. On the other hand, the developed FEM model also allowed accounting for reinforcement orientation. The calculated elastic properties are in excellent good agreement with the corresponding values calculated according to the classical lamination theory, indicating the accuracy of the developed FEM model for estimation of stiffness of composites and nanocomposites, also accounting for reinforcement orientation. Keywords: FEM, discontinuous fibers, orientation
Introduction Carbon nanotubes (CNTs) are regarded as one of the most promising reinforcement materials for the next generation of high-performance structural and multifunctional composites [1]. Some CNTs are stronger than steel, lighter than aluminum, and more conductive than copper [2]. Theoretical and experimental studies have shown that CNTs exhibit extremely high tensile modulus (1 TPa) and strength (150 GPa). Due to this unique combination of physical and mechanical properties, CNTs have emerged as excellent candidates for use in next generation nanocomposites. However, in order to attain the best properties for the bulk composite, several issues related to the alignment, dispersion, aspect ratio, orientation, and load transfer need to be optimized. Two approaches have generally been adopted in treating this class of materials, namely, experimental and theoretical. Although there exists plenty of notable experimental efforts, these are severely hindered by the technical difficulties encountered during the manipulation, fabrication, and processing of CNTs and their composites. As such, there exists a huge variability amongst experimental results. This has motivated the use of theoretical and computational approaches to predict the potential mechanical properties of nano-reinforced composite systems [3]. In order to simulate the mechanical behavior of nanocomposites, several peculiar effects should be taken into account, which are not available in common microscale particles reinforced composites. Among them, one should mention the peculiar relevance of interface/interphases, of nanoreinforcement distribution and clustering and of aspect ratio on nanocomposite properties [4]. Therefore, due to the wide ranges of scales associated with CNT, a multi-scale modeling is required to predict mechanical properties of investigated materials. All involved length scales, starting from atomic scale and lasting to bulk material level, are scanned in this approach. Explanation of the relationships between physical conditions governing on atomic and molecular scale and the behavior of materials at the macro level depends on the development of models, in which all the nano, micro, meso and macro scales are considered [5]. In this regard, several nano-mechanical models are developed to calculate the mechanical properties of CNT reinforced composites. Computational approaches like molecular dynamics (MD) and continuum mechanics play significant role in the development of the CNT/polymer composites. The atomistic modeling techniques like classical MD simulations and density functional theory have focused on small length scale calculation. Molecular modeling enables the simulation of single CNTs using different approaches with varying degrees of simplification. This allows to calculate the properties of the equivalent cylinders [6]; however, MD cannot be directly transferred to multi-scale models.
On the other hand, continuum mechanics based methods can be applied to parametrically evaluate the effective elastic properties of CNT-composites. These methods can be used to correlate nanoscale to macro-scale by incorporating the necessary information from the atomic scale such as properties of chemical bonds and the interphase between the different materials. The majority of the reported continuum-mechanics methods use finite element analysis (FEA) models, which are based on Representative Volume Element (RVE) concept [7]. Modeling of nanocomposites is performed by choice of proper RVE, formed by a single or few cylinders surrounded by the matrix, and numerical estimation of the nanocomposite mechanical properties [8, 9, 10, 11]; the phenomena of load transfer from the host matrix to CNT can be accomplished through an interface region between the two [12]. Therefore, simulation of interface region plays a critical role in understanding the behavior of CNT based nanocomposites under different loading conditions [13]. However, multi scale modeling have also been developed [14], which directly correlates the fully atomistic region with the surrounding region modeled using continuum mechanics approach and FEM. However, this approach is severely limited by the choice of a very small RVE, which is usually represented as single CNT surrounded by the matrix. In a third, higher scale step, an homogenization procedure is developed in order to consider the nanocomposite as homogeneous material, requiring the mathematical modeling of the overall composite formed by CNTs surrounded by the matrix and the numerical prediction of the mechanical properties of the overall composite [15]. A more detailed review of the numerical approaches used for simulation of CNT nanocomposites, is presented in [16]. Traditional continuum techniques, which neglect the details of atomic structure of the materials, were also proposed. Although the traditional continuum models are incapable of accurately describing the interactions at the CNT-polymer interface, they are able to provide valuable information concerning the influence of the microscale morphology on the properties of CNT-reinforced composites, such as curvature, aspect ratio and volume fraction of the dispersed CNTs [17]. Both micromechanics and finite element methods are widely used in modeling the effective bulk properties of nanocomposites at the microscale level [18]. However, the developed continuum techniques are suitable to study perfectly aligned (parallel to the direction of loading) or randomly oriented nanoparticles [19]. All the developed models recognize that fiber volume fraction and aspect ratio are the parameters that mainly influence the elastic properties. The effect of partial fiber orientation, which can be related to processing method, is still not well understood. On the other hand, in recent works, Greco et al. [20, 21, 22,] introduced a FEM model, based on a random array of lamellae, coupled with a non-interpenetrating condition, which is able to study at a
meso-scale the diffusion behavior of uniformly dispersed, randomly positioned nanocomposites based on montmorillonite (MMT) lamellae, also accounting for the orientation of nanoparticles. The model has also been adapted to study multi-scale diffusion in partially intercalated MMT nanocomposites, also accounting for intra-stack diffusion [23, 24, 25]. Therefore, this work is aimed at studying by continuum mechanics, the evolution of the mechanical properties of a generally discontinuous fiber reinforced composite. This analysis, therefore, includes both discontinuous glass fiber reinforced composite (DGFRC) as well as carbon nanotubes reinforced nanocomposites (CNTRNC). The main morphological parameters are accounted in this analysis, also including reinforcement orientation. The obtained results have been compared with closed form analytical models. Mechanical modeling In order to model the mechanical properties of DGFRC and CNTRNC, the finite element (FE) solutions were obtained by using the Comsol Multiphysics software (version 3.5, Comsol AB, Sweden). 3D geometries were adopted, modeling the composite as an array of high aspect ratio particles, uniformly and randomly dispersed in the integration domain. A parallelepiped domain was used, and reinforcing fibers were represented as cylinders of diameter D and length T. The position of fibers or CNTs inside the domain was randomized by means of a properly developed subroutine of Matlab 6.5. The program, based on a Monte Carlo stochastic approach, produces random numbers and provides as an output a file that contains the center of reinforcing particles, including non interpenetrating conditions between them. The composite/nanocomposite is characterized by three parameters:
Aspect ratio of the reinforcing particles AR
Volume fraction of the reinforcing phase, F
Angle of orientation, defined as the angle formed by the cylinder axis with the x axis; o x=0, corresponding to unidirectional composite; all the reinforcing particles are aligned along the x axis, which is therefore regarded as the principal axis of the composite; o 𝜃𝑥 ≠ 0, corresponding to oriented composite. In order to account for random orientation in the y-z plane, the azimuthal angle 𝜃𝑦𝑧 is randomized, which leads to a quasi-isotropic behavior in the y-z plane; o x=54°, corresponding to completely random orientation, accomplished by random orientation of the azimuthal angle 𝜃𝑦𝑧 , leads to a fully isotropic behavior.
A sketch of the simulation domain is reported in Figure 1 a). The domain is composed of “external” boundaries, which are made by the external contour of the domain, and “internal” boundaries, which
are made by the interface between the reinforcing particles and the polymer. The Matlab subroutine also includes a non intersecting condition between each cylinder and the boundaries of the external domain.
z y
x Figure 1. sketch of the FEM domain The problem is solved by means of the 3D structural mechanics module, solid stress strain. For each simulation, a force per unit area P was applied on two opposing faces of the external boundary. Therefore, by applying P in the x, y, and z direction, the corresponding values of the modulus, E xx, Eyy and Ezz can be obtained according to: eq. 1)
E𝑖𝑖 =
𝑃𝐿0𝑖 ∆𝐿𝑖
Where L0i is the initial length along the ii direction, and Li the corresponding elongation. When reinforcement is oriented along the x axis, x=0, the axis x coincides with the principal direction, and is therefore labeled as 1, whereas directions y and z are labeled as 2 and 3. For each loading condition along the ii axis, the boundaries parallel to the direction of the applied load are considered to be stress-free: eq. 2)
𝜎𝑖𝑖 = ±𝑃 indicates loading applied in the ii direction, perpendicular to the i axis
eq. 3)
𝜎𝑗𝑖 = 0 indicates shear stress free on the faces perpendicular to the i axis
eq. 4)
𝜎𝑗𝑗 = 0 indicates normal stress free on the faces perpendicular to the j axis;
The assumption of stress free boundaries parallel to the direction of loading has been verified and will be discussed in the next section. The stress free boundary condition eq. 3) and eq. 4) allow to calculate, for each loading condition, the corresponding Poisson modulus: eq. 5)
𝜀
𝜈𝑗𝑖 = − 𝜀𝑗𝑗 = 𝑖𝑖
∆𝐿𝑗 𝐿0𝑖 𝐿0𝑗 ∆𝐿𝑖
In addition, simulations were also run in order to evaluate the shear modulus of the composite. In this case, the boundary conditions are as follows: eq. 6)
𝜎𝑖𝑗 = ±𝑃 indicates loading applied in the j direction on a plane perpendicular to the i
axis; in order to attain equilibrium to rotation, 𝜎𝑗𝑖 = ±𝑃 is also applied. eq. 7)
𝜎𝑖𝑖 = 𝜎𝑗𝑗 = 0 indicates normal stress free on the face parallel to the direction of the
applied load eq. 8)
𝑢𝑘𝑘 = 0 indicates no displacement of the domain in the third, i.e. k direction; this
condition is necessary in order to make sure that displacement of the domain only occurs in the i and the j direction; Accordingly, the shear modulus is obtained as: eq. 9)
G𝑖𝑗 =
𝑃𝐿0𝑖 ∆𝐿𝑗
The developed model also allows to calculate the I and II order coefficients of mutual influence, as defined by Jones in [26]. For modeling DGFRC, the modulus of fibers was taken equal to that of glass fibers, with isotropic properties EF=72 GPa and F=0.22. The modulus of the matrix was taken as EM=1 GPa and F=0.33; it is worth observing that, according to all the proposed models, the modulus of the matrix is expected to have no influence on the relative modulus of the nanocomposite. The geometry representation of CNTs neglects the internal cavity modeling them as cylinders; the values of the Hill’s elastic moduli [27] were used for calculating the equivalent engineering constants reported in Table I. Table I: elastic constants for CNT
E11CNT=(GPa)
2100
E22CNT=E33CNT
G12CNT=G13CNT
G23CNT
(GPa)
(GPa)
(GPa)
420
790
130
13CNT=12CNT 23CNT
0.1716
0.595
The developed FEM model was initially tested by simulating a continuous unidirectional, transversely isotropic composite; the theoretical stiffness matrix was obtained using the materials elastic constant and widely recognized classical micromechanics theory [26]. Each property of the composite, PC, was therefore obtained from the corresponding property of the CNT, PCNT, and the property of the matrix, Pm, by means of the well known Halpin-Tsai equations: eq. 10)
𝑃𝐶 = 𝑃𝑚
𝑃𝐶𝑁𝑇 (1+𝜉𝑣𝐶𝑁𝑇 )−𝜉𝑣𝐶𝑁𝑇 +𝜉 𝑃𝑀 𝑃𝐶𝑁𝑇 (1−𝑣𝐶𝑁𝑇 )+𝑣𝐶𝑁𝑇 +𝜉 𝑃𝑀
𝜈31 = 𝜈21 = 𝜈13 𝜈23 = 𝜈32 = 𝜈13 𝐺23 =
𝐸33 𝐸11
1 − 𝜈31 1 − 𝜈13
𝐸22 1 − 𝜈23
In eq. 10), according to [26], is a function of aspect ratio of the reinforcement, AR; =2AR for calculation of E11, and 12=13; =2 for E22=E33; =1 for G12= G13. For a volume fraction of fibers F=0.2 and fibers aligned in the x direction, the stiffness matrix determined by FEM analysis was in excellent agreement to the stiffness matrix built with the material elastic constant, as reported in Table II (coefficient of regression R2=0.99). In order to further confirm the accuracy of the developed FEM model, simulations were also run on continuous fibers aligned with an orientation angle measured with respect to the x axis, x=10° (anisotropic behavior of the composite). The theoretical stiffness matrix was obtained by the elastic constant calculated according to eq. 10) and further rotation of the stiffness matrix according to the classical lamination theory. Again, a comparison between the two stiffness matrix, reported in Table III, shows an excellent agreement between the results from simulation and theoretical predictions (coefficient of regression was R2=0.985) which again indicates the suitableness of the developed approach for the prediction of the elastic constant of composites. In particular, the very good agreement found between the off-diagonal terms (Q12, Q13,Q23) from FEM and analytical modeling clearly indicates that the boundary conditions of eq. 2)-eq. 4) and eq. 6)-eq. 8) provide an affordable estimate of the elastic constants of multi-phase composites. Table II: stiffness matrix for continuous unidirectional composite. Black: Halpi-Tsai model. Red: simulation results
𝑄11 = 2.3
𝑄12 = 0.83
𝑄13 = 0.83
𝑄14 = 0
𝑄15 = 0
𝑄16 = 0
𝑄11 = 2.2
𝑄12 = 0.83
𝑄13 = 0.81
𝑄14 < 2 ∗ 10−6
𝑄15 < 2 ∗ 10−6
𝑄16 < 2 ∗ 10−6
𝑄21 = 0.83
𝑄22 = 1.8
𝑄23 = 0.83
𝑄24 = 0
𝑄25 = 0
𝑄26 = 0
𝑄21 = 0.83
𝑄22 = 1.7
𝑄23 = 0.81
𝑄31 = 0.83
𝑄32 = 0.83
𝑄33 = 1.8
𝑄31 = 0.81
𝑄32 = 0.81
𝑄33 = 1.7
𝑄41 = 0
𝑄42 = 0
𝑄43 = 0
𝑄41 < 2 ∗
10−6
𝑄51 = 0 𝑄51 < 2 ∗
10−6
𝑄61 = 0 𝑄61 < 2 ∗
10−6
𝑄42 < 2 ∗
10−6
𝑄52 = 0 𝑄52 < 2 ∗
10−6
𝑄62 = 0 𝑄62 < 2 ∗
10−6
𝑄43 < 2
∗ 10−6
𝑄53 = 0 𝑄53 < 2
∗ 10−6
𝑄63 = 0 𝑄63 < 2
∗ 10−6
𝑄24 < 2
∗ 10−6
𝑄25 < 2
𝑄34 = 0 𝑄34 < 2 ∗ 10
𝑄35 = 0 −6
𝑄35 < 2 ∗ 10
𝑄44 = 0.46 𝑄44 = 0.46
𝑄45 < 2
∗ 10−6
𝑄64 = 0 𝑄64 < 2
∗ 10−6
𝑄26 < 2 ∗ 10−6 𝑄36 = 0
−6
𝑄45 = 0
𝑄54 = 0 𝑄54 < 2
∗ 10−6
∗ 10−6
𝑄36 < 2 ∗ 10−6 𝑄46 = 0 𝑄46 < 2 ∗ 10−6
𝑄55 = 0.46
𝑄56 = 0
𝑄55 = 0.46
𝑄56 < 2 ∗ 10−6
𝑄65 = 0
𝑄66 = 0.49
𝑄65 < 2
∗ 10−6
𝑄66 = 0.48
Table III: stiffness matrix for continuous oriented composite. Black: Halpi-Tsai model and stiffness matrix rotaion according to lamination theory. Red: simulation results 𝑄̅11 = 2.3
𝑄̅12 = 0.83
𝑄̅13 = 0.83
𝑄̅14 = 0
𝑄̅15 = 0
𝑄̅16 = 0
𝑄̅11 = 2.2
𝑄̅12 = 0.84
𝑄̅13 = 0.84
𝑄̅14 < 3 ∗ 10−6
𝑄̅15 < 3 ∗ 10−6
𝑄̅16 < 3 ∗ 10−6
𝑄̅21 = 0.83
𝑄̅22 = 1.8
𝑄̅23 = 0.85
𝑄̅24 = 0
𝑄̅25 = 0
𝑄̅26 = 0.0083
𝑄̅21 = 0.84
𝑄̅22 = 1.7
𝑄̅23 = 0.84
𝑄̅24 < 3 ∗ 10−6
𝑄̅25 < 3 ∗ 10−6
𝑄̅26 = 0.0051
𝑄̅31 = 0.83
𝑄̅32 = 0.85
𝑄̅33 = 1.8
𝑄̅34 = 0
𝑄̅35 = 0
𝑄̅36 = −0.09
𝑄̅31 = 0.82
𝑄̅32 = 0.83
𝑄̅33 = 1.7
𝑄̅41 = 0
𝑄̅42 = 0
𝑄̅43 = 0
𝑄̅41 < 3 ∗
10−6
𝑄̅42 < 3 ∗
10−6
𝑄̅43 < 3 ∗
10−6
𝑄̅34 < 3 ∗
10−6
𝑄̅35 < 3 ∗
10−6
𝑄̅44 = 0.48
𝑄̅45 = −6 ∗
10−3
𝑄̅44 = 0.46
𝑄̅45 = −3 ∗
10−3
𝑄̅36 = −0.07 𝑄̅46 = 0 𝑄̅46 < 3 ∗ 10−6
𝑄̅51 = 0
𝑄̅52 = 0
𝑄̅53 = 0
𝑄̅54 = −6 ∗ 10−3
𝑄̅55 = 0.46
𝑄̅56 = 0
𝑄̅51 < 3 ∗ 10−6
𝑄̅52 < 3 ∗ 10−6
𝑄̅53 < 3 ∗ 10−6
𝑄̅54 = −3 ∗ 10−3
𝑄̅55 = 0.45
𝑄̅56 < 3 ∗ 10−6
𝑄̅61 = 0
𝑄̅62 = 0.0083
𝑄̅63 = −0.09
𝑄̅64 = 0
𝑄̅65 = 0
𝑄̅66 = 0.49
𝑄̅61 < 3 ∗ 10−6
𝑄̅62 = 0.0048
𝑄̅63 = −0.06
𝑄̅64 < 3 ∗ 10−6
𝑄̅65 < 3 ∗ 10−6
𝑄̅66 = 0.48
Preliminary tests Initially, in order to properly develop the FEM model, some optimization trials were performed. As demonstrated by Chen and Papathanasiou [28], the simulation data are affected by a high scattering due to the stochastic nature of the reinforcing particle distribution. In order to reduce this scattering, the number of random distributions (Nrd) and the number of reinforcing paticles (Nf) were optimized by keeping a constant vf= 0.01 and the aspect ratio equal to 20 for perfectly aligned reinforcing particles. At first, the attention was focused on the minimum number of particles to insert in the domain in order to obtain a stable solution. The approach is similar to those adopted for studying diffusion in nanocomposites; basically, since stress free boundary conditions are applied on the boundaries parallel to the direction of the applied stress, the number of particles must be adequately high in order to neglect the border effects [20]. Results reported in Figure 2 clearly show that, with increasing
number of reinforcing particles in the domain, an increasing trend of the modulus is observed the modulus increases. However, the FEM results tend to an asymptotic value; keeping in mind the necessity to reduce the computation times, the number of particles was chosen, in each simulation, equal to 25. 1.60
E11 (GPa)
1.56 1.52 1.48 1.44 1.40
10
15 20 25 30 number of reinforcing particles
35
Figure 2: effect of the number of reinforcing particles on simulation results Afterwards, the number of random distribution was optimized; as highlighted in Figure 3, a low number of Nrd involves a very high standard deviation of the simulated longitudinal modulus; as the number of random distributions increases, the standard deviation decreases, and also the average value tends to level off, reaching a stable value at Nrd=10. 1.60 1.55
E11 (GPa)
1.50 1.45 1.40 1.35 1.30
0
5 10 number of random distributions
15
Figure 3: effect of the number of random distributions on simulation results Therefore, for each simulation, results reported are obtained as the average of 10 simulations, each one run with 25 reinforcing particles.
Finally, we focused on the effect of the domain shape; in order to explain its effect on the FEM results, two different bidimensional simulation were performed on the domains reported in Figure 4: in the first case, Figure 4 a), the length of the domain is equal to the length of reinforcing particle LF, increased by an arbitrarily low amount of matrix along the x axis, LV (necessary in order to avoid continuous fiber condition). The additional matrix layer is also highlighted in Figure 4 a) in light red. In the second domain, the length of the domain is increased in order to allow for positioning of two reinforcing particles along the x axis. The interparticle distance along the x axis, for symmetry reasons, is equal to 2LV. Also, in Figure 4 b), the interparticle region is highlighted in light green. A plot of the deformation x along the red line of Figure 4 a) and blue and green lines of Figure 4 b) is reported in Figure 5, together with the average value of the deformation and the stress along the x axis. If one single fiber would represent a satisfactory RVE, then, for symmetry reasons, the three plots would be expected to be perfectly overlapped. Instead, in the interparticle region, a higher deformation is observed; in facts, due to stress continuity at the fiber/matrix interface, a higher stress is also found in the same region, compared to the region LV. However, the higher stress found in the interparticle region contributes to reduce the stress, and therefore the deformation, far away from the fibers; for example, a plot along the purple lines of Figure 4 a) and b), reported in Figure 5 b), clearly shows the reduced stresses and deformation, particularly in the area which is highlighted in light yellow in Figure 4 a). As a consequence of this, the modulus increases with increasing number of fibers along the x axis, Nx. Simulations run on domains correspondent to those reported in Figure 4 a) and b), with increasing length and increasing number of fibers, are reported in Figure 6 a), clearly showing that at least 5 fibers are required in the x direction in order to obtain a stable solution. No error bar is reported in Figure 6 a), since the simulations were obtained on a regular array of fibers arranged according to the scheme of Figure 4; therefore, one single simulation was performed for each value of Nx.
Ldom=LF+2LV
LV
a)
2Ldom=2LF+4LV
LV
2LV
b) x Figure 4: 2D simulation domain with one or two fibers along the principal direction
1 fiber along LV
-7
1.8x10
-8
-7
1.6x10
5.0x10
x,avg=76.5 Pa
1.4x10
2 fibers along 2LV
-7
1.2x10
-8
4.0x10
-8
x
x,avg=10.0*10
-7
1.0x10
x,avg=97 Pa
-8
3.5x10
-8
-8
3.0x10
-8
2.5x10
8.0x10
-8
6.0x10
-8
-8
4.0x10
1 fiber 2 fibers
-8
4.5x10
2 fibers along LV
-7
x
-8
x,avg=8.2*10
0.0
0.2
0.4
0.6
0.8
2.0x10
1.0
a)
0.0
0.2
0.4
0.6
0.8
1.0
x/Ldom
x/LV
b)
Figure 5: deformation along the x axis calculated from 2D simulations; a) along the axis of reinforcing fiber; b) far away from the axis of reinforcing fiber In analyzing the effect of the shape of the domain in 3D simulation, it must be initially emphasized that in any case, for oriented nanoparticles, the minimum length of the domain the x direction, Ltotx, must satisfy: eq. 11)
𝐿𝑡𝑜𝑡𝑥 > 𝐴𝑅 ∗ 𝐷 ∗ 𝑐𝑜𝑠𝜃𝑥
Where D is the diameter of the reinforcing particle. On the other hand, the total volume of the domain is given as: eq. 12)
𝜋𝐷 3
𝑁
𝑉𝑡𝑜𝑡 = 𝜙𝑓 ∗ 𝐴𝑅 ∗
4
𝑓
Where F is the volume fraction of reinforcing phase. Under the constraint of eq. 11), the domain will also satisfy eq. 13)
𝑁 ∗𝜋𝐷 2
𝑉
𝐿𝑡𝑜𝑡𝑦 = 𝐿𝑡𝑜𝑡𝑧 < √𝐿 𝑡𝑜𝑡 = √4∗𝜙𝑓
𝑓 ∗𝑐𝑜𝑠𝜃
𝑡𝑜𝑡𝑥
Therefore, in general, the simulation domain will be a squared base parallelepiped with its own minimum aspect ratio, given as: eq. 14)
𝐿
𝐴𝑅𝑚𝑖𝑛,𝑑𝑜𝑚𝑎𝑖𝑛 = 𝐿𝑡𝑜𝑡𝑥 = 2
𝐴𝑅∗(𝜙𝐹 )1/2 ∗(𝑐𝑜𝑠𝜃)3/2
𝑡𝑜𝑡𝑦
√𝑁𝑓 ∗𝜋
On the other hand, if, according to the previous discussion, Nx particles need to be aligned in the x direction, then eq. 11) is modified as: eq. 15)
𝐿𝑡𝑜𝑡𝑥 = 𝑁𝑥 𝐴𝑅 ∗ 𝐷 ∗ 𝑐𝑜𝑠𝜃𝑥
Which substituted into eq. 12) and eq. 13) yields an actual aspect ratio of the domain: 1
eq. 16)
𝐴𝑅𝑑𝑜𝑚𝑎𝑖𝑛 = 2
𝐴𝑅∗(𝜙𝐹 )2 ∗(𝑁𝑥 )3/2 ∗(𝑐𝑜𝑠𝜃)3/2 √𝑁𝑓 ∗𝜋
The results of Figure 6 b) are a plot of the modulus as a function of the ratio between the actual aspect ratio of the domain and the minimum aspect ratio of the domain with randomly dispersed discontinuous reinforcing cylinders in 3D simulations. As clearly observed, the modulus increases with increasing aspect ratio of the domain. Even in this case, the modulus reaches a plateau value for sufficiently high values of 𝐴𝑅𝑑𝑜𝑚𝑎𝑖𝑛 , which is roughly 10 times the minimum aspect ratio. It is worth observing that, according to eq. 16),
𝐴𝑅𝑑𝑜𝑚𝑎𝑖𝑛 𝐴𝑅𝑚𝑖𝑛,𝑑𝑜𝑚𝑎𝑖𝑛
= (𝑁𝑥 )3/2, and therefore a ratio of 10 roughly
corresponds to 𝑁𝑥 = 4.7, which is a value very close to that observed in Figure 6 a) for bi-dimensional simulation. Based on the results of Figure 6 -Figure 6, in each domain the minimum length of the domain along the x axis was taken according to eq. 15). Finally, it is worth observing that, for high values of the aspect ratio of the reinforcing particles and volume fractions, the aspect ratio of the domain can reach very high values (about 200). In order to test the stability of the developed approach towards such an elongated domain, the simulation on unidirectional continuous fiber reinforced composite, was also performed with an aspect ratio of the domain equal to 200. The stability of the solution is confirmed by the fact that the coefficient of regression between the simulated stiffness matrix and the stiffness matrix estimated by Halpin-Tsai, eq. 10), is R2=0.9999.
The results obtained indicate that with an adequately high number of nanoparticles, boundary effects do not significantly influence the FEM results; therefore, the choice of a “free” boundary condition on the faces parallel to the direction of the applied load is justified. This is a distinct advantage of the proposed approach. In most of the reported works, in facts, the boundary condition on the faces parallel to the direction of the applied load is considered to be symmetric, which, however, does not allow for its contraction, yielding a null Poisson modulus. Instead the proposed approach allows for estimating the whole stiffness matrix of the composite. 4.04
7.2
4.02
6.8 6.4 6.0
3.98
E11 (GPa)
E11 (GPa)
4.00
3.96 3.94
5.2 4.8
3.92 3.90
5.6
4.4 1
2
3
Nx
4
5
6
4.0
2
4
6 8 10 ARdomain/ARmin,domain
12
14
Figure 6: effect of domain shape on longitudical modulus; a) for 2D simualtions; b) for 3 D simulations Simulation in DGFRC In Figure 7, the longitudinal modulus along the principal direction 1 (E11) is reported for unidirectional discontinuous glass fibers. Also, model prediction according to the Halpin-Tsai theory, obtained with 𝜉 = 2𝐴𝑅 [29], clearly shows an excellent good agreement with the FEM simulation data, as highlighted by the very high coefficient of regression between FEM and HT, R2=0.92. On the other hand, as reported in Figure 7 a), the Tandon-Weng (TW) model was also applied. This approach, which represents a simplification of the Mori-Tanaka method, significantly overestimates the longitudinal tensile modulus. For the calculation of the TW model, the Eshleby’s tensor coefficients were calculated according to the procedure for prolate ellipsoid reported in [30]. From these, the elastic coefficients of the composite were obtained according to the formula reported in [31]. The principal Poisson ratio was estimated by the non-iterative approach proposed in [32]. However, the TW model shows a very good agreement with the FEM simulation results at low values of F, but significantly overpredicts the FEM simulation results at high values of F. The other elastic constants, E22, G12 and 12 from FEM modeling are reported in Figure 7 b), c) and d), together with Halpin-Tsai and Tandon-Weng models prediction. The HT values were initially obtained with =2 for E22, =1 for G12 and =2AR for 12.
The transversal modulus E22 from TW is in quite good agreement with FEM simulation data (R2=0.85), whereas the HT model, not shown in Figure 7 b), significantly underpredicts the E22 from FEM (R2=0.27). On the other hand, when using as a fitting parameter, the best fitting of the FEM simulation data was found for =4.4, (R2=0.94) as shown in Figure 7 b). As reported in Figure 7 c), the HT model (=1) significantly underpredicts the G12 from FEM simulation (R2=0.49). The TW model, not shown in Figure 7 c), provides lower estimate than the HT model (R2=-0.44). In this case, the best fitting with the HT model was found using =2.5 (R2=0.99). Finally, as reported Figure 7 d), the 12 from FEM is much higher than the corresponding HT (=2AR) and TW models. In both cases the poor agreement between analytical models and simulation is highlighted by the very low R2=-75. However, in this case, using as a fitting parameter did not allow to improve the accuracy of model prediction. FEM simulation Halpin-Tsai =2AR Tandon-Weng
4
1.18 1.14
F=0.019
2
F=0.0064 F=0.0032
1.12
F=0.019
1.10 1.08 1.06 1.04
F=0.0064
1.02
F=0.0032
1.00 20
40
60
80
100 AR
120
140
20
160
a)
40
60
80
100 AR
120
140
160
b) FEM simulation Halpin-Tsai =1 Halpin-Tsai =2.5
0.43
0.330
F=0.032
F=0.0032 F=0.0064
0.329 F=0.019
0.40
12
G12 (GPa)
0.41
0.39 F=0.0064
0.38 0.37
FEM simulation Halpin-Tsai =2AR Tandon-Weng
0.331
0.42
c)
F=0.032
1.16
F=0.032
E22 (GPa)
E11 (GPa)
3
1
FEM simulation Halpin-Tsai =4.4 Tandon-Weng
1.20
40
60
80
100 AR
120
140
F=0.019
0.327
F=0.0032
20
0.328
F=0.032
0.326
160
20
40
60
80
100 AR
120
140
160
d) Figure 7: elastic constants for unidirectional DGFRC; a) longitudinal modulus; b) transverse modulus; c) longitudinal shear modulus; d) principal Poisson ratio
Afterward, the effect of the orientation angle was taken into account. Therefore, the main orientation angle (measured with respect to the x axis, and therefore corresponding to a rotation around the y or the z axis) was changed between 0 and 54 °. However, the simple rotation along the z axis, would
lead to a composite with Eyy>Ezz, as clearly observed in Figure 8 a); though this is a condition which could be encountered in polymer composites or nanocomposites, in many cases the existence of a single direction of flow during processing (extrusion, injection molding) causes orientation of reinforcing particles along one single direction (the x axis in our example), but quasi-isotropic properties in the other two directions (y and z in the present case). In order to attain this condition, the azimuthal angle yz (measured with respect to the y or the z axis, and therefore corresponding to a rotation around the x axis, as highlighted in Figure 8), was randomly varied between 0 and 180°. The resulting domain reported in Figure 8 b) guarantees quasi-isotropic behavior in the y-z plane. x x z
a)
yz
y
b)
Figure 8: sketch of the rotation angles in oriented composites The results obtained for the different elastic constants are reported in Figure 9 a)-d). In this case, no comparison with closed form analytical model is possible. In facts, both Halpi-Tsai and TandonWeng are not suitable to accurately predict the evolution of the elastic constants in unidirectional composites, which is the basis for calculating the elastic constants in oriented composites; however, for each F and AR, the elastic constants of oriented composites were calculated from the elastic constant of unidirectional composites, obtained from FEM analysis, followed by rotating the compliance matrix by an angle x around the y axis and further rotating by an angle yz=45° around the x axis, which corresponds to a random orientation also in the y-z plane [33]. The azimuthal angle does not affect the longitudinal modulus Exx, but it affects the other elastic constants reported in Figure 9. In all cases, a very good agreement is observed between the FEM simulated constants and constants obtained from properties calculated in the principal direction followed by matrix rotation.
In particular, for Exx, a monotonic behavior with increasing orientation angle is found; for x=54 °, the value of Exx is lower than the corresponding value of E22; this is due to the fact that for all the simulated GFRP, the value of G12 is lower than
𝐸11 [26]. 𝐸 2( 11 +𝜈12 ) 𝐸22
However, the FEM simulation is able
to capture this behavior of the GFRP. For the same reason, in the range between 0 ° and 45 °, the value of Eyy decreases with increasing orientation angle. Even in this case, the FEM simulation is able to highlight the decrease of Eyy. The FEM model is also able to predict the evolution of Gxy, and of xy with their corresponding maximum. FEM simulation prediction from elastic constants along the principal axes and formula for rotation
1.8
AR=20
F=0.0064
AR=80
1.4
F=0.0064
AR=20 =0.0032 F
1.2
F=0.0019
1.08
Eyy (GPa)
Exx (GPa)
1.10
F=0.0019
1.6
FEM simulation prediction from elastic constants along the principal axes and formula for rotation
AR=40
F=0.0032
1.06
AR=20
F=0.0064
AR=80
F=0.0064
AR=20
1.04
F=0.00032
AR=40
F=0.0032
AR=20
1.02
AR=20 1.0 0
10
20 30 40 orientation angle (°)
50
a)
0
10
20 30 40 orientation angle (°)
50
60
b) FEM simulation prediction from elastic constants along the principal axes and formula for rotation
0.350
F=0.019
AR=20
0.435 0.420
FEM simulation prediction from elastic constants along the principal axes and formula for rotation
0.345 F=0.0064
F=0.0064
AR=80
AR=20 F=0.0032
0.405
AR=40
0.390
F=0.0032
0
10
20 30 40 orientation angle (°)
50
c)
F=0.0064
AR=20 0.335 F=0.0032 =0.0032 F
0.330 0.325
60
AR=80
AR=20
AR=40
AR=20
0.375
F=0.0064
F=0.019
0.340 xy
0.450
Gxy (GPa)
1.00
60
0
10
AR=20
20 30 40 orientation angle (°)
50
60
d)
Figure 9: elastic constants for GFRP at different orientation angles. a) longitudinal modulus; b) transverse modulus; c) longitudinal shear modulus; d) principal Poisson ratio. black F=0.0032, AR=20; red F=0.0032, AR=40; blue F=0.0064, AR=20; magenta F=0.0064, AR=80; olive F=0.019, AR=20 Simulation in CNT nanocomposites The longitudinal modulus for fully oriented nanoparticles, together with Halpin-Tsai prediction obtained with 𝜉 = 2𝐴𝑅, is reported in Figure 10 a). As in the case of glass fibers, the HT theory is in reasonable good agreement with the FEM simulation data. For CNTRNC, no comparison is provided
with Tandon-Weng, since closed form solutions for this approach are only available for isotropic reinforcement [34]. The results for the other elastic constant, reported in Figure 10 b), c) and d), substantially confirm the behavior already discussed for glass fibers: the HT model underpredicts both E22 (=2) and G12 (=1) from FEM; 12 from FEM is also much higher than the corresponding HT (=2AR). For E22, when considering as a fitting parameter, the best fit is obtained with =5.4; for G12, the best fit is obtained with =2.2. FEM simulation Halpin-Tsai =2AR
FEM simulation Halpin-Tsai =2 Halpin-Tsai =5.5
1.25 F=0.032
1.20
F=0.032
1.15
F=0.019
F=0.019
4
F=0.0064
3
E22 (GPa)
E11 (GPa)
10 9 8 7 6 5
F=0.0032
2
F=0.0064
1.10
F=0.0032
1.05 1
20
40
60
80
100 AR
120
140
a)
40
60
80
100 AR
120
140
0.330
F=0.032
160
0.405 0.400
F=0.0032
0.329
FEM simulation Halpin-Tsai =1 Halpin-Tsai =2.2
0.410
F=0.0064
FEM simulation Halpin-Tsai =2AR
0.328 F=0.019 12
G12 (GPa)
20
b) 0.415
0.395
0.327
F=0.019
0.326
0.390 0.385
F=0.0064
0.380
F=0.0032
20
c)
1.00
160
40
60
80
100 AR
120
140
0.325 0.324
160
F=0.032
20
40
60
80
100 AR
120
140
160
d)
Figure 10: elastic constants for unidirectional CNTRNC; a) transverse modulus; b) longitudinal shear modulus; c) principal Poisson ratio Also in the case of CNT nanocomposites, the effect of the orientation angle was taken into account. Following the same procedure used for GFRP, the results obtained for the different elastic constants are reported in Figure 11 a)-d). For comparison purposes, the elastic constant were also calculated from the elastic constant of unidirectional composite, followed by rotating the compliance matrix by an angle x around the y axis and further rotating by an angle yz=45 ° around the x axis. As in the case of GFRP, a very good agreement is observed between the FEM simulated constants and constants obtained from properties calculated in the principal direction followed by matrix rotation.
FEM simulation prediction from elastic constants along the principal axes and formula for rotation
3.6 3.2
F=0.019
1.14
AR=40
1.12
F=0.0064
F=0.019
1.10
AR=80
AR=40
2.8
F=0.0064
AR=80
2.4
F=0.0064
AR=40
2.0
Eyy (GPa)
Exx (GPa)
FEM simulation prediction from elastic constants along the principal axes and formula for rotation
1.16
F=0.0032
AR=40
1.6
F=0.0032
F=0.0064
1.06
AR=40
1.04
AR=20
1.2
1.08 F=0.0032
F=0.0032
AR=40
AR=20
1.02 0
10
20 30 40 orientation angle (°)
50
0
60
a)
10
20 30 40 orientation angle (°)
50
60
b)
0.52
FEM simulation prediction from elastic constants along the principal axes and formula for rotation F=0.019
0.50
AR=40
Gxy (GPa)
0.48 0.46
0.365
F=0.0064
0.360
AR=80
0.355
F=0.0064
0.350
AR=40
0.44
0.370
xy
0.54
F=0.0032
0.42
F=0.0032
0.38
AR=20 0
c)
10
20 30 40 orientation angle (°)
50
F=0.019
AR=40 F=0.0064
F=0.0064
AR=40
AR=80
0.340
AR=40
0.40
0.345
FEM simulation prediction from elastic constants along the principal axes and formula for rotation
0.335
F=0.0032
0.330
AR=40
0.325
60
0
10
20 30 40 orientation angle(°)
F=0.0032
AR=20 50
60
d)
Figure 11: elastic constants for CNTRNC at different orientation angles: a) longitudinal modulus; b) transverse modulus; c) longitudinal shear modulus; d) principal Poisson ratio. black F=0.005, AR=20; red F=0.005, AR=40; blue F=0.01, AR=40; magenta F=0.01, AR=80; olive F=0.03, AR=40 Conclusions In this work, a geometrical FEM model previously derived to study diffusion in polymer nanocomposites has been properly modified in order to study the effect of different morphological features on the stiffness evolution in short fiber reinforced composites and nanocomposites. Compared to similar works reported in literature, the developed model can account for different morphological features, including processing induced fiber orientation. The accuracy of the developed approach has been tested by comparing the simulation results with closed form solutions for continuous fiber reinforced composites, also accounting for orientation of the fibers. After, the FEM model has been used in order to calculated the stiffness matrix of discontinuous glass fiber reinforced composite (DGFRC) as well as carbon nanotubes reinforced nanocomposites (CNTRNC).
For DGFRC, comparison with closed form analytical models (Halpin-Tsai and Tandon-Weng models) showed that HT is able to predict the evolution of elastic constants, provided that the geometric factor is sued as a fitting parameter. On the other hand, the TW model does not provide an adequate estimate for the elastic constants of the DGFRC. For CNTRNC, comparison with Tandon-Weng model was not performed, since closed form solutions were only found for isotropic reinforcement; however, comparison with the Halpin-Tsai model substantially confirmed the results found for glass fibers. Finally, the developed FEM model was used for calculating the stiffness at different orientation of the reinforcement. In this case, results were not compared with classical micromechanical approaches, since these latter were proven to be unsuitable to predict the evolution of longitudinal modulus even for unidirectional reinforcement. Instead, for each value of volume fraction and aspect ratio of the reinforcement, the longitudinal, transversal, shear and Poisson modulus were calculated from the corresponding values obtained for uni-directional composites from FEM, followed by stiffness matrix rotation according to classical lamination theory. The elastic properties from FEM analysis are in excellent good agreement with the corresponding values calculated according to the classical lamination theory, indicating the accuracy of the developed FEM model for estimation of stiffness of composites and nanocomposites, also accounting for reinforcement orientation.
References 1 Endo M, Hayashi T, Kim YA, Terrones M, Dresselhaus MS. Applications of carbon nanotubes in the twenty-first century. Philos Trans R Soc London 2004; A 362: 2223–2238 2 Moniruzzaman M, Winey KI. Polymer nanocomposites containing carbon nanotubes. Macromolecules 2006; 39: 5194–5205. 3 Wernik JM, Meguid SA. Multiscale micromechanical modeling of the constitutive response of carbon nanotube-reinforced structural adhesives. Int J Solids Struct 2014; 51: 2575–2589 4 Wang HW, Zhou HW, Peng RD, Mishnaevsky L Jr. Nanoreinforced polymer composites: 3D FEM modeling with effective interface concept. Compos Sci Technol 2011; 71: 980–988 5 Rafiee R, Firouzbakht V. Multi-scale modeling of carbon nanotube reinforced polymers using irregular tessellation technique. Mech Mater 2014; 78: 74–84 6 Schiebold M, Schmidt H, Mehner J. A Finite Element Approach for Modeling and Simulation of CNT/Polymer Composites. Phys Status Solidi A 2019; 1800952 7 Kumar P, Srinivas J. Numerical evaluation of effective elastic properties of CNT-reinforced polymers for interphase effects. Comput Mater Sci 2014; 88: 139–144
8 Tserpes KI, Chanteli A. Parametric numerical evaluation of the effective elastic properties of carbon nanotube-reinforced polymers. Compos Struct 2013; 99: 366–374 9 Garcia-Macias E, Guzman CF, Saavedra Flores EI. Castro-Triguero R. Multiscale modeling of the elastic moduli of CNT-reinforced polymers and fitting of efficiency parameters for the use of the extended rule-of-mixtures. Compos Part B-Eng 2019; 159: 114–131 10 Wang JF, Liew KM. On the study of elastic properties of CNT-reinforced composites based on element-free MLS method with nanoscale cylindrical representative volume element. Compos Struct 2015; 124: 1–9 11 Joshi P, Upadhyay SH. Evaluation of elastic properties of multi walled carbon nanotube reinforced composite. Comput Mater Sci 2014; 81: 332–338 12 Shahzad Zuberi MJ, Esat V. Investigating the mechanical properties of single walled carbon nanotube reinforced epoxy composite through finite element modelling. Compos Part B-Eng 2015; 71: 1–9 13 Pantano A, Mantione P. A numerical-analytical model for the characterization of composites reinforced by carbon nanotubes. Appl Phys A 2010; 99: 895–902 14 Gupta AK, Harsh SP. Analysis of mechanical properties of carbon nanotube reinforced polymer composites using multi-scale finite element modeling approach. Compos Part B-Eng 2016; 95: 172178 15 Nie J, Jia Y, Qu P, Shi Q. Carbon Nanotube/Carbon Fiber Multiscale Composite: Influence of Interfacial Strength on Mechanical Properties. J Inorg Organomet Polym 2011; 21: 937–940 16 Yengejeh SI, Kazemi SA, Öchsner A. Carbon nanotubes as reinforcement in composites: A review of the analytical, numerical and experimental approaches. Comput Mater Sci 2017; 136: 85–101 17 Lusti HR, Gusev AA. Finite element predictions for the thermoelastic properties of nanotube reinforced polymers. Model Simul Mater Sci Eng 2004; 12: S107-S119 18 Chen X, Alian AR, Meguid SA. Modeling of CNT-reinforced nanocomposite with complex morphologies using modified embedded finite element technique. Compos Struct 2019; 227: 111329 19 Babu KP, Mohite PM, Upadhyay CS. Development of an RVE and its stiffness predictions based on mathematical homogenization theory for short fibre composites. Int J Solids Struct 2018; 130– 131: 80–104 20 Greco A, Maffezzoli A. Two-dimensional and three-dimensional simulation of diffusion in nanocomposite with arbitrarily oriented lamellae. J Membr Sci 2013; 442: 238-244 21 Greco A. Numerical simulation and mathematical modeling of 2D multi-scale diffusion in lamellar nanocomposite. Comput Mater Sci 2014; 90: 203–209
22 Greco A. Simulation and modeling of diffusion in oriented lamellar nanocomposites. Comput Mater Sci 2014; 83: 164–170 23 Greco A, Maffezzoli A. Finite Element Modeling of Multiscale Diffusion in Intercalated Nanocomposites. J Nanomater 2015; 482698 24 Greco A, Maffezzoli A. Finite element simulation and analytical modeling of 3D multi scale diffusion in nanocomposites with permeable stacks. Model Simul Mater Sci Eng 2016; 24 (1): 015003 25 Greco A, Esposito Corcione C, Maffezzoli A. Effect of multi scale diffusion on the permeability behaviour of intercalated nanocomposites. J Membr Sci 2016; 505: 92–99 26 Jones RM, Mechanics of Composite Materials, 2nd ed, 1999, Taylor & Francis, Philadelphia, US 27 Tornabene F, Bacciocchi M. Fantuzzi N, Reddy JN. Multiscale Approach for Three-Phase CNT/Polymer/Fiber Laminated Nanocomposite Structures. Polym Compos 2019; 40: E102-E126. 28 Chen X, Papathanasiou TD. Micro-scale modeling of axial flow through unidirectional disordered fiber arrays, Compos Sci Technol 2007; 67: 1286-1293 29 Halpin JC, Kardos JL. The Halpin-Tsai equations: a review. Polym Eng Sci 1976; 16(5): 344-352 30 David EC, Zimmerman RW, Compressibility and shear compliance of spheroidal pores: Exact derivation via the Eshelby tensor, and asymptotic expressions in limiting cases. Int J Solids Struct 2011; 48: 680–686 31 Zare Y, Development of simplified Tandon-Weng solutions of Mori-Tanaka theory for Young’s modulus of polymer nanocomposites considering the interphase. J Appl Polym Sci 2016: 43816 32 Tucker III CL, Liang E. Stiffness predictions for unidirectional short-fiber composites: Review and evaluation, Compos Sci Technol 1999; 59 (5): 655-671 33 Mortazavian S, Fatemi A. Effects of fiber orientation and anisotropy on tensile strength and elastic modulus of short fiber reinforced polymer composites. Compos Part B-Eng 2015; 72: 116–129 34 Tandon GP, Weng GJ. The effect of aspect ratio of inclusions on the elastic properties of unidirectionally aligned composites, Polym Compos 1984; 5(4): 327-333