Multiscale Isogeometric Topology Optimization with Unified
Structural Skeleton

This paper proposes a multiscale isogeometric topology optimization (ITO) method where the configuration and layout of microstructures are optimized simultaneously. At micro scale, a shape deformation method is presented to transform a prototype microstructure (PM) for obtaining a series of graded microstructures (GMs), where microstructural skeleton based on the level set framework is applied to retain more topology features and improve the connectability. For the macro scale calculation, the effective mechanical properties can be estimated by means of the numerical homogenization method. By adopting identical non-uniform rational basis splines (NURBS) as basis functions for both parameterized level set model and isogeometric calculation model, the isogeometric analysis (IGA) is integrated into the level set method, which contributes to improving the accuracy and efficiency. Numerical examples demonstrate that, the proposed method is effective in improving the performance and manufacturability.


Introduction
Research about artificial microstructures has already attracted great attention in engineering for the extraordinary properties [Chen, Ortiz and Huang (1998); Gibson and Ashby (1997)], such as impact energy absorption, thermal insulation. Inspired by the natural materials, cellular solids as promising materials have experienced an increasing study due to the superior performance [Valdevit, Jacobsen, Greer et al. (2011);Zheng, Lee, Weisgraber et al. (2014)]. In recent years, the advances of additive manufacturing technology have promoted the relative eases of fabrication for the complex structures. Thus, the cellular materials with periodically arranged microstructures are particularly brought into focus. In general, the constituent of periodic cellular materials is ordinary material such as mental.
Therefore, the key issue for gaining the superior performance of such materials is to design the configuration of the cellular microstructures at micro scale and assemble the cellular microstructures at macro scale, instead of adjusting the material composition. In spite of the large amounts of researches including analytical, numerical or experimental method devoted to cellular composite design, the design method for searching the optimal configuration and layout of cellular microstructures is still insufficient. Topology optimization (TO) is a mathematical method to calculate the optimal layout of materials within a prescribed design domain under boundary conditions. In recent years, a variety of methods have been proposed and developed to solve TO problems, such as the density-based homogenization methods [Bendsoe and Kikuchi (1988); Suzuki and Kikuchi (1991)], solid isotropic material with penalization (SIMP) methods [Zhou and Rozvany (1991); Bendsoe and Sigmund (1999)], evolutionary structural optimization (ESO) methods [Xie and Steven (1993) ;Huang, Radman and Xie (2011);Huang, Xie, Jia et al. (2012)], level set methods (LSMs) [Sethian and Wiegmann (2000); Wang, Wang and Guo (2003); Allaire, Jouve and Toader (2004)], and moving morphable components (MMC) [Guo, Zhang and Zhong (2014); Guo, Zhang, Zhang et al. (2016)]. Reference [Sigmund and Maute (2013)] gives a comparison and critical review to the various approaches. With respect to explicit material representation methods, the LSM employs an implicit high-dimensional function to describe the geometry. In the LSM, boundaries are represented as the zero contours of the level set function, by which the variation of topology structure can be achieved. Due to the characteristics like smooth boundaries and distinct interfaces [van Dijk, Maute, Langelaar et al. (2013)], the LSM has advantages on topology optimization structure design and fabrication. In conventional LSM, the evolution of structural boundary with respect to time is formulated as a Hamilton-Jacobi partial differential equation (PDE). The materials and shape derivatives from shape variational analysis [Wang, Wang and Guo (2003)] are applied to solve the problem. However, complicated numerical processes such as reinitializing the LSF to retain the feature of signed-distance, extending velocity field of boundary to the design domain decrease the computational efficiency. To eliminate these undesired numerical processes, the parameterized level set methods (PLSMs) which employed a variety of basis functions including finite element (FEM) basis functions [van Dijk, Langelaar and Keulen (2012)], radial basis functions (RBFs) [Wang and Wang (2006)] and compactly supported RBFs (CSRBFs) [Wendland (2006)] were proposed and developed. PLSM can avoid such redundant processes and combine the well-established optimization algorithms with their characteristics as well as inherit the advantages of LSM. For the past few years, isogeometric analysis (IGA) [Hughes, Cottrell and Bazilevs (2005) ;Cottrell, Hughes and Bazilevs (2009)], whose basis functions for analysis and computer aided design (CAD) processes are unified, has become an efficient alternative to conventional FEM in TO problems due to the advantages such as high continuity and high accuracy. The first study on isogeometric TO traces back to 2010, Seo et al. [Seo, Kim and Youn (2010)] employed the trimmed surface technique [Kim, Seo and Youn (2009)] to analysis structural response and calculate the sensitivity. Later, the optimality criteria (OC) and the method of moving asymptotes (MMA) were applied to isogeometric SIMP referred to Hassani et al. [Hassani, Khanzadi and Tavakkoli (2011); Tavakkoli, Hassani and Ghasemnejad (2013)]. Kumar et al. [Kumar and Parthasarathy (2011)] utilized Bspline finite elements to develop density-based topology optimization for avoiding the checkboard artifacts. Dede et al. [Dede, Borden and Hughes (2012)] applied a phasefield model to solve TO problem, and applied IGA to represent the design domain with the unique geometric exactness. Qian [Qian (2013)] proposed a density-based TO in Bspline space where the density field and the density for analysis are represented by tensorproduct B-spline. Thereafter, Wang et al. [Wang and Benson (2015)] presented an accurate and efficient ITO where the NURBS based IGA integrates with the PLSM by means of utilizing the NURBS basis functions of the CAD models for both the parameterization of level set function and the evaluation of structural analysis. Moreover, Hou et al. [Hou, Gai, Zhu et al. (2017)] applied MMC to explicit isogeometric topology optimization (ITO) to improve stability and robustness for the reason that the analysis calculation is replaced by IGA. ] presented an effective and efficient isogeometric topology design method by the adoption of enhanced density distribution function. Furthermore, for the systematic design of auxetic metamaterials, an isogeometric topology optimization method was proposed in Gao et al. [Gao, Xue, Gao et al. (2019)]. Wang et al. [Wang, Wang, Xia et al. (2018)] presented a comprehensive review on ITO with comparisons to traditional FEM design methods modestly. For improving the computational efficiency, Wang et al. [Wang, Liao, Ye et al. (2020)] proposed a highefficiency ITO with the collaboration of multilevel mesh, MGCG and local-update strategy. Referring to many applications of TO [van Dijk, Maute, Langelaar et al. (2013)], the design for engineering cellular materials attracted increasing attentions. Since the inverse homogenization method [Sigmund (1994)] has been adopted to predict effective properties of cellular microstructures, the design of microstructures is able to be formulated as an inverse problem for yielding the expected macro performance [Andreassen, Lazarov and Sigmund (2014); Takezawa, Kobashi and Kitamura (2015)] and prescribed properties [Xie, Yang, Shen et al. (2014)]. An extended design strategy of optimizing the functionally graded microstructures (GMs) is considered as alternative to offer more mechanical advantages for external stimuli. For the reason that both microstructures and macrostructures have effect on the performance of cellular composites, the multiscale concurrent design means that the microstructures are expected to be designed as well as assembled. Taking into account the linear [Fujii, Chen and Kikuchi (2001)] or nonlinear [Kato, Kyoya, Yachi et al. (2014)] macroscopic behavior, a periodic arranged microstructure was designed with the predefined distribution. In general, each element corresponds to an individual microstructure during multiscale optimization process, which results in varied optimized microstructures corresponding to different amounts of elements. With respect to such problem, Zhang et al. [Zhang and Sun (2006)] and Xie et al. [Xie, Zuo, Huang et al. (2012)] researched the scale effect and drew a conclusion that the configurations of microstructures would converge to a certain result along with the decrease of scale ratios. A reasonable design approach for hierarchical composite structure is supposed to optimize the configuration of cellular microstructures and assemble the cellular microstructures simultaneously. Therefore, Coelho et al. [Coelho, Fernandes, Guedes et al. (2008)] proposed a hierarchical optimization framework with double loops which correspond to macroscale materials distribution design and microscale unit cell structure design respectively. Nakshatrala et al. [Nakshatrala, Tortorelli and Nakshatrala (2013)] proposed another framework aimed to couple the nonlinear problems between macro and micro topology optimizations. Wang et al. [Wang, Xu and Pasini (2016)] presented an isogeometric topology optimization for the periodic cellular microstructures where the effective properties of cellular microstructure are formulated as functions with respect to the relative density. Referring to Wang et al. [Wang, Chen and Wang (2016)], a shape metamorphosis method was presented for obtaining a family of GMs with similar topology, thereby establishing the relation between effective properties and densities to achieve the multiscale concurrent design. Gao et al. [Gao, Luo, Xia et al. (2019)] gave the compact and efficient Matlab codes in both 2D and 3D for the concurrent TO of multiscale composite structures, in which the energy-based homogenization method was used to calculate the macroscale effective properties of microstructures. With a view to the high accuracy and efficiency of IGA and the distinct interfaces of LSM, we present a multiscale isogeometric topology optimization method where the configuration and layout of microstructures are optimized simultaneously. For the multiscale scheme, we apply the material distribution-based method [Cheng, Liu and Yan (2008)] to distribute cellular materials in macro scale, meanwhile, optimize the configuration of microstructures in micro scale by means of the isogeometric parameterized LSM combined with numerical homogenization method. Structural skeleton based on the level set framework [Xia and Shi (2015)] is employed to improve the topology similarity and connectability among GMs. The remainder of this paper is organized as follows: In Section 2, a shape deformation method is presented to yield a series of microstructures with the unified structural skeleton and similar topology. Section 3 briefly reviews NURBS-based IGA and introduces the integration of PLSM and IGA based on NURBS parametrization. The formulation of multiscale concurrent optimization and the sensitivity analysis are presented in Section 4. Thereafter in Section 5, some numerical examples are presented to demonstrate the effectiveness in improving the performance and manufacturability. Finally, a brief conclusion is given in Section 6.

Generation of microstructures
In this section, a series of microstructures with the similar topology are generated by means of a shape deformation technology. Therefore, these microstructures share a same prototype microstructure (PM) expressed by the level set function, where Ω pm denotes the structural regions, and D is the design domain. The structural boundary ∂Ω pm is represented as zero iso-contour of the implicit high-dimensional function Φ pm (x, t).
In order to obtain a series of GMs ranged from void to solid with similar structure to ensure the connectivity to each other, a shape deformation technique is proposed in this paper. Referring to Wang et al. [Wang, Chen and Wang (2016)], the GMs were generated along with parallelly moving the iso-contour of level set function as illustrated at the left side of Fig. 1. However, the GMs interpolated between void and the PM may yield structural fracture when structural component is too slender. To eliminate these unexpected fracture, the structural skeleton based on the level set framework is utilized to maintain the shape of microstructural components.
Prototype micorstructure with skeleton Figure 1: Schematic illustration of shape deformation technique. Left column are the GMs obtained by parallelly moving the iso-contour of level set function. Relatively, the GMs of right column are generated by scaling the components size of prototype microstructure based on the skeleton (marked as red line) As showed in Fig. 1, we can see the effects brought by the introduction of structural skeleton intuitively. Referring to Xia et al. [Xia and Shi (2015)], for the cellular materials with periodically arranged microstructures, the microstructural skeleton can be yielded as the marked red line of the figure. Taking the skeleton as a benchmark, level set values are recalculated by zooming in or out within cell domain to generate a series of microstructures with holding similar basic topology shape. Comparing the GMs of the left and right columns in Fig. 1, the microstructures with skeleton retain more topology features by means of adjusting the proportion of materials around skeleton when the structural components are too slender (when the relative densities in Fig. 1 is below 0.3).
The presented shape deformation technology works as scaling the components size of prototype microstructure corresponding to the different relative densities. Thus the GMs connected to each other naturally and continuously in theory. Because of the TO in macro scale just focusing on prototype microstructure, the weird-looking level set values of GMs do not have to be used for calculation.

Effective properties of graded microstructures
Assuming the microstructures are periodically arranged within the design domain which is much larger than the microstructure. The effective properties of microstructures can be estimated by means of the numerical homogenization method [Sigmund (1994)]. By imposing the periodic boundary conditions and referring to the small parameter perturbation, the effective elasticity tensor of a microstructure can be yielded as, where V N e and D klmn denote the volume of microstructure and the elasticity of base material respectively. For the 2D elasticity, the unit vectors (1, 0, 0) T , (0, 1, 0) T and (0, 0, 1) T are assembled to constitute strain field ε 0 kl . The period displacement field u pq and the strain field ε * kl (u pq ) produced by ε 0 kl can be solved in the following equation, whereŪ (V N e ) denotes the kinematically admissible displacement space. According to Zhou et al. [Zhou and Li (2008)], when the differences of spatial topology and effective property between neighboring microstructures are sufficiently small, the effective properties of microstructures can be applied to analysis the displacement field at macro scale. With regard to the GMs obtained in this paper, they possess similar topology structure and connect to each other naturally and continuously in theory. Exploring the effective properties of GMs mentioned in Fig. 1, the elasticity tensors of a few GMs corresponding to different relative densities have been calculated by Eq.
(2) and interpolated to estimate the properties of whole densities field as showed in Fig. 2. In order to compare the effective properties of microstructures to the property of base material, the regularized operations are implemented in this paper. In consideration of the zero element in the elasticity tensor of base material, the effective elasticity tensors are regularized by diving the factor E scale = E base /(1 − µ 2 ) which can also be written as D base (1, 1). Here E base and µ denote the Young's modulus and Poisson's radio of basis material respectively. In terms of Wang et al. [Wang, Chen and Wang (2016)], the number of symbol GMs selected as 21 will balance the accuracy of estimation and the computational efficiency. When the large-scale GMs are well-organized arranged in the design domain, the macroscopic displacement field can be obtained and utilized to analysis the microstructures.

Integration of PLSM and IGA based on NURBS parameterization
With a view to the smooth boundaries and distinct interfaces of LSM, the high continuity and accuracy of IGA, a method which integrates the LSM and IGA just meets the demands

Relative densities of GMs
Regularized effective elasticity tensors Figure 2: Example to illustrate the regularized effective elasticity tensors with respect to the relative densities. D(i, j) is the element of regularized elasticity matrix of the formulation for GMs mentioned in Section 2. As an evolution of conventional LSM, the parameterized LSM provides the bridge to combine with IGA by means of unified basis function for parametrization.

Review of NURBS-based IGA
For the IGA method, multiple NURBS patches are assembled to compose the geometry with complex topology and shape. Compared to the conventional FEM, the geometry representation and analysis of IGA are based on the same model. Beginning with the B-spline, a knot vector constituted by a sequence of non-decreasing real numbers is represented as: where p, n denote the order of B-spline and the number of basis functions respectively. Referring to the Cox-de Boor formula [De Boor (1972)], the B-spline basis functions are recursively generated as: Then, NURBS basis functions can be yielded by introducing the positive weight value ω i , Analogously, on the basis of the tensor product formulation, the 2D NURBS basis functions are represented as: where B i,p (ξ), B j,q (η) are B-spline basis functions constructed from the knot vectors Ξ= {ξ 1 , ξ 2 , · · · , ξ n+p+1 } and H= {η 1 , η 2 , · · · , η m+q+1 }. A NURBS surface can be represented as a tensor product of bivariate NURBS with p degree in ξ direction and q degree in η direction. where

Parameterized LSM based on NURBS basis function
To eliminate the undesired numerical processes such as reinitializing the LSF to retain the feature of signed-distance, extending velocity field for solving the PDE, the parameterized level set methods (PLSMs) are proposed. With the interpolation of basis functions, the Hamilton-Jacobi PDE is converted into ordinary differential equation (ODE) thereby simplifying the solving process. The parameterized level set function (LSF) with basis function is formulated as a summation.
where N (x) = [N 1 (x), N 2 (x), · · · , N n (x)] T is the basis function vector at the coordinate x. Correspondingly, A(t) = [α 1 (t), α 2 (t), · · · , α n (t)] T denotes the expansion coefficient at the time t. By the introduction of basis function, the representation of LSF is divided into two separate parts which are only associated with space and time respectively. Thus, the evolution of optimization is just performed by updating the expansion coefficient Substituting Eq. (9) to Eq. (10), it is converted into the form written as, thereby the normal velocity υ n = − ∂x ∂t ∇Φ |∇Φ| can be obtained as, Compared to the parametrization with basis functions such as finite element (FEM) basis functions [van Dijk, Langelaar and Keulen (2012)], radial basis functions (RBFs) [Wang and Wang (2006)] and compactly supported RBFs (CSRBFs) [Wendland (2006)], the utilization of NURBS [Wang and Benson (2015)] makes interpolation points not limited in the design domain no longer, that meet the demand of the isogeometric interpolation.
As one dimension problem with the parametric coordinate , the NURBS basis function is written as N i (ξ). Eq. (9) is expressed in the form as, In general, TO starts from an initial geometry with multiple holes evenly distributed in design domain. For the LSM, the initial level set values can be obtained by signed-distance function. The focus problem is to evaluate the initial expansion coefficients corresponded to the control points whose number is assumed as m. There will be m collocation points as well as their equations need to be established. The varied selection strategies of collocation points were presented such as the Gaussian quadrature scheme [De Boor (1973)] and the Greville abscissae scheme [Johnson (2005)]. On the basis of comparison in Qian [Qian (2011)], the Greville abscissae is favored for the advantages as the stability and accuracy during IGA. Greville abscissae are defined as, where ξ i denotes the ith knot of the vector Ξ= {ξ 1 , ξ 2 , · · · , ξ m+p+1 } from the NURBS which are p order and possess m control points. With respect to the collocation points located in the NURBS surface, two Greville abscissae for the coordinate expression in ξ and η direction are established and evaluated respectively. By means of substituting into Eq. (8), the physical coordinates corresponded to the collocation points are obtained. Afterwards, the corresponding level set values can be calculated. The linear equation set is established in a similar way to Eq. (9), where Φ denotes the vector containing all initial level set values of collocation points. As a matrix, A consists of NURBS basis functions corresponded to the collocation points. Similar to Eq. (9), the updates of expansion coefficients α are only related to time rather than the coordinates.

Multiscale isogeometric topology optimization
In this section, taking the compliance minimization problem as an example, the graded microstructures (GMs) based multiscale TO is presented and performed in a double-loop strategy. Although the macroscopic displacement U = U (X) is decoupled with the microscopic displacement u = u (x), the reasonable solution of optimization can be obtained in this decouple design scheme according to researches concerning in nonlinear coupling scheme [Xia and Breitkopf (2015); Alexandersen and Lazarov (2015)].

Formulation of multiscale design with GMs
The proposed multiscale optimization method is divided into two parts, the layout of GMs at macro scale and the configuration of prototype microstructure at micro scale.
Here, we utilize the density-based method to assemble cellular material due to its stability and the ability to produce the intermediate density element. Meanwhile, explore the configuration of microstructures by means of the isogeometric parameterized LSM. The numerical homogenization method is the bridge to establish the numerical relation between two scales.

Macroscopic topology optimization for the distribution of GMs
In order to generate the layout of GMs with material volume constraint, the formulation of TO based on material distribution method [Cheng, Liu and Yan (2008)] can be represented as, where the densities Θ of all the microstructure element are design variables and N is the total number of microstructures. G is the volume constraint for the global design domain when V 0 and V max denote the volume of element and the maximum volume in global respectively. J is the structural mean compliance calculated by the external load F and corresponding displacement field U . It is noted that the elasticity tensors D Ne of GMs are represented as a function of the element density ρ Ne . As discussed in Section 2.2, the GMs obtained in proposed method possess similar topology structure and connect to each other naturally and continuously in theory. Therefore, the mechanical properties of immediate density element can be estimated by interpolation, where ρ i , ρ i+1 are densities of the adjacent symbol GMs. Correspondingly, D i , D i+1 are the elasticity tensors which can be calculated by Eq.
(2). The density ρ min set in Eq. (16) is the density lower bound for the design variables. Mentioned in Section 2, although the microstructures with skeleton retain more topology features by means of adjusting the proportion of materials around the skeleton, structural fracture still occur inevitably when the relative density is too small. For avoiding the structural fractures of GMs, a density lower bound is set to constrain the sizes of microstructural components. The lower bound is not only depended on microstructural skeleton but also influenced by the precision of microstructural model to some extent. Considering the simplification of formulation for different structural skeletons and calculation efficiency related to DOF, ρ min = 0.1 is selected referred to the results of experiments.

Microscopic topology optimization with isogeometric PLSM
As mentioned above, the mechanical properties of immediate density microstructures can be estimated by Eq. (17) which is only influenced by the macroscopic design variable Θ and topology of the prototype microstructure (PM). Hence, the problem in micro scale focuses on the structure design of PM. Discretizing the PM into elements, the formulation of microscopic TO based on the isogeometric parameterized LSM can be represented as, where the expansion coefficient A P M of the PM are design variables and n is the total number of expansion coefficients. g P M is the volume constraint of the PM at unit cell when V P M and V P M,max denote the actual volume and the maximum volume of the PM respectively. J is the summation of all microstructural mean compliance calculated by the stiffness matrix k Ne of microstructure and the corresponding displacement field u Ne . It is noted that the elasticity tensors D H of GMs are represented as a function of the element density ρ Ne and the expansion coefficient A P M . At this point, ρ Ne is invariable during a microscopic optimization of double loop scheme. Therefore, the core of affecting the objective function J concentrate on the elasticity tensors D H P M of PM, which can be estimated by Eq. (2). From the perspective of the whole optimization process, macro design variables ρ Ne and micro design variables A P M are updated alternately by the double loop scheme. In this way, several merits are able to be gained, such as: (1) Integrated optimization problem is decoupled into the usual TO in single scale, where the more constrains can be concerned as well as the high efficiency of calculation: (2) The density based optimization in macroscopic level impels the densities of GMs to be distributed smoothly and continuously in the whole design domain; (3) With the help of the advantages as smooth boundaries and distinct interfaces of the LSM, the topology structure of GMs can be expressed accurately that benefits to the estimation of mechanical properties.

Sensitivity analysis
In the proposed design method, the problem will be settled by the gradient based algorithms. The optimality criteria (OC) [Rozvany and Kirsch (1995)] is introduced to solve the distribution of GMs, meanwhile the method of moving asymptotes (MMA) [Svanberg (1987)] is utilized to update the configuration of the PM. The first-order derivatives of the objective function and constraints with respect to the design variables are required in both the algorithms.

Sensitivity analysis for the distribution of GMs
At macro-scale, the derivatives of objective function in Eq. (16) with respect to the densities of GMs can be obtained as, As discussed in Section 2, when the prototype microstructure (PM) is certain, the topology structure and corresponding mechanical properties of the GMs are only affected by the relative densities. Considering the density lower bound, a projection method is adopted to allow for void microstructure. The elasticity tensors of GMs are calculated by is the Heaviside step function and D void = ρ min D denotes the pseudo elasticity tensor of the void regions. Thus, the derivatives of elasticity tensors can be approximately calculated by combining the interpolation function Eq. (17) from Fig. 2. Additionally, the derivatives of the volume constraint are calculated as,

Sensitivity analysis of the prototype microstructure
It is noted that all the GMs distributed within design domain possess the independent densities in theory. Therefore, the objective function represented as the summation of all microstructures in Eq. (18) requires the huge amount of computation. In order improve the efficiency of the design method, the objective function is simplified as, where (N 1 , N 2 , · · · , N m ) denotes the graded microstructure with the similar relative density ρ N which accounts for a considerable proportion in the global. Hence, the derivatives of Eq. (21) with respect to the expansion coefficients can be written as, The key issue focuses on the derivatives of the elasticity tensors with respect to the expansion coefficients. Referring to the shape derivative in Sokolowski et al. [Sokolowski and Zolesio (1992)], the derivatives with respect to t can be obtained as, is the Dirac function and Eq. (12) into Eq. (23), it is converted into, where N ne (x) is the NURBS basis function corresponding to the n e th collocation point and expansion coefficient, Considering the chain rules, the derivatives of the elasticity tensors with respect to the expansion coefficients will yield, In the same way, the derivatives of volume constraint can be calculated by

IGA for prototype microstructure
With respect to the equilibrium equations required in Eq. (18) during the microstructure optimization, the stiffness matrix K is composed of all the element stiffness matrix K e which is represented as, where B denotes the strain-displacement matrix while D is the elasticity tensors of material. Ω e is the physical domain of element, correspondingly, Ω e and Ω e are the domains in the NURBS parametric space and the integrated parametric space respectively. The Jacobian J 2 and J 1 as the bonds are utilized to build the relation from the integrated parametric space to the NURBS parametric space, then to the physical space. For the plane stress problem, B is represented as a matrix, where N i denotes the basis function of NURBS element corresponded to nc control points. The Jacobian J 1 is written as, Transforming the Gauss quadrature domain to the NURBS parametric domain [ξ i , ξ i+1 ] × [η j , η j+1 ], the relation can be established by the equation Thus, the Jacobian J 2 can be yield as, For the reason that the design domain of microstructure is constant, it is not necessary to remesh the element during optimization. In consideration of the singularity in analysis calculations, the "ersatz material" is employed. Thereby, the element stiffness matrix of the intermediate density can be expressed as a multiple of entity element stiffness matrix. (2) the GMs with unified structural skeleton are just generated by the PM and the relative densities Θ;

Numerical implementation
(3) the relation between the elasticity tensors yielded from Eq. (2) and the relative densities can be established as Fig. 2. On this basis, the layout of GMs and the configuration of PM are alternately updated to a locally optimal solution. Finally, the double-loop optimization will terminate when the number of total loop reaches the predefined maximum. It is noted that the generation of GMs and the corresponding estimation of mechanical properties are separate from each other. Therefore, the parallel computing strategy can be utilized to improve the efficiency. The same goes for the sensitivity analysis in Eq. (22) and Eq. (26) for the reason that the variables of the integration are almost same with respect to the different derivatives. The OC is introduced to solve the distribution of GMs, meanwhile the MMA is utilized to update the configuration of PM.

Numerical examples
In this section, numerical examples for the multiscale minimal compliance optimization problems are presented and illustrate the effectiveness in improving the performance and manufacturability. The programs are run on the desktop where the CPU is Intel core i7-2600 K of 3.4 GHz, the RAM is 12 GB, and the software environment is MATLAB 2015a. In macro scale, unified 4-node linear elements are adopted to represent the GMs one by one. Meanwhile in micro scale, the quadratic isogeometric elements with the 3 × 3 Gauss quadrature rule are established. All the microstructures consist of a base material whose Young's modulus and Poisson's radio are 1000 and 0.3 respectively. It is noted that the regularized effective elasticity tensors are obtained by diving the factor For the alternative one-scale optimization, the loop will terminate when the relative difference of the design variables between two iterations is under 10 −3 . The distribution of GMs and configuration of PM as the optimized solution are outputted as soon as the number of double-loop reaches the predefined maximum.

Layered beam design
As shown in Fig. 4, we test a beam problem which is fixed on the left edge. At the top edge, a linearly-increasing distributed load with the maximum value q max = 1.6 is applied. The length and height of macrostructure are 32 and 20 respectively when the thickness is set as 1.
The beam is discretized into 20 layers with 20 × 32 four-node elements, where each layer is composed by the identical GMs. Thus, the macro design variables ρ N e of each layers will keep consistent. 32 × 32 quadratic isogeometric elements constitutes the mesh of microstructure for generating the GMs and solving the corresponding displacement fields. The maximum volume of PM is set as V P M, max = 0.4 while the relative densities ρ N e of GMs are allowed to vary within the interval [0.1, 1]. A material usage constraint of 50% is applied for the macrostructure. to the GMs of the relative densities 0.28, 0.35 and 0.5 respectively. It is found that the upper and lower layers are fully-solids which will provide sufficient stiffness for the avoiding of bending deformation something like the solid face-sheets in sandwich. For the reason that the GMs of middle layers are generated from the unified PM as shown in the right side of Fig. 6, they possess identical structural skeleton and similar topology structure that benefits to the connectivity and manufacturability. In the theory, the connection among GMs will realize completely continuous with the infinite elements. The optimized solution of the PM under volume constrain 0.4 is plotted in Fig. 6 where the red line is marked as the structural skeleton. In order to support the linearly increasing loads rightwards, the PM contains the upper-right-aligned structural components to enhance the directional stiffness that will be mapped to all the GMs. The mechanical properties exhibited below the structure figure also illustrate the sense of the components distribution. Here, the elasticity tensors of the GMs with relative densities (2) and linearly interpolated as the chart in Fig. 6. It can be found that the GMs obtained from the proposed shape deformation technology almost possess similar effective properties and behave as a monotone function of the densities. Hence, the configuration of PM is representative for all microstructures and it is expected to generate the preferable solution.
As the double-loop scheme proceeding, the macro distribution design will convergence to a stable locally optimal solution fast with OC algorithm at a single scale. Here, the stages of the micro configuration design are gathered as shown in Fig. 7. It is seen that there are several dramatic changes around Stage 35, 50 and 70 of iterative curve. Comparing the configurations of prototype microstructures in corresponding Stages, the structural fractures as well as the holes mergence gradually occurred in some area. The GMs based on unified structural skeleton will be in the same situation. Thus, the temporary weakening of microstructural effective properties induced the sharp increase in the objective function. In preparation, the initial microstructure with the multiple holes is utilized to generate the GMs. Then one of them corresponded to V P M, max is selected as the PM, and the GMs corresponded to the material usage constraint are pre-distributed within the design domain. Thus, the volume ratio in Fig. 7 keep constant. The objective function of the final optimized solution is J GM = 14.32.
To illustrate the advantages, the identical problem is solved by the SIMP method in single macro scale. As shown in Fig. 8(a), the materials with the relative densities varying in interval [0.4, 1] are distributed in a similar pattern compared to the proposed method. The objective function finally convergence to J SIM P = 32.21, significantly greater than the proposed approach. It can be considered as that the multiscale design gives birth to broad design space, and the introduction of micro scale design releases the design freedoms to advance a remarkable improvement. Fig. 8 Figure 7: Objective function in the micro configuration design, and the PM of key iteration stage compared to the SIMP, the performance is inferior to the multiscale design with the proposed shape deformation technology. As discussed in Section 2, the microstructures with skeleton retain more topology features, thus the low-density GMs are able to be utilized with good connectivity. Therefore, the more plentiful candidates for GMs release the implicit design constrains modestly.

Cantilever beam
This example is to evaluate the result of the well-studied cantilever beam optimization problem as illustrated in Fig. 9. Here, a beam with 120 of length and 30 of height is fixed at the left edge, and a vertical concentrated force F = 5 is applied at the center of the right boundary. A mesh constituted by 32 × 32 quadratic isogeometric elements is adopted to generate the GMs and analysis the corresponding displacement fields. The maximum volume of PM is set as V P M, max = 0.5 while the relative densities ρ N e of GMs can vary between 0.1 and 1. A material usage constraint of 50% is applied for the macrostructure. As the material usage constraint, the 30 × 120 four-node elements in macro scale are set as the "gray materials" with the relative density of 0.5. On the foundation of such initial conditions, all the microstructures are divided to five groups as shown in Fig. 10 according to the orientations of major principal stresses. Each group corresponds to the GMs generated from the same PM. Considering the symmetry of the problem, the design is transformed to three parallel microscale designs which share the unified objective of optimization.
The final optimized solution is plotted in Fig. 11, as the microstructures defined by the PMs a b  Fig. 12 illustrates the configuration of three PMs, as well as the elasticity tensors of their corresponding GMs. Compared to the well-known solution of SIMP, there are a mass of the fullsolid microstructures distributed in a similar way. It can be considered for the reason that as the stiffest microstructure of GMs, the solid material should be placed in the mainly structural frame for supporting the load. The intermediate-densities GMs are distributed around the mainly components for stabilizing and strengthening stiffness for the efficient utilization. As the discussed characteristic of the identical structural skeleton, the advantages of connectivity and manufacturability are still maintained. The optimized compliance convergences to J GM = 4.356 in the end. Fig. 12 illustrates the configuration of three PMs with the marked structural skeletons. The elasticity tensors of the corresponding GMs are calculated and interpolated as the charts respectively. Although the trend of the interpolation curves belonged to the first PM behaves similar to the example 1, the interpolation curves in Figs. 12(b) and 12(c) take a sudden turn near ρ GM = 0.9. From the perspective of the structural shape and elasticity tensors of corresponding PMs, it can be aware that they appear the much weak bearing performance in vertical direction relative to horizontal direction even close to none. The GMs yielded from such PMs also maintain the directional features until the change of structural frame occurring caused by the relative densities close to 1. As shown in Fig. 11, the GMs yielded from the third PM mainly appear as a solid form or sandwiched between entity components that moderate the extreme bearing weakness. For comparison, the identical problem is solved by means of SIMP and the unified microstructure optimization. The optimized mean compliances are J SIM P = 4.437 and J homo_micro = 9.835 respectively. It once again confirms that the proposed multiscale design method is effective for improving structural stiffness compared with the single scale design. With respect to the little progress effect relative to SIMP, it is attributed to the release of the design freedom. In contrast to example 1, the microstructures of each layer are not restricted in an identical GM, thus the density optimization appears the polarization to 0 and 1.

Conclusions
This paper develops a multiscale isogeometric topology optimization method by the establishment of the structural similar microstructures. A shape deformation method is proposed to transform the PM for yielding a series of GMs where microstructural skeleton based on the level set framework is unified to retain more topology features and improve the connectability among the microstructures. Thus, the effective mechanical properties can be estimated by means of the numerical homogenization method. The global design problem is divided into two scale as the configuration design of PM in micro scale and the distributions optimization of GMs in macro scale. For the multiscale scheme, the material distribution-based method is applied to assemble the GMs while the isogeometric parameterized LSM is utilized for the structure optimization of PM. This method inherits the high accuracy and efficiency of IGA and the smooth boundaries and distinct interfaces of LSM. The introduction of structural skeleton contributes to maintaining the microstructural topology features especially for the low-density GMs. Thus, the more plentiful candidates with good connectivity can be assembled within the macrostructure for generating a preferable solution. The two numerical examples are presented to demonstrate the effectiveness in improving the performance and manufacturability.