Stochastic predictions of interfacial characteristic of polymeric nanocomposites (PNCs)

Composites: Part B 59 (2014) 80–95

Contents lists available at ScienceDirect

Composites: Part B journal homepage: www.elsevier.com/locate/compositesb

Stochastic predictions of interfacial characteristic of polymeric nanocomposites (PNCs) N. Vu-Bac a,b, T. Lahmer a, Y. Zhang a, X. Zhuang b,⇑, T. Rabczuk a,c,d,⇑ a

Institute of Structural Mechanics, Bauhaus-Universität Weimar, Marienstr. 15, D-99423 Weimar, Germany Department of Geotechnical Engineering, Tongji University, Shanghai 200092, China c Graduiertenkolleg 1462, Bauhaus-Universität Weimar, Berkaerstr. 9, D-99423 Weimar, Germany d School of Civil, Environmental and Architectural Engineering, Korea University, Republic of Korea b

a r t i c l e

i n f o

Article history: Received 7 July 2013 Received in revised form 30 October 2013 Accepted 12 November 2013 Available online 23 November 2013 Keywords: A. Nano-structures A. Polymer–matrix composites (PMCs) B. Interface/interphase C. Computational modeling Stochastic prediction

a b s t r a c t The effect of the single-walled carbon nanotube (SWCNT) radius, the temperature and the pulling velocity on interfacial shear stress (ISS) is studied by using the molecular dynamics (MD) simulations. Based on our MD results, the mechanical output (ISS) is best characterized by the statistical Weibull distribution. Further, we also quantify the influence of the uncertain input parameters on the predicted ISS via sensitivity analysis (SA). First, partial derivatives in the context of averaged local SA are computed. For computational efficiency, the SA is based on surrogate models (polynomial regression, moving least squares (MLS) and hybrid of quadratic polynomial and MLS regressions). Next, the elementary effects are determined on the mechanical model to identify the important parameters in the context of averaged local SA. Finally, the approaches for ranking of variables (SA based on coefficients of determination) and variancebased methods are carried out based on the surrogate model in order to quantify the global SA. All stochastic methods predict that the key parameters influencing the ISS is the SWCNT radius followed by the temperature and pulling velocity, respectively. Ó 2013 Elsevier Ltd. All rights reserved.

1. Introduction The Carbon nanotubes (CNTs) are considered as an ideal reinforcement for polymer composites due to their exceptional mechanical and electrical properties [1–3] and low density [4–6]. By adding CNTs into polymer, the properties of the resulting nanocomposite material such as strength and modulus are enhanced. However, the enhanced properties of these composites are strongly affected by the mechanics that govern the interface between CNT and polymer matrix [7]. The mechanical properties of composite materials greatly depend on the load transfer mechanism through the interface between the polymer matrix and CNTs and the strength of the interface. Therefore, the interface plays a significant role in the load transfer and the consequent improvements in modulus and strength. For conventional fiber-reinforced polymer composites, the single fiber pull-out test [8] is typically used to evaluate the interfacial shear strength. Great experimental efforts have been carried out in order to investigate interfacial properties.

⇑ Corresponding authors. Addresses: Department of Geotechnical Engineering, Tongji University, 1239 Siping Road, Shanghai, China (X. Zhuang), Institute of Structural Mechanics, Bauhaus-Universität Weimar, Marienstr. 15, D-99423 Weimar, Germany (T. Rabczuk). Tel.: +49 (0)3643 58 4511. E-mail addresses: [email protected] (X. Zhuang), timon.rabczuk@ uni-weimar.de (T. Rabczuk). 1359-8368/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compositesb.2013.11.014

For instance, the direct CNT pull-out experiments were made in transmission electron microscopy (TEM) [9,10] or by atomic force microscopy (AFM) [11,12], Raman spectroscopy [4], scanning probe microscope (SPM) [13,14]. Experiments in [4,5,9] demonstrated that since making atomic bonding between CNTs and polymer matrix is difficult, CNTs often interact mainly with the polymer matrix through van der Waals (vdW) forces. In the absence of atomic bonding (cross-links) between the CNT and the polymer matrix, the vdW forces especially govern the load transfer capability of the CNT/polymer interface. Due to difficulty of nanomanipulation in experiments, as mentioned by [15], numerical simulation is employed as a powerful alternative approach to study the CNT/polymer interface. Many molecular simulations have been carried out in order to investigate the CNT/polymer interfacial characteristics. In [16–20], the research presented molecular mechanics (MM) and molecular dynamics (MD) simulations, in which the authors examined the interfacial shear stress (ISS) during the entire pull-out process of CNT from polymer matrix. Wei et al. [21] used MD simulations to study the thermal expansion and diffusion characteristics of the single-walled carbon nanotube (SWCNT)/polyethylene (PE) composite. He found that insertion of SWCNT into a polymer matrix increases the glass transition temperature T g . For SWCNT/PE composites with long chained polymer molecules, the glass transition temperature T g

81

N. Vu-Bac et al. / Composites: Part B 59 (2014) 80–95

is around 400 K. He also studied the temperature dependent adhesion behavior of CNT/PE composites using MD simulations for tensile deformation of CNT/PE composites and reported that the interfacial shear stress between SWCNT and PE matrix decreases with the increase of temperature [22]. Frankland et al. [17] developed an interfacial friction model depending on the pull-out force, an effective viscosity and strain rate. Based on this model, a linear force-sliding velocity relation is used to estimate the average interfacial interaction. Zhang et al. [23] studied the rate-dependent interfacial behavior between SWCNT and PE in the absence of cross-links by using MD simulations. They indicated that the ISS increases with an increase in the sliding velocity of CNT. They also investigated the effects of the SWCNT radius on the ISS and reported that the ISS decreases when the SWCNT radius increases. Li et al. also showed that SWCNT radius have effect on the ISS, but the SWCNT length is independent of the ISS [24]. They have studied the interfacial properties between SWCNT and polymer matrix. However, the pull-out simulations based on MM simulations were carried out without considering the deformation of the polymer matrix during the pull-out process by fixing the matrix. While the influence of certain parameters on the ISS of CNT/ polymer composite has been obtained qualitatively in most of literatures, quantitative results in the context of stochastic analysis has not yet been investigated. In this paper, a SA is performed in order to quantify the influence of the SWCNT radius, temperature and pulling velocity on the ISS. Therefore, we compute first-order, total-effect sensitivity indices and SA based on coefficients of determination (COD) in the context of a global SA. Moreover, we approximate global SA by performing a series of local SA where we calculate partial derivatives, elementary effects and average the results. It should be noted that SA is implemented on the average ISS as presented in [16,18,19]. The variation in the distributions of ISS along the SWCNT length is not within the scope of this study. The paper is outlined as follows. In the next section, we briefly describe the MD model for SWCNT/polymer composites. The results of the MD model will be shown in Section 3. Section 4 describes the SA. In Section 5, the numerical regression models, the COD and the sensitivity indices will be presented before we discuss the numerical results. Finally, we close the manuscript with concluding remarks and future work.

2. Molecular dynamics model 2.1. Potential functions and parameters The adaptive intermolecular reactive empirical bond order (AIREBO) potential in [25,26] is adopted to consider the intra-carbon nanotube interactions. In this work, PE is chosen as a model matrix. A united-atom model is used in which each CH2 group is lumped into a single site to represent the molecular structure of PE [27]. The functional form

and the parameters of the intra-polymer potential [28] are outlined in Table 1. A non-bonded vdW force of the truncated Lennard–Jones (LJ) 6– 12 with  ¼ 0:1102 kcal/mol and r ¼ 3:65 Å is used to describe the interaction between the SWCNT and the PE matrix, see [21]. For the bulk PE, samples consist of 20 polymer chains with 500units (united atoms) in each chain. The SWCNT/polymer composite samples with periodic boundary condition in x, y directions are composed of a non-covalent SWCNT and an amorphous PE matrix (10,000 units of CH2 ) that was created by implementing a Monte Carlo random walk growth algorithm [29]. The simulation box after the equilibration process is a cubic space as shown in Fig. 1. LAMMPS [26] was used to perform the equilibration process. The equilibration involves four different steps; [30]. 1. Each generated structure is equilibrated for 100 ps (Dt ¼ 1  104 ps) so that the volume and energy of the system becomes stable. 2. The system is kept at a temperature of T ¼ 500 K and pressure P = 1 atm for 500 ps (Dt ¼ 5  104 ps) in the NPT ensemble 3. The system is cooled down to the desired temperature at a cooling rate of 0.8 K/ps. 4. The system is further equilibrated for (Dt ¼ 5  104 ps) at the desired temperature in the NPT ensemble. 2.2. Steered molecular dynamics Steered molecular dynamics (SMD) simulations are employed to simulate the pull-out of the SWCNT from the PE matrix. The potentials of mean force (PMF) (equilibrium property) are calculated from nonequilibrium processes (such as SMD simulations). By applying a moving spring force to the center of mass of the SWCNT’s atoms, the SWCNT is pulled along its axial (z-axis) direction at a constant speed while all longitudinal sides of the PE matrix is constrained as shown in Fig. 2. The restoring force of magnitude is given by [31]:

FðtÞ ¼ K spring





vspring ðtÞ  vpull ðtÞ

ð1Þ 2

where the elastic spring (Kspring ¼ 100 eV=Å ) is connected between the tether point and the center of mass of the SWCNT; vspring and vpull the spring and pulled group positions, respectively. For constant velocity pulling, the integral over time of the force projected on the pulling direction is accumulated and then used to compute the PMF by averaging over multiple independent trajectories along the same pulling path. More details can be found in [32]. 2.3. Pull-out simulations Fig. 3 illustrates snap-shots of the pull-out process. The PMF (pull-out energy) is computed and recorded during the pulling process. As shown in Fig. 4, the total potential energy of SWCNT/PE system and pull-out energy increases linearly with displacement

Table 1 Functional form and parameters of the force field of united atom polyethylene. Interaction

Form

Bond

Ebond ¼ 12 kb ðl  leq Þ

Angle

Eangle ¼ 12 kh ðcosðhÞ  cosðheq ÞÞ2 P Etorsional ¼ 12 3i¼0 kn cosn ð/Þ ( h  i 12 4 r  ðr Þ6 r 6 rc Enon-bond ¼ 0 r > rc

Torsional Non-bonded

Parameters 2

2

kb ¼ 350 kcal=mol Å ; leq ¼ 1:53 Å

2

2

kh ¼ 60 kcal=mol=rad ; heq ¼ 109:5 k0 ¼ 1:736; k1 ¼ 4:490; k2 ¼ 0:776; k3 ¼ 6:990ðkcal=molÞ

 ¼ 0:112 kcal=mol;r ¼ 4:01 Å; r c ¼ 10:5 Å

82

N. Vu-Bac et al. / Composites: Part B 59 (2014) 80–95

Fig. 10 illustrates the stress–strain response for glassy amorphous PE system for different temperatures and different strain rates. For the sake of comparision, we have simulated the PE system with 10 chains/1000-units. The results are in good agree with those in [30,28]. 3. Molecular dynamics simulation results

Fig. 1. Cross section of SWCNT/PE composite.

Fig. 2. Pull-out SWCNT from PE matrix using SMD simulations. The green boxes show the fix edges. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

during the pull-out process. This trend is a good agreement with previous results in [16,18,24]. The ISS can be calculated from the pull-out energy by [16,18]:

Epullout ¼

Z

L

2prðL  xÞsi dx ¼ pr si L2

Normal distribution is typically assumed to approximate the CNT’s radius distribution [33]. A mean diameter of 1.42 nm, and approximately 78% of the nanotubes ranging from 1.0 to 1.8 nm in diameter. As shown in Fig. 11, an increase in the SWCNT radius results in a decrease in the PE region of the interaction with a carbon atom on the SWCNT [23]. Therefore, the ISS decreases with increasing SWCNT radius. The interaction between a large SWCNT radius with PE matrix can be considered the interaction of graphite sheet with PE so that the SWCNT (16, 16) (radius = 10.848 Å) is chosen as upper bound and the SWCNT (4, 4) (radius = 2.712 Å) as lower bound of truncated normal distribution in our MD simulations. For the sake of estimating SA, the armchair and zigzag SWCNT are chosen to ensure the SWCNT length to be nearly the same for all of samples. The density of both the bulk PE and the composite is a function of temperature, see Wei et al. [21]. A decrease in these densities with an increase in the temperature indicate the thermal expansions of the materials. The ISS decreases with increasing temperature, as expected. The composite behaves as a glass solid at low temperatures T g < 400 K for long-chained PE, and its mechanical response is much weakened at high temperatures. Therefore, the temperature should be sufficient below the glass transition temperature and approximate the glass transition T g . Therefore, the temperature from 100 K to 400 K assuming uniform distribution is used in our study. Izrailevs et al. [34] indicated that the motion of the SWCNT proceeds in the strong friction with surrounding matrix. The average applied force measures the local slope of the binding potential plus a frictional contribution that depends linearly on the pulling velocity. Therefore, the pull-out energy and the ISS depend on the pulling velocity as the average applied force depends on the pulling velocity. More details can be found in [34]. An interfacial friction model developed by Frankland et al. [17] shows that the average SWCNT velocity is linearly related to the average applied force. Hence, the pulling velocity influences on the ISS. The constant pull-out velocities ranging from 0.1 ps/Å to 1.0 ps/Å are examined in our MD simulations to reduce the computational cost in pull-out simulations as well as avoid the influence of molecular thermal vibration as mentioned in [17,23]. We assume a uniform distribution. The random input variables are listed in Table 2.

ð2Þ

0

4. Sensitivity analysis

si ¼

Epullout

prL2

ð3Þ

where r and L are the radius and embedded length of the SWCNT, respectively; x is the displacement of the SWCNT and si is the ISS. Note that this formula gives average ISS. The pull-out energy of SWCNT from the polymer matrix with different radii, temperatures and pulling velocities are illustrated in Figs. 5–7, respectively. Fig. 8 shows that the ISS is independent of the simulation box size and fix-edge effect. When the chain length ranges from 100-units to 1000-units, the density of polymer matrix is 0:87—0:91 g=cm3 [30,28]. Fig. 9 shows that the chain length is independent of the ISS. In our simulations, the chain length of 500-units is used.

4.1. Stochastic modeling The purpose of quantitative uncertainty analysis is to use currently available information in order to quantify the degree of confidence in the existing simulations containing a set of assumptions and conceived models that are used to predict the real behavior response of material. Latin Hypercube Sampling (LHS) [35,36] is an improved sampling strategy that enables a reliable approximation of the stochastic properties even for a small number of samples N. LHS will be used to provide the design points which are spread throughout the design space. The base of the method is the subdivision of the design space of a variable X i into N classes Dm of equal probability

83

N. Vu-Bac et al. / Composites: Part B 59 (2014) 80–95

Fig. 3. Snap shots of CNT pull-out from SWCNT/PE composite.

SWCNT/PE system, 100 K, 0.1 Angs/ps −4

−8270

5

−8280

4.8

x 10

−8290 −8300 −8310 0

10

20

30

40

50

60

70

80

Displacement (Angs)

ISS (eV/Angs3)

Total potential energy (eV)

SWCNT(10,10)/PE system, 100 K, 0.1 Angs/ps −8260

Pull−out energy (eV)

40

4.6

4.4

4.2

4

30 3

4

5

6

7

8

9

10

SWCNT radius (Angs)

20

Fig. 5. The ISS for SWCNT (5, 5), SWCNT (10, 10) and SWCNT (15, 15) at 100 K and a pulling velocity of 0.1 Å/ps.

10

0 0

10

20

30

40

50

60

70

80

Displacement (Angs) Fig. 4. Energy variation during the pull-out of the SWCNT: (top) Total potential energy, and (bottom) Pull-out energy (PMF).

SWCNT(10,10)/PE system, 0.1 Angs/ps −4

4.15

P½X i 2 Dm  ¼

1 ; N

i ¼ 1; . . . ; k;

m ¼ 1; . . . ; N

ð4Þ

x 10

4.1

ISS (eV/Angs3)

4.05

where N and k are the number of samples and the number of input parameters, respectively. These methods can reach a sufficient accuracy by estimating statistical properties with only a few samples compared to plain Monte-Carlo sampling [37].

4 3.95 3.9 3.85

4.2. Sensitivity analysis According to [38], SA in general is the study of how much model output values are affected by changes in model input values. In the following, a general description of local and variance based SA is given, which is usually exploited to rank the models input parameters and their contribution to the model output.

3.8

100

150

200

250

300

Temperature (K) Fig. 6. The ISS at a pulling velocity of 0.1 Å/ps for the temperature of 100 K, 200 K and 300 K.

84

N. Vu-Bac et al. / Composites: Part B 59 (2014) 80–95

−4

fect of the input on the output value. A scatter plot can draw the correlation between Y and X i , or even by estimating variancebased measures by nonlinear regression [39].

SWCNT(10,10)/PE system, 100 K

x 10

4.5 4.45

4.4. Derivative-based approach

ISS (eV/Angs3)

4.4

Assume that the model output Y is a scalar function of a vector X ¼ ðX 1 ; X 2 ; . . . ; X k Þ of k input parameters. Partial derivatives can be seen as a mathematical definition of the sensitivity of the model output Y versus the input parameters X i , i.e. the partial derivatives will be shown for two following cases:

4.35 4.3 4.25 4.2

(1) if we fix the other parameters, the partial derivatives is written by [38]

4.15 4.1 0.1

0.2

0.3

0.4

0.5

0.6

0.7

Spd Xi

Velocity (Angs/ps)

N N 1X @Ys 1 X ¼ ¼ N s¼1 @X i N s¼1

Fig. 7. The ISS at 100 K for the pulling velocity of 0.1, 0.5 and 0.75 (Å/ps).

i ¼ 1; . . . ; k:

x 10−4

Spda Xi ¼

4.6

ISS (eV/Angs3)

¼ 4.4

N 1 X

Nk s1 ¼1 N 1 X

N

k s1 ¼1

...

N X @Ys

i

sk ¼1

...

@X i

N X sk ¼1

ðsi Þ

lim

½YðX 1 ; . .. ; X i

DX i !0

! ðs Þ þ DX i ;. . . ; X k Þ  YðX 1 ; . . .; X i i ;. . . ; X k Þ ; DX i

i

¼ 1; . . .; k; si ¼ 1; . . . ; N:

ð6Þ

4.2 4 3.8 3.6 1

1.5

2

2.5

3 x 10

Number of monomers (units)

4

Fig. 8. The ISS at 100 K for the simulation boxes with 104 ; 1:5  104 and 2  104 monomeric units.

4.8

x 10−4

4.6

ISS (eV/Angs3)

ð5Þ

(2) if we vary other parameters, the partial derivatives is written as follows

SWCNT/PE system, 100 K, 0.1 Angs/ps 4.8

! ðsÞ ðsÞ YðX i þ DX i Þ  YðX i Þ ; lim DX i !0 DX i

N is the number of samples. Saltelli [38] suggested computing these partial derivatives at a set of different points in the design space of input parameters so that an average response of Y s ; s ¼ 1; . . . ; N can be obtained when moving a parameter X i of a step DX i at different points. Using this concept, the elementary effects method was developed where now the choice of the parameter increment follows another strategy. 4.5. Elementary effects Assume that the mechanical model output Y is a scalar function of the vector X as mentioned in Section 4.4. The region denoted by X which is a regular k-dimensional p-level grid unit hypercube, where each X i ; i ¼ 1; . . . ; kmay be taken on values from 1=ðp  1Þ; . . . ; 1  1=ðp  1Þ. For a given value of X, the elementary effect of the ith input is defined as

4.4

EEi ðXÞ ¼ 4.2

ð7Þ

where X 2 X, except that X i  1  D; D is a preselected multiple of 1/(p-1) and ei vector of zeros but with a unit as its ith component. The performance of an efficient sampling scheme has been suggested by Morris [40] that constructs r trajectories of ðk þ 1Þ points in the input space and leads to n ¼ rðk þ 1Þ model executions. The reader is referred to [40,41] for more details.

4 3.8 3.6 100

½YðX þ ei DÞ  YðXÞ D

200

300

400

500

600

700

800

900

1000

Chain length (units) Fig. 9. The ISS at 200 K for the chain length of 100-units, 500-units and 1000-units.

4.3. Scatter plots A preliminary step in SA is to visualize the relationship of the input parameters and the output using scatter plots. We use the model as described in Section 2 and a set of data generated by the LHS method. The advantage of this approach is that we can get a first graphical understanding under acceptable efforts the ef-

4.5.1. The computation of the sensitivity measures Assume xðgÞ and xðgþ1Þ , with g in the set 1,. . .,k. The elementary effect for the ith input parameter can be estimated along the jth trajectory via two above sampling points as presented in [38]

EEji ðxðgÞ Þ ¼

½Yðxðgþ1Þ Þ  YðxðgÞ Þ D

ð8Þ

if the ith component of xðgÞ is increased by D, and

EEji ðxðgÞ Þ ¼

½YðxðgÞ Þ  Yðxðgþ1Þ Þ D

ð9Þ

85

N. Vu-Bac et al. / Composites: Part B 59 (2014) 80–95

Fig. 10. Stress–strain response at (a) different temperatures (b) different strain rates of pure PE system.

with respect to (w.r.t) X i . V Xi ½EX i ðYjX i Þ measures the first effect of X i on the model output. 4.7. Total effect sensitivity indices Since first order sensitivity indices measure the part of variance of model output due to model input X i , an extension for higher order coupling indices is necessary. The total effect index STi is used to estimate how much the input parameters X i contribute to the output, i.e the total value of first-order terms and all high-order terms. The total effect index is defined as the [42]

Fig. 11. The schematic of PE region of interaction with an a carbon atom on the SWCNT [23].

if the ith component of xðgÞ is decreased by D.

STi ¼ 1  4.6. First-order sensitivity indices

V X i ½EX i ðYjX i Þ ; VðYÞ

ð11Þ

where V Xi ½EX i ðYjXi Þ is the variance of the expected value of Y when conditioning w.r.t all parameters except for X i , denoted as X i . V Xi ½EX i ðYjXi Þ measures the first order effects of Xi on the model output which does not contain any effect corresponding to X i . The total effect sensitivity indices measure the residual influence of X i . Due to parameter interactions, the measured sensitivity by STi P of a parameter increases, therefore STi  1 will always hold. The difference STi  Si is a measure of how much X i interacts with other input parameters. In practice, the computation of Si and STi in Eqs. (10) and (11) require a high number of samples. The costs of applying these measures is of OðNðk þ 2ÞÞ, where N is the number of samples. Thus, a

An alternative to obtain global sensitivity measures is based on variances. The response Y of a model depends on a function, Y ¼ f ðX 1 ; X 2 ; . . . ; X k Þ. The first-order sensitivity indices are defined as [42]

Si ¼

V Xi ½EX i ðYjXi Þ ; VðYÞ

ð10Þ

where V(Y) is the unconditional variance of Y and V X i ½EX i ðYjX i Þ is named the variance of the expected value of Y when conditioning

Table 2 Model uncertainties. Input parameters

SWCNT radius

Type of distribution

Temperature

Pulling velocity

25

40

25

histogram PDF normal

35

histogram PDF uniform

histogram PDF uniform

20

20

25 20 15

Posibility density

Posibility density

Posibility density

30 15

10

15

10

10 5

5

5 0 2

4

6

8

10

12

0 100

150

SWCNT radius (Angs)

Mean value Standard deviation

7.1 2

Source

[33]

200

250

300

Temperature (K)

250 pffiffiffiffiffiffi 300= 12 Assumed

350

400

0

0

0.2

0.4

0.6

Velocity (Angs/ps)

0.55 pffiffiffiffiffiffi 0.9/ 12 Assumed

0.8

1

86

N. Vu-Bac et al. / Composites: Part B 59 (2014) 80–95

Fig. 12. Schematic diagram of all sensitivity assessment methods presented in this paper.

−4

−4

x 10

x 10 5.5

ISS (eV/Angs3)

ISS (eV/Angs3)

5.5 5 4.5 4 3.5

4

6

8

10

SWCNT radius (Angs)

5 4.5 4 3.5 100

200

300

400

Temperature (K)

−4

x 10

ISS (eV/Angs3)

5.5 5 4.5 4 3.5

0.2

0.4

0.6

0.8

Velocity (Angs/ps) Fig. 13. Scatter plots of the ISS (Y) versus the SWCNT radius (X 1 ), temperature (X 2 ) and pulling velocity (X 3 ), respectively.

87

N. Vu-Bac et al. / Composites: Part B 59 (2014) 80–95

b is used to estimate the real response that will be surrogate model Y presented in next section. 4.8. Surrogate model A surrogate model [43] is an approximate mathematical model used to evaluate the mechanical model. The mechanical model needs to be evaluated multiple times considering the specific range of input parameters in order to build the surrogate model. In the scope of this paper, the popular polynomial regression method, moving least squares (MLS) approximation and a hybrid model of quadratic polynomial regression and MLS are used for building the regression model. b has been chosen to approxBased on linear polynomial basis, Y imate the response of the mechanical model Y of N observed data ðsÞ ðsÞ points X ðsÞ ¼ ½X 1 ; . . . ; X k ; s ¼ 1; . . . ; N.

b ¼b þ Y 0

k X

bi X i þ e:

ð12Þ

i¼1

where k is the number of polynomial base functions, b are the unknown regression coefficients, e is the error term that is the differb and the mechanical model Y. ence between the approximation Y The mean square difference (residual) S is defined as

S¼

N  2 X bs : Ys  Y

ð13Þ

s¼1

b is determined by minimizing the mean The parameter vector b b and square difference S between the value of surrogate surface Y the observed response Y as presented in [44]

@S ¼ 0; @bj

j ¼ 1; . . . ; k

ð14Þ

An alternative to R2 is the adjusted COD R2adj . This measure additionally takes into account the number of supporting points N as well as the number of regression coefficients kR . If N is large compared to kR the R2 and R2adj provide in the limit case the same results.

R2adj ¼ 1 

 S=ðN  kR Þ N1  1  R2 : ¼1 Stot =ðN  1Þ N  kR

ð22Þ

Besides the regression based on linear polynomials, higher order approaches considering quadratic terms can be applied

b ¼ b þ b X 1 þ b X 2 þ þ b X kN þ b X 2 þ b X 2 þ Y 0 1 2 kN 11 1 22 2 þ bkk X 2k þ e:

ð23Þ

as well as quadratic and mixed terms

b ¼ b þ b X1 þ b X2 þ þ b Xk þ b X2 þ b X2 þ Y 0 1 2 k 11 1 22 2 þ bkk X 2k þ b12 X 1 X 2 þ þ bk1k X k1 X k þ e

ð24Þ

can be used. 4.9. Moving least squares (MLS) regression model The idea of MLS is to account for a set of unorganized point samples by introducing local weighting functions wðdÞ [45]. At a b MLS is written as follows point xm , the interpolated value Y

bm ¼ Y b MLS ðxm Þ ¼ pT ðxm Þb : Y Wm

ð25Þ

where pðxm Þ is base vector evaluated at the point xm

 pT ðxm Þ ¼ 1 x1m

x2m

x21m

x22m

x1m x2m



ð26Þ

which, together with Eq. (13) results in

( " #) N k X X ðsÞ ¼ 0; Xj Y s  bi X i s¼1

j ¼ 1; . . . ; k:

ð15Þ

i¼1

The parameter vector b b is obtained by solving the system of linear equations

Qb ¼ q

ð16Þ

in which the matrix Q and vector q are

Q ij ¼

N X ðsÞ ðsÞ Xi Xj ;

qj ¼

s¼1

N X ðsÞ Y sXj ;

i; j ¼ 1 . . . k

ð17Þ

s¼1

The coefficient of determination (COD) R2 allows an estimation of the quality of the approximation,

R2 ¼ 1 

S ; Stot

0  R2  1

Stot

 ¼ VðYÞ ¼ Y  Y

T 

 YY



N X ¼1 ðY s Þ: Y N s¼1

ð27Þ

w.r.t the unknown coefficients bW yielding

 1 @LW b W ¼ XT WðxÞX ¼0!b XT WðxÞY @bW

ð28Þ

where bW are the moving coefficients, and W(x) is the diagonal matrix

2

wðx  x1 Þ

6 6 WðxÞ ¼ 6 4

( ð19Þ

ð20Þ

The COD can be rewritten as:

 PN  b s¼1 Y s  Y s R ¼ 1  PN   : s¼1 Y s  Y

T

LW ¼ ðY  pT ðxÞbW Þ WðxÞðY  pT ðxÞbW Þ

0

0



0

wðx  x2 Þ

0









0

0

wðx  xN Þ

3 7 7 7 5

ð29Þ

We use the cubic polynomial weighting function [37]

 is being the mean value of all Y s ; Y

2

and bWm contains the coefficient of the polynomial. The coefficients are determined by a weighted least squares method minimizing the L2 -norm error LW

ð18Þ

where the total variation Stot is defined as





ð21Þ

A COD of R2 = 1 means that the regression model can reproduce the response of the mechanical model exactly. To construct surrogate model connecting an alternate input–output relation R2 should be at least larger equal 0.8.

wðdÞ ¼

2

3

1  3d þ 2d

d1

0

d>1

ð30Þ

where d ¼ kx  xm k=D is the normalized distance between the interpolation point and the supporting point; D is an influence radius which is constant or dependent on position x; in our examples, D is fixed. Substituting Eq. (28) into Eq. (25) yields the MLS approximation

h i1 b MLS ¼ pT ðxÞ XT WðxÞX XT WðxÞY: Y

ð31Þ

4.10. Hybrid algorithm A hybrid (quadratic polynomial regression/MLS regression) model is presented in this section with the purpose of capturing

88

N. Vu-Bac et al. / Composites: Part B 59 (2014) 80–95 40

0.99 0.96 0.90 0.75

histogram PDF Normal PDF Weibull PDF Lognormal

35

0.50

Probability density

30

Probability

25 20 15

0.10 0.05 0.02 0.01

10

0.003

5 0

−3.4

3

4

5

Data (ev/Angs )

x 10

Fig. 16. Weibull probability plot for the distribution of the ISS.

0.995 0.99

0.997 0.99 0.98

0.95 0.9

Probability

0.95 0.90 0.75 0.50 0.25

0.75 0.5 0.25 0.1 0.05

0.10 0.05

0.01 0.005

0.02 0.01

10

0.003 3.6

3.8

4

4.2

4.4

4.6

Data (ev/Angs3)

4.8

5

5.2

Table 3 Statistical results for the ISS. Output

ISS

Mode value

Mean value

(eV=Angs3 )

(eV=Angs3 )

Standard deviation

4:02  104

4:30  104

3:77  105

ð32Þ

4.11. Normalization of the input Since the input parameters in our simulation are at different scales, they need to be normalized before building the surrogate models. The sth realization of the ith input factor X i is normalized as follows:

Table 4 Uncertainties of mechanical output using various distributions. Type of assumed PDF

Parameter 1

Parameter 2

Error

Normal PDF

Mean value

Standard deviation

Sr

3:78  105 Shape parameter (B) 11:45

2.65

(eV=Angs3 ) Weibull PDF

4:30  104 Scale parameter (A)

Log–normal PDF

4:40  104 Mean value

Standard deviation

Sr

(eV=Angs3 ) 7:76

0:09

2.67

Sr 2.37

4.12. Sensitivity analysis based on coefficients of determination

ðsÞ

X i  minðX i Þ minðX i Þ  maxðX i Þ

−3.3

Fig. 17. Log–normal probability plot for the distribution of the ISS.

It is worth mentioning that the higher order of polynomial regression model results in a more accurate approximation of the input output relationship. However, the higher order the polynomial regression model, the more flexible is the regression model. It can lead to overfitting the output of the mechanical model, see [46].

¼

10

−4

x 10

the localities of global output without overfitting. This approximation model used the MLS in the dimension where the polynomial regression fails to account for. It can be written as:

b hybrid ¼ Y bþY b MLS ðRÞ Y h i b 1; . . . ; Y s  Y b s; . . . ; YN  Y bN R ¼ Y1  Y

−3.4

Data (ev/Angs3)

Fig. 15. Normal probability plot for the distribution of the ISS.

ðsÞ;norm

10 3

Fig. 14. Histogram of the mechanical output – ISS and assumed probability density functions (PDFs) for different distribution types.

Xi

−3.3

10

6 −4

ISS (ev/Angs3)

Probability

0.25

ð33Þ

where minðX i Þ and maxðX i Þ are the minimum and maximum of the samples of the ith input.

The COD R2 for the full model is compared to the value R2pi , which is evaluated when we remove the parameter pi from the full regression basis. This removal will result in a measurable drop

N. Vu-Bac et al. / Composites: Part B 59 (2014) 80–95

DR2pi ¼ R2  R2pi   ¼ R2 ½1X 1 X 2 . . . pi1 pi piþ1 . . . X 21 X 22 . . .    R2pi ½1X 1 X 2 . . . pi1 piþ1 . . . X 21 X 22 . . .

4.13. Sensitivity assessment procedure

ð34Þ

The values of 4R2pi show which parameter pi is important. The larger the value of 4R2pi , the more important the variable is and vice versa. We normalize the value of 4R2pi with the purpose of limiting the maximum value by 1. This defines a coefficient of importance

I pi ¼

R2pi DR p i ¼ 1  : R2 R2

ð35Þ

1

0.9

R

2

0.8

0.7 Linear Quadratic without mixed terms Full quadratic

0.6

MLS Hybrid 0.5

20

40

60

80

100

120

140

89

160

180

The schematic diagram of all sensitivity assessment methods are provided in Fig. 12. 5. Numerical results Scatter plots in Fig. 13 show the influence of the SWCNT radius, temperature and pulling velocity on the ISS, respectively. The highest influence seems to be the SWCNT radius. To quantify the affect of the input parameters, the sensitivity indices were studied using the surrogate model. The output (ISS) uncertainty of the investigated mechanical model can be described by a probability distribution function. Judging from the shape of the histogram shown in Fig. 14, either normal or Weibull, log–normal distributions of which the parameters are listed in Table 4 are proper candidates that can be used to account for the mechanical output uncertainty. The purpose of probability plot is to graphically assess if the ISS’s data can be characterized by a given distribution. If the probability plot produces the least deviations from a linear line, the exploited theoretical distribution is chosen as a goodness-of-fit of the data [47]. The normal, Weibull and log–normal probability plots are illustrated in Figs. 15–17, respectively. The difference between mean (4:30  104 ) and mode value (4:02  104 ) reveals that the output is skewed, see Fig. 14 and Table 3. Hence, the normal and log–normal distributions are not

200

Number of samples Fig. 18. The plot of R2 versus the number of samples when polynomial regression, MLS and hybrid models are used

Table 5 Regression coefficient, COD R2 and adjusted COD R2adj . Response surface methods

Regression coefficient

Linear regression R2 = 0.81

b0 = 5.27e4 b1 = 1.21e5 b2 = 1.80e7

R2adj = 0.81

b3 = 4.76e5

Quadratic without mixed terms

R2 = 0.85

b0 = 5.66e4 b1 = 3.34e5 b2 = 3.19e8 b3 = 8.59e5 b11 = 1.64e6 b22 = 3.14e10

R2adj = 0.85

b33 = 3.24e5

Full quadratic

R2 = 0.86

b0 = 5.76e4 b1 = 3.39e5 b2 = 1.07e7 b3 = 8.93e5 b11 = 1.48e6 b22 = 3.81e10 b33 = 3.47e5 b12 = 1.35e8 b13 = 1.25e6

R2adj = 0.86

b23 = 2.90e8

Moving least squares R2 = 0.88 R2adj = 0.88 Hybrid model of polynomial regression and MLS R2 = 0.90 R2adj

= 0.90

Fig. 19. Scatter points got from MD simulations (black points) and projection surface of surrogate model which express the ISS versus the SWCNT radius and temperature (a) and the ISS versus the SWCNT radius and pulling velocity (b). In this figure, the full quadratic regression is used.

90

N. Vu-Bac et al. / Composites: Part B 59 (2014) 80–95

−4

ISS (ev/Angs3)

x 10

5 4.5 4 3.5

300

Temperature (K)

200 100

2

4

6

8

10

SWCNT radius (Angs)

(a)

ISS (ev/Angs3)

x 10

−4

5

4.5

4

1 0.8

10

0.6

Velocity (Angs/ps)

8

0.4

6 4

0.2

SWCNT radius (Angs)

(b) Fig. 20. Scatter points got from MD simulations (black points) and projection surface of surrogate model which express the ISS versus the SWCNT radius and temperature (a) and the ISS versus the SWCNT radius and pulling velocity (b). In this figure, the MLS is used.

a good candidate to characterize the mechanical output. Note that in this example the log–normal distribution is nearly symmetric. It is more reasonable to describe the skewed data by a Weibull distribution. In order to quantitatively assess the best fits to the data, the sum of the squares of the estimated residuals (Sr ) are computed in Eq. (13) in which cross blue points of the probability plots are considered as Y s and the regression values (linear red lines) play the b s . Obviously, the Weibull probability plot provides same role as Y the best fit to characterize the ISS’s histogram due to the smallest Sr given in Table 4. Note that the Sr values are estimated within 95% confidence intervals. Fig. 18 shows that statistical convergence is achieved with 200 samples for all surrogate models the polynomial regression, MLS and hybrid models with appropriate chosen radius (normalized radius of 0.6). In the following, the surrogate models are built to approximate the output of the mechanical model (200 samples). The COD R2 and adjusted COD R2adj in accordance with the linear, quadratic (with and without mixed terms), MLS and hybrid regression models are calculated and it is demonstrated that the predicted response using the surrogate models are good approximations of the response of the mechanical model, see Table 5. As listed in Table 5, the COD is 0.90 for hybrid method, 0.88 for MLS with appropriate chosen radius (normalized radius of 0.6) and 0.86 (0.85) for quadratic (full and without mixed terms) regression,

Fig. 21. Scatter points got from MD simulations (black points) and projection surface of surrogate model which express the ISS versus the SWCNT radius and temperature (a) and the ISS versus the SWCNT radius and pulling velocity (b). In this figure, the hybrid model is used

indicating that the surrogate model reflects almost exactly the response of the mechanical model. In other words, 90%, 88%, 86 (85)% of the ISS variation are represented by the hybrid the MLS and quadratic (with and without mixed terms) regression models, respectively. The obtained R2 ; R2adj by using linear regression is 0.82 which is slightly lower than by the other methods. We have used the full quadratic regression model in this paper for further comparison. Three-dimensional scatter plots and the associated surrogate models are shown in dependence of the SWCNT radius and temperature, pulling velocity in Figs. 19 (full quadratic), Fig. 20 (MLS) and Fig. 21 (hybrid). The steepest gradient in the SWCNT radius direction indicates the SWCNT radius is a key parameter for the ISS. 5.1. Derivative-based approach Since long MD simulations are generating cumulative errors in the numerical integration that can be minimized with proper selection of algorithms and parameters, but not eliminated entirely. Therefore, the approximation of partial derivatives with finite different schemes is critical in particular when 4X i approaches to 0. This is the reason, why we apply the partial derivative approach solely on the generated response surfaces, which, by construction, are sufficiently smooth for approximating derivatives numerically. The full quadratic surrogate model is used for this method. The interval [0, 1] is split into 9 equal sized

91

N. Vu-Bac et al. / Composites: Part B 59 (2014) 80–95

(a)

(b)

Fig. 22. Flowchart of partial derivatives of the surrogate model output w.r.t X i when we fix the others (a) and when we vary the others (b); N is the number of samples.

sub-intervals and placing X i at each points. By applying the inverse cumulative function we obtain a 10-level grid for three input parameters (SWCNT radius (X 1 ), temperature (X 2 ) and pulling velocity (X 3 )) in the interval. We evaluate Y twice: first at X i and second at X i þ DX i with DX i ¼ 104 . The partial derivatives are computed by Eqs. (5) and (6) for two cases: (1) keeping the other parameters fixed; (2) considering the variation of other parameters, respectively. Flowcharts of these processes are provided in Fig. 22. Then, the partial derivatives are normalized w.r.t an input factor, e.g., X 1 in this example named SXX 1i ; i ¼ 1; 2; 3 as follows X ;ðpdÞ

SX 11

¼

X ;ðpdaÞ

SX 11

X ;ðpdÞ SX 12

¼

X ;ðpdaÞ

SX 12

X ;ðpdÞ

SX 13

¼

¼

X ;ðpdaÞ

SX 13

¼

¼

N 1X @Ys N s¼1 @X 1

!

N X N X N 1 X @Ys 3 N s¼1 j¼1 l¼1 @X 1

! N 1X @Ys X 2 N s¼1 @X 2 X 1 1

N X N X N X

N3

j¼1 s¼1 l¼1

! N 1X @Ys X 3 N s¼1 @X 3 X 1 1

N X N X N X

N3

l¼1 s¼1 j¼1

Normalized value of ARNPD calculated by Eqs. (39b), (40b), (41b)

bX1 ;ðpdÞ c X 1 ;ðpdÞ ¼ PS Xi NS Xi 3 bS XX1 ;ðpdÞ

bX1 ;ðpdaÞ c X 1 ;ðpdaÞ ¼ PS Xi NS Xi 3 bS XX1 ;ðpdaÞ

c X 1 ;ðpdÞ NS X1

c X 1 ;ðpdaÞ NS X1

i¼1

i¼1

i

= 0.50

i

= 0.51

c X 1 ;ðpdÞ = 0.35 NS X2

c X 1 ;ðpdaÞ = 0.34 NS X2

c X 1 ;ðpdÞ = 0.15 NS X3

c X 1 ;ðpdaÞ = 0.15 NS X3

These normalized partial derivatives (NPD) b S XX 1i ; i ¼ 1; 2; 3 are reduced by the COD R2

ð36bÞ

X ;ðpdÞ X ;ðpdÞ b ¼ R2 SX11 S X11

ð39aÞ

X ;ðpdaÞ X ;ðpdaÞ b ¼ R2 SX11 S X11

ð39bÞ

X ;ðpdÞ X ;ðpdÞ b ¼ R2 SX12 S X12

ð40aÞ

X ;ðpdaÞ b S X12

ð40bÞ

ð37aÞ !

ð37bÞ

ð38aÞ !

@Yl X 3 @X 3 X 1

Normalized value of ARNPD calculated by Eqs. (39a), (40a), (41a)

ð36aÞ !

@Yj X 2 @X 2 X 1

Table 6 Normalized average of reduced NPD for the ISS.

ð38bÞ

¼

X ;ðpdaÞ R2 SX12

X ;ðpdÞ X ;ðpdÞ b ¼ R2 SX13 S X13

ð41aÞ

X ;ðpdaÞ b S X13

ð41bÞ

¼

X ;ðpdaÞ R2 SX13

Table 6 shows the average of the reduced normalized partial derivatives (ARNPD) of the ISS versus the input parameters. Since X ;ðpdÞ the partial derivatives b are the same as the partial derivatives S 1 Xi

92

N. Vu-Bac et al. / Composites: Part B 59 (2014) 80–95

Fig. 23. Reduced pull-out energy of three trajectories (a) trajectory 1, (b) trajectory 2, and (c) trajectory 3.

X ;ðpdaÞ b ; i ¼ 1; 2; 3, there is no interaction between the input paramS X 1i

c 1 eters. The key input factor is the SWCNT radius ( NS X1

X ;ðpdÞ

followed by the temperature

c X 1 ;ðpdÞ ( NS X2

has smaller influence on the ISS

¼ 0:50)

¼ 0:35). The pulling velocity

c X 1 ;ðpdÞ ( NS X3

¼ 0:15).

Fig. 24. Reduced pull-out energy of three trajectories (a) trajectory 4, (b) trajectory 5, and (c) trajectory 6.

Table 7 shows the elementary effects relative to parameter i; i ¼ 1; 2; 3 computed along trajectory r; r ¼ 1; . . . ; 9, the average EEi and the normalized value of the average NEEi . The SWCNT radius is observed as the key factor (NEE1 ¼ 0:476) on the ISS followed by the temperature (NEE2 ¼ 0:329). The influence of the pulling velocity on the ISS is small (NEE3 ¼ 0:195). 5.3. Sensitivity analysis based on coefficients of determination

5.2. Elementary effects The elementary effects are calculated based on the mechanical model for 9 trajectories and 4 samples in each trajectory. So it is totally 36 samples tested. For the sake of comparison, the pullout energy is reduced by prL2 as shown from Fig. 23–25. Each of these figures corresponds to 4 points of a trajectory. Note that the legend rad (3.52) – temp (235) – vel (0.28) e.g., in the Fig. 23 expresses the SWCNT radius (X 1 ) = 3.52 Å, temperature (X 2 ) = 235 K and pulling velocity (X 3 ) = 0.28 Å/ps.

In order to evaluate the importance of variables, we use the SA based coefficients of determination approach by leaving out terms containing X 1 ; X 2 and X 3 from the full quadratic regression basis and compute the corresponding CODs, drop values 4R2i ; i ¼ 1; 2; 3 of reduced models, Eqs. (34) and (35). For example, carrying out a regression on the reduced full quadratic model containing X 2 ; X 3 , the reduced quadratic model is given by

b ¼ b þ b X2 þ b X3 þ b X2 þ b X2 þ b X2X3; Y 0 2 3 22 2 33 3 23

ð42Þ

93

N. Vu-Bac et al. / Composites: Part B 59 (2014) 80–95

The above model results in a COD (R21 ) and an adjusted COD The reduced MLS and hybrid regression models are constructed in the same fashion. The corresponding CODs, adjusted CODs, drop values 4R2i and normalized values Ii of 4R2i are listed in Table 8. As shown in Table 8, the big drop in the value of R2 and R2adj when removing terms X 1 ; X 21 ; X 1 X 2 ; X 1 X 3 from the regression shows the relative importance of the SWCNT radius (X 1 ) compared to the others. The insignificant decrease in R2 and R2adj when removing terms associated either the variable X 2 or the variable X 3 indicates that the temperature has less influence on the ISS compared to the SWCNT radius. The influence of the pulling velocity on the ISS is smallest. (R2adj1 ).

Table 7 Elementary effects for the ISS. Trajectory n Elementary effects

EE1;ðrÞ  103

EE2;ðrÞ  103

EE3;ðrÞ  103

Trajectory Trajectory Trajectory Trajectory Trajectory Trajectory Trajectory Trajectory Trajectory

0.08 0.11 0.12 0.11 0.13 0.16 0.13 0.13 0.05

0.05 0.07 0.08 0.09 0.09 0.11 0.08 0.10 0.04

0.03 0.05 0.04 0.03 0.07 0.08 0.04 0.06 0.03

Average elementary effects

EE1 ¼ 0:110

EE2 ¼ 0:076

EE3 ¼ 0:045

i Normalized NEEi ¼ PEE 3

NEE1 ¼ 0:476

NEE2 ¼ 0:329

NEE3 ¼ 0:195

1 2 3 4 5 6 7 8 9

i¼1

5.4. Global sensitivity analysis Based on the determined surrogate models, we generated 50,000 Latin hypercube samples. The regression based sensitivity

EEi

Note: EE1;ðrÞ ; EE2;ðrÞ ; EE3;ðrÞ indicate the elementary effects of the ISS relative to the SWCNT radius (X 1 ), temperature (X 2 ) and pulling velocity (X 3 ) computed along trajectory r, respectively.

indices by using those samples as approximation points are given in Table 9. The sensitivity indices Si and STi ; i ¼ 1; 2; 3, Eqs. (10) and (11), are calculated in accordance with the polynomial regression, the MLS and the hybrid approximation. These values are reduced by the COD inferring that only R2 % of response can be explained with surrogate model. The reduced sensitivity indices b S Ti are given by S i and b

b S i ¼ R 2 Si b S Ti ¼ R2 STi

ð43aÞ ð43bÞ

The results of the variance-based methods are summarized in Table 9 and illustrated in Fig. 26. They show: 1. The reduced first-order indices b S i and reduced total-effect indices b S Ti are nearly identical indicating the independence of the input-parameters. These results agree very well with results of the derivative-based approach. 2. All surrogate models are suitable to study SA. The highest COD is obtained by the hybrid model (0.90) followed by the MLSapproach (0.88), the quadratic (with and without mixed terms) regression model (0.86 and 0.85) and the linear regression model (0.82); Table 5. 3. The SWCNT radius is the most influential parameter followed by the temperature and pulling velocity, respectively. The ratio b S T2 =b S T1 is with a value of 0.64 slightly lower compared to the value of the derivative-based approach and elementary effects (0.69). While the ratio b S T3 =b S T1 ¼ 0:4 is identical to the ratio of the elementary effects, it is higher than the value of the derivative-based method (0.30), see Tables 9 + 7 + 6. 6. Conclusions and future work

Fig. 25. Reduced pull-out energy of three trajectories (a) trajectory 7, (b) trajectory 8, and (c) trajectory 9.

The mechanical properties of SWCNT/PE composites have been obtained by using MD simulations. The mechanical output (ISS) is best characterized by the Weibull distribution. Based on the mechanical and surrogate models, the influence of the SWCNT radius, temperature and pulling velocity on the ISS in the context of both averaged local and global SA is quantitatively studied. For the examples studied, the elementary-effects method requires less samples than the other methods. However, it fails to estimate the mutual interaction of the input parameters. The variance-based method is a global SA method that can be used for such purpose. The regression methods are robust methods because they do not require a large number of samples number of input parameters.

94

N. Vu-Bac et al. / Composites: Part B 59 (2014) 80–95

Table 8 The COD, adjusted COD, drop value of COD, and the coefficient of importance of the reduced quadratic polynomial, the reduced MLS and the reduced hybrid models. The reduced model containing X 2 ; X 3

The reduced model containing X 1 ; X 3

The reduced model containing X 1 ; X 2

R21 = 0.39

R22 = 0.70

R23 = 0.74

R2adj1 = 0.39 DR21 = 0.47

R2adj2 = 0.70 DR22 = 0.16

R2adj3 = 0.74

I1 = 0.54

I2 = 0.18

DR23 = 0.12 I3 = 0.13

R21 = 0.41

R22 = 0.70

R23 = 0.76

R2adj1 = 0.39 DR21 = 0.46

R2adj2 = 0.69 DR22 = 0.16

R2adj3 = 0.75

I1 = 0.53

I2 = 0.19

DR23 = 0.10 I3 = 0.12

R21 = 0.41

COD R22 = 0.70

R23 = 0.76

R2adj1 = 0.39 DR21 = 0.46

R2adj2 = 0.69 DR22 = 0.16

R2adj3 = 0.75

I1 = 0.53

I2 = 0.19

The reduced quadratic polynomial model

The reduced MLS model

The reduced hybrid model

DR23 = 0.10 I3 = 0.12

Note: R21 ; R2adj1 ; DR21 and I1 are the COD, the adjusted COD, the drop value of COD and the coefficient of importance when the SWCNT radius (X 1 ) is removed from the regression model, respectively. R22 ; R2adj2 ; DR22 and I2 are the COD, the adjusted COD, the drop value of COD and the coefficient of importance when the temperature (X 2 ) is removed from the regression model, respectively. R23 ; R2adj3 ; DR23 and I3 are the COD, the adjusted COD, the drop value of COD and the coefficient of importance when the pulling velocity (X 3 ) is removed from the regression model, respectively.

Table 9 First-order and total effects sensitivity indices computed on the surrogate model of the input parameters contributing to the ISS Linear regression model

Quadratic without mixed terms

Full quadratic

MLS

Hybrid

First-order indices b Si b S 1 = 0.44

First-order indices b Si b S 1 = 0.41

First-order indices b Si b S 1 = 0.41

First-order indices b Si b S 1 = 0.41

First-order indices b Si b S 1 = 0.40

b S 2 = 0.23 b S 3 = 0.14 P3 b i¼1 S i = 0.81

b S 2 = 0.26 b S 3 = 0.17 P3 b i¼1 S i = 0.84

b S 2 = 0.27 b S 3 = 0.17 P3 b i¼1 S i = 0.85

b S 2 = 0.30 b S 3 = 0.15 P3 b i¼1 S i = 0.86

b S 2 = 0.30 b S 3 = 0.16 P3 b i¼1 S i = 0.86

Total-effect indices b S Ti b S T1 = 0.44

Total-effect indices b S Ti b S T1 = 0.41

Total-effect indices b S Ti b S T1 = 0.42

Total-effect indices b S Ti b S T1 = 0.41

Total-effect indices b S Ti b S T1 = 0.44

b S T2 = 0.23 b S T3 = 0.15 P3 b i¼1 S Ti = 0.82

b S T2 = 0.27 b S T3 = 0.17 P3 b i¼1 S Ti = 0.85

b S T2 = 0.27 b S T3 = 0.17 P3 b i¼1 S Ti = 0.87

b S T2 = 0.29 b S T3 = 0.17 P3 b i¼1 S Ti = 0.88

b S T2 = 0.31 b S T3 = 0.17 P3 b i¼1 S Ti = 0.92

Note: b S 1 and b S T1 are the reduced first-order and total-effect indices of the ISS (Y) due to the SWCNT radius (X 1 ), respectively, b S 2 and b S T2 are the reduced first-order and total-effect indices of the ISS (Y) due to the temperature (X 2 ), respectively, b S 3 and b S T3 are the reduced first-order and total-effect indices of the ISS (Y) due to the pulling velocity (X 3 ), respectively.

First−order indices

Total−effect indices

0.45

0.45 Linear Quadratic w/o mixed terms Full quadratic MLS Hybrid

0.4 0.35

0.35

0.3

0.3

0.25

0.25

0.2

0.2

0.15

0.15

0.1

0.1

0.05

0.05

0

SWCNT radius

Temperature

Velocity

Linear Quadratic w/o mixed terms Full quadratic MLS Hybrid

0.4

0

SWCNT radius

Temperature

(a) Fig. 26. Sensitivity indices (a) first-order and (b) total effect.

(b)

Velocity

N. Vu-Bac et al. / Composites: Part B 59 (2014) 80–95

The hybrid and MLS surrogate models are good approximations of the mechanical model. However, the adjustment of the influence radius D used in hybrid and MLS surrogate models remains an open issue. All methods predict the same tendencies, i.e.: 1. The most influential parameter on the ISS is the SWCNT radius, followed by the temperature. 2. The pulling velocity has relatively small influence on the ISS. In the future, we will develop the work to stochastic multiscale modeling of CNT reinforced polymers and test a wider variety of input parameters on the predicted output. In particular, we intend to study the influence of the length, orientation, curvature and dispersion of CNTs on mechanical properties as presented by [48]. Graph theory could be employed to quantitatively assess the quality of those coupled partial models on different scales.

Acknowledgement We gratefully acknowledge the support by the Deutscher Akademischer Austausch Dienst (DAAD), IRSES-MULTIFRAC and the German Research Foundation (DFG) through the Research Training Group 1462. X. Zhuang acknowledges the supports from the NSFC (41130751), National Basic Research Program of China (973 Program: 2011CB013800), Changjiang Scholars and Innovative Research Team (PCSIRT, IRT1029). and Shanghai Pujiang Program (12PJ1409100).

References [1] Dai H. Carbon nanotubes: opportunities and challenges. Surf Sci 2002;500:218–41. [2] Salvetat-Delmotte JP, Rubio A. Mechanical properties of carbon nanotubes: a ber digest for beginners. Carbon 2002;40:1729–34. [3] Lau KT, Gu C, Hui D. A critical review on nanotube and nanotube/nanoclay related polymer composite materials. Compos Part B – Eng 2006;37:425–36. [4] Schadler LS, Giannaris SC, Ajayan PM. Load transfer in carbon nanotube epoxy composites. Appl Phys Lett 1998;73(26):3842–4. [5] Ajayan PM, Schadler LS, Giannaris SC, Rubio A. Single-walled carbon nanotube polymer composites: strength and weakness. Advan Mater 2000;12:750. [6] Shaffer MSP, Windle AH. Fabrication and characterization of carbon nanotube/ poly(vinyl alcohol) composites. Advan Mater 1999;11:937. [7] Desai AV, Haque MA. Review mechanics of the interface for carbon nanotube polymer composites. Thin-Walled Struct 2011;1854–60. [8] Piggott MR. A new model for interface failure in fibre-reinforced polymers. Compos Sci Technol 1995;55:269–76. [9] Qian D, Dickey EC, Andrews R, Rantell T. Load transfer and deformation mechanisms in carbon nanotube-polystyrene composites. Appl Phys Lett 2000;20:2868–70. [10] Deng F. Investigation of the interfacial bonding and deformation mechanism of the nano composites containing carbon nanotubes. Tokyo University, PhD Dissertation; 2008. [11] Barber AH, Cohen SR, Kenig S, Wagner HD. Measurement of carbon nanotubepolymer interfacial strength. Appl Phys Lett 2003;82(23):4140–2. [12] Barber AH, Cohen SR, Kenig S, Wagner HD. Interfacial fracture energy measurements for multi-walled carbon nanotubes pulled from a polymer matrix. Compos Sci Technol 2004;64:2283–9. [13] Cooper CA, Ravich D, Lips D, Mayer J, Wagner HD. Distribution and alignment of carbon nanotubes and nanofibrils in a polymer matrix. Compos Sci Technol 2002;62:1105–12. [14] Cooper CA, Ravich D, Lips D, Mayer J, Wagner HD. Load transfer in carbon nanotube epoxy composites. Appl Phys Lett 2002;81(20):3873–5. [15] Rafiee R, Rabczuk T, Pourazizi R, Zhao J, Zhang Y. Challenges of the modeling methods for investigating the interaction between the CNT and the surrounding polymer. Advan Mater Sci Eng; 2013, Article ID 183026, doi:http://dx.doi.org/10.1155/2013/183026. [16] Liao K, Li S. Interfacial characteristics of a carbon nanotube-polystyrene composite system. Appl Phys Lett 2001;79:4225–7. [17] Frankland SJV, Harik VM. Analysis of carbon nanotube pull-out from a polymer matrix. Surf Sci 2003;525(1-3):L103–8.

95

[18] Gou J, Minaei B, Wang B, Liang Z, Zhang C. Computational and experimental study of interfacial bonding of single-walled nanotube reinforced composites. Comp Mater Sci 2004(31):225–36. [19] Chowdhury SC, Okabe T. Computer simulation of carbon nanotube pull-out from polymer by the molecular dynamics method. Composites Part A 2007;3(38):747–54. [20] Zheng Q, Xia D, Xue Q, Yan K, Gao X, Li Q. Computational analysis of effect of modification on the interfacial characteristics of a carbon nanotubepolyethylene composites system. Appl Surf Sci 2009;255:3524–43. [21] Wei C, Srivastava D, Cho K. Thermal expansion and diffusion coefficients of carbon nanotube-polymer composites. Nano Lett 2002;2(6):647–50. [22] Wei C. Adhesion and reinforcement in carbon nanotube polymer composite. Appl Phys Lett 2006;88:093–108. [23] Zhang ZQ, Ward DK, Xue Y, Zhang HW, Horstemeyer MF. Interfacial characteristics of carbon nanotube-polyethylene composites using molecular dynamics simulations. ISRN Mater Sci; 2011, doi:10.5402/2011/145042. [24] Li Y, Liu Y, Peng X, Yan C, Liu S, Hu N. Pull-out simulations on interfacial properties of carbon nanotube-reinforced polymer nanocomposites. Comput Mater Sci 2011(50):1854–60. [25] Brenner DW, Shenderova OA, Harrison JA, Stuart SJ, Ni B, Sinnott SB. A secondgeneration reactive empirical bond order (REBO) potential energy expression for hydrocarbons. J Phys: Conden Matter 2002(14):783–802. [26] Plimpton SJ. Fast parallel algorithms for short-range molecular dynamics. J Comput Phys 1995;117:1–19. [27] Zhao J, Jiang JW, Wei N, Zhang Y, Rabczuk T. Thermal conductivity dependence on chain length in amorphous polymers. J Appl Phys 2013(113). http:// dx.doi.org/10.1063/1.4804237.. [28] Hossain D, Tschopp MA, Ward DK, Bouvard JL, Wang P, Horstemeyer MF. Molecular dynamics simulations of deformation mechanisms of amorphous polyethylene. Polymer 2010;51:6071–83. [29] Shepherd JE. Multiscale modeling of the deformation of semi-crystalline polymers. Ph.D thesis at Georgia Institute of Technology; 2006. [30] Vu-Bac N, Lahmer T, Keitel H, Zhao J, Zhuang X, Rabczuk T. Stochastic predictions of bulk properties of amorphous polyethylene based on molecular dynamics simulations. Mech Mater 2014;68:70–84. http://dx.doi.org/10.1016/ j.mechmat.2013.07.021. [31] Zhang Y, Zhao J, Wei N, Jiang JW, Rabczuk T. Effects of the dispersion of polymer wrapped two neighbouring single walled carbon nanotubes (swnts) on nanoengineering load transfer. Compos Part B: Eng 2013;45:1714–21. [32] Park S, Schulten K. Calculating potentials of mean force from steered molecular dynamics simulations. J Chem Phys 2004;120:5946–61. [33] Tian Y, Jiang H, Pfaler Jv, Zhu Z, Nasibulin AG, Nikitin T, et al. Analysis of the size distribution of single-walled carbon nanotubes using optical absorption spectroscopy. J Phys Chem Lett 2010;1(7):1143–8. [34] Izrailev S, Stepaniants S, Isralewitz B, Kosztin D, Lu H, Molnar F, et al. Steered molecular dynamics. In: Computational molecular dynamics: challenges, methods, ideas, Lecture notes in computational science and engineering, vol. 4. Berlin: Springer-Verlag; 1998. p. 39–65. [35] McKay MD, Conover WJ, Beckmann RJ. A comparison of three methods for selecting values of input variables in the analysis of output from a computer code. Technometrics 1979;21:239–45. [36] Iman RL, Conover WJ. A distribution-free approach to inducing rank correlation among input variables. Commun Stat – Simul Comput 1982;11:311–34. [37] Most T, Bucher C. A moving least squares weighting function for the elementfree Galerkin method which almost fulfills essential boundary conditions. Struct Eng Mech 2005;21(3):315–32. [38] Saltelli A, Ratto M, Andres T, Campolongo F, Cariboni J, Gatelli D, et al. Global sensitivity analysis. The Primer. John Wiley and Sons; 2008. [39] Paruolo P, Saisana M, Saltelli A. Ratings and rankings: voodoo or science? The Roy Stat Soc: J Ser A 2012;176(2):1–26. [40] Morris MD. Factorial sampling plans for preliminary computational experiments. Technometrics 1991;33(2):61–174. [41] Campolongo F, Braddock R. The use of graph theory in the sensitivity analysis of the model output: a second order screening method. Reliab Eng Syst Safety 1999:1–12. [42] Keitel H, Karaki G, Lahmer T, Nikulla S, Zabel V. Evaluation of coupled partial models in structural engineering using graph theory and sensitivity analysis. Eng Struct 2011:3726–36. [43] Myers RH, Montgomery DC. Response surface methodology process and product optimization using designed experiments. 2nd ed. New York: John Wiley; 2002. [44] Bucher C. Computational analysis of randomness in structural mechanics, Structures and infrastructures book series, vol. 3. CRC Press; 2009. [45] Lancaster P, Salkauskas K. An introduction: curve and surface fitting. Academic Press; 1986. [46] Karaki Ghada. Assessment of coupled models of bridges considering timedependent vehicular loading. Ph.D thesis, Bauhaus-Universität Weimar; 2011. [47] Gan F, Koehler KJ. Goodness-of-fit tests based on P–P probability plots. Technometrics 1990(32):289–303. [48] Shokrieh MM, Rafiee R. Stochastic multi-scale modeling of CNT/polymer composites. Comput Mater Sci 2010(50):437–46.