Concurrent optimization of structural topology and infill properties with a CBF-based level set method

Abstract

In this paper, a parametric level-set-based topology optimization framework is proposed to concurrently optimize the structural topology at the macroscale and the effective infill properties at the micro/meso scale. The concurrent optimization is achieved by a computational framework combining a new parametric level set approach with mathematical programming. Within the proposed framework, both the structural boundary evolution and the effective infill property optimization can be driven by mathematical programming, which is more advantageous compared with the conventional partial differential equation-driven level set approach. Moreover, the proposed approach will be more efficient in handling nonlinear problems with multiple constraints. Instead of using radial basis functions (RBF), in this paper, we propose to construct a new type of cardinal basis functions (CBF) for the level set function parameterization. The proposed CBF parameterization ensures an explicit impose of the lower and upper bounds of the design variables. This overcomes the intrinsic disadvantage of the conventional RBF-based parametric level set method, where the lower and upper bounds of the design variables oftentimes have to be set by trial and error. A variational distance regularization method is utilized in this research to regularize the level set function to be a desired distance-regularized shape. With the distance information embedded in the level set model, the wrapping boundary layer and the interior infill region can be naturally defined. The isotropic infill achieved via the mesoscale topology optimization is conformally fit into the wrapping boundary layer using the shape-preserving conformal mapping method, which leads to a hierarchical physical structure with optimized overall topology and effective infill properties. The proposed method is expected to provide a timely solution to the increasing demand for multiscale and multifunctional structure design.

Keywords

concurrent topology optimization / parametric level set method / cardinal basis function / shell-infill structure design / conformal mapping

1 Introduction

Structures containing cellular infills or micro/meso architectures can possess fine-tuned properties and extra functionalities with a low density [1–4]. With the help of modern additive manufacturing technologies [5–7], the multiscale structures are foreseeing a wide range of engineering applications [8–10]. However, how to identify the exact structural layout at different scales to maximize the performance is a challenging topic.

At micro/meso scale, topology optimization of metamaterials or lattice structures has been well studied since Sigmund’s pioneering work [11], where the material is optimally distributed inside the design domain to achieve desired effective material properties. The homogenization-based topology optimization method [12–14], the solid isotropic material with penalization (SIMP) method [15,16], the level set method (LSM) [17–19], and the evolutionary method [20] are most commonly used approaches in this field. With an inverse homogenization method, Sigmund designed truss-like base cells with different material properties [21,22]. With the SIMP method, Sigmund [23] elaborated the three-phased microstructure designs with extreme thermoelastic properties. As for the conventional LSM, Mei and Wang [24] designed multimaterial microstructures utilizing the “color” level sets. Vogiatzis et al. [25] employed a reconciled level set approach to design multimaterial microstructures with negative Poisson’s ratio (NPR). As an extension to the conventional level set approach, Wang et al. [26] employed a numerically robust parametric level set method (PLSM), in order to design a series of microstructures with different Poisson’s ratios. In this systematic computational design framework, the homogenization was used for calculating the material property and PLSM was used for updating the structural shape and topology. With the evolutionary approach, a series of papers have been published for designing microstructures with extreme properties [27–29]. A comprehensive review for the metamaterial design using evolutionary structural optimization is provided in Ref. [30].

Instead of merely designing the micro/meso structures, recent years have witnessed more and more researches on optimizing both the structure and its constructing material properties in a concurrent manner [31]. The SIMP-based concurrent structural design was reported by Rodrigues et al. [32] and Coelho et al. [33], for designing the overall structural geometry and its infill metamaterials. Liu et al. [34] uncoupled the macro and the micro densities in the SIMP method, to carry out a concurrent topology optimization of truss-like material. Deng and Chen [35] studied the concurrent topology optimization by designing structures with multiple porous metamaterials under radon field loading uncertainty. To achieve a high-performance shell-infill structure, Wu et al. [36,37] proposed a concurrent SIMP-based topology optimization with a local volume constraint to design bone-like structures. The SIMP-based concurrent approach has also been extended to designs utilizing thermal elastic materials [38]. Within the level set framework, Sivapuram et al. [39] concurrently optimized the overall topology and the metamaterial properties with a decomposition formulation for an easy algorithm parallelization. Wang et al. [40] reported a concurrent multiscale design with a shape metamorphosis technology, aiming at ensuring the connectivity of the adjacent microstructures. Li et al. [41] proposed an integrated multiscale structure design, where SIMP is used for macroscale structure design and a radial basis functions (RBF)-based PLSM is used for microscale structure design. Li et al. [42] further developed the multiscale structure design based on the concurrent PLSM approach, in which the connectivity between the adjacent microstructure patches were ensured by virtual kinematic connectors. Via the evolutionary approach, Xia et al. [10] designed the overall structural topology and local material property concurrently, using the multilevel finite element approach (FE²) to handle the material nonlinearity. Readers interested in evolutionary approaches to concurrent designs of multiscale structures can be referred to Refs. [43–46] for further details. In this paper we employ the level set approach and improve it to make it suitable for multiscale topology optimization.

With an implicit design representation, the designs represented by the level set model has clear design boundaries. Mathematically, the level set model can possess extra geometrical information such as curvatures, and can handle topological changes in a robust way. In the conventional level-set-based topology optimization, the Hamilton-Jacobi partial differential equation (PDE) describes the dynamics of the boundary motion [47]. The design velocity field derived from shape sensitivity analysis [19,48,49] drives the design boundary to its optimum location. However, the conventional LSM has several numerical limitations [50–52]. Besides, introducing multiple constraints to the design can be challenging when the Lagrange multiplier method is used. When utilized for multiscale structure design, introducing extra design variables can be difficult for the conventional level set approach. To increase the numerical efficiency and robustness of conventional LSM, Wang et al. [53] proposed PLSM with the RBF [54]. By utilizing the method of moving asymptotes (MMA) [55] as the optimizer, numerical efficiency can be improved and design constraints can be imposed in a straightforward way. The PLSM provides the flexibility in parameterizing the structural topology as well as the material properties, which enables the concurrent topology optimization for multiscale structures. Apart from using RBF, the R-function and B-spline are also used in the level set parameterization scheme. The readers can be referred to Refs. [56–58] for further details.

The RBF-based PLSM approximates the level set function as the weighted summation of RBF kernel functions. As design variables, those kernel function weights are updated during the optimization. However, those kernel function weights have no exact physical meaning, and their lower and upper bounds cannot be explicitly defined. In practice, those bounds are set by trial and error. To get the explicit bounds for the design variables, a cardinal basis function (CBF) was constructed based on the RBF partition of unity collocation method [59]. The CBF satisfies the Cardinal properties which equals to 1 at the selected node and equals to 0 elsewhere. This means when CBF is used as the kernel function, the corresponding weights have the physical meaning of being the heights of the level set function. Moreover, a distance regularization energy functional [60,61] can be introduced to regularize the level set function to be a distance function in a zone along the boundary and flat surfaces elsewhere. The regularized distance function can ensure an accurate material property interpolation from the level set model to the physics model. The flat surfaces can impose preferred numerical stability and can help the creation of new holes [62]. Besides, the signed distance function can preserve the distance information directly. This means with one level set function, a shell-infill multiscale structure can be designed [63–65]. The shell layer of this structure can function as a protective coating against loading or extreme exterior conditions [66] while the infills can provide the designed functionalities.

Another important issue of multiscale structure optimization is to consider its manufacturability. The proposed multiscale structure is generated by mapping the optimal infill into its wrapping shell. This process will be tricky on the boundary areas, since the infill cells will be rotated or deformed. A rotated or deformed anisotropic infill cell cannot preserve its designed property. To address this issue, an isotropic infill micro/meso structure is designed to be filled into the interior infill region. The isotropy of the infill structure is ensured by imposing the isotropy constraint [67–70] to the topology optimization process. As for the mapping, a local shape-preserving conformal mapping is employed in this research. With the conformal mapping, the material properties of the designed infill structure can be mathematically preserved. The final mapped structure can be sent to 3D printer directly to be manufactured with a high fidelity of its designed material property.

To sum up, in this research, the concurrent multiscale structure design is carried out via a variational CBF-based PLSM. The infill properties are treated as design variables to be updated together with the structural design simultaneously. A unified isotropic infill metamaterial is designed to lower the computational cost for the prefabrication of 3D printing. With a modified Heaviside material property interpolation, from macroscale, the void, the wrapping shell, and the infill material region will be naturally discriminated. By constructing a user-defined distance regularization energy functional, the distance regularization effect can be directly controlled. A second stage optimization is directly carried out to find out the isotropic infill metamaterial layout. With both the overall geometry and the infill structure design in hand, the local shape-preserving conformal mapping is performed to integrate them together. To make the mapping adaptive to complicated geometries, instead of using the conventional Ricci flow with four control points [71], in this paper, a multi-control-point conformal mapping is proposed so that it is much flexible in controlling area distortions, especially for irregular domains.

The remaining paper is organized as follows. In Section 2, the strain energy method is introduced to calculate the effective material properties of a unit cell. In Section 3, the level set representation, its parameterization, and the distance regularization are elaborated. The general problem setting for the proposed concurrent optimization is formulated in Section 4. The detailed concurrent optimization for given problems are detailed in Section 5 with numerical results. The conformal mapping method is introduced in Section 6, together with mapped multiscale structure results. At last, this work is summarized in Section 7.

2 Prediction of the effective material properties based on the strain energy method

For a given unit cell, noting its homogenized elasticity tensor as $C_{ijkl}^H$, a constitutive law can be formulated based on the Hooke’s law:

(1)${\overline{\sigma }}_{ij}=C_{ijkl}^H{\overline{\epsilon }}_{kl},$

where ${\bar{\sigma }}_{ij}$ and ${\bar{\epsilon }}_{kl}$ are the effective stress and strain tensor, respectively. Under a linear elastic problem setting, the stress and the strain tensors of the homogeneous medium should be equal to the average stress and strain in the microstructure, as ${\bar{\sigma }}_{ij}={\displaystyle \frac{\mathrm{1}}{V}}{\displaystyle \int {\sigma }_{ij}dV}$ and ${\bar{\epsilon }}_{kl}={\displaystyle \frac{\mathrm{1}}{V}}{\displaystyle \int {\epsilon }_{kl}dV}$, where V is the volume of the unit cell.

With a 2D plane-stress assumption, Eq. (1) can be expanded to the form of Eq. (2):

(2)$\left[\begin{array}{l}{\overline{\sigma }}_{\mathrm{11}}\\ {\overline{\sigma }}_{\mathrm{22}}\\ {\overline{\sigma }}_{\mathrm{12}}\end{array}\right]=\left[\begin{array}{lll}{C}_{\mathrm{1111}}^H& C_{\mathrm{1122}}^H& \mathrm{0}\\ C_{\mathrm{2211}}^H& C_{\mathrm{2222}}^H& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& C_{\mathrm{1212}}^H\end{array}\right]\left[\begin{array}{l}{\overline{\epsilon }}_{\mathrm{11}}\\ {\overline{\epsilon }}_{\mathrm{22}}\\ \mathrm{2}{\overline{\epsilon }}_{\mathrm{12}}\end{array}\right]。$

On the other hand, another equilibrium equation can be built based on the fact that the strain energy stored inside both the homogeneous medium, noted as U^H, and inside the unit cell, noted as U, is the same. With four special strain field as shown in Fig. 1:

(3)${\overline{\epsilon }}_{\_\mathrm{1}}=\left\{\begin{array}{l}\mathrm{1}\\ \mathrm{0}\\ \mathrm{0}\end{array}\right\},{\overline{\epsilon }}_{\_\mathrm{2}}=\left\{\begin{array}{l}\mathrm{0}\\ \mathrm{1}\\ \mathrm{0}\end{array}\right\},{\overline{\epsilon }}_{\_\mathrm{3}}=\left\{\begin{array}{l}\mathrm{0}\\ \mathrm{0}\\ \mathrm{1}\end{array}\right\},{\overline{\epsilon }}_{\_\mathrm{4}}=\left\{\begin{array}{l}\mathrm{1}\\ \mathrm{1}\\ \mathrm{0}\end{array}\right\}\mathrm{.}$

Four strain energies can be calculated accordingly as U₁, U₂, U₃ and U₄. For example, with the first loading case, as listed in Eq. (3):

(4)$U_{\mathrm{1111}}^H={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}V{\left[\begin{array}{l}\mathrm{1}\\ \mathrm{0}\\ \mathrm{0}\end{array}\right]}^T\left[\begin{array}{lll}{C}_{\mathrm{1111}}^H& C_{\mathrm{1122}}^H& \mathrm{0}\\ C_{\mathrm{2211}}^H& C_{\mathrm{2222}}^H& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& C_{\mathrm{1212}}^H\end{array}\right]\left[\begin{array}{l}\mathrm{1}\\ \mathrm{0}\\ \mathrm{0}\end{array}\right]={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}V{C}_{\mathrm{1111}}^H=U_{\mathrm{1}}\mathrm{.}$

For a unit cell with volume of 1, $C_{\mathrm{1111}}^H=\mathrm{2}{U}_{\mathrm{1}}$. Similarly, based on different boundary conditions, the elastic tensor can be re-written as

(5)$\left[\begin{array}{lll}{C}_{\mathrm{1111}}^H& C_{\mathrm{1122}}^H& \mathrm{0}\\ C_{\mathrm{2211}}^H& C_{\mathrm{2222}}^H& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& C_{\mathrm{1212}}^H\end{array}\right]=\left[\begin{array}{lll}\mathrm{2}{U}_{\mathrm{1}}& U_{\mathrm{4}}-U_{\mathrm{2}}-U_{\mathrm{1}}& \mathrm{0}\\ U_{\mathrm{4}}-U_{\mathrm{2}}-U_{\mathrm{1}}& \mathrm{2}{U}_{\mathrm{2}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{2}{U}_{\mathrm{3}}\end{array}\right]\mathrm{.}$

Fig.1 The different loading cases for a unit cell. (a) Load case 1; (b) Load case 2; (c) Load case 3; (d) Load case 4

Full size|PPT slide

3 Level set representation, parameterization and the distance-regularized level set function

3.1 Level set implicit boundary representation

The level-set based implicit boundary representation [72] can be expressed as

(6)$\{\begin{array}{l}\Phi (x,t)>\mathrm{0}\begin{array}{ll}& \end{array}(x\in \Omega )\\ \Phi (x,t)=\mathrm{0}\begin{array}{ll}& \end{array}(x\in \Gamma )\\ \Phi (x,t)<\mathrm{0}\begin{array}{ll}& \end{array}(x\in D\backslash \Omega )\end{array},$

where F represents the level set function. Within the design domain D, the structure region is denoted by W and the design boundary is represented as G. $x$ signifies an arbitrary point in the design domain D and t is a pseudo time for the dynamic shape optimization process. This implicit boundary representation is illustrated in Fig. 2. The dynamic boundary evolution is governed by the Hamilton-Jacobi PDE as:

(7)${\displaystyle \frac{\partial \Phi }{\partial t}}-V_n\left|\nabla \Phi \right|=\mathrm{0},$

where the V_n is the normal design velocity field.

Fig.2 Level set model for implicit boundary representation. (a) The three-dimensional level set function; (b) the embedded boundary (highlight in red) as the zero level set

Full size|PPT slide

3.2 Parameterization of the level set function using a kernel function

With a given kernel function, the level set function can be written as the weighted summation form:

(8)$\mathrm{\Phi }(x)={\displaystyle \sum _{j=\mathrm{1}}^n{{\displaystyle \Psi }}_j(x){{\displaystyle \mu }}_j}\mathrm{.}$

In this study, ${\mu }_j$, representing $\mu (x_j)$ for simplicity, is the weight on the jth node. ${\mathrm{\Psi }}_j(x)$ (j=1,2,..,n) is the newly constructed CBF, with the Kronecker delta property as

(9)${\Psi }_j(x_i)=\{\begin{array}{l}\mathrm{1},\begin{array}{ll}& \left(\mathrm{if}\begin{array}{ll}& i=j\end{array}\right)\end{array}\\ \mathrm{0},\begin{array}{ll}& \left(\mathrm{if}\begin{array}{ll}& i\ne j\end{array}\right)\end{array}\end{array}\begin{array}{ll}& j=\mathrm{1},\mathrm{2},\mathrm{...},n\end{array}\mathrm{.}$

The readers can be referred to Ref. [62] for the complete construction procedure of the CBF kernel function. A comparison is made in Fig. 3 between a regular Gaussian RBF and the CBF. One thing needs to be mentioned is that with the analytical parametric level set function expression, the optimization results generated via PLSM can be manufactured via the high-resolution layer-image-based continuous liquid interface production (CLIP) 3D printing technology [63,73]. By bridging the PLSM design methodology and the CLIP printing process, the optimal structure can be manufactured at a lower prefabrication computational cost.

Fig.3 Illustration of two kernel functions. (a) The Gaussian RBF; (b) the newly constructed CBF

Full size|PPT slide

3.3 Distance-regularized level set function

To maintain a distance-regularized level set evolution, a distance regularization energy functional is introduced in the current study. Generally, for a given level set function $\Phi$, the distance regularization energy functional $R(\left|\nabla \Phi \right|)$ can be expressed as

(10)$R(\left|\nabla \Phi \right|)={\displaystyle \int {}_{\Omega }}P(\left|\nabla \Phi \right|)d\Omega ,$

where $P\left|\nabla \Phi \right|$ is the regularization energy potential.

To obtain a distance-regularized level set function, Li et al. [60] proposed the “double-well” regularization potential which takes the form:

(11)$P(\left|\nabla \mathrm{\Phi }\right|)=\{\begin{array}{l}{\displaystyle \frac{\mathrm{1}}{{(\mathrm{2}\pi )}^{\mathrm{2}}}}[\mathrm{1}-\cos (\mathrm{2}\pi \left|\nabla \mathrm{\Phi }\right|)],\begin{array}{ll}& \end{array}\left|\nabla \mathrm{\Phi }\right|<\mathrm{1}\\ {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{(\left|\nabla \mathrm{\Phi }\right|-\mathrm{1})}^{\mathrm{2}},\begin{array}{lll}\begin{array}{ll}\begin{array}{ll}& \end{array}& \end{array}& & \end{array}\left|\nabla \mathrm{\Phi }\right|\ge \mathrm{1}\end{array}\mathrm{.}$

The “double-well” potential $P(\left|\nabla \mathrm{\Phi }\right|)$, and its diffusive rate $D=d_p(\left|\nabla \mathrm{\Phi }\right|)$ are plotted in Fig. 4, where the diffusive rate has the form:

(12)$D=d_p(\left|\nabla \mathrm{\Phi }\right|)={\displaystyle \frac{P^{\prime }(\left|\nabla \mathrm{\Phi }\right|)}{\left|\nabla \mathrm{\Phi }\right|}}\mathrm{.}$

It can be seen from Fig. 4(b) that this “double-well” distance regularization effect will be affected by the “selecting point” at 0.5. With the minimization of the energy R, when the $\left|\nabla \mathrm{\Phi }\right|$ is below the selecting point 0.5, with a positive diffusive rate, $\left|\nabla \mathrm{\Phi }\right|$ tends to be driven to 0. Oppositely, when $\left|\nabla \mathrm{\Phi }\right|$ is above the selecting point 0.5 but below 1, with a negative diffusive rate, $\left|\nabla \mathrm{\Phi }\right|$ tends to be driven to 1. When $\left|\nabla \mathrm{\Phi }\right|$ is above 1, with a positive diffusive rate, $\left|\nabla \mathrm{\Phi }\right|$ tends to be driven to 1. This will regulate the level set function $\mathrm{\Phi }$ to a desired distance-regularized shape.

Fig.4 (a) The distance regularization energy P of Eq. (11); (b) the corresponding diffusive rate D

Full size|PPT slide

Easily, this selecting point at 0.5 can be modified based on needs. One of the constructing procedures was recently reported in Ref. [61]. By following a similar procedure and choose 0.25 as the “selecting point”, a new energy potential P_n used in this research can be built as

(13)$P_n(\left|\nabla \mathrm{\Phi }\right|)={\displaystyle \frac{{(\left|\nabla \mathrm{\Phi }\right|)}^{\mathrm{2}}}{\mathrm{2}}}-{\displaystyle \frac{\mathrm{3}}{\mathrm{64}}}\left[\mathrm{3}\ln (\left|\mathrm{64}{(\left|\nabla \mathrm{\Phi }\right|)}^{\mathrm{2}}-\mathrm{80}\left|\nabla \mathrm{\Phi }\right|+\mathrm{34}\right|)+\mathrm{10}\arctan {\displaystyle \frac{\mathrm{8}\left|\nabla \mathrm{\Phi }\right|-\mathrm{5}}{\mathrm{3}}}\right]\mathrm{.}$

As can be seen from Fig. 5, when the $\left|\nabla \mathrm{\Phi }\right|$ is below the selecting point 0.25, with a positive diffusive rate, $\left|\nabla \mathrm{\Phi }\right|$ tends to be driven to 0 and when $\left|\nabla \mathrm{\Phi }\right|$ is above the selecting point 0.25 but below 1, with a negative diffusive rate, $\left|\nabla \mathrm{\Phi }\right|$ tends to be driven to 1. When $\left|\nabla \mathrm{\Phi }\right|$ is above 1, with a positive diffusive rate, $\left|\nabla \mathrm{\Phi }\right|$ tends to be driven to 1. With a wider range for $\left|\nabla \mathrm{\Phi }\right|$ to be driven to 1, the corresponding level set function will be relatively more likely to be driven to a signed distance function near the design boundary. In this case, the shell region for the multiscale structure will be well maintained.

Fig.5 (a) The new distance regularization energy Pn of Eq. (13); (b) the corresponding diffusive rate Dn

Full size|PPT slide

To provide an intuitive process of imposing the distance regularization energy functional, the whole process is visualized in Fig. 6. By importing a black/white image (TO is short for topology optimization) as a binary initial level set function and minimizing the distance regularization energy according to it, the final distance-regularized level set function can be achieved, as shown in Fig. 6(c). It can be noticed that the overall “TO” shape is preserved. As can be seen from the zoom-in view in Figs. 6(e) and 6(d), the transition area of the level set function is regularized into a signed distance function while the rest areas are kept flat. The signed distance function can be used to discriminate the wrapping shell region from the interior infill region, which can provide an accurate mapping from the level set model to the physics model.

Fig.6 The distance regularization effect. TO: Topology optimization. (a) The initial binary image; (b) the binary value level set function based on the binary image; (c) the distance regularized level set function; (d) the zoom-in view of the binary level set function; (e) the zoom-in view of the distance regularized level set function

Full size|PPT slide

4 General concurrent optimizing settings for multiscale structure design

To formulate the proposed multiscale concurrent structural topology optimization in a general way, the level set nodal height $\mu$, together with the infill material properties (Young’s modulus E and the Poisson’s ratio $\nu$) are set to be the design variables. Within the distance regularized parametric level set framework:

(14)$\begin{array}{l}\underset{\mu ,E,\nu }{\mathrm{Minimize}}:\begin{array}{l}\begin{array}{l}\end{array}{J}_{\mathrm{tot}}\end{array}=\mathrm{Objective}\hspace{0.17em}\mathrm{function}+wR\\ \mathrm{Subject}\hspace{0.17em}\mathrm{to:}\begin{array}{l}\end{array}\mathrm{Governing}\hspace{0.17em}\mathrm{equations,}\\ \mathrm{Design}\hspace{0.17em}\mathrm{constraints,}\end{array}$

where the total objective function J_tot contains the distance regularization energy functional R with weight w.

To discriminate a wrapping shell region and the infill region from the macroscale, the corresponding material interpolation from the level set model to the physics model needs to be modified. Conventionally, the material property interpolation is done by directly getting the Heaviside function value of the level set function. This process will generate a binary field in order to separate the material region from the void region, as shown in Fig. 7(a) and 7(b). In the proposed material interpolation scheme, a second Heaviside function is calculated within the material region at a given transition width $\Delta$. In this case, the design domain will be separated into three phases, which can be used to represent the void, the wrapping boundary layer, and the interior infill region of the multiscale structure, as shown in Fig. 7(d).

Fig.7 The material property interpolation. (a) The conventional level set implicit boundary representation by the level set function and the zero level-set; (b) the material region $\mathrm{\Omega }$, the design boundary $\mathrm{\Gamma }$ and the design domain D represented by the conventional level set model; (c) the proposed shell-infill representation by the level set function and the level $\mathrm{\Phi }=\mathrm{0}$ and $\mathrm{\Phi }=\Delta$ ($\Delta$ is the selected shell width); (d) within the material region $\mathrm{\Omega }$, the shell region ${\mathrm{\Omega }}_s$ and the infill region ${\mathrm{\Omega }}_i$ are bounded by the shell boundary ${\mathrm{\Gamma }}_s$ and the infill boundary ${\mathrm{\Gamma }}_i$ in the design domain D, respectively

Full size|PPT slide

With the gradient-based optimizer MMA, the optimization process requires the calculation of sensitivities. The sensitivities with respect to $\mu$ will be derived based on different optimization problems. The sensitivities for the material properties can be analytically calculated via a forward finite difference scheme. Generally, the derivative of a given function f at point x can be defined by the limit:

(15)$f^{\prime }(x)=\underset{h\to \mathrm{0}}{lim}{\displaystyle \frac{f(x+h)-f(x)}{h}}\mathrm{.}$

After the concurrent topology optimization process, the optimal macroscale structural topology and the optimal infill material properties can be obtained. To find out the actual layout of the infill structure, a second optimization can be set up in a least square manner. An extra isotropic condition is also included as an optimization constraint to ensure the isotropy of the infill structure.

(16)$\begin{array}{l}\underset{\mu }{\mathrm{M}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{m}\mathrm{i}\mathrm{z}\mathrm{e}}:\begin{array}{l}\end{array}{J}_{\mathrm{infill}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \sum _{ijkl}^n{(C_{ijkl}^H-C_{ijkl}^{*})}^{\mathrm{2}}}+w^{\prime }R\\ \mathrm{Subject}\hspace{0.17em}\mathrm{to}\begin{array}{l}\end{array}{\displaystyle {\int }_{D}H(\mathrm{\Phi })\mathrm{d}\mathrm{\Omega }}=V_{\mathrm{const}},\\ a(u,v)=l(v)\\ C_{\mathrm{1212}}^{*}=(C_{\mathrm{1111}}^H+C_{\mathrm{2222}}^H)/\mathrm{4}-C_{\mathrm{1122}}^H/\mathrm{2},\end{array}$

where $C_{ijkl}^H$ is the component of the homogenized elasticity tensor and their optimal value targets, got from the concurrent optimization, are noted as $C_{ijkl}^{*}$, $H(\mathrm{\Phi })$ is the Heaviside function of the level set $\mathrm{\Phi }$, V_const is the design volume constraint, u is the displacement field, and v is the test function. The weight of the distance regularization energy functional R is noted as $w^{\prime }$. The bilinear energy form a(u,v) and the linear load form l(u) inside the computational domain D, with p representing the body force and $\tau$ representing the boundary traction, are described as

(17)$\{\begin{array}{l}a(u,v)={\displaystyle {\int }_{\mathrm{\Omega }}{E}_{ijkl}{\epsilon }_{ij}(u){\epsilon }_{kl}(v)\mathrm{d}\mathrm{\Omega }}\\ l(v)={\displaystyle {\int }_{\mathrm{\Omega }}p\cdot v\mathrm{d}\mathrm{\Omega }}+{\displaystyle {\int }_{\Gamma }\tau \cdot v\mathrm{d}s}\end{array}\mathrm{.}$

The structural isotropy [68] is ensured when

(18)$C_{\mathrm{1212}}^{*}=(C_{\mathrm{1111}}^H+C_{\mathrm{2222}}^H)/\mathrm{4}-C_{\mathrm{1212}}^H/\mathrm{2.}$

Essentially, the isotropy is ensured when the current $C_{\mathrm{1212}}^H$ hits the targeted $C_{\mathrm{1212}}^{*}$ calculated from Eq. (18).

One more thing that needs to be mentioned is the monitoring of the isotropy during the optimization. Ideally, when this least square optimization comes to its convergence, the selected elastic tensor entries will hit the targeted values. However, in practice, this condition may not be satisfied 100%. But even without hitting all the selected targets, the isotropy may still be valid as long as the condition in Eq. (18) is hold. An isotropy polar plot index [70] is utilized in this research to visualize the isotropy of the infill structure throughout its optimization process, for an easy isotropy monitoring. The details are introduced as follows: For a unit cell with a rotation of angle $\theta$, as shown in Fig. 8, the original elastic tensor will be transformed as shown in Eq. (19) and the corresponding expressions are listed in Eq. (20).

(19)$\left[\begin{array}{lll}{C}_{\mathrm{1111}}^H& C_{\mathrm{1122}}^H& \mathrm{0}\\ C_{\mathrm{2211}}^H& C_{\mathrm{2222}}^H& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& C_{\mathrm{1212}}^H\end{array}\right]\Rightarrow \left[\begin{array}{lll}{\overline{C}}_{\mathrm{1111}}^H& {\overline{C}}_{\mathrm{1122}}^H& \mathrm{0}\\ {\overline{C}}_{\mathrm{2211}}^H& {\overline{C}}_{\mathrm{2222}}^H& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& {\overline{C}}_{\mathrm{1212}}^H\end{array}\right],$

(20)$\{\begin{array}{l}{\overline{C}}_{\mathrm{1111}}^H=C_{\mathrm{1111}}^H{\cos }^{\mathrm{4}}\theta +C_{\mathrm{2222}}^H{\sin }^{\mathrm{4}}\theta +\mathrm{2}(C_{\mathrm{1122}}^H+\mathrm{2}{C}_{\mathrm{1212}}^H){\sin }^{\mathrm{2}}\theta {\cos }^{\mathrm{2}}\theta \\ {\overline{C}}_{\mathrm{2222}}^H=C_{\mathrm{1111}}^H{\sin }^{\mathrm{4}}\theta +C_{\mathrm{2222}}^H{\cos }^{\mathrm{4}}\theta +\mathrm{2}(C_{\mathrm{1122}}^H+\mathrm{2}{C}_{\mathrm{1212}}^H){\sin }^{\mathrm{2}}\theta {\cos }^{\mathrm{2}}\theta \\ {\overline{C}}_{\mathrm{1122}}^H=\left({\overline{C}}_{\mathrm{1111}}^H+{\overline{C}}_{\mathrm{2222}}^H-\mathrm{4}{\overline{C}}_{\mathrm{1212}}^H\right){\sin }^{\mathrm{2}}\theta {\cos }^{\mathrm{2}}\theta +{\overline{C}}_{\mathrm{1122}}^H({\sin }^{\mathrm{4}}\theta +{\cos }^{\mathrm{4}}\theta )\\ {\overline{C}}_{\mathrm{1212}}^H=\left({\overline{C}}_{\mathrm{1111}}^H+{\overline{C}}_{\mathrm{2222}}^H-\mathrm{2}{\overline{C}}_{\mathrm{1122}}^H-\mathrm{2}{\overline{C}}_{\mathrm{1212}}^H\right){\sin }^{\mathrm{2}}\theta {\cos }^{\mathrm{2}}\theta +{\overline{C}}_{\mathrm{1212}}^H({\sin }^{\mathrm{4}}\theta +{\cos }^{\mathrm{4}}\theta )\end{array}\mathrm{.}$

Fig.8 The coordinates change of a unit cell

Full size|PPT slide

By plotting the normalized value of each expression in Eq. (20) with $\theta$ changing from 0° to 360°, an ideal isotropic material will generate a circle with the radius of 1. If the material is not isotropic, the plot will deviate away from the standard circle. This is illustrated in Fig. 9 where the red dashed line circle with a radius of 1 represents a plot of an isotropic material while the blue curve represents an anisotropic one. By plotting the current material with blue curve to be compared with the stand red circle throughout the optimization process, the isotropy of the structure can be monitored.

Fig.9 The isotropy polar plot

Full size|PPT slide

Another detail that needs to be mentioned is how to determine the bounds for the infill structure material properties. The bounds should be chosen reasonably to get a valid design. However, since the infill structure only occupies a portion of the total volume of the design domain, the overall Young’s modulus in each direction would be lower than the constructing material itself [74]. In the current research, the bounds are set by using Ref. [27] as a reference, where microstructures were designed to reach maximum material properties.

At last, a systematic flowchart is drawn in Fig. 10 for a more intuitive representation of the entire proposed concurrent topology optimization process.

Fig.10 The flow chart for the concurrent topology optimization process

Full size|PPT slide

5 Application examples and numerical results

5.1 The concurrent design of a Michell-type structure and its infill material property

In this section, a two-by-one design domain is used for the designing of a multiscale Michell-type structure with fixed-fixed supports and multiple loads. The design domain is shown in Fig. 11. The forces F₁=1 and F₂=0.5 are applied at the bottom edge. The overall structure optimization statement can be made as

(21)$\begin{array}{l}\underset{u,E,\nu }{\mathrm{M}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{m}\mathrm{i}\mathrm{z}\mathrm{e}}:\begin{array}{l}\end{array}OBJ={\displaystyle {\int }_D{E}_{ijkl}{\epsilon }_{ij}(u){\epsilon }_{kl}(v)}H(\mathrm{\Phi })\mathrm{d}\mathrm{\Omega }+w_{\mathrm{1}}R\\ \mathrm{Subject}\hspace{0.17em}\mathrm{to:}\begin{array}{l}\end{array}{\displaystyle {\int }_{D}H(\mathrm{\Phi })\mathrm{d}\mathrm{\Omega }}=V_{\mathrm{const}},\\ \begin{array}{llll}& & & \end{array}a(u,v)=l(v).\end{array}$

where $E_{ijkl}{\epsilon }_{ij}(u){\epsilon }_{kl}(v)$ represents the strain energy density inside the design domain. The sensitivity analysis can be taken from Ref. [62], which takes the form of Eq. (22). In these equations, a smoothed Dirac delta function is represented as $\delta (\mathrm{\Phi })$. As can be noted, the shape derivative is extended to the whole design domain naturally.

Fig.11 The boundary condition for designing the Michell-type structure with fixed-fixed supports and multiple loads

Full size|PPT slide

In this example, the design domain is discretized into 100×50 quadrilateral elements. The Young’s modulus for the shell material is set to be 1 and the void to be 0.001. As for the infill region, the range for the E varies from 0.01 to 0.165 and the range for the Poisson’s ratio varies from 0.3 to 0.4. The minimum mean compliance optimization is carried out under the total volume constraint of 60%. The optimization convergence history, together with the level set evolution and the design evolution are shown in Fig. 12. The optimal infill material properties E is 0.165 with infill material Poisson’s ratio 0.3. The overall objective function is minimized to 12.4199 and the final total volume is 0.60246.

(22)$\{\begin{array}{l}{\displaystyle \frac{\partial {\displaystyle {\int }_D{E}_{ijkl}{\epsilon }_{ij}(u){\epsilon }_{kl}(v)}H(\mathrm{\Phi })\mathrm{d}\mathrm{\Omega }}{\partial {\mu }_j}}=-{\displaystyle {\int }_D{E}_{ijkl}{\epsilon }_{ij}(u){\epsilon }_{kl}(v)\delta (\mathrm{\Phi }){\mathrm{\Psi }}_j(x)}\mathrm{d}\mathrm{\Omega }\\ {\displaystyle \frac{\partial R}{\partial {\mu }_j}}=-\mathrm{div}[d_p(\left|\nabla \mathrm{\Phi }\right|\nabla \mathrm{\Phi })]{\mathrm{\Psi }}_j(x)\\ {\displaystyle \frac{\partial V_{\mathrm{const}}}{\partial {\mu }_j}}={\displaystyle {\int }_D\delta (\mathrm{\Phi }){\mathrm{\Psi }}_j(x)}\mathrm{d}\mathrm{\Omega }\end{array},j=\mathrm{1},\mathrm{2},\mathrm{...},n\mathrm{.}$

Fig.12 Convergence history and the design evolution for the macroscale Michell-type structure with fixed-fixed supports and multiple loads. Upper: The level set function evolution; lower: The corresponding design evolution

Full size|PPT slide

With the optimal infill material property targets, the infill structure optimization can be carried out following Eq. (16). The optimization convergence history, the evolution-history for the level set function, the design evolution and the final design polar isotropy plot are shown in Fig. 13. It can be seen from Fig. 13(a) that with the least square objective function minimized, the selected targets are hit. The isotropy of the infill structure can be verified intuitively from Fig. 13(b) where the blue-color polar plot of the current infill structure is almost overlapping with the reference red unit circle at the end of the optimization process.

Fig.13 Convergence history and the isotropy of the infill structure. (a) The convergence history and design evolution for the infill structure: The level set function evolution (upper) and the corresponding design evolution (lower); (b) the isotropic polar plot of the corresponding infill structure throughout the optimization process: Reference standard circle (red) and isotropic polar plot of the current infill structure (blue); (c) the elastic tensor of the corresponding infill structure

Full size|PPT slide

5.2 Concurrent designing of an NPR structure and its infill material property

To achieve a structure with desired Poisson’s ratio, a common approach is to formulate a least square problem to design a structure with a specified elastic tensor. However, from the authors’ experience, trying to hit all targets simultaneously requires fine adjustment of parameters inside the optimization algorithm, which is tricky in most of the times. To decrease the number of adjustable weights, NPR structure design can be formulated directly as

(23)$\begin{array}{l}\underset{\mu ,E,\nu }{\mathrm{Minimize}}:\begin{array}{l}\end{array}{J}_{\mathrm{1}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left({\displaystyle \frac{C_{\mathrm{1122}}^H}{C_{\mathrm{1111}}^H}}-{\nu }_t\right)}^{\mathrm{2}}+w_{\mathrm{2}}R\\ \mathrm{Subject\; to}:\begin{array}{l}\end{array}{\displaystyle {\int }_{D}H(\Phi )d\Omega }=V_{\mathrm{const}},\\ \begin{array}{llll}& & & \end{array}a(u,v)=l(v),\\ \begin{array}{llll}& & & \end{array}{C}_{\mathrm{1111}}^H=C_{\mathrm{2222}}^H,\end{array}$

where ${\nu }_t$ is the Poisson’s ratio target. To ensure an overall orthotropic structure behavior, C₁₁₁₁ and C₂₂₂₂ are constrained to be equal. Considering that the homogenized elastic tensor entries can be isolated out by the aforementioned strain energy method, the design sensitivities can be borrowed from Eq. (22).

In this example, a 1×1 domain is discretized into 50×50 quadrilateral elements. The Young’s modulus for the void region is 0.001 and 1 for the wrapping shell, respectively. For the infill structure, the Poisson’s ratio is given a range from 0.4 to 0.5. The range of the Young’s modulus of the infill structure is from 0.01 to 0.13. The total volume is constrained within the range from 49% to 51%. The overall Poisson’s ratio target is set to be ‒0.3. The design evolution and the convergence history are shown in Fig. 14. The final total volume ratio is 0.5098 and the final overall Poisson’s ratio is ‒0.25841. For the infill structure, the final Young’s modulus is 0.13 and Poisson’s ratio is 0.4.

Fig.14 Convergence history and the design evolution for the macroscale NPR structure. Upper: The level set function evolution; Lower: The corresponding design evolution

Full size|PPT slide

With the infill material properties achieved, the second-stage topology optimization is formulated in the least square manner. The design domain and the material properties for void and material are kept the same. The design evolution and the convergence history for the infill structure are shown in Fig.15.

Furthermore, the NPR performance for the entire multiscale structure is verified in Fig. 16. A 0.3 downwards displacement is added on a 4×4 multiscale NPR structure array. The boundary conditions are shown in Fig. 16(a) and the pre/post-deformation simulation results for the multiscale structure are shown in Figs. 16(b) and 16(c), respectively. A shrink from the right side of the structure array can be observed from Fig. 16(c), which verifies the NPR performance.

Fig.15 Convergence history and the isotropy of the infill structure. (a) The convergence history and design evolution for the infill structure: The level set function evolution (Upper) and the corresponding design evolution (lower); (b) the isotropic polar plot of the infill structure throughout the optimization process: Reference standard circle (red) and isotropic polar plot of the final infill structure (blue); (c) the elastic tensor of the corresponding infill structure

Full size|PPT slide

Fig.16 The NPR effect verification. (a) A 3×3 array of the macroscale structure; (b) one unit cell of the macroscale structure (the outer shell is shown in blue and the inner infill is shown in green); (c) one unit cell of the infill structure

Full size|PPT slide

6 Multi-control-point conformal mapping, the mapped optimization results and its manufacturing

With the achieved optimal overall structure and the optimal infill structure, in this section, these two parts are integrated together via the local shape-preserving conformal mapping [75]. The shape preserving effect can be visualized in Fig. 17, as the small circles on a human face will remain their circular shape after the conformal mapping.

Fig.17 The local shape preserving effect of conformal mapping

Full size|PPT slide

Firstly, the concept of conformal mapping is introduced here in smooth settings. Let $\omega =f(z):ℂ\to ℂ$ be a complex function on the plane. Denote

(24)${\displaystyle \frac{\partial }{\partial \overline{z}}}:={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\displaystyle \frac{\partial }{\partial x}}+i{\displaystyle \frac{\partial }{\partial y}}\right),$

where i is the unit imaginary root. Then f is said to be conformal if

(25)${\displaystyle \frac{\partial f}{\partial \overline{z}}}=\mathrm{0.}$

In this research, the local shape-preserving character can ensure that the infill structure will maintain its designed material properties. This can be a great advantage compared with putting the infill structure into the design in a periodic manner and trim the ones on the design boundary. In that case, the elements cut by the design boundary might result in overhanging structure, isolated islands and some other unwanted results. Those defects will sabotage the designed infill structural properties and bring challenges to the manufacturing process.

Under discrete settings, conformal mappings can be computed by discrete Ricci flow method [76–78]. Some more discrete Ricci flow algorithms regarding efficiency and adaptivity improvements are reported in Refs. [79,80]. For further conformal mapping algorithms, the readers are referred to Ref. [81] for more information.

Given a triangular mesh ${\displaystyle \sum =(V,F,E)}$, a face element is denoted with corner vertex v_i, v_j, and v_k by f_ijk, and the angle between rays are denoted as $i\overrightarrow{j}$ and $\overrightarrow{i}k$ by ${\theta }_i^{jk}$. Then the discrete Gaussian curvature at vertex v_i is defined by

(26)$K_i=\{\begin{array}{l}\mathrm{2}\pi -{\sum }_{f_{ijk}\in F}{\theta }_i^{jk}\begin{array}{l}\end{array}if\hspace{0.17em}{v}_i\notin \partial \sum \\ \pi -{\sum }_{f_{ijk}\in F}{\theta }_i^{jk}\begin{array}{l}\end{array}if\hspace{0.17em}{v}_i\in \partial \sum \end{array},$

where $\partial \sum$ is the boundary of mesh $\Sigma$. Now the discrete Ricci flow is defined as follows. Given a circle packing metric to $\Sigma$, i.e., on each vertex v_i a positive real number is defined as ${{\displaystyle \gamma }}_i$, then the edge length between vertices v_i and v_j is $l_{ij}={\gamma }_i+{\gamma }_j$. With those parameters, all angles can be calculated in $\Sigma$. Denote u_i=logγ_i, then the discrete Ricci flow is defined as

(27)${\displaystyle \frac{d{u}_i(t)}{dt}}={\overline{K}}_i-K_i,$

where $\bar{k}={\left({\bar{K}}_{\mathrm{1}},{\bar{K}}_{\mathrm{2}},\mathrm{...},{\bar{K}}_n\right)}^T$ is the user defined target curvature.

In our case, a conformal mapping from an irregular planar region to a polygonal region needs to be found such that the inner angles are either $\pi /\mathrm{2}$ or $\mathrm{3}\pi /\mathrm{2}$. The polygonal region is filled with regular infill structures, and then they are mapped back with the inverse of the computed conformal mapping, which is also a conformal mapping that preserves local shapes. To realize this, on the boundary of input mesh $\partial \Sigma$, multiple control points

$W:=\{w_{\mathrm{1}},w_{\mathrm{2}},\mathrm{...},w_k\}\subseteq \partial \Sigma$

can be selected based on the need. Then the target curvature on each vertex are defined by

(28)${\overline{K}}_i=\{\begin{array}{l}\mathrm{0}\begin{array}{llll}& & & \end{array}\begin{array}{lll}& & \end{array}\mathrm{if}\begin{array}{l}\end{array}{v}_i\notin \partial \Sigma \\ \mathrm{0}\begin{array}{llll}& & & \end{array}\begin{array}{lll}& & \end{array}\mathrm{if}\begin{array}{l}\end{array}{v}_i\notin \partial \Sigma \backslash W\\ -\pi /\mathrm{2}\begin{array}{l}\end{array}\mathrm{or}\begin{array}{l}\end{array}\pi /\mathrm{2}\begin{array}{ll}& \end{array}\mathrm{if}\begin{array}{l}\end{array}{v}_i\notin W\end{array},$

where ${\bar{K}}_i=\pi /\mathrm{2}$ is chosen if the target polygonal region has an outward right angle and ${\bar{K}}_i=-\pi /\mathrm{2}$ is chosen if the target polygonal region has an inward right angle at point v_i. Compared to Ricci flow method with only four control points, as shown in Fig. 18, which will eventually map the input region to a rectangular region, the proposed multi-control-point method provides more flexibility with the benefit of lower area distortion. The proposed mapping process is illustrated in Fig. 19 using the mapping of the Mitchell-type structure as an example. Considering the symmetry, only the right half of the structure is taken into the mapping procedure and the final structure, as shown in Fig. 19(i) is achieved by a simple mirroring technique. Following the proposed mapping process, both the mapped Mitchell-type structure and the NPR structure are shown in Figs. 20 and 21, respectively. It is worth noting that after the mapping, the infill and the shell are integrated together. However, the actual material used for the final multiscale structure should be less than its designed total volume ratio since the mapped infill structure is porous. Practically, the actual volume ratio for the entire multiscale structure should be the macroscale structural volume ratio times the infill micro/meso scale structural volume ratio.

Fig.18 The conventional conformal mapping Ricci flow method with four control points

Full size|PPT slide

Fig.19 The multi-control-point conformal mapping process. (a) The shell-infill structure optimization result; (b) half of the original design; (c) the infill region; (d) the meshed infill region with multiple control points (red) and one central point (green); (e) the zoom-in view of the triangular mesh; (f) the isotropic infill unit cell; (g) the mapped infill structure; (h) half of the shell-infill mapped structure; (i) the final shell-infill multiscale structure

Full size|PPT slide

Fig.20 The conformal mapping result for the Mitchell-type structure

Full size|PPT slide

Fig.21 The conformal mapping result for the NPR structure

Full size|PPT slide

With the mapped structural design, the final manufacturing process can be carried out. The mapped result is sent to a fused deposition modeling (FDM) 3D printer and manufactured by using polylactic acid (PLA) filament. The printed result is shown in Fig. 22.

Fig.22 The 3D printed multiscale Mitchell-type structure. (a) The computational result; (b) the 3D printed result

Full size|PPT slide

7 Conclusions

In this paper, a concurrent CBF-based PLSM topology optimization is proposed for the designing of multiscale structures. By using the newly constructed CBF, the explicit design variable bounds can be passed to the optimizer, which is a significant advantage over the RBF-based PLSM. Driven by MMA, the proposed approach can efficiently handle multiple constraints. Coupled with mathematical programming, the multiscale structural topology optimization can be directly carried out in a concurrent manner. The distance-regularized level set function is used to discriminate the shell and the infill region from the macro scale. From the micro/meso scale, the isotropic infill structure will bring advantages to the structural mapping/manufacturing process. With the help of the local shape-preserving conformal mapping, the infill material properties can be mathematically preserved. By introducing multiple control points, the current proposed conformal mapping can be more flexible and adaptive in handling complex geometries. The mapped multiscale structure can be directly sent to 3D printers for manufacturing, to conclude the design-mapping-manufacturing multiscale structure topology optimization process.

References

[1] Gibson L J, Ashby M F. Cellular Solids: Structure and Properties. 2nd ed. Cambridge: Cambridge University Press, 1999

[2] Christensen R M. Mechanics of cellular and other low-density materials. International Journal of Solids and Structures, 2000, 37(1–2): 93–104

[3] Zheng X, Lee H, Weisgraber T H, Ultralight, ultrastiff mechanical metamaterials. Science, 2014, 344(6190): 1373–1377

[4] Valdevit L, Jacobsen A J, Greer J R, Protocols for the optimal design of multi-functional cellular structures: From hypersonics to micro-architected materials. Journal of the American Ceramic Society, 2011, 94(s1): s15–s34

[5] Clausen A, Wang F, Jensen J S, Topology optimized architectures with programmable Poisson’s ratio over large deformations. Advanced Materials, 2015, 27(37): 5523–5527

[6] Schwerdtfeger J, Wein F, Leugering G, Design of auxetic structures via mathematical optimization. Advanced Materials, 2011, 23(22–23): 2650–2654

[7] Murr L E, Gaytan S M, Ramirez D A, Metal fabrication by additive manufacturing using laser and electron beam melting technologies. Journal of Materials Science and Technology, 2012, 28(1): 1–14

[8] Han S C, Lee J W, Kang K. A new type of low density material: Shellular. Advanced Materials, 2015, 27(37): 5506–5511

[9] Xia Q, Wang M Y. Simultaneous optimization of the material properties and the topology of functionally graded structures. Computer Aided Design, 2008, 40(6): 660–675

[10] Xia L, Breitkopf P. Concurrent topology optimization design of material and structure within FE² nonlinear multiscale analysis framework. Computer Methods in Applied Mechanics and Engineering, 2014, 278: 524–542

[11] Sigmund O. Design of material structures using topology optimization. Dissertation for the Doctoral Degree. Lyngby: Technical University of Denmark, 1994

[12] 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

[13] Fujii D, Chen B, Kikuchi N. Composite material design of two-dimensional structures using the homogenization design method. International Journal for Numerical Methods in Engineering, 2001, 50(9): 2031–2051

[14] Neves M M, Sigmund O, Bendsøe M P. Topology optimization of periodic microstructures with a penalization of highly localized buckling modes. International Journal for Numerical Methods in Engineering, 2002, 54(6): 809–834

[15] Bendsøe M P. Optimal shape design as a material distribution problem. Structural Optimization, 1989, 1(4): 193–202

[16] Rozvany G I, Zhou M, Birker T. Generalized shape optimization without homogenization. Structural Optimization, 1992, 4(3–4): 250–252

[17] Sethian J A, Wiegmann A. Structural boundary design via level set and immersed interface methods. Journal of Computational Physics, 2000, 163(2): 489–528

[18] Wang M Y, Wang X, Guo D. A level set method for structural topology optimization. Computer Methods in Applied Mechanics and Engineering, 2003, 192(1–2): 227–246

[19] Allaire G, Jouve F, Toader A M. Structural optimization using sensitivity analysis and a level-set method. Journal of Computational Physics, 2004, 194(1): 363–393

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

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

[22] Sigmund O. Tailoring materials with prescribed elastic properties. Mechanics of Materials, 1995, 20(4): 351–368

[23] Sigmund O, Torquato S. Design of materials with extreme thermal expansion using a three-phase topology optimization method. Journal of the Mechanics and Physics of Solids, 1997, 45(6): 1037–1067

[24] Mei Y, Wang X. A level set method for structural topology optimization and its applications. Advances in Engineering Software, 2004, 35(7): 415–441

[25] Vogiatzis P, Chen S, Wang X, Topology optimization of multi-material negative Poisson’s ratio metamaterials using a reconciled level set method. Computer Aided Design, 2017, 83: 15–32

[26] Wang Y, Luo Z, Zhang N, Topological shape optimization of microstructural metamaterials using a level set method. Computational Materials Science, 2014, 87: 178–186

[27] Huang X, Radman A, Xie Y. Topological design of microstructures of cellular materials for maximum bulk or shear modulus. Computational Materials Science, 2011, 50(6): 1861–1870

[28] Huang X, Xie Y M, Jia B, Evolutionary topology optimization of periodic composites for extremal magnetic permeability and electrical permittivity. Structural and Multidisciplinary Optimization, 2012, 46(3): 385–398

[29] Huang X, Zhou S, Xie Y, Topology optimization of microstructures of cellular materials and composites for macrostructures. Computational Materials Science, 2013, 67: 397–407

[30] Xia L, Xia Q, Huang X, Bi-directional evolutionary structural optimization on advanced structures and materials: A comprehensive review. Archives of Computational Methods in Engineering, 2018, 25(2): 437–478

[31] Xia L, Breitkopf P. Recent advances on topology optimization of multiscale nonlinear structures. Archives of Computational Methods in Engineering, 2017, 24(2): 227–249

[32] Rodrigues H, Guedes J M, Bendsøe MP. Hierarchical optimization of material and structure. Structural and Multidisciplinary Optimization, 2002, 24(1): 1–10

[33] Coelho P, Fernandes P, Guedes J, A hierarchical model for concurrent material and topology optimisation of three-dimensional structures. Structural and Multidisciplinary Optimization, 2008, 35(2): 107–115

[34] Liu L, Yan J, Cheng G. Optimum structure with homogeneous optimum truss-like material. Computers & Structures, 2008, 86(13–14): 1417–1425

[35] Deng J, Chen W. Concurrent topology optimization of multiscale structures with multiple porous materials under random field loading uncertainty. Structural and Multidisciplinary Optimization, 2017, 56(1): 1–19

[36] Wu J, Clausen A, Sigmund O. Minimum compliance topology optimization of shell-infill composites for additive manufacturing. Computer Methods in Applied Mechanics and Engineering, 2017, 326: 358–375

[37] Wu J, Aage N, Westermann R, Infill optimization for additive manufacturing approaching bone-like porous structures. IEEE Transactions on Visualization and Computer Graphics, 2018, 24(2): 1127–1140

[38] Deng J, Yan J, Cheng G. Multi-objective concurrent topology optimization of thermoelastic structures composed of homogeneous porous material. Structural and Multidisciplinary Optimization, 2013, 47(4): 583–597

[39] Sivapuram R, Dunning P D, Kim H A. Simultaneous material and structural optimization by multiscale topology optimization. Structural and Multidisciplinary Optimization, 2016, 54(5): 1267–1281

[40] Wang Y, Chen F, Wang M Y. Concurrent design with connectable graded microstructures. Computer Methods in Applied Mechanics and Engineering, 2017, 317: 84–101

[41] Li H, Luo Z, Zhang N, Integrated design of cellular composites using a level-set topology optimization method. Computer Methods in Applied Mechanics and Engineering, 2016, 309: 453–475

[42] Li H, Luo Z, Gao L, Topology optimization for concurrent design of structures with multi-patch microstructures by level sets. Computer Methods in Applied Mechanics and Engineering, 2018, 331: 536–561

[43] Yan X, Huang X, Zha Y, Concurrent topology optimization of structures and their composite microstructures. Computers & Structures, 2014, 133: 103–110

[44] Yan X, Huang X, Sun G, Two-scale optimal design of structures with thermal insulation materials. Composite Structures, 2015, 120: 358–365

[45] Da D, Cui X, Long K, Concurrent topological design of composite structures and the underlying multi-phase materials. Computers & Structures, 2017, 179: 1–14

[46] Vicente W, Zuo Z, Pavanello R, Concurrent topology optimization for minimizing frequency responses of two-level hierarchical structures. Computer Methods in Applied Mechanics and Engineering, 2016, 301: 116–136

[47] Osher S, Sethian J A. Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations. Journal of Computational Physics, 1988, 79(1): 12–49

[48] Osher S J, Santosa F. Level set methods for optimization problems involving geometry and constraints: I. Frequencies of a two-density inhomogeneous drum. Journal of Computational Physics, 2001, (1): 272–288

[49] Allaire G, Jouve F, Toader A M. A level-set method for shape optimization. Comptes Rendus Mathematique, 2002, 334(12): 1125–1130

[50] Luo Z, Tong L, Wang M Y, Shape and topology optimization of compliant mechanisms using a parameterization level set method. Journal of Computational Physics, 2007, 227(1): 680–705

[51] Luo Z, Wang M Y, Wang S, A level set-based parameterization method for structural shape and topology optimization. International Journal for Numerical Methods in Engineering, 2008, 76(1): 1–26

[52] Jiang L, Chen S. Parametric structural shape & topology optimization with a variational distance-regularized level set method. Computer Methods in Applied Mechanics and Engineering, 2017, 321: 316–336

[53] Wang S, Wang M Y. Radial basis functions and level set method for structural topology optimization. International Journal for Numerical Methods in Engineering, 2006, 65(12): 2060–2090

[54] Buhmann M D. Radial Basis Functions: Theory and Implementations. Vol. 12. Cambridge: Cambridge University Press, 2003

[55] Svanberg K. MMA and GCMMA-two methods for nonlinear optimization. 2007. Retrieved from https://people.kth.se/~krille/mmagcmma.pdf.

[56] Chen J, Shapiro V, Suresh K, Shape optimization with topological changes and parametric control. International Journal for Numerical Methods in Engineering, 2007, 71(3): 313–346

[57] Qian X. Full analytical sensitivities in NURBS based isogeometric shape optimization. Computer Methods in Applied Mechanics and Engineering, 2010, 199(29–32): 2059–2071

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

[59] Safdari-Vaighani A, Heryudono A, Larsson E. A radial basis function partition of unity collocation method for convection-diffusion equations arising in financial applications. Journal of Scientific Computing, 2015, 64(2): 341–367

[60] Li C, Xu C, Gui C, Distance regularized level set evolution and its application to image segmentation. IEEE Transactions on Image Processing, 2010, 19(12): 3243–3254

[61] Zhu B, Wang R, Li H, A level set method with a bounded diffusion for structural topology optimization. Journal of Mechanical Design, 2018, 140(7): 071402

[62] Jiang L, Chen S, Jiao X. Parametric shape and topology optimization: A new level set approach based on cardinal basis functions. International Journal for Numerical Methods in Engineering, 2018, 114(1): 66–87

[63] Jiang L, Ye H, Zhou C, Parametric topology optimization toward rational design and efficient prefabrication for additive manufacturing. In: Proceedings of ASME 2017 12th International Manufacturing Science and Engineering Conference Collocated with the JSME/ASME 2017 6th International Conference on Materials and Processing. Los Angeles: ASME, 2017, V004T05A006

[64] Wang Y, Kang Z. A level set method for shape and topology optimization of coated structures. Computer Methods in Applied Mechanics and Engineering, 2018, 329: 553–574

[65] Fu J, Li H, Xiao M, Topology optimization of shell-infill structures using a distance regularized parametric level-set method. Structural and Multidisciplinary Optimization, 2019, 59(1):249–262 (in press)

[66] Tushinsky L, Kovensky I, Plokhov A, Coated Metal: Structure and Properties of Metal-Coating Compositions. Berlin: Springer, 2013

[67] Challis V, Roberts A, Wilkins A. Design of three dimensional isotropic microstructures for maximized stiffness and conductivity. International Journal of Solids and Structures, 2008, 45(14–15): 4130–4146

[68] Neves M, Rodrigues H, Guedes J M. Optimal design of periodic linear elastic microstructures. Computers & Structures, 2000, 76(1–3): 421–429

[69] Radman A, Huang X, Xie Y. Topological optimization for the design of microstructures of isotropic cellular materials. Engineering Optimization, 2013, 45(11): 1331–1348

[70] Guth D, Luersen M, Muñoz-Rojas P. Optimization of three-dimensional truss-like periodic materials considering isotropy constraints. Structural and Multidisciplinary Optimization, 2015, 52(5): 889–901

[71] Vogiatzis P, Ma M, Chen S, Computational design and additive manufacturing of periodic conformal metasurfaces by synthesizing topology optimization with conformal mapping. Computer Methods in Applied Mechanics and Engineering, 2018, 328: 477–497

[72] Sethian J A. Theory, algorithms, and applications of level set methods for propagating interfaces. Acta Numerica, 1996, 5: 309–395

[73] Tumbleston J R, Shirvanyants D, Ermoshkin N, Continuous liquid interface production of 3D objects. Science, 2015, 347(6228): 1349–1352

[74] Thomsen C R, Wang F, Sigmund O.Buckling strength topology optimization of 2D periodic materials based on linearized bifurcation analysis, Computer Methods in Applied Mechanics and Engineering, 2018, 339: 115–136

[75] Ahlfors L V. Conformal Invariants: Topics in Geometric Function Theory. Vol. 371. Providence: American Mathematical Society, 2010

[76] Jin M, Kim J, Luo F, Discrete surface Ricci flow. IEEE Transactions on Visualization and Computer Graphics, 2008, 14(5): 1030–1043

[77] Chow B, Luo F. Combinatorial Ricci flows on surfaces. Journal of Differential Geometry, 2003, 63(1): 97–129

[78] Zeng W, Gu X D. Ricci Flow for Shape Analysis and Surface Registration: Theories, Algorithms and Applications. New York: Springer, 2013

[79] Gu X, He Y, Jin M, Manifold splines with a single extraordinary point. Computer Aided Design, 2008, 40(6): 676–690

[80] Jin M, Luo F, Gu X. Computing surface hyperbolic structure and real projective structure. In: Proceedings of the 2006 ACM Symposium on Solid and Physical Modeling. New York: ACM, 2006, 105–116

[81] Gu X D, Zeng W, Luo F, Numerical computation of surface conformal mappings. Computational Methods and Function Theory, 2012, 11(2): 747–787

Acknowledgement

The authors acknowledge the support from the National Science Foundation of the United States (Grant Nos. CMMI1462270 and CMMI1762287), Ford University Research Program (URP), and the start-up fund from the State University of New York at Stony Brook.

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 you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made.

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/.

RIGHTS & PERMISSIONS

2019 The Author(s) 2019. This article is published with open access at link.springer.com and journal.hep.com.cn