polymer nanocomposites considering the interphase and interfacial imperfection

International Journal of Plasticity 41 (2013) 124–146

Contents lists available at SciVerse ScienceDirect

International Journal of Plasticity journal homepage: www.elsevier.com/locate/ijplas

Nonlinear multiscale modeling approach to characterize elastoplastic behavior of CNT/polymer nanocomposites considering the interphase and interfacial imperfection Seunghwa Yang a, Suyoung Yu a, Junghyun Ryu a, Jeong-Min Cho b, Woomin Kyoung b, Do-Suck Han b, Maenghyo Cho a,⇑ a WCU Multiscale Mechanical Design Division, School of Mechanical and Aerospace Engineering, Seoul National University, San 56-1, Shillim-Dong, Kwanak-Ku, Seoul 151-742, South Korea b Computer Aided Engineering and Materials Research Team, Central Advanced Research and Engineering Institute, Hyundai and Kia Corporate Research and Development Division, 460-30, Sam-Dong, Uiwang-Si, Gyeonggi-Do 437-040, South Korea

a r t i c l e

i n f o

Article history: Received 5 June 2012 Received in final revised form 10 September 2012 Available online 24 September 2012 Keywords: Nanocomposites Carbon nanotube Hierarchical multiscale model Interphase Weakened interface

a b s t r a c t A hierarchical multiscale modeling approach to characterize the elastic and plastic behavior of carbon nanotube (CNT)-reinforced polymer nanocomposites is proposed via molecular dynamics simulations and a continuum nonlinear micromechanics based on the secant moduli method. Even though the importance of the densified interphase zone formed between the CNT and polymer matrix has been demonstrated by many related studies for elastic properties, studies on how to identify the behavior and contribution of the interfacial condition and interphase zone in the overall elastoplastic behavior of nanocomposites is still an open issue. Different from conventional micromechanics approaches that homogenize overall elastoplastic behavior of heterogeneous structures from known behaviors of its constituent phases, the present study focuses on the identification of local elastoplastic behavior of the interphase region from the known elastoplastic behavior of nanocomposites through a hierarchical domain decomposition method. Firstly, the overall elastoplastic behavior of the CNT-reinforced nanocomposites is obtained from molecular dynamics (MD) simulations which are based on an ab initio force field. Due to a weak van der Waals interaction between the pristine CNT and the matrix polymer, the elastoplastic behavior of the nanocomposites clearly shows a weakened interface condition, while the matrix molecular structure in the vicinity of the CNT confirms the existence of the interphase zone. In upper level analysis, an effective matrix concept is adopted, and its elastoplastic behavior is inversely identified by equating the MD simulation result to a two-phase nonlinear micromechanics model that can consider imperfect interfacial condition. Then, the effective matrix domain is again decomposed into the interphase and pure matrix polymer regions in lower level analysis, and the elastoplastic behavior of the interphase is again identified through the same method. Using the constitutive relation of the interphase obtained from the proposed multiscale model, the overall elastoplastic behavior of the nanocomposites is obtained and compared with some available experimental results and an additional MD simulation result to validate the applicability and physical rigorousness of the proposed nonlinear multiscale approach. Crown Copyright Ó 2012 Published by Elsevier Ltd. All rights reserved.

⇑ Corresponding author. Address: WCU Multiscale Mechanical Design Division, School of Mechanical and Aerospace Engineering, Seoul National University, 599 Kwanak-Ro, Kwanak-Ku, Seoul, South Korea. Tel.: +82 2 880 1693; fax: +82 2 886 1693. E-mail address: [email protected] (M. Cho). 0749-6419/$ - see front matter Crown Copyright Ó 2012 Published by Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ijplas.2012.09.010

S. Yang et al. / International Journal of Plasticity 41 (2013) 124–146

125

1. Introduction Among the huge number of applications of carbon nanotubes (CNTs) (Ijjima, 1991), polymer nanocomposites that address CNTs as reinforcing fibers have received remarkable attention over the last two decades since their first development in 1994 (Ajayan et al.,1994; Thostenson et al., 2005; Spitalsky et al., 2010; Jia et al., 2011). Due to the excellent and attractive properties, carbon nanotubes have been widely applied as new rivaling alternatives to conventional reinforcing fillers in the form of independent reinforcing phases (Jin et al., 1998; Vaisman et al., 2006) or in the form of functional auxiliary phases that enhance the overall performance of conventional composites such as fuzzy fiber reinforced composites (Mathur et al., 2008). Despite their potential applicability to various fields to which multifunctional design of materials is required, there have been several challenges pointed out in regard to manifesting the efficient reinforcing effect of the CNTs, such as bundling of individual CNTs as a result of strong inter-tubular adhesion (Gojny et al., 2004; Barai and Weng, 2011a; Yang et al., 2012a; Shi et al., 2004), waviness of the CNTs in dispersion process (Qian et al., 2000; Fisher et al., 2003; Shi et al., 2004; Omidi et al., 2010), etc. These challenging tasks commonly aim at the realization of microstructure where straight CNTs are well dispersed and strongly bonded to the matrix to achieve superb load transfer between the matrix and CNTs. While these issues are important, attention has been focused on the identification, characterization and modeling of the interphase where the structural arrangement of matrix polymer in the vicinity of the CNTs is altered by the interaction with the CNTs, through observation of the densification and crystallization from molecular dynamics simulations (Wei et al., 2004; Wei, 2007; Han and Elliott, 2007), and experimental characterization (Bhattacharyya et al., 2003; Li et al., 2007). In computational modeling of CNT/polymer nanocomposites, therefore, various attempts based on diverse computational schemes have been carried out to account for the effect of interphase on overall multifunctional properties of CNT-reinforced polymer nanocomposites, for the elastic modulus (Odegard et al., 2003; Seidel and Lagoudas, 2006; Yang et al., 2012b), viscoelastic properties (Fisher and Brinson, 2001), thermal conductivity (Seidel and Lagoudas, 2008), and electrical conductivity (Seidel and Lagoudas, 2009). The importance of the interphase is not a restricted theme to the CNT-reinforced composites, and the issue is commonly shared by other types of composites that employ nano-sized fillers (Yang and Cho, 2008, 2009; Yu et al., 2009; Mesbah et al., 2009). Even if the consideration of the interphase in computational modeling of nanocomposites is currently straightforward, it is still an open issue as to how to describe its unique behaviors. Most of the approaches that consider the interphase have arbitrarily set the properties of interphase as an intermediate value between those of the matrix and filler or as a gradient from the properties of the filler to that of the matrix for qualitative parametric studies. Recently, inverse multiscale modeling approaches have been proposed and applied to various properties of nanocomposites by combining molecular dynamics (MD) simulations and continuum modeling approaches such as micromechanics and finite element (FE)-based homogenization methods for elastic, thermoelastic, and thermal conductivity (Yang et al., 2010; Cho et al., 2011; Yu et al., 2011). In these approaches, different from conventional micromechanics approaches that predict the overall properties of composites from known properties of individual phase, the overall properties of the nanocomposites are given by MD simulations, and the corresponding properties of the interphase are inversely obtained from micromechanics constitutive models or numerical calculation via finite element modeling. In addition to the elastic behavior, the yielding and plastic behavior can also be greatly enhanced by introducing CNTs into a polymer matrix (Granb et al., 2008; Bao and Tjong, 2008). Based on the Mori–Tanaka mean field approach, a nonlinear micromechanics approach that can consider the CNT size effect and clustering effect has been recently proposed (Barai and Weng, 2011a). By considering the CNT clustering effect via a two-scale homogenization scheme, the effect of the clustering on overall elastoplastic behavior of nanocomposites has been successfully identified and compared with experimental results. However, for relatively well dispersed polymer nanocomposites, where the interphase has a key responsibility in plastic behavior, the consideration and identification of the plastic behavior of interphase is critically required in equivalent continuum modeling. As has been recently pointed out, the effect of the interfacial elasticity and stress play important roles in the yielding and plastic behavior of nanocomposites (Zhang et al., 2010). Thus, the physical behavior of the interphase should be accounted for in the polymer-based nanocomposites that undergo plastic deformation. In contrast to the interphase, the inter-atomic strength between a pristine CNT and polymer matrix is disappointing and weak in nature unless the CNTs are chemically functionalized or covalently grafted to the matrix polymer (Sinnott, 2002; Qian et al., 2000; Watts and Hsu, 2003). Evidently, the weak interfacial strength guaranteed from experimental observation may be hindered by other candidate factors such as the dispersion state, which also affects the resultant properties of the nanocomposites. Recently, clearer evidence of the weakened interface condition has been demonstrated in molecular dynamics simulations, in which all possible noises that affect the degradation of the properties can be eliminated, and the inconsistencies of the perfect bonding condition in continuum micromechanics modeling of the pristine CNT reinforced composites has been pointed out (Yang et al., 2012b). For the sake of continuum level description of the imperfect interface condition, the relation between the interfacial compliance and the overall elastic and elastoplastic behaviors have been systematically demonstrated (Barai and Weng, 2011a; Esteva and Spanos, 2009) via linear spring model (Qu, 1993). Even if various types of surface treatments are currently being applied to the manufacture of CNT-reinforced composites for efficient load transfer and good dispersion of the fillers, interfacial imperfection is still inevitable and commonly occurs during the manufacturing process and operating conditions. Thus, there is truly a trade-off between the interphase zone and

126

S. Yang et al. / International Journal of Plasticity 41 (2013) 124–146

weak interfacial condition to decide the performance of the CNT-reinforced polymer nanocomposites in elastic and plastic deformation states, especially for reinforcement with untreated CNTs. Other factors aside, these two characteristics contribute to the overall plastic behavior of nanocomposites in a complicatedly coupled manner that makes it difficult to identify the effect of the two factors from the resultant elastic and plastic behavior of the nanocomposites. To assist the understanding of these two factors in linear elastic behavior, molecular dynamics (MD) simulations have ever been preferred (Yang and Cho, 2008). The MD simulation is a truly powerful computational experiment tool for qualitative and quantitative description of the physical behavior of polymer-based materials (Rottler and Robbins, 2001). Despite the computational capacity of current computers, however, there are still limitations in applying MD simulations to large length scale and long time scale problems such as plastic or viscoelastic behavior of polymeric systems. An alternative to resolve the computational efficiency and accuracy problems is the hierarchical multiscale modeling approach by reflecting nanoscale specific phenomena obtained from MD simulations to the continuum-based theory with some modifications (Yang and Cho, 2008). Motivated by the above requirements, the present study proposes a nonlinear multiscale modeling approach to predict the elastic and plastic behavior of polymer nanocomposites, and an inverse computational scheme to identify the elastoplastic behavior of the interphase by hierarchically combining the molecular dynamics simulations and a nonlinear micromechanics model based on the secant moduli method and field fluctuation method. The overall stress–strain relation of the CNT-reinforced composites in elastic and plastic states is basically obtained from MD simulations. In a nonlinear micromechanics model, the Mori–Tanaka mean field approach is used and modified to account for the imperfect interface condition adopting the effective fiber method. From the given stress–strain relation obtained from the MD simulations, the equivalent nanocomposites’ representative volume element (RVE) is sequentially decomposed into an effective matrix and the CNTs in upper level analysis, and the effective matrix domain is again decomposed into the interphase and the pure matrix phase in lower level analysis to identify the elastoplastic behavior of the interphase based on the framework of the nonlinear micromechanics model. Then the problem is turned back to the estimation of the overall elastoplastic behavior of the CNT-reinforced composites from the known behavior of the CNT, interphase, and matrix, and some validation studies for the accuracy and applicability of the present approach are accompanied.

2. Molecular model definition and the proposed overall analytical scheme Prior to the details of the MD simulation procedures and the nonlinear constitutive relations, the configuration of the molecular unit cell that represents the CNT-reinforced nanocomposites, and the overall multiscale analysis scheme that the present study focuses on is presented. Fig. 1(a) is the SEM (Scanning Electron Microscope) fractograph of polypropylene/multi-walled nanotube (MWNT) composites showing relatively well dispersed configuration. For computational resource and efficiency, consideration of such a large amount of CNTs in MD simulation is impossible. Thus, it has been preferred in MD simulation of polymer nanocomposites to model the RVE as a cubic or rectangular unit cell structure that consists of a single inclusion phase and surrounding matrix molecules with a periodic boundary condition. The molecular structure of the single-walled nanotube (SWNT) reinforced polypropylene composites considered in this study is shown in Fig. 2(b). Thanks to the periodic boundary condition, the arrangement of the SWNT in the present molecular model is periodic and regarded as a well dispersed one. In addition, the shape of the SWNT embedded into the matrix is straight. Thus, in MD simulations, the effects of the local aggregation and the waviness of the CNT are naturally excluded. These features of the molecular model are like two sides of a coin; such morphological noises are prevented while consideration of such factors is also limited unless huge numbers of atom are considered in the modeling process. As has been demonstrated in previous equivalent continuum modeling approaches (Barai and Weng, 2011a; Fisher et al., 2003), the

Fig. 1. (a) SEM image of the CNT/polypropylene nanocomposites (Bao and Tjong, 2008). The rectangular brackets represent the CNTs and the surrounding matrix considered in MD simulations. (b) RVE of nanocomposites molecular structure considered in the present MD simulations.

S. Yang et al. / International Journal of Plasticity 41 (2013) 124–146

127

(a)

(b) Fig. 2. Schematic diagrams of the proposed hierarchical two-step multiscale model to identify the elastoplastic behavior of interphase. (a) Upper level analysis to decompose the nanocomposites microstructure into the CNT and the effective matrix phase. (b) Lower level analysis to decompose the effective matrix into the interphase and pure matrix phase.

agglomeration and the waviness of the CNTs and the resultant degradation of the overall properties of nanocomposites can be effectively handled in the frame work of micromechanics approaches. In nonlinear micromechanics modeling approaches to identify the elastoplastic behavior of the interphase, the overall stress–strain relation of the nanocomposites obtained from MD simulations is firstly equated to the two-phase micromechanics model prediction. Here, we adopt the concept of an effective matrix that contains the pure matrix phase and the interphase. The basic idea of the effective matrix is a quite simple but reasonable approach to account for the structural change that forms the interphase. As the interphase region is uniquely formed by adding the CNTs into the matrix, it can be reasonably assumed that the CNT is embedded into the matrix which has already experienced structural change. Thus, the matrix phase that has several sites of structural change inside of its domain is defined as an effective matrix that has different physical behavior from pure matrix phase. When the CNTs are randomly oriented inside of the nanocomposites, the orientation of the interphase region is consequently random, thus, the effective matrix at this configuration becomes isotropic. In this process, the elastoplastic behavior of the effective matrix is set as an unknown variable and inversely obtained from the micromechanics model. For convenience, this process is simply referred to as the ‘‘Upper level analysis’’ as shown in Fig. 2(a). In upper level analysis, the volume concentration is defined as:

fp þ fem ¼ 1;

ð1Þ

where f is the volume fraction and the subscripts ‘p’ and ‘em’ respectively indicate the CNT and effective matrix phase. After obtaining the elastoplastic behavior of the effective matrix, the effective matrix phase is again decomposed into the pure matrix and interphase phase via the same nonlinear micromechanics modeling approach as shown in Fig. 2(b). In lower level analysis, the volume fraction of the pure matrix and the interphase are denoted as fmem and fiem , respectively. Respectively denoting the volume fractions of the interphase and pure PP matrix over the whole range of nanocomposites domain as fi and fm, the following relations are established:

128

S. Yang et al. / International Journal of Plasticity 41 (2013) 124–146

fem ¼ fi þ fm ;

fiem ¼

fi ; fem

fmem ¼

fm : fem

ð2Þ

At this stage, the unknown variable is the elastoplastic behavior of the interphase, and the known elastoplastic behavior of the effective matrix can be equivalently regarded as the overall elastoplastic behavior of two-phase composites where the interphase is embedded into the pure matrix as a reinforcement. As the effective matrix phase is a sub phase that constitutes the nanocomposites, this process is referred as the ‘‘Lower level analysis’’. Once the elastoplastic behavior of the interphase is obtained, the overall elastoplastic behavior at various volume fractionations can be easily estimated from the same nonlinear micromechanics model. Like the inverse process demonstrated above, the overall process is hierarchical as well, and the two-step homogenization process is required. The elastoplastic behavior of the effective matrix is firstly obtained from the known behavior of the matrix and the interphase in lower level analysis. Then, the overall elastoplastic behavior of the nanocomposites is obtained in upper level analysis using the elastoplastic behavior of the effective matrix calculated from the lower level analysis and the elastic constants of the CNT. 3. Molecular dynamics simulations 3.1. Construction of the unit cell In MD simulations, a zigzag type CNT with a chiral vector of (15, 0) is chosen as reinforcement and an amorphous isotatic polypropylene (PP) is used as base polymer. Neither surface functionalization, defects of CNT nor direct covalent grafting to the PP matrix are considered in MD simulations, thus, the CNT in PP matrix is pristine. One single chain of the amorphous PP consists of twenty unit monomers, and the number of PP chains in the unit cell varies according to the volume fraction of the CNT. In total, three different molecular unit cells are considered and have different volume fractions. The details of the cell composition at equilibrated states are arranged in Table 1. As a reference to compare the reinforcing effect of the CNTs, a pure amorphous PP unit cell is also constructed using the same PP chain. The procedures to construct the nanocomposites unit cell shown in Fig. 1(b) is as follows. At first, pure PP unit cells, whose corresponding number of unit monomers is listed in Table 1, are constructed as amorphous structures at a target density of 0.8 g/cm3. Then the total potential energy of the PP unit cell is minimized via the conjugate gradient method with a target convergence threshold for the specified maximum energy change of 0.001 kcal/mol. Then one CNT is embedded at the center of the unit cell aligned to the X direction as is depicted in Fig. 1(b). In this process, the position of some monomers at the center of the unit cell is modified to the Y and Z directions to make a vacancy to insert the CNT. Then the unit cells are minimized again using the same convergence threshold to settle the local energy instability possibly generated by the position modification process of the matrix molecules and to optimize the molecular arrangement surrounding the CNT. All the cell construction and the subsequent energy minimization processes are performed using a commercial molecular dynamics simulation package, Material Studio 5.0 (Accelrys Inc), with the polymer consistent force field (PCFF) (Sun et al., 1994) to describe inter- and intra-atomic interactions. As the PCFF force field does not describes the reactive bond order change, i.e. bond swap, bond break and bond nucleation is not allowed, such phenomena are not described using this force field. Thus, the mechanical behavior of the PP matrix and nanocomposites is simulated based on the intramolecular stretch, bent, distortion, and intermolecular slippage and disentanglement described by the PCFF forcefield. 3.2. Production run to predict the stress–strain curves of the nanocomposites After sufficiently minimized initial structures are prepared, an equilibrium process at a specified temperature and pressure is performed via statistical ensemble simulation methods. As the glass transition temperature of the PP is below room temperature (250 K), the target temperature is set to 200 K, where the PP molecules can be regarded as a glassy and rigid polymer. Using the canonical (NVT) ensemble simulation where the number of atoms N, the volume V, and the temperature T of the system are conserved, all the unit cells are equilibrated for 50 ps (pico second) with a time step of 1 fs (femto second) for time integration of the individual degrees of freedom (position and momentum) using the Nosé–Hoover extended Hamiltonian method (Hoover, 1985). Then, the unit cells are again equilibrated at 200 K and 1 atm for 3 ns (nano second) via the grand canonical (NPT) ensemble simulation where the pressure P of the systems is conserved using the Nosé–Hoover

Table 1 Composition of CNT-reinforced nanocomposites and pure PP matrix considered in MD simulations after equilibrated at 200 K and 1 atm. Volume fraction (%)

17.0 10.4 6.90 0

CNT diamter (Å)

11.58 11.58 11.58 –

CNT length (Å)

62.77 62.79 62.80 –

No. of PP chains

31 56 88 30

Cell density

0.95 0.91 0.88 0.83

Cell length X (Å)

Y (Å)

Z (Å)

62.77 62.79 62.80 35.81

31.24 41.60 49.62 37.20

32.82 40.72 51.40 38.04

S. Yang et al. / International Journal of Plasticity 41 (2013) 124–146

129

barostat (Hoover, 1986). For sufficient convergence of the pressure during the simulation, the pressure tensor of the unit cells during the NPT ensemble simulations are anisotropically controlled in each direction, thus, the Y and Z dimensions of the equilibrated unit cells listed in Table 1 differ each other. Due to the computational limitations in considering huge number of atoms in unit cell simulations of polymeric materials, such anisotropicity is usually inevitable in MD simulations. However, computational accuracy of the MD simulation results can be resolved from repeative calculations and averaging schemes from time to time. After the equilibration is finished, production runs are performed to obtain stress–strain curves of the unit cells by applying a continuous external strain to the unit cells. Due to the alignment of the CNT and the structure of the molecular unit cell in Fig. 1(b), it can be assumed that the elastic and secant stiffness of the nanocomposites in elastic and plastic deformation states are transversely isotropic. Thus, a total of six different loading conditions are applied: three uniaxial tension tests and three shear tests. From the three uniaxial tension tests, the stress–strain relation along the longitudinal direction and the two transverse directions can be obtained. As the XY and the XZ planes are the same symmetric planes, the two stress–strain curves of the transverse direction are averaged. Likewise, the shear stress–strain curves are obtained from the applied shear strain and the resultant shear stress, and the two in-plane shear stress–strain curves (in the XY and XZ planes) are averaged to resolve computational accuracy. In uniaxial tension simulation, according to the predefined strain rate of 0.0002/ps, a small amount of displacement is applied to the atoms at one edge of the unit cell and the system is relaxed via the NPT ensemble simulation. That is, a sudden change of the valence bond length, bending angle, dihedral angle and inter atomic distance, and the resultant increase of the potential energy, act as driving forces to increase the internal stress of the systems. After the relaxation process, the atoms at the edge are displaced again followed by the same relaxation process, and this process is repeated again until the strain of the unit cell reaches 20%. In order to mimic the real experimental test, the pressure tensor perpendicular to the loading direction is kept at 1 atm to reflect the Poisson effect, while the pressure along the loading direction is not controlled to maintain the strain state at each loading step. In each relaxation process, the time-averaged stress tensors calculated from the virial theorem are stored as the stress components corresponding to the strain state. The shear stress–strain curves are obtained in the same manner by sequentially changing the tilt angle of the unit cell at the same strain rate while the pressure tensor of the unit cell is not controlled, thus, the NVT ensemble simulation is used to relax the system after applying external shear strain. Again for computational accuracy, all the tensile and shear stress–strain curves are averaged over three different loading simulations by assigning different initial number seeds to determine the initial kinetic energy of the systems. As a result, uniaxial and shear stress–strain curves along the longitudinal direction are averaged over three different simulations, while the uniaxial stress–strain along the transverse direction and the in plan shear stress–strain curve are averaged over six different simulation results. For all the ensemble simulations including the equilibration process and production runs to drive stress–strain curves of nanocomposites, an open source code for large-scale atomic and molecular simulation tool LAMMPS (Plimpton, 1995) was used with the same force field (PCFF). 3.3. Prediction of the effective isotropic stress–strain relation From the MD simulation results, elastic and plastic behavior of a transversely isotropic unit cell can be obtained, thus, the transversely isotropic elastic stiffness and secant stiffness at each elastic and plastic strain state can be obtained. However, most of the manufactured nanocomposites have random orientation of the CNTs, thus, overall elastic and secant stiffness of the nanocomposites is almost isotropic. Moreover, due to the inherent computational inefficiency of the MD simulations, consideration of such random orientation by embedding many CNTs into the polymer matrix which leads to massive calculations to solve the equations of motion of individual atoms is almost impossible. As an alternative to address this computational limitation, we propose an orientation average of the transversely isotropic stiffness tensor at each strain state stored in MD simulations to generate an equivalent isotropic stiffness tensor, thus, to obtain the isotropic stress–strain relation. As the orientation average scheme is applied to every strain state, for convenience, we do not distinguish between the elastic and plastic regions in this chapter, and denote all the elastic constants as ‘secant’ components. From the stress–strain curves and the resultant cell dimensions obtained from the MD simulations, the longitudinal and transverse secant Young’s moduli, the corresponding secant Poisson’s ratios, and the secant shear moduli can be obtained at sec each strain state. Thus, the transversely isotropic secant compliance tensor Dsec comp and the stiffness tensor Ccomp of the nanocomposites can be constructed at each strain step as:

2

Dsec composites

1 Esec x

6 6 msec 6 secxy 6 Ex 6 6 msec 6 secxz E ¼6 6 x 6 0 6 6 6 0 6 4 0

msec yx Esec y

msec zx Esec z

1

msec zy

msec yz Esec y

0

0

0

0

1 Esec z

0

0

0

0

1 Gsec yz

0

0

0

0

1 Gsec xz

0

0

0

0

Esec y

Esec z

0

3

7 7 0 7 7 7 7 0 7 7; 7 0 7 7 7 0 7 7 5

1 Gsec xy

sec1 Csec composites ¼ Dcomposites ;

ð3Þ

130

S. Yang et al. / International Journal of Plasticity 41 (2013) 124–146

where the superscript sec indicates the secant components, and E, m, and G indicate the Young’s moduli, Poisson’s ratio, and the shear moduli, respectively. As the compliance and stiffness tensor obtained from MD simulations are not symmetric in nature since it is difficult to make fully homogeneous and symmetric polymeric structures, the off diagonal terms in Eq. (3) are again averaged to generate symmetric secant stiffness tensor components. Using Hill’s notation (Hill, 1964), the symmetric secant stiffness tensor in Eq. (3) can be simply written as: sec

sec

sec sec sec C sec comp ¼ ð2kcomp ; lcomp ; ncomp ; 2mcomp ; 2pcomp Þ;

ð4Þ

where k is the in-plane bulk modulus jyz, l is the cross modulus of C12, n is the axial modulus under an axial strain of C11, and m and p are respectively the in-plane and axial shear moduli lyz and lxy. Following Hill’s simple expression (Hill, 1965), the orientation-averaged isotropic stiffness tensor can be expressed using the two independent components of the bulk and shear moduli as: sec sec fC sec comp g ¼ ð3jcomp ; 2lcomp Þ;

ð5Þ

where the equivalent secant bulk modulus and the secant shear modulus are given as:

1 9 1 sec sec sec sec ¼  2lcomp þ nsec ½k comp þ 6ðmcomp þ pcomp Þ: 15 comp

sec sec ½4kcomp þ 4lcomp þ nsec jsec comp ¼ comp ;

lsec comp

ð6Þ

Once the equivalent isotropic secant stiffness tensors are obtained at each strain state, the equivalent isotropic stress– strain curve at uniaxial tension can be estimated from the secant Young’s moduli of the nanocomposites. 4. Nonlinear micromechanics constitutive model In the nonlinear micromechanics model to describe the elastoplastic behavior of the nanocomposites, the Mori–Tanaka mean field approach (Mori and Tanaka, 1973; Qiu and Weng, 1990) is used to calculate the effective elastic and secant constants of the nanocomposites. Under the assumption that the equivalent continuum structure of the present zigzag CNT is a solid cylinder, the elastic constants of the CNT are transversely isotropic (Yang et al. 2012b). Thus, the constitutive relation based on the Mori–Tanaka theory for two-phase composites involving transversely isotropic inclusion in isotropic media is adopted (Qiu and Weng, 1990), and its modification to account for the weakened bonding at the interface through a linear spring model is addressed (Qu, 1993; Barai and Weng, 2011a). In recent review of the Qu’s linear spring model, which reformulated the modified Eshelby tensor (Eshelby, 1957) component of a spherical inclusion, the singularity of the modified Eshelby tensor and the bounds of the admissible compliance of the linear spring have been demonstrated (Yanase and Ju, 2011). More recently, the same singularity problem of infinite cylinder type inclusion has been pointed out as well (Yang et al., 2012b). In order to make the modified Mori–Tanaka model for the weakened interface problem free from the singularity condition, the method of effective inclusion that replaces the weakened interface with a perfect one while the stiffness of the original inclusion is compatibly reduced (Yanase and Ju, 2010) is adopted in this study. For completeness of the present study, some important expressions in the constitutive model referred from other literature are selectively presented according to importance, while some omitted terms are described in detail for comprehension. 4.1. Method of effective fiber to account for the weakened interface In the linear spring model that describes the weakened interface via a finite value of spring compliance, it is assumed that the traction at the interface is continuous, whereas the displacement is discontinuous, given as:

Drij nj ¼ ½rij ðSþ Þ  rij ðS Þnj ¼ 0;

ð7Þ

Dui ¼ ui ðSþ Þ  ui ðS Þ ¼ gij rjk nk ;

ð8Þ

where S is the interface between the inclusion and matrix, with the superscripts + and – indicating positive (matrix) and negative (inclusion) sides, respectively, and nj is the outward unit normal at the interface to the positive side of S+. The second-order tensor gij is the linear compliance of the weakened interface and is defined as:

gij ¼ adij þ ðb  aÞni nj ;

ð9Þ

where a and b are the tangential and normal components of the compliance. On account of the displacement jump condition, there arises an additional term that contributes to the perturbed strain field inside the composites, and this perturbed strain field can be easily correlated to the eigenstrain field via the modified Eshelby tensor. The reformulated modified Eshelby tensor is written in an iterative form given below: M SM ðnþ1Þ ¼ S  S : R : Cm : ðI  SðnÞ Þ

or derived as a direct summation form as:

ð10Þ

S. Yang et al. / International Journal of Plasticity 41 (2013) 124–146

SM ¼ ðI  S : R : Cm Þ1 : S : ðI  R : Cm Þ;

131

ð11Þ

where the fourth order tensors S, R, and Cm are the original Eshelby tensor, interfacial compliance tensor, and the stiffness tensor of the matrix phase, respectively. The closed form solution of the original Eshelby tensor is available in the literature (Qiu and Weng, 1990), and the explicitly evaluated components of the interfacial compliance tensor (Barai and Weng, 2011a) are attached in Appendix A with one misprint fixed for the readers’ convenience. As the modified Eshelby tensor in Eqs. (10) and (11) has narrow admissible range of the interfacial compliance that describes the degree of separation and sliding at the interface, an alternative to resolve the singularity problem is required. The present study addresses the method of effective fiber (Yanase and Ju, 2010) that substitutes the modified Eshelby tensor with the original Eshelby tensor to prevent the singularity of the modified Eshelby tensor. Due to the displacement jump condition, the volume-averaged strain field of two-phase composites that are composed of filler (inclusion) and matrix is given as (Qu, 1993):

e ¼ fp ep þ fm em þ fp R : Cp : ep :

ð12Þ

By applying Hooke’s law to Eq. (12), we obtain:

e ¼ fm em þ fp ðI þ R : Cp Þ : C1  p : rp :

ð13Þ

Recalling the form of rule of mixtures, it can be intuitively found that the explicit expression of the stiffness tensor of the effective fiber can be obtained from the second term at the right hand side of Eq. (13) as:

~ p ¼ Cp : ðI þ R : Cp Þ1 : C

ð14Þ

When the interface is perfectly bonded, the compliance tensor becomes zero, thus, the stiffness of the effective fiber becomes that of the original fiber phase. 4.2. Modified Mori–Tanaka model with the weakened interface In the modified Mori–Tanaka model of two-phase composites with the effective fiber method, the closed form solution of the overall elastic stiffness is given in linear algebraic form:

C ¼ ðfm Cm þ fp Cp : Ap Þ : ðfm I þ fp Ap Þ1 ;

ð15Þ

where the fourth order tensor Ap is the dilute strain concentration tensor of the filler and is written as: 1 Ap ¼ ½I þ S : C1 m : ðCp  Cm Þ ;

ð16Þ

where I is the identity tensor. For randomly oriented fibers, the orientation-averaged isotropic stiffness tensor is given as:

fCg ¼ ðfm Cm þ fp fCp : Ap gÞ : ðfm I þ fp fAp gÞ1 ;

ð17Þ

where the bracket {} indicates the orientation average. In case of the weakly bonded interface, the modified Mori–Tanaka model involving randomly oriented fibers is written as:

fCg ¼ ðfm Cm þ fp fCp : Ap gÞ : ðfm I þ fp fAp g þ fp fR : Cp : Ap gÞ1 :

ð18Þ

In the present study, as the weakened interface is replaced with the perfect interface via the effective fiber method, the Eshelby tensor in the dilute strain concentration tensor of Eq. (16) is the original one, while the stiffness of the fiber in Eqs. (16)–(18) is the stiffness of the effective fiber given in Eq. (14) hereafter. As the stiffness tensor given in Eq. (18) is isotropic, it can be simply expressed in terms of the two independent elastic constants of the bulk and shear moduli following Hill’s notation as:

fCg ¼ ð3jc ; 2lc Þ;

ð19Þ

where the two elastic constants are given as:

jc ¼

fm jm þ fp nCA A

RCA

fm þ 3f p n þ 3f p n

;

lc ¼

fm lm þ fp gCA : fm þ 2f p gA þ 2f p gRCA

ð20Þ

Following the tensor manipulation process of the Mori–Tanaka model involving transversely isotropic fiber (Qiu and Weng, 1990) via Walpole’s mathematical scheme (Walpole, 1981), the complete expression of the components nA, nCA, nRCA and gA, gCA, gRCA can be obtained without difficulty. The complete expressions of these terms are given in Appendix B. In the plastic region, all the elastic constants of the matrix phase in Eq. (20) are replaced with the corresponding secant constants to draw stress–strain curves.

132

S. Yang et al. / International Journal of Plasticity 41 (2013) 124–146

5. Hierarchical multiscale modeling approach 5.1. Constitutive relation of PP matrix In order to utilize micromechanics modeling to predict the behavior of CNT reinforced nanocomposites, a suitable plastic constitutive relation should be established using the stress–strain curve of PP matrix obtained from MD simulations. Different from the experimental result of typical polymers that have smooth variation, MD simulation results generally show fluctuations in the stress–strain curve. In order to obtain best fit of the MD simulation results, we performed a nonlinear least square method. To express the yield surface of PP matrix, the Ludwik’s model is adopted as a yield function given as:

re ¼ rY þ hðep Þn ;

ð21Þ p

where rY, h and n are the yield strength, strength coefficient, and hardening exponent, and re and e are the von Mises’ effective stress and the effective plastic strain of the matrix, respectively, defined as:

re ¼

 1=2 3 0 0 rij rij ; 2

ep ¼

 1=2 2 p p eij eij ; 3

ð22Þ

where r0ij are the deviatoric components of the stress tensor. In order to use the uniaxial tensile stress–strain curve of pure polypropylene obtained from MD simulation in the micromechanics model, the coefficient of the Ludwik’s model that can best fit the fluctuating stress–strain curve obtained from MD simulation to a smooth and continuous curve should be decided. In this study, a non-linear least square method is introduced to determine the unknown coefficients of the Ludwik’s model and the details are as follows. First of all, it should be noted that in addition to the three unknown coefficients (ry, h, n) of Ludwick’s model, the Young’s modulus of PP matrix is also selected as an unknown coefficient, because the Young’s modulus obtained from MD simulation is sensitive to the number of data points used to fit the linear elastic stress–strain curve from MD simulation results. In elastoplastic behavior, the constitutive equation is expressed in terms of the Young’s modulus, total strain, and plastic strain as:

re ¼ Eðe  ep Þ;

ð23Þ

where E is the Young’s modulus and e is the total strain of pure PP matrix. The four unknown coefficients to reproduce the stress–strain relation of PP matrix can be numerically obtained by substituting Eq. (23) into Eq. (21) given as:

rY þ hðep Þn ¼ Eðe  ep Þ

ð24Þ

and the computational iteration procedures are as follows. At each strain state stored in the MD simulation, initial values of the four unknowns are assigned, and the plastic strain of PP matrix is obtained from Eq. (24). From the plastic strain, the effective stress of PP matrix is calculated from Eq. (21) or from Eq. (23) and compared with the corresponding stress values in MD simulation results. At this evaluation stage, we define the residue of the effective stress as a square sum over all strain state as:

R¼

N X

½ri  frY þ hðepi Þn g2 ;

ð25Þ

i¼1

where the subscript i indicates the ith strain, state and N is the total number of data points stored in MD simulations. Monitoring the magnitude of the residue, new rounds of iteration are continuously performed until the residue satisfies the predefined convergence criterion. In detail, the Newton–Raphson algorithm is used in the present non-linear least square fit process. The first and second derivatives (Hessian) of each unknown coefficient are used. Once the constitutive relation of the PP matrix is given from the MD simulations, all the secant moduli of the PP matrix at each plastic strain state can be calculated. In nonlinear micromechanics model using the Mori–Tanaka’s scheme, the secant moduli method (Talbot and Willis, 1985; Tandon and Weng, 1988; Qiu and Weng, 1992; Suquet, 1995) can be used to estimate the nonlinear behavior of nanocomposites. In the secant moduli method, the plastic strain and the corresponding stress above the yield strength of composites is described by defining a linear comparison material. Once the internal stress of the matrix phase in composites reaches its yield strength, the elastic moduli of the matrix phase is replaced with the secant moduli, and the corresponding secant moduli, stress, and strain of the comparison material is calculated from the secant moduli of the matrix and the elastic moduli of the filler phase step by step to update the plastic behavior of the composites. Besides the secant moduli method adopted in this study, other nonlinear models (Hill, 1965; Hutchinson, 1976; Berveiller and Zaoui, 1979; Fotiu and Nemat-Nasser, 1996; Pierard and Doghri, 2006; Mercier and Molinari, 2009) can be selectively applied to describe the nonlinear behavior of nanocomposites according to the advantageous features of each models, and to the specialty of the extraordinary behavior uniquely revealed by nanoscale microstructures. From the constitutive relation of the PP matrix, the secant Young’s and shear modulus of the matrix phase at a given plastic strain can be obtained as:

S. Yang et al. / International Journal of Plasticity 41 (2013) 124–146

Esm ¼

1 1 Em

þr

;

ep

Y þhðe

lsm ¼

p Þn

Esm ; 2ð1 þ msm Þ

133

ð26Þ

where the super script s indicates the secant modulus, and m is Poisson’s ratio. For von Mises solids, the bulk modulus is kept constant throughout the elastic and plastic deformation, thus, the following relations are satisfied: s

km ¼

Esm ¼ km ; 3ð1  2msm Þ

1 2

msm ¼ 

  s 1 E  mm m : 2 Em

ð27Þ

5.2. Calculation of stress in individual phase When the secant moduli method is used to describe the plastic behavior of composites, the most important part is to calculate the stress of matrix phase with accuracy because yielding of the composites is evaluated from the local yielding of matrix phase that surrounds the rigid inclusion phase. In the present study, the field fluctuation method based on the average second order stress moment (Bobeth and Dienner, 1986; Kreher and Pompe, 1989; Suquet, 1995; Hu, 1996, 1997) is adopted. In the field fluctuation method, the fluctuation of the local compliance tensor of an individual phase gives rise to a change in overall stored energy, thus, the perturbation of the macroscopic compliance is correlated to the local field fluctuation given as:

r : dD : r ¼

X hrr : ddr : rr iV r ;

ð28Þ

r

 and where D and dr are the overall compliance tensor of composites and the local compliance tensor of the rth phase, and r rr are the uniform macroscopic stress and the local stress tensor of the rth phase, respectively. The brackets hiVr indicate the volume-averaged quantity. When the shear modulus of the rth phase only undergoes a small variation, by splitting the stress tensor into deviatoric and spherical parts, Eq. (28) can be reduced to:

      Vr 0 1 1 1  m2 d 0 : r 0d ¼r þr ; hrr : r0r iV r d 2l r 2lc V jc

ð29Þ

 m is the average stress of composites. The average stress of the where the superscript 0 indicates the deviatoric part, and r matrix phase obtained from the second order momentum is more accurate than the method based on the first order momentum, i.e. via the global stress concentration tensor, and can prevent underestimation of the overall stress–strain relation of composites (Qiu and Weng, 1992). Using the field fluctuation method, the average effective stress of matrix phase is defined as:

(

r em ¼

1 3 fm



lm jc

2

@ jc m2 r þ @ lm



lm lc

2

@ lc e2 r @ lm

!)1=2 ;

ð30Þ

 e and r  em are the effective stress of composites and matrix phase, respectively. In the plastic region, all the elastic where r constants in Eq. (30) are replaced with the corresponding secant constants. In order to utilize Eq. (30) to calculate the effective stress of the matrix phase in plastic deformation state, the secant shear and bulk modulus of the matrix phase are required. These properties are also functions of the effective stress of matrix phase as defined in Eqs. (26) and (27). Thus, in general, a numerical iteration scheme has been applied to solve Eq. (28) (Barai and Weng, 2011a) in the plastic region. Meanwhile, a more recent application of the field fluctuation method has suggested a direct calculation method (Zhang et al., 2010). Different from the iterative method, where the stress (or strain) state of composites is firstly defined in Eq. (30), one can prescribe the effective strain and resultant secant moduli of matrix phase and can solve for the macroscopic stress of composites in a direct calculation method. In this study, both methods have been selectively chosen according to the situation, and the details of the computational procedures applied at each analysis step will be demonstrated later. 5.3. Upper level analysis to characterize effective matrix There are two main objectives in upper level analysis: one is to find out the lower bound of the interfacial compliance between the embedded CNT and surrounding phase, and the other is to obtain an effective stress–strain relation of the effective matrix which contains both the densified PP molecules that constitute the interphase and the rest of the PP matrix. In order to extract important physical information on the interfacial strength of the CNT nanocomposites, it is essential to compare the MD simulation results with the conventional micromechanics model prediction that assumes a perfect bonding condition between the fiber and matrix phase. At the same volume fraction conditions considered in MD simulations, the elastoplastic stress–strain curves of nanocomposites with a random orientation of CNTs are obtained from the Mori–Tanaka model with the secant moduli method by setting the interfacial compliance as zero, and comparing with the MD simulations to identify the difference of the stress–strain curves obtained from the two different methods. Then, by changing the interfacial compliance from zero to a finite value, the lower bound of the interfacial compliance that can best reproduce the MD

134

S. Yang et al. / International Journal of Plasticity 41 (2013) 124–146

simulation results from the Mori–Tanaka model can be easily estimated. As the elastoplastic behavior of pure PP matrix is used in this micromechanical model prediction without considering the polymer densification effect for the evaluation of the interfacial strength, the corresponding compliance found from this stage is the smallest value to reproduce the MD simulation results from the micromechanics model. In the real situation, however, the interphase region plays an important role to enhance the mechanical behavior of the surrounding matrix. If the elastoplastic behavior of the matrix considered in the micromechanics model is stiffer than the pure PP matrix, the resultant interfacial compliance that can best fit the same MD simulation in the micromechanics model increases correspondingly. Thus, the interfacial compliance that can reproduce the MD simulation from the micromechanics model prediction is referred as the lower bound. Even though the exact value of the interfacial compliance is not known at this stage, knowing the possible lower bound of the compliance is truly valuable for further processes. Together with the interfacial compliance that degrades the overall elastoplastic behavior of nanocomposites, the elastoplastic behavior of the interphase should be considered as well. As has been mentioned, however, it is quite difficult to exactly estimate the elastoplastic behavior of the interphase and the corresponding interfacial compliance from a given stress– strain relation of nanocomposites because they rival each other. Alternatively, the combination of the interfacial compliance and the relevant elastoplastic behavior of the interphase can be obtained by combining MD simulation results and micromechanics model prediction. In order to identify the elastoplastic behavior of interphase, the nanocomposites microstructure is firstly divided into CNTs and effective matrix phase in upper level analysis. Different from conventional micromechanics model that predicts the overall behavior of composites from known behavior of an individual phase, the elastoplastic behavior of the effective matrix that includes both the interphase and pure matrix phase is yet unknown, thus, a different route is required. However, as the overall elastoplastic behavior of nanocomposites is given by the MD simulation, the overall elastoplastic behavior of effective matrix can be inversely estimated from the micromechanics constitutive equations. At each strain step of nanocomposites, the elastic constants of effective matrix are numerically obtained from Mori–Tanaka solution given as:

jc ¼

fem jem þ fp nCA A

fem þ 3f p n þ 3f p n

RCA

;

lc ¼

fem lem þ fp gCA : fem þ 2f p gA þ 2f p gRCA

ð31Þ

In the present inverse estimation method, the terms at the left hand side of the above equations are given values. While the terms at the right had side with the subscript ‘em’ are unknowns to be numerically determined. In this inverse process, as the numerical calculation process to obtain the properties of the interphase proceeds without defining the effective plastic strain of nanocomposites, it is not required to distinguish elastic and plastic regions in this part, thus, we replace the elastic constant with the secant constant for simplicity. Once the secant constants of the effective matrix at a certain strain state of the nanocomposites are decided, the corresponding effective stress of the effective matrix is obtained from the field fluctuation method as:

(

r eem ¼



1

3

fem

lem jc

2

@ jc m2 r þ @ lem



lem lc

2

@ lc e2 r @ lem

!)1=2

:

ð32Þ

In order to obtain several possible sets of interfacial compliance and corresponding elastoplastic behavior of interphase, the interfacial compliance was gradually increased from the lower bound. In this study, five different interfacial compliance values including the lower bound and corresponding elastoplastic behavior of interphase are obtained. The elastoplastic behavior of the effective matrix obtained in upper level analysis is again used as a given elastoplastic behavior to finally obtain the elastoplastic behavior of interphase in lower level analysis. The numerical iteration scheme and the convergence criterion of the upper level analysis will be demonstrated later. 5.4. Lower level analysis to characterize effective interphase Because the effective matrix consists of the interphase and matrix, it can be again regarded as a heterogeneous structure where the interphase is embedded into a pure PP matrix as a reinforcement. Thus, the inverse micromechanics method used in Section 5.3 can be applied to identify elastoplastic behavior of the interphase again. While the behavior of base material was an unknown value to be determined in upper level analysis, the behavior of reinforcement now becomes an unknown value in lower level analysis. In lower level analysis, however, the interface between interphase and matrix is perfectly bonded because there is no obvious structural discontinuity. Similar to the process used in upper level analysis, the elastic constants of the interphase are numerically obtained from the following Mori–Tanaka solution using the known behavior of pure PP matrix and the effective matrix obtained in part 5.5:

jem ¼

fmem jm þ fiem nCA fmem

þ

3fiem nA

;

lem ¼

fmem lm þ fiem gCA : fmem þ 2fiem gA

ð33Þ

In Eq. (33), the orientation averaged terms nA, nCA and gA, gCA are obtained from the stiffness tensor and the resultant dilute strain concentration tensor of the interphase which is regarded as a reinforcement embedded into the pure PP matrix. Once the elastic constants of the interphase at each strain step of the effective matrix are obtained, the corresponding stress of the interphase is again obtained from the field fluctuation method as:

S. Yang et al. / International Journal of Plasticity 41 (2013) 124–146

(  ei

r ¼

1

fiem



li 3 jem

2

@ jem m2 r þ @ li em



li lem

2

@ lem e2 r @ li em

135

!)1=2 ;

ð34Þ

where the subscript ‘i’ indicates the interphase. After calculating the stress, the corresponding strain of the interphase is obtained from Hooke’s law, and the final stress–strain of the interphase at the given interfacial compliance considered in upper level analysis is deducted. Different from the upper level analysis where two elastic constants of the unknown phase are assumed in iteration process, in our experience, the numerical iteration was not stable in lower level analysis by setting two unknown variables. At the same time, as the pure matrix phase was assumed as von Mises solid, it is reasonable that the interphase also follows the von Mises yield criterion. That is, the bulk moduli of the interphase should be constant throughout the elastic and plastic deformation. In order to prevent the numerical instability and to satisfy the von Mises criterion, we set only one elastic constant as an unknown value while the bulk modulus of the interphase is predefined. For this, the Poisson’s ratio of the interphase was set as 0.32 which is identical to the pure PP matrix. The way to predefine the bulk modulus of the interphase in this assumption is as follows. In the small strain range, the shear modulus of the interphase is set as an unknown while a finite value of the bulk modulus is assigned in iteration steps. By gradually increasing the predefined bulk modulus of the interphase at each iteration process while the shear modulus is still set as an unknown value, and by monitoring the resultant Poisson’s ratio of the interphase that satisfies the convergence criterion, we can easily find the bulk modulus where the Poisson’s ratio of the interphase becomes 0.32. Then, using the bulk modulus of the interphase found from this small strain range as a predefined value while the shear modulus of the interphase is still set as an unknown variable, the same iteration process is performed throughout the whole strain range of the effective matrix to obtain other elastic constants of the interphase. The convergence criterion applied in lower level analysis is the same as in upper level analysis, and the details will be demonstrated in next part. Once the elastoplastic behavior of the interphase is obtained with the corresponding interfacial compliance, the overall elastoplastic behavior of nanocomposites at various volume fractions can be efficiently predicted from the micromechanics model itself. In this process, the overall computational process is identical to those of the conventional micromechanical approach. However, a two-step homogenization is required. Using the elastoplastic behavior of the interphase obtained from lower level analysis, the elastoplastic behavior of the effective matrix at any volume fraction of CNT is firstly obtained. Then, using the elastic properties of CNT and the elastoplastic behavior of the effective matrix, the overall elastoplastic behavior of nanocomposites is finally predicted. 6. Computational procedures 6.1. Estimation of the overall elastoplastic behavior of composites: iterative method The conventional computation process to obtain the plastic behavior of composites from secant moduli method combined with field fluctuation method is well known, and various schemes have been suggested. In our case, macroscopic strain of the composites is firstly defined. After the matrix phase in nanocomposites reaches the yield condition at a predefined strain state of the nanocomposites, a small amount of its plastic strain is assumed. Then, the secant moduli of the matrix phase are obtained from Eqs. (26) and (27) and substituted into Eq. (20) to calculate the bulk and shear modulus of the nanocomposites. After substituting the moduli of the nanocomposites into Eq. (30), the effective stress of the matrix phase is obtained and compared with the effective stress obtained from the constitutive relation in Eq. (21). If the difference satisfies the convergence criterion, the secant moduli of nanocomposites at such strain state is decided. In the iteration process, we used the following convergence criterion:

  em field error ¼ r

fluctuation

  em constitutiv e r

relation

:

ð35Þ

When the error becomes smaller than 106, we decided that the results are converged. This iteration scheme was used in Section 5.3 to find the lower bound of the interfacial compliance. 6.2. Estimation of the overall elastoplastic behavior of composites: direct method Different from the iterative method that prescribes the macroscopic stress (or strain) of nanocomposites to obtain the effective stress of matrix phase, the direct method prescribes the effective stress of the matrix phase which is described via the constitutive relation in Eq. (21). At the given effective stress and the resultant secant moduli of the matrix obtained from Eqs. (26) and (27), the overall secant moduli of nanocomposites can be simply calculated, thus, the macroscopic stress of nanocomposites can be simply solved for at the predefined stress of the matrix phase from the field fluctuation method. The direct method was used to predict the overall stress–strain relation of nanocomposites at various volume fractions using the stress–strain relation of effective matrix and the elastic constants of CNT. Once the stress–strain relation of the interphase is obtained from the proposed multiscale model, the stress–strain relation of effective matrix at any volume fraction of CNT can be simply predicted from lower level analysis. Because the typical output format of this calculation is the set of strain, stress and secant moduli of the effective matrix, the direct method is more convenient in the hierarchical approach

136

S. Yang et al. / International Journal of Plasticity 41 (2013) 124–146

than the iterative method, which requires the nonlinear least square fitting process that was introduced in Section 5.1 to express the stress–strain curve of the effective matrix in terms of the effective plastic strain and hardening parameters given as in Eq. (20). Importantly, referring to the Eqs. (30), (32), and (34), the plastic behavior of the nanocomposites and the effective matrix are no longer independent of the hydrostatic pressure. It is essential in the iterative method that the plastic behavior of the matrix phase (soft phase) is assumed to be plastically incompressible to calculate the secant bulk moduli of the matrix phase from Eq. (27), while the direct method does not apply to this assumption.

6.3. Estimation of the overall elastoplastic behavior of effective matrix: combination of the iterative and the direct method The most complicated computational procedure of the proposed hierarchical multiscale modeling approach is to estimate the elastoplastic behavior of effective matrix from the constitutive relation of pure matrix phase and the stress, strain, and secant moduli of interphase obtained from the proposed inverse identification process. As in the case of the secant moduli method of polycrystalline structure (Barai and Weng, 2011b) where the constitutive relation of the matrix (grain boundary zone) and the filler (grain zone) are elastoplastic, the constitutive relation of the pure PP matrix as well as the interphase are elastoplastic. In this case, the iterative method to calculate the plastic behavior of effective matrix defines two unknown variables, and this situation can make the iteration process unstable. Alternatively, we combined the iterative method and the direct method to obtain the stress–strain curve and the corresponding secant moduli of the effective matrix. As the elastoplastic behavior of the interphase is obtained as a set of strain, stress, and elastic (secant) moduli from lower level analysis, the computation process begins by defining the stress of the interphase in the effective matrix. Using the given stress and moduli of the interphase and the elastic constants of pure matrix, the elastic constants of the effective matrix can be obtained from Eq. (33). From Eq. (34), the macroscopic stress of the effective matrix can be directly solved for from the predefined effective stress of interphase. The effective stress of pure matrix can also be calculated from Eq. (34) by replacing all the subscripts ‘i’ with ‘m’. After the pure PP matrix phase reaches the yield condition, at predefined stress and secant moduli of the interphase and the resultant macroscopic stress of the effective matrix calculated from Eq. (34), a small amount of effective plastic strain of pure PP matrix is assumed, and the rest of the iteration scheme follows the same process introduced in Section 6.1. When the iteration converges, the stress and corresponding moduli of the effective matrix at the predefined effective stress of interphase is completely decided. Then, the next effective stress of the interphase is defined in turn and the same iteration is repeated until the stress of the effective matrix corresponds to the final set of the effective stress, strain, secant moduli of interphase is completely obtained.

6.4. Inverse estimation of the elastic constants of effective matrix and interphase In the inverse estimation process for the elastoplastic behavior of the effective matrix introduced in Section 5.3, the elastic constants and stress–strain curves of the nanocomposites are given by the MD simulation results. As the effective isotropic stress–strain curve and secant moduli of nanocomposites obtained from MD simulation are orientation-averaged quantities from individual stress–strain curves that show fluctuations, the resultant stress–strain curve of the effective matrix inversely obtained from the MD simulation results also have the same fluctuation in its evolution. Alternatively, one can use the reproduced data from micromechanics model with the lower bound of the interfacial compliance to obtain the same stress–strain curve while the results are smooth without fluctuations. We simply used the data reproduced from the micromechanics model in this process. At given strain, stress, and moduli of the nanocomposites, one can choose any two elastic constants of the effective matrix among the Young’s modulus, shear modulus, bulk modulus, and Poisson’s ratio as unknown variables to be defined in the iteration process. We set the bulk modulus and the shear modulus of the effective matrix as unknown variables, and an initial guess starts from the moduli of the pure PP matrix. Using the assumed moduli of the effective matrix and elastic constants of the CNT, the bulk and shear moduli of the nanocomposites are obtained from Eq. (31) and compared with the given bulk and shear moduli of the nanocomposites. For evaluation of the iteration, we used the following convergence criterion using the norm of the difference between predicted and given moduli of the nanocomposites:

error ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðjcgiv en  jpredict Þ2 þ ðlcgiv en  lpredict Þ2 : c c

ð36Þ

The same tolerance of 106 used in Section 6.1 is applied for convergence evaluation. The numerical procedures for inverse estimation of the interphase from the given elastoplastic behavior of the effective matrix are almost identical, while only one elastic constant is defined as an unknown variable (we set the shear moduli of the interphase) in the iteration step as has been briefly elucidated in Section 5.4. Even though only one elastic constant is defined, more than one elastic constant can be chosen to evaluate the convergence of the iteration. In our experience, the following criterion generated the best solution:

error ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðEcgiv en  Epredict Þ2 þ ðlcgiv en  lpredict Þ2 : c c

The same tolerance was used.

ð37Þ

137

S. Yang et al. / International Journal of Plasticity 41 (2013) 124–146

7. Results and discussions 7.1. MD simulation results of pure PP matrix and nanocomposites The stress–strain curve of pure PP matrix obtained from MD simulation is shown in Fig. 3(a). At the strain range less than 5%, the stress of PP gradually increases according to the imposed external strain. After 5% of strain, the increase of the stress is saturated and the slope of the stress–strain curve reaches a plateau showing a gentle hardening behavior. In addition, there is no obvious overshoot or softening behavior which are commonly observed in the plastic behavior of typical polymeric materials. Despite the fluctuations in stress evolution of the PP matrix, the plastic behavior of polypropylene can be effectively described via the Ludwik’s model without controversy. The nonlinear least square fitted stress–strain curve of PP matrix is compared with the result from MD simulation in Fig. 3(b). The resultant variables defined in Ludwik’s model are ry = 41 MPa, h = 15.91 MPa, n = 0.092 respectively. The elastic moduli and the parameters obtained to describe the plastic deformation of pure matrix are then used in the micromechanics model. The elastic constants of (15, 0) CNT for micromechanical approaches are referred from molecular mechanics simulation results (Yang et al., 2012b) and arranged in Table 2. In order to utilize the linear elastic constitutive relation to calculate the elastic constants of CNT, it has been assumed that the equivalent continuum model of CNT atomic structure is a solid cylinder with the effective wall thickness of 3.4 Å. The uniaxial and shear stress–strain curves of nanocomposites obtained from MD simulations are depicted in Fig. 4 and compared with those of pure PP polymer. In uniaxial tension, due to the exceptionally high elastic modulus of CNTs, the stress increases above 10GPa and no clear yielding behavior is observed. Moreover, the evolution of the stress–strain relation gradually becomes concave, which is not in accordance with the typical stress–strain relation of polymeric materials. This is due to the limitation of the inter-atomic potential used in the present MD simulation that the covalent bond breaking or bond order change at high strain of CNT cannot be described. Thus, the resultant plastic deformation of nanocomposites in uniaxial loading does not occur within the description of the present inter-atomic potential. According to the MD simulation result of CNTs with a reactive bond order potential such as REBO potential (Liew et al., 2004), CNTs gradually show plastic deformation above 10% of strain and failure at 20% of strain, showing a clear bond order transition. Thus, the MD simulation results of less than 12% of strain are considered further in the process of the multiscale modeling. On the other hand, the stress–strain relation in transverse loading shows similar plastic behavior to the pure PP polymer as shown in Fig. 4(b). One important feature obtained from the present MD simulation is that the stress–strain relations of 6.8% and 10% nanocomposites are almost identical to that of the PP polymer, indicating that the interface between CNTs and surrounding PP molecules is imperfect. However, the stress–strain relation of 17% reinforced nanocomposites clearly shows enhanced reinforcing effect. The reason for this discrepancy does not mean that the interfacial strength is enhanced by the volume fraction, because the same CNT and matrix material is used to construct the unit cell. From our own analysis on the in-plain Poisson’s ratio of nanocomposites in transverse loading, it has been found that the in-plain Poission’s ratio myz of the nanocomposites increases as the volume fraction of CNT increases. As a result, even at the same external strain state, CNTs in smaller unit cells (having larger volume fraction) are more likely to be compressed perpendicular to the loading direction, thus, the internal stress of the CNTs along the loading direction increases. Similar to the transverse loading, the stress–strain curves of nanocomposites in shear loadings also show the same weakened interfacial bonding as depicted in Fig. 4(c) and (d). As the longitudinal and in-plane shear moduli of CNT is much greater than the shear modulus of pure PP polymer, and as the non-bond interaction between the CNT and PP matrix is weak, the

70

(b) 70

60

60

50

50

Stress(MPa)

Stress(MPa)

(a)

40 30 20

30 20 MD simulation result Least square fitted

10

10 0

40

0

0.05

0.1

Strain

0.15

0

0

0.05

0.1

0.15

Strain

Fig. 3. Stress–strain curve of pure PP matrix obtained from MD simulations at 200 K and 1 atm. (a) MD simulation result of PP matrix and (b) Least square fitted result of PP matrix. Elastic constants of PP matrix obtained from 3% of strain range are E = 1.379 GPa, G = 0.52 GPa.

138

S. Yang et al. / International Journal of Plasticity 41 (2013) 124–146

Table 2 Transversely isotropic elastic constants of CNT/PP nanocomposites obtained from molecular dynamics simulations (GPa). Chirality

EL

ET

GL

GT

KT

(15, 0)

775

64

378

17

281

120

30 Vol.f: 6.8% Vol.f: 10% Vol.f: 17%

100

20

Stress(MPa)

Stress(GPa)

25

15 10 5 0

80 60 40 20

0

0.05

0.1

0

0.15

0.05

0.1

Strain

(a) Longitudinal tension

(b) Transverse tension

0.15

80 Vol.f: 6.8% Vol.f: 10% Vol.f: 17% Pure PP

Vol.f: 6.8% Vol.f: 10% Vol.f: 17% Pure PP

70 60

Stress(M Pa)

50

0

Strain

60

Stress(MPa)

Vol.f: 6.8% Vol.f: 10% Vol.f: 17% Pure PP

40 30 20

50 40 30 20

10

10 0 0

0.05

0.1

0.15

0

0

0.05

Strain

0.1

0.15

Strain

(c) Longitudinal shear

(d) In-plane shear

Fig. 4. Stress–strain curves of CNT/PP nanocomposites obtained from MD simulations. Due to the critical anisotropicity of the CNT, longitudinal tensile stress reaches 10 GPa.

Radial Density Distribution

2 PP in nanocomposites Pure PP 1.5

1

0.5

0 5

10

15

20

Radial Distance ( ) Fig. 5. Radial density distribution of surrounding PP matrix in nanocomposites.

139

S. Yang et al. / International Journal of Plasticity 41 (2013) 124–146

CNTs are not distorted in shear loading of nanocomposites while it slips at the interface maintaining its original shape. From the stress–strain curves of nanocomposites, it can be reasonably concluded that the interface between a pristine CNT and pure PP matrix is imperfect, thus, it should be rigorously considered in equivalent continuum modeling of the nanocomposites via the micromechanics model. The uniqueness of the interphase surrounding the CNT is confirmed from the radial density distribution of the PP matrix in nanocomposites. Fig. 5 shows the radial density distribution of PP matrix with the reference density of homogeneous pure PP structure. As the radius of (15, 0) including half of the wall thickness of 3.4 Å is 7.49 Å, the non-zero density profile of PP matrix starts above this radial distance and abruptly increases near the surface of CNT, showing a highly densified region. Right after the density decreases to the level of homogeneous PP matrix, it fluctuates near the horizon of homogeneous density, indicating that the rest of the matrix region is regarded as pure PP matrix phase. As the density profile in the vicinity of the CNT is conical, from the radial distance where the peak of the density appears, we defined the outer thickness of the interphase as 4 Å. Thus, the thickness of the interphase is about 6.8 Å, and we simply set it as 7 Å for further process in the micromechanics constitutive model. Even though we only demonstrate the radial density distribution of surrounding polymer, other related simulations on nanoparticulate composites have revealed that the radial stress distribution of surrounding polymer shows similar transition to the radial density distribution of surrounding polymer (Papakonstantopoulos et al., 2005). 7.2. Orientation averaged stress–strain relation and lower bound of interfacial compliance From the transversely isotropic stress–strain relation obtained from MD simulations as demonstrated in the previous part, the equivalent isotropic stress–strain curves for randomly oriented nanocomposites can be generated from the orientation-averaged scheme introduced in Section 3.3. At the same time, the orientation-averaged stress–strain behavior of nanocomposites can be estimated from the Mori–Tanaka nonlinear micromechanics model using the properties of pure PP polymer and CNT with a perfect bonding condition to evaluate the interfacial strength. The orientation-averaged stress–strain curves of nanocomposites obtained from the two different methods are compared in Fig. 6. Because the longitudinal stress–strain relation dominates the overall stress–strain behavior of nanocomposites, the resultant isotropic stress– strain curves linearly increase according to the external strain. As the CNT waviness and the local CNT clustering effect could be naturally excluded in our molecular model, the MD simulation results are the possible upper bound that we can expect by embedding pristine CNTs into PP matrix. As has been demonstrated from the previous two-scale micromechanics model that considered the CNT clustering effect, the plastic deformation of nanocomposites is rather severely degraded by the CNT agglomeration (Barai and Weng, 2011a). Moreover, a small waviness of CNT has been found to greatly affect the overall elastic moduli of CNT-reinforced nanocomposites (Fisher et al., 2003). If both parameters were considered in MD simulation, the resultant stress–strain relation would be more realistic than the present MD simulation results. However, it is desirable to exclude such parameters in identifying the effects of interphase and interfacial strength. After identifying them, some parametric studies can be performed considering all possible factors together via a fully continuum-based analysis, and we leave it for future work. One important feature that we can find from the compared results in Fig. 6 is that the difference between the MD simulation results and Mori–Tanaka model prediction increases as the volume fraction of the CNT increases. Again from this result, we can definitely conclude that the interfacial strength is imperfect. However, at the volume fraction of 6.9%, the MD simulation results are almost identical to the Mori–Tanaka model prediction. This variation implies that one can lose important physical information that should be rigorously considered in equivalent continuum modeling of nanocomposites unless

Stress(MPa)

3000 2500 2000

Vol.f: Vol.f: Vol.f: Vol.f: Vol.f: Vol.f:

(b)

17.17% M-T 17.17% MD 10.04% M-T 10.04% MD 6.90% M-T 6.90% MD

1500 1000 500 0 0

4000 3500

Stress(M Pa)

(a) 3500

3000 2500

Vol.f: Vol.f: Vol.f: Vol.f: Vol.f: Vol.f:

17.17% M-T 17.17% MD 10.04% M-T 10.04% MD 6.90% M-T 6.90% MD

2000 1500 1000

0.02

0.04

0.06

0.08

0.1

0.12

500 0.06

0.07

0.08

0.09

0.1

0.11

0.12

Strain Fig. 6. CNT orientation-averaged stress–strain curves of nanocomposites. (a) Comparison of MD simulation results with the Mori–Tanaka model prediction with the perfect bonding condition, (b) magnified view to emphasize the difference between the two different methods.

140

S. Yang et al. / International Journal of Plasticity 41 (2013) 124–146

various volume fraction conditions are considered in MD simulations. Moreover, as the Mori–Tanaka model used to plot the stress–strain curves in Fig. 6 did not consider the effect of interphase, it can be reasonably conjectured that the difference in stress–strain curves can be addressed by considering the effect of interfacial imperfection and the interphase. Even if the exact values of the interfacial compliance and the corresponding properties of the interphase are yet unknown in detail, the possible lower bound of the interfacial compliance and the possible combinations of the finite interfacial compliance and interphase elastoplastic behavior can be systematically investigated via the proposed hierarchical multiscale scheme. First of all, by gradually increasing the interfacial compliance of a defined in Eq. (9), the lower bound of the interfacial compliance of the three nanocomposites could be easily obtained as 0.5, 1.45, and 1.2 nm/GPa for 6.9%, 10.4% and 17.17% of volume fractions respectively. Because all the nanocomposite molecular structures considered in MD simulations have been found to have similar lower bounds of the interfacial compliance near 1 nm/GPa, nanocomposites with 6.9% volume fraction are considered in the rest of the present multiscale approach to obtain the combinations of interfacial compliance and the elastoplastic behavior. The stress–strain curves of the effective matrix in 6.9% nanocomposites inversely obtained from the upper level analysis using the elastoplastic behavior of nanocomposites are depicted in Fig. 7 at various interfacial compliances. When the lower bound of the interfacial compliance is used, the resultant stress–strain curve of the effective matrix exactly coincides with that of the pure PP matrix. As the interfacial compliance increases, the corresponding elastic and secant moduli and the resultant stress of the effective matrix increase. At 1.0 nm/GPa of interfacial compliance, the resultant elastic modulus of the effective matrix is 1.75 GPa, while at 2.5 nm/GPa compliance it becomes 2.88 GPa, which is twice that of the pure PP matrix. The stress–strain curves of the interphase corresponding to the given elastoplastic behavior of the effective matrix are shown in Fig. 8 according to the interfacial compliance. Similar to the results in Fig. 7, the slope of the stress–strain curves increases as the corresponding interfacial compliance increases. Compared with the stress–strain behavior of pure PP matrix, the behaviors of interphase at the interfacial compliance above the lower bound are significantly stiffened and even become almost linearly elastic at high interfacial compliance. Even though the elastoplastic behaviors of the interphase are inversely obtained through the nonlinear micromechanics model and MD simulation results of nanocomposites, the key structural property relationship of the stiffened behavior of densely packed interphase can be explained via molecular level description using the concept of ‘free volume’. As the density of the interphase is much higher than that of pure amorphous PP polymer, the interphase region has less unoccupied space (free volume) to allow molecular motion of individual subunits of PP molecules according to the external loading than pure PP polymer, thus, the interphase can have higher elastic moduli and yield strength. A similar phenomenon of this density-related mechanical behavior of polymeric materials is the hydrostatic pressure-dependent yielding in compressive loading to which the same theoretical elucidation is applied. Once the elastoplastic behaviors of the interphase are obtained, the stress–strain behavior of nanocomposites at various volume fractions can be simply obtained from the micromechanics model itself. However, there is still no obvious answer to the question of which combination of the interfacial compliance and elastoplastic behavior is more realistic to predict the elastoplastic behavior of nanocomposites. In order to narrow down the possible combinations of interfacial compliance and interphase properties and to validate the applicability of the proposed multiscale model, we predicted the stress–strain curves of nanocomposites at small volume fractions and compared the result with some available experimental results and additional MD simulation results. The stress–strain curve of nanocomposites at 0.5 wt.% (weight fraction) are predicted from the proposed multiscale model and compared with the real experimental results of PP/MWNT (multi-walled nanotube) nanocomposites (Bao and Tjong, 2008) in Fig 9. The SEM image of the experimental nanocomposites microstructure is shown in Fig. 1(a) and the MWNTs are relatively well dispersed into the PP matrix without large clustering and waviness of CNTs. For the proposed multiscale

140 120

Stress(MPa)

100 80 60 α: α: α: α: α:

40 20 0

0

0.02

0.04

0.06

2.5nm/GPa 2.0nm/GPa 1.5nm/GPa 1.0nm/GPa 0.5nm/GPa

0.08

0.1

Strain Fig. 7. Stress–strain curves of the effective matrix obtained from the given elastoplastic behavior of 6.9% nanocomposites at various interfacial compliances.

141

S. Yang et al. / International Journal of Plasticity 41 (2013) 124–146

900 α: α: α: α: α:

800

Stress(MPa)

700 600

2.5nm/GPa 2.0nm/GPa 1.5nm/GPa 1.0nm/GPa 0.5nm/GPa

500 400 300 200 100 0 0

0.02

0.04

0.06

0.08

Strain Fig. 8. Stress–strain curves of interphase obtained from the given elastoplastic behavior of effective matrix and pure PP matrix.

45 40

(b) 14

Nanocomposites: Experiment W interphase, α=3nm/GPa W /O interphase, perfect bonding PP matrix: Experiment PP matrix: MD

σ composites - σ matrix (MPa)

(a) 50

Stress(MPa)

35 30 25 20 15 10

12

Experiment W interphase, α=3nm/GPa W /O interphase, perfect bonding

10 8 6 4 2

5 0 0

0.005

0.01

Strain

0.015

0.02

0 0

0.005

0.01

0.015

0.02

Strain

Fig. 9. Comparison of the proposed multiscale model with the experimental stress–strain curve of PP/MWNT(multi-walled nanotube) and original Mori– Tanaka model with perfect bonding condition while without interphase effect. (a) Direct comparison of stress–strain curves, (b) comparison of stressdifference between nanocomposites and PP matrix.

model, the interfacial compliance was set as 3 nm/GPa, and the resultant stress–strain curve of the interphase was considered in estimating the overall stress–strain curve of nanocomposites. As a reference, we predicted the stress-train curve at the same volume fraction using the original Mori–Tanaka model in which neither the interphase nor the weakened bonding effect are considered. Even though the stress–strain curve of PP matrix obtained from the present MD simulation differs slightly from the experimental results, the proposed model can accurately describe the reinforcing effect of CNT as shown in Fig. 9(a), while the original Mori–Tanaka model predicts less stiff behavior even while assuming perfect bonding condition. For a more clarified comparison, the differences of the stress between the nanocomposites and the pure matrix are compared in Fig. 9(b), and the proposed multiscale model prediction with the rigorous consideration of the interfacial imperfection and interphase are closer to the experimental results than the conventional Mori–Tanaka model predictions. Going back to the results in Fig. 6, where the effect of interfacial imperfection becomes critical at high volume fraction condition, it can be concluded that the relative dominance of the interfacial imperfection and the interphase varies according to the volume fractions. At small volume fractions, the interphase plays a more important role than the imperfect interface, while the trend is reversed at high volume fraction condition, thus, both effects balance each other. As an another validation to support the reinforcing mechanism, an additional MD simulation was performed using a 3% nanocomposites unit cell with the same procedure demonstrated in Section 3 to obtain isotropic stress–strain relation of nanocomposites. The stress–strain curves of nanocomposites are compared in Fig. 10. In this comparison, the same combination of the interfacial compliance and interphase behavior used in Fig. 9 is used for the present multiscale model. At a small strain range of less than 4%, which can be reasonably regarded as in the elastic region, the difference in stress–strain curves obtained from MD simulation, the present two-step nonlinear multiscale model, and the original Mori–Tanaka model

142

S. Yang et al. / International Journal of Plasticity 41 (2013) 124–146

600 MD simulation W interphase, α=3nm/GPa W /O interphase, perfect bonding

Stress(MPa)

500 400 300 200 100

0

0

0.02

0.04

0.06

0.08

0.1

Strain Fig. 10. Comparison of the stress–strain curve of nanocomposites obtained from the proposed two-step multiscale model with those from MD simulations and the Mori–Tanaka model with a perfect bonding condition.

are almost negligible. Above 4% strain, however, the gap between the MD simulation result and the prediction from the Mori–Tanaka model, which assumes the perfect bonding condition, becomes prominent. Even if the present model assigned a finite interfacial compliance of 3 nm/GPa, the result is closer to the MD simulation as a result of rigorous considerations of the interphase as an individual phase that contributes to the plastic behavior of nanocomposites. Integrating the results in Figs. 10 and 6, we can again deduce very important physics of the elastoplastic behavior of CNT-reinforced nanocomposites: that the interfacial imperfection and the formation of densely packed interphase plays important roll canceling each other. Because the relative dominance of the interfacial imperfection and interphase depends on the volume fractions, one can miss the important structure–property relationship of the reinforcing mechanism of CNTs unless the overall behavior of CNTreinforced nanocomposites is estimated at a wide range of volume fractions and evaluated via the comparison with wellestablished conventional micromechanics model predictions. 8. Concluding remarks In this study, a nonlinear two-step multiscale modeling approach to characterize elastic and plastic behavior of CNTreinforced PP nanocomposites was proposed with a rigorous consideration of the interfacial imperfection between the CNT and surrounding medium. and the interphase formed by the adsorbed polymer molecules of PP matrix to the surface of CNT. One important finding in the elastoplastic behavior of CNT/PP nanocomposites is that both the weakened interface and polymer-densified interphase play important roles, and they cancel each other out according to the volume fraction of CNT. From molecular dynamics simulations of well dispersed nanocomposite RVEs at several volume fractions, the transversely isotropic stress–strain curves are obtained and converted to isotropic stress–strain curves via a CNT orientation averaging scheme. As a nonlinear micromechanics model to describe the elastoplastic behavior of nanocomposites, we adopted the secant moduli method combined with the Mori–Tanaka mean field approach and the field fluctuation method. From the evaluation of the stress–strain relation by comparing the MD simulation results with the Mori–Tanaka model prediction with a perfect bonding condition, it has been found that the interface between CNT and PP is weak. At the same time, the existence of a dense interphase region made of PP in the vicinity of the CNT has been confirmed from the radial density distribution of the surrounding PP polymer. In order to consider the weakened interface and to characterize the contribution of the interphase to the overall elastoplastic behavior of nanocomposites, a two-step inverse analysis was suggested by sequentially decomposing nanocomposites microstructure incorporating the concept of effective matrix. In upper level analysis, the stress–strain relation of effective matrix at a given stress–strain relation of nanocomposites was numerically obtained in the frame work of Mori–Tanaka method with a finite value of interfacial compliance, followed by a lower level analysis that decomposes the effective matrix into the pure matrix and interphase using the same numerical scheme. According to the interfacial compliance considered in upper level analysis, several combinations of the interfacial compliance and corresponding stress–strain relation of the interphase could be successfully obtained. Using the elastoplastic behavior of the interphase with a finite interfacial compliance, the stress–strain relation of nanocomposites at very small volume fractions are obtained and compared with available experimental results and additional MD simulation results. With a rigorous consideration of the weakened interface and interphase, it has been proved that the present multiscale modeling approach can generate a quite reasonable estimation of the stress–strain relation of nanocomposites, while the direct application of original Mori–Tanaka method has some limitations in describing the balancing effect of the weakened interface and densely packed interphase.

S. Yang et al. / International Journal of Plasticity 41 (2013) 124–146

143

For design of nanocomposites to achieve excellent load transfer from matrix to CNTs or to other types of nanofillers, many kinds of surface treatments have become key to functionalizing the surface of nanocomposites. For evaluation of the newly developed surface functionalization, the present two-step multiscale modeling approach based on molecular level simulation and a continuum micromechanics model can be efficiently applied to characterize the enhancement of interfacial strength. In addition, the inverse scheme proposed in this study to identify the interfacial condition and interphase behavior can be extended to other types of rate-dependent mechanical behaviors such as viscoelastic (Jia et al., 2011; Tehrani et al., 2011) behavior of composites that employ CNTs or other types of nanofillers to take advantage of structural change induced by such fillers to maximize the reinforcing effect. As the elastic and plastic behaviors of polymeric materials have a hydrostatic pressure effect in compression (Ghorbel, 2008), consideration of such effect in the constitutive relation of polymer matrix is required for more generalized description of the plastic deformation of polymer nanocomposites. We will consider such effect in future work by performing compression simulation and defining the constitutive relation of polymer matrix as a function of the hydrostatic pressure to identify the anisotropic yield contour of polymer matrix, nanocomposites, and interphase. Acknowledgments This work was supported by the Hyundai and Kia Corporate Research and Development Division, WCU (Grant No. R312009-000-10083-0), funded by the MEST of the Korean government, and National Research Foundation of Korea (NRF) Grant funded by the Korea government (MEST) (No. 2012R1A3A2048841). Also, This research was kindly supported by a grant from the Seed Collaborative R&D Program funded by the Korean Research Council of Fundamental Science and Technology (KRCF), Republic of Korea. The authors gratefully acknowledge this support. Appendix A. Explicit formulae of the nonzero components of the compliance tensor Rijkl in Eq. (11) Following Qu’s original formulae, the fourth order interfacial compliance tensor Rijkl is expressed as a function of the interfacial compliance a and b as:

Rijkl ¼ aðPijkl  Q ijkl Þ þ bQ ijkl ;

ðA:1Þ

where the two fourth order tensors Pijkl and Qijkl are of integral form with respect to the surface normal vector of the spherical inclusion. The explicit form of the nonzero components of the tensor Pijkl and Qijkl, and their symmetry are as follows (Barai and Weng, 2011a):

) pffiffiffiffiffiffiffiffiffiffiffiffiffiffi a a2  1 1 1 ffiffiffiffiffiffiffiffiffiffiffiffiffiffi p ;  sin 2aða2  1Þ a a a2  1 ( ) pffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3 a a2  2 a2  1 a 1 ffiffiffiffiffiffiffiffiffiffiffiffiffiffi p ; ¼ þ sin 2a a 2ða2  1Þ a2  1 2ða2  1Þ ( ( )) pffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3 a a2  1 a2  2 a a2  1 1 1 1 pffiffiffiffiffiffiffiffiffiffiffiffiffiffi sin pffiffiffiffiffiffiffiffiffiffiffiffiffiffi sin ¼   ; 8a 2ða2  1Þ a a a a2  1 a2  1

P1111 ¼ P2222 P1212

3

P1313 ¼ P1212 ;

(

P3333 ¼ P 2222 ;

P2323 ¼

ðA:2Þ

P2222 ; 2

where a and a respectively indicate the diameter and aspect ratio (length to the diameter) of embedded inclusion.

Q 1111 Q 2222 Q 1122

( pffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffi) 3 a a2  1 3 a2  1 1 1 pffiffiffiffiffiffiffiffiffiffiffiffiffiffi sin ; ¼  sin 5=2 2 2 2a a a a 1 ða  1Þ ( ) pffiffiffiffiffiffiffiffiffiffiffiffiffiffi 9 a3 ð2a2  5Þ a3 ða2  4Þ a2  1 1 ; ¼ a sin þ 5=2 2 8a 2ða2  1Þ a 2ða  1Þ ( ) pffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3 aða2 þ 2Þ a2  1 aða2  4Þ 1 ¼ sin þ ; 4a ða2  1Þ5=2 a ða2  1Þ2

Q 3333 ¼ Q 2222 ;

Q 1212 ¼ Q 1313 ¼ Q 1122 ¼ Q 1133 ;

Q 2233 ¼ Q 2323 ¼

ðA:3Þ

Q 2222 : 3

Appendix B. Complete expression of constants nA ; nCA ; nRCA and gA ; gCA ; gRCA for randomly oriented nanocomposites Following the tensor manipulation of transversely isotropic fiber, such as carbon nanotube-reinforced composites with Hill’s notation of the transversely isotropic stiffness tensor (Hill, 1964; Qiu and Weng, 1990), the complete expressions of scalar components in Eq. (20) are given as follows:

144

S. Yang et al. / International Journal of Plasticity 41 (2013) 124–146

nA ¼ nCA ¼ nRCA

1

ðpÞ

ðpÞ

9l 1

½2d

ðpÞ

 2ðg ðpÞ þ h Þ þ cðpÞ ; ðpÞ

ðpÞ

ðpÞ

ðpÞ

½4ðkp d  lp g ðpÞ Þ þ 2ðlp d  np g ðpÞ þ lp cðpÞ  2kp h Þ þ nr cðpÞ  2lp h ; ðpÞ 9l" ðpÞ ðpÞ 1 4ðR2222 þ R2233 Þðkp d  lp g ðpÞ Þ 4R2211 ðlp d  np g ðpÞ Þ ¼ þ ðpÞ ðpÞ 9 l l ðpÞ

2R1122 ðkp d

þ2

 lp g ðpÞ Þ

ðpÞ

l

ðpÞ

þ

R1111 ðnp cðpÞ  2lp h Þ ðpÞ

l

ðpÞ

þ

R1111 ðlp d

l

 np g ðpÞ Þ

ðpÞ ðpÞ

þ

2R1122 ðlp cðpÞ  2kp h Þ l

ðpÞ

þ

R2211 ðnp cðpÞ  2lp h Þ l

#

ðpÞ

ðpÞ

þ

! ðpÞ ðR2222 þ R2233 Þðlp cðpÞ  2kp h Þ l

ðpÞ

ðB:1Þ

;

and

gA ¼ CA

g ¼ gRCA

1

ðpÞ

½d ðpÞ

30l 1

h

ðpÞ

ðpÞ

þ 2ðg ðpÞ þ h Þ þ 2cðpÞ  þ ðpÞ

ðkp d

ðpÞ

 lp g ðpÞ Þ  ðlp d

  1 1 1 ; þ 5 eðpÞ f ðpÞ

 i 2 m pp p ðpÞ ðpÞ ;  np g ðpÞ þ lp cðpÞ  2kp h Þ þ ðnp cðpÞ  2lp h Þ þ þ 5 eðpÞ f ðpÞ

15l" ðpÞ ðpÞ 1 ðR2222 þ R2233 Þðkp d  lp g ðpÞ Þ R2211 ðlp d  np g ðpÞ Þ ¼ þ ðpÞ ðpÞ 15 l l ðpÞ

2R1122 ðkp d



l

 lp g ðpÞ Þ

ðpÞ

ðpÞ

þ

R1111 ðlp d

l

 np g ðpÞ Þ

ðpÞ

ðpÞ

þ

R2211 ðnp cðpÞ  2lp h Þ l

ðpÞ

þ

! ðpÞ ðR2222 þ R2233 Þðlp cðpÞ  2kp h Þ

 # R1111 ðnp cðpÞ  2lp h Þ 2R1122 ðlp cðpÞ  2kp h Þ R2323 mp R1212 pp ; þ þ 12 þ þ ðpÞ ðpÞ eðpÞ f ðpÞ l l ðpÞ

l

ðpÞ

ðpÞ

ðB:2Þ

where k; l; n; m; p, with subscripts ‘m’ and ‘p’ corresponding to the matrix phase and inclusion phase (CNT in this study), are the five elastic constants of transversely isotropic medium given in Eq. (4). In applying Hill’s notation for transversely isotropic stiffness to an isotropic matrix, the five elastic constants are correlated with the isotropic shear and bulk moduli as km = jm + lm/3, lm = jm  2lm/3, nm = jm + 4lm/3, and mm = pm = lm. Together with the transversely isotropic stiffness components of individual phase, six transversely isotropic parameters with a bracket superscript ‘(p)’are defined as:

2ðkp  km Þ 2ðlp  lm Þ ½ð1  mm ÞðS2222 þ S2233 Þ  2mm S2211  þ ½S2211  mm ðS2222 þ S2233 Þ; Em Em ðkp  km Þ ðlp  lm Þ ½ð1  mm ÞS1122  mm S1111  þ ½S1111  2mm S1122 ; g ðpÞ ¼ 2 Em Em ðnp  nm Þ ðlp  lm Þ ðpÞ ½S2211  mm ðS2222 þ S2233 Þ þ ½ð1  mm ÞðS2222 þ S2233 Þ  2mm S2211 ; h ¼ Em Em ðnp  nm Þ 2ðlp  lm Þ ðpÞ ½S1111  2mm S1122  þ ½ð1  mm ÞS1122  mm S1111 ; d ¼1þ Em Em 2ðmp  mm Þ S2323 ; eðpÞ ¼ 1 þ mm 2ðpp  pm Þ f ðpÞ ¼ 1 þ S1212 ; pm

cðpÞ ¼ 1 þ

l

ðpÞ

ðpÞ

¼ cðpÞ d

 2g ðpÞ h

ðB:3Þ

ðpÞ

where Em and mm are the Young’s modulus and Poisson’s ratio of matrix phase, respectively. For multi-particle-reinforced cases, the subscript ‘p’ and the superscript ‘(p)’ are expanded over the whole number of inclusions having different elastic stiffness. One can refer to (Qiu and Weng, 1990; Barai and Weng, 2011a) for detailed derivation of Eqs. (B.1)–(B.3). References Accelrys Inc. San Francisco. Ajayan, P.M., Stephan, O., Colliex, C., Trauth, D., 1994. Aligned carbon nanotube arrays formed by cutting a polymer resin-nanotube composite. Science 265, 1212–1214. Bao, S.P., Tjong, S.C., 2008. Mechanical behaviors of polypropylene/carbon nanotube nanocomposites: the effects of loading rate and temperature. Mater. Sci. Eng. A 485, 508–516. Barai, P., Weng, G.J., 2011a. A theory of plasticity for carbon nanotube reinforced composites. Int. J. Plasticity 27, 539–559. Barai, P., Weng, G.J., 2011b. Mechanics of nanocrystalline coating and grain-size dependence of its plastic strength. J. Mech. Mater. Struct. 43 (9), 496–504. Berveiller, M., Zaoui, A., 1979. An extension of the self-consistent scheme to plastically-flowing polycrystals. J. Mech. Phys. Solids 26, 325–344.

S. Yang et al. / International Journal of Plasticity 41 (2013) 124–146

145

Bhattacharyya, A.R., Sreekumar, T.V., Liu, T., Kumar, S., Ericson, L.M., Hauge, R.H., Smalley, R.E., 2003. Crystallization and orientation studies in polypropylene/single wall carbon nanotube composite. Polymer 44, 2373–2377. Bobeth, M., Dienner, G., 1986. Field fluctuation in multicomponent mixtures. J. Mech. Phys. Solids 34, 1–17. Cho, M., Yang, S., Chang, S., Yu, S., 2011. 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. 85 (12), 1564–1583. Eshelby, J.D., 1957. The determination of the elastic field of an ellipsoidal inclusion and related problems. Proc. R. Soc. Lond. A 241, 376–396. Esteva, M., Spanos, P.D., 2009. Effective elastic properties of nanotube reinforced composites with slightly weakened interfaces. J. Mech. Mater. Struct. 4, 887–900. Fisher, F.T., Brinson, L.C.C., 2001. Viscoelastic interphases in polymer-matrix composites: theoretical models and finite-element analysis. Compos. Sci. Technol. 61, 731–748. Fisher, F.T., Bradshaw, R.D., Brinson, L.C., 2003. Fiber waviness in nanotube-reinforced polymer composites-I: modulus predictions using effective nanotube properties. Compos. Sci. Technol. 63, 1689–1703. Fotiu, P.A., Nemat-Nasser, S., 1996. Overall properties of elastic-viscoplastic periodic composites. Int. J. Plasticity 12, 163–190. Ghorbel, E., 2008. A viscoplastic constitutive model for polymeric materials. Int. J. Plasticity 24, 2032–2058. Gojny, F.H., Wichmann, M.H.G., Kopke, U., Fiedler, B., Shulte, K., 2004. Carbon nanotube-reinforced epoxy-composites: enhanced stiffness and fracture toughness at low nanotube content. Compos. Sci. Technol. 64 (15), 2363–2371. Granb, M., Satapathy, B.K., Thunga, M., Weidisch, R., Potschke, P., Jehnichen, D., 2008. Structural interpretations of deformation and fracture behavior of polypropylene/multi-walled carbon nanotube composites. Acta Mater. 56, 2247–2261. Han, Y., Elliott, J., 2007. Molecular dynamics simulations of the elastic properties of polymer/carbon nanotube composites. Compos. Mater. Sci. 39, 315–323. Hill, R., 1964. Theory of mechanical properties of fiber-strengthened materials: I. Elastic behavior. J. Mech. Phys. Solids 12, 199–212. Hill, R., 1965. Continuum micro-mechanics of elastoplastic polycrystals. J. Mech. Phys. Solids 13, 89–101. Hoover, W.G., 1985. Canonical dynamics: equilibrium phase-space distributions. Phys. Rev. A 31, 1695–1697. Hoover, W.G., 1986. Constant-pressure equations of motion. Phys. Rev. A 34, 2499–2500. Hu, G.K., 1996. A method of plasticity for general aligned spheroidal void or fiber-reinforced composites. Int. J. Plasticity 12, 439–449. Hu, G.K., 1997. Composites Plasticity based on matrix average second order stress moment. Int. J. Solids. Struct. 34, 1007–1015. Hutchinson, J.W., 1976. Bounds and self-consistent estimates for creep of polycrystalline materials. Proc. R. Soc. Lond. A 348, 101–127. Jia, Y., Peng, K., Gong, X., Zhang, Z., 2011. Creep and recovery of polypropylene/carbon nanotube composites. Int. J. Plasticity 27, 1239–1251. Jin, L., Bower, C., Zhou, O., 1998. Alignment of carbon nanotubes in a polymer matrix by mechanical stretching. Appl. Phys. Lett. 73, 1197–1199. Kreher, W., Pompe, W., 1989. Internal Stress in Heterogeneous Solids. Akademie, Berlin. Li, L., Li, C.Y., Ni, C., Rong, L., Hsiao, B., 2007. Structure and crystallization behaviors of nylon 66/multi-walled carbon nanotube nanocomposites at low carbon nanotube contents. Polymer 48, 3452–3460. Liew, K.M., He, X.Q., Wong, C.H., 2004. On the study of elastic and plastic properties of multi-walled carbon nanotubes under axial tension using molecular dynamics simulation. Acta Mater. 52, 2521–2527. Mathur, R.B., Chatterjee, S., Singh, B.P., 2008. Growth of carbon nanotubes on carbon fiber substrates to produce hybrid/phenolic composites with improved mechanical properties. Compos. Sci. Technol. 68, 1608–1615. Mercier, S., Molinari, A., 2009. Homogenization of elastic–viscoplastic heterogeneous materials: self-consistent and Mori–Tanaka schemes. Int. J. Plasticity 25, 1024–1048. Mesbah, A., Zairi, F., Boutaleb, S., Gloaguen, J.M., Nait-Abdelaziz, M., Xie, S., Boukharouba, T., Lefebvre, J.M., 2009. Experimental characterization and modeling stiffness of polymer/clay nanocomposites within a hierarchical multiscale framework. J. Appl. Poly. Sci. 114, 3274–3291. Mori, T., Tanaka, K., 1973. Average stress in matrix and average elastic energy of materials with misfitting inclusions. Acta. Metall. 21, 571–574. Odegard, G.M., Gates, T.S., Wise, K.E., Park, C., Siochi, E.J., 2003. Constitutive modeling of nanotube-reinforced polymer composites. Compos. Sci. Technol. 63, 1671–1687. Omidi, M., Rokni, D.T.H., Milani, A.S., Seethaler, R.J., Arasteh, R., 2010. Prediction of the mechanical characteristics of multi-walled carbon nanotube/epoxy composites using a new form of the rule of mixtures. Carbon 48, 3218–3228. Papakonstantopoulos, G.J., Yoshimoto, K., Doxastakis, M., Nealey, P.F., de Pablo, J.J., 2005. Local mechanical properties of polymeric nanocomposites. Phys. Rev. E 72, 031801-1–031801-6. Pierard, O., Doghri, I., 2006. An enhanced affine formulation and the corresponding numerical algorithms for the mean-field homogenization of elasto– viscoplastic composites. Int. J. Plasticity 22, 131–157. Plimpton, S., 1995. Fast parallel algorithms for short-range molecular dynamics. J. Comput. Phys. 117, 1–19. Qian, D., Dickey, E.C., Andrews, R., Rantell, T., 2000. Load transfer and deformation mechanisms in carbon nanotube-polystyrene composites. Appl. Phys. Lett. 76 (20), 2868–2870. Qiu, Y.P., Weng, G.J., 1990. On the application of Mori–Tanaka theory involving transversely isotropic spheroidal inclusions. Int. J. Eng. Sci. 28 (11), 1121– 1137. Qiu, Y.P., Weng, G.J., 1992. A theory of plasticity for porous materials and particle-reinforced composites. J. Appl. Mech. 59, 261–268. Qu, J., 1993. Eshelby tensor for and elastic inclusion with slightly weakened interface. J. Appl. Mech. 60, 1048–1050. Rottler, J., Robbins, M.D., 2001. Yield conditions for deformation of amorphous polymer glasses. Phys. Rev. E 64 (5), 051801-1–051801-8. Seidel, G.D., Lagoudas, D.C., 2006. Micromechanical analysis of the effective elastic properties of carbon nanotube reinforced composites. J. Mech. Mater. Struct. 38 (8–10), 884–907. Seidel, G.D., Lagoudas, D.C., 2008. A micromechanics model for the thermal conductivity of nanotube-polymer nanocomposites. J. Appl. Mech. 75 (4), 041025-1–041025-9. Seidel, G.D., Lagoudas, D.C., 2009. A micromechanics model for the electrical conductivity of nanotube-polymer nanocomposites. J. Compos. Mater. 43 (9), 917–941. Shi, D.L., Feng, X.Q., Huang, Y.Y., Hwang, K.C., Gao, H., 2004. The effect of nanotube waviness and agglomeration on the elastic property of carbon nanotubereinforced composites. J. Eng. Mater. Tech. 126, 250–257. Sinnott, Susan B., 2002. Chemical functionalization of carbon nanotubes. J. Nanosci. Nanotechnol. 2 (2), 113–123. Spitalsky, Z., Tasis, D., Papagelis, K., Galiotis, C., 2010. Carbon nanotube-polymer composites: chemistry, processing, mechanical and elelctrical properties. Progress. Poly. Sci. 35, 357–401. Sun, H., Mumby, S.J., Maple, J.R., Hagler, A.T., 1994. An abi initio CFF90 all-atom force field for polycarbonates. J. Am. Chem. Soc. 116, 2978–2987. Suquet, P., 1995. Overall properties of nonlinear composites: a modified secant moduli theory and its link with Ponte Castañeda’s nonlinear variational procedure. C. R. Acad. Des. Sci. 320, 563–571. Talbot, D.R.S., Willis, J.R., 1985. Variational principles for inhomogeneous nonlinear media. IMA J. Appl. Math. 35, 39–54. Tandon, G.P., Weng, G.J., 1988. A theory of particle-reinforced plasticity. J. Appl. Mech. 55, 126–135. Tehrani, M., Safdari, M., Al-Haik, M.S., 2011. Nanocharacterization of creep behavior of multiwall carbon nanotubes/epoxy nanocomposite. Int. J. Plasticity 27, 887–901. Thostenson, E.T., Li, C., Chou, T.W., 2005. Nanocomposites in context. Compos. Sci. Technol. 65, 491–516. Vaisman, L., Marom, G., Wagner, H.D., 2006. Dispersions of surface-modified carbon nanotubes in water-soluble and water-insoluble polymers. Adv. Funct. Mater. 16, 357–363. Walpole, L.J., 1981. Elastic behavior of composite materials: theoretical foundations. Adv. Appl. Mech. 21, 169–241. Watts, P.C.P., Hsu, W.K., 2003. Behaviors of embedded carbon nanotubes during film cracking. Nanotechnology 14, L7–L10.

146

S. Yang et al. / International Journal of Plasticity 41 (2013) 124–146

Wei, C.Y., 2007. Structural phase transition of alkane molecules in nanotube composites. Phys. Rev. B 76, 134104-1–134104-10. Wei, C.Y., Srivastava, D., Cho, K.J., 2004. Structural ordering in nanotube polymer composites. Nano Lett. 4 (10), 1494–1952. Yanase, K., Ju, J.W., 2010. The modified Eshelby tensor and effective properties of spherical particle reinforced composites. In: Proceedings of the Ninth World Congress on Computational Mechanics and the Fourth Asian Pacific Congress on Computational Mechanics. Yanase, K., Ju, J.W., 2011. Effective elastic moduli of spherical particle reinforced composites containing imperfect interfaces. Int J. Damage Mech. 21 (1), 97– 127. Yang, S., Cho, M., 2008. Scale bridging method to characterize mechanical properties of nanoparticle/polymer nanocomposites. Appl. Phys. Lett. 93 (4), 043111-1–043111-3. Yang, S., Cho, M., 2009. A scale bridging method for nanoparticulate polymer nanocomposites aned their non-dilute concentration effect. Appl. Phys. Lett. 94, 223104-1–223104-3. Yang, S., Yu, S., Cho, M., 2010. Sequential thermoelastic multiscale analysis of nanoparticulate composites. J. Appl. Phys. 108 (5), 056102-1–056102-3. Yang, Q.S., He, X.Q., Liu, X., Leng, F.F., Mai, Y.W., 2012a. The effective properties and local aggregation effect of CNT/SMP composites. Compos. Part B: Eng. 43 (1), 33–38. Yang, S., Yu, S., Kyoung, W., Hah, D.S., Cho, M., 2012b. Multiscale modeling of size-dependent elastic properties of carbon nanotube/polymer nanocomposites with interfacial imperfections. Polymer 53 (2), 623–633. Yu, S., Yang, S., Cho, M., 2009. Multi-scale modeling of cross-linked epoxy nanocomposites. Polymer 50 (3), 945–952. Yu, S., Yang, S., Cho, M., 2011. Multiscale modeling of cross-linked epoxy nanocomposites to characterize the effect of particle size on thermal conductivity. J Appl. Phys. 110, 124302-1–124302-9. Zhang, W.X., Wang, T.J., Chen, X., 2010. Effect of surface/interface stress on the plastic deformation of nanoporous materials and nanocomposites. Int. J. Plasticity 26, 957–975.