Optimal Tailoring of CNT Distribution in Functionally Graded Porous CNTRC Beams

This paper is concerned with the multi-objective optimization of thickness-wise CNT distribution in functionally graded porous CNT-reinforced composite (FG-porous CNTRC) beams. The mechanical behaviors of FG-porous CNTRC structures are strongly influenced by the thickness-wise distributions of CNTs and porosity. Nevertheless, several linear functions were simply adopted to represent the thickness-wise CNT distribution without considering the porosity distribution, so these assumed linear primitive CNT distribution patterns are not sufficient to respond to arbitrary loading and boundary conditions. In this context, this study presents the multi-objective optimization of thickness-wise CNT distribution in FG-CNTRC porous beams to simultaneously minimize the peak effective stress and the peak deflection. The multi-objective function is defined by the larger value between two normalized quantities and the design variable vector is composed of the layer-wise CNT volume fractions. The constrained multi-objective optimization problem is formulated by making use of the exterior penalty-function method and the aspiration-level adjustment. The proposed optimization method is demonstrated through the numerical experiments, and the optimization solutions are investigated with respect to the porosity distribution and the combination of aspiration levels for two single-objective functions. It is found from the numerical results that the optimum CNT distribution is significantly affected by the porosity distribution. Furthermore, the proposed method can be successfully used to seek an optimum CNT distribution within FG-porous CNTRC structures which simultaneously enhances the multi-objective functions.


Introduction
Functionally graded carbon nanotube-reinforced composites have been spotlighted as a state-of-the-art composite due to their excellent mechanical properties [1,2]. These advanced composites were developed based on the excellence of carbon nanotubes (CNTs) and the notion of functionally graded materials (FGMs) [3]. The mechanical strength of polymer composites increases dramatically when only a small amount of CNTs are inserted [4], so that CNT-reinforced composites (CNTRCs) can provide higher bending stiffness and superior free vibration and buckling behaviors than the conventional polymer composites [5][6][7][8]. In mechanical applications, CNTRCs have been produced in the form of beams, plates, and shells which exhibit the thickness-wise variation in their mechanical behaviors. Hence, the uniform distribution of CNTs through the thickness may not be appropriate to respond to such thickness-wise variation. According to Seidel and Lagoudas [9] and Qian et al. [10], it has been reported that the enhancement of mechanical properties of CNTRCs is limited when CNTs are uniformly distributed through the thickness. To overcome this limitation, Shen [11] and Ke et al. [12] proposed the purposeful thickness-wise CNT distributions by utilizing the notion of FGM, which is characterized by the continuity and functional tailoring of thickness-wise material composition distribution [13]. a multi-objective optimization method for seeking a best suitable thickness-wise CNT distribution for FG-porous CNTRC structures.
As an extension of our previous work [33], the optimization problem is defined by the minimization of the peak deflection and the peak effective stress of FG-CNTRC beams at the same time by considering the porosity. The relatively larger values between these two peak quantities are defined by the multi-objective function, and the thickness-wise porosity distribution is modeled by adopting cosine functions. An FG-porous CNTRC beam is uniformly divided into a finite number of sub-layers with layer-wise uniform CNT volume fractions, in order to suppress the increase of design variable number and maintain the CNT distribution flexibility at the same time. The CNT volume fractions of sub-layers are defined as the design variables and subjected to the inequality and equality constraints. The sensitivity analysis of the multi-objective function is explicitly carried out by employing the central difference method (CDM), while the solution of the constrained multi-objective optimization problem is sought by making use of the exterior penaltyfunction method and the golden section method. The numerical experiments are carried out to demonstrate the proposed optimization method and to investigate the optimization results with respect to the porosity distribution pattern and the weighting factors for two single-objective functions. The numerical results inform that the optimum CNT distribution can be successfully sought by the proposed optimization method and the optimization results are remarkably affected by the porosity distribution and the weighting factor.

Modeling of FG-Porous CNTRC Plate
A CNT-reinforced composite plate is shown in Figure 1a, where length, depth and thickness of the plate are denoted by L, D, and h. Single-walled CNTs (SWCNTs) are aligned along the x−axis and uniformly distributed in the thickness direction. The thickness-wise distribution of CNTs may have a functional gradation as represented in Figure 1b, where four primitive CNT distribution patterns showing different functional gradations are denoted by FG-U, FG-V, FG-O, and FG-X, respectively. In FG-U, the CNT distribution is uniform through the thickness, while the other three have linear variations such that zero at the bottom and the peak at the top in FG-V, zero at the bottom and top and the peak at the mid-surface in FG-O, and zero at the mid-surface and the peak at the bottom and top in FG-X. Being considered as a sort of dual-phase composite material, the effective material properties of FG-CNTRC plates can be predicted using either the modified linear rule of mixture (LRM) or the Mori-Tanaka method [11,21] in which the CNT efficiency parameters η j (j = 1, 2, 3) are introduced. Usually, the modified LRM is widely adopted due to its accuracy and easy application.
The scale effect of CNTs on the effective material properties of CNTRCs is taken into consideration by the CNT efficiency parameters j  [34], which were determined by matching the effective CNTRC properties obtained by the molecular dynamics (MD) simulation with those predicted by the LRM.  The volume fraction distributions of CNTs in the above four FG-CNTRCs through the thickness are mathematically expressed by: The CNTRC plates are usually modeled as an orthotropic material, and their effective elastic and shear modui are estimated by [11]: according to the modified LRM. Here, V cnt and V m = 1 − V cnt indicate the thickness-wise volume fractions of CNTs and polymer. The material properties of the SWCNT and polymer matrix are labeled by the scripts cnt and m, and it is further assumed that E 3 = E 2 and G 23 = G 31 = G 12 . The scale effect of CNTs on the effective material properties of CNTRCs is taken into consideration by the CNT efficiency parameters η j [34], which were determined by matching the effective CNTRC properties obtained by the molecular dynamics (MD) simulation with those predicted by the LRM. Table 1 presents η j for the CNTRCs composed of Poly (methyl methacrylate) (PMMA) matrix and CNTs at room temperature. The volume fraction distributions of CNTs in the above four FG-CNTRCs through the thickness are mathematically expressed by: with V * cnt being the total volume fraction of CNTs contained within the polymer matrix. Meanwhile, the effective Poisson's ratios ν αβ and the effective density ρ of CNTRCs are determined in a similar manner: Figure 2a represents an FG-porous CNTRC plate in which the porosity density varies in the z−direction only. In the current study, three different porosity distributions are considered: center-biased, lower-biased, and upper-biased, as depicted in Figure 2b. These distributions are called sym, unsym-1, and unsym-2 in this paper, and those are mathematically repressed by [35]: by referring to Phani and Niyogi [36].
Viewing FG-porous CNTRC plates as a 3-D orthotropic body occupying a bounded domain under the action of body force f and external load tˆ is governed by the 3-D elasticity theory: (10) with the boundary conditions: Here,  Viewing FG-porous CNTRC plates 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 the 3-D elasticity theory: with the boundary conditions: Here, Γ D and Γ N stand for the displacement and force boundary regions, and σ ij and n j indicate the stress components and the unit normal vector, respectively.

Analysis of Bending Deformation
The previous static equilibrium in Equation (10) is converted to the following variational form for solving the bending deformation field: for every admissible displacement v. Using iso-parametric finite elements, both the trial and test displacement fields u and v are approximated as: where [Φ] is a (3 × 3N) matrix expressed in terms of FE basis functions and u and v denote the (3N × 1) vectors of nodal displacements.
where C ij indicate the orthotropic material constants [37] expressed in terms of the effective elastic moduli E i in Equation (1) and the effective Poisson's ratios ν ij in Equation (4). Then, both the strain and stress components are approximated as: and Plugging Equations (17) and (18) into the variational form (13), one can derive the simultaneous equation system for solving the plate bending deformation: where the stiffness matrix [K] and the load vector {F} are calculated as:

Multi-Objective Optimization of CNT Distribution
For the effective numerical optimization with a reasonable number of design variables, an FG-porous CNTRC beam is divided into (NDV) numbers of uniform homogenized sub-layers, as shown in Figure 3. Then, the CNT volume fractions (V cnt ) I of each sub-layer constitutes the design variable vector X such that: ..., 2 1 = X (22) Owing to the physical constraint, each sub-layer-wise CNT volume fraction ( ) with  being the entire material domain of the FG-porous CNTRC beam. Then, the multi-objective (MO) function ( ) (26) in terms of the ideal levels ( ) 2 1, and the aspiration levels I f of two so functions. The ideal level is the optimum solution obtained by the SO optimization, and the aspiration level indicates the level to be reached by the MO optimization. In the current study, the aspiration level is set by: once both the ideal and aspiration levels are given.
For a detailed explanation, refer to reference [38]. Then, the constrained MO optimization problem of thickness-wise CNT distribution is formulated as follows: Owing to the physical constraint, each sub-layer-wise CNT volume fraction (V cnt ) I should obey the following lower and upper bounds: and their arithmetic average should be equal to the preset volume fraction V * cnt of CNTs: Letting the vertical displacement u 3 be w, two single-objective (SO) functions f 1 (x, X) are f 2 (x, X) are defined by: with Ω being the entire material domain of the FG-porous CNTRC beam. Then, the multi-objective (MO) function F(X) defined by: in terms of the ideal levels f * I (I = 1, 2) and the aspiration levelsf I of two so functions. The ideal level is the optimum solution obtained by the SO optimization, and the aspiration level indicates the level to be reached by the MO optimization. In the current study, the aspiration level is set by:f with a 0 = f 0 I / f * I , where f 0 I denote the values obtained by the initial design variable vector X 0 . Note that the relative weights w f I of each SO function is automatically determined such that w f I = 1/ f I − f * I once both the ideal and aspiration levels are given. For a detailed explanation, refer to reference [38].
Then, the constrained MO optimization problem of thickness-wise CNT distribution is formulated as follows: It is worth noting that Equation (31) indicates the equality constraint in Equation (24), while Equations (32) and (33) present the inequality constraints in Equation (23), respectively. The above weighted constrained optimization problem is solved by the exterior penalty-function method (EPFM) [39], which converts the objective function F(X) subject to the constraints to an unconstrained pseudo-objective function Φ X; r p given by: by introducing the exterior penalty parameters r p . Here, c 1 and c 2 are the normalization factors to keep the balance between the magnitudes of the objective function F(X) and the constraints, which are calculated according to: Here, the inequality constraints given in Equations (32) and (33) leads to: So, Equation (33) ends up with: The optimization iteration starts with a guessed initial design variable X 0 . At each iteration step, the sensitivity analysis is performed, and the convergence is checked using the convergence criterion defined by: The iteration is terminated when the convergence tolerance ε T is satisfied; otherwise, it moves to the next iteration by updating the exterior penalty parameter: where an iteration-independent update constant γ (γ > 1) uniformly increases the penalty parameter r p along the optimization iteration. Figure 4 represents the flowchart of the present optimization process.

Sensitivity Analysis
The searching of optimization direction is essential in the mathematical optimization, and it is accomplished by computing the direction vector S of the design variables through the sensitivity analysis. The sensitivity of the pseudo-objective function Φ X; c, r p with respect to the I − th design variable X I is mathematically expressed by: according to Equations (25) and (31)- (33). This direct mathematical derivation may be adopted when the thickness-wise CNT distribution is assumed such that the objective function F(X) can be explicitly expressed in terms of design variables (V cnt ) I . However, even though it can be assumed, the mathematical derivation of the first term on the righthand side of Equation (40) is highly painstaking.

Sensitivity Analysis
The searching of optimization direction is essential in the mathematical optimiza- according to Equations (25) and (31)- (33). This direct mathematical derivation may be adopted when the thickness-wise CNT distribution is assumed such that the objective function ( ) X F can be explicitly expressed in terms of design variables ( ) . However, even though it can be assumed, the mathematical derivation of the first term on the right-hand side of Equation (40) is highly painstaking.
This direct method can be effectively replaced with the finite difference method (FDM) when the problem size is not so large. In the finite difference method, the direction vector at the th k iteration is calculated as: with the initial exterior penalty parameter 0 p r . Once the direction vector is calculated, the design variable vector is updated according to: where k  is the iteration-dependent parameter to determine the magnitude of k S , and it is calculated by the golden section method [39], which always guarantees the local minimum with respect to the design variable vector X . This direct method can be effectively replaced with the finite difference method (FDM) when the problem size is not so large. In the finite difference method, the direction vector at the k th iteration is calculated as: with the initial exterior penalty parameter r 0 p . Once the direction vector is calculated, the design variable vector is updated according to: where β k is the iteration-dependent parameter to determine the magnitude of S k , and it is calculated by the golden section method [39], which always guarantees the local minimum with respect to the design variable vector X.

Results and Discussion
Figure 5a shows a simply supported FG-porous CNTRC beam under a uniform vertical distributed load q = 0.1 MPa. The beam length L is 0.1 m and the width a and the thickness h are equally 0.01 m, respectively. The porosity parameter e of three cosine-type porosity distributions shown in Figure 2b is set by 0.2. The beam is manufactured with isotropic matrix of Poly (methyl methacrylate) (PMMA) and orthotropic (10,10) single-walled CNTs (SWCNTs). The material properties of constituent materials are given in Table 2, where 1, 2, and 3 indicate x, y and z, respectively. Figure 5b represents ten uniform sub-layers which are divided to suppress the increase of total design variable number by maintaining the flexibility of thickness-wise CNT distribution at the same time. Here, V I cnt indicates the volume fraction of the I−th sub-layer so that the total number of design variables becomes ten. Each sub-layer is uniformly discretized using 100 × 10 4-node cubic finite elements such that 100 along the x-axis and 10 along the y-axis. So, the whole FG-porous CNTRC beam is uniformly discretized by 10,000 finite elements, and the FE static analysis was carried out by commercial FEM code midas NFX [40]. same time. Here, I cnt V indicates the volume fraction of the − I th sub-layer so that the total number of design variables becomes ten. Each sub-layer is uniformly discretized using 10 100  4-node cubic finite elements such that 100 along the x -axis and 10 along the y -axis. So, the whole FG-porous CNTRC beam is uniformly discretized by 10,000 finite elements, and the FE static analysis was carried out by commercial FEM code midas NFX [40].  First, the single-objective optimization for minimizing the peak effective stress is performed with the symmetric porous distribution, for which the initial CNT distribution pattern and the initial CNT volume fractions  Figure 5a exhibits the edge effect in the stress field near the left and right ends of the beam. Thus, the peak value of stress is extracted from the thickness-wise stress distribution at the mid-span of the beam. The optimization process is terminated in five iterations, as given in Table 3, where the objective function shows a uniform convergence. The peak effective stress occurs simultaneously at the beam top and bottom, and a total of 232 FEM analyses were performed, mostly for the sensitivity analyses. The initial peak effective stress   First, the single-objective optimization for minimizing the peak effective stress is performed with the symmetric porous distribution, for which the initial CNT distribution pattern and the initial CNT volume fractions V I cnt of ten sub-layers are set by FG-U and 0.12. The values of simulation parameters ε T , r 0 p and γ are set by 1.0 × 10 −3 , 1.0 and 2.0, respectively. The beam bending shown in Figure 5a exhibits the edge effect in the stress field near the left and right ends of the beam. Thus, the peak value of stress is extracted from the thickness-wise stress distribution at the mid-span of the beam. The optimization process is terminated in five iterations, as given in Table 3, where the objective function shows a uniform convergence. The peak effective stress occurs simultaneously at the beam top and bottom, and a total of 232 FEM analyses were performed, mostly for the sensitivity analyses. The initial peak effective stress 7.859 MPa is reduced to 5.985 MPa at the final stage, and the reduction amount is 1.874 MPa, which corresponds to 23.8% of the initial peak effective stress.  Figure 6a comparatively represents the initial and optimum CNT distributions, where the optimum CNT distribution for a non-porous (i.e., dense) CNTRC beam is added for comparison purposes. The arithmetic average of layer-wise CNT volume fractions V I cnt is found to be 0.1208 so that the equality constraint in Equation (31) is strictly satisfied. One can see the difference in CNT distributions between porous and non-porous cases, and the difference becomes more apparent at the central region. This is because the porosity in the symmetric porosity distribution is dominated in the central region, as depicted in Figure 2b. Figure 6b comparatively represents the thickness-wise effective stress distributions of initial and optimum CNT distributions. The initial one is symmetric and shows a linear variation from zero at the mid-surface to the peak value at the beam bottom and top. In the bending of a homogeneous beam, the linearly varying bending strain produces the linear bending stress distribution through the thickness. Meanwhile, the optimal CNT distribution also leads to a symmetric effective stress distribution with zero at the midsurface, but its distribution is not linear anymore and its peak is smaller than that of the initial CNT distribution. This is because the parabolic-type distribution of CNTs and porosity leads to the parabolic distribution of elastic modulus, which leads to the non-linear stress distribution with smaller peak effective stress. The difference in the effective stress distribution between porous and non-porous cases is not so remarkable. satisfied. One can see the difference in CNT distributions between porous and non-porous cases, and the difference becomes more apparent at the central region. This is because the porosity in the symmetric porosity distribution is dominated in the central region, as depicted in Figure 2b. Figure 6b comparatively represents the thickness-wise effective stress distributions of initial and optimum CNT distributions. The initial one is symmetric and shows a linear variation from zero at the mid-surface to the peak value at the beam bottom and top. In the bending of a homogeneous beam, the linearly varying bending strain produces the linear bending stress distribution through the thickness. Meanwhile, the optimal CNT distribution also leads to a symmetric effective stress distribution with zero at the mid-surface, but its distribution is not linear anymore and its peak is smaller than that of the initial CNT distribution. This is because the parabolic-type distribution of CNTs and porosity leads to the parabolic distribution of elastic modulus, which leads to the non-linear stress distribution with smaller peak effective stress. The difference in the effective stress distribution between porous and non-porous cases is not so remarkable. The numerical optimization was also performed for two different porosity distributions, unsym-1 and unsym-2, and the optimization results are summarized in Table 4. The total numbers of iteration and FEM analyses of unsym-1 and 2 are smaller than those The numerical optimization was also performed for two different porosity distributions, unsym-1 and unsym-2, and the optimization results are summarized in Table 4. The total numbers of iteration and FEM analyses of unsym-1 and 2 are smaller than those of sym, but unsym-1 and 2 show the same total numbers in iteration and FEM analyses. The peak effective stress occurs at the top of the beam for unsym-1 while at the bottom for unsym-2, but its magnitude is shown to be the same for both cases. When compared with the case of sym, unsym-1 and 2 lead to smaller initial and final effective stresses. Thus, it is found that there exists a remarkable difference in the CNT distribution between symmetric and unsymmetric porous distributions, but two unsymmetric porous distributions show the difference only in the location of peak effective stress.  Figure 7a,b comparatively represent the optimum CNT distributions between unsym-1 and unsym-2 porosity distributions. It is observed that both distributions of unsym-1 and 2 are exactly anti-symmetric with respect to the mid-surface. This anti-symmetry is solely owing to the anti-symmetry in their porosity distributions. When compared with the optimum CNT distribution of sym, unsym-1 shows a slightly larger difference in the lower beam region while unsym-2 leads a slightly larger difference in the upper beam region. This is because the porosity in unsym-1 is dominated in the lower region but that in unsym-2 is dominated in the upper region, as depicted in Figure 2b. of sym, but unsym-1 and 2 show the same total numbers in iteration and FEM analyses. The peak effective stress occurs at the top of the beam for unsym-1 while at the bottom for unsym-2, but its magnitude is shown to be the same for both cases. When compared with the case of sym, unsym-1 and 2 lead to smaller initial and final effective stresses. Thus, it is found that there exists a remarkable difference in the CNT distribution between symmetric and unsymmetric porous distributions, but two unsymmetric porous distributions show the difference only in the location of peak effective stress.  Figure 7a,b comparatively represent the optimum CNT distributions between unsym-1 and unsym-2 porosity distributions. It is observed that both distributions of unsym-1 and 2 are exactly anti-symmetric with respect to the mid-surface. This anti-symmetry is solely owing to the anti-symmetry in their porosity distributions. When compared with the optimum CNT distribution of sym, unsym-1 shows a slightly larger difference in the lower beam region while unsym-2 leads a slightly larger difference in the upper beam region. This is because the porosity in unsym-1 is dominated in the lower region but that in unsym-2 is dominated in the upper region, as depicted in Figure 2b.  Figure 8a compares the effective stress distributions of three porosity distributions for the optimal CNT distribution. The zero line in the effective stress distribution of unsym-1 moves slightly upwards, and vice versa for unsym-2. This is because the porosity dominance in the lower beam region for unsym-1 gives rise to higher elastic modulus in the upper beam region, and vice versa for unsym-2. However, it is seen that the relative difference in the magnitudes of effective stress is not remarkable. Figure 8b represents the dependence of iteration history of objective function on the porosity distribution. Un-  Figure 8a compares the effective stress distributions of three porosity distributions for the optimal CNT distribution. The zero line in the effective stress distribution of unsym-1 moves slightly upwards, and vice versa for unsym-2. This is because the porosity dominance in the lower beam region for unsym-1 gives rise to higher elastic modulus in the upper beam region, and vice versa for unsym-2. However, it is seen that the relative difference in the magnitudes of effective stress is not remarkable. Figure 8b represents the dependence of iteration history of objective function on the porosity distribution. unsym-1 and 2 show exactly the same iteration history because both porosity distributions produce exactly anti-symmetric stress distributions through the thickness. Except for the relative , respectively. As will be presented later in detail, the SO optimization for minimizing the peak deflection was terminated in five iterations, and the initial peak deflection of mm .05229 0 was reduced to mm .03174 0 . Table 5 presents the variation of MO function ( ) to the iteration, where the objective function decreases with small fluctuation and the optimization terminates in six iterations with a total of 271 FEM analyses. The peak effective stress occurs initially at the beam top and bottom, but the peak location moves to the inside of the beam during the next three iterations, and thereafter the peak effective stress occurs again at the beam top and bottom. The optimal CNT distribution is represented in Figure 9a, where those obtained by two SO optimizations are also given for comparison purposes. Here, denote the SO optimizations for the peak effective stress and for the peak deflection, respectively. It is clearly observed that SO . This result was obtained because MO tries to minimize the peak effective stress and the peak deflection at the same time. This trend of MO is also observed from Figure 9b, where the thickness-wise effective stress distribu- Next, the multi-objective optimization for minimizing the peak effective stress and the peak deflection was performed for the sym porosity distribution using the same FGporous CNTRC beam shown in Figure 5a. The aspiration levels in Equation (27) for two SO functions f 1 and f 2 were set by s 1 = s 2 = 1.2, while the other simulation parameters were kept the same with the previous single-objective optimization. For the sym porosity distribution, the ideal levels f * 1 and f * 2 were 5.98457 MPa and 0.03174 mm which were obtained by SO optimization for the peak effective stress and for the peak deflection, respectively. Thus, the values taken for two aspiration levels becomef 1 = s 1 f * 1 = 7.18148 MPa andf 2 = s 2 f * 2 = 0.03809 mm, respectively. As will be presented later in detail, the SO optimization for minimizing the peak deflection was terminated in five iterations, and the initial peak deflection of 0.05229 mm was reduced to 0.03174 mm. Table 5 presents the variation of MO function F(X) to the iteration, where the objective function decreases with small fluctuation and the optimization terminates in six iterations with a total of 271 FEM analyses. The peak effective stress occurs initially at the beam top and bottom, but the peak location moves to the inside of the beam during the next three iterations, and thereafter the peak effective stress occurs again at the beam top and bottom. The optimal CNT distribution is represented in Figure 9a, where those obtained by two SO optimizations are also given for comparison purposes. Here, SO −σ and SO −w denote the SO optimizations for the peak effective stress and for the peak deflection, respectively. It is clearly observed that SO −σ and SO −w lead to the optimal CNT distributions showing the opposite thickness-wise variations. The former leads to smaller CNT volume fraction at the beam top and bottom to decrease the peak effective stress by decreasing the elastic modulus, while the latter leads to larger CNT volume fraction at the beam top and bottom to decrease the peak deflection by increasing the beam bending stiffness. Meanwhile, MO leads to the optimal CNT distribution which is placed between those of SO −σ and SO −w . This result was obtained because MO tries to minimize the peak effective stress and the peak deflection at the same time. This trend of MO is also observed from Figure 9b, where the thickness-wise effective stress distribution of MO is placed between those of SO −σ and SO −w . Note that the order in the magnitude of peak effective stress is SO −w > MO > SO −σ , which is consistent with the minimization goals of SO −w , MO, and SO −σ . Table 5. Variation of MO function to the iteration (sym, s 1 = s 2 = 1.2 ).    Table 6 and represented in Figure 10. First of all, it is found that max eff  and max w obtained by the MO optimization are larger than those obtained by the SO optimizations. This is consistent with the fact that the optimum solution of MO optimization is always larger than the optimum solution of SO optimization. The total number of iterations of MO is longer than one of SO because the MO optimization exhibits larger fluctuation than the SO optimization, as can be realized from the comparison of Tables 3 and 5.  The iteration histories of σ max e f f and w max between MO and SO optimizations are compared in Table 6 and represented in Figure 10. First of all, it is found that σ max e f f and w max obtained by the MO optimization are larger than those obtained by the SO optimizations. This is consistent with the fact that the optimum solution of MO optimization is always larger than the optimum solution of SO optimization. The total number of iterations of MO is longer than one of SO because the MO optimization exhibits larger fluctuation than the SO optimization, as can be realized from the comparison of Tables 3 and 5.   Table 7, where the results of the previous case (1.2:1.2) are also given for the purpose of comparison. The case of (1.3:1.1) shows the total numbers of iteration and FEM analyses which are similar to those of the case of (1.2:1.2), but the case of (1.1:1.3) leads to smaller total numbers in the optimization iteration and FEM analyses. Meanwhile, it is found that the case of (1.3:1.1) leads to the optimum solution equal to F(X) = 0.66486, which is smaller than those of (1.2:1.2) and (1.1:1.3). This is because the case of (1.3:1.1) places more weight on the reduction of peak deflection, and the optimum value of F(X) was determined from the relative reduction in the peak deflection, as given in Equation (26). This can be confirmed from the fact that the case of (1.3:1.1) leads to smaller peak deflection w max = 0.03248 mm but larger peak effective stress σ max e f f = 7.17380 MPa compared with the previous case of (1.2:1.2). On the other hand, the case of (1.1:1.3) places more weight on the reduction in the peak effective stress, so that it leads to smaller peak effective stress σ max e f f = 6.08067 MPa but larger peak deflection w max = 0.04219 mm compared with the case of (1.2:1.2). Thus, it has been justified that the trade-off between the peak effective stress and the peak deflection is successfully accomplished by adjusting the aspiration levels (i.e., the combination of (s 1 : s 2 )).  Table 7, where the results of the previous case (1.2:1.2) are also given for the purpose of comparison. The case of (1.3:1.1) shows the total numbers of iteration and FEM analyses which are similar to those of the case of (1.2:1.2), but the case of (1.1:1.3) leads to smaller total numbers in the optimization iteration and FEM analyses. Meanwhile, it is found that the case of (1.3:1.1) leads to the optimum solution equal to ( )

Multi-Objective Function
, which is smaller than those of (1.2:1.2) and (1.1:1.3). This is because the case of (1.3:1.1) places more weight on the reduction of peak deflection, and the optimum value of ( ) X F was determined from the relative reduction in the peak deflection, as given in Equation (26). This can be confirmed from the fact that the case of (1 compared with the case of (1.2:1.2). Thus, it has been justified that the trade-off between the peak effective stress and the peak deflection is successfully accomplished by adjusting the aspiration levels (i.e., the combination of (

Items
Combinations of Aspiration Levels (   Figure 11a comparatively represents the optimal CNT distributions for three different combinations of aspiration levels. When compared with the case of (1.2:1.2), the case of (1.1:1.3) leads to the distribution in which the CNT volume fraction V cnt is high in the central region while low in the vicinity of top and bottom. This is because the peak effective stress occurring at the beam top and bottom can be effectively reduced by reducing the elastic modulus of the beam in the vicinity of the beam top and bottom, and the reduction of elastic modulus can be made by decreasing the CNT volume fraction. Meanwhile, the case of (1.3:1.1) leads to the optimum CNT distribution which is almost opposite to that of the case of (1.1:1.3). This is because the beam deflection can be reduced by increasing the bending flexural rigidity, and the flexural rigidity can be effectively increased by placing more CNTs in the vicinity of the beam top and bottom. Figure 11b compares the thicknesswise effective stress distributions, where the case of (1.3:1.1) shows the highest level while the case of (1.1:1.3) leads to the lowest level. Meanwhile, the effective stress distribution of the case of (1.2:1.2) passes through between those of the cases of (1.1:1.3) and (1.3:1.1). This comparison justified again the trade-off between σ max e f f and w max according to the adjustment of aspiration levels.
the reduction of elastic modulus can be made by decreasing the CNT volume fraction. Meanwhile, the case of (1.3:1.1) leads to the optimum CNT distribution which is almost opposite to that of the case of (1.1:1.3). This is because the beam deflection can be reduced by increasing the bending flexural rigidity, and the flexural rigidity can be effectively increased by placing more CNTs in the vicinity of the beam top and bottom. Figure 11b compares the thickness-wise effective stress distributions, where the case of (1.3:1.1) shows the highest level while the case of (1.1:1.3) leads to the lowest level. Meanwhile, the effective stress distribution of the case of (1.2:1.2) passes through between those of the cases of (1.

Conclusions
The multi-objective optimization of FG-CNTRC porous beams was presented to simultaneously minimize the peak effective stress and the peak deflection. An FG-CNTRC porous beam was divided into uniform sub-layers, and the layer-wise uniform CNT volume fractions were chosen as the design variables. The peak effective stress and the peak deflection were normalized using their ideal and aspiration levels which were obtained by the single-objective optimization and set by the designer's decision making. The maximum value between two normalized quantities was defined as the multi-objective function, and it was minimized by making use of the exterior penalty-function method and the trade-off process. The proposed method was demonstrated and verified through the benchmark experiment, and the optimization solutions were investigated with respect to the porosity distribution pattern and the combination of aspiration levels. The numerical results provide the following major observations: Figure 11. Trade-off between σ max e f f and w max : (a) optimal CNT distribution, (b) thickness-wise effective stress distribution.

Conclusions
The multi-objective optimization of FG-CNTRC porous beams was presented to simultaneously minimize the peak effective stress and the peak deflection. An FG-CNTRC porous beam was divided into uniform sub-layers, and the layer-wise uniform CNT volume fractions were chosen as the design variables. The peak effective stress and the peak deflection were normalized using their ideal and aspiration levels which were obtained by the single-objective optimization and set by the designer's decision making. The maximum value between two normalized quantities was defined as the multi-objective function, and it was minimized by making use of the exterior penalty-function method and the trade-off process. The proposed method was demonstrated and verified through the benchmark experiment, and the optimization solutions were investigated with respect to the porosity distribution pattern and the combination of aspiration levels. The numerical results provide the following major observations:

•
The proposed MO optimization method successfully seeks the optimum CNT distribution, which simultaneously reduces the peak effective stress and the peak deflection with the stable convergence.

•
The optimal CNT distribution of a porous beam is different from that of a non-porous beam, and the difference is apparent in the region where the porosity is dominant. However, the effect of porosity on the optimum effective stress distribution is not so remarkable.

•
The porosity distribution pattern significantly affects the optimum CNT distribution but its effect on the optimum effective stress distribution is insignificant. unsym-1 and 2 porosity distributions lead to the optimum CNT distributions, which are antisymmetric to each other.

•
The MO optimization leads to the optimal CNT distribution and the optimal effective stress distribution, which are placed between those of SO −σ and SO −w . The order in the magnitude of σ max e f f is SO −w > MO > SO −σ , and that of w max is SO −σ > MO > SO −w .

•
The trade-off between the peak effective stress and the peak deflection can be successfully accomplished according to the adjustment of aspiration levels such that the decrease of s I places more weight on the corresponding SO function f I .

•
Regarding the trade-off between two SO functions, the case of (1.1:1.3) leads to the lowest peak effective stress while the case of (1.3:1.1) provides the smallest peak deflection.