Multiscale homogenization model for thermoelastic behavior of epoxy-based composites with polydisperse SiC nanoparticles

Composite Structures 128 (2015) 342–353

Contents lists available at ScienceDirect

Composite Structures journal homepage: www.elsevier.com/locate/compstruct

Multiscale homogenization model for thermoelastic behavior of epoxy-based composites with polydisperse SiC nanoparticles Seongmin Chang a, Seunghwa Yang b, Hyunseong Shin a, Maenghyo Cho a,⇑ a b

Department of Mechanical and Aerospace Engineering, Seoul National University, San 56-1, Shillim-Dong, Kwanak-Ku, Seoul 151-742, South Korea Department of Mechanical Engineering, Dong-A University, Nakdong-daero 550 beon-gil, Saha-Ku, Busan 604-714, South Korea

a r t i c l e

i n f o

Article history: Available online 26 March 2015 Keywords: Multiscale Homogenization Nanocomposites

a b s t r a c t A multiscale statistical homogenization method based on a finite element analysis is proposed to predict the embedded filler size-dependent thermoelastic properties of nanoparticulate polymer composites. A molecular dynamics simulation is used to predict the particle size-dependent elastic constants and coefficient of thermal expansion (CTE) of nanocomposites. Because of the densified interphase zones formed in the vicinity of nanoparticles, and their relative dominance according to particle size, the thermoelastic properties are more prominently increased by smaller nanoparticles. The equivalent continuum microstructure of a nanocomposite is modeled as a three-phase periodic unit cell consisting of the nanoparticle, matrix, and effective interphase zone. An inverse numerical scheme is proposed to predict the thermoelastic properties of the effective interphase zone from the known properties of nanocomposites. In order to account for more realistic variation of the nanoparticle radius and spatial distribution in nanocomposites, statistical homogenization is performed by assigning randomness to the particle size using a beta distribution and to the spatial distribution through larger many-particle embedded representative volume elements (RVEs). Compared with the result from a regular homogenization method using a mono-particle RVE, the mean value of the elastic moduli increases, while that of the CTE decreases in the many-particle RVEs whose mean particle radius corresponds to that of the mono-particle RVE at the same volume fraction of nanoparticles. Ó 2015 Elsevier Ltd. All rights reserved.

1. Introduction Polymer nanocomposites have gained considerable attention in both industrial and academic circles studying advanced functional materials because of their exceptional properties, which allow the use of nanostructured materials as reinforcements [1–4]. In nanocomposites, the structural arrangement and packing of matrix molecules near the surface of a nanofiller are known to be altered, and a unique phase is formed, referred to as the ‘‘interphase,’’ whose material properties are different from those of the bulk matrix region [2,5–11]. With smaller nanoparticles, therefore, the relative volume fraction of such an interphase zone increases and a better reinforcing effect can be anticipated by reducing the size of the nanoparticles. In most of the studies that have reported the filler size-dependent properties of nanocomposites, these filler size-dependent properties appeared in nano-sized particles and disappeared in micron or larger scale particles [2]. In this respect, controlling the size of the embedded nanofiller is a very important ⇑ Corresponding author. Tel.: +82 2 880 1693; fax: +82 2 886 1693. E-mail address: [email protected] (M. Cho). http://dx.doi.org/10.1016/j.compstruct.2015.03.041 0263-8223/Ó 2015 Elsevier Ltd. All rights reserved.

design variable to maximize the weight advantage of nanocomposites compared with micron or larger scale composites. Among the various applications of polymer nanocomposites, the underfill for an electronic packaging system requires an excellent elastic moduli and thermal expansion coefficient of the polymer to prevent interfacial delamination from a rigid substrate. In this application, tailoring the thermoelastic properties by adding nanofillers is a promising way to guarantee the reliability of the electronic packaging system. Regarding the thermoelastic properties of nanocomposites, Goyal et al. [12] compared the thermal expansion coefficient of nano- and micro-alumina-reinforced poly-ether-ether-ketone (PEEK) composites and revealed the superiority of nanoparticles to decrease the thermal expansion coefficient of pure PEEK. Sun et al. [13] compared the dimensional stabilities of micro- and nano-silica-reinforced epoxy composites and demonstrated that the nanofillers could decrease the coefficient of thermal expansion (CTE) of a base polymer better than micron-sized fillers. Goertzen and Kessler [14] reported the same experimental result for the CTE of nanosilica/cyanate ester nanocomposites. More recently, using the same silica/epoxy composites, Jang et al. [15] reported a clear filler size effect in the

S. Chang et al. / Composite Structures 128 (2015) 342–353

elastic moduli and CTE of nanocomposites, where smaller nanoparticles could increase (decrease) the elastic moduli (CTE) of an epoxy better than larger ones. Using molecular dynamics (MD) simulations, Choi et al. [9] demonstrated the same tendency in nano-SiC/epoxy composites and correlated the size-dependent thermoelastic properties of the nanocomposites with the interphase zone formed right in the vicinity of a SiC nanoparticle. To improve the design of polymer nanocomposites, however, the direct application of the conventional continuum theory has limitations in describing the particle size effect. To resolve this limitation, proposals have been made for micromechanics-based multiscale approaches and a noble inverse estimation method to characterize the properties of the interphase zone, which plays a key role in tailoring the properties of nanocomposites [5–7,16]. However, micromechanics models based on Eshelby’s solution [17] have some limitations in the rigorous consideration of various particle sizes and spatial distributions, as well as a high volume fraction condition. As pointed out by Yang and Cho [6], the equivalent continuum model of a nanoparticle surrounded by an interphase zone in a nanocomposite can be interpreted as a coated rigid inclusion. Thus, because of the interphase zone, the effective volume fraction of the nanoparticle when including the interphase zone is higher than the real volume fraction of the nanoparticle itself. Moreover, the precise control of the nanoparticle size in nanocomposites is still challenging, and there is some variation in the nanoparticle size in real nanocomposites. Therefore, most of experimental studies have adopted the mean particle size rather than a definite value for it [18,19]. This situation makes the analytical estimation of the effective properties of nanocomposites with a definite value of the particle size further away from real experimental observations [20]. Motivated by the above requirement to establish an efficient and reliable analysis and design scheme for nanocomposites, the present study proposed a stochastic homogenization method to characterize the thermoelastic properties of polydisperse SiC/ epoxy nanocomposites across multiple spatial scales. To capture the nanoparticle size-dependent elastic moduli and CTEs of nanocomposites, we utilized MD simulations of the molecular structures of monoparticulate nanocomposites with different embedded nano-SiC particles in the lowest scale. To reproduce the MD simulation results from the homogenization method, an inverse numerical scheme was proposed to estimate the thermoelastic properties of the interphase zone. By defining the properties of the interphase as a function of the particle size, a full continuum scale stochastic analysis was implemented to account for the random variation of the particle size and spatial distribution inside the many-particle embedded nanocomposite representative volume element (RVE). Together with the randomness of the size and distribution of nanoparticles, randomness in cross link density and mechanical properties of epoxy matrix, imperfect interfacial bonding between SiC and epoxy can also generate stochastic variation of the properties of composites. To understand the effect of polydisperse nanoparticle more clearly, however, we let the randomness of other remaining uncertain factors and properties in composites microstructures aside. To evaluate the validity of the present stochastic analysis, the mean and probability distribution of the thermoelastic properties of the nanocomposite were compared with those from regular homogenization, MD-based stochastic analysis, and micromechanics-based stochastic analysis. 2. Molecular dynamics simulation In all the molecular modeling and simulations of the nanocomposites and other constituents, Material StudioÒ [21] was used, which is a commercially available molecular modeling and simulation package. To describe the inter- and intra-atomic potentials,

343

the condensed-phase optimized molecular potential for atomistic simulation studies (COMPASS) force field was used, which was parameterized using ab initio calculations [22]. 2.1. Preparation of nanocomposite unit cell The nanocomposite molecular structure consisted of a spherical SiC nanoparticle embedded at the center of a unit cell and crosslinked epoxy molecules. The cross-linked epoxy consisted of a diglycidyl ether of bisphenol F(EPON862Ò) and a curing agent, triethylenetetramine (TETA). In contrast to thermoplastic composites, the modeling of a thermoset polymer requires an additional process to form cross-linked networks. Regarding the modeling of an epoxy, three major approaches have been applied: a Monte Carlo method [23], dynamic in situ cross-linking method based on continuous close contact detection [24–27], and representative molecule method [7,9,10,28]. With the goal of chemical rigorousness and the realistic modeling of the cross-linked network structure of the epoxy, close contact monitoring between candidate carbons in the epoxide groups and nitrogen or other atoms in the curing agent molecules has been preferably chosen in molecular dynamics simulations. However, this method requires a longer simulation time and careful monitoring of ring catenation and spearing to prevent an unrealistic molecular configuration after the cross-linking process. In modeling epoxy-based polymer composites, because of the existence of the embedded filler molecule inside the unit cell, reaching the target conversion ratio also becomes more difficult than in the case of a pure epoxy matrix. On the other hand, even if the randomness of the real cross linking process is not fully taken into account, the representative molecule method can be a useful alternative to study a cross-linked epoxy and an epoxybased nanocomposite structure. According to the previous studies [7,29], the representative molecule method can provide the comparable mechanical properties and thermal stability of a pure epoxy matrix with simplicity, as well as accuracy. In this respect, we used the representative molecule method with a curing ratio of 50%. As the name of this method indicates, the representative molecule method was used to construct a bulk periodic epoxy unit cell using several identical and pre-cured molecules composed of EPON862Ò and TETA molecules. Because the TETA molecule has six nitrogen atoms, which can be bonded to the carbon atoms in the epoxide groups of EPON862Ò, various cross-linked structures can be made by selectively choosing the cross-linking sites. The representative cross-linked molecule considered in this study consisted of nine EPON862Ò molecules and three TETA molecules, as shown in Fig. 1(a). Using equal numbers of each molecule, other types of cross-linked structures such as ladder-like networks or triangular structures are also possible. In our simulation, however, the linear linkage structure shown in Fig. 1(a) showed the best performance to obtain thermoelastic properties for the epoxy that were very comparable to the experimental values from molecular dynamics simulations. The spherical SiC nanoparticle was made from a bulk SiC crystal by trimming the carbon and silicon atoms out of the predefined radius, and is shown in Fig. 1(b). To investigate the particle size dependency of the thermoelastic properties, five different nanoparticles were prepared, with radii in the range of 5.18–11.63 Å. The nanoparticle volume fraction was fixed at 5.8%. No surface functionalization or valence grafting between the nanoparticle and epoxy matrix were considered in this study. Once the nanoparticle molecular structures were prepared, they were mixed into the representative epoxy molecules to make nanocomposite unit cells. The nanocomposite unit cells were constructed using the Amorphous cellÒ module of Material StudioÒ [21]. Since EPON862Ò has two hexagonal ring structures, ring spearing and catenation may occur when packing the representative molecules

344

S. Chang et al. / Composite Structures 128 (2015) 342–353

(a) Cross-linked epoxy representative molecule

(b) SiC nanoparticle

(c) Nanocomposite unit cell Fig. 1. Molecular models of cross-linked epoxy representative molecule, SiC nanoparticle, and nanocomposite unit cell for MD simulation.

into the unit cell. Thus, we defined the initial and final target densities of 0.1 g/Cm3 and 1.23 g/Cm3 for the nanocomposites, and 0.1 g/Cm3 and 1.0 g/Cm3 for pure epoxy molecular structures, respectively. Once the initial molecular configurations were determined at the initial density, the volume of the unit cell was gradually decreased in steps. At each step of cell densification, 1000 iterations of the conjugate gradient energy minimization process were implemented, followed by 1000 steps of isothermal ensemble simulation (the NVT ensemble simulation, N: total number of atoms in a unit cell, V: volume of the unit cell, T: temperature of the unit cell) at 500 K via the Andersen thermostat, with a time step of 1 fs [30]. In the unit cell construction process, periodic boundary conditions were imposed in all directions to prevent the finite size effect in the nanocomposites. The constructed nanocomposite periodic unit cell is depicted in Fig. 1(c), and the details of the unit cell composition and dimensions are summarized in Table 1.

simulation was implemented, followed by a 1,000-ps isothermalisobaric ensemble simulation (the NPT ensemble simulation, P: internal pressure tensor of the unit cell) adopting the Berendsen barostat [31]. During the NPT ensemble simulation process, the convergence of the unit cell dimension was stable, and the volume fraction of all the unit cells could be successfully conserved. The elastic constants were calculated using the Parinello– Rahmann strain fluctuation method [32]. This was combined with the constant stress ensemble simulation (the N rT ensemble simulation, r: internal stress tensor of the unit cell), which allowed anisotropic cell swelling [33], as well as distortion. The elastic stiffness tensor can be calculated based on the variation of the unit cell dimension and thermodynamic properties of the unit cell during the N rT ensemble simulation, as follows:

2.2. Calculation of thermoelastic properties

where k is the Boltzmann constant, T is the temperature of the unit cell, V is the volume of the unit cell, and the bracket hi indicates that the scalar or vector quantities are ensemble averages during the N rT simulation. Before calculating the elastic stiffness tensor, the unit cell was further equilibrated for 600 ps using the NrT ensemble simulation at 300 K and 1 atm of internal stress. Then, 100 ps of the N rT simulation was again implemented to collect the strain fluctuations. In this study, we used the final 100 ps of the N rT simulation to collect 10,000 strain fluctuations. For computational accuracy, the final production run to collect the strain fluctuations was repeated three times using different random number seeds to allocate the initial temperatures of individual atoms. Then, the elastic stiffness tensor was averaged over the three different simulation results. In addition to the strain fluctuation method adopted in this study, a uniaxial tension simulation using the MD simulations can be used to predict the elastic constants of the unit cell structures of interest. By applying a constant rate for the external strain or stress to the unit cell with a suitable boundary

After preparing the initial unit cell structures, the elastic constants and CTEs of the nanocomposites, pure epoxy, and SiC bulk crystal were calculated using ensemble simulations. To equilibrate the unit cell structures at 300 K and 1 atm, before a production run to calculate the thermoelastic properties, a 100-ps NVT ensemble

Table 1 Details of unit cell composition and dimensions of SiC/epoxy nanocomposites. Radius (Å)

No. of cross-linked epoxies

Cell length (Å)

Volume Fraction (%)

Cell density (g/cm3)

7.54 9.00 10.00 10.90 11.63

6 10 14 18 22

31.40 37.32 41.73 45.38 48.41

5.8

1.23

C ijkl ¼

1 kT  deij dekl hV i

ð1Þ

345

S. Chang et al. / Composite Structures 128 (2015) 342–353

condition to mimic a real experimental test, the elastic constants of the nanocomposites and other periodic structures can be calculated from the resultant stress or strain of the unit cell corresponding to the external strain or stress. The CTEs of the nanocomposites were calculated based on the results of the cooling down simulation from 340 K to 280 K, which was the glassy region for the pure epoxy matrix. After equilibrating at 300 K for 1,000 ps, the unit cells were heated up to 340 K and again equilibrated for 500 ps using the NPT ensemble simulation. Then, the unit cells were gradually cooled down with a cooling rate of 20 K/500 ps. At each cooling down stage, the density or volume of the unit cells could be monitored to determine the temperature– volume (or density) relation and obtain the CTEs of the nanocomposites and other molecular structures. The linear CTE could then be obtained from the slope of the temperature–volume relation obtained from the cooling-down simulation as follows:

aL ¼

aV 3

¼

1 @V 3V 0 @T

ð2Þ

where aL and aV are the linear and volumetric CTEs, and V 0 is the initial unit cell volume. The slope of the temperature–volume relation was obtained from the least square approximations. In a similar manner, the linear CTEs of the pure epoxy and SiC crystal were calculated using the multiscale homogenization methods. The linear CTEs of the nanocomposites were averaged over five different cooling-down simulations. 3. Homogenization theory 3.1. Fundamental formulation The homogenization method [34–37] and mean field theory have been the two main methods used to describe the overall behavior of heterogeneous materials. While the mean field theory is a physics-based theory, the homogenization method establishes a rigorous mathematical foundation. One major purpose of the homogenization method is to calculate the equivalent macroscopic homogeneous properties of heterogeneous structures such as porous materials or composites with periodic microstructures. To describe both the microscopic and macroscopic behaviors of nanocomposites, it is convenient to distinguish micro- and macro-coordinate systems as follows:

X ¼ Xðx; yÞ;

x ¼ X;

y ¼ X=e

ð3Þ

where the scale ratio e is the ratio of the microscopic representative length x to the macroscopic length y. In this study, we defined a microscopic-level unit cell of the nanocomposite as a three-phase structure consisting of the nanoparticle, interphase, and matrix zones, as shown on the right-hand side of Fig. 2. Then, the overall

macroscopic domain of the nanocomposite could be defined as periodic simple cubic (SC) aggregates, as shown on the left-hand side of Fig. 2. In Fig. 2, V X is the overall macroscopic domain of the nanocomposite, except for the nanoparticle and matrix; At is the  y is the boundary where the external traction is prescribed; and V unit cell domain, except for the nanoparticle and interphase zones. In the homogenization method, the macro-scale displacement field is asymptotically expanded with respect to the scale ratio as follows:

uðXÞ ¼ u0 ðx; yÞ þ eu1 ðx; yÞ þ e2 u2 ðx; yÞ þ    :

ð4Þ

For a thermoelastic problem, the virtual work principles of a nanocomposite subjected to traction and thermal loads are given as follows:

Z

Z Z @uk @ v i @v i dV ¼ t i v i dA þ C ijkl akl ðT  T 0 Þ dV @X l @X j @X j VX At VX @wðx; yÞ @w 1 @w ¼ þ ; wðx; yÞ ¼ uðx; yÞ or v ðx; yÞ @X i @xi e @yi C ijkl

ð5Þ

where uk and v i are the real and virtual displacements; t i is the traction; T and T 0 are the current and reference temperatures, respectively; and C ijkl and akl are the fourth- and second-order elastic stiffness and thermal expansion tensor of the nanocomposites, respectively. Since local coordinate y can be expressed in terms of the scale ratio e, all the partial derivatives of the displacement fields, which are a function of the coordinates x and y, can be expanded as depicted in Eq. (5). By substituting Eq. (4) into the governing equation, Eq. (5), the governing equation can be arranged by the power of the parameter e, as shown in Eq. (6):

Oð1=e2 Þ

1

e2

1 þ Oð1=eÞ þ Oð1Þ þ    ¼ 0

ð6Þ

e

Eq. (6) is always satisfied regardless of the value of parameter e. Thus, the following equations can be derived:

Z y V

Z y V

Z VX

C ijml

@ vkl m @dv i dV y ¼ @yn @yj

C ijkl

@/k @dv i dV y ¼ @yl @yj

CHijkl

@uk @ v i dV ¼ @xl @xj

Z y V

Z

Z At

y V

C ijml

@v i dV y @yj

C ijkl akl

ð7Þ

@dv i dV y @yj

t i du0i dA þ

Z VX

CHijkl aHkl ðT  T 0 Þ

ð8Þ @v i dV @xj

ð9Þ

In Eqs. (7) and (8), vkl m and /k describe the elastic and thermal expansion behaviors of the microscopic base cell, respectively, and can be numerically solved using finite element discretization. The homogenized elastic stiffness CHijkl and thermal expansion coefficient aHkl are obtained as follows:

1 C Hijkl ¼   Vy

 Z  @ vkl C ijkl  C ijmn m dV y @yn y V

1 1 aHij ¼ ðC Hijkl Þ   Vy

  @/ C pqkl akl  k dV y @yl V y

Z

ð10Þ

ð11Þ

  where V y  is the volume of a homogenized microscopic unit cell. As

Fig. 2. Macroscopic and periodic microscopic structures of the nanocomposite and RVE.

shown in Eqs. (10) and (11), both vkl m and /k should be provided to calculate the homogenized thermoelastic properties of nanocomposites. More details on the homogenization method and the derivation of the homogenized properties can be found in Ref. [38]. In the next section, we demonstrate the finite element discretization process used to obtain these two quantities and the homogenized thermoelastic properties.

346

S. Chang et al. / Composite Structures 128 (2015) 342–353

3.2. Finite-element discretization The displacement field and strain can be expressed in terms of the shape functions and nodal displacement for finite-element discretization as follows:

u ¼ ½Nfue g

ð12Þ

ru ¼ ½@ ½Nfue g ¼ ½Bfue g

ð13Þ

where ½N and ½B are the matrices of the shape functions and their first order partial derivatives with respect to the local coordinate system. The tensor vðx; yÞ has the following symmetry due to the periodicity of the base-cell:

vijk ¼ vikj

ð14Þ



v ¼ vi11 ; vi22 ; vi33 ; vi12 ; vi23 ; vi13



ð15Þ

and the tensor /k is expressed as,

/k ¼ ½U:

ð16Þ

By substituting Eqs. (13)–(16) into Eqs. (7) and (8), the following equations are derived:

Z

½BT ½C½B½vdV y ¼

Z

Vy

Z

½BT ½CdV y

ð17Þ

½BT ½CfagdV y

ð18Þ

Vy

½BT ½C½B½UdV y ¼

Z

Vy

Vy

Thus, v and U can be obtained from the following equations:

½K½v ¼ ½Q 

ð19Þ

½K½U ¼ ½P

ð20Þ

where,

½K ¼

Z

½BT ½C½BdV y

Vy

½Q  ¼ ½P ¼

Z Z

½BT ½CdV y

ð21Þ

Vy

½BT ½CfagdV y

Vy

Eqs. (10) and (11) can be expressed in the finite-element discretized form as follows:

h

Z i 1 CH ¼   ½C  ½C½B½vdV y V y Vy





h

aH ¼ CH

i1 1 Z   V y  V ½Cðfag  ½B½UÞdV y y

ð22Þ

ð23Þ

The homogenized elastic stiffness tensor CH is obtained via Eq. (22) using the tensor v obtained from Eq. (19). Then, the homogenized CTE aHkl is obtained via Eq. (23), with the homogenized elastic stiffness tensor CHijkl and tensor u, which is obtained from Eq. (20). 3.3. Numerical procedures to obtain thermoelastic properties of interphase To fully take advantage of the homogenization method to predict the overall thermoelastic properties of a nanocomposite, both the volume fractions and thermoelastic properties of the individual phases that constitute the microscopic RVE of the nanocomposite should be provided. Among the three phases that constitute the present nanocomposite, the volume fraction of the nanoparticle and thermoelastic properties of the nanoparticle and pure matrix

phase can be provided from MD simulations. Thus, the remaining quantities are the volume fractions of the interphase and matrix, and the thermoelastic properties of the interphase. Once the volume fraction of the interphase is determined, the volume fraction of the matrix phase is naturally fixed. Thus, the key of the present multiscale homogenization method for nanocomposites is the determination of the volume fraction and thermoelastic properties of the interphase zone. One acceptable suggestion for determining the volume fraction of the interphase is to use the radial density distribution of the surrounding polymer [5–7]. As a result of the non-bond interaction between the nanoparticle and the surrounding polymer, the structural packing and orientation of the polymer right in the vicinity of the nanoparticle are altered [39,40]. Because the density and crystallinity of the polymer truly affect its properties, it can be reasonably assumed that the highly densified zone in the radial density distribution can be modeled as an independent phase whose thermoelastic properties are different from those of the pure matrix phase. In our previous study, the peak radial density of an epoxy was found to be 3 Å away from the surface of the nanoparticles [10]. Thus, the thickness of the interphase was determined to be 7 Å by adding the outer thickness of the interphase zone of 4 Å to account for the gradual decrease in the density distribution after the peak point. In this study, however, we set the outer thickness of the interphase zone as 2 Å because of the modeling consideration. Because the spatial distribution of the nanoparticles in the microscopic RVE has a SC periodicity with a maximum closed packing ratio of 52.3%, the interphase zone can be overlapped at smaller particle sizes with an outer thickness of 4 Å at a volume fraction of 5.8%. Even if the thickness of the interphase varies, the interphase zone can be modeled as thick and less stiff or thin and stiffer as one wishes. However, the combination of the thickness and thermoelastic properties of the interphase does not affect the estimation of the overall properties of the nanocomposite [38]. Once the thickness of the interphase is determined, the volume fractions of the interphase and matrix can be calculated. Then, the only unknown to be obtained is the thermoelastic properties of the interphase. Because the overall thermoelastic properties of the nanocomposite have already been provided from MD simulations, we can estimate the properties of an interphase zone that can generate the same thermoelastic properties for the nanocomposite as the MD simulation results through the homogenization method. In the case of a micromechanics-based multiscale model, the properties of the interphase can be inversely estimated by manipulating the closed form solution of the linear algebraic micromechanics solution with respect to the properties of the interphase. In the case of the homogenization method, however, such a manipulation is meaningless because the mathematical forms of the homogenized elastic stiffness and CTE given in Eqs. (10) and (11) are general forms regardless of the number, shape, and properties of each subdomain in the microscopic RVE. Therefore, we propose a simple numerical iteration scheme for the inverse estimation of the thermoelasticity of the interphase. At the given thermoelastic properties of the nanocomposite, initial estimations of the elastic moduli and CTE of the interphase are assigned, and the overall thermoelastic properties of the interphase are calculated using the homogenization method. If the resultant properties from the homogenization method are close to the MD simulation results at the same particle radius, then the solution is obtained. Otherwise, a new round of iteration is followed until the convergence criterion is satisfied. The numerical error used in the convergence criterion is the Euclidean norm of the difference between the homogenized properties and the MD simulation results, normalized by the MD simulation results given as,

347

S. Chang et al. / Composite Structures 128 (2015) 342–353

Table 2 Thermoelastic properties of nanocomposites and other constituents obtained from MD simulations.

Fig. 3. Probability density of nanoparticle size following beta distribution. The mean, upper bound, and lower bound of the nanoparticle size were 7.54. 12.39, and 2.42 Å, respectively.

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi !2 !2 u MD H H u E GMD comp  Gcomp ðEint ; Gint Þ comp  Ecomp ðEint ; Gint Þ : error ¼ t þ EMD GMD comp comp ð24Þ MD H H In the above, EMD comp and Ecomp are Young’s moduli; Gcomp and Gcomp are the shear moduli of the nanocomposite obtained from MD simulations and the homogenization method, respectively; and Eint and Gint are the Young’s modulus and shear modulus of the effective interphase – the solution of the iteration. Likewise, the error for the CTE of the effective interphase is defined as,

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi !ffi u u CTEMD  CTEH ðCTEint Þ 2 comp comp t error ¼ CTEMD comp

ð25Þ

H where CTEMD comp and CTEcomp are the CTEs of the nanocomposite obtained from the MD simulations and homogenization method, respectively, and CTEint is the solution of this iteration process, which is the CTE of the effective interphase zone. The error toler-

ance for the convergence was set as 1010 . In estimating the thermoelastic properties of the effective interphase, as has already been stated, the spatial distribution of the nanoparticles was an SC array. As a result, the overall elastic constants obtained from the homogenization method have cubic anisotropy. However, the MD simulation results for the nanocomposite used to obtain the properties of the effective interphase using the homogenization

Radius (Å)

E (GPa)

G (GPa)

CTE (106/K)

ti (Å)

7.54 9.00 10.00 10.90 11.63 Epoxy SiC

4.59 4.65 3.69 3.48 3.56 3.16 451.6

1.71 1.70 1.33 1.29 1.31 1.14 182.5

86.7 91.0 84.2 91.8 96.0 103.6 2.48

5.11 5.85 4.85 4.85 5.72 – –

method are isotropic. Therefore, in the iterative steps using Eq. (23), a direct comparison of the overall elastic moduli of the nanocomposite, as obtained from the homogenization method, with those from the MD simulation is not straightforward. To resolve this problem, we used an additional treatment to generate equivalent isotropic elastic constants from the homogenization method, that is, Hill’s model [41] for polycrystalline aggregates. The basic idea for this special manipulation can be found in Refs. [38,42]. 3.4. RVE modeling for stochastic analysis To account for the uncertainties that may affect the overall properties of nanocomposites, both the particle radius and spatial distribution were considered. Therefore, many-particle RVEs embedding randomly distributed nanoparticles were constructed for the stochastic analysis. In this study, we assumed that all the nanoparticles were uniformly distributed without any noticeable agglomeration that could lead to a merging of the interphase zone. In contrast to the interphase characterization process, which deals with only one nanoparticle in a microscopic domain with periodic boundary conditions, we considered three different sets of manyparticle-embedded RVEs, which contained 8, 16, and 24 nanoparticles in the matrix domain. In each set, a total of 70 different microscopic RVEs were constructed to derive the stochastic variation of the thermoelastic properties from the homogenization method. The volume fraction of the three different sets was fixed at 5.8%. To assign a random distribution of the particle radii inside the RVEs, we adopted the beta-distribution function. In contrast to a normal distribution function, the beta distribution function can allow lower and upper bounds for the possible random numbers to prevent the assignment of negative numbers, which is not desirable when assigning a particle radius. To compare the thermoelastic properties of the many-particle RVEs with the reference values obtained from the MD simulation or micromechanics model at a given nanoparticle size, the mean and maximum radii in the beta distribution function were set at 7.54 and 12.3988 Å, respectively.

Fig. 4. RVE of nanocomposite containing 24 nanoparticles. Left: geometry of nanoparticles and surrounding interphase zone generated using DIGIMATÒ FE, except for the nanoparticles and interphase in the matrix zone. Right: Discretized model with tetrahedral elements.

348

S. Chang et al. / Composite Structures 128 (2015) 342–353

18 Micromechanics Homogenization

Young's modulus(GPa)

16 14 12 10 8 6 4 2 5

6

7

8

9

Particle Radius(

10

11

12

)

(a) Young’s modulus 9 Micromechanics Homogenization

Shear modulus(GPa)

8 7 6

4. Results and discussion

5

4.1. Thermoelastic properties of nanocomposites from MD simulation 4 3 2 1

5

6

7

8

9

Particle Radius(

10

11

12

)

(b) Shear modulus 120 Micromechanics Homogenization 100

CTE(10-6 /K)

with coated effective interphase zones were constructed using a composite RVE modeling software, DIGIMATÒ FE [43]. Fig. 4 shows one of the RVEs of the 24-nanoparticle-embedded nanocomposite before and after the mesh generation. After constructing each RVE, we imported each set of geometrical data into our own inhouse finite element homogenization program and used further numerical steps to determine the stochastic variation of the thermoplastic properties of the nanocomposite. Even if we modeled many-particle RVEs for stochastic homogenization, the maximum of 24 nanoparticles in a unit cell was still not enough to guarantee the statistical isotropy of the unit cell in a mechanical test. Therefore, a huge number of inclusion phases should be modeled in a microscopic unit cell. However, such a situation is not desirable in a stochastic analysis because of the enormous number of computations needed using different microstructures for the nanocomposite to determine the distribution of its overall properties according to the radius of the nanoparticles. Thus, we again used Hill’s model [41] in calculating the overall isotropic elastic moduli of the nanocomposites for the stochastic analysis. Since the cubic arrangement of nanoparticles in a periodic macroscopic domain is broken by embedding many particles in an RVE, the Zener anisotropy ratios for a regular SC arrangement and randomly distributed arrangement of nanoparticles were compared.

80

60

40

20 5

6

7

8

9

Particle Radius(

10

11

12

)

(c) Coefficient of linear thermal expansion Fig. 5. Thermoelastic properties of interphase zone obtained from proposed multiscale homogenization method and micromechanics-based approach.

The elastic moduli and CTEs of the nanocomposites, epoxy, and SiC are listed in Table 2, along with the thickness of the effective interphase zone. As the radius of the embedded SiC nanoparticle decreases, both the Young’s modulus and shear modulus of the nanocomposite increases, while the CTE decreases, demonstrating a clear filler size dependency. This positive reinforcing effect of the smaller nanoparticles can again be elucidated from the effective interphase zone. As listed in Table 2, the thickness of the interphase has finite values that are close to 5 Å. If the thickness of the effective interphase is set as 5 Å, regardless of the nanoparticle size, the volume fraction of the effective interphase zone increases as the size of the nanoparticle decreases. Thus, the relative contribution of the interphase zone to the overall properties of the nanocomposite increases. Another important mechanism for this size effect is the non-bond interaction energy between the nanoparticle and surrounding polymer. According to Yang et al. [40], the magnitude of the non-bond interaction energy between the nanoparticle and matrix normalized by the embedded nanoparticle volume increases as the size of the nanoparticle decreases. The positive reinforcing effect of the nanocomposite can also be interpreted in terms of the hydrostatic pressure effect. When the density of the interphase zone is much higher than the pure matrix zone, as in the case of asymmetric yielding of an amorphous polymer in tension and compression, the atoms in the interphase zone have less space to allow the motion of individual atoms according to the external mechanical or thermal loads. As a result, an interphase zone that has experienced structural densification and ordering by embedding nanoparticles into a polymer matrix naturally has a different elastic stiffness and CTE compared to the pure matrix region. 4.2. Thermoelastic properties of effective interphase zone

Fig. 3 shows the probability density of the particle radius considered in the present stochastic homogenization analysis. Based upon the distribution of the nanoparticle size shown in Fig. 3, the RVEs of nanocomposites containing 8, 16, and 24 nanoparticles

This section presents the elastic moduli and CTE of the effective interphase zone obtained using the present homogenization method with the inverse numerical scheme demonstrated in Section 3.3. For comparison with the micromechanics-based

349

S. Chang et al. / Composite Structures 128 (2015) 342–353 Table 3 Coefficient of thermoelastic properties of effective interphase zone and nanocomposites fitted into function of nanoparticle radius. Property

Homogenization (interphase)

a

b

c

a

b

c

a

b

c

E (GPa) G (GPa)

3.16 1.14 103.6

158 151 58:6

0:48 0:62 0:36

3.16 1.14 103.6

122.8 119.0 56:1

0:41 0:53 0:35

3.57 1.29 96

6.05 2.53 29:96

0:86 0:92 0:46

CTEð106 =KÞ

Micromechanics (interphase)

MD (nanocomposites)

Fig. 6. Comparison of mean elastic moduli of nanocomposites obtained from normal homogenization and stochastic analyses.

multiscale approach, the elastic moduli and CTE of the effective interphase zone were inversely obtained from the multi-inclusion model [44,45] following the linear algebraic form provided in the previous work [46]. While the thickness of the effective interphase zone was set at 7 Å in Ref. [46], we reproduced the thermoelastic properties of an effective interphase zone with a thickness of 5 Å for consistency with the present homogenization approach. The Young’s modulus, shear modulus, and CTE of the interphase zone at each particle radius obtained from the present multiscale homogenization method and the micromechanics-based method are depicted in Fig. 5. As the radius of the nanoparticle increased, both Young’s modulus and the shear modulus of the effective interphase zone gradually decreased. Even though larger nanoparticles were not considered in the MD simulations, it seems reasonable to conjecture that the moduli of the effective interphase zone would gradually converge to the Young’s and shear moduli of the epoxy matrix zone, 3.16 and 1.14 GPa, respectively. Compared with the micromechanics-based bridging model, the elastic moduli obtained from the multiscale homogenization method were smaller, regardless of the nanoparticle size, and the gap between the two interphase properties decreased as the size of the nanoparticle increased. In contrast to the homogenization method, which is based on a rigorous mathematical foundation, a micromechanics model based on the average of the strain or stress has limitations in a rigorous description of the internal field fluctuations, which can affect the overall properties of heterogeneous structures. As revealed by Yang and Cho [6], the description of such a local field fluctuation becomes critical in nanocomposites because of the interphase zone. For example, at a high-volume-fraction condition, the Mori–Tanaka (M–T) model [47,48] can underestimate the effective properties of composites, and the self-consistent model [49,50] may show better performance because it assumes that

the average stiffness of the matrix zone is the effective stiffness of the composites itself, which is stiffer than that of the matrix phase in general. Referring to the conclusion on the non-dilute concentration effect of a nanocomposite by Yang and Cho [9], the application of micromechanics-based approaches to predict the properties of the effective interphase zone can encounter an overestimation problem, which is analogous to the underestimation problem in predicting the overall properties of composites. However, the difference between the present homogenization method and the micromechanics-based approach is not critical in this case. In the case of a high-volume-fraction condition for nanoparticles above 10%, however, such a difference can be magnified and result in a non-negligible error in reproducing the size dependent properties of nanocomposites using continuum-based approaches such as a finite element analysis to validate the accuracy of the interphase properties. Similar to the variation in the CTE of nanocomposites, the CTE of the effective interphase zone increases as the radius of the nanoparticle increases, as shown in Fig. 5. In this case, however, the difference between the micromechanics-based model and the present homogenization method is almost negligible. As the variation in the CTE of the effective interphase is also asymptotic, we may fit the thermoelastic properties of the interphase as a continuous function of the particle radius to utilize the homogenization method at various particle sizes. Even though all the particle radii were not considered in the present MD simulations, such an approach can guarantee an enormous computational efficiency, especially for a stochastic analysis and repeated calculations for the optimum design. In the next section, by taking advantage of the description of the interphase properties as a function of the particle size, the mean and variation of the effective thermoelastic properties of a nanocomposite

350

S. Chang et al. / Composite Structures 128 (2015) 342–353

determine using the stochastic homogenization method will be demonstrated. 4.3. Stochastic variation of thermoelastic properties of nanocomposites Once the elastic moduli and CTE of the interphase zone are obtained, to account for the size effect induced by the interphase zone in the stochastic homogenization method or micromechanics model, they can be fitted into a continuous function of the nanoparticle size [10,16]. The thermoelastic properties of the effective interphase zone surrounding each nanoparticle are assigned according to the nanoparticle size in many-particle RVEs. Using a simple least square approximation, the thermoelastic properties of the effective interphase zone shown in Fig. 5 are fitted into a function of the nanoparticle size, as PropertyðE; G; CTEÞ ¼ a þ bexpcrp , where rp is the nanoparticle radius. The three coefficients are listed in Table 3. Among the three coefficients, a is the converging point and coincides with the thermoelastic properties of the pure epoxy. Therefore, all the curve-fitted properties of the interphase zone gradually converge to those of the epoxy zone as the size of particle increases to describe the gradual decrease in the size dependent properties of the nanocomposite. In addition to the properties of the effective interphase, the moduli and CTEs of nanocomposites obtained from MD simulations can also be expressed in terms of the nanoparticle radius using the same least square approximation. These coefficients are also listed in Table 3. Therefore, for the nanoparticle size distribution shown in Fig. 3, a simple stochastic analysis directly using the least square-fitted thermoelastic properties of the nanocomposites is also an available approach. For the micromechanicsbased and MD simulation-based stochastic approaches, the elastic moduli and CTEs of 560 monoparticle unit cell structures were considered. Since the micromechanics model represents random heterogeneous structures through an equivalent mono-inclusion and deals with the mean internal field quantities in a constitutive relation, the spatial distribution of the nanoparticles was generally not considered. Such a restriction was also applied to the MD simulation-based stochastic analysis. A comparison of the mean elastic properties of the nanocomposites obtained from regular homogenization with a mono-particle periodic RVE and the stochastic analysis based on micromechanics, MD, and many-particle periodic RVEs is shown in Fig. 6. For the regular homogenization method, the nanoparticle radius was 7.54 Å. Since the Zener anisotropy ratio of the homogenization results for the mono-particle and many-particle RVEs was still not zero, the moduli shown in Fig. 6 are orientation-averaged values using Hill’s model. One thing noteworthy in Fig. 6 is that both the Young’s modulus and shear modulus obtained from the stochastic homogenization method are larger than those obtained from regular homogenization and other stochastic analyses based on the micromechanics and MD simulation results. Since the latter three methods do not consider the random spatial distribution of nanoparticles, the consequent local non-dilute interaction between individual particles and the effective interphases are not taken into account. Therefore, the elastic moduli are underestimated compared with the stochastic homogenization method. Even though the maximum of 24 nanoparticles in an RVE was still not fully representative of real nanocomposites microstructure, the effect of a random spatial distribution of embedded nanoparticles on the overall elastic properties could be successfully identified. Fig. 7 shows the probability density of the elastic moduli and CTEs of the nanocomposites obtained using each stochastic approach. In the MD- and micromechanics-based results, the probability density distribution has a wider deviation, while the deviation of the stochastic homogenization method is narrow and gradually decreases as the number of nanoparticles in the

Fig. 7. Probability densities of thermoelastic properties of the nanocomposites obtained using stochastic approaches.

RVE increases. Moreover, the mean of the properties is not sensitive to the number of nanoparticles, while the probability density becomes clearer as the number of nanoparticles increases. To study how the anisotropy ratio is changed by the random distribution of nanoparticles, Young’s modulus in the X direction ðE1 Þ and the shear modulus in the XY plane ðG12 Þ were obtained from 8-particle and mono-particle RVEs, and are compared in

351

S. Chang et al. / Composite Structures 128 (2015) 342–353

Fig. 8. To trace the change in the elastic moduli of the nanocomposites by the random distribution of (1) only the particle radius, (2) only the spatial location, and (3) both parameters, the probability densities of the moduli were compared with those for the moduli of the mono-particle RVE with SC symmetry. When a random distribution of the nanoparticle radius was considered, the mean of the elastic modulus was not prominently changed, and the probability density showed a wide deviation. On the other hand, the Young’s modulus of the 8-particle RVEs decreased and the shear moduli increased when the random spatial distribution of nanoparticles was considered, with a fixed particle radius. When random distributions of both the particle radius and spatial distribution were considered, both moduli were further changed by more than 1 GPa. This result indicates that a random spatial distribution of nanoparticles in a nanocomposite is an important parameter in determining its effective properties. Since an SC array of periodic microstructures in the typical homogenization method is a simplified assumption and has an intrinsic cubic anisotropy, a consideration of many random particles in the size and distribution is more desirable to determine the realistic behavior of a nanocomposite.

The necessity of the stochastic approach is again confirmed from the probability density of the Zener anisotropy ratio defined as,

A¼

2C 44 C 11  C 12

ð26Þ

where C ij is the stiffness matrix components. When the stiffness matrix of a system is fully isotropic, the Zener anisotropy becomes one. The probability density of the Zener anisotropy ratio of the 8particle RVE is shown in Fig. 9, for random distributions of the nanoparticle size and spatial location. For the mono-particle RVE, the Zener anisotropy ratio is 0.7816. When only a random

(a)

(a) (b)

(c) (b) Fig. 8. Changes in elastic modulus of nanocomposite by random particle size and spatial distribution.

Fig. 9. Changes in probability density of Zener anisotropy ration according to random microstructural parameters: (a) only nanoparticle radius varies, (b) only spatial distribution of nanoparticles varies, and (c) both radius and spatial distribution of nanoparticles vary.

352

S. Chang et al. / Composite Structures 128 (2015) 342–353

nanoparticle size distribution is considered, the most probable anisotropy ratio is close to 0.8. On the other hand, the most probable anisotropy ratio becomes larger than one when a random spatial distribution of nanoparticles is considered. The change in the Zener anisotropy ratio shown in Fig. 9 correlates well with the change in the elastic moduli according to the random microstructural parameters shown in Fig. 8. Even if the Zener anisotropy ratio of the 8-particle RVE is still not one, it is obvious that RVEs with eight or more nanoparticles are more realistic than a mono-particle RVE with a periodic SC array. Since the difference in the average overall thermoelastic properties of nanocomposites between mono-particle and many-particle RVEs is not negligible, a random spatial distribution of nanoparticles should be rigorously taken into account to predict the behavior of nanocomposites. 5. Conclusion In this study, a stochastic multiscale homogenization method was provided to determine the thermoelastic properties of nanoparticle composites. Based on MD simulations, the nanoparticle size-dependent elastic moduli and CTEs of nanocomposites were obtained. To describe the nanoparticle size effect in an equivalent continuum model, a three-phase microstructure incorporating the effective interphase zone was considered as a periodic microstructure in the finite element homogenization method. By equating the MD simulation results to the homogenization solution, the thermoelastic properties of the effective interphase zone were numerically obtained and compared with the results using a multi-inclusion model. To describe the nanoparticle sizedependent behavior over a wide range of nanoparticle sizes, the thermoelastic properties of the effective interphase zone were described as a function of the particle radius. Using the properties of the effective interphase zone and a beta distribution of the particle radius, a stochastic analysis was performed to observe the variation of the overall thermoelastic properties of nanocomposites using many-particle RVEs. Compared with the regular homogenization method with a mono-particle unit cell, the stochastic homogenization approach could better describe the local non-dilute concentration effect arising from locally close-packed nanoparticles, and thus predicted enhanced thermoelastic properties. In view of the computational cost and availability, a micromechanics-based approach or regular homogenization approach with a mono-particle periodic unit cell is truly superior to the stochastic homogenization approach. Unless the spatial distribution of embedded nanoparticles was considered, however, it was found that nonnegligible error in the predicted properties could occur. Since the random spatial distribution is an important microstructural parameter, a consideration of severe particle agglomeration and the consequent overlap of the effective interphase zone is an open issue in stochastic approaches. For a rigorous consideration of the degradation of the effective properties of nanocomposites by the agglomeration, molecular-level observation of the alteration in the density and crystallization of the effective interphase zone is of primary concern. To provide a more robust and versatile microstructural design tool, we will extend the current multiscale stochastic homogenization approach to the nanoparticle agglomeration problem. Acknowledgement This work was supported by a National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIP) (No. 2012R1A3A2048841).

References [1] Thostenson ET, Li C, Chou TW. Nanocomposites in context. Compos Sci Technol 2005;65(3-4):491–516. [2] Fu SY, Feng XQ, Lauke B, Mai YW. Effects of particle size, particle/matrix interface adhesion and particle loading on mechanical properties of particulate-polymer composites. Compos Part B 2008;39:933–61. [3] Koo J. Polymer nanocomposites: processing, characterization, and applications. In: McGraw-Hill nanoscience and technology series. McGraw-Hill; 2006. [4] Yan LT, Xie XM. Computational modeling and simulation of nanoparticle selfassembly in polymeric systems: structures, properties and external field effects. Prog Polym Sci 2013;38(2):369–405. [5] Yang S, Cho M. Scale bridging method to characterize mechanical properties of nanoparticle/polymer nanocomposites. Appl Phys Lett 2008;93(4):043111. [6] Yang S, Cho M. A scale-bridging method for nanoparticulate polymer nanocomposites and their nondilute concentration effect. Appl Phys Lett 2009;94(22):223104. [7] Yu S, Yang S, Cho M. Multi-scale modeling of cross-linked epoxy nanocomposites. Polymer 2009;50(3):945–52. [8] Adnan A, Sun CT, Mahfuz H. A molecular dynamics simulation study to investigate the effect of filler size on elastic properties of polymer nanocomposites. Compos Sci Technol 2007;67(3):348–56. [9] Choi J, Yu S, Yang S, Cho M. The glass transition and thermoelastic behavior of epoxy-based nanocomposites: a molecular dynamics study. Polymer 2011;52(22):5197–203. [10] Yu S, Yang S, Cho M. Multiscale modeling of cross-linked epoxy nanocomposites to characterize the effect of particle size on thermal conductivity. J Appl Phys 2011;110(12):124302. [11] Yang S, Choi J, Cho M. Intrinsic defect-induced tailoring of interfacial shear strength in CNT/polymer nanocomposites. Compos Struct 2015;127:108–19. [12] Goyal RK, Tiwari AN, Mulik UP, Negi YS. Thermal expansion behavior of high performance PEEK matrix composites. J Phys D Appl Phys 2008;41:085403. [13] Sun Y, Zhang Z, Wong CP. Study and characterization on the nanocomposite underfill for flip chip applications. IEEE Trans Compon Packag Technol 2006;29(1):190–7. [14] Goertzen WK, Kessler MR. Thermal expansion of fumed silica/cyanate ester nanocomposites. J Appl Polym Sci 2008;109(1):647–53. [15] Jang JS, Bouveret B, Suhr J, Gibson RF. Combined numerical/experimental investigation of particle diameter and interphase effects on coefficient of thermal expansion and Young’s modulus of SiO2/epoxy nanocomposites. Polym Compos 2012;33(8):1415–23. [16] Choi J, Yang S, Yu S, Shin H, Cho M. Method of scale bridging for thermoelasticity of cross-linked epoxy/SiC nanocomposites at a wide range of temperatures. Polymer 2012;53(22):5178–89. [17] Eshelby JD. The determination of the elastic field of an ellipsoidal inclusion, and related problems. Proc R Soc London Ser A 1957;241(1226):376–96. [18] Longo A, Carotenuto G, Palomba M, Nicolar SD. Dependence of optical and microstructure properties of thiol-capped silver nanoparticles embedded in polymeric matrix. Polymer 2011;3(4):1794–804. [19] Gam S, Meth JS, Zane SG, Chi C, Wood BA, Winey KI, et al. Polymer diffusion in a polymer nanocomposite: effect of nanoparticle size and polydispersity. Soft Matter 2012;8:6512–20. [20] Silani M, Talebi H, Ziaei-Rad S, Kerfriden P, Bordas SP, Rabczuk T. Stochastic modelling of clay/epoxy nanocomposites. Compos Struct 2014;183:241–9. [21] Accelrys Inc. San Francisco. . [22] Sun H. COMPASS: an ab initio force-field optimized for condensed-phase applications-overview with details on alkane and benzene compounds. J Phys Chem B 1998;102:7338–64. [23] Komarov PV, Yu-Tsung C, Shih-Ming C, Khalatur PG, Reineker P. Highly crosslinked epoxy resins: an atomistic molecular dynamics simulation combined with a mapping/reverse mapping procedure. Macromolecules 2007;40:8104–13. [24] Yarovsky I, Evans E. Computer simulation of structure and properties of crosslinked polymers: application to epoxy resins. Polymer 2002;43:963–9. [25] Varshney V, Patnaik SS, Roy A, Farmer BL. A molecular dynamics study of epoxy-based networks: cross-linking procedure and prediction of molecular and material properties. Macromolecules 2008;41(18):6837–42. [26] Li C, Strachan A. Molecular simulations of crosslinking process of thermosetting polymers. Polymer 2010;51(25):6058–70. [27] Bandyopadhyay A, Valavala PK, Clancy TC, Wise KE, Odegard GM. Molecular modeling of crosslinked epoxy polymers: the effect of crosslink density on thermomechanical properties. Polymer 2011;52(11):2445–52. [28] Li C, Browning A, Christensen S, Strachan A. Atomistic simulations on multilayer graphene reinforced epoxy composites. Compos Part A 2012;43:1293–300. [29] Shin H, Yang S, Chang S, Yu S, Cho M. Multiscale homogenization modeling for thermal transport properties of polymer nanocomposites with Kapitza thermal resistance. Polymer 2013;5(28):1543–54. [30] Andersen HC. Molecular dynamics simulations at constant pressure and/or temperature. J Chem Phys 1980;72:2384–93. [31] Berendsen HJC, Postma JPM, van Gunsteren WF, DiNola A, Haak JR. Molecular dynamics with coupling to an external bath. J Chem Phys 1984;81:3684–90. [32] Parrinello M, Rahman A. Strain fluctuations and elastic constants. J Chem Phys 1982;76(5):2662–6.

S. Chang et al. / Composite Structures 128 (2015) 342–353 [33] Parrinello M, Rahman A. Crystal structure and pair potentials: a molecular dynamics study. Phys Rev Lett 1980;45:1196–9. [34] Sanchez-Palencia E. Non-homogeneous media and vibration theory. Lecture notes in physics, vol. 127. Berlin: Springer; 1981. [35] Bakhvalov N, Panasenko G. Homogenization: averaging processes in periodic media. New York: Kluwer; 1984. [36] Bendsøe MP, Kikuchi M. Generating optimal topologies in structural design using a homogenization method. Comput Methods Appl Mech Eng 1988;71:197–224. [37] Guedes JM, Kikuchi N. Preprocessing and postprocessing for materials based on the homogenization method with adaptive finite element methods. Comput Methods Appl Mech Eng 1990;83:143–98. [38] Cho M, Yang S, Chang S, Yu S. A study on the prediction of the mechanical properties of nanoparticulate composites using the homogenization method with the effective interface concept. Int J Numer Methods Eng 2011;85(12):1564–83. [39] Harton SE, Kumar SK, Yang H, Koga T, Hicks K, Lee H, et al. Immobilized polymer layers on spherical nanoparticles. Macromolecules 2010;43(7):3415–21. [40] Yang S, Choi J, Cho M. Elastic stiffness and filler size effect of covalently grafted nanosilica polyimide composites: molecular dynamics study. ACS Appl Mater Interfaces 2012;4(9):4792–9.

353

[41] Hill R. The elastic behavior of a crystalline aggregate. Proc R Soc London A 1952;65(5):349–54. [42] Friebel C, Doghri I, Legat V. General mean-field homogenization schemes for viscoelastic composites containing multiple phases of coated inclusions. Int J Solids Struct 2006;43:2513–41. [43] e-Xstream engineering SA. MSC Software Corporation. . [44] Hori M, Nemat-Nasser S. Double-inclusion model and overall moduli of multiphase composites. Mech Mater 1993;14:189–206. [45] Li JY. Thermoelastic behavior of composites with functionally graded interphase: a multi-inclusion model. Int J Solids Struct 2000;37(39):5579–97. [46] Yang S, Cho SYuM. Sequential thermoelastic multiscale analysis of nanoparticulate composites. J Appl Phys 2010;108(5):056102. [47] Benveniste Y. A new approach to the application of Mori–Tanaka’s theory in composites materials. Mech Mater 1987;6:147–57. [48] Weng GJ. The theoretical connection between Mori–Tanaka’s theory and the Hashin–Shtrikman–Walpole bounds. Int J Eng Sci 1990;28:1111–20. [49] Hill R. Elastic properties of reinforced solids: some theoretical principles. J Mech Phys Solids 1963;11:357–72. [50] Hill R. A self-consistent mechanics of composites materials. J Mech Phys Solids 1965;13:213–22.