Numerical Optimization of CNT Distribution in Functionally Graded CNT-Reinforced Composite Beams

This paper is concerned with the numerical optimization of the thickness-wise CNT (carbon nanotube) distribution in functionally graded CNT-reinforced composite (FG-CNTRC) beams to secure the structural safety. The FG-CNTRC in which CNTs are inserted according to the specific thickness-wise distribution pattern are extensively investigated for high-performance engineering applications. The mechanical behaviors of FG-CNTRC structures are definitely affected by the distribution pattern of CNTs through the thickness. Hence, the tailoring of suitable CNT distribution pattern is an essential subject in the design of FG-CNTRC structure for a given boundary and loading conditions. Nevertheless, the thickness-wise CNT distribution pattern has been assumed by several linear functions so that these assumed primitive patterns cannot appropriately respond to arbitrary loading and boundary conditions. In this context, this paper aims to introduce a numerical method for optimally tailoring the CNT distribution pattern of FG-CNTRC beams. As a preliminary stage, the effective stress is defined as the objective function and the layer-wise CNT volume fractions are chosen as the design variables. The exterior penalty-function method and golden section method are adopted for the optimization formulation, together with finite difference scheme for the design sensitivity analysis. The proposed optimization method is illustrated and validated through the benchmark experiments, such that it successfully provides an optimum CNT distribution which can significantly minimize the effective stress, with a stable and rapid convergence in the iterative optimization process.


Introduction
Carbon nanotube (CNT) has been spotlighted as a state-of-art material of the 21st century due to its excellent physical, electrical and chemical properties [1,2]. A representative quality is the high weight-stiffness ratio, so that CNTs can be used as a next-generation pillar for polymer-matrix composites. In fact, the mechanical strength of the CNT-reinforced polymer composite dramatically increases when only a small amount of CNTs are inserted [3]. In mechanical applications, CNTRCs have been developed in the form of beams, plates and shells, because these are used as ingredient components in a variety of engineering applications. When compared with the conventional composites, CNTRCs provide superior mechanical behaviors, such as bending deformation, free vibration and buckling [4][5][6].
Meanwhile, the above-mentioned structural components exhibit variations in their mechanical behaviors through their thickness, so that the uniform thickness-wise distribution of CNTs may not be sufficient to respond to such variations. It has been reported that the improvement in the mechanical behavior of CNTRCs is limited when CNTs are uniformly distributed [7,8]. To resolve this limitation, the functional gradient in the thickness-wise CNT distribution was introduced according to the concept of functionally graded material (FGM) [9]. Shen [10] and Ke et al. [11] proposed a purposeful thickness-wise CNT become the design variables and their arithmetic sum should satisfy the preset entire CNT volume fraction. The constrained numerical optimization is formulated by employing the exterior penalty-function method, and the sensitivity of objective function to the design variable vector is evaluated using the central difference scheme. The proposed numerical optimization method is validated through the benchmark examples, and the optimum CNT distributions and the associated stress distributions are investigated with respect to the loading condition and compared with those of conventional CNT distribution patterns. The numerical experiments confirmed that the present optimization method successfully seeks the optimum CNT distribution, which provides the minimum effective stress.
The paper is organized as follows: FG-CNTRC structures are introduced in Section 2, together with the CNT volume fraction distribution and the effective material properties. The numerical formulations for the finite-element static analysis and the CNT distribution optimization are addressed in Section 3. The optimization results are presented in Section 4, together with the discussion, and the final conclusion is given in Section 5. Figure 1a shows a typical CNT-reinforced rectangular composite in which singlewalled CNTs (SWCNTs) are functionally distributed through the thickness within a polymer matrix. The carbon nanotubes are inserted such that their axes direct to the x-direction, and the composite has the length L, the depth D and the thickness h. The distribution of CNT through the thickness has a functional gradient as shown in Figure 1b, where four types of functionally graded CNTRCs are depicted: FG-U, FG-V, FG-O and FG-X. Since FG-CNTRCs are sort of dual-phase composite, their effective mechanical properties can be determined using the approximate method [9,28] such as the rule of mixture or the Mori-Tanaka estimate. For the current study, the modified linear rule of mixtures (LRM) is used by introducing the CNT efficiency parameters η j (j = 1, 2, 3). The volume fraction functions of CNTs in the above four FG-CNTRCs through the thickness are expressed by

Modeling of Functionally Graded CNT-Reinforced Composites
Here, the total volume fraction * cnt V of CNTs contained within the polymer matrix is calculated from the CNT mass fraction cnt w such that with cnt ρ and m ρ being the densities of CNTs and the polymer matrix.  According to the modified LRM, the effective elastic moduli and the effective shear modulus are calculated as [10] in which V cnt and V m = 1 − V cnt denote the volume fractions of CNTs and polymer which are in function of z. In Equations (1)-(3), the scripts cnt and m are used to denote the properties of SWCNT and polymer matrix. The polymer matrix is assumed to be homogeneous and isotropic while SWCTs are modeled as an orthotropic material, and, furthermore, it is assumed that E 3 = E 2 and G 23 = G 31 = G 12 . The scale dependence of the effective material properties of CNTRCs is reflected through the CNT efficiency parameters η j [29]. These parameters were obtained by matching the effective CNTRC properties predicted by the molecular dynamics (MD) simulation with those estimated by LRM. Table 1 presents the CNT efficiency parameters for a poly(methyl methacrylate) (PMMA) matrix with CNTs reinforcement at room temperature. The volume fraction functions of CNTs in the above four FG-CNTRCs through the thickness are expressed by Here, the total volume fraction V * cnt of CNTs contained within the polymer matrix is calculated from the CNT mass fraction w cnt such that with ρ cnt and ρ m being the densities of CNTs and the polymer matrix. Meanwhile, the effective Poisson's ratio ν 12 and the effective density ρ are determined in a similar manner Considering the FG-CNTRCs as a 3-D orthotropic body occupying a bounded domain Ω ∈ 3 , its displacement field u(x) = u(x, y, z) under the action of body force f and external loadt is governed by with the boundary conditions: Here, Γ D and Γ N denote the displacement and force boundaries, and σ ij and n j are the stress components and the unit normal vector, respectively.

Analysis of Bending Deformation
By virtue of the virtual work principle, the previous static equilibrium equations in Section 2 are converted to the following weak form for the elastic deformation field The components in [E] denote the effective orthotropic material properties of CNTRC structures, which will be given later. Using these two matrices, together with the FE approximation (12), both strains and stresses are approximated as and Introducing Equations (15) and (16) into the weak form (11), one can obtain the simultaneous linear equations to solve the static deformation problem: (17) where the stiffness matrix [K] and the load vector {f} are defined by

Optimization of CNT Distribution
For the present numerical optimization, an FG-CNTRC beam with the specific thicknesswise CNT gradient is divided into (ND) numbers of uniform homogenized sub-layers, as depicted in Figure 2. Then, the CNT volume fractions (V cnt ) i of each sub-layer define the design variable vector X: ..., 2 1 = X (20) From the physical constraint, each layer-wise CNT volume fraction ( ) i cnt V must satisfy the following upper and lower bounds ( ) and their sum should be equal to the preset volume fraction * cnt V of CNTs: Next, the peak effective stress is defined as the objective function ( ) with Ω being the entire material domain of FG-CNTRC beam. Then, the constrained optimization of thickness-wise CNT distribution is formulated as follows, together with the FE approximation (17): From the physical constraint, each layer-wise CNT volume fraction (V cnt ) i must satisfy the following upper and lower bounds and their sum should be equal to the preset volume fraction V * cnt of CNTs: Next, the peak effective stress is defined as the objective function F(X), such that with Ω being the entire material domain of FG-CNTRC beam. Then, the constrained optimization of thickness-wise CNT distribution is formulated as follows, together with the FE approximation (17): Note that Equation (27) is the equality constraint, while Equations (28) and (29) are the inequality constraints.
The constrained optimization problem is solved by employing the exterior penaltyfunction method (EPFM) [30]. This method converts the constrained objective function F(X) to an unconstrained pseudo-objective function Φ X; r p , by adopting the exterior penalty parameters r p , such that Meanwhile, the maximum value of objective function is in the order of 10 6 ∼ 10 7 while those of the constraints are 1.0, as will be shown later. Therefore, the normalization factors c 1 and c 2 are inserted to maintain the balance between the magnitudes of constraints and objective function. These two factors are calculated using the usual normalization given by The definition of inequality constraints given in Equations (28) and (29) leads to Plugging Equation (32) into Equation (31) results in The iterative optimization process starts with an initial design variable X 0 . At each iteration step k (k = 1, 2, . . .), the sensitivity analysis presented in Section 3.3 is performed and the convergence is examined as follows with ε T being the convergence tolerance. The iterative optimization is terminated when the convergence criterion is satisfied; otherwise, the optimization process proceeds to the next iteration by updating the current exterior penalty parameter Here, γ (γ > 1) indicates an iteration-independent update constant, which successively increases the penalty parameter during the iterative optimization. The flowchart of the proposed optimization process is represented in Figure 3.
The definition of inequality constraints given in Equations (28) and (29) leads to ( ) The iterative optimization process starts with an initial design variable 0 X . At each iteration step ( ) , the sensitivity analysis presented in Section 3.3 is performed and the convergence is examined as follows with T ε being the convergence tolerance. The iterative optimization is terminated when the convergence criterion is satisfied; otherwise, the optimization process proceeds to the next iteration by updating the current exterior penalty parameter indicates an iteration-independent update constant, which successively increases the penalty parameter during the iterative optimization. The flowchart of the proposed optimization process is represented in Figure 3.

Sensitivity Analysis
The sensitivity analysis calculates the direction vector S of the design variable, which is essential for searching the optimization direction. From Equations (23) and

Sensitivity Analysis
The sensitivity analysis calculates the direction vector S of the design variable, which is essential for searching the optimization direction. From Equations (23) and (27)-(29), the sensitivity of the pseudo-objective function Φ X; r p to the ith design variable X i is expressed by This direct method may be considered when the thickness-wise CNT distribution is assumed a priori such that the objective function F(X) is expressed in an explicit form of CNT volume fractions (V cnt ) i . However, it can be assumed that the first term on the right-hand side of Equation (36) requires a complex and painstaking analytical derivation. An effective and attractive alternative is to employ the finite difference method because the current optimization problem is not large-scale. When the finite difference is used, the direction vector S k = S k 1 , S k 2 , . . . , S k ND at the kth stage is computed by with r 0 p being the initial exterior penalty parameter. When the direction vector is obtained, the design variable is updated such that with the iteration-dependent constants β k for determining the direction vector magnitude. Two representative methods are widely used to decide the magnitude of direction vector: Lagrange interpolation and golden section methods [30]. Both methods commonly seek the critical value of β that minimizes the pseudo-objective function Φ X; c, r p , but the latter is preferable because it always secures the local minimum to the design variable vector. The detailed description for this method is skipped because the numerical implementation of this method is standardized.

Results and Discussion
Figure 4a represents the first model problem, a simpply supported FG-CNTRC beam subject to uniform distributed load q = 0.1 MPa. The length L is 0.1m and the depth D and the thickness h are equally set by 0.01m, respectively. The matrix is manufactured with poly(methyl methacrylate) (PMMA) and its isotropic material properties are given in Table 2, where the orthotropic material properties of the (10,10) SWCNTs are also presented. The FG-CNTRC beam is divided into 10 uniform sub-layers with layer-wise CNT volume fractions (V cnt ) i (i = 1, 2, · · · , 10), as represented in Figure 4b, in order to suppress the increase in the total number of design variables by considering the flexibility of thicknesswise CNT distribution at the same time. For the finite element structural analysis, each sub-layer is discretized into 100 × 10 uniform 8-node cubic elements, with 100 in the x-direction and 10 in the y-direction. Therefore, the total number of finite elements for the whole FG-CNTRC beam reaches 10,000.
This direct method may be considered when the thickness-wise CNT distribution is assumed a priori such that the objective function ( ) is expressed in an explicit form of CNT volume fractions ( ) i cnt V . However, it can be assumed that the first term on the right-hand side of Eq. (36) requires a complex and painstaking analytical derivation.
An effective and attractive alternative is to employ the finite difference method because the current optimization problem is not large-scale. When the finite difference is used, the direction vector with 0 p r being the initial exterior penalty parameter. When the direction vector is obtained, the design variable is updated such that with the iteration-dependent constants k β for determining the direction vector magnitude.
Two representative methods are widely used to decide the magnitude of direction vector: Lagrange interpolation and golden section methods [30]. Both methods commonly seek the critical value of β that minimizes the pseudo-objective function , but the latter is preferable because it always secures the local minimum to the design variable vector. The detailed description for this method is skipped because the numerical implementation of this method is standardized. , respectively. The matrix is manufactured with poly(methyl methacrylate) (PMMA) and its isotropic material properties are given in Table 2, where the orthotropic material properties of the (10,10) SWCNTs are also presented. The FG-CNTRC beam is divided into 10 uniform sub-layers with layer-wise CNT volume fractions ( ) ( )

Results and Discussion
, as represented in Figure 4b, in order to suppress the increase in the total number of design variables by considering the flexibility of thickness-wise CNT distribution at the same time. For the finite element structural analysis, each sub-layer is discretized into 10 100× uniform 8-node cubic elements, with 100 in the x -direction and 10 in the y -direction. Therefore, the total number of finite elements for the whole FG-CNTRC beam reaches 10,000.    The initial CNT volume fractions (V cnt ) i of ten sub-layers are set by the preset total CNT volume fraction V * cnt . For example, the initial layer-wise CNT volume fractions (V cnt ) i are equally designated by 0.12 when the preset value of V * cnt is 0.12. The convergence tolerance ε T in Equation (34) is set by 1.0 × 10 −3 and the initial penalty parameter r 0 p and the update constant γ in Equation (35) are taken by 1.0 and 2.0, respectively. Since this problem exhibits a remarkable edge effect in the bending stress field near the left and right ends of beam, the peak stress value is examined from the thickness-wise stress distribution at the beam mid-span (i.e., at x = L/2). The finite element analysis was carried out by midas NFX [31], a commercial FEM software.
For the preset CNT volume fraction V * cnt = 0.12, the optimization process terminates after five iterations, as shown in Table 3. The peak effective stress occurs at the top and bottom of beam, and 311 FEM analyses were performed, mostly for the sensitivity analysis. Table 3 shows that the objective function uniformly decreases proportionally to the iteration number. The peak effective stress σ max e f f is 7.504 MPa at the initial stage and 5.947 MPa at the final stage, so that it is reduced by 1.557 MPa, which corresponds to 20.7% of the initial peak effective stress. Figure 5a compares the initial and optimum CNT distributions, where the layer-wise CNT volume fractions (V cnt ) i are larger than 0 and less than 1.0 and their sum is found to be 11.95. Hence, it is confirmed that the equality constraint in Equation (27) and the inequality constraints in Equations (28) and (29) are strictly enforced. Figure 5b compares the thickness-wise effective stress distributions between the initial and optimum CNT distributions, where the initial one varies linearly from zero at the mid-surface to the peak value at the top and bottom. The optimum one is also zero at the mid-surface and symmetric with respect to the mid-surface, but is not linear any longer and its peak is smaller than that of the initial one. In the bending deformation, the bending strain exhibits a linear thickness-wise variation, which produces the linear bending stress distribution through the thickness when the elastic modulus is uniform. However, the parabolic-type optimum CNT distribution leads to the parabolic-type distribution of elastic modulus of FG-CNTRC beam, so that the non-linear stress distribution with the smaller peak effective stress is revealed.  The optimization is also performed for two different CNT volume fractio 17 0. V * cnt = and 0.28, and the optimization results are compared in Table 4. The total i ation numbers are similar, but the total numbers of FE analyses are different. It is p sumed that the sensitivity analysis is influenced by the value of * V . The peak effec The optimization is also performed for two different CNT volume fractions, V * cnt = 0.17 and 0.28, and the optimization results are compared in Table 4. The total iteration numbers are similar, but the total numbers of FE analyses are different. It is presumed that the sensitivity analysis is influenced by the value of V * cnt . The peak effective stresses occur at the top and bottom of beam for all the three cases, and the initial and optimum peak effective stresses are similar. Thus, it is found that the influence of the total CNT volume fraction V * cnt on the CNT distribution optimization for minimizing the peak effective stress is not significant.  Figure 6a comparatively represents the optimum CNT distributions for three different values of V * cnt . It can be seen that three optimum CNT distributions are commonly parabolic-type and show different layer-wise CNT volume fractions (V cnt ) i to fulfill the equality constraint imposed upon V * cnt . Figure 6b compares the effective stress distributions between three different values of V * cnt , where it is observed that the difference between three distributions is not significant. Thus, it is again confirmedthat the total CNT volume fraction V * cnt does not impose any significant influence on the stress results.  Table 5, where σ Δ indicates the difference in max eff σ between the optimum and FG CNT distributions a the values within parenthesis denote the relative differences with respect to the optim stress. It is found that max eff σ is lowest for the optimum case, while this is the highest FG-V. Thus, it has been justified that the present optimization method successfully se the optimum CNT distribution, which can provide the maximum effective stress, wh is smaller than those of standard FG CNT distributions.  Next, the thickness-wise stress distribution for the optimum CNT distribution when V * cnt is 0.12 is compared with those of four standard functionally graded CNT distributions, FG-U, FG-V, FG-O and FG-X. The results are presented in Table 5, where ∆σ max e f f indicates the difference in σ max e f f between the optimum and FG CNT distributions and the values within parenthesis denote the relative differences with respect to the optimum stress. It is found that σ max e f f is lowest for the optimum case, while this is the highest for FG-V. Thus, it has been justified that the present optimization method successfully seeks the optimum CNT distribution, which can provide the maximum effective stress, which is smaller than those of standard FG CNT distributions. This fact is well represented in Figure 7a, where FG-V shows an unsymmetric effective stress distribution with the highest σ max e f f at the top of the beam because its CNT distribution is unsymmetric, as shown in Figure 1b. Meanwhile, FG-O produces the peak effective stress at z = ±0.3 in accordance with its CNT distribution. Furthermore, its peak effective stress is higher than that of FG-U even though its CNT distribution is similar to the optimum CNT distribution. For the sake of comparison, the axial stress distributions are comparatively represented in Figure 7b, to examine whether there is a remarkable difference between the effective stress, which is calculated by all six stress components, and a single bending stress component. Since the axial stress is the dominant stress component in the current bending deformation problem, the axial stress distributions are observed to be similar to the effective stress distributions, except for the minus sign in the upper-half region.  , which is 23.7% of the initial peak effective stress.

Iteration
Objective  Next, the optimization was performed again for the previous simply supported beam by inclining the distributed load. Referring to Figure 4a, the vertical distributed load was inclined by α = 45 • with respect to the z-axis. However, except for this inclined distributed load, the beam geometry and the simulation parameters were kept the same as the previous case. The layer-wise CNT volume fractions (V cnt ) i were initially set by 0.12, and the numerical optimization terminates in five iterations, as presented in Table 6. When compared with the previous case, the total number of FEM analyses was smaller and the location of peak effective stress was fixed at the bottom. The peak effective stress σ max e f f ws 5.666 MPa at the initial stage and 4.321 MPa at the final stage, being reduced by 1.345 MPa, which is 23.7% of the initial peak effective stress.  Figure 8a compares the iteration histories of objective function between the previous and current cases, where both cases produce a difference in effective stress magnitude but commonly show rapid and stable convergence. Figure 8b comparatively represents the initial and optimum CNT distributions, where the arithmetic sum of layer-wise CNT volume fractions (V cnt ) i was found to be 12.08. Thus, the preset CNT volume fraction V * cnt = 0.12 is strictly satisfied. Meanwhile, the current optimum CNT distribution was seen to be asymmetric, differing from the previous one shown in Figure 5a. This is because the inclined distributed load produces the asymmetric effective stress distribution, as will be seen below. Thus, it has been found that the optimum thickness-wise CNT distribution is significantly influenced by the loading condition. Additionally, this asymmetric CNT distribution cannot be tailored, even though the conventional primitive CNT distribution patterns shown in Figure 1b are combined.   Figure 7b, the axial stress in the present case does not vanish at the mid-surface. In other words, the zero-line moves slightly upwards because the inclined load not only produces the bending deformation but also the axial extension. Thus, the initial axial stress distribution is not symmetric but asymmetric with respect to the mid-surface. The optimum axial stress distribution becomes more asymmetric because the optimum CNT distribution is additionally asymmetric, as shown in Figure 8b. Figure 9b compares the thickness-wise effective stress distributions between the initial and optimum CNT distributions, where both the initial and optimum CNT distributions show higher peak effective stresses at the bottom because the inclined distributed load produces additional tensile stress. When compared with the previous case shown in Figure 5b, the present inclined load leads to a remarkably different effective stress distribution. The effective stress does not vanish near the mid-surface because the shear stress caused by the inclined distributed load is not negligible. Hence, it is found that not only the optimum CNT distribution, but also the resulting stress distribution, is strongly influenced by the loading condition. Figure 9a comparatively represents the axial stress distributions of the initial and optimum CNT distributions. Differing from the previous case shown in Figure 7b, the axial stress in the present case does not vanish at the mid-surface. In other words, the zero-line moves slightly upwards because the inclined load not only produces the bending deformation but also the axial extension. Thus, the initial axial stress distribution is not symmetric but asymmetric with respect to the mid-surface. The optimum axial stress distribution becomes more asymmetric because the optimum CNT distribution is additionally asymmetric, as shown in Figure 8b. Figure 9b compares the thickness-wise effective stress distributions between the initial and optimum CNT distributions, where both the initial and optimum CNT distributions show higher peak effective stresses at the bottom because the inclined distributed load produces additional tensile stress. When compared with the previous case shown in Figure 5b, the present inclined load leads to a remarkably different effective stress distribution. The effective stress does not vanish near the mid-surface because the shear stress caused by the inclined distributed load is not negligible. Hence, it is found that not only the optimum CNT distribution, but also the resulting stress distribution, is strongly influenced by the loading condition.
is not symmetric but asymmetric with respect to the mid-surface. The optimum axial stress distribution becomes more asymmetric because the optimum CNT distribution is additionally asymmetric, as shown in Figure 8b. Figure 9b compares the thickness-wise effective stress distributions between the initial and optimum CNT distributions, where both the initial and optimum CNT distributions show higher peak effective stresses at the bottom because the inclined distributed load produces additional tensile stress. When compared with the previous case shown in Figure 5b, the present inclined load leads to a remarkably different effective stress distribution. The effective stress does not vanish near the mid-surface because the shear stress caused by the inclined distributed load is not negligible. Hence, it is found that not only the optimum CNT distribution, but also the resulting stress distribution, is strongly influenced by the loading condition.

Conclusions
A numerical optimization method was introduced to seek the optimum thickness-wise CNT distribution that minimizes the effective stress in FG-CNTRC beams. An FG-CNTRC beam with a continuous CNT distribution was divided into several sub-layers with uniform CNT distributions, and the optimum CNT volume fractions of each sub-layer were sought by the exterior penalty-function method. This approach not only reduced the total number of design variables but simplified the numerical formulation. The benchmark numerical experiments were conducted to illustrate and validate the proposed numerical method and investigate the optimization results. Through the numerical result, the following major observations can be drawn:

•
The proposed method successfully seeks the optimum thickness-wise CNT distribution, which minimizes the effective stress, with rapid and stable convergence.

•
For the vertical distributed load, the optimum CNT distribution is symmetric and parabolic, and the effective stress distribution is also symmetric.

•
The peak effective stress occurred at the top and bottom, and was reduced by 20.7% after the optimization.

•
The optimum CNT distribution is different from the four primitive CNT distributions and leads to the peak effective stress, which is reduced by at least 26.2% compared to those of four primitive patterns. • Both the optimum CNT distribution and the associated stress distributions are significantly influenced by the loading condition.

•
For the inclined distributed load, both distributions become un-symmetric, but the initial peak stress that occurred at the bottom is reduced by 23.7% after optimization.
The proposed optimization method can be practically used to tailor the thickness-wise CNT distribution and secure the structural safety of FG-CNTRC beams against the external load. The current work is limited to the minimization of effective stress, so the optimized CNT distribution pattern may lead to a reduced structural sfiffness and, consequently, a larger bending deflection. Hence, the trade-off between the effective stress and the bending deflection would be worthwhile, and is a deserving topic for future work.