Maximizing natural frequencies of inhomogeneous cellular structures by Kriging-assisted multiscale topology optimization

Computers and Structures 230 (2020) 106197

Contents lists available at ScienceDirect

Computers and Structures journal homepage: www.elsevier.com/locate/compstruc

Maximizing natural frequencies of inhomogeneous cellular structures by Kriging-assisted multiscale topology optimization Yan Zhang, Liang Gao, Mi Xiao ⇑ State Key Laboratory of Digital Manufacturing Equipment and Technology, Huazhong University of Science and Technology, 1037 Luoyu Road, Wuhan, Hubei 430074, China

a r t i c l e

i n f o

Article history: Received 17 July 2019 Accepted 20 December 2019

Keywords: Natural frequency Multiscale topology optimization Cellular structure Inhomogeneous microstructures Kriging metamodel Shape interpolation

a b s t r a c t This paper proposes a Kriging-assisted multiscale topology optimization method for maximizing natural frequencies of inhomogeneous cellular structures, where spatially-varying microstructural configurations and their macroscopic distribution are simultaneously optimized. At the beginning, under macroscopic boundary conditions for a cellular structure, the configurations of multiple prototype microstructures are topologically optimized by the parametric level set method (PLSM) combined with the numerical homogenization approach. A kinematical connective constraint is considered to ensure the connectivity between adjacent prototype microstructures. Then, a shape interpolation method is adopted to interpolate shapes of the prototype microstructures, so as to generate a series of sample microstructures. Based on these samples, Kriging metamodels are constructed to predict the effective property of each microstructure within the macrostructure. Finally, the variable thickness sheet (VTS) method is applied to optimize the material distribution pattern at macroscale for maximizing the natural frequency of the cellular structure, where an efficient mode-tracking strategy based on modal assurance criterion (MAC) is employed to track the target mode accurately. Numerical examples are provided to test the performance of the proposed method in natural frequency optimization of cellular structures. The results indicate that the multiscale cellular structures obtained by the proposed method show higher natural frequency compared with the monoscale macrostructural and microstructural designs. Ó 2019 Elsevier Ltd. All rights reserved.

1. Introduction Dynamic optimization problems are of major importance to suppress the undesirable vibration of structures with the demand for structural stability, durability and noise reduction. Amongst various approaches developed to achieve this goal, structural topology optimization has been extensively investigated in the last three decades. Essentially, topology optimization is a numerical iterative procedure that aims to find the optimal material layout inside a fixed design domain to maximize structural performance while satisfying a set of constraints, such as a given amount of material. Starting from the pioneering work of employing the homogenization theory by Bendsøe and Kikuchi [1], various popular methods have been invented and applied to topology optimization, such as solid isotropic material with penalization (SIMP) [2,3], bi-directional evolutionary structural optimization (BESO) [4,5], and level set method (LSM) [6,7].

⇑ Corresponding author. E-mail address: [email protected] (M. Xiao). https://doi.org/10.1016/j.compstruc.2019.106197 0045-7949/Ó 2019 Elsevier Ltd. All rights reserved.

Currently, topology optimization is widely employed at macroscale to improve the macrostructural performance. Meanwhile, topology optimization methods have broaden their applications at microscale to design new material microstructures with prescribed or extreme properties, where the effective physical properties of the material microstructures are calculated by a numerical homogenization method [8,9]. For instance, Sigmund [10] formulated an inverse homogenization problem to find the simplest microstructure with the prescribed elastic properties. Gao et al. [11], Da et al. [12,13], and Zhang et al. [14] designed periodic microstructures of cellular materials with the maximum bulk or shear modulus by various topology optimization methods. Wang et al. [15] applied the topology optimization method in material design to achieve a type of mechanical metamaterials with negative Poisson’s ratios (auxetic materials). Recently, the optimal material microstructures with functionally graded properties [16] and multifunctional properties [17] have been investigated by topology optimization methods. The aforementioned researches on topology optimization were restricted at a single scale, namely either at macroscale or at microscale. To fully explore the design space and further improve the

2

Y. Zhang et al. / Computers and Structures 230 (2020) 106197

macrostructural performance, multiscale topology optimization of macrostructures has been studied. Generally, multiscale topology optimization aims to simultaneously optimize the topology of a macrostructure and the geometrical configurations of material microstructures. In this sense, current multiscale topology optimization methods can be roughly categorized into three branches. The first is under the assumption that material microstructures are identical throughout the entire material region of the macrostructure [18–23]. The second has been studied under another assumption that the macrostructure comprises multiple types of material microstructures [24–28]. Some heuristic criteria for distributing multiple material microstructures in the entire macrostructure are necessary to implement multiscale topology optimization [27], such as trial-and-error criterion to divide the geometrical domain [25] and the principal stress distribution [24]. The third is a general multiscale design strategy inspired by the pioneering work of Rodrigues et al. [29], where material microstructures vary from point to point throughout the entire macrostructure [30–33]. This kind of multiscale design strategy can provide the sufficient design freedom at both scales for pursuing the superior macrostructural performance. However, there are two main obstacles in academic and industrial applications of this multiscale design strategy, i.e., the prohibitive computational burdens and the connectivity between two adjacent microstructures. Recently, a variable-density lattice design method for additive manufactured cellular structures is innovatively proposed [34–36] to achieve spatially-varying mechanical properties, where the explicit cellular structure is finally reconstructed by mapping the optimized continuous parameters (e.g., density) to cell structural parameters (e.g., strut diameter). Subsequently, the work is extended to cooling channel system design [37,38]. Frequency optimization is one of the most important topics in topology optimization, which is desirable to shift the natural frequency of macrostructures in order to avoid destructive responses caused by external excitations [39–42]. At macroscale, various material interpolation schemes were proposed to improve the classical SIMP model, which is unsuitable for frequency optimization due to artificial localized modes [41–43]. At microscale, Huang et al. [44] extended the bi-directional evolutionary structural optimization (BESO) method to design the composite microstructure with the optimal viscoelastic characteristics for the natural frequency of a macrostructure. Youn et al. [45] designed the material microstructures with high damp properties to reduce unwanted vibrations under dynamic loads using density-based topology optimization. Andreassen and Jensen [46] presented a topology optimization to design periodic composites with enhanced dynamic properties for maximizing loss/attenuation of propagating waves. To sum up, the aforementioned studies are confined to monoscale structural design for frequency optimization. On the other hand, at multiscale, Niu et al. [47] presented a two-scale topology optimization method for maximum structural fundamental frequency using a porous anisotropic material with penalization (PAMP) model. Zuo et al. [48] introduced a hierarchical concurrent design method to maximize the natural frequency of a structure. Based on the similar work, Liu et al. [49] realized the automatic allocation of material usage at both macro and micro scales to develop a genuine concurrent topology optimization algorithm for maximizing natural frequency. However, these aforementioned works at multiscale is under the assumption that the macrostructure is composed of identical material microstructure. Very little attention has been devoted to the multiscale design of macrostructures with spatially-varying connectable inhomogeneous microstructures for maximizing the natural frequency.

This paper concentrates on maximizing the natural frequency of cellular structures with spatially-varying inhomogeneous microstructures by a Kriging-assisted multiscale topology optimization method, which can fully explore the design space at both scales with an affordable computation burden. Firstly, at microscale, the optimal design of multiple prototype microstructures is gained by the parametric level set method (PLSM) integrated with the numerical homogenization approach. The connectivity between microstructures is ensured by a kinematical connective constraint. Then, a shape interpolation method is employed to interpolate level set functions of the optimized prototype microstructures and generate a family of key inhomogeneous microstructures, which are considered as sample points and used to construct a Kriging metamodel. The effective properties of all inhomogeneous microstructures within macrostructure are predicted by the built Kriging metamodel and integrated into the analysis of the macroscopic structure. Finally, the material distribution pattern at macroscale is topologically optimized for maximizing the natural frequency of the cellular structure, where the variable thickness sheet (VTS) method combined with an efficient modetracking strategy based on modal assurance criterion (MAC) is employed. In this way, the macrostructural topology, configurations of spatially-varying inhomogeneous microstructures and their global distribution in macrostructure can be simultaneously optimized, which can fully explore the design space for maximizing the natural frequency of the cellular structure. Some numerical examples are provided to test the applications of the Krigingassisted multiscale topology optimization method in maximizing the natural frequency of cellular structures. 2. The core framework of Kriging-assisted multiscale topology optimization A cantilever beam with lumped mass at the bottom right corner is considered as an illustrative example for maximizing the natural frequency subject to a global volume constraint. The core framework of the multiscale topology optimization method assisted by Kriging metamodel is schematically illustrated in Fig. 1. The design domain and boundary conditions of the cantilever beam as well as the initial design of three prototype microstructures are defined in Fig. 1(a). The multiscale optimization method involves three stages. At the stage 1, at microscale, three prototype microstructures are topologically optimized by integrating the numerical homogenization approach into PLSM, as illustrated in Fig. 1(b). At the stage 2, a shape interpolation method is employed to map the optimized prototype microstructures and generate a family of key inhomogeneous microstructures shown in Fig. 1(c), which are considered as sample points and used to construct a Kriging metamodel. The built Kriging metamodel shown in Fig. 1(d) is then employed to predict the effective properties of all the inhomogeneous microstructures within macrostructure. At the stage 3, based on the predicted effective properties of all the inhomogeneous microstructures, the material distribution patterns at macroscale is optimized by VTS method combined with MAC-based mode-tracking strategy, as shown in Fig. 1(e). The optimized cantilever beam is displayed in Fig. 1(f), which is composed of spatially-varying inhomogeneous microstructures obtained by mapping the shape of the corresponding prototype microstructures. To sum up, in this framework of the proposed multiscale topology optimization method, designs of the macrostructure and spatially-varying inhomogeneous microstructures are coupled, so that the most compatible topology of the macrostructure and spatially-varying inhomogeneous microstructures as well as their global distribution in macrostructure can be simultaneously achieved for maximizing the natural frequency of the cellular structure.

Y. Zhang et al. / Computers and Structures 230 (2020) 106197

3

Fig. 1. Schematic of cellular structure design with spatially-varying inhomogeneous microstructures.

3. Kriging-assisted multiscale topology optimization for maximizing natural frequency In this section, to predict the effective property of each inhomogeneous microstructure within the macrostructure, the construction of Kriging metamodel is introduced. Then, the topology optimization problems at macro and micro scales are respectively formulated for maximizing the natural frequency of cellular structure. 3.1. Property prediction of inhomogeneous microstructures Kriging metamodel is originally developed and employed for predictions in mining engineering and geostatistics [50,51]. Recently, it has been migrated and used in the field of structural topology optimization to reduce the computational burdens [52].

In this section, the effective property of each microstructure within the entire macrostructure is predicted by the Kriging metamodel in order to fully explore the design space for macrostructure optimization at a low computational effort. Note that the prototype microstructures refer to ‘‘mother” microstructures that give birth to a set of key microstructures by the shape interpolation technology. These key microstructures will be regarded as sample points in the construction of the Kriging metamodel. 3.1.1. Configuration mapping for key microstructures The key concept of level set-based methods is to implicitly represent the structural boundary as the zero level set of a higherdimensional level set function U, as shown in Fig. 2. One attractive merit of the level set-based boundary representation lies in that some different connectable structural domains can be obtained by setting different level set values of U [30]. Based on this idea,

4

Y. Zhang et al. / Computers and Structures 230 (2020) 106197

Fig. 2. Level-set boundary representation of a prototype microstructure.

an efficient shape mapping method is introduced to produce a series of inhomogeneous microstructures by setting different level set function values of a prototype microstructure UPM , which is represented as follows:

8 PM PM PM > < U ðxÞ > 0; 8x 2 X @ X PM PM 8x 2 @ X U ðxÞ ¼ 0; > : PM U ðxÞ < 0; 8x 2 DPM XPM

ðSolidÞ ðBoundaryÞ

ð1Þ

ðVoidÞ

where DPM is a fixed Eulerian reference domain containing all admissible shapes XPM of microstructures. @ XPM represents the structural boundary and x denotes the point coordinates in DPM . Then, a series of inhomogeneous microstructures can be obtained by interpolating the level set function of the prototype microstructure as follows: PM UMM ðxÞ  ue ; e ðxÞ ¼ U

ðe ¼ 1; 2; :::; neÞ

ð2Þ



    where ue 2 max UPM ; min UPM is a shape-mapping coefficient, i.e., the level set value of UPM , ne is the number of mapped key is the level set function of the eth mapped microstructures. UMM e key microstructure. From Eq. (2), it can be found that a unique characteristic level set function is taken to describe the shapes of all the microstructures, and their corresponding volume fractions (effective densities at macroscale) range from 0 to 1, namely qMM 2 ½0; 1. For the microstructure with special effective density e qMM e , the corresponding shape-mapping coefficient ue can be found by a bi-sectioning algorithm. The zero level set of the eth microstructure, i.e., UMM e ðxÞ ¼ 0, is considered as the structural MM are defined boundary @ XMM e . The corresponding solid regions Xe n o PM MM . 8x 2 D by Ue ðxÞ > 0;

Different inhomogeneous microstructures are featured with different layers of the level set function of a prototype microstructure. In this sense, all the microstructures would offer similar configuration features, especially at their edges. Thus, they are connected well with each other. For example, as shown in Fig. 3, four microstructures with effective densities from 0.3 to 0.7 are naturally connected owing to their highly similar configurations, which are extracted from different layers of the level set function of the prototype microstructure with effective density qPM ¼ 0:5. However, it is difficult to avoid yielding unexpected member part breaks when the of mapped microstructures cover the range effective densities qMM e from 0 to 1, such as the mapped microstructures with effective densities 0.1 and 0.9 shown in Fig. 3(b). Even if the geometrical reinitialization algorithm is employed to preserve the gradient property of the signed-distance level set function [30], member part breaks are not avoided. In this study, the predefined density interval

mapped by a prototype microstructure is considered to avoid yielding member part breaks. Meanwhile, to obtain mapped microstructures with effective densities covering from 0 to 1, multiple prototype microstructures would be used, such as  PM  PM PM qMM m;e 2 qm  a; qm þ a , where qm is the effective density of the mth prototype microstructure, qMM m;e is the effective density of the eth mapped microstructure of the mth prototype microstructure and a is the interval coefficient to ensure that all the mapped density intervals of prototype microstructures cover the range from 0 to 1. In addition, as demonstrated in [30], the unexpected member part breaks would be yielded when the mapped relative density approaches the zero. Thus, a relatively small mapped density value f pm will be set. To allow the occurrence of void regions in optimized macrostructure, a Heaviside projection method will be utilized (see Section 3.2). The aforementioned shape mapping method is implemented without any prior knowledge about the prototype microstructural configuration. Thus, it can achieve the mapped configurations of any shapes without explicit geometric parameters. Meanwhile, it allows for topological changes of the obtained microstructures. Additionally, multiple prototype microstructures can be used to expand the available property space of microstructures, so as to provide more design freedoms for macrostructure optimization.

3.1.2. Property estimation of key microstructures In this paper, the numerical homogenization method [8,9] is employed to estimate the effective properties of the mapped key microstructures. During the numerical implementation, an appropriate number of key microstructures with effective density values qMM m;e would be selected firstly as sample points to construct the MM Kriging metamodel. qMM m;e is an arithmetic sequence and qm;e 2  PM  PM PM PM qm  a; qm þ a ; f pm 6 qm  a 6 qm þ a 6 1. In order to find a level set value ue corresponding to the key microstructure with PM the effective density qMM , the bim;e from the level set function U sectioning algorithm can be utilized. Firstly, it is assumed that the level set value ue corresponding to qMM m;e is located between   the level sets with ul and uu , where the initial ul ¼ max UPM ,  PM  uu ¼ min U , and the corresponding initial effective density qMM ¼ 0 and qMM ¼ 1. Then the effective densities qMM , qMM and u u l l qMM of these microstructures corresponding to u , u and the midl u m dle level set um ¼ ðul þ uu Þ=2 can be calculated by qMM m;e ¼   R  MM  MM MM 1 dXm;e , where Xm;e and Xm;e  are the design domain and XMM m;e jXMM m;e j the area or the volume of the eth mapped key microstructure of the mth prototype microstructure, respectively. By comparing MM qMM ,qMM and qMM m;e , ql m u , the search interval can be narrowed by

5

Y. Zhang et al. / Computers and Structures 230 (2020) 106197

Fig. 3. Schematic example to illustrate the shape mapping method.

the bi-sectioning algorithm, and ul , uu and um can be updated. This process is repeated until the final level set value ue corresponding to qMM m;e is found. The effective elasticity tensor of the mapped key microstructure represented by the level set function UMM m;e is calculated via the numerical homogenization method [8,9] as follows:

  DHijkl

m;e

Z    1 0ðijÞ  MMðijÞ  ¼  e  e u Dpqrs pq pq m;e MM  Xm;e  XMM m;e      MM ðklÞ  e0rsðklÞ  ers uMM H UMM m;e m;e dXm;e

ð3Þ

MM ðijÞ

domain. The displacement um;e corresponding to locally varying   MMðijÞ can be calculated by solving the following strain fields epq um;e equation at microscale with the given macroscopic strain





where

e

0ðijÞ pq :



ðklÞ vMM m;e



ð4Þ

kinematically admissible displacement field satisfying the periodic condition for the eth mapped key microstructure. 3.1.3. Property prediction of microstructures based on the Kriging metamodel The Kriging metamodel corresponding to the mth prototype microstructure assumes the combination of a global model and local deviations:





PM qMA 2 qPM m  a; qm þ a

ð5Þ 



is the where q denotes the prediction point, and D q observed response. The first term on the right side of Eq. (5) pro  vides the global response estimation. f qMA is an ne  1 vector of regression functions and b is a vector of regression coefficients. In this paper, the first order polynomial regression model is used   [53,54]. Z qMA is a realization function of a stochastic process with mean zero and variance r2Z , and its nonzero covariance can be calculated by MA

H

MA

MA MA MA where qMA and i;mk , qj;mk and hmk are the mk-th components of qi , qj





   T  T lD^ H qMA ¼ f qMA ^b þ r qMA R1 Y  F^b 





r^ 2D^ H qMA ¼ r^ 2Z 1 þ uT FT R1 F

   is the is the virtual displacement field, and U XMM m;e

   T   DH qMA ¼ f qMA b þ Z qMA ;

ð7Þ

h, respectively. MK represents the dimensions of training points. The rationale of the Kriging metamodel is to provide the unbiased predictor as Eq. (8). The mean squared error of the Kriging metamodel at the untried point qMA can be written as Eq. (9).



      MM ðijÞ ðklÞ e0pqðijÞ  epq uMM Dpqrs ers v MM H UMM m;e m;e m;e dXm;e ¼ 0;

ðklÞ 8v MM 2 U XMM m;e m;e

where qMA and qMA are two sample points, and R is the stochastic i j process correlation function with respect to unknown correlation parameter vector h. The Gaussian correlation function is the most frequently used, which has the following expression:

mk¼1

the Heaviside function to indicate different parts of the design

XMM m;e

ð6Þ

MK    2

Y MA MA R qMA exp hmk qMA i ; qj ; h ¼ i;mk  qj;mk

where i, j, k, and l are equal to 1, 2, . . ., d, and d is the spatial dimen  is sion. Dpqrs is the elasticity tensor of the solid material. H UMM m;e

Z

  h    MA i MA Cov Z qMA ¼ r2Z R qMA ; Z qj i i ; qj ; h

1

 T   u  r qMA R1 r qMA

ð8Þ

ð9Þ

    where u ¼ FT R1 r qMA  f qMA . ^ b is the generalized least square  MA  is the correlation vector between the untried estimator of b. r q MM MM point qMA and each of the ne training points (i.e. qMM m;1 ; qm;2 ; :::; qm;ne ). R is an ne  ne correlation matrix. Y is the column vector of true responses at the training points. F is an ne  ne matrix with each  T row f qMA .

 1 ^b ¼ FT R1 F FT R1 Y

ð10Þ

  MA MA   T     r qMA ¼ R qMA ; qMA ; q2 ; h ; :::; R qMA ; qMA 1 ;h ;R q ne ; h

ð11Þ

It is worth noting that all of the left terms in Eqs. (7)–(11) are functions of the correlation parameter vector h. The maximum likelihood estimation method in statistics is usually employed to select the value of h, which can be obtained by solving the maximization problem as

    ^ 2Z þ ln½jRj max : LðhÞ ¼  neln r

ð12Þ

where

r^ 2Z ¼

  T 1 Y  F^b R1 Y  F^b ne

ð13Þ

6

Y. Zhang et al. / Computers and Structures 230 (2020) 106197

Finally, the construction of Kriging metamodel is transformed into solving an mk-dimensional unconstrained optimization problem. According to the Gaussian process regression theory, the response of qMA is subject to a normal distribution:

 H  2  MA  ^H  N D ^ qMA ; r ^ ^H q D D

ð14Þ

Therefore, for any microstructure with the effective density qMA in the following optimizations, its effective elasticity tensor can be predicted by the corresponding Kriging metamodel:

    ^ H qMA DHijkl qMA ¼ D

ð15Þ

Obviously, based on the constructed Kriging metamodels, the effective properties of all the microstructures within the whole macrostructure can be predicted at a low computation cost. In this study, unless otherwise specified, ne ¼ 10 key microstructures for the Kriging metamodel of each prototype microstructure are selected. 3.2. Macrostructure topology optimization The dynamic behavior of the macrostructure without damping can be expressed by



 K  x2k M uk ¼ 0

ð16Þ

where K and M denote the global stiffness and mass matrix of the macrostructure, respectively. xk denotes the kth natural frequency and uk is the eigenvector corresponding to xk . Based on the Rayleigh quotient, the natural frequency can be written as

x2k ¼

uTk Kuk uTk Muk

ð17Þ

Here, the VTS method is used to topologically optimize the material distribution at macroscale for maximizing the kth natural frequency of a macrostructure. It is stated as

  MA MA Find : qMA ¼ qMA ðNe ¼ 1; 2; :::; NEÞ 1 ; q2 ; :::; qNe ;  12   uT KðqMA Þu Maximize : xk qMA ¼ uTkM qMA uk Þ k k ð  MA  R MA MA Subject to : G q ¼ MA qNe V 0 dXMA  V MA max 6 0;   MA   MA  X 2 K q  xk M q uk ¼ 0;

ð18Þ

th where qMA Ne denotes the relative density of the Ne macro element as design variable, and NE denotes the total number of the finite elements over the whole macrostructure. xk is the kth natural frequency as the objective function. Note that any eigenfrequency within the above optimization formulation can be regarded as the objective function. The first natural frequency is considered in this

study. G denotes the global volume constraint. V MA is the volume 0 of element, and

K¼

N Z X Ne¼1

M¼

XMA Ne

N Z X Ne¼1

XMA Ne

MA BT DMA ijkl BdXNe

ð19Þ

MA NT qMA Ne NdXNe

ð20Þ

where B and N denote the strain-displacement matrix and shape th function matrix, respectively. XMA Ne denotes the domain of the Ne th macro element. DMA macro eleijkl is the elasticity tensor of the Ne

ment with the density qMA Ne , which can be interpolated by

  MA ~0 DMA ijkl ¼ qNe  DNe

ð21Þ

ijkl

~ 0 is a variable used to identify the equivalent base material where D Ne property for the Neth macro element. During the multiscale optimization process, the elasticity tensor  MA  DMA of the Neth macro element is equal to the homogenized ijkl qNe   of one microstructure with the same elasticity tensor DHijkl qMA Ne density, i.e.,

 MA    DMA ¼ DHijkl qMA Ne ijkl qNe 

ð22Þ

 MA

where DHijkl qNe can be predicted by the aforementioned built Kriging metamodel.   ~0 According to Eqs. (21) and (22), D can be calculated by Ne

  ~0 D Ne

ijkl

¼ DHijkl



denotes the allowable maximum volume of

the macrostructure. XMA denotes the entire macro design domain. qmin ¼ 0:001 and qmax ¼ 1 are the lower and upper bounds of the macro design variable, respectively. It is worth noting that the simple eigenvalue problems are considered in this study. Repeated eigenvalues may occur in the field of eigenfrequencies optimization involving structures with geometric symmetries [36,55]. Unlike simple eigenvalue problems, sensitivity information for repeated eigenvalues cannot be obtained because they are nondifferentiable with respect to the design variable in the common mathematical sense [36]. According to Seyranian [55], this difficulty can be overcome by utilizing the perturbation analysis in the sensitivity analysis [36].

ijkl

 MA qMA Ne =qNe

  ~0 Note that D Ne

ijkl

ð23Þ

is no longer the natural property of the solid

element Ne, and it practically acts as a temporal value varying in the multiscale optimization. Before macrostructure optimization   ~0 of each microstructure needs at each iteration, the term D Ne

qmin < f pm 6 qMA Ne 6 qmax :

V MA max

The global stiffness matrix K and mass matrix M of the macrostructure can be computed by the effective properties of inhomogeneous microstructures within macrostructure.

ijkl

  to be calculated by Eq. (23), where the DHijkl qMA can be and qMA Ne Ne  0 ~ obtained from the previous iteration. Then, the new DNe of ijkl

each microstructure will be used to update qMA Ne . In Eq. (18), f pm is a mapped minimum density value that determines the lower bound of the design variable qMA Ne . In order to allow for the occurrence of void regions during the optimization process, a Heaviside projection function Hpm is utilized as follows: MA q~ MA Ne ¼ H pm qNe

ð24Þ

      ~ H qMA ¼ Hpm DH qMA þ 1  Hpm Dv oid D ijkl Ne Ne ijkl 

ð25Þ



where Hpm ¼ H qMA is equal to 1 for the elements whose Ne  f pm effective densities are larger than f pm ; otherwise, Hpm = 0.

Dv oid ¼ qmin Dsolid

ð26Þ

where Dv oid is the pseudo elasticity matrix of the void element to avoid the singularity when solving Eq. (16). By assigning f pm a relatively small value close to qmin , only a few microstructures will be artificially taken as voids, which guarantees that the numerical convergence has a small effect on the solving accuracy [30]. In this study, an empirical value f pm ¼ 0:05 is used to produce reasonable designs. It is worth nothing that the localized mode [41] is overcome naturally in the proposed method since the penalization between of mass and stiffness is restricted to a fixed value 1 in low density areas, as stated Eq. (26).

7

Y. Zhang et al. / Computers and Structures 230 (2020) 106197

To track the target mode accurately during optimization process, an efficient MAC-based mode-tracking strategy is employed for maximizing the eigenfrequencies of desired modes. The MACbased mode-tracking strategy is very effective and accurate in keeping tracking each of the desired modes, even when the structural topology and configuration change substantially from the initial configuration [56]. In the proposed method, a desired mode is defined as the target mode in the initial design stage. The definition of MAC is stated as

MACðua ; ub Þ ¼ 

 T 2 u ub  a   T ua ua uTb ub

ð27Þ

where ua and ub represent two mode shape vectors, i.e., desired mode shape and extracted mode shape. The value of MAC varies between 0 and 1. When MAC is equal to 1, ua and ub exactly represent the same mode shape. When MAC is equal to 0, ua and ub exactly denote the different mode shapes. The advantage of MACbased mode-tracking strategy is the substantial reduction in computational cost since the mass orthogonalization is eliminated. Meanwhile, it is independent of the order of the eigenfrequency function and thereby has high numerical efficiency.

prototype microstructure. UMI m is the level set function in the design th domain XMI m of the m prototype microstructure, which are interpolated by CSRBFs as follows:

MI MI UMI m ðx; t Þ ¼ um ðxÞam ðt Þ ¼

In the present work, PLSM is employed at microscale to evolve the shape and topology of multiple prototype microstructures. These prototype microstructures will be used to generate a series of inhomogeneous microstructures by the aforementioned shape mapping method, which are considered as sample points for constructing Kriging metamodel. Before the topological optimization of multiple prototype microstructures, the density clustering is implemented based on the material distribution obtained by solving Eq. (18), in order to efficiently implement the macro finite element analysis by using the effective properties of multiple prototype microstructures. Here, all the macro elements with the  PM  PM density qMA Ne 2 qm  a; qm þ a will be represented by a unique ~ MA Ne

after prototype microstructure with the density qPM m . The term q the density clustering operation can be defined by the following heuristic scheme [33]:

(

q~ MA Ne ¼

PM m1 PM m

q q

if if

MA PM PM MA Ne  m1 6 m  Ne MA PM PM MA m  Ne Ne  m1 >

q

q

q

q

q

q

q

q

ð30Þ

h i MI MI MI where uMI is CSRBFs vector. m ðxÞ ¼ um;1 ðxÞ; um;2 ðxÞ; :::; um;N ðxÞ h iT MI MI MI aMI is the vector of actual design m ðt Þ ¼ am;1 ðt Þ; am;2 ðt Þ; :::; am;N ðt Þ variables of the mth prototype microstructure. aMI m;n denotes expansion coefficient of the CSRBF interpolation of the nth knot in the micro level set grid of the mth prototype microstructure. In Eq. (29), uk is the eigenvector corresponding to the kth natural frequency xk , which is calculated by substituting the effective properties DHijkl of all the prototype microstructures into the corresponding      state equation K qMA ; DHijkl  x2k M qMA uk ¼ 0 at macroscale. V MA denotes the volume of a macroscale finite element with solid 0 material. Gm is the local volume fraction constraint for the mth pro~ MI ~ MI totype microstructure. a min ¼ 0:001 and amax ¼ 1 are the lower and ~ m;n is regularized design variable, which ~ m;n . And a upper bounds of a will be used in the optimization algorithm to facilitate the numerical implementation. The moving of a prototype microstructural boundary towards to its optimum is equivalent to solving the following first-order Hamilton-Jacobi partial differential equation (H-J PDE) [6,7]: MI

  @ UMI MI m ðx; t Þ    tMI m  rUm ðx; t Þ ¼ 0; @t

MI UMI m ðx; 0Þ ¼ Um0 ðxÞ

ð31Þ

where t is an artificial pseudo-time to enable the dynamic motion of the shape deformations, and tMI m denotes the normal velocity field. Substituting Eq. (30) into Eq. (31), the H-J PDE is transformed into a new form of an ordinary differential equation (ODE) [6,7]:

uMI m ðxÞ

daMI MI MI m ðt Þ  tMI m  jrum ðxÞam ðt Þj ¼ 0 dt

ð32Þ

Thus, the normal velocity field tMI m can be denoted as follows:

 tMI m ¼ 

ru

u

MI m ðxÞ  MI ðxÞ MI ðt Þ m m

a

_ MI a_ MI m ðt Þ; where am ðt Þ ¼

daMI m ðt Þ dt

ð33Þ

ð28Þ 4. Sensitivity analysis

Thus, a limited number of prototype microstructures will be topologically optimized to represent all microstructures with the   PM density qMA of each prototype Ne 2 f pm ; 1 . The density qm microstructure is regarded as their volume fraction constraint during optimization process. Using the PLSM with the compactly supported radial basis functions (CSRBFs) [26,57–59], the optimization formulation of multiple prototype microstructures can be mathematically expressed as

Find : aMI n ¼ 1; 2; :::; NÞ m;n ðm ¼ 1; 2; :::; M;  112 0 P R M uTk BT DH UMI aMI ÞÞBdXMA uk  MI  m ijkl ð m ð m¼1 XMA  A Maximize : xk a ¼ @ PM R m MA T MI MI uTk N H U a Nd X u ð m ð ÞÞ m k m¼1 XMA m   MI  R  MI  MI MI MA Subject to : Gm a ¼ XMI H Um a dXm  qPM m V 0 6 0; m       K qMA ; DHijkl  x2k M qMA uk ¼ 0;

MI uMI n ðxÞan ðt Þ

n¼1

MI

3.3. Prototype microstructure topology optimization

N X

In the proposed multiscale optimization method, a gradientbased algorithm, such as the method of moving asymptotes (MMA) [60], is employed to update the design variables at both scales. Hence, it is necessary to calculate the first-order derivatives of the objective and constraint functions with respect to the design variables. In this section, the sensitivity analysis at both scales is carried out. 4.1. Sensitivity analysis at macroscale

ð29Þ

~ MI ~ MI a~ MI min 6 am;n 6 amax : where M is the type number of prototype microstructures and N is the total number of the knots in the micro level set grid for a

With regard to the Rayleigh quotient in Eq. (17), the sensitivity of the objective function xk with respect to the design variable qMA Ne can be written as

  @ xk qMA 1 Ne ¼ 2xk uTk Muk @ qMA Ne 

 2@uTk  @K @M  K  x2k M uk þ uTk  x2k MA uk MA MA @ qNe @ qNe @ qNe ð34Þ

8

Y. Zhang et al. / Computers and Structures 230 (2020) 106197

With the help of Eq. (16), Eq. (34) can be simplified as

  

@ xk qMA 1 @K Ne T 2 @M uk ¼ u  x k k 2xk uTk Muk @ qMA @ qMA @ qMA Ne Ne Ne

ð35Þ

Recalling Eqs. (19)–(21), the derivatives of the global stiffness K and global mass matrix M with respect to the design variable qMA Ne can be calculated as follows: NE Z   X @K ~0 ¼ BT D BdXMA Ne Ne MA MA ijkl @ qNe Ne¼1 XNe

ð36Þ

NE Z X @M ¼ NT NdXMA Ne MA @ qMA Ne Ne¼1 XNe

ð37Þ

As a result, the sensitivity of the kth eigenfrequency xk at macroscale can be stated as

"   NE Z   X @ xk qMA 1 Ne T ~0 ¼ u BT D k Ne T MA ijkl 2 x u Mu @ qMA k k k Ne Ne¼1 XNe ! # Z NE X 2 uk NT NdXMA  BdXMA Ne  xk Ne Ne¼1

where



~0 D Ne

 ijkl

XMA Ne

ð38Þ

Eqs. (23)–(25).

  ~0 D Ne

ijkl

¼

 MA





Hpm DHijkl qNe þ 1  Hpm Dv oid Hpm qMA Ne

ð39Þ

The sensitivity of the volume constraint with respect to the design variable qMA Ne can be obtained by

  @G qMA Ne ¼ V MA 0 @ qMA Ne

th XMI prototype microstructure, respectively. m of the m

Eqs. (42) and (43) can be rewritten as the following form by substituting the normal velocity tMI m defined in Eq. (33).

  MI  @DHijkl UMI m a @t

0 Z  N X  ðijÞ  B 1  e0pqðijÞ  epq uMI ¼ Dpqrs @ MI  m MI Xm  Xm n¼1 











 @ aMI m;n @t ð44Þ

MI MI ðklÞ e0rsðklÞ  ers uMI uMI  m m;n ðxÞ  d Um dXm

0 1   MI  Z MI N   X @H UMI a 1 B C @ am;n MI MI m uMI ðxÞ  d Um dXm A  ¼ @ MI  m;n @t @t Xm  XMI m n¼1 ð45Þ

can be written as the following form considering



  1 where d UMI ¼ p  MI n 2 2 is the derivative of the Heaviside funcm ðUm Þ þn  MI  tion H Um , and n is chosen as 2–4 times the mesh size from MI numerical experience [26,59]. uMI m and tm denote the microscale displacement filed and the normal velocity in the design domain

ð40Þ

On the other hand, the derivatives of the effective elasticity ten     MI  and H UMI with respect to t can be sor DHijkl UMI aMI m a expressed by the chain rule:

  MI    MI  N X @DHijkl UMI @DHijkl UMI @ aMI m a m a m;n  ¼ MI a @t @t @ m;n n¼1

ð46Þ

  MI    MI  N X @ aMI @H UMI @H UMI m;n m a m a ¼  @t @ aMI @t m;n n¼1

ð47Þ

Comparing the corresponding terms in Eqs. (44) and (46), as well as Eqs. (45) and (47), the derivatives of the effective elasticity      MI  tensor DHijkl UMI aMI and H UMI with respect to the design m a variable aMI m;n can be stated as

4.2. Sensitivity analysis at microscale In this section, the shape derivative [6,7] is introduced to conduct the sensitivity of boundary perturbations with respect to the time variable t at the microscale. The shape derivatives of   xk aMI can be calculated by

"  MI  MI    H M Z X @ xk aMI 1 T @Dijkl Um a T ¼ u B BdXMA m @t @t 2xk uTk Muk k m¼1 XMA m ! #     M Z MI X @H UMI m a uk ð41Þ NdXMA NT x2k m MA @t X m m¼1 Based on [15,26,32], the derivative of the effective elasticity     tensor DHijkl UMI aMI of the prototype microstructure and H UMI m  MI  a Þ with respect to t can be stated as Eqs. (42) and (43), respectively.

  MI  @DHijkl UMI m a @t

Z   ðijÞ  1  ¼  e0pqðijÞ  epq uMI Dpqrs m MI  MI Xm  Xm       MI  MI 0ðklÞ MIðklÞ MI    ers  ers um tMI m  rUm d Um dXm ð42Þ

  MI  Z    MI  MI @H UMI 1 MI  m a   tMI ¼  m  rUm d Um dXm  MI MI @t Xm  Xm

ð43Þ

  MI  Z   ðijÞ  @DHijkl UMI 1 m a   ¼ e0pqðijÞ  epq uMI Dpqrs m MI   MI MI @ am;n Xm  Xm     ðklÞ  MI MI  e0rsðklÞ  ers uMI um;n ðxÞ  d UMI m dXm m ð48Þ   MI  Z   @H UMI 1 m a MI  dXMI ¼  u ðxÞ  d UMI m m m;n MI  MI MI @ am;n X Xm  m

ð49Þ

Substituting Eqs. (46) and (47) into Eq. (41) yields   @ xk aMI 1 ¼ @t 2xk uTk Muk 2 0 M 1 3 N H MI MI

MI P R MA T P @Dijkl ðUm ða ÞÞ @ am;n Bd  X MA B MI m 6 B C 7 X @t @ am;n n¼1 6 B m¼1 m C 7  6uTk B Cuk 7 N

M R MI MI MI 4 @ A 5 P P @HðUm ða ÞÞ @ am;n MA T 2 Nd xk  X MA N MI m X @t @a m¼1

m

n¼1

m;n

ð50Þ 

The derivative of the objective function xk a t can be expressed by the chain rule as

    M X N X @ xk aMI @ xk aMI @ aMI m;n  ¼ @t @t @ aMI m;n m¼1 n¼1

MI



with respect to

ð51Þ

Y. Zhang et al. / Computers and Structures 230 (2020) 106197

Again, comparing the corresponding terms in Eqs. (50) and (51),   the derivative of the objective function xk aMI with respect to the design variable aMI m;n can be obtained as



 MI

@ xk a @ aMI m;n

"

  MI  Z @DHijkl UMI 1 m a T T ¼ u B BdXMA m @ aMI 2xk uTk Muk k XMA m;n m ! #   MI  Z @H UMI m a NT NdXMA uk x2k m @ aMI XMA m;n m

9

Similarly, the derivative of the local volume constraints   Gm aMI with respect to the design variable aMI m;n can be given by

  Z   @Gm aMI MI ¼ uMI ðxÞ  d UMI m dXm m;n MI @ am;n XMI m

ð53Þ

5. Numerical implementation

ð52Þ

Fig. 4 shows the flowchart of numerical implementation of the Kriging-assisted multiscale topology optimization method, which

Fig. 4. Flowchart of the proposed method.

Y. Zhang et al. / Computers and Structures 230 (2020) 106197

contains three stages, i.e., topology optimization of multiple prototype microstructures at microscale, effective property prediction of microstructures within macrostructure by Kriging metamodel, and material distribution optimization at macroscale. Firstly, the effective property DH of each prototype microstructure is calculated by the homogenization method (Eq. (3)). Based on the material distri~ MA at macroscale obtained by clustering operation (Eq. bution q (28)), the macro finite element analysis (FEA) is implemented via Eq. (16). The sensitivity information at microscale is calculated via Eqs. (52) and (53). The MMA is employed to update the micro design variable aMI . The new shape of each prototype microstructure UMI m is output. Secondly, an appropriate number of key microstructures are obtained by mapping the shape of the corresponding prototype microstructure UMI m (Eq. (2)), and their effective properties DHe are calculated by the homogenization method (Eq. (3)). These key microstructures are considered as sample points and used to construct a Kriging metamodel in the corresponding density interval via Eqs. (5)–(14). The built Kriging metamodel is then employed to predict the effective properties of all the microstructures within the corresponding density interval. Finally, ^ H. the macro FEA is implemented via Eq. (16) using the predicted D

terminate when the difference of the objective function values between two successive iterations is lower than 104, or the maximum 200 iteration steps of the multiscale optimization is reached. 6.1. Double-clamped structure The first example investigates the multiscale optimization of the double-clamped structure with length L ¼ 1:4m and height H ¼ 0:2m, which is clamped at both ends as shown in Fig. 5. A concentrated lump mass of the magnitude 2000 is placed at the center of the design domain. 210  30 ¼ 6300 four-node quadrilateral 1 1 elements with length l ¼ 150 m and height h ¼ 150 m are applied to discretize the macro design domain. The objective function is to maximize the first natural frequency x1 of the double-clamped structure under a global volume constraint of 40%.

1.4m

0.2m

10

e

The sensitivity information for macroscale optimization can be calculated via Eqs. (38) and (40). The MMA is employed to update the macro design variables qMA . The multiscale optimization process repeats until the convergent criterion is satisfied. With the help of the shape interpolation method, the microstructures generated by mapping a prototype microstructure in the predefined density interval will be featured with similar configuration features. Thus, they can be well connected with each other. However, when multiple prototype microstructures are employed, the connectivity between the different prototype microstructures is also important to ensure the connectivity of all the inhomogeneous microstructures within macrostructure. In this study, the kinematical connective constraint approach [61] is employed to guarantee the connectivity of the adjacent prototype microstructures. Note that these predefined connectors will more or less limit the design space in multiscale optimization. Hence, they may slightly compromise the performance of the optimized structure. However, in engineering, it is acceptable to reasonably sacrifice the structural performance in order to achieve a manufacturable design [26,32].

Fig. 5. Sketch of double-clamped structure.

(a) Optimized macrostructural topology

(b) The first mode shape of optimized macrostructure

(c) The second mode shape of optimized macrostructure

6. Numerical examples In this section, three numerical examples are presented to test the performance of the Kriging-assisted multiscale optimization method for maximizing the first natural frequency of a cellular structure. Specifically, the first example is investigated for the validity and advantage of the proposed method, the second example for its effectiveness under different structural geometries, and the third example for its effectiveness under various lumped masses. For all examples, the materials are subject to plane stress conditions. The base material has Young’s modulus E0 ¼ 201 GPa, Poisson’s ratio l ¼ 0:3, and density q0 ¼ 7:8  103 kg=m3 . For simplicity, four-node quadrilateral elements are applied in the finite element discretization at both scales, and each finite element at macroscale is represented by an individual microstructure at microscale. Each microstructure is discretized into 50  50 ¼ 2500 four-node quadrilateral elements. In this study, unless otherwise specified, four prototype microstructures (i.e., qPM m ¼ 0:2; 0:4; 0:6; 0:8) are employed to generate the inhomogeneous microstructures. In the mapping density interval of each prototype microstructure, a total number of ne ¼ 10 key microstructures are used to construct the Kriging metamodel. The optimization will

(d) The three mode shape of optimized macrostructure

(e) Optimized multiscale structure with details of local regions

(f) The first mode shape of optimized multiscale structure Fig. 6. Multiscale design of double-clamped structure with first nature frequency 136.6298 rad/s.

11

Y. Zhang et al. / Computers and Structures 230 (2020) 106197 Table 1 Microstructural configurations, level set functions, first mode shapes, and effective elastic properties of four optimized prototype microstructures. Microstructural configurations

Level set functions

(a) Optimized prototype microstructure with effective density 0.2

Effective elastic properties DH (GPa) 2 3 12:0132 11:4228 0:0639 4 11:4228 12:0423 0:0880 5 0:0639 0:0880 10:5642

2

29:8611 4 24:6585 0:7344

(b) Optimized prototype microstructure with effective density 0.4

2

88:6822 4 26:9643 0:5098

(c) Optimized prototype microstructure with effective density 0.6

2

149:4142 4 35:4836 0:0744

24:6585 29:8131 0:8576

3 0:7344 0:8576 5 21:8949

26:9643 46:1609 0:5479

3 0:5098 0:5479 5 28:5998

35:4836 82:7998 0:0497

3 0:0744 0:0497 5 43:0271

(d) Optimized prototype microstructure with effective density 0.8

Fig. 7. Effective property versus effective density for the mapped microstructures and their configurations.

6.1.1. Multiscale design of double-clamped structure The multiscale design of double-clamped structure is shown in Fig. 6, including the optimized macrostructural topology, the first three mode shapes of the optimized macrostructure, the optimized multiscale structure with spatially-varying microstructures and its

first mode shape. The detailed microstructural member parts in some local regions are magnified in Fig. 6(e). From Fig. 6, it can be seen that the multiscale design contains a number of inhomogeneous microstructures with intermediate densities and they are well connected with each other. The external frame region of

12

Y. Zhang et al. / Computers and Structures 230 (2020) 106197

optimized macrostructural topology is occupied by a number of completely solids and high-density microstructures to effectively resist the deformation of the double-clamped structure in its first vibration mode, which accords with the widely accepted design for the double-clamped structure. To verify the possibility of local vibration of each inhomogeneous microstructure, the modal analysis of optimized cellular structure is implemented directly. It can be observed from Fig. 6(f) that the local vibration mode does not occur in any microstructure in the optimized cellular structure. The microstructural configurations, level set functions, first mode shapes, and effective elastic properties of four optimized prototype microstructures with effective densities 0.2, 0.4, 0.6, and 0.8 are presented in Table 1. The connectivity of the four prototype microstructures is guaranteed by setting predefined connectors between adjacent prototype microstructures. From

Table 1, it can be observed that the four prototype microstructural geometries are symmetrical due to the symmetrical geometry of the double-clamped structure, even without symmetrical constraints imposed on prototype microstructures during optimization. For each of key microstructures obtained by mapping each prototype microstructure, its effective property versus the effective density, together with its configuration are plotted in Fig. 7, where the terms DH ð1; 3Þ and DH ð2; 3Þ are eliminated since the effective properties of microstructures are almost orthotropic. Based on the shape interpolation technology, the microstructures within the same density interval are featured with similar configuration features, especially at their edges, and thus they are well connected with each other. From Fig. 7, it can be noticed that microstructures obtained by mapping prototype microstructures behave in an orthotropic form to flexibly offer directional stiffness, and their

Fig. 8. Iterative histories of the objective function and global volume function, together with intermediate optimized macrostructural topologies.

Fig. 9. Iterative histories of local volume constraints of four prototype microstructures, together with their initial designs and intermediate optimized configurations.

13

Y. Zhang et al. / Computers and Structures 230 (2020) 106197

Fig. 8 shows iterative histories of the objective function and global volume constraint, together with some intermediate optimized results. It is illustrated that the macrostructural topology, spatially-varying microstructural configurations and their global distribution in macrostructure are simultaneously optimized. As seen in Fig. 8, the first natural frequency gradually decreases as the global volume fraction decreases at the first few iterations

effective properties are bounded by those of the microstructure with the minimum density f pm ¼ 0:05 and the solid one. The microstructures with effective density from f pm ¼ 0:05 to 1 can offer a large design space; meanwhile, the topological variations of four prototype microstructures can further extend the design space to sufficiently maximize the first natural frequency of the cellular structure.

Fig. 10. Iteration histories of the first three natural frequencies of the double-clamp structure.

Table 2 Comparisons between the effective properties by homogenization method and Kriging metamodel. Microstructural configurations

Effective densities

Effective elastic property DH(GPa)

0.13

Homogenization method 2 3 7:4001 7:0906 0:0587 4 7:0906 7:4164 0:0715 5 0:0587 0:0715 6:7395

Kriging metamodel 2 7:3972 7:0881 4 7:0881 7:4138 0:0605 0:0763

2

2

0.45

35:8513 4 28:1967 0:6744 2

0.66

104:0611 4 29:3629 0:4916 2

0.88

178:1283 4 40:3635 0:0928

Maximum error

3 0:6744 0:8147 5 24:9440

28:1967 35:8023 0:8147

29:3629 54:7926 0:5545

40:3635 104:5275 0:0048

35:8389 4 28:2012 0:6781

3 0:4916 0:5545 5 32:2050

2

104:0540 4 29:3656 0:4886

3 0:0928 0:0048 5 52:2707

2

178:1587 4 40:3821 0:0782

3 0:0605 0:0763 5 6:7371 3 0:6781 0:8173 5 24:9457

28:2012 35:7957 0:8173

29:3656 54:7894 0:5540

40:3821 104:4238 0:0006

3 0:4886 0:5540 5 32:2043 3 0:0782 0:0006 5 52:2985

0.039%

0.034%

0.009%

0.046%

Table 3 Optimized the first natural frequency with different numbers of prototype microstructures. M

2

3

4

9

18

First natural frequency

131.9879

133.8084

136.6298

136.7988

137.3736

Table 4 Optimized the first natural frequency with different numbers of key microstructures. ne

5

10

20

40

First natural frequency

136.2816

136.6298

136.6056

136.3488

14

Y. Zhang et al. / Computers and Structures 230 (2020) 106197

(a) Monoscale microstructural design with the first natural frequency 97.6321 rad/s

(b) Monoscale macrostructural design with the first natural frequency 128.2650 rad/s

(a) Optimized macrostructural topology (c) Multiscale design with the first natural frequency 136.6298 rad/s Fig. 11. Comparison of optimized results between the multiscale design and the monoscale designs.

and then increases. Finally, it converges to 136.6298 rad/s when the global volume fraction constraint is satisfied. It takes 115 iterative steps for the multiscale optimization to achieve a convergent solution. This reveals the efficiency and validity of the proposed multiscale method for maximizing the first natural frequency of a cellular structure. Fig. 9 shows convergence histories of local volume fractions of four prototype microstructures along with some intermediate optimized microstructural configurations. It can be seen that the convergence configurations of the four prototype microstructures are quickly obtained after about 20 iterations steps. Fig. 10 shows the evolution histories of the first three natural frequencies, x1, x2 and x3. It is noticed that the first three natural frequencies have the similar variation tendency. Specifically, due to the use of MACbased mode-tracking strategy, there is no mode switch between the first and second natural frequencies. To test the prediction accuracy of Kriging metamodels for effective properties of microstructures, four microstructures with effective densities of 0.13, 0.45, 0.66 and 0.88 are selected, which are obtained by mapping the four optimized prototype microstructures with effective densities of 0.2, 0.4, 0.6 and 0.8, respectively. As shown in Table 2, the effective properties of selected four microstructures predicted by constructed Kriging metamodels (Eq. (15)) are slightly different from those calculated directly by the homogenization method (Eq.(3)). The maximum error of Kriging metamodels for the four microstructures is 0.046%. Moreover, the effective properties of all mapped microstructures within the entire macrostructure are calculated directly by the homogenization method. The obtained first natural frequency of the cellular structure is xHM ¼ 136:5803 rad/s, which is very close to that 1 obtained by the use of the Kriging metamodels, i.e., x1 ¼ 136:6298 rad/s. Thus, it is demonstrated that the Kriging metamodel has high accuracy in prediction of the effective properties of microstructures. In addition, the comparison of optimized results is conducted under different numbers of prototype microstructures and selected key microstructures for each prototype microstructure. Table 3

H

L

Fig. 12. Sketch of pinned beam.

(b) The first mode shape of optimized macrostructure

(c) Optimized multiscale design Fig. 13. Multiscale design of the pinned beam with geometry H  L = 0.2 m  0.2 m.

lists the optimized first natural frequency of the double-clamped structure under different numbers of prototype microstructures, in which ten key microstructures are used for each prototype microstructure. On the other hand, Table 4 presents the optimized first natural frequency of the double-clamped structure under different numbers of selected key microstructures for each prototype microstructure, where four prototype microstructures are employed. Because the optimized configurations of the microstructures and their distribution within macrostructure are almost the same as those in Fig. 7 and Table 1, they are not presented here. From Table 3, it can be seen that a larger first natural frequency can be achieved when more prototype microstructures are used. This is because more prototype microstructures can obtain a larger available material property space. However, topology optimization of more prototype microstructures will remarkably increase the computational cost. From Table 4, it can be observed that the optimized first natural frequency has a very small variation as the increase of the number of selected key microstructures for each prototype microstructure. From the optimized results in Tables 3 and 4, it can be preliminarily found that four prototype microstructures and ten key microstructures for each prototype microstructure would be a reasonable choice to achieve a balance among

15

Y. Zhang et al. / Computers and Structures 230 (2020) 106197

(a) Optimized macrostructural topology

(a) Optimized macrostructural topology (b) The first mode shape of optimized macrostructure

(c) Optimized multiscale design Fig. 15. Multiscale design of the pinned beam with geometry H  L = 0.2 m  0.8 m.

Find : aMI n ð n ¼ 1; 2; :::; N Þ   Maximize : xk uk ; aMI    P  MI  MI R MI a Subject to : G aMI ¼ NE dX  V max 6 0; Ne¼1 XMI H U       K DHijkl UMI  x2k M UMI uk ¼ 0;

(b) The first mode shape of optimized macrostructure

~ MI ~ MI a~ MI min 6 an 6 amax : ð55Þ where aMA and aMI n are design variables of macrostructural optik mization and microstructural optimization, respectively. K and N are the total number of knots in the macro and micro level set grid, respectively. uk is the eigenvector corresponding to the kth natural frequency xk , which is the objective function. G denotes the global structural volume constraint. V max denotes the allowable maximum

(c) Optimized multiscale design Fig. 14. Multiscale design of the pinned beam with geometry H  L = 0.2 m  0.4 m.

the structural performance, the computational burdens and accuracy. 6.1.2. Comparison with single-scale designs To illustrate the advantage of the Kriging-assisted multiscale optimization method for maximizing the first natural frequency of the cellular structure, the optimized multiscale design is compared with the monoscale designs, i.e., monoscale macrostructural and microstructural designs, under the same global volume constraint of 40%. For fair and meaningful comparison, the PLSM is used to optimize the designs of monoscale macrostructure and microstructure. The corresponding optimization models are formulated in Eqs. (54) and (55), respectively.

Find : aMA ð k ¼ 1; 2; :::; K Þ k   Maximize : xk uk ; aMA    R   Subject to : G aMA ¼ XMA H UMA aMA dXMA  V max 6 0;      K UMA  x2k M UMA uk ¼ 0;

volume. UMA and UMI are the corresponding level set function in the macro and micro design domains XMA and XMI . The meshes of 210  30 ¼ 6300 and 50  50 ¼ 2500 with four-node quadrilateral elements are respectively used in monoscale macrostructural and microstructural designs, which correspond to the two mesh schemes used at macroscale and microscale in the proposed multiscale optimization method. For all optimization designs, the base materials are the same and subject to plane stress conditions. The optimized designs by the monoscale macrostructural and microstructural design method are shown in Fig. 11(a) and (b), respectively. For expedient comparison, the optimized multiscale design with global volume constraint 40% is presented here again, as shown in Fig. 11(c). It is obvious that the multiscale design has a larger first natural frequency than the two monoscale designs, especially compared with the monoscale microstructural design. This is mainly because the Kriging-assisted multiscale optimization method can sufficiently explore the design space by simultaneously optimizing the macrostructural topology, the configurations of spatially-varying inhomogeneous microstructures and their global distribution within macrostructure. 6.2. Pinned beam with different geometries

~ MA ~ MA a~ MA min 6 ak 6 amax : ð54Þ

In this example, the Kriging-assisted multiscale design method is applied to optimize a pinned beam as depicted in Fig. 12. The

16

Y. Zhang et al. / Computers and Structures 230 (2020) 106197

beam is simply point-pinned at the middle of both ends. A concentrated lump mass of the magnitude 2000 is located at the center of the design domain. Three different geometries of the macro design domain are considered, i.e., the height H = 0.2 m and the length L = 0.2 m, 0.4 m, 0.8 m. And three kinds of finite element mesh are used to discretize the macro design domain for the three cases, i.e., 40  40 elements, 40  80 elements and 40  160 elements, where each finite element is characterized with four-node quadrilateral elements with length l ¼ 0:005 m and height h ¼ 0:005 m. The goal is to maximize the first natural frequency of the pinned beam under a global volume constraint of 40%. Optimization results of the pinned beam under three different geometries by the Kriging-assisted multiscale optimization

method are respectively presented in Figs. 13–15, including the optimized macrostructural topologies, their first vibration modes, and the optimized macrostructures with spatiallyvarying inhomogeneous microstructures. Fig. 16 present the effective property of the mapped key microstructures versus the effective densities from f pm ¼ 0:05 to 1, and four optimized prototype microstructures together with their corresponding effective properties for the pinned beam under three different geometries. From Figs. 13(c), 14(c) and 15(c), it can be seen that the completely solids and high-density microstructures are mainly distributed in the external frame regions of the multiscale design, and the shear resistant microstructures with low densities are

(a) The geometry of the pinned beam H×L=0.2m×0.2m

(b) The geometry of the pinned beam H×L=0.2m×0.4m Fig. 16. Effective properties of the mapped microstructures versus effective densities, and four prototype microstructures together with their effective properties under different geometries.

Y. Zhang et al. / Computers and Structures 230 (2020) 106197

(c) The geometry of the pinned beam H×L=0.2m×0.8m Fig. 16 (continued)

(a) The geometry of the pinned beam H×L=0.2m×0.2m

(b) The geometry of the pinned beam H×L=0.2m×0.4m

(c) The geometry of the pinned beam H×L=0.2m×0.8m Fig. 17. Iterative histories of the objective function and global volume constraint under different geometries.

17

18

Y. Zhang et al. / Computers and Structures 230 (2020) 106197

Fig. 18. Sketch of cantilever beams with various lumped masses.

(a) Optimized macrostructural topology

(b) The first mode shape of optimized macrostructure

(c) Optimized multiscale design

(d) Effective properties of the mapped microstructures versus effective densities, and four prototype microstructures together with their effective properties Fig. 19. Multiscale design under case Ⅰ, and the first natural frequency is 125.3763 rad/s.

Y. Zhang et al. / Computers and Structures 230 (2020) 106197

placed in the interval regions, which can effectively improve the macrostructural capacity to resist the dominant vertical deformation of the pinned beams in their first vibration modes as shown in Figs. 13(b), 14(b) and 15(b). By mapping four optimized prototype microstructures, all the microstructures within the macrostructures are well connected with each other. From Fig. 16, it is noticed that four optimized prototype microstructures have remarkable orthotropic properties to flexibly offer directional stiffness, and the effective properties of all the mapped microstructures are bounded by those of the microstructure with the mini-

(a) Optimized macrostructural topology

19

mum density f pm ¼ 0:05 and the solid one, which enable to supply enough large design space. Also, from Fig. 16, it can be observed that the geometry of the macrostructure affects substantially the final configurations of the prototype microstructures. With the increment of the length of the pinned beam, the microstructural stiffness in x-direction DH ð1; 1Þ (line with blue solid square marks) increases and the microstructural stiffness in y-direction DH ð2; 2Þ (line with green solid rhombus marks) decreases to make the pinned beam, in order to well resist bending deformation in its first vibration mode.

(b) The first mode shape of optimized macrostructure

(c) Optimized multiscale design

(d) Effective properties of the mapped microstructures versus effective densities, and four prototype microstructures together with their effective properties Fig. 20. Multiscale design under case II, and the first natural frequency is 90.5848 rad/s.

20

Y. Zhang et al. / Computers and Structures 230 (2020) 106197

Fig. 17 show iterative histories of the first natural frequency and the global volume fraction under three geometries. During the optimization process, the macrostructural topology, configurations of all the microstructures and their distributions within the macrostructure are simultaneously optimized. The first natural frequency of the pinned beam under three geometries converges to the final values of 495.5710 rad/s, 333.2769 rad/s and 165.0631 rad/s after 89, 106 and 113 iterative steps, respectively. The first natural frequency decreases gradually as the length of the pinned beam increases. Thus, it is illustrated that the proposed multiscale method is effective for maximizing the first natural fre-

(a) Optimized macrostructural topology

quency of the cellular structure under different structural geometries. 6.3. Cantilever beam with various lumped masses To validate the effectiveness of the Kriging-assisted multiscale optimization method, the cantilever beam with various lumped masses is considered in this example. The cantilever beam with length L ¼ 0:4 m and height H ¼ 0:2 m is shown in Fig. 18. Various lumped masses are allocated at the macro design domain, and each lumped mass is featured with the magnitude 2000. In this

(b) The first mode shape of optimized macrostructure

(c) Optimized multiscale design

(d) Effective properties of the mapped microstructures versus effective densities, and four prototype microstructures together with their effective properties Fig. 21. Multiscale design under case III, and the first natural frequency is 80.4713 rad/s.

Y. Zhang et al. / Computers and Structures 230 (2020) 106197

example, three cases are considered. In case Ⅰ, a lumped mass is placed at the middle point of the right free edges; in case Ⅱ, a lumped mass is placed at the top right corner; in case III, two equal-quality lumped masses are placed at the top and bottom right corners, respectively. The macro design domain is discretized into 80  40 four-node quadrilateral elements with length l ¼ 0:005 m and height h ¼ 0:005 m. The objective function is to maximize the first natural frequency of the cantilever beam under a global volume constraint of 40%. The multiscale design results of the cantilever beam under case Ⅰ is presented in Fig. 19, and the first natural frequency converges to 125.3763 rad/s. Fig. 19(a) and (b) show the optimized

21

macrostructural topology and its first vibration mode, respectively. The optimized multiscale design with spatially-varying inhomogeneous microstructures is presented in Fig. 19(c). Fig. 19(d) shows the effective property of the mapped microstructures versus effective densities from f pm ¼ 0:05 to 1, and four optimized prototype microstructures together with their effective properties. Similarly, Figs. 20 and 21 present the multiscale design results of the cantilever beam under the other two cases. The first natural frequency converges to 90.5848 rad/s and 80.4713 rad/s, respectively. Iteration histories of the first three natural frequencies and the global volume constraint for the cantilever beam under three cases are plotted in Fig. 22.

Fig. 22. Iteration histories of the first three natural frequencies of the cantilever beam under different cases.

22

Y. Zhang et al. / Computers and Structures 230 (2020) 106197

Fig. 22 (continued)

As shown in Figs. 19(c), 20(c) and 21(c), it can be seen that the external frame regions of the multiscale design are enhanced by the completely solids and high-density microstructures. The interval regions are occupied by a number of shear resistant microstructures with low densities. These designs are effective to resist the dominant deformation in their corresponding first vibration modes as shown in Figs. 19(b), 20(b) and 21(b). By mapping four optimized prototype microstructures, all microstructures within the macrostructures are well connected with each other. From Figs. 19 (d), 20(d) and 21(d), it can be observed that the effective properties of the microstructures are bounded by those of the microstructure with minimum density f pm ¼ 0:05 and the solid one, which can provide sufficient exploration for the multiscale design space and four optimized prototype microstructures with orthotropic properties further broaden the advantage. As for four optimized prototype microstructures under different cases, they behave different orthotropic properties. It is interestingly noted that the four resulting prototype microstructural topologies are highly symmetric about x and y axes (as shown in Figs. 19(d), 20(d) and 21(d)) when the boundary conditions of the macrostructures are symmetric, and highly asymmetric (as shown in Fig. 20(d)) when the boundary conditions of the macrostructure are asymmetric. From Fig. 22, it can be observed that under the three cases, the first natural frequency and global volume fraction have the same iterative tendency. The multiscale designs take 97, 114 and 112 iterative steps to achieve the convergent solutions, respectively. With MAC-based mode-tracking strategy during the optimization process, there is no mode switch between the first and second frequencies under three cases, although the first three natural frequencies have the similar variation tendencies as shown in Fig. 22. Therefore, it is illustrated that the proposed multiscale method is effective for maximizing the first natural frequency of the cellular structure under various lumped masses.

of inhomogeneous cellular structures. In the proposed method, both spatially-varying microstructural configurations and their macroscopic distribution are simultaneously optimized. Firstly, the PLSM integrated with the numerical homogenization approach is used to topologically optimize multiple prototype microstructures. Then, a shape interpolation method is employed to map these optimized prototype microstructures to generate a series of key microstructures with different volume fractions. Kriging metamodels are constructed based on the key microstructures and used to predict the effective properties of each microstructure with any volume fraction within the macrostructure. Finally, with the predicted effective properties of all microstructures, an optimized material distribution pattern at macroscale is achieved for maximizing the natural frequency of the cellular structure by the VTS method combined with MAC-based mode-tracking strategy. By mapping multiple optimized prototype microstructures, all the microstructures within macrostructure are well connected with each other. The constructed Kriging metamodels for predicting the effective properties of all the microstructures have a great contribution in considerably reducing the computational burden in the multiscale topology optimization. Compared with the traditional monoscale macrostructural and microstructural designs, the Kriging-assisted multiscale topology optimization method enables to fully explore the design space, thereby achieving the cellular structure with higher natural frequency. Numerical examples are provided to demonstrate the validity and advantage of the proposed method. Although the first natural frequency of a cellular structure is maximized in the test examples, the proposed method can be equally extended for the maximization of the natural frequency of other orders.

7. Conclusions

The authors declared that they have no conflicts of interest to this work. We declare that we do not have any commercial or associative interest that represents a conflict of interest in connection with the submitted manuscript entitled ‘‘Maximizing natural fre-

This paper proposes an effective Kriging-assisted multiscale topology optimization method for maximizing natural frequencies

Declaration of Competing Interest

Y. Zhang et al. / Computers and Structures 230 (2020) 106197

quencies of inhomogeneous cellular structures by Kriging-assisted multiscale topology optimization”. Acknowledgements This research was supported by the National Natural Science Foundation of China [grant numbers 51675196 and 51721092], the Natural Science Foundation of Hubei Province [grant number 2019CFA059], and the Program for HUST Academic Frontier Youth Team [grant number 2017QYTD04]. References [1] Bendsøe MP, Kikuchi N. Generating optimal topologies in structural design using a homogenization method. Comput Methods Appl Mech Eng 1988;71:197–224. [2] Bendsøe MP. Optimal shape design as a material distribution problem. Struct Optimization 1989;1:193–202. [3] Bendsøe MP, Sigmund O. Topology optimization: theory, methods, and applications. Berlin, Heidelberg: Springer; 2003. [4] Xie YM, Steven GP. A simple evolutionary procedure for structural optimization. Comput Struct 1993;49:885–96. [5] Da DC, Yvonnet J, Xia L, Li GY. Topology optimization of particle-matrix composites for optimal fracture resistance taking into account interfacial damage. Int J Numer Meth Eng 2018;115:604–26. [6] Wang MY, Wang X, Guo D. A level set method for structural topology optimization. Comput Methods Appl Mech Eng 2003;192:227–46. [7] Allaire G, Jouve F, Toader A-M. Structural optimization using sensitivity analysis and a level-set method. J Comput Phys 2004;194:363–93. [8] Sigmund O. Materials with prescribed constitutive parameters: an inverse homogenization problem. Int J Solids Struct 1994;31:2313–29. [9] Andreassen E, Andreasen CS. How to determine composite material properties using numerical homogenization. Comput Mater Sci 2014;83:488–95. [10] Sigmund O. Tailoring materials with prescribed elastic properties. Mech Mater 1995;20:351–68. [11] Gao J, Li H, Gao L, Xiao M. Topological shape optimization of 3D microstructured materials using energy-based homogenization method. Adv Eng Softw 2018;116:89–102. [12] Da DC, Cui XY, Long K, Huang GX, Li GY. Design of material microstructures for maximum effective elastic modulus and macrostructures. Eng Comput 2018;35:622–40. [13] Da DC, Chen JH, Cui XY, Li GY. Design of materials using hybrid cellular automata. Struct Multidisciplinary Optimization 2017;56:131–7. [14] Zhang Y, Xiao M, Li H, Gao L. Topology optimization of material microstructures using energy-based homogenization method under specified initial material layout. J Mech Sci Technol 2019;33:677–93. [15] Wang Y, Luo Z, Zhang N, Kang Z. Topological shape optimization of microstructural metamaterials using a level set method. Comput Mater Sci 2014;87:178–86. [16] Radman A, Huang X, Xie YM. Topology optimization of functionally graded cellular materials. J Mater Sci 2012;48:1503–10. [17] de Kruijf N, Zhou S, Li Q, Mai Y-W. Topological design of structures and composite materials with multiobjectives. Int J Solids Struct 2007;44:7092–109. [18] Liu L, Yan J, Cheng G. Optimum structure with homogeneous optimum trusslike material. Comput Struct 2008;86:1417–25. [19] Yan X, Huang X, Zha Y, Xie YM. Concurrent topology optimization of structures and their composite microstructures. Comput Struct 2014;133:103–10. [20] Da DC, Cui XY, Long K, Li GY. Concurrent topological design of composite structures and the underlying multi-phase materials. Comput Struct 2017;179:1–14. [21] Da DC, Yvonnet J, Xia L, Le MV, Li GY. Topology optimization of periodic lattice structures taking into account strain gradient. Comput Struct 2018;210:28–40. [22] Chu S, Gao L, Xiao M, Li H. Design of sandwich panels with truss cores using explicit topology optimization. Compos Struct 2019;210:892–905. [23] Gao J, Luo Z, Xia L, Gao L. Concurrent topology optimization of multiscale composite structures in Matlab. Struct Multidiscip Optim 2019;60 (6):2621–51. [24] Xu L, Cheng G. Two-scale concurrent topology optimization with multiple micro materials based on principal stress orientation. Struct Multidiscip Optim 2018;57:2093–107. [25] Sivapuram R, Dunning PD, Kim HA. Simultaneous material and structural optimization by multiscale topology optimization. Struct Multidiscip Optim 2016;54:1267–81. [26] Zhang Y, Xiao M, Li H, Gao L, Chu S. Multiscale concurrent topology optimization for cellular structures with multiple microstructures based on ordered SIMP interpolation. Comput Mater Sci 2018;155:74–91. [27] Gao J, Luo Z, Li H, Gao L. Topology optimization for multiscale design of porous composites with multi-domain microstructures. Comput Methods Appl Mech Eng 2019;344:451–76.

23

[28] Wang YG, Kang Z. Concurrent two-scale topological design of multiple unit cells and structure using combined velocity field level set and density model. Comput Methods Appl Mech Eng 2019;347:340–64. [29] Rodrigues H, Guedes JM, Bendsøe MP. Hierarchical optimization of material and structure. Struct Multidiscip Optim 2002;24:1–10. [30] Wang Y, Chen F, Wang MY. Concurrent design with connectable graded microstructures. Comput Methods Appl Mech Eng 2017;317:84–101. [31] Xia L, Breitkopf P. Concurrent topology optimization design of material and structure within FE2 nonlinear multiscale analysis framework. Comput Methods Appl Mech Eng 2014;278:524–42. [32] Zhang Y, Xiao M, Zhang X, Gao L. Topological design of sandwich structures with graded cellular cores by multiscale optimization. Comput Methods Appl Mech Eng 2019. https://doi.org/10.1016/j.cma.2019.112749. [33] Zhang Y, Xiao M, Gao L, Gao J, Li H. Multiscale topology optimization for minimizing frequency responses of cellular composites with connectable graded microstructures. Mech Syst Signal Process 2020;135(106369). https:// doi.org/10.1016/j.ymssp.2019.106369. [34] Zhang P, Toman J, Yu Y, Biyikli E, Kirca M, Chmielus M, et al. Efficient designoptimization of variable-density hexagonal cellular structure by additive manufacturing: theory and validation. J Manuf Sci Eng 2014;137. 021004. [35] Cheng L, Bai JX, To AC. Functionally graded lattice structure topology optimization for the design of additive manufactured components with stress constraints. Comput Methods Appl Mech Eng 2019;344:334–59. [36] Wang X, Zhang P, Ludwick S, Belski E, To AC. Natural frequency optimization of 3D printed variable-density honeycomb structure via a homogenization-based approach. Addit Manuf 2018;20:189–98. [37] Cheng L, Liu JK, Liang X, To AC. Coupling lattice structure topology optimization with design-dependent feature evolution for additive manufactured heat conduction design. Comput Methods Appl Mech Eng 2018;332:408–39. [38] Takezawa A, Zhang X, Kato M, Kitamura M. Method to optimize an additivelymanufactured functionally-graded lattice structure for effective liquid cooling. Addit Manuf 2019;28:285–98. [39] Yoon GH. Structural topology optimization for frequency response problem using model reduction schemes. Comput Methods Appl Mech Eng 2010;199:1744–63. [40] Zhao J, Wang C. Topology optimization for minimizing the maximum dynamic response in the time domain using aggregation functional method. Comput Struct 2017;190:41–60. [41] Pedersen NL. Maximization of eigenvalues using topology optimization. Struct Multidiscip Optim 2000;20:2–11. [42] Huang X, Zuo ZH, Xie YM. Evolutionary topological optimization of vibrating continuum structures for natural frequencies. Comput Struct 2010;88:357–64. [43] Yoon GH. Maximizing the fundamental eigenfrequency of geometrically nonlinear structures by topology optimization based on element connectivity parameterization. Comput Struct 2010;88:120–33. [44] Liu Q, Ruan D, Huang X. Topology optimization of viscoelastic materials on damping and frequency of macrostructures. Comput Methods Appl Mech Eng 2018;337:305–23. [45] Yun KS, Youn SK. Microstructural topology optimization of viscoelastic materials of damped structures subjected to dynamic loads. Int J Solids Struct 2018. [46] Andreassen E, Jensen JS. Topology optimization of periodic microstructures for enhanced dynamic properties of viscoelastic composite materials. Struct Multidiscip Optim 2014;49:695–705. [47] Niu B, Yan J, Cheng G. Optimum structure with homogeneous optimum cellular material for maximum fundamental frequency. Struct Multidiscip Optim 2008;39:115–32. [48] Zuo ZH, Huang X, Rong JH, Xie YM. Multi-scale design of composite materials and structures for maximum natural frequencies. Mater Des 2013;51:1023–34. [49] Liu Q, Chan R, Huang X. Concurrent topology optimization of macrostructures and material microstructures for natural frequency. Mater Des 2016;106:380–90. [50] Xiao M, Zhang J, Gao L, Lee S, Eshghi AT. An efficient Kriging-based subset simulation method for hybrid reliability analysis under random and interval variables with small failure probability. Struct Multidiscip Optim 2019;59:2077–92. [51] Zhang J, Xiao M, Gao L, Chu S. Probability and interval hybrid reliability analysis based on adaptive local approximation of projection outlines using support vector machine. Comput-Aided Civ Infrastruct Eng 2019;34:991–1009. [52] Yoshimura M, Shimoyama K, Misaka T, Obayashi S. Topology optimization of fluid problems using genetic algorithm assisted by the Kriging model. Int J Numer Meth Eng 2017;109:514–32. [53] Zhang J, Xiao M, Gao L, Fu J. A novel projection outline based active learning method and its combination with Kriging metamodel for hybrid reliability analysis with random and interval variables. Comput Methods Appl Mech Eng 2018;341:32–52. [54] Zhang J, Xiao M, Gao L, Chu S. A combined projection-outline-based active learning Kriging and adaptive importance sampling method for hybrid reliability analysis with small failure probabilities. Comput Methods Appl Mech Eng 2019;344:13–33. [55] Seyranian AP, Lund E, Olhoff N. Multiple eigenvalues in structural optimization problems. Struct Optimization 1994;8:207–27.

24

Y. Zhang et al. / Computers and Structures 230 (2020) 106197

[56] Kim TS, Kim YY. Mac-based mode-tracking in structural topology optimization. Comput Struct 2000;74:375–83. [57] Chu S, Gao L, Xiao M, Luo Z, Li H, Gui X. A new method based on adaptive volume constraint and stress penalty for stress-constrained topology optimization. Struct Multidiscip Optim 2018;57:1163–85. [58] Chu S, Gao L, Xiao M, Zhang Y. Multiscale topology optimization for coated structures with multifarious-microstructural infill. Struct Multidiscip Optim 2019. https://doi.org/10.1007/s00158-019-02428-y.

[59] Luo Z, Tong L, Kang Z. A level set method for structural shape and topology optimization using radial basis functions. Comput Struct 2009;87:425–34. [60] Svanberg K. The method of moving asymptotes—a new method for structural optimization. Int J Numer Meth Eng 1987;24:359–73. [61] Cadman JE, Zhou S, Chen Y, Li Q. On design of multi-functional microstructural materials. J Mater Sci 2013;48:51–66.