Stochastic multiscale modeling of heat conductivity of Polymeric clay nanocomposites

Stochastic multiscale modeling of heat conductivity of Polymeric clay nanocomposites

Journal Pre-proof

Stochastic multiscale modeling of heat conductivity of Polymeric clay nanocomposites Bokai Liu, Nam Vu-Bac, Xiaoying Zhuang, Timon Rabczuk PII: DOI: Reference:

S0167-6636(19)30857-9 https://doi.org/10.1016/j.mechmat.2019.103280 MECMAT 103280

To appear in:

Mechanics of Materials

Received date: Revised date: Accepted date:

3 October 2019 9 December 2019 10 December 2019

Please cite this article as: Bokai Liu, Nam Vu-Bac, Xiaoying Zhuang, Timon Rabczuk, Stochastic multiscale modeling of heat conductivity of Polymeric clay nanocomposites, Mechanics of Materials (2019), doi: https://doi.org/10.1016/j.mechmat.2019.103280

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier Ltd.

Highlights • We propose a stochastic multiscale method to quantify the heat conductivity of polymeric nanocomposites (PNCs) with clay reinforcement • We systematically studied the effects of different parameters on the thermal conductivity of composites. • The most influential input parameters are the aspect ratio followed by the volume fraction. • The Kapitza Resistance has no significant effect on the thermal conductivity of the PNCs. • We show that the most accurate surrogate model is the moving least square (MLS).

1

Stochastic multiscale modeling of heat conductivity of Polymeric clay nanocomposites Bokai Liud , Nam Vu-Bacc , Xiaoying Zhuangc , Timon Rabczuka,b,∗ a

Division of Computational Mechanics, Ton Duc Thang University, Ho Chi Minh City, Viet Nam b Faculty of Civil Engineering, Ton Duc Thang University, Ho Chi Minh City, Viet Nam c Institute of Continuum Mechanics, Leibniz University of Hannover, 30167 Hannover, Germany d Institute of Structural Mechanics, Bauhaus-Universit¨ at Weimar, Marienstr. 15, D-99423 Weimar, Germany

Abstract We propose a stochastic multi-scale method to quantify the most significant input parameters influencing the heat conductivity of polymeric nanocomposites (PNCs) with clay reinforcement. Therefore, a surrogate based global sensitivity analysis is coupled with a hierarchical multi-scale method employing computational homogenization. The effect of the conductivity of the fibers and the matrix, the Kapitza resistance, volume fraction and aspect ratio on the ’macroscopic’ conductivity of the composite is systematically studied. We show that all selected surrogate models yield consistently the conclusions that the most influential input parameters are the aspect ratio followed by the volume fraction. The Kapitza Resistance has no significant effect on the thermal conductivity of the PNCs. The most accurate surrogate model in terms of the R2 value is the moving least square (MLS). Keywords: Multi-scale modeling, Uncertainty quantification, Polymeric nano-composites(PNCs), Heat conductivity, Stochastic modeling 1. Introduction Polymeric (epoxy) clay nanocomposites (PCNs) have been researched extensively recent years because of their excellent thermal and mechanical ∗

Corresponding Author. E-mail address: [email protected]

Preprint submitted to Elsevier

December 14, 2019

properties. Compared to other materials like metallic and ceramic composites, polymers are commonly employed as matrix since they are easy to fabricate and have reasonable strength and durability, superior resistance to chemicals and corrosion, and attractive thermal/mechanical properties[1]. The thermal conductivity is important in many applications, especially in aerospace engineering, for electronic devices as well as in civil engineering. It is also essential in high-energy density systems since electronic components need heat dissipation functionality. Or in other words, in electronic devices the generated heat should ideally be dissipated by light and small heat sinks. Since most polymer exhibit a low thermal conductivity, the polymer matrix is filled with a reinforcement to increase the macroscopic thermal conductivity. Promising candidates for fillers are nano-materials such as graphene or CNTs with extraordinary thermal properties. Another good candidate, which is much cheaper, is clay. The macroscopic thermal conductivity of composites can be influenced by many factors, for instance the type and properties of fillers, their dispersion and orientation, the property of the polymer matrix or in general the microstructures of the composite. Although the nanofillers improve the properties of the composites, predicting reliably and accurately the properties of the composite still remains a challenging task. Since experiments are time-consuming, expensive and also sometimes unfeasible, computational methods have been developed in order to support the design of new composite materials and to gain a better understanding of nanocomposites [2][3]. Computational models are useful and powerful tools in order to study the influence of numerous aspects on the effective thermal conductivity (ETC) of composites. Mean field mircomechanics methods for instance are simple and efficient approaches to predict the effective properties of different composites. The Mori-Tanaka(MT) method [4] is probably among the most popular mean field approaches. Strong contrast method [5] is a statistical continuum theory based on correlation functions to statistically describe the microstructure. Liu [6] devised an analytical model, the step-by-step method of the rule-of-mixture, to predict and analyze the properties of fiber-reinforced composites. Sheng et al. defined an ’effective particle’ and modelled a nanoclay structure [7]. Zhai et al. proposed a novel theoretical model to estimate the effective thermal conductivity of nanoparticle reinforced polymeric composites [8]. Pourasghar et al. estimated the density, thermal relaxation time, heat capacity and thermal conductivity for CNT based polymeric composites of different volume fractions [9]. Mortazavi et al. built a two-phase random 3

composite using the Finite Element Method (FEM) to compute its thermal and mechanical properties [10]. Yang et al. carried out molecular dynamics (MD) simulations to investigate the effect of carbon nanobud (CNB) on the mechanical properties and thermal conductivity of cross-linked epoxy resins based on the DGEBA and MTHPA [11]. The studies in [12] [13] [14] indicate that there is a specific interphase layer between the fillers and polymer matrix which might affect the homogenized thermal conductivity of the composite. This so-called Kapitza resistance is widely used to analyse such interfacial effects without modelling the interphase layer explicitly [15].In general, there are many effective methods to deal with interface problems, such as theoretical analysis or numerical simulation[16][17][18]. While most studies are based on deterministic approaches, there is a comparatively lower number of stochastic methods accounting for uncertainties in the input parameters. In deterministic models, the output of the model is fully determined by the parameter values and the initial conditions. However, uncertainties in the input parameters such as aspect ratio, volume fraction, thermal properties of fiber and matrix need to be taken into account for reliable predictions. In this manuscript, we provide a stochastic multi-scale method in order to study the influence of numerous uncertain input parameters on the thermal conductivity of the composite. Therefore, we present a hierarchical multi-scale method based on computational homogenization in order to predict the macroscopic thermal conductivity based on the fine-scale structure. We use the finite element method and employ surrogate models to conduct a Global Sensitivity Analysis (GSA). The SA is performed in order to quantify the influence of the conductivity of the fiber, matrix, Kapitza resistance, volume fraction and aspect ratio on the macroscopic conductivity. Therefore, we compute first-order and total-effect sensitivity indices with different surrogate models. This paper is organized as follows. The RVE generation algorithm, the details of the models are presented in Section 2. The FE model is briefly described subsequently. We present the stochastic approach in Section 3 before the sensitivity analysis based on different surrogate models is discussed in Section 4. Finally, we discuss the numerical results in Section 5. The conclusion is presented in Section 6.

4

2. Polymeric Nanocomposite Models The focus of this study is on PNCs consisting of an epoxy resin matrix reinforced either with diglycidyl of bisphenol A(D.E.R.TM 332) or sodium montmorillonite (MMT), respectively. The curing agent is diethyltoluenediamine (Ethacur 100-LC,Albemarle)[19][20][21][22]. Various periodic representative volume elements (RVEs) of the PNCs where the carbon nanotube (CNT) are randomly dispersed and oriented within the polymer matrix are considered. The aspect ratio p = H/D is used to represent the shape of the clay [23], H denoting the height of the clay cylinder and D denotes its diameter. We assume a representative size of the clay as 1µm at the microscale and account for the interface resistance (Kapitza resistance) between inclusions and the matrix. We also study the effect of the volume fraction (λ) 2.1. RVE definition and generation algorithm The essence of the RVE method is to homogenize the original heterogeneous medium and achieve equivalent homogeneous properties in the case of homogenization. It is of major importance to determine the appropriate RVE size and a detailed discussion can be found in [24]. The fillers are randomly dispersed in the matrix according to a specific pre-defined statistical distribution. Fig.1 shows the flowchart. Firstly, we set an initial governing adjustment parameters which need to include the numbers of RVEs for each value of volume fraction, the size of RVE, the aspect ratio, the radius of the fillers and the orientation distribution. We subsequently compute the number of fillers according to the cubic length which is the RVE’s size, the aspect ratio and the randomly generated volume fraction. Exploiting the random sequential addition algorithm (RSA) [25][26], the fillers are placed randomly and sequentially into the cubic matrix. According to their probability density functions (PDFs), the fillers corner positions are determined by different generators and we use a uniform probability distribution function to place the fillers inside the RVE [27]. Fig 2 shows the 3D cubic RVE, where the fillers are different and randomly distributed. In order to enforce periodicity, we check if the filler distribution is suitable for non-overlapping and non-crossing conditions, which means that the portion that is missing on one side of the RVE needs to form a complete filler with the portion of the filler on the other side of the RVE boundary. Typically, there are three types of boundary conditions that can be applied to the RVE, i.e. linear displacement, uniform traction, and hybrid 5

BCs (including periodic BCs). In all our simulations, we employed periodic boundary conditions (PBCs). Therefore, we developed our own program written in C++, which generates these different randomly distributed fillers. The algorithm is based on the Non-collision theory in 3D geometry packing problems. We also developed a program based on Python script to automatically generated the associated discretization within the commercial software package ABAQUS. The code can be downloaded from: https://github.com/jackylbk/Generation-Algorithm-of-PNCs.git. 2.2. Output of the model The output of this model is the effective (or macroscopic) thermal conductivity (ETC). Therefore, we apply two heat fluxes on the X-axis direction surface of the RVE cube, one of which is (+q), applying another heat flux (-q) on its opposite surface. This introduces a heat flow transferring a constant heat flux. The thermal energy is transferred from the hot surface to which is evaluated by the cold surface creating a temperature gradient dT dx computing the differences in the temperature between the hot and the cold surface. The effective thermal conductivity obtained from Fourier Law [28] is given as: dT (1) K =q× dx where q denotes the heat flux and dT is the temperature gradient. Only if dx we average the values of all three axes, the effective thermal conductivity of the RVE can be obtained.

6

Figure 1: Flowchart of RVE generation algorithm

7

Figure 2: 3D cubic RVE

3. Stochastic modelling The main goal of an uncertainty analysis is to quantify the uncertainties in the input parameters w.r.t to an output of interest. The input parameters in our studies will be related ’only’ to material properties while boundary conditions are assumed to be deterministic. This is a reasonable assumptions for ’laboratory experiments’ where boundary conditions can be controlled well. The output of interest is the effective thermal conductivity of the composite (ETC). Iman et al. [29] developed an efficient Latin Hypercube Sampling (LHS) method. The method can be applied to generate multiple input parameters randomly so that the response of the mechanical system can be studied. The variance and the standard deviation of the output parameter in this model can be evaluated once we know the probability density function (PDF) of a set of input parameters. We consider a model with k input parameters. If the cumulative probability curve of each input parameter is partitioned into 8

N equal intervals, a N × k design matrix is obtained with independently distributed random values. Subsequently, we map the generated values to the parameter’s physical distribution to obtain representative values. As shown in [30], the number of requested samples is reduced significantly compared to Monte-Carlo sampling (MCS). Furthermore, since samples are stratified densely across the parameter’s range, the computational cost is reduced considerably. Based on the samples generated by the LHS, surrogate models are employed for the uncertainty analysis to reduce the number of simulations of the mechanical model, which is the most time-consuming part of the entire stochastic analysis. It is well known that the thermal conductivity of single-layer graphene is related to the size of graphene sheets and ranges from K=1000 - 4000 W/m×k at room temperature (RT) [31]. According to Suchismita, et al.[32], the Kapitza resistance is around Kp=2000W/m×k.Usually, the thermal conductivity of graphene is related to the number of atomic planes. When the number of atomic planes is gradually increased, the thermal conductivity of graphene will decrease. In our simulation, 15002500(W/m×k) are chosen as lower and upper bounds based on a truncated normal distribution. According to A. Moisala et al.[33], the thermal conductivity of an epoxy matrix is 0.20(W/m × k). We assume this input as uniform distribution from 0.1 to 0.3 with the mean value of 0.20(W/m×k). Interface problems are also involved in this study, and the interface resistance (RBd ) (Kapitza resistance) between inclusions and the matrix is also considered. According to M.Freitag, et al.[34], the interface resistance is around 25 MW/m2 × K while the computational studies in [35][31] suggest values between 11 and 32 MW/m2 × K. We assume again a uniform distribution from 10 to 30 MW/m2 × K. Mahmood M. Shokrieh, et al.[36] studied CNTs volume fraction of 0.03 to 0.05 employing a normal distribution ranging from 0 to 0.05. We will also study the effect of the geometry of the fillers characterized by the aspect ratio p=H/D. According to Khoa Bui, et al. [37] and Hai M. Duong, et al.[38], 20-70 are considered as lower and upper bounds. They assumed a truncated normal distribution. The random input variables are listed in Table 1. 4. Sensitivity Analysis We use local and variance based sensitivity analysis (SA) to rank the model input parameters and their contribution to the model output. Fig. 9

Table 1: Model uncertainties mean

standard deviation

Thermal conductivity of Fiber (X1 )

2006

227.2

Suchismita, et al.[32]

Thermal conductivity of Matrix(X2 )

0.20

0.013

A. Moisala et al.[33]

20

0.15

Assumed

Volume fraction (X4 )

0.0235

0.01

M.Shokrieh,et al.[36]

Aspect ratio (X5 )

41.96

1.4

Assumed

Inputs

Interface (X3 )

resistance

Type of distribution

Sources

3 shows the structure of all sensitivity analysis methods, which will be described subsequently. 4.1. Scatter plot Scatter plots visualize the relationship between the input and output parameters and provide a first and fast (graphical) though rather qualitative overview. We can also use quantitative measures, by measuring e.g. the correlation between Y and Xi , or variance-based approach in non-linear regression [39]. 4.2. First-order sensitivity indices The function Y = f (X1 , X2 , , Xk ) is the response of the mechanical model. The first-order sensitivity indices are presented as [40] Si =

Vxi [EX∼i (Y |Xi )] V (Y )

(2)

where V (Y ) in the denominator position is regarded as an unconditional variance of Y and Vxi [EX∼i (Y |Xi )] denotes the variance of the expected value of Y conditional on Xi , i.e. Xi is related and fixed to a value Xji , j = 1, . . . , N ; with N being a natural number whose value depends on the number of samples. The Vxi [EX∼i (Y |Xi )] at the molecular position marks the main influence of Xi on the output variable corresponding to the model. 10

START

SAMPLING Conductivity of Fiber(X1) Conductivity of Matrix(X2) Kapitza Resistance (interfacial Resistance) (X3) Volume fraction (X4) Aspect Ratio(X5)

FEM simulations Conductivity(Y)

No

Surrogate Model: Polynomial Regression MLS approximation

Mechanical Model

Statistical Convergence?

Scatter plot

Yes

End

Linear Regression Variancebased methods Si Sti i=1,2,3,4,5

Quadratic without mixed terms Regression Variancebased methods Si Sti i=1,2,3,4,5

Full Quadratic Regression Variancebased methods Si Sti i=1,2,3,4,5

MLS Regression Variancebased methods Si Sti i=1,2,3,4,5

End

End

End

End

Figure 3: Schematic diagram of all sensitivity assessment methods

11

4.3. Total effect sensitivity indices The first-order sensitivity index has certain limitations, mainly in the partial measurement of the variance of the model.Since that, an extension to higher-order coupling indices is essential. Therefore, the total-effect indices are commonly employed. A definition of the total effect index is made, the formula is as follows [40] STi =

Vx∼i [EXi (Y |X∼i )] Ex∼i [VXi (Y |X∼i )] =1− V (Y ) V (Y )

(3)

where Vx∼i [EXi (Y |X∼i )] at the molecular position is the variance which it is derived from of the expected value of all parameters Y we need in this mode except for Xi , which subsequently denoted as X∼i . On the model output the first-order effect of X∼i is represented by Vx∼i [EXi (Y |X∼i )], which the effect part of Xi is not included in this item. The overall effect of input parameter Xi on the output is represented by an index of total effect STi , i.e. the total effect is expressed by the sum of the first-order term and all high-order terms: STi = Si + Si ,∼i = 1 − S∼i

(4)

where S∼i is a sum, including all parameters’ sensitivity indices except of i. The items Si and STi in Eqs [2] and [3] require numerous samples, which can be generated for instance by LHS. The cost of applying these measures is usually large. Thus, a model Y which is called surrogate model is used to estimate and represent the response of the ’real’ model. 4.4. Variance-based method using Sobol‘ quasi-random sequences for computing Si and STi Let us summarize the Monte Carlo estimation based on the Sobol’ indices [41]. The matrix A and B are two independent samples with corresponding entries aji and bji,both have the same dimension for the input parameters (i) (i) X, and the matrix is represented as (N × k) . Matrix AB BA are defined with entries from A(B) except the i-th column, which is taken from B(A)

12

(i)

(i)

and compute Si from A, BA or B, AB [42]: N 1 X (i) Vxi [EX∼i (Y |Xi )] = f (A)j f (BA )j − f02 N j=1

(5)

where (B)j denotes the j-th row of matrix B. We rebuilt STi in Eqs [3] with Vx∼i [EXi (Y |X∼i )] as follows: Vx∼i [EXi (Y |X∼i )] =

N 1 X (i) f (A)j f (AB )j − f02 N j=1

(6)

Similarly, the estimation of Si from Eqs[5] can be derived [43] [44]: Vxi [EX∼i (Y |Xi )] =

N 1 X (i) f (A)j (f (BA )j − f (B)j ) N j=1

(7)

Saltelli et al.[41]suggested a more appropriate improvement under the framework of the variance-based approach based on Quasi-Monte Carlo sam(i) (i) ples where the triplet A,B and AB are used instead of B, A and BA in Eqs[7]. It can be shown that the STi estimator in Eq [6] is obtained by [45] Vx∼i [EXi (Y |X∼i )] =

N 1 X (i) f (A)j (f (A)j − f (AB )j ) N j=1

(8)

An alternative formula for the estimators of Si and STi has been proposed by Jansen [46] where Vxi [EX∼i (Y |Xi )] is obtained from Vxi [EX∼i (Y |Xi )] =

N 1 X (i) (f (B)j − f (AB )j )2 2N j=1

(9)

After using Ex∼i [VXi (Y |X∼i )] in Eqs[3] STi is given by Vx∼i [EXi (Y |X∼i )] =

N 1 X (i) (f (A)j − f (AB )j )2 2N j=1

(10)

According to the theory from Saltelli, he suggested that using the quasi(i) random numbers (QRN) to generate A and AB , and then compute STi estimates. The quasi-random sequences algorithm(QRSA) i.e.LP τ [47] has been 13

implemented in MATLAB. It is applied to generate the samples X1 , X2 , ......, Xk as uniformly as possible over the unit hypercube ω. There are three main requirements on LP τ sequences [40]: 1. The uniformity of the distribution is optimal only if the sequence length is increased. 2. This distribution can also be applied to a relatively small initial set, when the scale is small, that is, when the conditional N value is relatively small. 3. Computational cost is not very expensive with this algorithm. After we first sample the quasi-random sequence with matrix size (N × 2k) by QRSA method, we divide it into matrix A (left half) and B (right (i) half) with matrix size (N × k).The triplet A, B and AB include a higher amount of ’ideal’ points, which is better than others. In order to compute Si (i) and STi ,they are ’assigned’ to the Triplet A, B and BA . 4.5. Surrogate Model The surrogate model [23][48][49] [50] [51] is an effective method to approximate the method of the ’real’ model. In this study, the popular polynomial regression method(with and without mixed terms) and moving least squares (MLS) approximation are applied. 4.5.1. Polynomial Regression model There is a basic function which can approximate the ’real’ response Y in the polynomial regression model (both linear and quadratic with mixed terms): T PX (X) = [X1

. . . X5 . . . X12

. . . X52 . . . X1 X2

. . . X1 X5 . . . ]

(11)

Then, the SA methods yield the sensitivity indices for different input parameters; based on the basis of linear polynomial, which are integrated T into vector a vector PX (X) = [X1 . . . Xk ].We use Yˆ to approximate the response of the ’real’ model Y that with k-parameter. Yˆ = β0 +

K X

T βi PX (Xi ) + e

i=1

14

(12)

The vector β can be computed in a kind of way that we minimize the mean squared difference (MSD) S between the measured value at N training points in ’real’ model xj = [xj1 , . . . , xjk ], j = 1, . . . , N and the surrogate model’s value. There is a least-square estimate (LSE) βˆ which is given by T T βˆ = (PX PX )−1 PX Y

(13)

There is a versatile metric to estimate model’s accuracy which is defined as coefficient of determination (COD) R2 and it is given by: PN ˆ 2 j=1 (Yj − Yj ) 2 (14) R = 1 − PN 2 j=1 (Yj − Y )

The main purpose of R2 is to describe the similarity between those two models, the the ’real’ model and the the surrogate model. In addition, the R2 2 improvement is the adjusted COD Radj . The value of purchasing coefficients kR and the value of supporting points N are also considered into this measure method.

(N − 1) (1 − R2 ) (15) (N − kR ) Not only the regression in linear polynomials, we also use the quadratic without mixed terms basis function in this research: 2 Radj =1−

PXT (X) = [X1

X2

. . . Xk

X12

X22

. . . Xk2 ]

(16)

The corresponding regression model is given as follows: Yˆ = β0 +

K X (βi Xi + βii Xi2 ) + e

(17)

i=1

4.5.2. Moving Least Squares (MLS) regression model In the MLS regression model, a local weighting functions w(d) (Lancaster and Salkauskas, 1986 [52]) is introduced in order to calculate a set of unordered point samples. We obtain the unknown coefficient by minimizing a difference value which is from the interpolation point x in local approximation and the supporting point xi in input part of the real model collection. The function YˆM LS is given as below: YˆM LS = P T (X)βw 15

(18)

Linear or quadratic terms are commonly employed in the polynomial basis p: T PX (X) = [1 X1

X2

. . . Xk

X12

X22

. . . Xk2 ]

(19)

After minimizing a discrete weighted L2 -norm LW , we can obtain the weighted least squares coefficient βˆW LW = (Y − pT(X) βW )T W (X)(Y − pT (x)βW )

(20)

with the unknown coefficients βw yielding ∂LW ∂βw

= 0 → βˆW = (X T W (x)X)−1 X T W (x)Y

where W (x) is a diagonal matrix  w(x − x1 ) 0  0 w(x − x2 )   ... ... 0 0

 ... 0  ... 0   ... ... . . . w(x − xN )

(21)

(22)

the cubic polynomial weighting function is applied in this paper [53]: ( 1 − 3s2 + 2s3 s 6 1 w(s) = (23) 0 s>1

where s = ||x − xi ||/D is a distance that is normalized between the supporting point and the interpolation point; D is the influential radius that is assumed to be a certain value (fixed) in our examples. Substituting Eq. 21 into Eq. 18 the MLS approximation can be obtained shown as: YˆM LS = pT (X)[X T W (x)X]−1 X T W (x)Y

(24)

4.6. Normalization of the input Shrinking the data to a certain interval [0, 1] is called data normalization. It eliminates the unit limit of the data and converts or changes it to a dimensionless value, which allows for better comparison and weighting of different units or different sized indicators. Normalization can improve the convergence speed of the model and the accuracy of the model. 16

We need to normalize our input before building the surrogate models since the input parameters have different scales in our computation. The i-th input factor Xi , the normalization produced by the s-th is shown as: (s)

(S),norm

Xi

=

Xi − min(Xi ) minX(Xi ) − max(Xi )

(25)

where max(Xi ) and min(Xi ) represent the maximum and minimum values of the i-th input sample, respectively. 5. Numerical Results In this section, we first show the convergence of our Finite Element RVE model. As the size increases, the thermal conductivity of the model also fluctuates due to the uncertainty of the stochastic model. Fig 4 shows the resulting thermal conductivity versus RVE size. Let us name the output of interest, i.e. the thermal conductivity, Y . We consider five uncertain input parameters, i.e. the thermal conductivity of the fiber (X1 ), the thermal conductivity of the matrix (X2 ), the interface resistance (X3 ), the volume fraction (X4 ) and the aspect ratio (X5 ). Fig 5 shows the temperature distribution inside the RVE. The output uncertainty can be represented by a probability distribution function illustrated in Figs. 6, 7 and 8. According to the shape of the curve fitting in these figures, the Log-normal distribution, Normal distribution and Weibull distribution are listed in Table 2, which are proper candidates for the output uncertainty. We use probability plots which are based on a best fit to get the parameters in this table. The probability plot is a useful tool to assess in visual if the output data can be represented by a certain distribution. When the minimum deviations from the linear line in the probability plot is generated, the assumed distribution is taken as a goodness-of-fit of those data. Figs. 9, 10 and 11 show the Log-normal, Normal and Weibull distribution probability plots. The Lognormal distribution is better suitable than the other distributions. Since we need a quantitative assessment of the best fit of these data, the AndersonDarling test is applied [54]. The Anderson Darling statistic determines how well the data follows a particular distribution. For a specified data set and distribution, the better the distribution fits the data, the smaller this statistic will be [55][56]. It is characterized by the AD and P-value summarized in table 2. According to the values in this table, the smallest AD appears for 17

the Log-normal distribution indicating the Log-normal distribution indeed fits those data best. Also we should note that the AD values are estimated in a certain area which is within 95% confidence intervals (CIs). Table 2: Uncertainties of mechanical output using different distribution

Type of as- Parameters Parameters sumed PDF 1 2

Anderson- P-Value Darling(AD)

Log-normal PDF

Location -0.6177

0.378

0.403

Normal PDF

Mean value Standard 0.6007 deviation 0.2835

3.757

<0.005

Weibull PDF

shape 2.275

1.994

<0.005

Scale 0.4712

Scale 0.6810

Figure 4: Thermal conductivity versus RVE size

18

Figure 5: Result of 3D cubic RVE

19

Figure 6: Histogram of mechanical output and assumed probability density functions(PDFs) for Log-normal distribution

20

Figure 7: Histogram of mechanical output and assumed probability density functions(PDFs) for Normal distribution

21

Figure 8: Histogram of mechanical output and assumed probability density functions(PDFs) for Weibull distribution

22

Figure 9: Probability plot for Log-normal distribution of output data

23

Figure 10: Probability plot for Normal distribution of output data

24

Figure 11: Probability plot for Weibull distribution of output data

5.1. Scatter plots The scatter plots in Figs. 12, 13, 14, 15 and 16 show the influence of the thermal conductivity of the fiber (X1 ), the thermal conductivity of the matrix (X2 ), interface resistance (X3 ), volume fraction (X4 ) and aspect ratio(X5 ) on the homogenized thermal conductivity (Y) of the composite. Obviously X4 and X5 significantly increase when the X-value is growing; X1 and X2 have upward trends, and there is almost no trend in X3 . A quantitative assessment is subsequently made through a sensitivity analysis.

25

Figure 12: Scatter plots of X1

Figure 13: Scatter plots of X2

26

Figure 14: Scatter plots of X3

Figure 15: Scatter plots of X4

27

Figure 16: Scatter plots of X5

5.2. Regression Model In the following, we use 200 samples to build surrogate models to simulate and approximate the output results of the ’real’ model. According to the theory of linear, quadratic, and MLS regression, the value of COD R2 2 and adjusted COD Radj are calculated. We find that those results are good approximations of the response of the ’real’ model, see Table 3. The COD is 0.8 for the linear regression, 0.82 for the quadratic one without mixed terms, 0.83 for the full quadratic one and 0.92 for the MLS regression indicating that the MLS regression is the best though the other models are within general accepted tolerances [40]. 5.3. Global sensitivity analysis After we identified the specific surrogate model, we generated 100,000 Latin hypercube samples in this surrogate model. Table 3 shows the regression models and detailed coefficient in every model.According to the Eqs.2 and 3, the Si and STi , i = 1, 2, 3, 4, 5, are calculated in MLS approximation method and polynomial regression. R% reflects the approximation of the surrogate model to the ’real’ model to a certain extent, so the value of the sensitivity index is reduced. We define those reduced sensitivity indices Sˆi and SˆTi as: Sˆi = R2 Si 28

(26)

2 Table 3: Regression coefficient COD R2 and adjusted COD Radj

Regression Model

Regression Coefficient

Linear regression

β0 β1 β2 β3 β4 β5 β0

R2 = 0.8 2 Radj = 0.8 Quadratic without mixed terms regression

= 1.1515 = 2.23e − 5 = −2.85e − 4 = 7.02e − 3 = 1.56e − 2 = 6.24e + 5 = −1.6338

β1 = 5.32e − 4 β2 = 5.23e − 4 β3 = 8.37e − 4 β4 = 6.08e − 2 β5 = −1.50e − 2 β11 = −1.35e − 07 β22 = 2.85e − 06 β33 = −7.e − 4 β44 = 2.43e − 2 β55 = 2.87e − 3 β0 = −1.2526 β1 = 6.21e − 4 β2 = 2.63e + 5 β3 = 4.73e − 2 β4 = 7.24e + 5 β5 = 6.12e − 4 β11 = −3.57e − 6 β22 = 9.34e − 4 β33 = 5.90e − 5 β44 = −3.75e − 7 β55 = 6.45e − 2 β12 = −5.72e − 5 β13 = 5.35e + 3 β14 = 2.56e − 2 β15 = 4.76e − 3 β23 = 3.75e − 4 β24 = 1.25e + 4 β25 = 5.27e − 2 β34 = 1.11e − 5 β35 = 6.31e − 6 β45 = 2.65e + 2

R2 = 0.82 2 Radj = 0.82 Full quadratic regression

R2 = 0.83 2 Radj = 0.83 Moving Least Square R2 = 0.92 2 Radj = 0.92

29

SˆTi = R2 STi

(27)

Table 4 shows the results of the variance-based methods. The key results of the SA can be summarized as follows: 1. The total-effect indices SˆTi and the reduced first-order indices Sˆi exhibit only subtle differences. Or in other words: There are almost no correlations between the input-parameters. 2. All surrogate models are suitable and predict the same trend. The highest COD is obtained by the MLS approach (0.92), followed by the quadratic (with and without mixed terms) regression model (0.83 and 0.82) and the linear regression model (0.80). 3. The aspect ratio and volume fraction have the largest influence on the overall thermal conductivity of the composite followed by the thermal conductivity of the matrix and conductivity of the fiber. The kapitza resistance (=interfacial resistance) has almost no effect on the overall thermal conductivity of the material indicating that it might be neglected in the model. This is an interesting finding also considering that it is difficult to extract the interface resistance experimentally.

30

Table 4: First-order and total effects sensitivity indices computed on different surrogate models

Linear Regression Model

Quadratic with- Full Quadratic out mixed terms Regression

MLS

First-order indices Sˆi Sˆ1 =0.02 Sˆ2 =0.1 Sˆ3 =0 Sˆ4 =0.34 Sˆ5 =0.42 P5 ˆ i=1 Si =0.88 Total-effect indices SˆT i SˆT 1 = 0.02 SˆT 2 = 0.12 SˆT 3 = 0 SˆT 4 = 0.34 SˆT 5 = 0.43 P5 ˆ i=1 ST i = 0.91

First-order indices Sˆi Sˆ1 =0.02 Sˆ2 =0.12 Sˆ3 =0 Sˆ4 = 0.34 Sˆ5 = 0.42 P5 ˆ i=1 Si =0.9 Total-effect indices SˆT i SˆT 1 =0.02 SˆT 2 =0.12 SˆT 3 =0 SˆT 4 =0.34 SˆT 5 =0.43 P5 ˆ i=1 ST i =0.91

First-order indices Sˆi Sˆ1 =0.02 Sˆ2 =0.12 Sˆ3 =0 Sˆ4 =0.35 Sˆ5 =0.45 P5 ˆ i=1 Si =0.94 Total-effect indices SˆT i SˆT 1 =0.02 SˆT 2 =0.13 SˆT 3 =0 SˆT 4 =0.36 SˆT 5 =0.44 P5 ˆ i=1 ST i =0.95

First-order indices Sˆi Sˆ1 =0.02 Sˆ2 =0.12 Sˆ3 =0 Sˆ4 =0.34 Sˆ5 =0.43 P5 ˆ i=1 Si =0.91 Total-effect indices SˆT i SˆT 1 =0.02 SˆT 2 =0.12 SˆT 3 =0 SˆT 4 =0.34 SˆT 5 =0.43 P5 ˆ i=1 ST i =0.91

6. Conclusions and future work 1.According to the data of ST i and Si , the aspect ratio and volume fraction have the largest influence on the total thermal conductivity of the composites. For this result, we separately explain several aspects from geometry and material mechanics and model generation. First, from a geometric point of view, increasing the aspect ratio can significantly improve the effective part of the fillers, which is equivalent to indirectly increasing the unit volume[57]. The same effect can be achieved by increasing the volume fraction, but the degree of influence is less compared to the aspect ratio. Although the thermal conductivity of the matrix itself is quite low, its volume is relatively large and hence it will also have a certain degree of impact on the overall composite material. The thermal conductivity of carbon fibers is very high, but the content of carbon fibers in composite materials is very low so that this parameter does not have a significant influence on the whole 31

material for the scenarios studied here. Previous research had the same trend in their study[31][58]. On the other hand, from the material physical theory, the substance is prone to phase change at a certain voltage, pressure or temperature[59]. Previous studies by scholars have shown that carbon nano-fibers may have a phase change at a certain temperature, and will produce a liquid crystal state[60][61][62]. In this case, a large amount of heat is absorbed but the temperature is kept constant, which will greatly improve the overall heat conductivity of the composite material. In addition, previous scholars have shown that elongated carbon nano-fiber tubes produce nematic phase transitions even at relatively small volume fractions[63][64]. This phase transition can significantly improve the physical properties of carbon nano-fibers, making thermal conductivity greatly improved,that has a obvious impact on the overall thermal conductivity of the material[65]. Finally, regarding the model generation, the model filler we used in this study were generated by RSA. According to this algorithm, the fibers in the model are randomly added. The structural statistical parameters such as fiber orientation distribution and spatial correlation function were not explicitly constrained during the generation. That makes the actual fibers selected were inevitably related to existing ones, which lead to orientation correlation[66]. The aspect ratio in the input parameters is particularly sensitive to orientation correlation, which is most likely a reason for observing that the aspect ratio is the most sensitive parameter for thermal conductivity. 2. In any regression model, there is almost no difference between ST i and Si . According to the theory of global sensitivity analysis, this shows that there will be no interaction between the input parameters. 3. The Kapitza resistance effect is sufficient to be considered at the nanoscale, and the effect is more pronounced, but the effect on the meso-scale is extremely weak and almost negligible. The results shows it has almost no effect on the overall thermal conductivity of the material at meso-scale. At nano-scale the interface resistance might have a lager influence, but in this RVE model, it might be neglected. 4. Comparing different regression models, the MLS has the highest precision followed by the full quadratic polynomial regression and the quadratic polynomial regression without mixed term and finally the linear polynomial regression. However, the same conclusions can be derived from all surrogate models. 5. The thermal-mechanical properties of PNCs have been obtained by 32

using FEM simulations. The output (thermal conductivity) is best characterized by a Log-normal distribution. In the future,we will extend our work to Polymeric clay nano-composites, especially in model generation. In our study, the RSA method was used to generate fillers. For structural statistics, such as fiber orientation distribution, spatially related functions did not receive explicit constraints during the generation process. One possible way to rigorous control the orientation order is to use inverse Monte Carlo method to further optimize the RSA packing[66].So in future work, we need to further consider the constraints of statistical data in the algorithm of filler generation. Acknowledgement We gratefully acknowledge the support of the China Scholarship Council (CSC). References References [1] B. Mortazavi, J. Bardon, S. Ahzi, Interphase effect on the elastic and thermal conductivity response of polymer nanocomposite materials: 3d finite element study, Computational Materials Science 69 (2013) 100 – 106. [2] B. Li, H. Fang, H. He, K. Yang, C. Chen, F. Wang, Numerical simulation and full-scale test on dynamic response of corroded concrete pipelines under multi-field coupling, Construction and Building Materials 200 (2019) 368–386. [3] H. Fang, J. Lei, M. Yang, Z. Li, Analysis of gpr wave propagation using cuda-implemented conformal symplectic partitioned runge-kutta method, Complexity 2019 (2019). [4] T. Mori, K. Tanaka, Average stress in matrix and average elastic energy of materials with misfitting inclusions, Acta metallurgica 21 (1973) 571– 574. [5] A. K. Sen, S. Torquato, Effective conductivity of anisotropic two-phase composite media, Physical Review B 39 (1989) 4504. 33

[6] G. Liu, A step-by-step method of rule-of-mixture of fiber-and particlereinforced composite materials, Composite structures 40 (1997) 313–322. [7] N. Sheng, M. C. Boyce, D. M. Parks, G. Rutledge, J. Abes, R. Cohen, Multiscale micromechanical modeling of polymer/clay nanocomposites and the effective clay particle, Polymer 45 (2004) 487–506. [8] S. Zhai, P. Zhang, Y. Xian, P. Yuan, D. Yang, Modelling and analysis of effective thermal conductivity for polymer composites with sheet-like nanoparticles, Journal of Materials Science 54 (2019) 356–369. [9] A. Pourasghar, Z. Chen, Hyperbolic heat conduction and thermoelastic solution of functionally graded cnt reinforced cylindrical panel subjected to heat pulse, International Journal of Solids and Structures 163 (2019) 117–129. [10] B. Mortazavi, M. Baniassadi, J. Bardon, S. Ahzi, Modeling of two-phase random composite materials by finite element, mori–tanaka and strong contrast methods, Composites Part B: Engineering 45 (2013) 1117 – 1125. [11] X. Yang, Y. Wan, X. Wang, Y. Fu, Z. Huang, Q. Xie, Molecular dynamics studies of the mechanical behaviors and thermal conductivity of the dgeba/mthpa/cnb composites, Composites Part B: Engineering (2019). [12] Y. Liu, A.-L. Hamon, P. Haghi-Ashtiani, T. Reiss, B. Fan, D. He, J. Bai, Quantitative study of interface/interphase in epoxy/graphenebased nanocomposites by combining stem and eels, ACS applied materials & interfaces 8 (2016) 34151–34158. [13] W. Ding, A. Eitan, F. Fisher, X. Chen, D. Dikin, R. Andrews, L. C. Brinson, L. Schadler, R. Ruoff, Direct observation of polymer sheathing in carbon nanotube- polycarbonate composites, Nano letters 3 (2003) 1593–1597. [14] C. Zhang, J. Zhao, T. Rabczuk, The interface strength and delamination of fiber-reinforced composites using a continuum modeling approach, Composites Part B: Engineering 137 (2018) 225–234.

34

[15] Y. Xian, P. Zhang, S. Zhai, P. Yuan, D. Yang, Experimental characterization methods for thermal contact resistance: a review, Applied Thermal Engineering 130 (2018) 1530–1548. [16] S. Zhou, T. Rabczuk, X. Zhuang, Phase field modeling of quasi-static and dynamic crack propagation: Comsol implementation and case studies, Advances in Engineering Software 122 (2018) 31–49. [17] S. Zhou, X. Zhuang, T. Rabczuk, Phase-field modeling of fluid-driven dynamic cracking in porous media, Computer Methods in Applied Mechanics and Engineering 350 (2019) 169–198. [18] S. Zhou, X. Zhuang, T. Rabczuk, Phase field modeling of brittle compressive-shear fractures in rock-like materials: A new driving force and a hybrid formulation, Computer Methods in Applied Mechanics and Engineering 355 (2019) 729–752. [19] K. Wang, L. Chen, J. Wu, M. L. Toh, C. He, A. F. Yee, Epoxy nanocomposites with highly exfoliated clay: mechanical properties and fracture mechanisms, Macromolecules 38 (2005) 788–800. [20] M. Silani, S. Ziaei-Rad, M. Esfahanian, V. Tan, On the experimental and numerical investigation of clay/epoxy nanocomposites, Composite Structures 94 (2012) 3142–3148. [21] L. Wang, K. Wang, L. Chen, C. He, Y. Zhang, Hydrothermal effects on the thermomechanical properties of high performance epoxy/clay nanocomposites, Polymer Engineering & Science 46 (2006) 215–221. [22] K. Wang, L. Wang, J. Wu, L. Chen, C. He, Preparation of highly exfoliated epoxy/clay nanocomposites by slurry compounding: process and mechanisms, Langmuir 21 (2005) 3613–3618. [23] N. Vu-Bac, T. Lahmer, Y. Zhang, X. Zhuang, T. Rabczuk, Stochastic predictions of interfacial characteristic of polymeric nanocomposites (pncs), Composites Part B: Engineering 59 (2014) 80–95. [24] N. Vu-Bac, T. Rabczuk, X. Zhuang, Continuum/finite element modeling of carbon nanotube–reinforced polymers, in: Carbon NanotubeReinforced Polymers, Elsevier, 2018, pp. 385–409. 35

[25] G. Zhang, S. Torquato, Precise algorithm to generate random sequential addition of hard hyperspheres at saturation, Physical Review E 88 (2013) 053312. [26] Y. Pan, L. Iorga, A. A. Pelegri, Numerical generation of a random chopped fiber composite rve and its elastic properties, Composites Science and Technology 68 (2008) 2792–2798. [27] M. Silani, H. Talebi, S. Ziaei-Rad, P. Kerfriden, S. P. Bordas, T. Rabczuk, Stochastic modelling of clay/epoxy nanocomposites, Composite Structures 118 (2014) 241–249. [28] B. Mortazavi, J. Bardon, S. Ahzi, Interphase effect on the elastic and thermal conductivity response of polymer nanocomposite materials: 3d finite element study, Computational Materials Science 69 (2013) 100– 106. [29] R. L. Iman, W. Conover, Small sample sensitivity analysis techniques for computer models. with an application to risk assessment, Communications in statistics-theory and methods 9 (1980) 1749–1842. [30] D. Nov´ak, B. Tepl` y, Z. Kerˇsner, The role of latin hypercube sampling method in reliability engineering, in: Proc. of ICOSSAR, volume 97, pp. 403–409. [31] B. He, B. Mortazavi, X. Zhuang, T. Rabczuk, Modeling kapitza resistance of two-phase composite material, Composite Structures 152 (2016) 939–946. [32] S. Ghosh, W. Bao, D. L. Nika, S. Subrina, E. P. Pokatilov, C. N. Lau, A. A. Balandin, Dimensional crossover of thermal transport in few-layer graphene, Nature materials 9 (2010) 555. [33] A. Moisala, Q. Li, I. Kinloch, A. Windle, Thermal and electrical conductivity of single-and multi-walled carbon nanotube-epoxy composites, Composites science and technology 66 (2006) 1285–1288. [34] M. Freitag, M. Steiner, Y. Martin, V. Perebeinos, Z. Chen, J. C. Tsang, P. Avouris, Energy dissipation in graphene field-effect transistors, Nano letters 9 (2009) 1883–1888. 36

[35] B. Mortazavi, O. Benzerara, H. Meyer, J. Bardon, S. Ahzi, Combined molecular dynamics-finite element multiscale modeling of thermal conduction in graphene epoxy nanocomposites, Carbon 60 (2013) 356–365. [36] M. M. Shokrieh, R. Rafiee, Stochastic multi-scale modeling of cnt/polymer composites, Computational Materials Science 50 (2010) 437–446. [37] K. Bui, B. P. Grady, D. V. Papavassiliou, Heat transfer in high volume fraction cnt nanocomposites: Effects of inter-nanotube thermal resistance, Chemical Physics Letters 508 (2011) 248–251. [38] H. M. Duong, D. V. Papavassiliou, K. J. Mullen, B. L. Wardle, S. Maruyama, Calculated thermal properties of single-walled carbon nanotube suspensions, The Journal of Physical Chemistry C 112 (2008) 19860–19865. [39] P. Paruolo, M. Saisana, A. Saltelli, Ratings and rankings: voodoo or science?, Journal of the Royal Statistical Society: Series A (Statistics in Society) 176 (2013) 609–634. [40] A. Saltelli, M. Ratto, T. Andres, F. Campolongo, J. Cariboni, D. Gatelli, M. Saisana, S. Tarantola, Global sensitivity analysis: the primer, John Wiley & Sons, 2008. [41] A. Saltelli, P. Annoni, I. Azzini, F. Campolongo, M. Ratto, S. Tarantola, Variance based sensitivity analysis of model output. design and estimator for the total sensitivity index, Computer Physics Communications 181 (2010) 259–270. [42] I. M. Sobol, Sensitivity estimates for nonlinear mathematical models, Mathematical modelling and computational experiments 1 (1993) 407– 414. [43] A. Saltelli, Making best use of model evaluations to compute sensitivity indices, Computer physics communications 145 (2002) 280–297. [44] S. Tarantola, D. Gatelli, S. Kucherenko, W. Mauntz, et al., Estimating the approximation error when fixing unessential factors in global sensitivity analysis, Reliability Engineering & System Safety 92 (2007) 957–960. 37

[45] I. M. Sobol’, Global sensitivity indices for the investigation of nonlinear mathematical models, Matematicheskoe Modelirovanie 17 (2005) 43–52. [46] M. J. Jansen, Analysis of variance designs for model output, Computer Physics Communications 117 (1999) 35–43. [47] S. Kucherenko, High dimensional Sobols sequences and their application, Technical Report, Technical Report, 2008. [48] N. Vu-Bac, T. Lahmer, H. Keitel, J. Zhao, X. Zhuang, T. Rabczuk, Stochastic predictions of bulk properties of amorphous polyethylene based on molecular dynamics simulations, Mechanics of Materials 68 (2014) 70 – 84. [49] N. Vu-Bac, M. Silani, T. Lahmer, X. Zhuang, T. Rabczuk, A unified framework for stochastic predictions of mechanical properties of polymeric nanocomposites, Computational Materials Science 96 (2015) 520 – 535. Special Issue Polymeric Composites. [50] N. Vu-Bac, R. Rafiee, X. Zhuang, T. Lahmer, T. Rabczuk, Uncertainty quantification for multiscale modeling of polymer nanocomposites with correlated parameters, Composites Part B: Engineering 68 (2015) 446 – 464. [51] N. Vu-Bac, T. Lahmer, X. Zhuang, T. Nguyen-Thoi, T. Rabczuk, A software framework for probabilistic sensitivity analysis for computationally expensive models, Advances in Engineering Software 100 (2016) 19 – 31. [52] P. Lancaster, K. Salkauskas, An introduction: curve and surface fitting, Academic, San Diego (1986). [53] T. Most, C. Bucher, A moving least squares weighting function for the element-free galerkin method which almost fulfills essential boundary conditions, Structural Engineering and Mechanics 21 (2005) 315–332. [54] T. W. Anderson, D. A. Darling, A test of goodness of fit, Journal of the American Statistical Association 49 (1954) 765–769. [55] N. M. Razali, Y. B. Wah, et al., Power comparisons of shapiro-wilk, kolmogorov-smirnov, lilliefors and anderson-darling tests, Journal of statistical modeling and analytics 2 (2011) 21–33. 38

[56] F. W. Scholz, M. A. Stephens, K-sample anderson–darling tests, Journal of the American Statistical Association 82 (1987) 918–924. [57] F. Deng, Q.-S. Zheng, L.-F. Wang, C.-W. Nan, Effects of anisotropy, aspect ratio, and nonstraightness of carbon nanotubes on thermal conductivity of carbon nanotube composites, Applied Physics Letters 90 (2007) 021914. [58] F. H. Gojny, M. H. Wichmann, B. Fiedler, I. A. Kinloch, W. Bauhofer, A. H. Windle, K. Schulte, Evaluation and identification of electrical and thermal conduction mechanisms in carbon nanotube/epoxy composites, Polymer 47 (2006) 2036–2045. [59] A. Goncharuk, N. Lebovka, L. Lisetski, S. Minenko, Aggregation, percolation and phase transitions in nematic liquid crystal ebba doped with carbon nanotubes, Journal of Physics D: Applied Physics 42 (2009) 165411. [60] S. Zhang, S. Kumar, Carbon nanotubes as liquid crystals, Small 4 (2008) 1270–1283. [61] I. A. Kinloch, J. Suhr, J. Lou, R. J. Young, P. M. Ajayan, Composites with carbon nanotubes and graphene: An outlook, Science 362 (2018) 547–553. [62] N. Lebovka, T. Dadakova, L. Lysetskiy, O. Melezhyk, G. Puchkovska, T. Gavrilko, J. Baran, M. Drozd, Phase transitions, intermolecular interactions and electrical conductivity behavior in carbon multiwalled nanotubes/nematic liquid crystal composites, Journal of Molecular Structure 887 (2008) 135–143. [63] H. Duran, B. Gazdecki, A. Yamashita, T. Kyu, Effect of carbon nanotubes on phase transitions of nematic liquid crystals, Liquid crystals 32 (2005) 815–821. [64] W. Song, A. H. Windle, Isotropic- nematic phase transition of dispersions of multiwall carbon nanotubes, Macromolecules 38 (2005) 6181– 6188.

39

[65] E. Lafuente, M. Pi˜ nol, L. Oriol, E. Mu˜ noz, A. M. Benito, W. K. Maser, A. Dalton, J. L. Serrano, M. T. Martinez, Polyazomethine/carbon nanotube composites, Materials Science and Engineering: C 26 (2006) 1198– 1201. [66] L. Liu, Z. Li, Y. Jiao, S. Li, Maximally dense random packings of cubes and cuboids via a novel inverse packing method, Soft matter 13 (2017) 748–757.

40

Author Statement Bokai Liu: Formal analysis, Validation, Visulaization Nam Vu-Bac: Methodology Xiaoying Zhuang: Conceptualization, Supervision Timon Rabczuk: Conceptualization, Supervision, Project Administration

41

Declaration of Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

42