Isogeometric topology optimization based on energy penalization for symmetric structure

Abstract

We present an energy penalization method for isogeometric topology optimization using moving morphable components (ITO–MMC), propose an ITO–MMC with an additional bilateral or periodic symmetric constraint for symmetric structures, and then extend the proposed energy penalization method to an ITO–MMC with a symmetric constraint. The energy penalization method can solve the problems of numerical instability and convergence for the ITO–MMC and the ITO–MMC subjected to the structural symmetric constraint with asymmetric loads. Topology optimization problems of asymmetric, bilateral symmetric, and periodic symmetric structures are discussed to validate the effectiveness of the proposed energy penalization approach. Compared with the conventional ITO–MMC, the energy penalization method for the ITO–MMC can improve the convergence rate from 18.6% to 44.5% for the optimization of the asymmetric structure. For the ITO–MMC under a bilateral symmetric constraint, the proposed method can reduce the objective value by 5.6% and obtain a final optimized topology that has a clear boundary with decreased iterations. For the ITO–MMC under a periodic symmetric constraint, the proposed energy penalization method can dramatically reduce the number of iterations and obtain a speedup of more than 2.

Key words： topology optimization; moving morphable component; isogeometric analysis; energy penalization method; symmetric constraint

Introduction

Topology optimization (TO) from the work of Bendsøe and Kikuchi [¹] has been used extensively in the academic community due to its effectiveness in obtaining efficient and lightweight structures. With the endeavors of many excellent scientific researchers over past three decades [²–⁵], many methods, such as homogenization approach [¹], solid isotropic material with penalization (SIMP) [⁶–¹¹], level set method (LSM) [¹²–²²], and the evolutionary approach (ESO) [²³,²⁴], have been presented to solve TO problems. From the viewpoint of application, TO has been applied to many industrial fields owing to its tremendous superiorities in designing optimal structures, such as aircraft [²⁵,²⁶], additive manufacturing [²⁷,²⁸], photonics design [²⁹,³⁰], and hip implant design [³¹]. A remarkable work of a full-scale airplane wing with the optimal design of the internal structure has been elaborated in the journal Nature [³²].

Large engineering parts, such as airplane wing and car wheel, are often composed of small symmetric components for the benefit of stress equilibrium, manufacture convenience, structure stability, transportation, installation, maintenance, and recycling [¹⁸]. Many researchers have investigated the TO of symmetric structures. Sigmund [³³] proposed an inverse homogenization method based on periodic boundary conditions to tailor the effective properties of cellular materials possessing periodic symmetric microstructures. With this method, the optimal cellular material distribution is obtained considering extreme macroscopic mechanical properties. Afterward, the layout of macrostructures considering periodic symmetric geometries is investigated through the macroscopic design to identify the preliminary layout of materials, where the material distribution of microstructure is locally determined using a refined design [³⁴]. Using bi-directional ESO (BESO) and with the aim of minimization of the pressure frequency response, the optimal material distribution within the design domain can be generated for an acoustic-structure coupled system with periodic symmetric geometry constraint [³⁵]. Unit cells having different geometries and discretized into irregular finite element meshes for periodic structures are optimized by a TO method proposed in Ref. [³⁶] based on the BESO and Shepard interpolation functions.

However, the lack of geometric information about the boundary of the optimized results is a serious drawback in the optimal design of symmetric structures using the traditional TO approaches. This drawback can be resolved by the explicit TO method based on moving morphable component (MMC) [³⁷]. Compared with the conventional TO approaches, TO approaches based on MMC and moving morphable void (MMV) can directly obtain the explicit geometric information from the final result. Some imperfections existing in conventional TO methods based on MMC and MMV may cause failure of structural optimization. First, unstable results may occur due to the numerical instabilities resulting from the use of linear elements in finite element analysis (FEA) [³⁸]. Second, the computational burden is huge when the order of elements used in the FEA is no less than quadratic [¹⁶]. Last, low-order elements lead to the checkerboard pattern in the final result, which can be circumvented using high-order elements [⁶]. Therefore, high-efficient MMC–TO must be developed for symmetric structures to maximize the aforementioned effectiveness of MMC–TO in the structural design for the symmetric structure and to overcome the drawbacks of the traditional TO method in structural optimization.

The problems mentioned above in TO can be dealt with using the isogeometric analysis (IGA) proposed by Hughes et al. [³⁹], which features high accuracy and computational efficiency when high-order elements are used in the analysis. Given the aforementioned advantages, IGA has attracted more and more academic researchers [⁴⁰–⁴⁸]. In recent years, using IGA to replace finite element method (FEM) in structural TO has received increasing attention. An isogeometric TO (ITO) combining with a phase-field model is presented in Ref. [⁴⁹], and higher accuracy is obtained with the help of the embedded CAD geometric expression in optimization and analysis. Shape irregularities, such as “islanding” and “layering” are solved by taking B-spline as the shape function to obtain the corresponding shape boundaries exploiting the high smoothness of B-spline [⁵⁰]. An accurate and efficient ITO is proposed by utilizing the advantages of IGA and parameterized LSM (PLSM) [¹⁶], and an arbitrary geometric constraint is added in such ITO scheme to satisfy some practical requirements [¹⁵]. Furthermore, the TO problem of isotropic and anisotropic lattice materials is settled by an ITO scheme with the aid of asymptotic homogenization theory [⁵¹]. A parallel computation approach based on GPU has been used to overcome the tremendous computational burden of the level set-based ITO [¹⁷]. Xu et al. [⁵²] developed an isogeometric approach to the TO of spatially graded hierarchical structures. FEM is replaced by IGA in MMC–TO (ITO–MMC), which inherits the merits of MMC–TO and IGA, to promote the stability and robustness of MMC–TO [³⁸]. Xie et al. [⁵³] introduced a new ITO–MMC based on R-functions and Greville abscissae collocation scheme to accelerate the convergence rate of MMC–TO and obtain an effective constitutive model. However, the ITO–MMC may encounter slow convergence for some TO problems, such as the short beam subjected to an asymmetric load shown in Section 6.1.

To solve this problem, we propose an energy penalization method for ITO–MMC. With this method, the sensitivities of MMC geometric parameters with respect to the objective function are modified by the penalization of the strain energy of elements. In addition, an ITO–MMC is applied to the design of symmetric structures to obtain explicit structural boundary information. Then, the energy penalization method is extended to symmetric structures under asymmetric loads optimized by an ITO–MMC with a symmetric constraint. In the remainder of our work, the basic theory of MMC–TO, the non-uniform rational B-splines (NURBS) basis function used in IGA and the mathematic model of ITO–MMC are given in Section 2. Then, an energy penalization method for ITO–MMC is proposed in Section 3. With the consideration of two types of symmetric constraints, an optimization model based on ITO–MMC is displayed in Section 4. Section 5 describes the sensitivities of MMC geometric parameters under the symmetric constraint based on the extension of the energy penalization method. Section 6 expounds on the numerical results of several symmetric TO problems in terms of the effectiveness of the ITO–MMC with a symmetric constraint and the energy penalization method. Finally, conclusions are provided in Section 7.

Fundamentals

A summary of the MMC–TO

In MMC–TO, topology description functions (TDFs) [⁵⁴] constructed from the MMC geometric parameters describe the relationship between the element Young’s modulus and the geometric variables of MMCs. For the sake of simplicity, MMCs with straight skeleton and variable depth are used in this paper, which have the same geometric parameters as those illustrated in Fig. 1.

Fig.1 Illustration of the ith component, which has a straight skeleton and quadratically varies along the direction of depth.

Full size|PPT slide

The TDF of the component with the same shape as Fig. 1 can be expressed by a set of geometric parameters explicitly, including the position of the center (x_Oi, y_Oi) in the global coordinate system, the half-length (L_i) of the component, the inclined angle with respect to the x-axis (q_i), and the variable depth (t_1i, t_2i, t_3i). The formula of TDF corresponding to the ith structural component is shown in Eq. (1):

(1){\varphi }_i\left(x,y\right)={\left({\displaystyle \frac{x^{\prime }}{L_i}}\right)}^s\hspace{0.17em}+\hspace{0.17em}{\left({\displaystyle \frac{y^{\prime }}{f\left(x^{\prime }\right)}}\right)}^s\hspace{0.17em}-\mathrm{1},

with

(2)\left\{\begin{array}{l}{x}^{\prime }\hfill \\ y^{\prime }\hfill \end{array}\right\}\hspace{0.17em}=\hspace{0.17em}\left[\begin{array}{cc}\hfill \cos {\theta }_i\hfill & \sin {\theta }_i\hfill \\ -\sin {\theta }_i\hfill & \cos {\theta }_i\hfill \end{array}\right]\left\{\begin{array}{c}x-x_{Oi}\hfill \\ y-y_{Oi}\hfill \end{array}\right\},

and

(3)f\left(x^{\prime }\right)={\displaystyle \frac{t_{\mathrm{1}i}+t_{\mathrm{2}i}-\mathrm{2}{t}_{\mathrm{3}i}}{\mathrm{2}{L}_i^{\mathrm{2}}}}{\left(x^{\prime }\right)}^{\mathrm{2}}+{\displaystyle \frac{t_{\mathrm{2}i}-t_{\mathrm{1}i}}{\mathrm{2}{L}_i}}{x}^{\prime }+t_{\mathrm{3}i}.

In Eq. (1), s is an even integer (s=\mathrm{6} herein), and (x,y) represents one of the nodal coordinates of linear FEM elements used to discretize the design domain. Each individual component has a unique TDF such that the number of TDF values of all the element nodal points is equal to the that of components. If the structural topology of the design domain is composed of n structural components, then the implicit expression of structural topology is formulated as Eq. (4):

(4)

where

(5){\varphi }^s\left(x\right)=\max \left({\varphi }_{\mathrm{1}},{\varphi }_{\mathrm{2}},\mathrm{...},{\varphi }_n\right).

In Eq. (4), D denotes the design domain, and {\Omega }^s represents a proportion of the designable domain occupied with n solid components.

The general theoretical model of the framework of TO using MMC can be written as follows [³⁷]:

(6)\begin{array}{l}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{Find}\hspace{0.17em}\hspace{0.17em}d={\left({\mathit{d}}_{\mathrm{1}}^{\mathrm{T}},{\mathit{d}}_{\mathrm{2}}^{\mathrm{T}},\mathrm{...},{\mathit{d}}_n^{\mathrm{T}}\right)}^{\mathrm{T}},\hfill \\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{Minimize}\hspace{0.17em}I=I\left(\mathit{d}\right),\hspace{0.17em}\hfill \\ \mathrm{s}\mathrm{.t}.\hspace{0.17em}\hspace{0.17em}\hfill \\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{g}_j\left(\mathit{d}\right)\hspace{0.17em}\le \mathrm{0},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}j=\mathrm{1},\mathrm{2},\mathrm{...},m,\hfill \\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathit{d}\subset {\mathit{U}}_{\mathit{d}},\hfill \end{array}

where {\mathit{d}}_i=(x_{Oi},y_{Oi},L_i,t_{\mathrm{1}i},t_{\mathrm{2}i},t_{\mathrm{3}i},{\theta }_i) is the union of geometric variables of ith components, the geometric variables of all components d determine the objective functional I, m represents the number of inequality constraints, and the admissible set of d is expressed by {\mathit{U}}_{\mathit{d}}.

With the aid of Eulerian mesh description and fixed finite element mesh as well as the values of TDF as defined in Eq. (5), the problem formulation shown in Eq. (6) can be reformulated as Eq. (7) when the minimization of structural compliance is the objective of TO under a prescribed volume constraint.

(7)\begin{array}{l}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{Find}\hspace{0.17em}\hspace{0.17em}\mathit{d}={\left({\mathit{d}}_{\mathrm{1}}^{\mathrm{T}},{\mathit{d}}_{\mathrm{2}}^{\mathrm{T}},\mathrm{...},{\mathit{d}}_n^{\mathrm{T}}\right)}^{\mathrm{T}},\hfill \\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{Minimize}\hspace{0.17em}I=I\left(\mathit{d},\mathit{u}\right),\hspace{0.17em}\hfill \\ \mathrm{s}\mathrm{.t}.\hspace{0.17em}\hspace{0.17em}\hfill \\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\displaystyle {\int }_{D}H\left({\varphi }^s\right)\hspace{0.17em}}\mathbb{E}:\mathit{\epsilon }\left(\mathit{u}\right):\mathit{\epsilon }\left(\mathit{v}\right)\mathrm{d}V\hfill \\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}={\displaystyle {\int }_{D}H\left({\varphi }^s\right)\hspace{0.17em}}\mathit{f}·\mathit{v}\mathrm{d}V+\hspace{0.17em}{\displaystyle {\int }_{{\mathrm{\Gamma }}_t}\mathit{t}·\mathit{v}\mathrm{d}S},\hspace{0.17em}\forall \mathit{v}\in {\mathit{U}}_{\mathrm{ad}},\hfill \\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\displaystyle {\int }_{D}H\left({\varphi }^s\right)\hspace{0.17em}\mathrm{d}V}\le \overline{V},\hspace{0.17em}\hfill \\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathit{d}\subset {\mathit{U}}_{\mathit{d}},\hfill \\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathit{u}=\overline{\mathit{u}},\hspace{0.17em}\hspace{0.17em}\mathrm{on}\hspace{0.17em}{\mathrm{\Gamma }}_{\mathit{u}},\hfill \end{array}

where H denotes the Heaviside function, \mathbb{E}=\mathbb{E}[\mathrm{I/m:mi>} (\mathbb{I} is the fourth-order identity tensor, \mathit{\delta } represents the second-order identity tensor, and E and v are the Young’s modulus and Poisson’s ratio, respectively) is the fourth-order isotropic elasticity tensor of the material, \mathit{f} is the external force, and t represents the surface traction on Neumann boundary {\mathrm{\Gamma }}_{\mathit{t}}. u is the displacement field, and v denotes the test function corresponding to u with {\mathit{U}}_{\mathrm{a}\mathrm{d}}=\{\mathit{v}|\mathit{v}\in {\mathrm{H}}^{\mathrm{1}}(\mathrm{\Omega }),\hspace{0.17em}\mathit{v}=\mathrm{0}\mathrm{o}\mathrm{n}{\mathrm{\Gamma }}_u\}. ε represents the second-order linear strain tensor. Moreover, \overline{V} is the prescribed solid material volume, and \overline{\mathit{u}} is the fixed displacement on Dirichlet boundary {\mathrm{\Gamma }}_u.

Basic theory of NURBS

NURBS [⁵⁵,⁵⁶] is a vital ingredient of numerical discretization in IGA. Having n spline basis functions of p-degree, the (ordered) knot vector \mathbf{\Xi }=\{{\xi }_{\mathrm{1}},{\xi }_{\mathrm{2}},\mathrm{...},{\xi }_{n+p+\mathrm{1}}\} denotes a non-decreasing sequence of real numbers representing parametric coordinates in 1D parameter space. The interval [{\xi }_{\mathrm{1}},{\xi }_{n+p+\mathrm{1}}]\hspace{0.17em}represents a patch, and the subinterval [{\xi }_i,{\xi }_{i+\mathrm{1}}] denotes a span. The B-spline basis functions of a knot vector \mathbf{\Xi } are recursively derived from the celebrated Cox-de Boor recursion formula [⁵⁷]:

B_{i,p}\left(\xi \right)={\displaystyle \frac{\xi -{\xi }_i}{{\xi }_{i+p}-{\xi }_i}}{B}_{i,p-\mathrm{1}}\left(\xi \right)+{\displaystyle \frac{{\xi }_{i+p+\mathrm{1}}-\xi }{{\xi }_{i+p+\mathrm{1}}-{\xi }_{i+\mathrm{1}}}}{B}_{i+\mathrm{1},p-\mathrm{1}}\left(\xi \right),

(8)

where the convention \mathrm{0}/\mathrm{0}=\mathrm{0} is used. Figure 2 depicts seven cubic basis functions constructed from an open knot vector.

Fig.2 B-spline basis functions for open knot vector

\mathbf{\Xi }=\left\{\mathrm{0},\mathrm{0},\mathrm{0},\mathrm{0},\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{3},\mathrm{4},\mathrm{4},\mathrm{4},\mathrm{4}\right\}

.

Full size|PPT slide

A tensor product is used to construct the basis functions for multi-dimensional parameter space, formulated as

(9)N_{i,p}^{j,q}(\xi ,\eta )=N_{i,p}\left(\xi \right)N_{j,q}\left(\eta \right),

where N_{i,p}\left(\xi \right) is the B-spline basis function of the knot vector \mathbf{\Xi }=\left\{{\xi }_{\mathrm{1}},{\xi }_{\mathrm{2}},\mathrm{...},{\xi }_{n+p+\mathrm{1}}\right\}, and N_{j,q}\left(\eta \right) represents the B-spline basis function for \Pi =\left\{{\eta }_{\mathrm{1}},{\eta }_{\mathrm{2}},\mathrm{...},{\eta }_{m+q+\mathrm{1}}\right\}. A NURBS basis function is obtained from scaling a B-spline basis function with a corresponding positive weight w_i:

(10)N_{i,p}\left(\xi \right)={\displaystyle \frac{B_{i,p}\left(\xi \right)w_i}{{\displaystyle \mathrm{\sum }_{j=\mathrm{1}}^n{B}_{j,p}\left(\xi \right)w_j}}}.

By means of the tensor product, NURBS basis functions in 2D parameter space are as follows:

(11)N_{i,p}^{j,q}\left(\xi ,\eta \right)={\displaystyle \frac{B_{i,p}\left(\xi \right)B_{j,q}\left(\eta \right)w_{i,j}}{{\displaystyle \mathrm{\sum }_{k=\mathrm{1}}^n{\displaystyle \mathrm{\sum }_{t=\mathrm{1}}^m{B}_{k,p}\left(\xi \right)B_{t,q}\left(\eta \right)w_{k,t}}}}},

where w_{i,j} denotes the weight of the tensor product B_{i,p}\left(\xi \right)B_{j,q}\left(\eta \right).

R-functions and Greville abscissae-based ITO–MMC

In the MMC–TO method and ITO–MMC by Hou et al. [³⁸], the values of TDF are calculated by max function which leads to C¹ discontinuity. However, the C¹ discontinuity caused by max function obviously deteriorates the convergence of TO because the method of moving asymptotes (MMA) is used in the MMC–TO and ITO–MMC as the optimizer, where the variables are optimized by the gradient. The TDF values of overlapping regions of components are computed by R-functions [⁵⁸] in Greville abscissae-based ITO–MMC to solve the problem [⁵³]. Thus, the formula of TDF illustrated in Eq. (5) is replaced by Eq. (12):

{\varphi }^s=\{\begin{array}{ll}{R}^{\alpha }\left(\cup {\varphi }_{i=\mathrm{1}}^m\right),\hfill & \mathrm{if}\hspace{0.17em}{\varphi }_i\in {\varphi }_{+}\hspace{0.17em}\mathrm{with}\hspace{0.17em}\hspace{0.17em}\mathrm{2}\le m\le n,\hfill \\ \hspace{0.17em}\mathrm{m}\mathrm{a}\mathrm{x}\left({\varphi }_{\mathrm{1}},{\varphi }_{\mathrm{2}},\mathrm{...},{\varphi }_n\right),\hfill & \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e},\hfill \end{array}

\hspace{0.17em}\mathrm{with}\hspace{0.17em}\hspace{0.17em}{\varphi }_{+}=\left\{{\varphi }_i\in \{{\varphi }_{\mathrm{1}},{\varphi }_{\mathrm{2}},\mathrm{...},{\varphi }_n\}|{\varphi }_i>\mathrm{0}\right\},

$$\begin{gathered} R^{\alpha }\left(\cup {\varphi }_{i=\mathrm{1}}^m\right)\hspace{0.17em}={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+\alpha }}(R^{\alpha }\left(\cup {\varphi }_{i=\mathrm{1}}^{m-\mathrm{1}}\right)+{\varphi }_m \\ +\sqrt{{\left(R^{\alpha }\left(\cup {\varphi }_{i=\mathrm{1}}^{m-\mathrm{1}}\right)\right)}^{\mathrm{2}}+{\varphi }_m^{\mathrm{2}}-\mathrm{2}\alpha R^{\alpha }\left(\cup {\varphi }_{i=\mathrm{1}}^{m-\mathrm{1}}\right){\varphi }_m})\hspace{0.17em},\hspace{0.17em}\hspace{0.17em} \end{gathered}$$

(12)\alpha \in \left(-\mathrm{1},\mathrm{1}\right).

In Eq. (12), {\varphi }_{+} represents a subset of \left\{{\varphi }_{\mathrm{1}},{\varphi }_{\mathrm{2}},\mathrm{...},{\varphi }_n\right\} and the elements of the set {\varphi }_{+} (the number of elements is m) are positive values. In addition, R^{\alpha }\left(\cup {\varphi }_{i=\mathrm{1}}^m\right) is a straightforward extension of Eq. (13):

R^{\alpha }\left({\varphi }_{\mathrm{1}}\cup {\varphi }_{\mathrm{2}}\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+\alpha }}\left({\varphi }_{\mathrm{1}}+{\varphi }_{\mathrm{2}}+\sqrt{{\varphi }_{\mathrm{1}}^{\mathrm{2}}+{\phi }_{\mathrm{2}}^{\mathrm{2}}-\mathrm{2}\alpha {\varphi }_{\mathrm{1}}{\varphi }_{\mathrm{2}}}\right),\hspace{0.17em}

(13)\alpha \in (-\mathrm{1},\mathrm{1}).

To solve the problem of numerical instability resulting from linear FEM elements used in traditional MMC–TO and the high computation cost of high-order FEM elements illustrated in Fig. 3, Xie et al. [⁵³] used IGA in ITO–MMC to replace FEM. The Greville abscissae collocation scheme is adopted to establish the ersatz material model in the new ITO–MMC, which is different from the indirect and complicated ersatz material model used in ITO–MMC by Hou et al. [³⁸]. Hence, the ITO–MMC formula on the basis of R-functions and Greville abscissae for the minimization compliance problem under a specified volume constraint is

(14)\begin{array}{l}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{Find}\hspace{0.17em}\mathit{d}={\left({\mathit{d}}_{\mathrm{1}}^{\mathrm{T}},{\mathit{d}}_{\mathrm{2}}^{\mathrm{T}},\mathrm{...},{\mathit{d}}_n^{\mathrm{T}}\right)}^{\mathrm{T}},\hspace{0.17em}\hfill \\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{Minimize}\hspace{0.17em}C={\displaystyle \mathrm{\sum }_{e=\mathrm{1}}^{N_e}{\mathit{\gamma }}_e{\mathit{u}}_e^{\mathrm{T}}}{\mathit{K}}_e{\mathit{u}}_e,\hspace{0.17em}\hfill \\ \mathrm{s}\mathrm{.t}.\hspace{0.17em}\hspace{0.17em}\hfill \\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\left({\displaystyle \mathrm{\sum }_{e=\mathrm{1}}^{N_e}{\mathit{\gamma }}_e}{\mathit{K}}_e\right)·\mathit{u}=\mathit{F},\hfill \\ \\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\displaystyle \mathrm{\sum }_{e=\mathrm{1}}^{N_e}{\rho }_e}{V}_e\hspace{0.17em}\le \overline{V},\hfill \\ \\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathit{d}\subset {\mathit{U}}_{\mathit{d}},\hfill \end{array}

where \mathit{d} is the union of the geometric variables, \mathit{u} denotes the displacement field, C represents the structural mean compliance, {\mathit{\gamma }}_e and {\rho }_e denote the eth element Young’s modulus coefficient and density calculated by Eq. (15), respectively. {\mathit{u}}_e is the displacement vector of the eth element, {\mathit{K}}_e is the standard stiffness matrix of the eth element, \mathit{F} represents the external load vector, V_e is the eth element volume, \overline{V} represents the specified volume ratio, and {\mathit{U}}_{\mathit{d}} denotes the admissible set of \mathit{d}.

(15)\{\begin{array}{c}\mathit{\gamma }=\frac{1}{n_{col}}\left(\mathrm{\sum }_{icol=1}^{n_{col}}{H}_{\epsilon }^q\left({\varphi }_{icol}\right)\right),\hfill \\ \rho =\frac{1}{n_{col}}\left(\mathrm{\sum }_{icol=1}^{n_{col}}{H}_{\epsilon }\left({\varphi }_{icol}\right)\right),\hfill \end{array}

where \mathit{\gamma } and \rho represent Young’s modulus coefficient and density based on Greville abscissae collocation points, n_{col} denotes the number of element collocation points, {\varphi }_{icol} denotes the value of TDF of the icolth collocation point calculated by Eq. (12), q=\mathrm{2}, H_{\epsilon }^q\left({\varphi }_{icol}\right) is the component of the Young’s modulus coefficient at the icolth collocation point, and H_{\epsilon }\left({\varphi }_{icol}\right) obtained from Eq. (16) represents the regularized TDF value of icolth collocation point:

(16)

where \beta denotes a positive parameter to assure the non-singularity of the global stiffness matrix, which is very small, and \epsilon controls the level of regularization.

Fig.3 A comparison: Total number of DOFs is 162 in the left for FEA, which are much larger than 72 for IGA.

Full size|PPT slide

Energy penalization method for ITO--MMC

In MMC–TO, the sensitivity of the minimization compliance problem is obtained from the first-order variation of the mean compliance to the variation of component geometric parameters. When an arbitrary design variable a_j is changed, the corresponding sensitivity of the mean compliance in ITO–MMC is approximately defined as [⁵³]

(17){\displaystyle \frac{\partial C}{\partial a_j}}=-\mathrm{\sum }_{e=1}^{N_e}[\frac{{{\mathit{u}}_e}^{\mathrm{T}}q}{n_{col}}{\left(\mathrm{\sum }_{icol=1}^{n_{col}}{H}_{\epsilon }{\left({\varphi }_{icol}\right)}^{q-1}\frac{\partial H_{\epsilon }\left({\varphi }_{icol}\right)}{\partial a_j}\right)}_e·{\mathit{K}}_e{\mathit{u}}_e],

where N_e is the number of NURBS elements used to obtain the displacement field of the design domain, n_{col} denotes the number of Greville abscissae collocation points of a NURBS element, and it takes q=\mathrm{2}, and {\displaystyle \frac{\partial H_{\epsilon }\left({\varphi }_{icol}^i\right)}{\partial a_j}} represents the finite difference quotient of {\varphi }_{icol}^i and the design variable a_j.

The sensitivity of the volume constraint function in ITO–MMC is illustrated in Eq. (18):

(18){\displaystyle \frac{\partial V}{\partial a_j}}={\displaystyle \mathrm{\sum }_{e=\mathrm{1}}^{N_e}\left[{\displaystyle \frac{\mathrm{1}}{n_{col}}}{\left({\displaystyle \mathrm{\sum }_{icol=\mathrm{1}}^{n_{col}}\frac{\partial H\left({\varphi }_{icol}\right)}{\partial a_j}}\right)}_e{V}_e\right]}\hspace{0.17em}.

The original ITO–MMC procedure is straightforwardly proceeded with gradually optimizing the geometric parameters of components by the MMA [⁵⁹].

During the above objective functional sensitivity analysis, the sensitivities for geometric design variables depend on the strain energy of elements. If the differences of strain energy of elements are overly large, which usually occurs in the optimization of a structure subjected to an asymmetric load, it may lead to the very limited effect of the elements with small strain energy on the optimization, which will deteriorate the convergence and numerical stability of the optimization process. To solve the aforementioned problem, we propose an energy penalization method based on the strain energy obtained by the p-power of strain energy of each element in this part. The sensitivity with respect to the mean compliance of component geometric parameters is reformulated as follows:

(19)

The basic idea of the penalization method is originated from the variation curves of a power function y=x^p taking different p values, and the curves are shown in the left of Fig. 4. The curves depicted in Fig. 4(a) show that the dependent variable y is decreased and a penalization effect occurs in the dependent variable y, when p is decreased and the variable x is set to a fixed value. In addition, the difference between 10^p and 2^p is decreased with the decrease in p-value on the basis of the variation curve shown in Fig. 4(b). Although the strain energy of elements in the void region will be enlarged, the void elements’ magnifying effect on the sensitivity calculated by Eq. (19) is very limited because of the small Heaviside values of Greville abscissae of the void elements. Hence, the large difference in the strain energy of solid elements can be eased by reducing the value of p used as the power exponent of the strain energy, for the structure optimized by ITO–MMC. With the modified sensitivity, the variations of the design variables are reduced over the optimization process, so a converged solution is easily obtained.

Fig.4 Variation curves of (a) y=x^p with different p values and (b) the difference between 10^p and 2^p with respect to p.

Full size|PPT slide

ITO–MMC with a symmetric constraint

The target of the present TO is to generate an optimal material distribution within the design domain simultaneously subjected to volume and symmetric constraints. For the sake of simplicity, minimization of the mean compliance is only presented in this paper. The symmetric constraint has two concrete forms, namely, bilateral and periodic symmetric constraints. Two TO models based on ITO–MMC briefly introduced in Section 2.3 are developed in this work. One is for TO problems with bilateral symmetric constraint, and the other is for TO problems with periodic symmetric constraint.

Bilateral symmetric TO problems using ITO–MMC

Huang and Xie [⁶⁰] proposed a method of applying the periodic constraint to a structure based on BESO, where the densities of periodic symmetric elements are equal and only the densities of elements in a unit cell are the design variables. Similarly, two bilateral symmetric unit cells construct the design domain as elaborated in Fig. 5 to apply the bilateral symmetric constraint in a TO problem. Moreover, the basic structural components in the first and the second unit cells are bilateral symmetric with respect to the dashed symmetry axis presented in Fig. 5. Hence, the basic structural components in any unit cell can be used as the design variables, which reduce the design variables by half than the ITO–MMC without a bilateral symmetric constraint. According to the ersatz material model shown in Eq. (15), the Young’s modulus coefficient and density of each element in the first cell are obtained, whereas that information in the second cell can be easily gained from the first cell according to the bilateral symmetric constraint. The total number of basic structural components and the initial parameters corresponding to those components in the second unit cell can be specified from the first unit cell. Thus, to take the bilateral symmetric constraint into consideration, we can reformulate the minimization problem of mean compliance in terms of the ersatz material model generated from the geometric design variables in the first cell (Fig. 5), as

(20)\begin{array}{l}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{Find}\hspace{0.17em}\mathit{d}={\left({\mathit{d}}_{\mathrm{1}}^{\mathrm{T}},{\mathit{d}}_{\mathrm{2}}^{\mathrm{T}},\mathrm{...},{\mathit{d}}_{n/\mathrm{2}}^{\mathrm{T}}\right)}^{\mathrm{T}},\hfill \\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{Minimize}\hspace{0.17em}C={\displaystyle \mathrm{\sum }_{i=\mathrm{1}}^{nelx}{\displaystyle \mathrm{\sum }_{j=\mathrm{1}}^{nely}{\mathit{\gamma }}_{i,j}}{\mathit{u}}_{i,j}^{\mathrm{T}}}{\mathit{K}}_{i,j}{\mathit{u}}_{i,j},\hfill \\ \mathrm{s}\mathrm{.t}.\hspace{0.17em}\hspace{0.17em}\hfill \\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\left({\displaystyle \mathrm{\sum }_{i=\mathrm{1}}^{nelx}{\displaystyle \mathrm{\sum }_{j=\mathrm{1}}^{nely}{\mathit{\gamma }}_{i,j}}}{\mathit{K}}_{i,j}\right)·\mathit{u}=\mathit{F},\hfill \\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\displaystyle \mathrm{\sum }_{i=\mathrm{1}}^{nelx}{\displaystyle \mathrm{\sum }_{j=\mathrm{1}}^{nely}{\rho }_{i,j}}}{V}_{i,j}\hspace{0.17em}\le \overline{V},\hfill \\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\mathit{\gamma }}_{nelx-i+\mathrm{1},j}={\mathit{\gamma }}_{i,j},\hspace{0.17em}\left(i=\mathrm{1},\hspace{0.17em}\mathrm{2},\hspace{0.17em}\mathrm{...},\hspace{0.17em}nelx/\mathrm{2}\mathrm{;}\hspace{0.17em}j=\mathrm{1},\hspace{0.17em}\mathrm{2},\hspace{0.17em}\mathrm{...},\hspace{0.17em}nely\right),\hfill \\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\rho }_{nelx-i+\mathrm{1},j}={\rho }_{i,j},\hspace{0.17em}\left(i=\mathrm{1},\hspace{0.17em}\mathrm{2},\mathrm{...},\hspace{0.17em}nelx/\mathrm{2}\mathrm{;}\hspace{0.17em}j=\mathrm{1},\hspace{0.17em}\mathrm{2},\hspace{0.17em}\mathrm{...},\hspace{0.17em}nely\right),\hfill \\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathit{d}\subset {\mathit{U}}_{\mathit{d}},\hspace{0.17em}\hfill \end{array}

where \mathit{d} is the set of design variables associated with the components in the first unit cell, \mathit{u} and \mathit{F} are the displacement vector of all control points of NURBS mesh and the applied external force, and C denotes the mean compliance. {\mathit{\gamma }}_{i\mathit{,}j}, {\rho }_{i,j}, and V_{i,j} correspond to the Young’s modulus coefficient, density, and the element volume (i and j represent the index of the element in the x and y directions). The bilateral symmetric constraint is satisfied by {\mathit{\gamma }}_{nelx-i+\mathrm{1},j}={\mathit{\gamma }}_{i,j} and {\rho }_{nelx-i+\mathrm{1},j}={\rho }_{i,j}.

Using the geometric parameters of basic structural components in the first or the second unit cell as the design variables shows no considerable difference.

Fig.5 2D design domain with bilateral symmetric constraint, where 12 basic structural components are only placed in the first unit cell, the design domain is divided into nelx and nely quadratic NURBS elements along the x and y directions, respectively, and the blue and yellow elements are bilateral symmetric, about the dark red symmetric axis.

Full size|PPT slide

Periodic symmetric TO problems using ITO–MMC

Similar to the problems under bilateral symmetric constraint in ITO–MMC, the TO problems subjected to periodic symmetric constraint can be solved by ITO–MMC. A 2D designable domain with m=\mathrm{8} unit cells is shown in Fig. 6 for the periodic symmetric case. For the simplicity of illustration, the parameters of components placed in the first unit cell are used as the design variables. The conventional TO problem appears when m=\mathrm{1}\times \mathrm{1}. Therefore, the TO problem subjected to periodic symmetric using R-functions and Greville abscissae-based ITO–MMC can be formulated as

(21)\begin{array}{l}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{Find}\hspace{0.17em}{\mathit{d}}_{\mathrm{1}}={\left({\mathit{d}}_{\mathrm{1},\mathrm{1}}^{\mathrm{T}},{\mathit{d}}_{\mathrm{1},\mathrm{2}}^{\mathrm{T}},\mathrm{...},{\mathit{d}}_{\mathrm{1},n/m}^{\mathrm{T}}\right)}^{\mathrm{T}},\hfill \\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{Minimize}C={\displaystyle \mathrm{\sum }_{i=\mathrm{1}}^m{\displaystyle \mathrm{\sum }_{j=\mathrm{1}}^{n_e}{\mathit{\gamma }}_{i,j}}{\mathit{u}}_{i,j}^{\mathrm{T}}}{\mathit{K}}_{i,j}{\mathit{u}}_{i,j},\hspace{0.17em}\hfill \\ \mathrm{s}\mathrm{.t}.\hspace{0.17em}\hspace{0.17em}\hfill \\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\left({\displaystyle \mathrm{\sum }_{i=\mathrm{1}}^m{\displaystyle \mathrm{\sum }_{j=\mathrm{1}}^{n_e}{\mathit{\gamma }}_{i,j}}}{\mathit{K}}_{i,j}\right)·\mathit{u}=\mathit{F},\hfill \\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\displaystyle \mathrm{\sum }_{i=\mathrm{1}}^m{\displaystyle \mathrm{\sum }_{j=\mathrm{1}}^{n_e}{\rho }_{i,j}}}{V}_{i,j}\hspace{0.17em}\le \overline{V},\hfill \\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\mathit{\gamma }}_{\mathrm{1},j}=\hspace{0.17em}{\mathit{\gamma }}_{\mathrm{2},j}=\cdots \hspace{0.17em}={\mathit{\gamma }}_{m,j},\hspace{0.17em}\left(j=\mathrm{1},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{2},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{...},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{n}_e\right),\hfill \\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\rho }_{\mathrm{1},j}=\hspace{0.17em}{\rho }_{\mathrm{2},j}=\cdots \hspace{0.17em}={\rho }_{m,j},\hspace{0.17em}\left(j=\mathrm{1},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{2},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{...},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{n}_e\right),\hfill \\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathit{d}\subset {\mathit{U}}_d,\hfill \end{array}

where {\mathit{d}}_{\mathrm{1}} denotes the set of all design variables that represent the geometric parameters of the structural components placed in the first unit cell, n/m denotes the number of components placed in the first unit cell, m is the number of unit cells used to decompose the design domain, and n_e represents the number of NURBS elements in a unit cell. {\mathit{\gamma }}_{i,j}, {\rho }_{i,j} and V_{i,j} represent the element Young’s modulus coefficient, density, and volume, respectively (i and j denote the index of a unit cell and the local element index within a unit cell, respectively).

Fig.6 2D design domain with m=8 unit cells, m₁ is the number of unit cells along the x direction, and m₂ is the number of unit cells along the y direction.

Full size|PPT slide

Numerical implementation for ITO–MMC with a symmetric constraint

Sensitivities of ITO–MMC with a symmetric constraint based on the energy penalization method

When an arbitrary geometric design variable a_m of components located in the first unit cell is changed, the corresponding sensitivity of the structural mean compliance for the ITO–MMC under a bilateral symmetric constraint illustrated in Fig. 5 is obtained by reformulating Eq. (17) as Eq. (22), and the sensitivity of the volume constraint Eq. (18) is transformed into Eq. (23), exploiting the equality of Young’s modulus coefficient defined by {\mathit{\gamma }}_{nelx-i+\mathrm{1},j}=\hspace{0.17em}{\mathit{\gamma }}_{i,j} and density formulated as {\rho }_{nelx-i+\mathrm{1},j}={\rho }_{i,j}.

(22){\displaystyle \frac{\partial C}{\partial a_m}}=-{\displaystyle \mathrm{\sum }_{i=\mathrm{1}}^{nelx/\mathrm{2}}{\displaystyle \mathrm{\sum }_{j=\mathrm{1}}^{nely}\left\{{\displaystyle \frac{q}{n_{col}}}{\left({\displaystyle \mathrm{\sum }_{icol=\mathrm{1}}^{n_{col}}{H}_{\epsilon }{\left({\varphi }_{icol}\right)}^{q-\mathrm{1}}{\displaystyle \frac{\partial H_{\epsilon }\left({\varphi }_{icol}\right)}{\partial a_m}}}\right)}_{i,j}\right({\mathit{u}}_{i,j}^{\mathrm{T}}{\mathit{K}}_{i,j}{\mathit{u}}_{i,j}+{\mathit{u}}_{k,j}^{\mathrm{T}}{\mathit{K}}_{k,j}{\mathit{u}}_{k,j}\left)\right\}\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}k=nelx-i+\mathrm{1},}}

(23){\displaystyle \frac{\partial V}{\partial a_m}}=-{\displaystyle \mathrm{\sum }_{i=\mathrm{1}}^{nelx/\mathrm{2}}{\displaystyle \mathrm{\sum }_{j=\mathrm{1}}^{nely}\left\{{\displaystyle \frac{\mathrm{1}}{n_{col}}}{\left({\displaystyle \mathrm{\sum }_{icol=\mathrm{1}}^{n_{col}}\frac{\partial H_{\epsilon }\left({\varphi }_{icol}\right)}{\partial a_m}}\right)}_{i,j}\left(V_{i,j}+V_{k,j}\right)\right\}}}\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}k=nelx-i+\mathrm{1.}

Due to the periodicity of the cells as those depicted in Fig. 6, the basic structural components constituting the topology of an arbitrary unit cell should be the same. In other words, the geometric parameters of the components that are periodic symmetric should have the same amount of change simultaneously. Therefore, we can choose the geometric parameters of components in any unit cell as the design variables of ITO–MMC. Thus, optimization can be implemented in a representative unit cell, which can be arbitrarily selected from all of the unit cells due to periodicity. The first unit cell is selected in this work. The computational domain of IGA is still the whole design domain of the TO problem.

According to Eq. (17), the sensitivity of the objective function caused by the variation of the arbitrary geometric design variable is approximately by the total strain energy of the Greville abscissae collocation points located in the corresponding components. The strain energy of the Greville abscissae collocation points is derived from the strain energy of the corresponding elements. In the ITO–MMC subjected to periodic symmetric constraint, the sensitivity of the representative cell is formulated as Eq. (24), and the sensitivity of the volume constraint can be obtained from Eq. (25), due to the repeatability of the components constituting the topology of an arbitrary unit cell.

(24){\displaystyle \frac{\partial C}{\partial a_j}}=-{\displaystyle \mathrm{\sum }_{e=\mathrm{1}}^{n_e}\left\{{\displaystyle \frac{q}{n_{col}}}{\left({\displaystyle \mathrm{\sum }_{icol=\mathrm{1}}^{n_{col}}{H}_{\epsilon }{\left({\varphi }_{icol}\right)}^{q-\mathrm{1}}{\displaystyle \frac{\partial H_{\epsilon }\left({\varphi }_{icol}\right)}{\partial a_j}}}\right)}_e{\displaystyle \mathrm{\sum }_{i=\mathrm{1}}^m{\mathit{u}}_{i,e}^{\mathrm{T}}{\mathit{K}}_{i,e}{\mathit{u}}_{i,e}}\right\},}

(25){\displaystyle \frac{\partial V}{\partial a_j}}=-{\displaystyle \mathrm{\sum }_{e=\mathrm{1}}^{n_e}\left\{{\left({\displaystyle \mathrm{\sum }_{icol=\mathrm{1}}^{n_{col}}\frac{\partial H_{\epsilon }\left({\varphi }_{icol}\right)}{\partial a_j}}\right)}_e{\displaystyle \mathrm{\sum }_{i=\mathrm{1}}^m{V}_{i,e}}\right\}}.

Similar to the optimization process of ITO–MMC without a symmetric constraint, the convergence and numerical stability of the optimization process are deteriorated by the large differences in strain energy of the symmetric elements for the ITO–MMC under a symmetric constraint. To tackle those problems, an extension of the energy penalization method elaborated in Section 3 is performed for the ITO–MMC under a symmetric constraint. Accordingly, Eqs. (22) and (24) are changed into Eqs. (26) and (27), respectively. Using the p-norm calculation with p<1, the strain energy of elements with small strain energy is enlarged, and the influence of elements with small strain energy can be magnified during the optimization of symmetric structures.

(26){\displaystyle \frac{\partial C}{\partial a_m}}=-{\displaystyle \mathrm{\sum }_{i=\mathrm{1}}^{nelx/\mathrm{2}}{\displaystyle \mathrm{\sum }_{j=\mathrm{1}}^{nely}\{{\displaystyle \frac{q}{n_{col}}}{\left({\displaystyle \mathrm{\sum }_{icol=\mathrm{1}}^{n_{col}}{H}_{\epsilon }{\left({\varphi }_{icol}\right)}^{q-\mathrm{1}}{\displaystyle \frac{\partial H_{\epsilon }\left({\varphi }_{icol}\right)}{\partial a_m}}}\right)}_{i,j}·{\left({\left({\mathit{u}}_{i,j}^{\mathrm{T}}{\mathit{K}}_{i,j}{\mathit{u}}_{i,j}\right)}^p+{\left({\mathit{u}}_{k,j}^{\mathrm{T}}{\mathit{K}}_{k,j}{\mathit{u}}_{k,j}\right)}^p\right)}^{\frac{\mathrm{1}}{p}}\}}}\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}k=nelx-i+\mathrm{1}\mathrm{a}\mathrm{n}\mathrm{d}p\le \mathrm{1},

(27){\displaystyle \frac{\partial C}{\partial a_j}}=-{\displaystyle \mathrm{\sum }_{e=\mathrm{1}}^{n_e}\{{\displaystyle \frac{q}{n_{col}}}{\left({\displaystyle \mathrm{\sum }_{icol=\mathrm{1}}^{n_{col}}{H}_{\epsilon }{\left({\varphi }_{icol}\right)}^{q-\mathrm{1}}{\displaystyle \frac{\partial H_{\epsilon }\left({\varphi }_{icol}\right)}{\partial a_j}}}\right)}_e·{\left({\displaystyle \mathrm{\sum }_{i=\mathrm{1}}^m{\left({\mathit{u}}_{i,e}^{\mathrm{T}}{\mathit{K}}_{i,e}{\mathit{u}}_{i,e}\right)}^p}\right)}^{\frac{\mathrm{1}}{p}}\}}\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}p\le \mathrm{1.}

The extension of the energy penalization method to the ITO–MMC with a symmetric constraint is illustrated by a simple case, where the design domain is decomposed into two symmetric unit cells, and the relationship between different strain energy are established in Eq. (28):

(28)\{\begin{array}{l}{C}_{\mathrm{total}}=C_{\mathrm{1}}^p+C_{\mathrm{2}}^p=C_{\mathrm{2}}^p\left({\alpha }^p+\mathrm{1}\right),\hfill \\ \alpha =C_{\mathrm{1}}/C_{\mathrm{2}},\hfill \\ {\theta }_{\mathrm{1}}=\hspace{0.17em}{C}_{\mathrm{1}}^p/C_{\mathrm{total}}={\alpha }^p/\left({\alpha }^p+\mathrm{1}\right),\hfill \\ {\theta }_{\mathrm{2}}=\hspace{0.17em}{C}_{\mathrm{2}}^p/C_{\mathrm{total}}=\mathrm{1}/\left({\alpha }^p+\mathrm{1}\right),\hfill \end{array}

where C₁ is the larger strain energy of element in the cell close to the load, C₂ is the smaller strain energy of element in the other cell, α represents the ratio of strain energy of the symmetric elements in different unit cells, q₁ and q₂ are respectively the proportion of the p power of the larger and smaller strain energy in the total strain energy C_total, and \alpha \ge 4 and p\le \mathrm{1}.

Basing from the variation curves illustrated in Fig. 7, we can conclude that for p<1 cases, q₁ is less than p=1 and q₂ is larger than p=1; q₁ is smaller if a smaller p value is used in Eq. (28); q₂ is larger if a smaller p value is used in Eq. (28). When α is fixed, the ratio of q₂ to q₁ is increased with the decrease in p. Therefore, for p<1 cases, the elements with smaller strain energy play a more important role in the optimization than the p=1 case due to the increase in q₂, the shrinking effect for the element in the cell close to load becomes more obvious with the decrease in p values, and the magnifying effect for the element in the other cell becomes more obvious with the decrease in p values.

Fig.7 Energy ratio variation curve for the design domain divided into two symmetric unit cells with different p values used in Eq. (28): Variation curve of the ratio (a) q₁=C₁/C_total, (b) q₂=C₂/C_total, and (c) q₂ to q₁.

Full size|PPT slide

Optimization for the ITO–MMC with a symmetric constraint

On the basis of the ITO–MMC by Xie et al. [⁵³], a new ITO–MMC is developed to solve TO problems with bilateral and periodic symmetric constraints. The flowchart for the ITO–MMC under a symmetric constraint is illustrated in Fig. 8. The meshing algorithm for the conventional FEA in MMC–TO is replaced by patch refinement [³⁹] in ITO–MMC to eliminate geometric discretization errors. MMC is used to obtain explicit boundary information from the TO results, whose geometric parameters are the design variables in ITO–MMC. The formula of the ITO–MMC with a symmetric constraint is developed in Section 4 to satisfy the symmetric constraint. The sensitivity of the mean compliance and the volume, as well as the sensitivities based on the energy penalization method for the ITO–MMC with a symmetric constraint, is shown in Section 5.1. Finally, the MMA algorithm is used as the optimizer in the ITO–MMC with a symmetric constraint.

Fig.8 Flowchart of ITO–MMC subjected to bilateral symmetric constraint (left) and periodic symmetric constraint (right).

Full size|PPT slide

Numerical examples

Five benchmark TO problems are studied in this part to verify the effectiveness of the proposed energy penalization method for ITO–MMC and the ITO–MMC for symmetric problems as well as the energy penalization method for the ITO–MMC with a symmetric constraint. The first one is intended to demonstrate the performance of the energy penalization method for asymmetric structure under ITO–MMC framework. The second and the third ones are used to testify the effectiveness of the proposed ITO–MMC for symmetric problems. Finally, the energy penalization method for the ITO–MMC considering a symmetric constraint is verified in the last two benchmarks. All physical quantities in the problems are dimensionless parameters due to the major concern of numerical performances. In any benchmark of our work, the Young’s modulus and the Poisson’s ratio used in standard elasticity matrix {\mathit{K}}_{i,j} formulated as Eqs. (20) and (21) are respectively set to 1 and 0.3 without special explanation. MMA [⁵⁹] is the optimizer used in the optimization process.

Asymmetric structure optimized by ITO–MMC based on the energy penalization method

The problem setting and the initial design of components for an asymmetric structure are illustrated in Fig. 9. The asymmetric structure is fixed at the left edge and subjected to a vertical load at the right-bottom corner, whereas the prescribed volume ratio is set to 0.4. The maximum number of iterations is 1000. Given that the mesh size and the parameters of MMA influence the convergence rate of MMC–TO, the ITO–MMC based on the energy penalization method is investigated for four cases (Table 1).

Tab.1 Parameters of MMA and the mesh size for the four cases

Case	Parameters of MMA						Mesh size of the NURBS mesh
	epsimin	raa0	albefa	asyinit	asyincr	asydecr
1	10–10	0.01	0.3	0.1	0.4	0.2	30 × 30
2	10–10	0.01	0.3	0.1	0.4	0.2	40 × 40
3	10–10	0.01	0.4	0.1	0.8	0.6	30 × 30
4	10–10	0.01	0.4	0.1	0.8	0.6	40 × 40

Fig.9 Problem setting and initial design of components for the asymmetric structure. (a) The problem setting for asymmetric structure; (b) the initial design of components for all cases.

Full size|PPT slide

The parameters of MMA used in Cases 1 and 2 are the same to those elaborated in Ref. [⁵³], and the parameters of MMA for Cases 3 and 4 are identical to those used in Ref. [⁶¹]. Then, the comparison results between the ITO–MMC based on the energy penalization method and the ITO–MMC for the four cases are shown in Fig. 10. p=1 means the ITO–MMC is not penalized by the energy penalization method.

As shown in Fig. 10, the improvement in convergence rate for all cases listed in Table 1 ranges from 18.6% to 44.5% through the energy penalization method used in ITO–MMC, without deteriorating the stiffness of the optimized structure. A large difference in the improvement effect of the energy penalization method on the convergence rate for different mesh sizes is observed. In specific, the improvement is no less than 40% for Cases 1 and 3 where a 30 × 30 mesh size is used while 20.2% for Cases 2 and 4 where a 40 × 40 mesh size is used. Therefore, the energy penalization method exerts a significance to the convergence rate.

Fig.10 Improvement in convergence rate by the energy penalization method for ITO–MMC.

Full size|PPT slide

Bilateral symmetric structure by the ITO–MMC with a symmetric constraint

The structure shown in Fig. 11 is a bilateral symmetric beam. The TO of the bilateral symmetric beam can be viewed as a TO problem subjected to a bilateral symmetric constraint. Hence, a beam having a bilateral symmetric boundary can be optimized by the process illustrated in the left of Fig. 8. The designable domain, unit cells, the dimensions of the design domain and the prescribed material volume ratio, as well as the boundary condition, are illustrated in Fig. 11. In Fig. 11, the green dashed line is the symmetric axis of the design domain, the beam is fixed at the right and left bottoms, and a vertical downward force is applied to the middle of the top boundary of the design domain.

Fig.11 Design domain of the bilateral symmetric beam example.

Full size|PPT slide

The beam is first optimized without a bilateral symmetric constraint to testify the feasibility of ITO–MMC in the optimal design of the beam illustrated in Fig. 11. The design domain is discretized into a quadratic NURBS mesh with the size of \mathrm{180}\times \mathrm{30}, and the initial topology of the structure consisting of 48 basic structural components is depicted in Fig. 12. The final optimized structure of the beam using ITO–MMC is shown in Fig. 13(a). Then, the optimal results of the beam using the traditional SIMP method and SIMP with a density filter are presented in Figs. 13(b) and 13(c), respectively. Apparently, the final result optimized by ITO–MMC has a similar optimal structure obtained by either SIMP or SIMP with a Heaviside density filter, which indicates the effectiveness of ITO–MMC.

Fig.12 Initial design of the bilateral symmetric beam.

Full size|PPT slide

Fig.13 Optimal results of the bilateral symmetric beam (a) obtained by the ITO–MMC without a bilateral constraint and the number of optimization cycle is ste{p}_{\mathrm{i}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}}=\mathrm{8}\mathrm{3}\mathrm{7}, C=57.54; (b) generated by SIMP with filter radius r_{\min }=\mathrm{4} based on the 88-line code using sensitivity filter, C=66.58; and (c) gained by SIMP with r_{\min }=\mathrm{4} based on the variant of 88-line code using density filter, C=58.32.

Full size|PPT slide

Afterward, the bilateral symmetric beam is optimized by the ITO–MMC with a bilateral symmetric constraint, and the optimization process is presented in the left of Fig. 8. Quadratic NURBS elements with the size of \mathrm{180}\times \mathrm{30} are adopted, which are the same as those used in the case of the ITO–MMC without a bilateral symmetric constraint. Then, the initial design of the first cell illustrated in Fig. 11 consists of 24 basic structural components, which is half of that used in the ITO–MMC without a bilateral symmetric constraint, and those components are as depicted in Fig. 14(a). The optimal design of the beam is obtained through 348 optimization iterations by using the model of the ITO–MMC subjected to a bilateral symmetric constraint presented in Eq. (20), and the final optimal layout is presented in Fig. 14(b). Finally, the convergent history of the optimization process using the ITO–MMC with a bilateral symmetric constraint is illustrated in Fig. 15. Basing from the convergence history and the corresponding results, we can conclude that the ITO–MMC with a bilateral symmetric constraint can handle the TO problem subjected to a bilateral symmetric constraint effectively. In addition, the number of iterations is reduced from 837 to 348, and the time consumption is decreased from 49.6 to 11.9 s for each iteration, when the ITO–MMC is replaced by the ITO–MMC under a bilateral symmetric constraint. This result indicates that the convergence rate is largely improved by the new ITO–MMC model considering the bilateral symmetric constraint due to the reduction of half of the design variables. The convergence rate of the ITO–MMC under a bilateral symmetric constraint is improved up to 58.4% and the computational efficiency is accelerated by 76% for each iteration compared with the ITO–MMC without a bilateral symmetric constraint. Hence, the speedup is up to 10, which is defined by the time of the ITO–MMC with a bilateral symmetric constraint versus to the time of the ITO–MMC.

Fig.14 (a) Initial design of components located in the first cell of which the geometric parameters are the design variables used in the ITO–MMC model formulated as Eq. (20); (b) the final result generated by Eq. (20).

Full size|PPT slide

Fig.15 Convergence histories of volume fraction, objective function, and topology for the beam depicted in Fig. 11 using the model of the ITO–MMC subjected to a bilateral symmetric constraint.

Full size|PPT slide

Periodic symmetric structure by the ITO–MMC with a symmetric constraint

The proposed method for optimizing the periodic structure as shown in the right of the flowchart illustrated in Fig. 8 can be used to generate the optimal design of a sandwich structure with a periodic symmetric constraint as depicted in Fig. 16 where skins are attached to the sandwich structure and the ends of skins are fixed. The design domain is a rectangular region with the size of \mathrm{16}\times \mathrm{4}, which is discretized into 160 × 40 quadratic NURBS elements, and the top and bottom layers have 0.1 thickness and are discretized into 160 quadratic NURBS beam elements. Only the mean compliance of the structure is considered, and the upper volume fraction is set at \overline{V}=\mathrm{0.4}. The Young’s modulus and Poisson’s ratio of the material of skins are assumed to be 100 and 0.3, respectively. The middle of the top skin is subjected to a vertical force F=\mathrm{1}.

Fig.16 Design domain of 2D sandwich structure.

Full size|PPT slide

The convergent histories of volume fraction, objective function, topology for m=\mathrm{2}\times \mathrm{1}, m=\mathrm{4}\times \mathrm{1}, and m=\mathrm{8}\times \mathrm{2} are respectively illustrated in Fig. 17. All topologies shown in Fig. 17 satisfy the prescribed periodic symmetric constraint. Figure 18 depicts the initial designs of components for m=\mathrm{2}\times \mathrm{1}, m=\mathrm{4}\times \mathrm{1}, and m=\mathrm{8}\times \mathrm{2}. The final topology and the corresponding mean compliance for the three periodic symmetric constraints are presented in Fig. 19, where the material layout inside dash lines is the topology of a representative unit cell for each specific case. Those topologies depicted in Fig. 19 are similar to those generated by BESO [⁶⁰]. In addition, the final optimal topologies are different for different numbers of unit cells, indicating that the optimal topology in ITO–MMC under a periodic symmetric constraint depends on the number of unit cells for decomposing the design domain. The objective function is increased with increasing cell number, considering that the design space is decreased with increasing cell number and the solution in the small space is inferior to that obtained in the large design space. Furthermore, compared with the BESO with a periodic symmetric constraint, the ITO–MMC with a periodic symmetric constraint can more effectively obtain the explicit geometric parameters from the final optimal results, which are vital to manufacturing the final optimized structure. In addition, the number of design variables is much less under the same mesh. In specific, the number of design variables in Fig. 19(c) is 70, whereas that in the BESO is up to 6400 when a finite element mesh with the size of 160 × 40 is used.

Fig.17 Curves of volume fraction, compliance, and topology for 2D sandwich structure with (a) m=\mathrm{2}\times \mathrm{1}, (b) m=\mathrm{4}\times \mathrm{1}, and (c) m=\mathrm{8}\times \mathrm{2} using the model of ITO–MMC under a periodic symmetric constraint.

Full size|PPT slide

Fig.18 Initial design of components: (a) m=\mathrm{2}\times \mathrm{1}; (b) m=\mathrm{4}\times \mathrm{1}; and (c) m=\mathrm{8}\times \mathrm{2}.

Full size|PPT slide

Fig.19 Optimized results for 2D sandwich structures under three periodic constraints: (a) m=\mathrm{2}\times \mathrm{1}, C=10.54; (b) m=\mathrm{4}\times \mathrm{1}, C=12.94; and (c) m=\mathrm{8}\times \mathrm{2}, C=15.09.

Full size|PPT slide

Optimization of symmetric structure by the energy penalization method

Two benchmarks considering the bilateral symmetric constraint and the periodic symmetric constraint are considered to testify the validity of the energy penalization method.

Bilateral symmetric structure

A 2D beam subjected to a lateral load is depicted in Fig. 20(a), where the dimension of the design domain is set to \mathrm{2}\times \mathrm{1}. The corner of the left and right bottom of the asymmetric loaded beam is fixed, the asymmetric load F_{\mathrm{1}} is applied at the upper right corner, the specified volume ratio is 0.5, and a bilateral symmetric constraint is considered. The beam is divided into 80 × 40 quadratic NURBS elements for analysis. In addition, the max number of iterations is set to 1000. Figure 20(b) depicts the initial design of components. The results are listed in Fig. 21 obtained by the energy penalization method with p taking different values used in Eq. (26). Equation (26) is equivalent to Eq. (22), and the penalization effect is disappeared when p=\mathrm{1}.

Fig.20 Problem setting and initial design of components for beam subjected to symmetric and volume constraints. (a) Problem setting; (b) initial design of components.

Full size|PPT slide

According to the final optimized structure presented in Fig. 21, the objective values represented by C are reduced by 5.6% and the optimized structures have a clear boundary in the final topology when p takes 0.6, 0.7, and 0.8 than when p takes 1 equally to the ITO–MMC considering a symmetric constraint without energy penalization. The overall iteration number for p=0.6 is 491, which is less than the iterations for p=0.7, p=0.8, and p=1. Moreover, from both the objective value and the clearness of boundary of the final topology standpoint, the results cannot be improved for p=0.4, p=0.5, and p=0.9 than p=1. Hence, the optimal penalization parameter used in Eq. (26) is p=0.6, and the energy penalization method is valid for the ITO–MMC optimizing the symmetric structure because the strain energy of the cell with a smaller strain energy is magnified when p=0.6, p=0.7, and p=0.8, which enlarge the influence of the element having a small strain energy on the design variables. Obviously, overly small penalization parameter causes excessive punishment, which is indicated by the results obtained when p=0.4 and p=0.5.

Fig.21 Final optimized results for p taking different values in Eq. (26).

Full size|PPT slide

Periodic symmetric structure

A 2D beam subjected to multiple loads (Fig. 22) is used as the benchmark to verify the validity of energy penalization in the optimization of the periodic symmetric structure. In addition, the size of the design domain is set to \mathrm{2}\times \mathrm{1}. The left edge of the beam is fixed, the prescribed external load F_{\mathrm{1}} and F_{\mathrm{2}} are respectively applied at the upper and bottom right corners, and the upper-bounded volume is specified to 0.5, and a periodic symmetric constraint is considered. As shown in Section 6.3, the optimal results depend on the cell number. The design domain is divided into m=\mathrm{2}\times \mathrm{1}, m=\mathrm{4}\times \mathrm{1}, and m=\mathrm{4}\times \mathrm{2} unit cells to validate the versatility of the proposed energy penalization method for the optimization of a periodic structure. Meshes with the sizes of 80×40, 120×60, and 120×60 are used for the three cases where m=\mathrm{2}\times \mathrm{1}, m=\mathrm{4}\times \mathrm{1}, and m=\mathrm{4}\times \mathrm{2} unit cells are the specific periodic symmetric constraints. The maximum iteration of the optimization process is set to 2000. The initial designs of components for those three cases are shown in Fig. 23. Similar to the bilateral symmetric case, the penalization effect is disappeared when p=\mathrm{1} is used in Eq. (27).

Fig.22 Problem setting for multiple load beam considering periodic symmetric constraint.

Full size|PPT slide

Fig.23 Three initial designs of components for different numbers of unit cells as the periodic symmetric constraint of the multiple load beam.

Full size|PPT slide

The optimized structures and the corresponding objective values, as well as the iterations for the three specific periodic symmetric constraints, are illustrated in Figs. 24–26. For the m=\mathrm{2}\times \mathrm{1} periodic symmetric constraint, the iterations for all p<1 cases are less than that for the p=1 case. The objective values for the p=0.1, p=0.2, and p=0.3 cases are considerably larger than those for the other cases, and the objective values are decreased with the increase in p for p=0.1, p=0.2, and p=0.3, as shown in Fig. 24. For p=0.1, p=0.2, and p=0.3, the objective values are decreased with the increase in p because the excessive penalization effect of the strain energy of elements in the cells close to the load is alleviated with the increase in p. Then, the excessive penalization effect disappears when p is no less than 0.4 (Fig. 24). Furthermore, the optimal factor used in Eq. (27) is p=0.5 in terms of convergence and objective value as well as the boundary clearness of topology. In addition, the iteration number of p=0.5 is only 235, which is lower by 51% than that of the optimization process without energy penalization, whereas the objective values for p=0.5 are slightly less than that when p takes 1 for a m=\mathrm{2}\times \mathrm{1} periodic symmetric constraint. The comparison in convergence history between p=0.5 and p=1 is shown in Fig. 27. Hence, the energy penalization method accelerates the optimization of the ITO–MMC considering a m=\mathrm{2}\times \mathrm{1} periodic symmetric constraint.

Fig.24 Final results are presented, for the beam divided into m=\mathrm{2}\times \mathrm{1} unit cells optimized by the energy penalization method, where different penalization values are used in Eq. (27).

Full size|PPT slide

As shown in Fig. 25 for a m=\mathrm{4}\times \mathrm{1} periodic symmetric constraint applied to the multiple load beam, the iterations for p=0.4, p=0.6, p=0.686, and p=0.7 are less than the iterations for the p=1 case, whereas the optimization becomes suboptimal when p takes 0.4 due to the excessive penalization effect of the strain energy of elements close to the load. Moreover, p=0.686 is the optimal coefficient used in the energy penalization formula from the viewpoint of objective function and convergence as well as the boundary clearness of topology. The iteration number of p=0.686 is only 880, which is lower by 56% than that when p=1, whereas the objective value for p=0.686 is slightly larger than that obtained from the optimization process without energy penalization for a m=\mathrm{4}\times \mathrm{1} periodic symmetric constraint. Thus, similar to the m=\mathrm{2}\times \mathrm{1} periodic symmetric constraint case, the energy penalization method is highly effective for the ITO–MMC taking a m=\mathrm{4}\times \mathrm{1} periodic symmetric constraint into consideration.

Fig.25 Final results obtained by the energy penalization method where different penalization values used in Eq. (27) are presented for the beam subjected to a m=\mathrm{4}\times \mathrm{1} periodic symmetric constraint.

Full size|PPT slide

As demonstrated in Fig. 26, the iterations for p=0.83, p=0.85, p=0.9, p=0.92, and p=0.95 in a m=\mathrm{4}\times \mathrm{2} periodic symmetric constraint are less than the iterations for the p=1 case, and the differences in objective value among the p=0.83, p=0.85, p=0.9, p=0.92, and p=0.95 cases are negligible. According to the principle of the least iterations, p=0.92 outperforms the other cases listed in Fig. 26, and the iteration number of p=0.92 is only 100, which is 88% lower than that when p=1. Meanwhile, the stiffness of the optimized structure for p=0.92 is slightly higher than that obtained in the optimization process without energy penalization. Thus, similar to the m=\mathrm{2}\times \mathrm{1} and m=\mathrm{4}\times \mathrm{1} periodic symmetric constraint cases, the validity of the energy penalization method for the ITO–MMC considering a m=\mathrm{4}\times \mathrm{2} periodic symmetric constraint is verified.

Fig.26 Optimized structures obtained by the energy penalization method where the coefficient used in Eq. (27) takes different values for a m=\mathrm{4}\times \mathrm{2} periodic symmetric constraint.

Full size|PPT slide

Fig.27 Convergence histories for p=0.5 and p=1 under a m=\mathrm{2}\times \mathrm{1} periodic symmetric constraint.

Full size|PPT slide

Therefore, the proposed energy penalization method is very valid for enhancing the convergence and the numerical stability of the ITO–MMC under a periodic constraint. Moreover, the optimal factor used in Eq. (27) increases with cell number, considering that the optimal factors for m=\mathrm{2}\times \mathrm{1}, m=\mathrm{4}\times \mathrm{1}, and m=\mathrm{4}\times \mathrm{2} unit cells are 0.5, 0.686, and 0.92, respectively. The reason for the aforementioned conclusion is that the enlargement effect of small-strain-energy elements increases with the increase in the number of unit cells for the same factor used in Eq. (27).

Conclusions

An energy penalization method is proposed for ITO–MMC to solve the low convergence problem of asymmetric structure by ITO–MMC. The improvement in convergence rate for all cases listed in Table 1 ranges from 18.6% to 44.5% with the energy penalization method within the optimization framework of ITO–MMC. Hence, the ITO–MMC based on the energy penalization method is more effective for the optimization of asymmetric structure than ITO–MMC. Then, to solve the TO problems of symmetric structures, including bilateral and periodic structures, we develop two formulae of ITO–MMC and extend the energy penalization method to solve the optimization of structures subjected to asymmetric loads. An additional bilateral or periodic constraint has been considered in ITO–MMC, which ensures the symmetry of the structure. The numerical stability and convergence rate of the optimization process of asymmetric loaded structure is largely improved by the energy penalization method for the ITO–MMC with a symmetric constraint.

The optimal material distribution within a unit cell is obtained from gradually optimizing the basic structural components located in the unit cell using the ITO–MMC method. Four benchmarks are studied in this work to demonstrate the validity of the new ITO–MMC formulae and the energy penalization method. Numerical results show that the improvement from the view of the iteration number and time consumption for each iteration are up to 58.4% and 76%, respectively, compared with the ITO–MMC without a bilateral symmetric constraint. The optimal topology in the ITO–MMC under a periodic symmetric constraint depends on the number of unit cells. In addition, the speedup of the overall optimization process is up to 10 when ITO–MMC considers the bilateral symmetric constraint for the bilateral symmetric structures.

With the aid of the energy penalization method, the numerical stability for the ITO–MMC under a bilateral symmetric constraint is largely improved, which can be observed from the final optimized topology that has a clear boundary for p=0.6, whereas the boundary of the final topology is not clear for p=1, and the improvement in objective value is up to 5.6% with less iterations when p takes 0.6 compared with the case without the energy penalization effect used in the optimization. Meanwhile, the proposed energy penalization method can largely speedup the convergence rate and improve the numerical stability of the ITO–MMC under a periodic symmetric constraint, and the speedups are more than 2 (defined by the time of ITO–MMC with energy penalization effect versus to the time of ITO–MMC without energy penalization effect).

An ITO–MMC subjected to a symmetric constraint with an automatic reconstruction of boundary condition [⁶²,⁶³] is under construction. In addition, the computational efficiency of the proposed optimization framework can be largely improved by a GPU algorithm [⁶⁴], and it needs further research. Furthermore, the proposed method will be extended to the problem of frequency–stiffness optimization.

Acknowledgements

This research was supported by the National Natural Science Foundation of China (Grant Nos. 51675197 and 51705158), the National Engineering Research Center of Near-Net-Shape Forming for Metallic Materials, Ministry of Education Key Laboratory of High Efficient Near-Net-Shape Forming Technology and Equipment for Metallic Materials (Category B) Opening Foundation (Grant No. 2018005). The support is gratefully acknowledged.ƒ

Open Access

This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution, and reproduction in any medium or format, as long as credit is given to the original author(s) and the source, a link to the Creative Commons license is provided, and any changes made are indicated.ƒ

The images or other third-party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.ƒ

To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

References

1 Bendsøe M P, Kikuchi N. Generating optimal topologies in structural design using a homogenization method. Computer Methods in Applied Mechanics and Engineering, 1988, 71(2): 197–224

2 Gholizadeh S, Ebadijalal M. Performance based discrete topology optimization of steel braced frames by a new metaheuristic. Advances in Engineering Software, 2018, 123: 77–92

3 Zhang Y, Ge W, Zhang Y, et al. Topology optimization of hyperelastic structure based on a directly coupled finite element and element-free Galerkin method. Advances in Engineering Software, 2018, 123: 25–37

4 Csébfalvi A, Lógó J. A critical analysis of expected-compliance model in volume-constrained robust topology optimization with normally distributed loading directions, using a minimax-compliance approach alternatively. Advances in Engineering Software, 2018, 120: 107–115

5 Gao J, Li H, Gao L, et al. Topological shape optimization of 3D micro-structured materials using energy-based homogenization method. Advances in Engineering Software, 2018, 116: 89–102

6 Sigmund O, Petersson J. Numerical instabilities in topology optimization: A survey on procedures dealing with checkerboards, mesh-dependencies and local minima. Structural Optimization, 1998, 16(1): 68–75

7 Sigmund O. A 99 line topology optimization code written in Matlab. Structural and Multidisciplinary Optimization, 2001, 21(2): 120–127

8 Andreassen E, Clausen A, Schevenels M, et al. Efficient topology optimization in MATLAB using 88 lines of code. Structural and Multidisciplinary Optimization, 2011, 43(1): 1–16

9 Liao Z, Zhang Y, Wang Y, et al. A triple acceleration method for topology optimization. Structural and Multidisciplinary Optimization, 2019, 60(2): 727–724

10 Zhou M, Liu Y, Lin Z. Topology optimization of thermal conductive support structures for laser additive manufacturing. Computer Methods in Applied Mechanics and Engineering, 2019, 353: 24–43

11 Zhou M, Lian H, Sigmund O, et al. Shape morphing and topology optimization of fluid channels by explicit boundary tracking. International Journal for Numerical Methods in Fluids, 2018, 88(6): 296–313

12 Liu T, Wang S, Li B, et al. A level-set-based topology and shape optimization method for continuum structure under geometric constraints. Structural and Multidisciplinary Optimization, 2014, 50(2): 253–273

13 Liu T, Li B, Wang S, et al. Eigenvalue topology optimization of structures using a parameterized level set method. Structural and Multidisciplinary Optimization, 2014, 50(4): 573–591

14 Mei Y, Wang X. A level set method for structural topology optimization and its applications. Computer Methods in Applied Mechanics and Engineering, 2004, 35(7): 415–441

15 Wang Y, Benson D J. Geometrically constrained isogeometric parameterized level-set based topology optimization via trimmed elements. Frontiers of Mechanical Engineering, 2016, 11(4): 328–343

16 Wang Y, Benson D J. Isogeometric analysis for parameterized LSM-based structural topology optimization. Computational Mechanics, 2016, 57(1): 19–35

17 Xia Z, Wang Y, Wang Q, et al. GPU parallel strategy for parameterized LSM-based topology optimization using isogeometric analysis. Structural and Multidisciplinary Optimization, 2017, 56(2): 413–434

18 Liu Y, Li Z, Wei P, et al. Parameterized level-set based topology optimization method considering symmetry and pattern repetition constraints. Computer Methods in Applied Mechanics and Engineering, 2018, 340: 1079–1101

19 Li Z, Shi T, Xia Q. Eliminate localized eigenmodes in level set based topology optimization for the maximization of the first eigenfrequency of vibration. Advances in Engineering Software, 2017, 107: 59–70

20 Liu J, Li L, Ma Y. Uniform thickness control without pre-specifying the length scale target under the level set topology optimization framework. Advances in Engineering Software, 2017, 115: 204–216

21 Xia Q, Shi T, Xia L. Stable hole nucleation in level set based topology optimization by using the material removal scheme of BESO. Computer Methods in Applied Mechanics and Engineering, 2019, 343: 438–452

22 Li Z, Shi T, Xia L, et al. Maximizing the first eigenfrequency of structures subjected to uniform boundary erosion through the level set method. Engineering with Computers, 2019, 35(1): 21–33

23 Xie Y M, Steven G P. A simple evolutionary procedure for structural optimization. Computers & Structures, 1993, 49(5): 885–896

24 Huang X, Xie Y M. Convergent and mesh-independent solutions for the bi-directional evolutionary structural optimization method. Finite Elements in Analysis and Design, 2007, 43(14): 1039–1049

25 Stanford B, Ifju P. Aeroelastic topology optimization of membrane structures for micro air vehicles. Structural and Multidisciplinary Optimization, 2009, 38(3): 301–316

26 Zhu J H, Zhang W H, Xia L. Topology optimization in aircraft and aerospace structures design. Archives of Computational Methods in Engineering, 2016, 23(4): 595–622

27 Zegard T, Paulino G H. Bridging topology optimization and additive manufacturing. Structural and Multidisciplinary Optimization, 2016, 53(1): 175–192

28 Guo X, Zhou J, Zhang W, et al. Self-supporting structure design in additive manufacturing through explicit topology optimization. Computer Methods in Applied Mechanics and Engineering, 2017, 323: 27–63

29 Li Y F, Huang X, Meng F, et al. Evolutionary topological design for phononic band gap crystals. Structural and Multidisciplinary Optimization, 2016, 54(3): 595–617

30 Li Y F, Huang X, Zhou S. Topological design of cellular phononic band gap crystals. Materials, 2016, 9(3): 186–208

31 Wang Y, Arabnejad S, Tanzer M, et al. Hip implant design with three-dimensional porous architecture of optimized graded density. Journal of Mechanical Design, 2018, 140(11): 111406

32 Aage N, Andreassen E, Lazarov B S, et al. Giga-voxel computational morphogenesis for structural design. Nature, 2017, 550(7674): 84–86

33 Sigmund O. Materials with prescribed constitutive parameters: An inverse homogenization problem. International Journal of Solids and Structures, 1994, 31(17): 2313–2329

34 Zhang W, Sun S. Scale-related topology optimization of cellular materials and structures. International Journal for Numerical Methods in Engineering, 2006, 68(9): 993–1011

35 Vicente W M, Picelli R, Pavanello R, Xie Y M. Topology optimization of periodic structures for coupled acoustic-structure systems. In: Proceedings of the VII European Congress on Computational Methods in Applied Sciences and Engineering. Crete Island, 2016

36 He G, Huang X, Wang H, et al. Topology optimization of periodic structures using BESO based on unstructured design points. Structural and Multidisciplinary Optimization, 2016, 53(2): 271–275

37 Guo X, Zhang W, Zhong W. Doing topology optimization explicitly and geometrically—A new moving morphable components based framework. Frontiers in Applied Mechanics, 2014, 81(8): 081009

38 Hou W, Gai Y, Zhu X, et al. Explicit isogeometric topology optimization using moving morphable components. Computer Methods in Applied Mechanics and Engineering, 2017, 326: 694–712

39 Hughes T J R, Cottrell J A, Bazilevs Y. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Computer Methods in Applied Mechanics and Engineering, 2005, 194(39–41): 4135–4195

40 Cottrell J A, Hughes T J R, Bazilevs Y. Isogeometric Analysis: Toward Integration of CAD and FEA. Hoboken: John Wiley & Sons, 2009

41 Wang X, Zhu X, Hu P. Isogeometric finite element method for buckling analysis of generally laminated composite beams with different boundary conditions. International Journal of Mechanical Sciences, 2015, 104: 190–199

42 Marussig B, Hughes T J R. A review of trimming in isogeometric analysis: Challenges, data exchange and simulation aspects. Archives of Computational Methods in Engineering, 2018, 25(4): 1059–1127

43 Pan Q, Rabczuk T, Chen C, et al. Isogeometric analysis of minimal surfaces on the basis of extended Catmull–Clark subdivision. Computer Methods in Applied Mechanics and Engineering, 2018, 337: 128–149

44 An Z, Yu T, Bui T Q, et al. Implementation of isogeometric boundary element method for 2-D steady heat transfer analysis. Advances in Engineering Software, 2018, 116: 36–49

45 Rypl D, Patzák B. Construction of weighted dual graphs of NURBS-based isogeometric meshes. Advances in Engineering Software, 2013, 60–61: 31–41

46 Rypl D, Patzák B. Object oriented implementation of the T-spline based isogeometric analysis. Advances in Engineering Software, 2012, 50(1): 137–149

47 Lai W, Yu T, Bui T Q, et al. 3-D elasto-plastic large deformations: IGA simulation by Bézier extraction of NURBS. Advances in Engineering Software, 2017, 108: 68–82

48 Rypl D, Patzák B. From the finite element analysis to the isogeometric analysis in an object oriented computing environment. Advances in Engineering Software, 2012, 44(1): 116–125

49 Dedè L, Borden M J, Hughes T J R. Isogeometric analysis for topology optimization with a phase field model. Archives of Computational Methods in Engineering, 2012, 19(3): 427–465

50 Kumar A V, Parthasarathy A. Topology optimization using B-spline finite elements. Structural and Multidisciplinary Optimization, 2011, 44(4): 471–481

51 Wang Y, Xu H, Pasini D. Multiscale isogeometric topology optimization for lattice materials. Computer Methods in Applied Mechanics and Engineering, 2017, 316: 568–585

52 Xu M, Xia L, Wang S, et al. An isogeometric approach to topology optimization of spatially graded hierarchical structures. Composite Structures, 2019, 225: 111171

53 Xie X, Wang S, Xu M, et al. A new isogeometric topology optimization using moving morphable components based on R-functions and collocation schemes. Computer Methods in Applied Mechanics and Engineering, 2018, 339: 61–90

54 Zhang W, Zhou J, Zhu Y, et al. Structural complexity control in topology optimization via moving morphable component (MMC) approach. Structural and Multidisciplinary Optimization, 2017, 56(3): 535–552

55 Piegl L, Tiller W. The NURBS Book. Berlin: Springer, 1997

56 de Falco C, Reali A, Vázquez R. GeoPDEs: A research tool for isogeometric analysis of PDEs. Advances in Engineering Software, 2011, 42(12): 1020–1034

57 de Boor C D. On calculating with B-splines. Journal of Approximation Theory, 1972, 6(1): 50–62

58 Fougerolle Y D, Gribok A, Foufou S, et al. Boolean operations with implicit and parametric representation of primitives using R-functions. IEEE Transactions on Visualization and Computer Graphics, 2005, 11(5): 529–539

59 Svanberg K. The method of moving asymptotes—A new method for structural optimization. International Journal for Numerical Methods in Engineering, 1987, 24(2): 359–373

60 Huang X, Xie Y M. Optimal design of periodic structures using evolutionary topology optimization. Structural and Multidisciplinary Optimization, 2008, 36(6): 597–606

61 Zhang W, Yuan J, Zhang J, et al. A new topology optimization approach based on moving morphable components (MMC) and the ersatz material model. Structural and Multidisciplinary Optimization, 2016, 53(6): 1243–1260

62 Xia Z, Wang Q, Liu Q, et al. A novel approach for automatic reconstruction of boundary condition in structure analysis. Advances in Engineering Software, 2016, 96: 38–57

63 Xia Z, Wang Q, Wang Y, et al. A CAD/CAE incorporate software framework using a unified representation architecture. Advances in Engineering Software, 2015, 87(C): 68–85

64 Wang Y, Wang Q, Deng X, et al. Graphics processing unit (GPU) accelerated fast multipole BEM with level-skip M2L for 3D elasticity problems. Advances in Engineering Software, 2015, 82(2): 105–118