polymer nanocomposites, using complex and time-dependent homogenizations

Computational Materials Science 142 (2018) 395–409

Contents lists available at ScienceDirect

Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci

Multiscale modeling of the viscoelastic properties of CNT/polymer nanocomposites, using complex and time-dependent homogenizations Ali R. Shajari a, Rahmatollah Ghajar a,⇑, Mahmood M. Shokrieh b a

Mechanical Properties Research Lab (MPRL), Faculty of Mechanical Engineering, K.N. Toosi University of Technology, No. 17, Pardis St., Mollasadra Ave., Vanak Square, Tehran, Iran Composites Research Laboratory, Center of Excellence in Experimental Solid Mechanics and Dynamics, School of Mechanical Engineering, Iran University of Science and Technology, Narmak, Tehran 16846-13114, Iran b

a r t i c l e

i n f o

Article history: Received 16 May 2017 Received in revised form 10 September 2017 Accepted 1 October 2017

Keywords: CNT reinforced polymer composite Viscoelastic properties Multiscale model Equivalent inclusion Equivalent fiber

a b s t r a c t In the present research, a combined molecular structural/micromechanics multiscale model was developed to predict the viscoelastic properties of carbon nanotube (CNT) reinforced polymer composites. At the nanoscale, transversely isotropic time-dependent properties of an embedded CNT in the matrix were obtained using a finite element (FE) based molecular structural model in which the beam, solid continuum and nonlinear spring elements were used for carbon bonds, surrounding matrix and van der Waals interactions, respectively. To homogenize the heterogeneous nature of nanocomposites at microscale, Mori-Tanaka method was applied and the corresponding overall viscoelastic properties were obtained. To this end, nano-structural equivalent structures (equivalent fiber (EF) and equivalent inclusion (EI)) were developed to consider the nano-structural interactions. To apply micromechanics formulations, using the Fourier transform viscoelastic properties were transformed to frequency domain in algebraic forms. Additionally, an iterative-based solution for the Mori-Tanaka method was introduced and performed to use the time-dependent properties in the micromechanics model. The overall time dependent, storage and loss moduli of nanocomposites obtained by different methods were presented. A comparison of the creep response of the two models with available experimental data reveals that the EI model can predict more realistic properties compared to those obtained by the EF model. Ó 2017 Elsevier B.V. All rights reserved.

1. Introduction Applications of the CNT reinforced polymer composites have been spread in several industries such as automotive and aerospace in recent years. To the end of mechanical design of the nanocomposite structures, modeling and simulation of the properties of these materials and structures are needed. Various modeling techniques have been developed for prediction of material properties e.g. quantum, atomistic, continuum and multiscale modeling. In some studies of continuum modeling approaches, [1–4], the CNTs were assumed to be long cylindrical inclusions and micromechanical methods with considering perfect or imperfect bonding conditions were directly used. Due to the importance of interfacial bonding, some researches [5–7], considered interfacial sliding at the interface of the CNT and matrix and the modified forms of the micromechanical formulations were used. In the mentioned investigations, due to the continuum modeling of the CNTs, the discrete nature of the nanostructures was not taken into account. ⇑ Corresponding author. E-mail address: [email protected] (R. Ghajar). https://doi.org/10.1016/j.commatsci.2017.10.006 0927-0256/Ó 2017 Elsevier B.V. All rights reserved.

Considering discrete nature of the CNT, molecular dynamic (MD) simulations for an embedded CNT in the matrix were implemented in several studies, e.g. [8–10]. These methods of simulation are not applicable for the structures contain large number of atoms for long time. So, in some investigations, it was tried to correlate atomistic simulation with continuum modeling of CNT reinforced nanocomposites. In the pioneering work of Odegard et al. [11], the CNT and interfacial region were modeled by an equivalent fiber. The transversely isotropic properties of the fiber were obtained by equating the results of molecular dynamic simulations and micromechanics for an embedded CNT in the matrix. Shokrieh and Rafiee [12] presented equivalent long fiber properties for the CNT and surrounding interface by developing a FE-based molecular structural multiscale model. They used beam, spring and solid continuum elements for Carbon bonds, van de Waals interfacial interactions and polymer matrix, respectively [12–14]. The nanocomposite materials naturally contain constituents in various scales, such that multiscale modeling methods have been widely implemented to account different scales. The results of the multiscale modeling methods are more efficient than atomistic simulation approaches and more accurate than those of continuum

396

A.R. Shajari et al. / Computational Materials Science 142 (2018) 395–409

modeling [15]. Savvas et al. [16] used an equivalent beam structure for the CNT, in which the corresponding properties were evaluated and considered using the results of the MD simulations. At the macro scale, they applied a FE model for prediction of the properties and investigation of the effects of interfacial shear strength on the properties. As separate studies, the solid continuum elements were used to model the CNT, matrix and interphase, respectively and the elastic properties of the composites were evaluated [17– 19]. To study the nonlinear elastic response of the nanocomposites reinforced with CNTs, Wernik and Meguid [20] developed a continuum multiscale representative volume element (RVE) based on the atomistic structure of CNT, interface and matrix. In their continuum RVE, the beam, brick and spring elements were used for Carbon bonds, surrounding matrix, and interfacial interactions, respectively. They applied a nonlinear form of micromechanical formulations to predict overall elastic properties. To consider interfacial interactions in the framework of multiscale modeling, Liu et al. [21] derived cohesive law for van der Waals interactions of CNT/polymer composites. They evaluated tensile cohesive strength and cohesive energy in terms of volume densities of polymer and nanotube. For the viscoelastic properties, the atomistic-based simulation methods have not been used due to the long term of creep and relaxation responses of polymers. Therefore, the micromechanical approaches have been widely applied ([22–26]). In these studies, the constitutive time-dependent equations of the viscoelastic materials were converted to an algebraic form, applying the Laplace or Fourier transforms. The transformed form of the constitutive equations can be used in formulations of micromechanical methods. The interface conditions can be considered as either perfect ([23,24]) or imperfect ([25,26]) bonds for viscoelastic properties. To consider the non-bonded interface for viscoelastic properties, recently a viscoelastic time-dependent formulation for sliding and normal separation of the CNT in the matrix was developed by the present authors [26]. Using a multiscale model, the corresponding timedependent interfacial parameters were estimated and implemented to predict overall properties. Applying time dependent form of the constitutive equations of viscoelastic materials in computational micromechanics were carried out by El Kouri et al. [27]. They used dynamic Green’s functions, integral equations and Volterra product to estimate properties through Mori-Tanaka approach. It is worth mentioning that the researches dealing with modeling of the viscoelastic properties of the CNT/polymer composites, the CNTs were considered as elastic reinforcement and only time-dependent structural properties of the polymer matrix were taken into account. In the present study, two multiscale modeling approaches were considered. In the first approach, the CNT and the surrounding interface were considered as an equivalent fiber and its properties were evaluated using the properties of the FE-based molecular structural mechanics modeling of a CNT embedded in the matrix. To predict overall properties, the micromechanical method of Mori-Tanaka was used with the assumption of randomly oriented fibers. In the second method, developed in the present research, the embedded CNT in the matrix was considered as an equivalent inclusion with same geometries and viscoelastic properties. Then, the equivalent inclusions were assumed to be dispersed throughout the matrix and the Mori-Tanaka method was utilized. The effect of the CNT length on the transversely isotropic viscoelastic properties was investigated and it was found that the upper bonds of properties belong to the case of the long CNT model. On the other hand, a comparison between the predictions of the two models with experimental data shows that the creep properties of the models were underestimated compared to experimental data. However, using the second model leads to more satisfactory predictions. Moreover, a micromechanical modeling for relaxation properties of composites using the time integral constitutive

equations without considering a transformed form of these equations was proposed. The transversely time dependent properties of composites obtained by this modeling method were compared and verified by those obtained by the conventional complex formulations of micromechanics. 2. Modeling strategy Multiscale models of the present study were introduced as: Model I) Equivalent fiber (EF) and Model II) Equivalent inclusion (EI). In model I, at first, a representative volume element (RVE) containing a CNT in the matrix with the same volume fraction of the reinforcement was modeled based on the molecular structural mechanics and corresponding viscoelastic properties were obtained. Then, the CNT and surrounding interface were equated by an EF and its transversely isotropic properties were estimated, using the rule of mixture (ROM) method. The obtained properties of the EFs were applied in Mori-Tanaka approach to predict the overall viscoelastic properties of the composites. Fig. 1 shows the steps of the model I. The random distribution of the CNTs in the resin was taken into account by micromechanics modeling. Hence, the properties of EFs were inserted in micromechanical formulations and the overall viscoelastic properties of the composite were estimated with the assumption of randomly distributed CNTs. In model II, three following steps were considered: At the first step, the embedded CNT in the matrix were modeled using the molecular structural mechanics. Next, the whole of the mentioned model were considered as an equivalent inclusion (EI) and the corresponding viscoelastic properties were obtained. At the last step, EIs with different orientations were considered to be dispersed in the matrix and the Mori-Tanaka method was used to predict overall properties. The schematic of the above steps are shown in Fig. 2. It is worth mentioning that the EI’s size can be expressed by two volume fractions: volume fractions of CNTs in the EIs (cCNT ) and the volume fraction of EIs in the matrix (cEI ). Therefore the overall volume fraction of the CNTs in the composite (cr ) can be stated as:

cr ¼ cCNT  cEI

ð1Þ

The main difference between two modeling strategies is related to considering time-dependent properties of the reinforcements. The properties of the reinforcement at the microscale homogenization were considered elastic and viscoelastic in EF and EI models, respectively. In addition, EF model was applied in prior studies (e.g., [11,12]) for modeling of elastic properties of the CNTreinforced composites. In the present investigation, both EF and EI models were utilized in order to account viscoelastic properties of the reinforcement containing CNT, interface and surrounding polymer. 3. Molecular structural mechanics According to the time order of creep and relaxation responses of the polymers about 103 – 106 sec [28], the atomistic-based simulations cannot be used for the simulation of the structures containing large number of atoms in a large range of time. Hence, the FE-based model which was proposed by Shokrieh and Rafiee [12–14] was implemented here for viscoelastic properties using ABAQUS [29] software. 3.1. CNT modeling Based on the model developed by Li and Chou [30], a spaceframe structure containing beam elements as Carbon bonds and joint nodes as Carbon atoms, was used to model the CNT structure. In this model, the strain energies of the structure were correlated

397

A.R. Shajari et al. / Computational Materials Science 142 (2018) 395–409

Fig. 1. Schematic of model I steps: (a) CNT Polymer nanocomposite with volume fraction c r , (b) molecular structural mechanics model for an embedded CNT in the matrix with the same volume fraction of (a), (c) Obtaining properties of equivalent fibers which represent the CNT and interface, d) Model of the nanocomposite containing EFs.

Fig. 2. Schematic of model II steps: (a) CNT/Polymer nanocomposite with volume fraction cr , (b) molecular structural mechanics model for CNT and matrix with volume fraction cCNT , (c) Equivalent inclusion with same volume, geometry and viscoelastic properties as (b), (d) Model of the nanocomposite containing viscoelastic EIs with volume fraction c EI .

with the interatomic molecular potential energies to obtain structural properties. The energies related to stretching, bending and torsion of the bonds can be expressed as:

1 U r ¼ kr ðDrÞ2 ; 2

1 U h ¼ kh ðDhÞ2 ; 2

1 U s ¼ ks ðDsÞ2 2

ð2Þ

where Dr, Dh and D/ denote bond stretching increment, bond angle change and angle change of bond twisting, respectively and kr , kh and ks are the bond stretching force constant, bond angle bending force constant and torsional resistance, respectively. The relations between sectional stiffness parameters and constants of force fields can be written as:

kr ¼

EA ; a

kh ¼

EI ; a

ks ¼

GJ a

ð3Þ

where EA, EI, GJ are sectional rigidities under tension, bending and torsion of the beam elements, respectively. Moreover, a is the Carbon bond length, chosen as 0:142 nm. According to the AMBER force field, the constants of force fields are given by [30]:

kr ¼ 6:52  107 Nnm1 kh ¼ 8:76  1010 N nm1 rad 10

ks ¼ 2:78  10

1

N nm

2

rad

2

ð4Þ

398

A.R. Shajari et al. / Computational Materials Science 142 (2018) 395–409

Considering a circular cross section with a diameter d for Carbon bonds and substituting Eq. (4) in Eq. (3), the elastic properties of the bonds are obtained as:

E ¼ 5:780 TPa; G ¼ 0:871 TPa; d ¼ 0:146 nm

ð5Þ

These values are inserted in the software for the properties of the beam elements.

p1 ¼ 2:9164  103 sec; q0 ¼ 1:8175  109 Pa; q1 ¼ 6:9308  1012 Pa:sec

ð10Þ

The result of the viscoelastic model with above constants and the experimental data are illustrated in Fig. 3. In this figure, maximum deviations between the model and experimental data are up to 1%.

3.2. Matrix modeling 3.3. Interfacial interactions Generally, CNT reinforced composites have volume fractions lower than 5%, so that at the molecular scale, it can be stated that the volume of polymer is much higher than that of CNTs [12–14]. Hence, modeling of the composites at a molecular scale is a difficult task due to the large amount of polymer chains. Therefore, in the model of the present study, continuum solid elements were used for simulation of the surrounding polymer. In the ABAQUS software, the second order 3-D solid continuum brick elements with 20 nodes (C3D20R) were utilized. The surrounding matrix was modeled as a hollow cylinder; and the inner and outer radii were determined according to radius and the volume fraction of the CNT (cCNT ), respectively [26]. The size of the element was selected as small as the length of hexagonal rings of CNTs in axial direction and the depth of the region was calculated in a manner to reflect the proper volume fraction for an embedded CNT [12–14]. The constitutive equations of the viscoelastic material can be expressed as:

rm ðtÞ ¼

Z 0

t

Cm ðt  sÞ :

dem ðsÞ ds ds

ð6Þ

The matrix was assumed to be isotropic, so the relaxation tensor is:

2

C 11

6C 6 12 6 6 C 13 Cm ðtÞ ¼ 6 6 0 6 6 4 0 0

C 12

C 13

0

0

C 22

C 23

0

0

C 23

C 33

0

0

0

0

C 44

0

0

0

0

C 55

0

0

0

0

C 11 ¼ C 22 ¼ C 33 ¼

Em ðtÞð1  mm Þ ð1 þ mm Þð1  2mm Þ

C 12 ¼ C 23 ¼ C 13 ¼

Em ðtÞmm ð1 þ mm Þð1  2mm Þ

C 44 ¼ C 55 ¼ C 66 ¼

Em ðtÞ 2ð1 þ mm Þ

0

3

0 7 7 7 0 7 7 0 7 7 7 0 5

ð7Þ

C 66

The CNT-Polymer interactions, are naturally occurred via the vdW forces which are modeled by the Lennard-Jones ‘‘6–12” potential as follows [33]:

 12  6  r r V LJ ¼ 4e  r r

ð11Þ

where e ¼ 0:0566 kcal=mol and r ¼ 0:34 nm, are Lennard-Jones parameters and r is the atomic distance. According to the Lennard-Jones potential, the vdW force can be stated as:

F vdW ¼ rV LJ ¼ 4

  12 r6  r 6 12 r r r

e

ð12Þ

In order to model the vdW force in the finite element, nonlinear spring elements were implemented, due to the nonlinear behavior of vdW force with atomic distance as can be seen in Eq. (12). The spring elements were located between each Carbon atoms and inner nodes of the polymer elements with the distance less than rcutoff ¼ 0:85 nm [26]. The FE model of the embedded CNT with the chirality of (10,10) and the volume fraction of 8% in the matrix is illustrated in Fig. 4. The dimensions of the RVE were calculated based on the CNT length and volume fraction. In the present study, for 8% volume fraction of CNT, for several CNT lengths the diameters of the RVEs were presented in Table 1. Noting that, the CNT in the molecular structural mechanics model was considered to be straight. Shi et al. [34] showed that for the case of randomly oriented CNT, the influence of CNT waviness is negligible on the effective properties of the composite. Hence, based on the assumption of random orientation which will be applied in Section 6, the CNT waviness was not modeled in the molecular structural mechanics.

where subscript ‘‘m” denotes ‘‘matrix”, Em and mm are uniaxial relaxation modulus and the Poisson’s ratio, respectively. To simulate viscoelastic properties of the matrix, the three-parameter standard linear solid model were used. In this model, the creep compliance and uniaxial relaxation modulus versus time can be expressed as [31]:

JðtÞ ¼

  1 q1  p1 q0 q  exp  0 t q0 q0 q1 q1

Em ðtÞ ¼ q0 þ

  q1  p1 q0 1 exp  t p1 p1

ð8Þ

ð9Þ

The creep compliance equation of the model (Eq. (8)) was fitted on the experimental data which was reported by Starkova et al. [32] and the constitutive parameters of the model were evaluated as:

Fig. 3. The result of the viscoelastic model and the experimental creep data [32] for epoxy resin.

399

A.R. Shajari et al. / Computational Materials Science 142 (2018) 395–409

Fig. 4. (a) Section view of the model, without showing vdW springs, (b) section view of cap region of the CNT, with showing vdW springs, (c) section view of cylindrical region of the CNT, with showing vdW springs.

Table 1 RVE dimensions. CNT length (nm) 10 20 RVE diameter (nm)

30

40

50

60

70

80

90

100

200

300

400

500

1

4.359

4.665

4.718

4.752

4.775

4.793

4.806

4.816

4.825

4.865

4.879

4.886

4.890

4.908

4.572

C 44 ¼ C 55 ¼ pEI ðtÞ;

4. Viscoelastic properties of the model The viscoelastic properties of the embedded CNT in the matrix which was modeled in the previous section, were evaluated in the present section. According to the transversely isotropic elastic properties of the CNTs, the effective viscoelastic properties of the model are transversely isotropic with the following form for the relaxation tensor:

2

C 11 6C 6 12 6 6 C 13 CEI ðtÞ ¼ 6 6 0 6 6 4 0 0

C 12

C 13

0

0

C 22

C 23

0

0

C 23

C 33

0

0

0

0

C 44

0

0

0

0

C 55

0

0

0

0

0

where nðtÞ, lðtÞ, kðtÞ, mðtÞ and pðtÞ are the axial relaxation modulus, cross relaxation modulus, transverse bulk relaxation modulus, transverse shear relaxation modulus, and axial shear relaxation modulus of the composites, respectively. 4.1. Determination of the viscoelastic properties

3

0 7 7 7 0 7 7 0 7 7 7 0 5 C 66

C 66 ¼ mEI ðtÞ

ð13Þ

To determine the properties of the embedded CNT in the matrix, five loading cases were applied on the model [11]. Considering time-dependent strain energy of the system as KðtÞ, volume and outer boundary of the model as V and B, the prescribed displacement/traction on the bounders of the model for each cases were presented in Table 2. 4.2. Effect of CNT length on the viscoelastic properties

C 11 ¼ nEI ðtÞ;

C 22 ¼ C 33 ¼ kEI ðtÞ þ mEI ðtÞ

C 12 ¼ C 13 ¼ lEI ðtÞ;

C 23 ¼ kEI ðtÞ  mEI ðtÞ

To model a single CNT with surrounding matrix, two assumptions may be used; the CNT length-dependent model (short CNT) and CNT length-independent model (long CNT). In the first one,

400

A.R. Shajari et al. / Computational Materials Science 142 (2018) 395–409

Table 2 Different loading cases of the model. Case No.

Displacement/Traction

1

u1 ðBÞ ¼ 0

u2 ðBÞ ¼ 2e x3

u3 ðBÞ ¼ 2e x2

2

u1 ðBÞ ¼ 0

u2 ðBÞ ¼ ex2

u3 ðBÞ ¼ ex3

kEI ðtÞ ¼ KðtÞ=2Ve2

3

u1 ðBÞ ¼ 2e x2

u2 ðBÞ ¼ 2e x1

u3 ðBÞ ¼ 0

pEI ðtÞ ¼ 2KðtÞ=Ve2

4

u1 ðBÞ ¼ ex1

u2 ðBÞ ¼ 0

u3 ðBÞ ¼ 0

5

u1 ðBÞ ¼ ex1

T 2 ðBÞ ¼ 0

T 3 ðBÞ ¼ 0

Relaxation modulus mEI ðtÞ ¼ 2KðtÞ=Ve2

2

nEI ðtÞ ¼ 2KðtÞ=Ve qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi lEI ðtÞ ¼ kEI ðtÞ½nEI ðtÞ  2KðtÞ=Ve2 

matrix were inserted in the ROM and properties of the EF with the assumption of a prefect bonding condition were obtained. In order to use formulations of the micromechanics for viscoelastic materials, the half-sided Fourier transform according to the following formula must be applied on the constitutive equations of the matrix and RVE.

~f ðxÞ ¼

Z

0

1

f ðtÞ expðixtÞdt

ð14Þ

The constitutive equations of the matrix and RVE were obtained as

CNTs do not have end capes and the matrix was modeled by a hollow cylinder. The loading cases with their boundary conditions at the ends were applied on both CNT and matrix nodes. In the one, the CNTs have end caps and fully embedded in the matrix. In this model, all loads and boundary conditions at the ends were applied on the matrix nodes. For a RVE contacting 8% CNT in the matrix, the results of the axial, cross, transverse bulk, transverse shear and axial shear relaxation moduli for different CNT length were presented in Figs. 5–9. From these figures, it can be found that, the transverse bulk and axial shear relaxation moduli were fairly independent of the CNT length and the results for long CNT can be used for different lengths. For the axial, cross and transverse shear relaxation moduli the effect of CNT length was considerable, especially for axial relaxation modulus. It was observed that the relaxation properties of the embedded CNT in the matrix were increased by increasing the length of CNT and converged to results of the long CNT. So, it can be stated that the upper bonds of the relaxation properties belong to the long CNT. For the case of long CNT, the viscoelastic properties of the model with CNT volume fractions 2%, 4%, 6%, 8% and 10% were obtained and listed in Table 3, based on the standard linear solid viscoelastic model. These values will be applied in Section 6 as inputs of the EI model. 5. Developing equivalent fiber In this section, equivalent fibers representing the CNT and its surrounding interface were conducted and corresponding properties were evaluated, using the results of structural mechanics simulation. To this end, the viscoelastic properties of the RVE and the

r~ k ðxÞ ¼ ixC~ k ðxÞ : ~ek ðxÞ; k ¼ m or RVE

ð15Þ

~ k ðxÞ are the stress, strain and relaxation ~ k ðxÞ, ~ek ðxÞ and C where r stiffness tensors of the matrix or RVE in a frequency domain, respectively. Defining the tensor of complex moduli as:

~ k ð xÞ Ck ðxÞ  ixC

ð16Þ

The constitutive equations of the matrix and RVE were expressed as:

r~ k ðxÞ ¼ Ck ðxÞ : ~ek ðxÞ

ð17Þ

The tensors of complex moduli contain real and imaginary parts where denote tensors of the storage and loss moduli, respectively as the following:

Ck ðxÞ ¼ C0k ðxÞ þ iC00k ðxÞ

ð18Þ

Here and throughout, over tilde, asterisk, prime and double prime indicate quantities at the frequency domain, complex form, storage and loss moduli, respectively.According to the ROM method, the stiffness tensor of the EFs can be presented as:

CEF ¼

1  ½C ðxÞ  ð1  cr ÞCm ðxÞ cr RVE

ð19Þ

Applying half-sided Fourier transform to the relaxation tensors of the matrix and the RVE, and using the magnitudes which were reported in Eq. (10) and Table 3, the transversely isotropic properties of EFs can be calculated as:

nEF ¼ 655:3; lEF ¼ 12:29; kEF ¼ 37:64; mEF ¼ 3:69; pEF ¼ 16:21 GPa

Fig. 5. Transverse bulk relaxation modulus for different CNT lengths.

ð20Þ

401

A.R. Shajari et al. / Computational Materials Science 142 (2018) 395–409

Fig. 6. Cross relaxation modulus for different CNT lengths.

Fig. 7. Axial relaxation modulus for different CNT lengths.

6. Micromechanical Mori-Tanaka approach To evaluate the overall viscoelastic properties of nanocomposites, the Mori-Tanaka method was applied in which EFs and EIs were considered as reinforcement phases for models I and II, respectively. Therefore, the stiffness tensor of EFs with the volume fraction cr and the tensor of complex moduli of EIs with the volume fraction cEI must be inserted in the Mori-Tanaka formulations for models I and II, respectively. According to Eqs. (15)–(18), the relaxation stiffness tensors in the frequency domain, tensors of complex, storage and loss moduli of the EIs are respectively denoted ~ EI ðxÞ, C ðxÞ, C0 ðxÞ and C00 ðxÞ. by C EI

EI

EI

Based on the Mori-Tanaka method, the overall properties of the composite in terms of the properties of its constituents can be stated as [35]:

CI ðxÞ ¼ ½ð1  cr ÞCm ðxÞ þ cr CEF : AI  : ½ð1  cr ÞI þ cr AI 1 CII ðxÞ ¼ ½ð1  cEI ÞCm ðxÞ þ cEI CEI ðxÞ : AII  : ½ð1  cEI ÞI þ cEI AII 1 ð21Þ where xÞ and xÞ are the effective tensors of complex moduli of the composite with aligned reinforcements for models I and II, respectively, and AI and AII are the dilute strain concentration tensors which are given by [35]: CI ð

CII ð

AI=II ¼ ½I þ S : ½Cm ðxÞ

1

1

: ½CEF=EI  Cm ðxÞ

ð22Þ

Where I and S are the fourth order identity and the Eshelby’s tensors. For cylindrical inclusions which were considered for the EFs and EIs, the non-zero components of tensor S can be stated as [35]:

402

A.R. Shajari et al. / Computational Materials Science 142 (2018) 395–409

Fig. 8. Transverse shear relaxation modulus for different CNT lengths.

Fig. 9. Axial shear relaxation modulus for different CNT lengths.

Table 3 The constants of the standard linear solid model for different volume fractions. cCNT

10%

8%

6%

4%

2%

Constants

p1 ð103 secÞ q0 ðGPaÞ q1 ðTPa:secÞ p1 ð103 secÞ q0 ðGPaÞ q1 ðTPa:secÞ p1 ð103 secÞ q0 ðGPaÞ q1 ðTPa:secÞ p1 ð103 secÞ q0 ðGPaÞ q1 ðTPa:secÞ p1 ð103 secÞ q0 ðGPaÞ q1 ðTPa:secÞ

Relaxation modulus ðGPaÞ kEI ðtÞ

lEI ðtÞ

mEI ðtÞ

nEI ðtÞ

pEI ðtÞ

3.312

3.813

3.185

3.489

3.227

6.272 24.00 3.312

3.642 16.77 3.813

1.911 6.614 3.185

68.84 243.7 3.489

2.171 7.548 3.227

5.691 21.90 3.334

3.445 15.72 3.813

1.450 5.252 3.185

55.86 198.8 3.489

1.874 6.596 3.204

5.066 20.22 3.312

3.233 14.84 3.813

1.098 4.140 3.185

42.84 155.8 3.489

1.565 5.628 3.227

4.488 17.88 3.312

3.034 14.24 3.813

0.851 3.347 3.185

29.83 108.3 3.489

1.261 4.715 3.227

3.844 16.21

2.819 13.23

0.704 2.866

16.87 62.87

0.959 3.701

403

A.R. Shajari et al. / Computational Materials Science 142 (2018) 395–409

S2222 ¼ S3333 ¼

5  4mm ; 8ð1  mm Þ

mm

S2211 ¼ S3311 ¼ S2233 S1212 S2323

Eq. (27) can be rewritten as following:

rk ðtÞ ¼ Dk1 : ek ðtÞ þ

2ð1  mm Þ 4mm  1 ¼ S3322 ¼ ; 8ð1  mm Þ 1 ¼ S1313 ¼ ; 4 3  4mm ¼ 8ð1  mm Þ

xÞ ¼ xÞ ¼

ð23Þ

Z

Z

1

ð24Þ where the angle bracket denotes the average value of the term for all orientations. The viscoelastic properties of the composite with randomly oriented reinforcements are isotropic with two structural constants: bulk and shear moduli, obtsained as follows [34]:

j x

x

xÞ

Þ  3 m ð ÞnI2 ð 3½ð1  cr Þ þ cr nI2 ð Þ cr ½nI3 ð Þ  2 m ð ÞnI4 ð 2½ð1  cr Þ þ cr nI4 ð Þ

cr ½nI1 ð

l x

x

x

x

ð25Þ

xÞ

 mð

7. Time incremental micromechanics For a direct applying of micromechanics formulations in the time domain to obtain the overall time dependent properties of composites, an iterative-recursive method was proposed in this section. At first, stress-strain relation of a viscoelastic material was presented by a time recursive equation. Next, an iterative based solution procedure was proposed to determine the overall properties. 7.1. Time recursive constitutive equation The relaxation tensor of a viscoelastic material obtained by experiments can be approximated by a well-known Prony series as:

Ck ðtÞ ¼ Dk1 þ

N X

Dkn expðkkn tÞ

ð26Þ

n¼1

where subscript/superscript ‘‘k” can be ‘‘m”, ‘‘r” or ‘‘EI”; and Dk1 is the tensor of the steady modulus. Substituting Eq. (26) into Eq. (6) yields:

rk ðtÞ ¼ Dk1 : ek ðtÞ þ

Z N X Dkn :

0

n¼1

t

exp½kkn ðt  sÞ

dek ðsÞ ds ds

ð27Þ

Defining tensor qkn ðtÞ as:

Z

qkn ðtÞ ¼

0

t

exp½kkn ðt  sÞ

dek ðsÞ ds ds

dek ðsÞ ds ds dek ðsÞ exp½kkn ðt  sÞ ds ds

exp½kkn ðt  sÞ

t

tDt

ð30Þ

Using Eq. (28), the first term of the right hand side of Eq. (30) was estimated as:

Z

tDt

0

exp½kkn ðt  sÞ

dek ðsÞ ds ¼ expðkkn tÞqkn ðt  DtÞ ds

ð28Þ

ð31Þ

Moreover, the second term of the right hand side of Eq. (30) was calculated as:

Z

t

tDt

exp½kkn ðt  sÞ

dek ðsÞ 1  expðkkn tÞ ½ek ðtÞ  ek ðt ds ¼ ds kkn t  DtÞ

ð32Þ

Based on Eqs. (31) and (32), the tensor qkn ðtÞ can be obtained as following:

qkn ðtÞ ¼ expðkkn tÞqkn ðt  DtÞ þ

1  expðkkn tÞ kkn t

½ek ðtÞ  ek ðt

 DtÞ

where j xÞ and l xÞ are the bulk and shear moduli of the matrix in the frequency domain. Similarly, bulk and shear moduli of model II were obtained by replacing subscript ‘I’ with ‘II’ in Eq. (25). The parameters nI1 ðxÞ to nI4 ðxÞ and nII1 ðxÞ to nII4 ðxÞ will be presented in Appendix A.  mð

tDt

0

þ

xÞ þ cr hCEF : Ai : ½ð1  cr ÞI þ cr hAI i xÞ þ cEI hCEI ðxÞ : Ai : ½ð1  cEI ÞI þ cEI hAi1

lI ðxÞ ¼ lm ðxÞ þ

The integral of Eq. (28) was separated by two integrals with different domains as:

qkn ðtÞ ¼

½ð1  cr ÞCm ð ½ð1  cEI ÞCm ð

jI ðxÞ ¼ jm ðxÞ þ

ð29Þ

n¼1

Eq. (21) gives the overall properties of the composite with the assumption of aligned EIs or EFs. For randomly oriented reinforcements, the average values of the stress and strain tensors over all orientations must be considered. So the Mori-Tanaka formulations were modified as [11]:

CI ð CII ð

N X Dkn : qkn ðtÞ

ð33Þ

Applying Eq. (33) in Eq. (29), the stress-strain relation for a viscoelastic material was reduced as the following recursive equation:

"

# N X 1  expðkkn tÞ Dkn kkn t n¼1 " # N k X k 1  expðkn tÞ : ek ðtÞ  Dn kkn t n¼1

rk ðtÞ ¼ Dk1 þ

: ek ðt  DtÞ þ

N X

expðkkn tÞDkn : qkn ðt  DtÞ

ð34Þ

n¼1

7.2. Solution procedure To use the Mori-Tanaka formulation in the time domain for a composite with viscoelastic constituents, it was needed to evaluate relaxation tensors at a given time with considering the history of load and/or deformation in the material. So, a time incremental iterative method was developed in which by applying a constant strain tensor and using the strain concentration tensor A, the strain tensors in the matrix and reinforcement were calculated. Due to the dependency of the tensor A to the relaxation tensor of the matrix and reinforcement (see Eq. (22)), a trial and error procedure was applied at the given time increment. Using this procedure, using Eq. (34), stress tensors of constituents and composites were obtained. Consequently, the overall relaxation of composites was evaluated. In order to obtain transversely isotropic relaxation properties of composites, three characteristic unit strain tensors were introduced:

e#1 ¼ ½ 1 0 0 0 0 0 T e#2 ¼ ½ 0 1 0 0 0 0 T e#3 ¼ ½ 0 0 0 1 0 0 T

ð35Þ

404

A.R. Shajari et al. / Computational Materials Science 142 (2018) 395–409

Fig. 10. The flowchart of computational solution procedure.

Fig. 11. Schematically view of chosing c CNT and cEI in model II: (a) large EI, (b) small EI.

These tensors were applied and corresponding stress tensors were obtained and the several relaxation moduli were calculated as follows:

r#1 ðtÞ ¼ ½ nðtÞ lðtÞ lðtÞ 0 0 0 T r#2 ðtÞ ¼ ½ lðtÞ kðtÞ þ mðtÞ kðtÞ  mðtÞ 0 0 0 T r#3 ðtÞ ¼ ½ 0 0 0 pðtÞ 0 0 T

ð36Þ

The proposed solution procedure was presented in the flowchart shown in Fig. 10. 8. Results and discussions In this section, the results of EF and EI models were compared in the frequency and time domains. It is to be mentioned, in model II, for a constant volume fraction of the CNTs in a composite (cr ), the results of the model must be independent of the EI’s size. Hence, in this section, for a constant cr , different sizes were considered for

EIs by considering different values for volume fractions of CNTs in the EIs (cCNT ) and the volume fraction of EIs in the matrix (cEI ). Therefore, the dependency of the results on the size of EI was investigated. Fig. 11 shows the model II for a constant cr and different sizes of EIs (i.e. magnitudes of cCNT and cEI ) schematically. For a 0.5 wt% CNT-reinforced composite, in which CNT volume fraction is cr ¼ 0:38%, the magnitudes of cCNT and cEI are taken in different cases of the EI’s size as presented in Table 4. Fig. 12(a) and (b) show the bulk storage and loss moduli, respectively for models I and II based on Eq. (25). It can be observed that the storage modulus of the composite is increased at least 8%, compared to the pure resin at different frequencies. However, for the loss modulus of the resin and composites, no considerable difference was found. The CNTs were assumed to be elastic, so the loss properties have not been considered for them. Consequently, it was concluded that the CNT only affects on the storage properties of composites. The difference between results of model II for various sizes of EIs is ignorable, so it can be stated that these results are independent of EI’s size. It is also found that,

405

A.R. Shajari et al. / Computational Materials Science 142 (2018) 395–409 Table 4 Different volume fractions cases of c CNT and c EI considered in model II for c r ¼ 0:38%. Case

1

2

3

4

5

cr (%) cCNT (%) cEI (%)

0.38 2 19.00

0.38 4 9.50

0.38 6 6.33

0.38 8 4.75

0.38 10 3.80

Fig. 12. (a) Bulk storage modulus, (b) Bulk loss modulus versus frequency for models I and II.

the predicted storage properties of model I are overestimated in comparison with those of model II. For the shear storage and loss moduli, the above statements are valid as can be observed in Fig. 13. In the following, the viscoelastic properties of composites were presented in a time-dependent form, applying the inverse Fourier transform to the bulk and shear moduli which was mentioned in Eq. (25). Fig. 14(a) and (b) present the bulk and shear relaxation moduli for models I and II. Independency of the results respect to the EI’s size for the results of the model II, can be found in these figures. In addition, it is observed that for the relaxation properties, results of model I are overestimated compared to those of model II. To compare the predicted properties of the models with experimental data, the results of uniaxial creep compliance tests reported by Starkova et al. [32] for 0.5 wt% CNT reinforced epoxy were used. This comparison was shown in Fig. 15. As can be seen in this figure, both models underestimate the creep behavior in comparison with the experimental data. Therefore, it was found that the relaxation properties, predicted by the models are overestimated in comparison with experimental

results. Moreover, the prediction of properties can be proportionally improved when equivalent inclusion (model II) model is utilized. The better prediction of model II is the consequence of considering viscoelastic properties for the reinforcement in the micromechanical modeling. Applying viscoelastic properties of EIs to estimate the overall viscoelastic properties of the nanocomposite leads to increase creep response compared to the case of applying elastic properties of EFs. Therefore, the improved prediction of model II in comparison with experimental data was observed. The difference between the results of both models with experimental data can be justified by nanoscale phenomena such as the CNTs waviness and agglomeration that reduce the effective stiffness properties of nanocomposites. To compare the results of the molecular structural mechanics modeling of an embedded CNT in the polymer matrix with available results of MD simulation, the study of Zhu et al. [36] was considered. In their study, a MD simulation for a singlewalled CNT reinforced Epon 862 was carried out and the stressstrain curve of unidirectional composite was obtained. To calculate

406

A.R. Shajari et al. / Computational Materials Science 142 (2018) 395–409

Fig. 13. (a) Shear storage modulus, (b) Shear loss modulus versus frequency for models I and II.

the stress-strain behavior of the RVE in FE model, all nodes at one end were axially constrained and a uniform axial displacement was applied at the other end. The mean value of the reaction forces was obtained and implemented for evaluation of stress. The stressstrain curves of MD simulation and present model were shown in Fig. 16. From this figure, it can be found that the results of the structural mechanics have satisfactory agreement with those obtained by MD simulation. The properties of the matrix in the present model, was an input of the model. But in MD, these properties were determined from the results of the simulation. Hence, the difference between results of the two approaches can be justified. To evaluate results of the Mori-Tanaka time incremental method, an aligned unidirectional composite with 0.5 wt% CNT reinforcement was considered and five independent relaxation moduli were obtained. Likewise, these moduli were estimated using the complex Mori-Tanaka formulations and compared with the results of the time incremental method. 4The results were presented in Fig. 17. For the present solution, the time increment was considered to be: Dt ¼ 0:01sec, and the convergence condition was:

max

ði1Þ ði1Þ kCm  CðiÞ  CrðiÞ k2 m k2 kCr ; ðiÞ ðiÞ kCm k2 kCr k2

! < 0:01

ð37Þ

where kBk2 returns the 2-norm of the B matrix. As can be observed from this figure, the predicted properties of time incremental method were satisfactory agreed with those

obtained by the complex micromechanics. Therefore, the computational iterative solution was verified.

9. Conclusion In the present study, a multiscale modeling for simulation of viscoelastic properties of CNT reinforced polymer nanocomposites was presented. At the nano scale, a pristine CNT, embedded in the matrix was modeled using a FE-based molecular structural mechanics and corresponding transversely isotropic relaxation properties were obtained. At the micro scale, two modeling approaches were implemented. In model I, the CNT and the interface were modeled by an equivalent fiber. Moreover, elastic properties of the equivalent fiber were estimated using the rule of mixture method. In model II, which was developed in the present research, the whole CNT, interface and the matrix were modeled by an equivalent inclusion with the same geometries and properties. Finally, the Mori-Tanaka method was used to predict the overall viscoelastic properties of composites. The FE-based model of the embedded CNT in the matrix was developed using the assumptions of short and long CNTs. It was found that the axial relaxation modulus was strongly dependent on the length of the CNT. This modulus is increased by increasing the CNT length and tends to the result of the long CNT model. Other relaxation moduli, such as cross and transverse shear moduli have lesser dependency on the length of the CNT compared to the axial modulus. The relaxation moduli of the axial shear and trans-

A.R. Shajari et al. / Computational Materials Science 142 (2018) 395–409

407

Fig. 14. (a) Bulk relaxation modulus, (b) Shear relaxation modulus for models I and II.

Fig. 15. Comparison between creep compliance of the models with the available experimental dada.

Fig. 16. Comparison between stress strain curve of MD results [36] with present model.

verse bulk can be considered approximately independent of the CNT length. The investigation of the results dependency of model II on the size of the equivalent inclusion showed that the size did not have significant effects on the predicted properties by the model.

Therefore, it can be stated that this multiscale model is capable of prediction of the viscoelastic properties of CNT/polymer composites. A comparison between results of models I and II indicated that for relaxation and storage properties, the results of model I are

408

A.R. Shajari et al. / Computational Materials Science 142 (2018) 395–409

Fig. 17. A comparison of independent relaxation moduli of the time incremental and complex methods: (a) Transverse bulk, (b) Cross, (c) Axial, (d) Transverse shear, (e) Axial shear.

higher than those of the model II. While for the creep properties, the prediction of model I is lesser than those of obtained by model II. On the other hand, using available experimental data and comparing with the results of the two models, it was revealed a more realistic prediction of model II. The proposed time incremental formulations of micromechanics and the corresponding solution procedure are efficient tools for prediction of the viscoelastic properties of composites in the time domain without using any transformation.

Appendix A Considering the stiffness tensor of EFs as the following form:

2 6 6 6 6 CEF ðxÞ ¼ 6 6 6 6 4

nEF

lEF

lEF

0

0

lEF

kEF þ mEF

kEF  mEF

0

0

lEF

kEF  mEF

kEF þ mEF

0

0

0

0

0

pEF

0

0

0

0

0

pEF

0

0

0

0

0

0

3

0 7 7 7 0 7 7 0 7 7 7 0 5 mEF

A.R. Shajari et al. / Computational Materials Science 142 (2018) 395–409

and the complex moduli tensor of EIs as the following form:

2

nEI 6 l 6 EI 6  6 lEI 6

CEI ðxÞ ¼ 6 6 0 6 4 0

 lEI

0

0

kEI þ mEI

kEI  mEI

0

0

kEI  mEI

kEI þ mEI

0

0

0

0

pEI

0

0

0

0

pEI



 lEI



0

 

0

The parameters

0

xÞ to

nI1 ð

nI4 ð

0

0

0

3

0 7 7 7 0 7 7 0 7 7 7 0 5 mEI

xÞ in Eq. (25) are obtained as

nI1 ðxÞ ¼

  1 ð2kEF þ lEF Þð3jm þ 2lm  lEF Þ nEF þ 2lEF þ 3 lm þ kEF

nI2 ðxÞ ¼

3ðj þ l 3ðl  m

 m Þ þ kEF  þk Þ EF m

 lEF

 1 2 8l  p 8mEF lm ð3jm þ 4lm Þ ðnEF  lEF Þ þ  m EF þ  5 3 lm þ pEF 3jm ðmEF þ lm Þ þ lm ð7mEF þ lm Þ   2ðkEF  lEF Þð2lm þ lEF Þ þ 3ðlm þ kEF Þ

nI3 ðxÞ ¼

[12]

[13]

[14]

[15] [16]

[17]

[18] [19]

[20]



1 4lm þ 2kEF þ lEF 4l þ  m  5 3ðlm þ kEF Þ lm þ pEF  lm ð3jm þ lm Þ þ lm ð3jm þ 7lm Þ þ2  lm ð3jm þ lm Þ þ mEF ð3jm þ 7lm Þ

nI4 ðxÞ ¼

[11]

Likewise, the parameters nII1 ðxÞ to nII4 ðxÞ can be calculated by   replacing nEF , lEF , kEF , mEF and pEF with nEI , lEI , kEI , mEI and pEI , respectively in above formulas. References [1] G.D. Seidel, D.C. Lagoudas, Micromechanical analysis of the effective elastic properties of carbon nanotube reinforced composites, Mech. Mater. 38 (8) (2006) 884–907. [2] B. Arash, Q. Wang, V.K.Varadan, Mechanical properties of carbon nanotube/ polymer composites, Scientific reports. 2014, 4. [3] A. Alva, A. Bhagat, S. Raja, Effective moduli evaluation of carbon nanotube reinforced polymers using micromechanics, Mech. Adv. Mater. Struct. 22 (10) (2015) 819–828. [4] J. Pan, L. Bian, H. Zhao, Y. Zhao, A new micromechanics model and effective elastic modulus of nanotube reinforced composites, Comput. Mater. Sci. 113 (2016) 21–26. [5] M. Esteva, P. Spanos, Effective elastic properties of nanotube reinforced composites with slightly weakened interfaces, J. Mech. Mater. Struct. 4 (5) (2009) 887–900. [6] S. Yang, S. Yu, J. Ryu, J.M. Cho, W. Kyoung, D.S. Han, M. Cho, Nonlinear multiscale modeling approach to characterize elastoplastic behavior of CNT/ polymer nanocomposites considering the interphase and interfacial imperfection, Int. J. Plast 41 (2013) 124–146. [7] P. Barai, G.J. Weng, A theory of plasticity for carbon nanotube reinforced composites, Int. J. Plast 27 (4) (2011) 539–559. [8] B. Coto, I. Antia, J. Barriga, M. Blanco, J.R. Sarasua, Influence of the geometrical properties of the carbon nanotubes on the interfacial behavior of epoxy/CNT composites: a molecular modelling approach, Comput. Mater. Sci. 79 (2013) 99–104. [9] K. Sharma, M. Shukla, Molecular modeling of the mechanical behavior of carbon fiber-amine functionalized multiwall carbon nanotube/epoxy composites, New Carbon Mater. 29 (2) (2014) 132–142. [10] Y. Zhang, X. Zhuang, J. Muthu, T. Mabrouki, M. Fontaine, Y. Gong, T. Rabczuk, Load transfer of graphene/carbon nanotube/polyethylene hybrid

[21]

[22]

[23]

[24]

[25] [26]

[27]

[28] [29] [30] [31] [32]

[33]

[34]

[35] [36]

409

nanocomposite by molecular dynamics simulation, Compos. B Eng. 63 (2014) 27–33. G.M. Odegard, T.S. Gates, K.E. Wise, C. Park, E.J. Siochi, Constitutive modeling of nanotube–reinforced polymer composites, Comp. Sci. Technol. 63 (11) (2003) 1671–1687. M.M. Shokrieh, R. Rafiee, Prediction of mechanical properties of an embedded carbon nanotube in polymer matrix based on developing an equivalent long fiber, Mech. Res. Commun. 37 (2) (2010) 235–240. M.M. Shokrieh, R. Rafiee, Investigation of nanotube length effect on the reinforcement efficiency in carbon nanotube based composites, Compos. Struct. 92 (10) (2010) 2415–2420. M.M. Shokrieh, R. Rafiee, On the tensile behavior of an embedded carbon nanotube in polymer matrix with non-bonded interphase region, Compos. Struct. 92 (3) (2010) 647–652. Pavliotis GA, Stuart A. Multiscale methods: averaging and homogenization. Springer Science & Business Media; 2008. D.N. Savvas, V. Papadopoulos, M. Papadrakakis, The effect of interfacial shear strength on damping behavior of carbon nanotube reinforced composites, Int. J. Solids Struct. 49 (26) (2012) 3823–3837. E. Mohammadpour, M. Awang, S. Kakooei, H.M. Akil, Modeling the tensile stress–strain response of carbon nanotube/polypropylene nanocomposites using nonlinear representative volume element, Mater. Des. 58 (2014) 36–42. P. Joshi, S.H. Upadhyay, Evaluation of elastic properties of multi walled carbon nanotube reinforced composite, Comput. Mater. Sci. 81 (2014) 332–338. K.N. Spanos, S.K. Georgantzinos, N.K. Anifantis, Investigation of stress transfer in carbon nanotube reinforced composites using a multi-scale finite element approach, Compos. B Eng. 31 (63) (2014 Jul) 85–93. J.M. Wernik, S.A. Meguid, Multiscale micromechanical modeling of the constitutive response of carbon nanotube-reinforced structural adhesives, Int. J. Solids Struct. 51 (14) (2014) 2575–2589. X. Liu, Q.S. Yang, X.Q. He, K.M. Liew, Cohesive laws for van der Waals interactions of super carbon nanotube/polymer composites, Mech. Res. Commun. 72 (2016) 33–40. K. Li, X.L. Gao, A.K. Roy, Micromechanical modeling of viscoelastic properties of carbon nanotube-reinforced polymer composites, Mech. Adv. Mater. Struct. 13 (4) (2006) 317–328. X.Y. Zhu, X. Wang, Y. Yu, Micromechanical creep models for asphalt-based multi-phase particle-reinforced composites with viscoelastic imperfect interface, Int. J. Eng. Sci. 76 (2014) 34–46. Y. Pan, G.J. Weng, S.A. Meguid, W.S. Bao, Z.H. Zhu, A.M. Hamouda, Interface effects on the viscoelastic characteristics of carbon nanotube polymer matrix composites, Mech. Mater. 58 (2013), 1 1. R. Li, L.Z. Sun, A micromechanics-based viscoelastic model for nanocomposites with imperfect interface, Int. J. Damage Mech 22 (7) (2013) 967–981. M.M. Shokrieh, R. Ghajar, A.R. Shajari, The effect of time-dependent slightly weakened interface on the viscoelastic properties of CNT/polymer nanocomposites, Compos. Struct. 146 (2016) 122–131. M. El Kouri, A. Bakkali, L. Azrar, Mathematical modeling of the overall timedependent behavior of non-ageing viscoelastic reinforced composites, Appl. Math. Model. (2015), https://doi.org/10.1016/j.apm.2015.11.031. H.F. Brinson, L.C. Brinson, Polymer engineering science and viscoelasticity, Springer, Berlin, 2008. Hibbitt, Karlsson, Sorensen, ABAQUS/standard user’s Manual (vol. 1). Hibbitt, Karlsson & Sorensen, 2001. C. Li, T.W. Chou, A structural mechanics approach for the analysis of carbon nanotubes, Int. J. Solids Struct. 40 (10) (2003) 2487–2499. R.S. Lakes, Viscoelastic materials, Cambridge University Press, 2009. O. Starkova, S.T. Buschhorn, E. Mannov, K. Schulte, A. Aniskevich, Creep and recovery of epoxy/MWCNT nanocomposites, Compos. A Appl. Sci. Manuf. 43 (8) (2012) 1212–1218. K.I. Tserpes, P. Papanikos, G. Labeas, S.G. Pantelakis, Multi-scale modeling of tensile behavior of carbon nanotube-reinforced composites, Theoret. Appl. Fract. Mech. 49 (1) (2008) 51–60. D.L. Shi, X.Q. Feng, Y.Y. Huang, K.C. Hwang, H. Gao, The effect of nanotube waviness and agglomeration on the elastic property of carbon nanotubereinforced composites, J. Eng. Mater. Technol. 126 (3) (2004) 250–257. S. Nemat-Nasser, M. Hori, Micromechanics: overall properties of heterogeneous materials, Elsevier, 2013. R. Zhu, E. Pana, A.K. Roy, Molecular dynamics study of the stress–strain behavior of carbon-nanotube reinforced Epon 862 composites, Mater. Sci. Eng., A 447 (2007) 51–57.