Bézier extraction based isogeometric approach to multi‐objective topology optimization of periodic microstructures

This article presents the Bézier extraction based isogeometric approach to multi‐objective topology optimization of periodic microstructures. The approach utilizes the B‐splines based Bézier elements as the finite element representation at both macro and micro levels. The equivalent elastic properties of the representative volume element (RVE) are computed by the numerical homogenization method with the periodic boundary conditions on the control points (CPs) of the B‐spline mesh. The multi‐objective optimization problems including the material bulk modulus maximization, negative Poisson's ratio, and concurrent topology optimization of the composite structures are formulated by using the Bézier elements based isogeometric analysis. The distance‐based density distribution on the CPs of the B‐spline mesh is proposed as the initialization for designing the RVE, which is more robust than the uniform density distribution towards the optimal results. Several numerical examples are presented to illustrate the effectiveness of the proposed approach, and a variety of isotropic and anisotropic RVEs and composite structures are obtained. Meanwhile, various influences on the optimal design of the RVE and macrostructure are also discussed.

behavior of the RVE. 4,5 In recent years, the additive manufacturing techniques enable the possibility of successfully fabricating these microstructures that are designed by TO with the prescribed elasticity tensor. 6 Due to the special physical properties of the maximum stiffness and thermal conductivity, 7,8 negative Poisson's ratio (NPR), 9 and prescribed elasticity tensor, 10 these composite materials of man-made microstructures have attracted attention of many engineering applications for structural stiffness, thermal transfer and energy absorption among a broad range of fields, such as aerospace, transportation, civil and mechanical engineering. 11 An efficient method for designing the composites of periodic microstructures is the TO method. TO using the homogenization method was originally developed by Bendsøe and Kikuchi. 12 The homogenized elasticity tensor was computed in terms of a mutual energy concept after solving the characteristic deformation equation with the periodic boundary conditions (PBCs). The finite element (FE) formulation in terms of element mutual energy was used to design the material RVE with the prescribed constitutive properties by using the solid isotropic material with penalization (SIMP) model. The optimization problems of the constitutive tensor and NPR were formulated. A hierarchical computational procedure was proposed by Rodrigues et al., 13 which integrated the macrostructure with local microstructures using the SIMP model. A modified SIMP approach was employed to implement the concurrent topological design of the micro and macro structures. 14 The energy-based homogenization method was used to evaluate the elasticity matrix of the RVE. TO of the thermoelastic composites was formulated to maximize the stiffness and heat transfer, which was solved by the method of moving asymptotes (MMA). 15 In recent years, isogeometric analysis (IGA) as one of the most active research topics on the numerical approximation of partial differential equations was originally developed by Hughes and his co-workers. [16][17][18][19][20] It can be considered as a generalization of the classical FE method, and the NURBS basis functions were taken as the shape functions. They demonstrated that the B-splines or NURBS functions can be used to describe the geometry as well as construct the FE approximation. With IGA, the interoperability and integration between CAD and FE analysis are enhanced, since the NURBS representation is a part of numerous industry-wide standards in CAD.
IGA has attracted the attention of researchers in the structural optimization field, and shows the most potential application for the density-based optimization problems. In the framework of the IGA-based approach, the density field or design variables were represented by a NURBS scalar function. For the minimization of structural compliance, the material distribution function was approximated by the NURBS basis functions over the whole design domain. The design variables were defined on the CPs of the NURBS surface. 21 The weights were also considered as the design variables of the NURBS-based SIMP method. A thorough study aiming at investigating the influences of the B-spline and NURBS characteristic parameters on the optimized topology was systematically performed. [22][23][24][25][26][27][28] For 2D problems, the z-coordinate of each CP was replaced by the value of the optimized density. The postprocessing was provided in detail and the NURBS-based density surface can be imported in a CAD software for the material boundaries and checking analysis. 22 The similar postprocessing can be implemented for 3D problems. 23 Considering manufacturability, the minimum length scale of geometrical requirement was controlled in a NURBS-based algorithm for TO. 24,25 An enhanced density distribution function using the Shepard function was presented for the isogeometric TO. 29 The optimal results were featured with the smooth boundaries and distinct interfaces between the solids and voids. The multi-material isogeometric TO was also investigated for the optimization of multiple materials with the NURBS-based multi-material interpolation model. 30 An intrinsic filter for TO was analyzed with the IGA-based density distribution model. 31 The filter size was determined by the B-spline degrees and knot spans. According to the density distribution model, the IGA-based filter belonged to the sensitivity filter. A phase field model for the minimum compliance was developed with IGA. 32 IGA was combined with the level set function for TO and the NURBS basis functions were used to parameterize the level set function, which extended the applications of the IGA and level set method. [33][34][35] Based upon the above framework, the isogeometric TO is considered as an alternative approach to design the periodic microstructures with the desired properties. The NURBS hyper-surface based density field was combined with the strain energy-based homogenization method for TO of periodic microstructures by taking into account for the scale separation constraints and manufacturing requirements. 28 A rigorous discussion about the influence of the PBCs of the Dirichlet's type on the sensitivity analysis was carried out, and a rigorous proof of the effectiveness of the homogenization procedure based on the element strain energy was further provided. The IGA-based homogenization method was implemented with the smooth and continuous density distribution function. [36][37][38] Different deformation mechanisms were achieved with sufficiently smooth 2D or 3D material boundaries. The multi-material TO was also investigated for the design of the periodic microstructures. For the optimal design of the sandwich panels in real-world engineering problems, the NURBS basis functions were used to describe the RVE for the least-weight design subject to constraints of different properties. [39][40][41][42] The other applications of the IGA related to the displacement, 26 eigenfrequencies, 27 and structural stress 43 were also integrated into TO making use of the NURBS-based material interpolation model. Gao et al. provided a comprehensive discussion about the isogeometric TO in methods and applications. 44 The main directions for the development of the isogeometric TO were discussed in detail.
As discussed above, the IGA-based TO has the intrinsic filter which avoids the checkerboards and generates optimal results towards the black-and-white solutions. Thus, the present work proposes an IGA-based two-scale structural TO of composites. To facilitate the B-splines into the existing FE and TO procedures, the Bézier extraction based IGA 45,46 is incorporated into structural and sensitivity analyses. The Bézier elements can be applied to the standard FE programs through the Bézier extraction operator, and utilize the similar data processing structure. Thus, the introduction of the Bézier extraction allows the IGA to be performed with minor alterations to an existing FE program. 47 IGA based on the Bézier extraction preserves the features and advantages of the traditional IGA. In the Bézier extraction framework, the basis functions, which describe the geometry, are also used as the basis functions for FE analysis. Three multi-objective optimization problems, including the material bulk modulus maximization, NPR, and concurrent topology optimization (CTO) of composite structures, are formulated by using the isogeometric Bézier elements. For TO of seeking minimum compliance with a prescribed volume constraint, the initialization of the design variables usually consists of a uniformly distributed density field for the density-based optimization model. However, the uniformly distributed initialization for designing the RVE generates the uniformly distributed derivatives of the objective function, which fails to obtain the optimal solution by the gradient-based methods. Thus, the distance-based density distribution on the CPs of the B-spline mesh is proposed as the initial guess, which shows better performance towards the optimal design of the RVE. To compare with previous works, the distance-based density distribution is applied to the optimal designs of the material bulk modulus maximization, NPR, and composite structures, respectively. In summary, this article has the following main contributions: (1) The Bézier extraction based IGA is formulated to the numerical homogenization of the RVE. (2) The distance-based density distribution on the CPs as the initialization presents better performance towards the optimal design of the RVE. (3) The multi-objective CTO of composites with periodic microstructures is formulated on the basis of the Bézier elements. (4) The sensitivity analyses of the optimization problems are derived with respect to the design variables defined on the CPs. (5) The Bézier elements based IGA facilitates the incorporation of the existing TO programs.
The remainder of this article is organized as follows: Section 2 presents the Bézier extraction based IGA and the numerical homogenization with the PBCs on the CPs. In Section 3, the multi-objective optimization problems are formulated by using the B-spline shape functions based on the Bézier extraction operator. In Section 4, the sensitivity analyses and the numerical implementation with the distance-based density initialization are detailed. Section 5 presents several numerical examples to demonstrate the effectiveness of the proposed method. Finally, some concluding remarks are drawn in Section 6.

B-spline surfaces based on Bézier extraction
Given m × n CPs P i,j , and two scalar dimensionless parameters and in the parametric space, the NURBS surface is parametrically represented by 48 : where R p,q i,j ( , ) is the rational basis function, and it is defined by the tensor product functions with the following relationship: where i,j is the weight. It should be noted that the NURBS basis functions will become the B-spline basis functions when the weights are identical. The weights are set to be identical in this work. For the sake of completeness of the Bézier extraction, 45,46 the weights are preserved in the following derivation process. N i,p ( ) and N j,q ( ) represent the ith and jth B-spline basis functions with degrees p and q along the and directions, corresponding to the knot vectors = { 1 , 2 , … , m+p+1 } and H = { 1 2 , … , n+q+1 }, respectively. For a p-degree B-spline, the boundary knots can be repeated at most p + 1 times. The B-spline basis function N i,p ( ) is defined by using the Cox-de Boor recursive formula: wherep = 1, 2, … , p, and i is the ith component of the knot vector . Correspondingly, the B-spline basis function N j,q ( ) can be obtained on the knot vector H. For an in-depth presentation of the B-splines, the reader is referred to the excellent textbook. 48 The Bézier extraction operator provides an element structure that can be easily incorporated into the existing FE programs, which was proposed by Borden et al. 45 To derive the FE formulation, the matrix form of Equation (2) with the Bézier extraction operators is written as follows: where W is the diagonal matrix of weights, and C is the element Bézier extraction operator that is defined by the extraction operators C i and C j along the and directions. The element Bézier extraction operator can be computed as: The algorithm of the element Bézier extraction operators C i and C j is detailed in Reference 45.B( , ) in (5) is defined by: where B( ) and B( ) are the Bernstein polynomial vectors that can be defined recursively. In one dimension, the ith degree-p Bernstein polynomial function is defined as: where p is the degree of the Bernstein polynomials, and i = 0, 1, … , p. For a quadratic Bézier curve, the vector B( ) is defined by the Bernstein polynomials: Analogously, B( ) can also be defined. W( , ) in (5) has the following form: Therefore, the function R( , ) in vector form is the shape function following the traditional concept of the standard FE method.

Numerical homogenization using Bézier elements based IGA
The homogenization theory has been widely used to determine the mechanical properties of composite material, and TO using the optimality criteria (OC) method was further reviewed. [49][50][51] The region of the RVE is discretized and analyzed using the traditional FE method, which was usually referred as the numerical homogenization. 11,28,52 Using the homogenization theory, the macroscopic elasticity tensor E H ijkl of an RVE can be computed as [53][54][55][56][57] : where |Y | stands for the volume of the RVE, E pqrs is the elasticity tensor of the linear elasticity, and 0(ij) pq are the prescribed macroscopic unit strain fields. For 2D problems, the prescribed macroscopic unit strain fields take the following forms: The displacement fields are obtained by solving the following equilibrium [52][53][54][55][56][57] : where is a virtual displacement field. The strain fields (kl) rs in Equation (11) can be computed by the characteristic deformations obtained by solving Equation (13). Then, the homogenized elasticity tensor E H ijkl can be obtained by solving Equation (11) with respect to the principle coordinate axes of the RVE. Consider the macroscopic behavior of 2D composite material, the relationship between stress and strain is given in matrix form: where ij is the stress, ij denotes the strain, and the elastic tensor is written in matrix form. Thus, the homogenized elasticity tensor has the following components E H 1111 , E H 2222 , E H 1122 , and E H 1212 . To illustrate the nodes and degrees of freedom on the CPs of the B-spline mesh, a coarse example consisting of six Bézier elements is presented in Figure 1. The geometry can be exactly described by a quadratic B-spline surface. When the degree p = 2 of the Bernstein polynomial basis functions is selected, there are 2 + 1 Bernstein base polynomials generated by Equation (9). The knot vectors are separately set to = {0, 0, 0, 1, 2, 3, 3, 3} and H = {0, 0, 0, 1, 2, 2, 2} with 2 + 1 times repeats on the boundary knots. There are 8 and 7 knots along 2 directions. Therefore, the numbers of the basis functions are 8 -2 − 1 = 5 and 7 -2 − 1 = 4, respectively, corresponding to the CPs marked as the blue solid circles, as shown in Figure 1A. One can see that the displacement fields are solved on the CPs, and that the displacements of arbitrary points, as shown in Figure 1B, can be obtained by interpolation using the shape function Equation (5). The nodes distribution is non-uniform compared with the uniform distribution of the quadrilateral elements for TO of the RVE. 58,59 The PBCs are illustrated using Figure 1A,B. The nodes are defined on the CPs numbered as 1 , 2 , … , 12 , as presented in Figure 1A. Then, by mirroring the left nodes 1 , 5 , and 9 to the right surface, and then mirroring the bottom nodes 1 , 2 , 3 , 4 , and 1 to the top surface, the imposed boundary conditions defined on the CPs are generated using the same nodes numberings on two opposite surfaces, as shown in Figure 1B. The information of Figure 1B is used to construct the connection between these global and their local numberings on the elements (IEN array), as presented in Table 1. The IEN array with the PBCs is obviously different when compared to the IEN array of the traditional IGA presented in Table 2. The first columns of Tables 1 and 2 are the local numberings of elements, and other columns are the global numberings of elements. For the traditional IGA without the PBCs, the nodes are defined in order on the CPs including the upper and right surfaces from 1 to 20 , as shown in Figure 1C. Then, the IEN array is generated to show the connection between these global and their local numberings on the elements, as presented in Table 2. The PBCs are imposed in terms of nodal displacements which are the same on two opposite surfaces, as it can be inferred from the Figure 1B and Table 1. Figure 1D presents the Bézier elements computed by the extraction operator C. For example, the first Bézier element in Figure 1D is constructed by the CPs 1, 2, 3, 6, 7, 8, 11, 12, and 13 marked in blue circles, as shown in Figure 1A. Figure 1D  shows the mesh points of the first Bézier element 1, 2, 3, 8,9,10,15,16, and 17 marked in red circles. Another example is the third Bézier element with the PBCs constructed by the CPs 3,4,5,8,9,10,13,14,and 15, and its nodes numberings are 3 , 4 , 1 , 7 , 8 , 5 , 11 , 12 , and 9 , as presented in Figure 1B and the fourth column of Table 1.
Once the linear elastic problem (13), characterized by the PBCs of Figure 1B and Table 1, is solved the displacement fields can be obtained over the whole RVE. Then, the homogenized elasticity tensor E H ijkl can be calculated by using Equation (11). It is noteworthy that the assemblies of the global stiffness matrix and load matrix should use the information provided by the PBCs. The more detailed implementations can refer to the related works using the traditional FE analysis. 11,52,60,61 From a theoretical viewpoint, the work of Montemurro et al. 28 gave a deeper insight into this matter. The equilibrium Equation (13) for the REV was subject to the PBCs that were imposed in terms of displacements. Thus, the equilibrium problem (13) is of the Dirichlet's type. 28 The non-reduced (singular) stiffness matrix, displacement and force vectors of the FE model were separated according to the PBCs. Then, the equilibrium problem of the REV belonging to the Dirichlet's type was directly solved. Furthermore, the computational efficiency of the strain energy-based homogenization based on elements strain energy over that of the strain energy-based homogenization based on elements averaged stresses was rigorously proven.

MULTI-OBJECTIVE OPTIMIZATION PROBLEMS
The density-based approach to TO is used and each Bézier element is assigned a fictitious density x e with the following modified SIMP model 62 : where E 0 is the Young's modulus of solid material, E min is the Young's modulus of void region to prevent the stiffness matrix from becoming singular, and ℏ is a penalization factor to suppress the intermediate density and makes the optimal result towards the black-and-white solution. The fictitious densities at all CPs are the design variables. Applying the material interpolation scheme, the density x e of each Bézier element is determines by the shape function R e ( , ): where X e is a density vector constructed by the CPs related to the Bézier element e. For example, X e of element 1 is constructed by the density values at CPs 1, 2, 3, 6, 7, 8, 11, 12, and 13, as shown in Figure 1A. The Bézier element 1 consists of element CPs 1, 2, 3, 8, 9, 10, 15, 16, and 17, as shown in Figure 1D. The density of point 9 in the center of the Bézier element 1 can be computed by using Equation (16) with the parameters = 0.5 and = 0.5. The density in the center of each Bézier element is used to represent the element density during the optimization process. From Equations (11) and (14), the composite material with different microstructures and volume fractions holds different extremal bulk modulus, which affects the mechanical behavior of composites. Thus, a density-based TO, where the objective is to maximize the material bulk modulus, can be stated as follows: find ∶ x m , where 1 and 2 are the weighting factors from 0 to 1 satisfying 1 + 2 ≤ 1, V m , and V m0 are the material volume and design domain volume, m is the prescribed volume fraction, K m and are the global stiffness matrix and displacement fields derived from the left side of Equation (13), F m is the force matrix generated by the unit strain fields that take the following forms: 0(11) A second choice of the objective function is the material design of NPR, which is also a typical material design problem. Material, that will expand transversely when stretched longitudinally, is composed of special microstructures with NPR. TO was applied to the material design by Sigmund, 1 where the periodic microstructures were modeled using the homogenization method. A relaxed optimization form of NPR is adopted, 59,62 and the mathematical formulation is expressed as follows: where is a fixed parameter at each iteration, 1 and 2 are two constants, and is the iteration number. The other parameters are same as those of the optimization problem (17). If the NPR E H 1122 ∕E H 1111 or E H 1122 ∕E H 2222 is taken as the objective function, the material design with NPR requires imposing additional constraints for instance on isotropy or on bulk modulus. 1,59 Thus, the relaxed linear combination was proposed for designing the material with NPR. 59,63 For 2D problems, the relaxed form was originated by Xia and Breitkopf 59 for constructing the NPR materials. The objective function J 2 denotes a multi-objective design with the weighting factors 1 and ( ). As the iteration number increases, ( ) gradually approaches towards the prescribed value 1e−8 according to the setting. The multiple objectives of the optimal result have a tradeoff relation and are affected by the weighting factors. Therefore, the multi-objective function J 2 is formulated to obtain the maximization of the absolute value of E H 1122 and minimization of E H 1111 , which has an equivalent performance when using E H 1122 ∕E H 1111 as the objective function for the material design. The IGA based TO has the intrinsic filter to optimal result with structural compliance as the optimization objective. Therefore, both macrostructure and its periodic microstructures are concurrently optimized to verify the effectiveness of the Bézier elements based optimization at two scales. In the CTO model, the macro and micro densities are taken as the design variables for the macrostructure and RVE independently. The above concurrent multi-objective TO problem can be formulated as: , where the weighting factors , 1 , and 2 range from 0 to 1 with 1 + 2 ≤ 1, f 0 , and e 0 are the given values for the normalizations of structural compliance and extremal bulk modulus, f 0 can be estimated using structural compliance with the initial density at the first optimization step, e 0 can be computed by e 0 = 1 E 0 1111 , and E 0 1111 , E 0 2222 , and E 0 1122 are the material properties from the elasticity tensor, V and V 0 are the material volume and design domain volume of macrostructure, is the prescribed volume fraction for macrostructural TO, K is the global stiffness matrix, U and F are the global displacement and force vectors of macrostructure, and X is the vector of design variables for macrostructure. The other parameters are described as those of the problem (17). In the framework of the CTO, the overall mass of macrostructure is a non-linear constraint coupling the design variables defined at both lower and upper F I G U R E 2 Multi-objective CTO of macrostructure with its periodic microstructures scales, 28 and this constitutes a further complexity in the formulation of the problem (19). Montemurro et al. also derived the derivatives of the macrostructural mass with respect to the design variables. 28 The CTO problem (19) is an unrealistic and simplistic version when compared to real-world engineering problems as those presented in works of Montemurro and his co-works. [39][40][41][42] A two-level optimization strategy was proposed for designing the lattice structures made of periodic microstructures, which did not make any simplifying assumption to obtain a true global optimum configuration of the practical system. The least-weight of the periodic microstructure was the optimization objective subject to a given set of requirements including the first buckling load, the positive-definiteness of the stiffness tensor of the core, the ratio between skins and core thickness, and the admissible moduli for the laminated skins. Macrostructure is composed of the uniform periodic arrangement of RVEs, as shown in Figure 2. Macrostructural TO uses the homogenized elasticity tensor E H ijkl that is computed in terms of the material distribution of the RVE. For CTO, the RVE and macrostructure are simultaneously discretized with the B-spline meshes. The CPs of the B-spline mesh are non-uniform, as shown in Figure 1A, but the optimal results can be presented by the Bézier elements, that are computed in the physical space using the shape function R( , ), as shown in Figures 1D and 2.

Bézier elements based finite element formulation
The FE analysis should be conducted for structural analysis and TO at the macro and micro scales. The FE formulation of Equation (13) can be written as: where K m denotes the global stiffness matrix of the RVE, is the global displacement vector, and F m indicates the load matrix. The global stiffness matrix is assembled by the element stiffness matrix: where the summation denotes the assembly of N m FEs, and the element stiffness matrix k e is computed by: where B e ( , ) is the element strain-displacement matrix, and C e is the elasticity matrix. Using the shape function R e ( , ) with degrees 2, the element strain-displacement matrix B e is then written as: It must be noted that the assemblies of the global stiffness matrix and load matrix for Equation (20) should use the IEN array with the PBCs, as presented in Figure 1B and Table 1. Then, the other programs of the Bézier elements based IGA can fully use the existing traditional FE codes. 11,52,60,61 Therefore, the Bézier extraction based isogeometric approach facilitates the incorporation of the smooth B-splines into the existing FE codes. As shown in Figure 3, the Bézier element 2 as an example is extracted from Figure 1. The Bézier element consists of nine CPs that are computed by using the Bézier extraction operator and the CPs of the B-spline mesh, as shown in Figure 3. Therefore, when the polynomial orders of element are set to 2, the number of CPs of each Bézier element are nine.
In Equation (23), the derivatives R x and R y of the shape function R e ( , ) with respect to x and y in the physical space are computed by: where the transformation is defined as the Jacobian matrix: Transformation of a Bézier element between the physical space Ω e and the parametric domainΩ e It is computed by the transformation of a Bézier element from the parametric spacẽe to the physical domain e with the following geometrical mapping: where ⌢ x and ⌢ y are the Bézier element CPs, as shown in Figure 3. Thus, the determinant |J| of the Jacobian matrix in Equation (22) can be computed by Equation (25). For the plane stress problem, the elasticity modulus matrix C e in Equation (22) can be expressed as: where E(x e ) is the element Young's modulus, and is the Poisson's ratio. Using the PBCs, the global load matrix is obtained in the following form: Substituting the unit strain fields to the element load matrix yields: where the prescribed unit strain field 0 is the transpose of 0(11) Applying Equation (22), the FE formulation of Equation (11) takes the following form: where 0(ij) e are the displacement fields corresponding to three unit strains, the characteristic deformations kl are the global solutions by solving Equation (20) with the PBCs.
For the macrostructural analysis in Figure 2, the FE formulation is written as: where K and U are the global stiffness matrix and displacement field of macrostructure composed of periodic microstructures, and F is the force vector of macrostructure. The assembly of the global stiffness matrix of macrostructure can be completed by using the IEN array without the PBCs, as presented in Figure 1C and Table 2. In Equation (22), the elasticity matrix C e is replaced by the homogenized elasticity matrix for computing the element stiffness matrix, which is used to assemble the global stiffness matrix K. Thus, the element stiffness matrix for macrostructure is given as: where the element strain-displacement matrix B e ( , ) and the Jacobian matrix J are computed according to the CPs of the B-spline mesh for macrostructure.

The B-splines based sensitivity analysis
To implement TO of the RVE and CTO of the composites, the sensitivity results are required. The sensitivity of E H ijkl with respect to the Bézier element density x e is calculated from Equation (30) as follows 28,[53][54][55][56][57] : As indicated in work of Montemurro et al., 28 under the PBCs of the Dirichlet's type, the gradient of the elasticity tensor components at the macroscopic scale can be expressed in terms of the element strain energy. The detailed derivation process was provided by separating the global displacement vector into two parts: the known displacement values on the periodic boundaries and the unknowns corresponding to the interior nodes. Therefore, under the hypothesis that the external nodal forces (no external forces) and the known displacements on periodic boundaries do not depend on the design variables, the gradients of the elasticity tensor components at the macroscopic scale can be obtained in term of the element strain energy.
Using Equations (15) and (27), the derivative of Equation (22) with respect to the design variable x c defined on the CP reads: The sensitivity of E H ijkl with respect to the design variable x c can be calculated, using Equations (33) and (34), as, where I e represents the Bézier elements related to the CP with the design variable x c . Using the interpolation model (16), the gradient of x e with respect to x c is given as: where R c ( e , e ), which belongs to the shape function R e ( e , e ), is the term corresponding to the CP with the design variable x c , x c ∈ x m . Now, substituting Equations (16) and (36) into (35), the sensitivity of E H ijkl can be expressed as: Thus, the sensitivity of the objective function J 1 (E H ijkl (x m )) is found as: Similarly, the sensitivity of the objective function J 2 (E H ijkl (x m )) is: For CTO, the element stiffness matrix (32) of macrostructure is defined by: where X mac can be computed using Equation (16) at the macro-scale level. The sensitivity result of J 3 (x) with respect to X c ∈ X is given as: where the CP vector X e belongs to the B-spline mesh of macrostructure, and R e ( e , e ) is computed according to the CPs of the macrostructural B-spline mesh. At the micro-scale level, the sensitivity result of J 3 (x) with respect to x c ∈ x m is given by: where N is the macrostructural FE number. Using Equations (32) and (40), the derivative of K e (X mac ) can be obtained: Therefore, substituting the results (33) and (43) into (42) yields the derivative of J 3 (x) with respect to x c . Using Equation (16) for the concurrent design, the derivatives of the material volume are given as: where R c ( e , e ), which belongs to the shape function R e ( e , e ), is the component related to the CP with the design variable x c of the RVE, x c ∈ x m . Similarly, the sensitivity of the material volume V(X) for macrostructural TO can also be obtained using the CPs of the macrostructural B-spline mesh.

Numerical implementation
The MMA algorithm is used to solve the multi-objective TO problems. For TO of the RVE, the initial guess has effect on the optimal results. 1,59,64 In other words, different initializations may generate different solutions. The design variables are initialized as a uniform density distribution for macrostructural TO when the density-based model is adopted. However, the uniform initialization cannot be applied to the design of the RVE when the PBCs are imposed to the optimization problem, which yields a uniformly distributed sensitivity result and fails to obtain the solution. Thus, the density distributions with the circular regions of different radii were proposed, as shown in Figure 4B for material design, which generates satisfactory results. 59 The distance-based density distribution as initialization is proposed to design the RVE, which obtains the optimum components of the elasticity tensor for the following macrostructural optimization. The distance-based density distribution is defined as: where is a constant that is usually defined as the prescribed volume fraction, which means that the maximum value in Figure 4A is the prescribed volume fraction, is a minimum relative density, and (P i,j ) is defined according to a distance function: where H = 0 or H = 1, P i,j are the m × n CPs of the B-spline mesh, and P c is the center point of the B-spline mesh. Figure 4A shows the results when = 0.5, H = 1, and H = 0, respectively. The values of the distance function Equation (46) range from 0 to 1. The distance-based density distribution of Equation (45) is constructed according to CPs, which is applied to the initialization of the design of RVE. In the macro scale design, the design variables are initialized as a uniform density distribution for macrostructural TO when the density-based model is adopted. The optimization problems (17) and (18) can be taken as the special cases of the CTO problem (19). Thus, the CTO procedure is outlined as follows: Step 1 Initializations for the CTO of the RVE and macrostructure at two scales.
(1) Initialize the design domains of the RVE and macrostructure, give the mechanical properties: Young's moduli E 0 , E min and Poisson's ratio . (2) Set up the CPs of the B-spline mesh, knot vectors, and degrees to describe the design domains at two scales.
Construct the IEN arrays and compute the Bézier extraction operators. (3) Initialize the parameters of the MMA algorithm. The distance-based density distribution defined on the CPs as the topological guess is adopted for designing the RVE. The uniform density distribution is used for the macrostructural TO.
Step 2 Implement the FE analysis of the RVE using Equation (20). The global stiffness and load matrices are computed according to Equations (21) and (28) with the imposed PBCs, as shown in Figure 1B. Step 3 Using the solution of Step 2, compute the homogenized elasticity tensor E H ijkl by Equation (30). Substitute the resulting E H ijkl into Equation (32) and implement macrostructural FE analysis with Equation (31).
Step 4 Calculate the objective function J 3 (x), volume constraints V(X) and V m (x m ) of the problem (19), and the corresponding sensitivity results using Equations (41)- (44). Implement the MMA algorithm to update the design variablesx in the CTO design.
Step 5 Repeat Steps 2-4 until the termination condition is satisfied or the maximum iteration number is reached.
Once the IEN arrays and Bézier extraction operators have been initialized, the results do not require recomputation during the optimization process. In the MATLAB code, the matrix edofMat and vector edofVec 58 schemes can be explored to efficiently assemble the global FE matrices K m and K for the concurrent designs. According to the above optimization procedure, the Bézier elements and shape function (5) can be processed in the similar way as those in the existing TO program. The MMA algorithm is used to solve the TO problem, the parameters tuning the performance of the MMA algorithm are listed in Table 3. The optimization iteration is terminated when the convergence criterion or the maximum iteration number is attained. The convergence criterion is that the maximum absolute value (MAV) of densities changes between two consecutive iterations is below a specified error limit. The specified error value is set to 0.01 for 2D problems. To reduce the total computational time, the convergence criterion is set to 0.02 for 3D TO problems.

NUMERICAL RESULTS
In this section, several examples for designing the RVE and macrostructure are presented. The square region is defined as the design domain of the RVE. The distance-based density distribution and the uniform density distribution with the specified circular regions, as shown in Figure 4, are taken as the initial guesses for designing the RVE. Different optimal results are compared to show the effectiveness of the Bézier elements based IGA and the distance-based density initialization for designing the RVE. The TO codes are written in MATLAB R2020b Version. All computations are performed on a PC with Intel(R) Core(TM) i7-8850H CPU @ 2.60GHz 2.59GHz processor and 64.0GB RAM. The iteration number to achieve convergence, the values of the constraint functions, and the CPU time are all provided for each optimal result.

Periodic microstructures with bulk modulus maximization
The design domain of the RVE is a square area with size 1 × 1 at the macro-scale level for Figures 5 and 6. This case corresponds to the optimization problem (17). The square base cell is discretized into 80 × 80 Bézier elements. The B-spline degrees along the and directions are set to p = 2 and q = 2. The knot vectors and H with repeated 3 times on the boundaries are defined as = {0, 0, 0, 1, … , 79, 80, 80, 80} 1×85 , and H = {0, 0, 0, 1, … , 79, 80, 80, 80} 1×85 , respectively. Therefore, the number of the basis functions is 85-2-1 = 82 along two directions. According to the CPs of the B-spline mesh as shown in Figure 1A, the CP vectors denoted by P are set to P x = {0, 0.0063, 0.0188, … , 0.9813, 0.9938, 1} 6742×1 and P y = {0, 0, 0, … , 1, 1, 1} 6742×1 . The weights in the shape functions are set to i,j = 1. The prescribed volume fraction of the RVE is set to m = 0.5. The solid material is assumed to be isotropic with the Young's modulus E 0 = 1, E min = 1e − 9 and Poisson's ratio = 0.3. The penalization factor ℏ is set to 5. Figures 5 and 6 show the optimal results of the RVEs and the corresponding elasticity matrices starting from different initializations. With the distance-based density distribution, the parameter in Equation (45) is equal to the prescribed volume fraction 0.5. H = 1 and H = 0 in Equation (45) are used for the first rows of Figures 5 and 6, respectively. For the uniform density distribution with the circular regions, the diameters are set to 1∕3 and 2∕3 of the height of the design domain. When 1 = 0.25 and 2 = 0.25 are adopted for the problem (17), the final solution presented in Figure 5A with the initialization of the distance-based density distribution shows better performances than those in Figure 5B,C according to the material distributions and optimal objectives. When the parameters 1 = 0.1 and 2 = 0.1 are set, the results are given in Figure 6. The optimal result in Figure 6A presents better performance than those of Figure 6B,C in terms of the objective values.
For Figure 7, the design domain is a rectangular area with size 1.2 × 0.6 and discretized into 120 × 60 Bézier elements. The distance-based initialization is more robust than the other initial guesses. Furthermore, the classical filter for the density-based TO 58 does not apply to the optimization process to prevent the checker-boards and ensure the existence of solution. The optimized results are almost black-and-white, as shown in the second columns of Figures 5-7. The structural topology of the RVE with the distance-based initialization is simpler than the optimal results starting from the circular initializations. Moreover, the objective values starting from the distance-based initialization are slightly smaller than the values obtained from the circular initializations. The optimal results are presented with the element densities, as shown in Figures 5-7. For the exact description of geometric boundaries, the smooth nodal densities using the Shepard function were incorporated into the NURBS-based TO framework. 29,30,36 The enhanced density distribution function generated the desired smoothness and continuity for the design of continuum structures. According to the numerical results, the final topological distributions of material were featured with smooth boundaries and distinct interfaces between the solids and voids. The promising results showed the potential performance for the TO of continuum structures. The future work of this article would explore the Shepard function for the smooth nodal densities due to its excellent characteristic of exact geometry. Costa et al. presented the postprocessing of the optimal result in a CAD software. 22 The fictitious densities were taken as the z-coordinates of the CPs, and the description of the NURBS-based density surface was obtained. The NURBS surface can be imported in a CAD software for extracting the 2D structural boundaries. Using the B-splines based Bézier elements, the density surface can be obtained using Equation (1), where the CPs P i,j are replaced by the fictitious densities on the CPs. For the representation between density and geometry, a post-processing scheme to define the topology is the same as the works of Costa et al. 22 and Gao et al., 29 as: boundary, where x b are the density values defined on the mesh points of Bézier elements, and x c is a threshold value to cut the density surface. The material topology is obtained by cutting the density surface to obtain the isocontour of the density surface with a threshold value x c = 0.5. As shown in Figures 8-10, the geometrical boundaries tend to become smooth as a convergent solution is found. In particular, the convergent solutions obtained by means of the NURBS surfaces exhibit smoother boundaries than those obtained through the B-spline ones, with better performances for structural compliance. 22,23 The weights are considered as the design variables in their works. With the isocontour of the density surface based representation, the volume fractions of the optimal results of Figures 8-10, which are close to the prescribed volume fractions, are 0.5, 0.4981, and 0.5039, respectively.
The first-order basis functions are used for the density-based TO model, and the optimal results of the periodic microstructures with bulk modulus maximization are given in Figure 11. The optimal results of Figure 11 are compared with Figures 5A, 6A, and 7A. It is observed that the optimal shape of the RVE is a little different when comparing results of Figures 5-7 and 11. Furthermore, when the first-order basis functions are adopted, the intermediate densities in the central region appear as the convergence solution is achieved, as shown in Figure 11C. The reason underlying the intermediate densities is the filter zone that is affected by the degrees of the B-spline basis functions. 22,31 The optimization objective increases with the increasement of the order of the B-spline basis functions. The computational performances of optimal results are given in Table 4. The mean computational time decreases as the first-order basis functions are adopted. Figure 5A. Isocontour of the density surface based representation: the first row is the density surfaces in the physical space, and the second row is the corresponding isocontours and material topologies. The final volume fraction is 0.5 with the isocontour of the density surface based representation F I G U R E 9 Iteration process of Figure 6A. Isocontour of the density surface based representation: the first row is the density surfaces, and the second row is the corresponding isocontours and material topologies. The final volume fraction is 0.4981 with the isocontour of the density surface based representation F I G U R E 10 Iteration process of Figure 7A. Isocontour of the density surface based representation: the first row is the density surfaces, and the second row is the corresponding isocontours and material topologies. The final volume fraction is 0.5039 with the isocontour of the density surface based representation Figures 5-7 Table 4 shows the computational performances for designing the periodic microstructures with the bulk modulus maximization. The mean CPU time (s) for each iteration is very close. Thus, the total CPU time increases with the number of iterations, as presented in Table 4. Figures 5A, 6A, and 7A using the distance-based density distributions as initializations achieve the optimal results with less iterations and better performances in terms of the objective values. When the initial density distribution is approximate to the final solution, the iteration number can be reduced, as presented in Figure 6C of Table 4.

TA B L E 4 Computational performances of optimal results in
The optimization iteration is terminated when the convergence criterion or the 1000 design iterations are attained. The convergence criterion is that the MAV of densities changes is below 0.01. The objective value and KKT norm can also be considered as the stopping criteria according to the MMA algorithm. Table 5 gives the final parameters when the MAV of densities changes and the maximum number of iterations are taken as the termination criteria, respectively. The

F I G U R E 11
Optimal results of the RVEs using the first-order basis functions, and distance-based density distribution as initialization: objective and elasticity matrix (left), RVE (middle) and 3 × 3 microstructures (right) objective change and KKT norm have different magnitudes and they are difficult to be selected as the termination criteria for TO, as presented in Table 5. Therefore, a prescribed value about the objective change or KKT norm cannot be easily determined for the termination condition since the ranges of their values cannot be predicted and change with the TO problems, while the densities variables hold the fixed range. For case 1 in Table 5, the MAV of densities changes, objective change and KKT norm are normalized to values between 0 and 1 after 1000 iterations, as shown in Figure 12. The curve of the MAV of densities changes drops more slowly than the other two curves. Therefore, the MAV of densities changes as the termination condition is reasonably, which can achieve a steady state after those of the objective change and KKT norm. The OC based TO method also takes the MAV of densities changes as the termination condition, 2,58,59 which further illustrates the selection of the termination condition is effective. For TO of the bulk modulus maximization, the optimal results have the same shape and topology compared with those of two types of termination conditions. According to the MAV of densities changes in Table 5, the objective value has slight difference due to the different iteration numbers for two types of termination conditions. For the following material designs of NRP and CTO, the similar conclusions regarding the selection of the termination condition can be found.
For 3D problems, the design domain is a cube with size 1 × 1 × 1 at the macro-scale level, as shown in Figure 13. The 3D RVE is discretized into 20 × 20 × 20 Bézier elements. The B-spline degrees are set to 2 along three directions. The knot vectors with repeated 3 times on the boundaries are defined as {0, 0, 0, 1, … , 19, 20, 20, 20} 1×25 . The CP vectors of the B-spline hyper-surfaces denoted by P are set to P x = {0, 0.025, 0.075, … , 0.925, 0.975, 1}, P y = {0, 0, 0, … , 1, 1, 1}, and P z = {0, 0, 0, … , 1, 1, 1} of dimension 10,648 × 1. The prescribed volume fraction is set to m = 0.5. The solid material is assumed to be isotropic with the Young's modulus E 0 = 1, E min = 1e − 9, and Poisson's ratio = 0.3. The penalization factor ℏ is set to 5. The optimal results are presented in Figures 14 and 15 using different initializations. The distance-based densities values are defined in descending order from the center to the outside of the cube ranged from = m to 1e-3, as shown in Figure 13A, while the densities values are arranged in ascending order, as shown in Figure 13B. The objective function for 3D problems is defined as J 1 (x m ) = − (

F I G U R E 12
Iteration processes of the MAV of densities changes, objective change and KKT norm after 1000 iterations for case 1 in Table 5 F I G U R E 13 Two initial designs for RVE Homogenized elasticity matrix: 3D cases, please refer to References 57, 61, 65, and for more details. The PBCs and numerical implements were provided using the traditional FE method for 3D problems. 11 The PBCs based on the CPs of the B-spline surface can be deduced from the results of Section 2.2 and the work of Dong et al. 11 The convergence criterion for 3D problems is that the MAV of densities changes is below 0.02. With different initializations, different topological distributions of material are obtained while the objective values are close with −0.2199 and − 0.2211, as presented in Table 6. The results are provided as the optional cases for designer when the manufacturing requirements are considered. Furthermore, the mean CPU times for two cases are 8.60 and 8.78, respectively, which largely exceeds those of 2D problems, as presented in Table 4.

Periodic microstructures with negative Poisson's ratio
To show the effectiveness of the Bézier extraction based IGA and the distance-based density distribution for TO of the RVE, the proposed method is further used to design the RVE of NPR. The design domain, B-spline mesh and physical parameters for TO are the same as those for Figure 5, unless otherwise stated. Xia and Breitkopf 59 pointed out that the optimized solutions were more sensitive to the initial guess and other parameters as compared to the material design of the bulk modulus maximization. The Bézier extraction based IGA for TO, that is applied to the material design of NPR, also fails in obtaining the optimal results without filter. Thus, the density filter is used with the filter radius r min = 3. The penalization factor ℏ is set to 2 for Equation (15). The optimized solutions are presented in Table 7 with different initial guesses. The objectives and NPRs are given in the second column, which uses the distance-based density distribution with different in Equation (45) as initializations. It can be seen that the satisfactory result is obtained with = 0.5 ( m = 0.5) in the fourth row. Therefore, it is effective and robust for TO of the RVE when the parameter is equal to the prescribed volume fraction m according to the previous examples and Table 7. The last two rows with the initializations of circular regions in Table 7  designs, which illustrates the dependence of the optimal results on the initial guesses. Meanwhile, the sequences of iterations for TO of NPR including the objective and volume fraction are shown in Figure 16, where the intermediate designs are also presented corresponding to the case of = 0.5 in Table 7. The number of the intermediate densities is reduced. The elements with the intermediate densities distribute along the interfaces between the solids and voids, as shown in Table 7 with = 0.5. Therefore, the intrinsic filter provided by the B-splines has effect on the optimal result. Figure 17 is the iteration process of the density surfaces for case 3 of Table 7 and their corresponding material topologies in the physical space. The similar conclusion can be drawn as those of Figures 8-10. With the isocontour of the density surface based representation, the final volume fraction 0.5016 is close to the prescribed volume fraction 0.5.
The optimal result with the first-order basis functions is given in Figure 18. Comparing with the first-order basis functions for the periodic microstructures with bulk modulus maximization, the similar conclusion can be drawn. A higher order value of the B-spline basis functions can achieve better convergence solution according to the intermediate densities and optimization objective. The computational performances are presented in Table 8 including the iteration number and CPU time. When the iteration step reaches 968, the termination criteria is satisfied, as presented in Table 8.
The material design of NPR is also implemented using the traditional four-node quadrilateral elements 58,59 to compare the computational cost and efficiency with those of the Bézier elements based IGA method. The density filter is used to ensure the existence of solution. The MMA algorithm is used to solve the optimization problem. The other parameters used for Figure 19 are the same as those for Table 7. The MAV of densities changes is 0.3242 after 1000 iterations, as shown in Figure 20B. There are oscillations along the curves of the objective function, volume fraction and densities changes, F I G U R E 16 Iterative processes of designing NPR material with = 0.5 in Table 7 F I G U R E 17 Iteration process of case 3 in Table 7. Isocontour of the density surface based representation: the first row is the density surfaces, and the second row is the corresponding isocontours and material topologies. The final volume fraction is 0.5016 with the isocontour of the density surface based representation F I G U R E 18 Optimal result of designing material with NPR using the first-order basis functions and distance-based density distribution as initialization ( = 0.5) Table 7   as shown in Figure 20. For the IGA-based TO of NPR, the iteration number is 414 when the termination condition is satisfied, as shown in Figure 16. The total CPU time is 138.52 corresponding to the case 3 presented in Table 8. For the quadrilateral element based optimization, the total CPU time is 3.27e2 given in the last column of Table 8. However, the mean CPU time of each iteration is very close for each case presented in Table 8. The iteration number, total CPU time, mean CPU time and volume fraction are presented in Table 8, when the convergence criterion or the maximum iteration number is satisfied. The mean CPU time for each iteration is very close in terms of the results presented in Table 8. For the material design of NPR, the iteration number is affected by the initial guess of the density distribution. By comparing with Tables 4 and 8, the mean CPU times are also close for these two types of optimization problems.

TA B L E 8 Computational performances of optimal results in
For 3D material design of NPR, the distance-based density distribution and multiple holes initializations are both considered, as shown in Figure 21. The other parameters are the same as those used in Section 5.1 for 3D problems. The optimal result cannot be easily achieved for the optimization model (18) when the numerical examples are tested. Thus, the following objective function 1,65,66 is used to replace that in the problem (18): F I G U R E 20 Iterative processes of TO for NPR using the quadrilateral elements and MMA algorithm

F I G U R E 21 Two initial designs for NPR
where is the weighting factor, and P J 2 is the penalization term that is added to the objective function for making the components of the elasticity tensor equal, and defined as: In this example, is set to 0.01 and 0.02 for two initial cases. The prescribed components of the elasticity tensor are set to E * 1111 = E * 2222 = E * 3333 = 0.1 and E * 1122 = E * 1133 = E * 2233 = −0.01. According to the optimal results presented in Table 9, the components of the elasticity tensor approach towards the prescribed constants. The final topologies are largely different when starting from different initializations, as presented in Figures 22 and 23, which demonstrates the dependence of the optimal result on the initial guess. However, the result obtained from the distance-based density distribution holds better performance in terms of the objective value.

Concurrent topology optimization at two scales
The first CTO of macrostructure with periodic microstructures is a half MBB beam, as shown in Figure Homogenized elasticity matrix:  According to the B-spline mesh, as shown in Figure 1C without the PBCs, the CP vectors P are set to P x = {0, 0.5, 1.5, … , 98.5, 99.5, 100} 5304×1 and P y = {0, 0, 0, … , 50, 50, 50} 5304×1 . The settings for the B-spline mesh of the RVE are like those of example of Figure 5. As shown in Tables 10 and 11, the optimal results are given with different volume fractions. The weighting factors , 1 , and 2 are set to 0.5, 0.25, and 0.25. The penalization factor ℏ in Equation (15) is set to 5 and 3 for designing the RVE and macrostructure. The other physical parameters are the same as those of above example. Table 10 shows the final solutions with the prescribed volume fractions = 0.4 and m = 0.2. The optimized solution cannot be reached with the initial guess of small circular region, so it is not listed in Table 10. Table 11 lists  Figure 25 shows the sequences of iterations including the macro and micro volume fractions for the half MBB beam with = 0.5, m = 0.5, and = 0.5 in Table 11, where the intermediate designs of the RVE and macrostructure are also presented at two scales.
According to the optimal results, as presented Tables 10 and 11, the thickness of the smallest topological branches of the macrostructure is approximate to one or two RVEs. However, when the PBCs are imposed for the homogenized elasticity tensor, the infinite number of the RVEs for this application of the PBCs should be satisfied. 28,60 Thus, the macrostructural result does not satisfy the hypothesis. Therefore, the minimum thickness scale could be considered valuable from a practical viewpoint. Furthermore, the manufacturing and geometrical requirements should be further addressed for real-world engineering applications. 24,25,[39][40][41][42] In this work, the hypothesis of an infinite number of RVEs, that allows for the application of the PBCs on a single RVE, is neglected for the thin topological branches that could appear during the optimization process of CTO.

F I G U R E 25
Iterative processes of CTO of half MBB beam with = 0.5 in Table 11 F I G U R E 26 A cantilever beam for CTO A cantilever beam for CTO at two scales is shown in Figure 26, where a force F = −1 is applied at the right side of the cantilever beam. The prescribed volume fractions are = 0.5 and m = 0.5 for the RVE and macrostructure. The other parameters are the same as those of the above half MBB beam. Table 12 presents the resulting RVEs and macrostructures starting from the initializations of the distance and circular region based density distributions. The final solutions with different are compared in Table 12 when the distance-based density distribution is used as the initial guess. It can be seen that the objective, shape, and topology with = m are similar to the optimal results when the big circular region as the initialization is adopted in the last row of Table 12. Figure 27 shows the CTO process of case = 0.5 with the distance-based density distribution.
For the CTO, the computational performances related to the iteration number, total CPU time, mean CPU time, and volume fractions of the RVE and macrostructure are presented in Table 13. m_vol fraction and M_vol fraction represent the RVE and macrostructural volume fractions, respectively. Compared with the computational time of the bulk modulus maximization in Table 4, the mean CPU time of CTO increases with the combination of macrostructural TO. The number of iterations also becomes larger than that of the material design of the bulk modulus maximization due to more design variables.
In macro scale design, the design variables are initialized as a uniform density distribution for macrostructural TO when the density-based model is adopted. The isocontour of the density surface based representation for CTO of case 3 in Table 12 can also be obtained with Equation (1), as shown in Figure 28. The final volume fractions of the RVE and macrostructure, which are close to the prescribed volume fractions, are 0.5 and 0.5024, respectively.
For optimization problem with several objective functions, the most common way to obtain an optimal solution, though not the most rigorous, is to form a scalar objective function by using the weighting factors. Then, the optimization problem can be solved by the single-objective optimizer with specified weighting factor for each objective. The result of the multi-objective optimization depends on the choice of the weighting factors , 1 , and 2 . The optimal results are presented in Table 14 with different combinations of the weighting factors. The RVE and macrostructure obtained from parameters = 0.01, 1 = 0.25, and 2 = 0.25 present obviously different shape and topology when compared to those TA B L E 12 The optimal results of CTO of cantilever beam with = 0.5 and m = 0.5

F I G U R E 27
Iterative processes of CTO of cantilever beam with = 0.5 in Table 12 TA B L E 13 Computational performances of optimal results in Tables 10-12

F I G U R E 28
Iteration process of the CTO of case 3 in Table 12. Isocontour of the density surface based representation: the first and third rows are the RVEs, and the second and fourth rows are the corresponding macrostructural topologies. The final volume fractions of the RVE and macrostructure are 0.5 and 0.5024, respectively with = 1.0, 1 = 0.25, and 2 = 0.25. With the decrease of , the objective function J 3 becomes smaller, which means the influence related to the RVE on the objective function becomes larger. The weighting factors 1 and 2 affect the material distribution of the RVE, as shown in the last two rows of Table 14. When the weighting factor 2 is set to 0.9, there is a change corresponding to E H 2222 . It is expected that with different weighting factors, more optimal solutions are provided. The designer can observe the changes of different objectives related to the RVEs and macrostructures, and choose a final design more rationally according to their practical requirement. By varying the weighting factors , 1 , and 2 , different optimal solutions are obtained. It should be remarked, however, that in general not every optimal result may be obtained TA B L E 14 The optimal results of CTO of cantilever beam with different weighting factors , 1 Table 15 corresponding to the results in Table 14. The weighting factors have effect on the final results, so the iteration number and total CPU time are also different. Figure 29 is the optimal results using the first-order basis functions. The similar conclusion can be drawn as previous designs of periodic microstructures with bulk modulus maximization and NPR when using the first-order basis functions. The intermediate densities appear in the optimal results of CTO of the half MBB beam. Table 16 gives the computational performances in comparison with those of degrees p = 2 and q = 2 in Table 13. Figure 30 is the optimized RVE and macrostructure using the quadrilateral elements and MMA algorithm. The parameters are the same as those in Table 12. There are also oscillations along the curves of the objective function, volume fractions and densities changes, as shown in Figure 30D,E. When the iteration step reaches 1000, the MAV of densities changes is 0.3203. The RVE and macrostructural volume fractions cannot be satisfied due to the volume oscillations. For the Bézier elements based method with = 0.5, the convergence iteration is 455, as presented in Table 13. The total CPU time is 1.99e3, and the mean CPU time is 4.37 presented as case 3 in Table 13. Though the Bézier elements based method takes a little more time than the quadrilateral element based TO by comparing the mean CPU times in Tables 13 and 17, it reaches the termination condition with less iterations than the quadrilateral elements based TO with the same parameter settings.

F I G U R E 29
Optimal results of CTO using the first-order basis functions, and distance-based density distribution as initialization for RVE. Three columns are the initial design, optimal RVE, and optimal macrostructure, respectively The CTO developed by Gao et al. 14 is also implemented. The OC method is used to solve the optimization problem. The optimization objective is the structural compliance. The input parameters are given as: Macro_struct = [100, 50, 100, 50, 0.5], Micro_struct = [1, 1, 60, 60, 0.5], penal = 3 and rmin = 2. The optimal result is shown in Figure 31. The performances of the concurrent designs of the cantilever beam are presented in Table 17 using the MMA and OC algorithms, respectively. The computational efficiency with the quadrilateral elements and OC algorithm is higher than that of the Bézier elements and MMA algorithm in terms of the mean CPU time by comparing Tables 13 and 17, since IGA takes more time than the traditional FE analysis. The convergence iteration is 504 when the termination condition is reached by the OC method. The convergence iterations are close by comparing the quadrilateral elements and OC algorithm with the Bézier elements and MMA algorithm.

F I G U R E 31
Iterative processes of CTO of cantilever beam using the quadrilateral elements and OC algorithm F I G U R E 32 Two initial designs for 3D RVE initialization, as presented in Figure 32A. The iteration process takes more time than that of Figure 34. Thus, when the initialization of the RVE is approximate to the optimal result, the iteration number can be reduced, as presented in Table 18.
It should be noted that the design variables are defined only on the CPs of the B-spline mesh in this work, while Montemurro and his co-workers further included the weights among the design variables for the NURBS-based TO. [22][23][24][25][26][27][28] Structural compliance, the displacements, the eigen-value problems, and the dynamic responses were comprehensively investigated for 2D and 3D problems. The influences of the main NURBS parameters (weights, degrees, and number of CPs) on the optimal results were presented in detail. Compared to the B-splines based solution, a high-quality solution in terms of the objective value was achieved with a smoother boundary when the weights were included among the design variables. Furthermore, when the B-spline entities were adopted to describe the density distribution of the RVE, the greater the number of CPs (for a given degree) or the bigger the degree (for a given number of CPs) the better convergence solution can be obtained in terms of the objective value. Meanwhile, when the NURBS entities were used to describe the RVE topological variables, similar conclusion can be drawn and the shape and topology seem to approximately satisfy the global result observed from their works. The higher the number of the CPs, the smaller the filter size and the improved performances can be achieved. In this work, the issues not investigated here are the influences of the main NURBS parameters, which will be a further topic in future work, aiming at the application and extension of the works of Montemurro and his co-workers. [22][23][24][25][26][27][28] F I G U R E 33 Concurrent optimal results of 3D cantilever beam starting from Figure 32A

CONCLUSIONS
This article presents the Bézier extraction based IGA framework for multi-objective TO of the RVE and CTO of the composite structure with periodic microstructures. The density-based optimization model defines the design variables on the CPs of the B-spline mesh, and the PBCs are imposed on the B-spline mesh for the numerical homogenization to determine the effective macroscopic properties. The sensitivity analyses of the mentioned problems are derived with respect to the design variables on the CPs of the B-spline mesh. The Bézier elements based IGA facilitates the incorporation of the existing FE and TO procedures. 58,59 The proposed Bézier elements based IGA for TO also preserves the local support property of the traditional IGA, which generates optimal results towards the black-and-white solutions. Furthermore, the distance-based density distribution is proposed as the design initialization, which avoids to guess the initialization with special shape and parameters, such as the circular region and its radius. Numerical examples, including the bulk modulus maximization, NPR and CTO of macrostructure with periodic microstructures, show the robustness and effectiveness of the distance-based density distribution as the initialization for designing the composites with periodic microstructures.