Multi-objective concurrent isogeometric topology optimization of multiscale structures

Abstract

Multiscale structures require excellent multiphysical properties to withstand the loads in various complex engineering fields. In this study, a concurrent isogeometric topology optimization method is proposed to design multiscale structures with high thermal conductivity and low mechanical compliance. First, the mathematical description model of multi-objective topology optimization for multiscale structures is constructed, and a single-objective concurrent isogeometric topology optimization formulation for mechanical and thermal compliance is proposed. Then, by combining the isogeometric analysis method, the material interpolation model and decoupled sensitivity analysis scheme of the objective function are established on macro and micro scales. The solid isotropic material with penalization method is employed to update iteratively the macro and microstructure topologies simultaneously. Finally, the feasibility and advantages of the proposed approach are illustrated by several 2D and 3D numerical examples with different volume fractions, while the effects of volume fraction and different boundary conditions on the final configuration and multi-objective performance of the multiscale structure are explored. Results show that the isogeometric concurrent design of multiscale structures through multi-objective optimization can produce better multi-objective performance compared with a single-scale one.

Graphical abstract

Keywords

isogeometric topology optimization / multiscale structure / multi-objective optimization / thermal conductivity / mechanical compliance

1 Introduction

Multiscale structures have different extraordinary physical properties, such as ultralightness and ultrastiffness [1,2], high ionic conductivity [3], heat-transfer capability [4], and thermal insulation property [5], which make multiscale structures exhibited in a wide variety of complex engineering applications, including mechanical engineering, aerospace, optoelectronics, and biomedical engineering [6]. Multi-objective optimization allows the design of multiscale structure with several different properties, including improved thermal and mechanical properties [7,8]. As a research hotspot in the field of structural optimization, topology optimization is a numerical iterative procedure to find the optimal material distribution by following the specified objective functions and constraints in a predefined design domain, which has the advantage of more remarkable material savings and improved structural performance. Currently, the main methods for topology optimization are homogenization method [9], variable density method [10,11], evolutionary structural optimization [12], level set method [13], and moving morphable components [14]. In addition, topology optimization methods based on finite difference method, finite volume method, and finite element method (FEM) can hardly guarantee the smoothness of the discrete model on the curved boundary [15,16]. Integrating the isogeometric analysis (IGA) method [17], which uses nonuniform rational B-splines (NURBS) as basis functions, with topology optimization is expected to solve the above problems.

In recent years, concurrent topology optimization for multiscale structures has been investigated by numerous scholars. Kato et al. [18] assumed that the topology of the macrostructure is invariant, while the topology of material microstructure is unique with respect to the macrostructure; moreover, they optimized the topology of material microstructure to reflect the mechanical behavior of the macrostructure. In contrast to other independent designs of macrostructures or material microstructures, a two-scale topology optimization algorithm, where the design of the macrostructure also changes with the evolution of the material microstructure, is proposed [19,20]. These concurrent topology optimization methods for one material microstructure are simple to formulate and computationally efficient but have relatively limited design freedom. Another concurrent topology optimization approach that considers multiple material microstructures [21,22] effectively expands the design space but is computationally expensive due to the large number of different microstructures that must be optimized simultaneously, thus posing a serious challenge of microstructural connectivity. To address the above issues, Coelho et al. [23] developed a computational procedure for two-scale topology optimization problems using parallel computing techniques, and the procedure demonstrated the relevance of parallel computing in solving such problems. The problem of high computational cost can also be solved by some alternative design formulations [24], such as conformal mapping [25], preparameterized micro unit cells [26], and the projection method [27]. Numerous studies have been conducted to address connectivity in multiscale topology optimization [28–30], which can be broadly classified into three subtypes: morphology deformation, strong control, and weak adjustment [24]. Most of the above concurrent topology optimization studies are based on the mechanical properties to obtain the final designs of macro and microstructure, but the concurrent topology optimization can also be extended to the design of thermal materials [31], electromagnetic materials [32], and so on.

In many aerospace applications, thermal stress from high local temperatures may trigger unexpected failure mechanisms [33]. Furthermore, with the advancement in integrated circuits downscaling of micro and nanoelectronics, the demand is increased for reducing the size and the weight of the heat passive components and sustaining appropriated structural stiffness [34]. The shell in limited space must perform the structural role of supporting the electric part from external loads and the role of a heat sink to release heat generated by the electric device into the open air simultaneously [35]. In above design tasks where the devices have to dissipate heat and withstand an external load, design for minimum compliance does not adequately guarantee structural integrity, and thus, structural thermal performance must be considered. Therefore, multi-objective topology optimization for multiscale structures has a wide range of application scenarios. Yan et al. [36] obtained a stiff but lightweight structure consisting of thermal insulating materials by simultaneously optimizing the topology of the macrostructure and its material microstructure driven by the multi-objective functions, macrostructure stiffness, and material thermal conductivity. The structures with excellent multi-objective properties of mechanical and heat have also been obtained in fields, such as multi-materials and metamaterials [37,38]. For other multi-objective properties, Al Ali and Shimoda [39] used concurrent multiscale multiphysics topology optimization to design lightweight porous composite structures that have resilience toward mechanical and hygral and thermal loads. Bao et al. [40] proposed a two-stage optimization design method to obtain mechanical metamaterials with customized low-frequency bandgap and negative Poisson’s ratio. In most cases, the selection of weighting coefficients among different objectives becomes critical because performance conflicts often occur among multiple objectives, and further improvement of one part of the objective leads to considerable weakening of another part of the objective. Zheng et al. [5] defined a novel multi-objective optimization function in which the weighting coefficients can be adaptively adjusted to obtain topological designs of hierarchical materials with multiple properties in terms of thermal insulation and mechanics. These results prove the excellent performance of multiscale structures and fully demonstrate the necessity of topology optimization under multi-objective. However, all of the above finite element mesh-based topology optimization methods rely on discrete mesh models, where the continuity of the shape function is usually of low order and obtaining high-quality meshes is difficult [41], leading to final results prone to numerical instability and checkerboarding when filtering techniques are not used [16,42].

IGA has been applied in the field of topology optimization with its advantages of high accuracy and continuity, and it has been employed in topology optimization to replace the FEM, in what is termed as isogeometric topology optimization (ITO) [43]. Seo et al. [44] first proposed the use of ITO for structural response analysis and sensitivity calculation in topology optimization by treating topologically complex spline surfaces using trimmed surface analysis. Wang and Benson [45] proposed a geometrically constrained isogeometric level-set based TO, which provides higher accuracy and efficiency than previous methods underpinned on FEM because of the properties of IGA. In addition, the framework of ITO can likewise address multimaterial, multi-objective, and multiscale issues. Gao et al. [41] introduced a NURBS-based multimaterial interpolation model and developed a multimaterial ITO method. Zhang et al. [16] proposed a multi-objective topology optimization model of thermal–mechanical coupling anisotropic structures based on the IGA method, which was verified to be correct and superior to the finite element model. Gai et al. [46] applied IGA-based topology optimization to the material optimization design with extreme shear modulus and negative Poisson’s ratio.

In this study, a new concurrent topology optimization of multiscale structures is proposed using isogeometric mesh as the basis for multi-objective performance considering structural flexibility minimization and thermal flexibility minimization. The energy-based homogenization method [47,48] is used to calculate multiple effective properties of material microstructures and solid isotropic material with penalization (SIMP) method is used as a material interpolation scheme to optimize the macrostructure and material microstructure simultaneously. Compared with previous work, this study introduces IGA into concurrent topology optimization of multiscale structures and combines it with multi-objective optimization to determine the model, formulation, and flow of multi-objective concurrent ITO problems. The results show that the proposed method remarkably improves the thermal conductivity and stiffness of multiscale structure.

The rest of this paper is organized as follows: Section 2 introduces the underlying theory of concurrent ITO methods, mainly the theory of ITO and homogenization theory. Section 3 develops a multi-objective concurrent ITO model, showing the details of its formulation and sensitivity analysis and the algorithmic flow. Section 4 demonstrates several 2D and 3D numerical examples. Finally, Section 5 presents the concluding remarks.

2 Theoretical basis

Concurrent ITO introduces IGA into the concurrent topology optimization framework so that the macro and micro design domains are discretized into two independent IGA grids, based on which the topologies of the macro and micro structures are optimized respectively. Therefore, ITO theory and homogenization theory for concurrent topology optimization are summarized in this section as basic theories. Detailed discussion on ITO and homogenization can be found in Refs. [17,47,49,50].

2.1 Isogeometric topology optimization

ITO is based on NURBS, which are the most common representation of curves and surfaces in computer-aided design and computer graphics. Given a knot vector $\mathbf{\Xi }=\{{\xi }_1,{\xi }_2,\dots ,{\xi }_{m+q+1}\}$, where q is the degree of the spline curve and m is the number of basis functions (also the number of control points). The B-spline basis functions used in the construction of NURBS can be defined recursively via the Cox-de Boor recursion formula:

(1)$\begin{array}{c}{B}_{i,0}(\xi )=\{\begin{array}{ll}1& \text{if}{\xi }_i⩽\xi ⩽{\xi }_{i+1},\\ 0& \text{otherwise,}\end{array}\\ B_{i,q}(\xi )=\frac{\xi -{\xi }_i}{{\xi }_{i+q}-{\xi }_i}{B}_{i,q-1}(\xi )+\frac{{\xi }_{i+q+1}-\xi }{{\xi }_{i+q+1}-{\xi }_{i+1}}{B}_{i+1,q-1}(\xi ),\end{array}$

where the convention 0/0 = 0 is adpoted. By introducing weights ${\omega }_i$ to the B-spline basis functions, the NURBS basis function is obtained as

(2)$N_{i,q}(\xi )=\frac{B_{i,q}(\xi ){\omega }_i}{{\displaystyle \sum _{j=1}^m{B}_{j,q}(\xi ){\omega }_j}}.$

Through the tensor product property, Eq. (2) can be further extended to two-dimensional or high-dimensional NURBS basis functions.

The numerical computation of a given physical field in ITO is computed at the control points so that the variable $X$ (e.g., coordinate, force, or displacement) at the position of $\xi$ can be obtained from the control point values:

(3)$X\left(\xi \right)=\sum _j{N}_j\left(\xi \right)X_j,$

where $N_j$ is the basis function for the jth control point, and $X_j$ is the value corresponding to the control point. When using IGA, replacing the design variable from the element density to the control point density is advantageous [50]. The element density ${\rho }_i^{\mathrm{e}}$ is expressed as the density ${\rho }_i^c$ at the center of the element, which can be calculated by the following equation:

(4)${\rho }_i^{\mathrm{e}}={\rho }_i^{\mathrm{c}}=\sum _{j\in {\mathrm{S}}_{\mathrm{c}}}^{}{N}_{ij}^{\mathrm{c}}{\rho }_{ij},$

where $\mathrm{S}{}_{\mathrm{c}}$ denotes the set of control points affecting the ith element, $N_{ij}^{\mathrm{c}}$ is the basis function of the jth control point corresponding to the center of the ith element, and ${\rho }_{ij}$ denotes the density of the jth control point. The sensitivity of the objective function $c$ to ${\rho }_k$ at any control point can be expressed as a 1st-order derivative of the density of the elements that the objective function affects at that control point:

(5)$\frac{\mathrm{\partial }c}{\mathrm{\partial }{\rho }_k}=\sum _{h\in {\mathrm{S}}_{\mathrm{e}}}^{}\frac{\mathrm{\partial }c}{\mathrm{\partial }{\rho }_{hk}^{\mathrm{e}}}\frac{\mathrm{\partial }{\rho }_{hk}^{\mathrm{e}}}{\mathrm{\partial }{\rho }_k},$

where ${\mathrm{S}}_{\mathrm{e}}$ represents the set of elements affected by the kth control point. ${\rho }_{hk}^{\mathrm{e}}$ then represents the density of the hth element in ${\mathrm{S}}_{\mathrm{e}}$.

According to Eq. (4), the 1st-order derivative of the density of any element with respect to the density of the control point of the influencing element is

(6)$\frac{\mathrm{\partial }{\rho }_i^{\mathrm{e}}}{\mathrm{\partial }{\rho }_{ij}}=N_{ij}^{\mathrm{c}}.$

On the basis of the above equation, the $\mathrm{\partial }{\rho }_{hk}^{\mathrm{e}}/\mathrm{\partial }{\rho }_k$ in Eq. (5) can be obtained. Given that ITO is based on control points, NURBS filters are preferred over the conventional distance-based filters in this study. In Eq. (5), the NURBS basis functions are used as weighting functions, i.e., the sensitivity values are calculated by the spatial average of the neighboring NURBS basis functions, where the influence region is a rectangular region, as shown in Fig.1(b). The filtering region can be adjusted by increasing/decreasing the NURBS degree of the filter, and the specific details about NURBS filtering can be found in Ref. [50].

Fig.1 Illustration of control points and elements controlling and influencing each other. (a) Element i controlled by a set of control points. (b) Control point k and the set of elements affected by it.

Full size|PPT slide

2.2 Homogenization theory

The equivalent properties of microstructures in the linear elastic range can be assessed by homogenization methods. As shown in Fig.2, the microstructures are periodically arranged in the macro design domain and can be considered as homogenized materials in the macro scale. The displacement field at the macrostructure can be expressed in terms of asymptotic expansion theory:

Fig.2 Homogenization for multiscale structures.

Full size|PPT slide

(7)${\mathit{u}}^{\zeta }\left(\mathit{x}\right)={\mathit{u}}_0\left(x,y\right)+\zeta {\mathit{u}}_1\left(x,y\right)+{\zeta }^2{\mathit{u}}_2\left(x,y\right)+\cdots ,$

where $\zeta =x/y\text{(}0<\zeta <<1)$ is the scale factor, i.e., the aspect ratio between the macro and micro scale. When only the 1st-order terms of the asymptotic expansion in Eq. (7) are considered, the homogenized elastic tensor ${\mathit{E}}^{\mathrm{h}\mathrm{o}\mathrm{m}}$ can be obtained by integrating over the periodic unit cells (PUC), according to the Einstein summation convention, given by

(8)$E_{ijmn}^{\mathrm{h}\mathrm{o}\mathrm{m}}=\frac{1}{\left|{\mathrm{\Omega }}_{\mathrm{m}}\right|}{\int }_{{\mathrm{\Omega }}_{\mathrm{m}}}^{}{E}_{ijpq}({\epsilon }_{pq}^{0\left(mn\right)}-{\epsilon }_{pq}^{\ast \left(mn\right)})\mathrm{d}{\mathrm{\Omega }}_{\mathrm{m}},$

where $E_{ijmn}^{\mathrm{hom}}$ is the homogenized elastic tensor in index notation $ijmn$; $\left|{\mathrm{\Omega }}_{\mathrm{m}}\right|$ is the area (2D) or volume (3D) of the periodic cells; $E_{ijpq}$ is the locally varying elastic tensor, which is equal to the elasticity matrix of the base material in the solid region and is approximated as a zero matrix in the blank region. In the 2D case, ${\epsilon }_{pq}^{0\left(mn\right)}$are the prescribed unit strain fields corresponding to three (mn = 11,22,12) linearly independent unit test strain fields, i.e., unit strain in horizontal direction ${\epsilon }_1$, unit strain in vertical direction ${\epsilon }_2$, and unit shear strain ${\epsilon }_3$. Fig.3(a) illustrates the components of the strain field vectors correspond to pq = 11,22,12. Two unit test temperature gradients $g_1$ and $g_2$ are imposed on the PUC with respect to thermal conductivity problem, as shown in Fig.3(b). For the 3D case in Section 4.3, the unit test strain contains six ($mn=11,22,33,12,23,13$) linearly independent unit strain fields corresponding to ${\left(1,0,0,0,0,0\right)}^{\mathrm{T}}$, ${\left(0,1,0,0,0,0\right)}^{\mathrm{T}}$, ${\left(0,0,1,0,0,0\right)}^{\mathrm{T}}$, ${\left(0,0,0,1,0,0\right)}^{\mathrm{T}}$, ${\left(0,0,0,0,1,0\right)}^{\mathrm{T}}$, and ${\left(0,0,0,0,0,1\right)}^{\mathrm{T}}$, while three unit test temperature gradients are imposed on the PUC in the case of the thermal conductivity problem. The periodic characteristic strain ${\epsilon }_{pq}^{\ast \left(mn\right)}$ is obtained by solving the following equilibrium equation [51]:

Fig.3 Two kinds of test fields imposed on the boundaries of isogeometric PUC. (a) Three unit test strain fields. (b) Two unit test temperature gradients.

Full size|PPT slide

(9)${\int }_{{\mathrm{\Omega }}_{\mathrm{m}}}^{}{E}_{ijpq}{\epsilon }_{pq}^{\ast \left(mn\right)}\frac{\mathrm{\partial }{v}_i}{\mathrm{\partial }{y}_i}\mathrm{d}{\mathrm{\Omega }}_{\mathrm{m}}={\int }_{{\mathrm{\Omega }}_{\mathrm{m}}}^{}{E}_{ijpq}{\epsilon }_{pq}^{0\left(mn\right)}\frac{\mathrm{\partial }{v}_i}{\mathrm{\partial }{y}_i}\mathrm{d}{\mathrm{\Omega }}_{\mathrm{m}},$

where $v_i$ is periodic admissible displacement field in the micro scale.

Efficient existing algorithms used in topology optimization are supported by rewriting Eq. (8) as an equivalent form of elemental mutual energies [47,51].

(10)$E_{ijmn}^{\mathrm{h}\mathrm{o}\mathrm{m}}=\frac{1}{\left|{\mathrm{\Omega }}_{\mathrm{m}}\right|}{\int }_{{\mathrm{\Omega }}_{\mathrm{m}}}^{}{E}_{klpq}\left({\epsilon }_{kl}^{0\left(ij\right)}-{\epsilon }_{kl}^{\ast \left(ij\right)}\right)\left({\epsilon }_{pq}^{0\left(mn\right)}-{\epsilon }_{pq}^{\ast \left(mn\right)}\right)\mathrm{d}{\mathrm{\Omega }}_{\mathrm{m}}.$

The homogenized elastic tensor in Eq. (10) can be approximated if the PUC is divided into n isogeometric elements:

(11)$E_{ijmn}^{\mathrm{h}\mathrm{o}\mathrm{m}}=\frac{1}{\left|{\mathrm{\Omega }}_{\mathrm{m}}\right|}\sum _{e=1}^n{\left({\mathit{u}}_e^{0\left(ij\right)}-{\mathit{u}}_e^{\ast \left(ij\right)}\right)}^{\mathrm{T}}{\mathit{k}}_e\left({\mathit{u}}_e^{0\left(mn\right)}-{\mathit{u}}_e^{\ast \left(mn\right)}\right),$

where ${\mathit{u}}_e^{\ast \left(mn\right)}$ and ${\mathit{u}}_e^{\ast \left(ij\right)}$ are control point unknown displacements of eth element, and ${\mathit{k}}_e$ is the element stiffness matrix. In the energy-based isogeometric homogenization method (EBIGHM), the unit test strains are directly imposed to the boundaries of PUC, inducing ${\epsilon }_{kl}^{A\left(ij\right)}$ corresponding to the superimposed strain fields [47] $\left({\epsilon }_{kl}^{0\left(ij\right)}-{\epsilon }_{kl}^{\ast \left(ij\right)}\right)$ in Eq. (10); thus, Eq. (11) can be written in the following equivalent form:

(12)$\begin{array}{rl}{E}_{ijmn}^{\mathrm{h}\mathrm{o}\mathrm{m}}& =\frac{1}{\left|{\mathrm{\Omega }}_{\mathrm{m}}\right|}\sum _{e=1}^n{Q}_{ijmn}^e\\ & =\frac{1}{\left|{\mathrm{\Omega }}_{\mathrm{m}}\right|}\sum _{e=1}^N{\left({\mathit{u}}_e^{A\left(ij\right)}\right)}^{\mathrm{T}}{\mathit{k}}_e{\mathit{u}}_e^{A\left(mn\right)},\end{array}$

where ${\mathit{u}}_e^{A\left(ij\right)}$ denotes the element’s control point displacement solution corresponding to ${\epsilon }_{kl}^{A\left(ij\right)}$, and $Q_{ijmn}^e$ denotes element mutual energy. Under the assumption of periodicity [47], the difference between the displacements on a pair of opposite parallel boundaries of PUC shown in Fig.4(a) can be expressed as

Fig.4 Periodic boundary condition of PUC and its mesh in 2D case. (a) PUC; (b) bilinear finite elements; (c) biquadratic isogeometric elements.

Full size|PPT slide

(13)${\omega }_i^k=u_i^{k+}-u_i^{k-}={\epsilon }_{ij}^0\left(y_j^{k+}-y_j^{k-}\right)={\epsilon }_{ij}^0\mathrm{\Delta }{y}_j^k,$

where $y_j^{k+}$ and $y_j^{k-}$ represent the coordinates of the two opposite boundaries “k+” and “k−” of the PUC, respectively; the subscript j represents the coordinate direction.

For any given parallel hexagonal PUC, $\mathrm{\Delta }{y}_j^k$ is constant (note that in Fig.4(a), we have $\mathrm{\Delta }{y}_1^1=y_1^0$, $\mathrm{\Delta }{y}_2^2=y_2^0$, $\mathrm{\Delta }{y}_2^1=0$, and $\mathrm{\Delta }{y}_1^2=0$). Thus, with a specified ${\epsilon }_{ij}^0$, ${\omega }_i^k={\epsilon }_{ij}^0\mathrm{\Delta }{y}_j^k$ becomes a constant.

This form of boundary condition converts the unit test strain into prescribed displacement applied to the PUC boundary. After dividing the PUC into a finite element mesh (Fig.4(b)) or isogeometric mesh (Fig.4(c)), the boundary condition can be applied directly to the mesh model by constraining the corresponding pairs of nodal or control point displacements. The boundary control points of the parallel hexagonal isogeometric model are located exactly on the boundary, which is the same as the boundary nodes for finite element model. Thus, even though the basis functions in IGA are not interpolating functions [52], the boundary conditions can still be enforced directly as in the finite element method, without leading to major errors. The study [46] shows that the elasticity matrices obtained from isogeometric homogenization are almost identical compared with finite element homogenization. For the PUC shown in Fig.4(c), the equilibrium equations for IGA are simplified due to the periodic boundary conditions and the specified boundary displacements.

(14)$\begin{array}{rl}& \left[\begin{array}{cc}{\mathit{K}}_{22}& {\mathit{K}}_{23}+{\mathit{K}}_{24}\\ \mathrm{S}\mathrm{y}\mathrm{m}& {\mathit{K}}_{33}+{\mathit{K}}_{34}+{\mathit{K}}_{43}+{\mathit{K}}_{44}\end{array}\right]\left[\begin{array}{c}{\mathit{U}}_2\\ {\mathit{U}}_3\end{array}\right]\\ =& -\left[\begin{array}{c}{\mathit{K}}_{21}\\ {\mathit{K}}_{31}+{\mathit{K}}_{41}\end{array}\right]{\overline{\mathit{U}}}_1-\left[\begin{array}{c}{\mathit{K}}_{24}\\ {\mathit{K}}_{34}+{\mathit{K}}_{44}\end{array}\right]\overline{\mathit{W}},\end{array}$

where the global displacement vector is divided into four parts: ${\overline{\mathit{U}}}_1$ is the prescribed displacements of the four corner control points (blue control points in Fig.4(c)); ${\mathit{U}}_2$ is the unknown displacements corresponding to the interior control points (red control points); ${\mathit{U}}_3$ and ${\mathit{U}}_4$ denote the displacements of the control points on the two pairs of parallel boundaries (corresponding to the green and yellow control points, respectively), which should satisfy ${\mathit{U}}_4={\mathit{U}}_3+\overline{\mathit{W}}$. $\overline{\mathit{W}}$ is the prescribed value computed in accordance with Eq. (13). The details of the periodic boundary conditions for the 3D case can be found in Ref. [53].

In the 2D case, the homogenized elastic tensor can be expanded as

(15)${\mathit{E}}^{\mathrm{hom}}=\left[\begin{array}{ccc}{E}_{1111}^{\mathrm{hom}}& E_{1122}^{\mathrm{hom}}& E_{1112}^{\mathrm{hom}}\\ E_{2211}^{\mathrm{hom}}& E_{2222}^{\mathrm{hom}}& E_{2212}^{\mathrm{hom}}\\ E_{1211}^{\mathrm{hom}}& E_{1222}^{\mathrm{hom}}& E_{1212}^{\mathrm{hom}}\end{array}\right].$

The homogenized thermal conductivity tensor can be given similarly as Eq. (8):

(16)${\mathrm{\Lambda }}_{ijmn}^{\mathrm{h}\mathrm{o}\mathrm{m}}=\frac{1}{\left|{\mathrm{\Omega }}_{\mathrm{m}}\right|}{\int }_{{\mathrm{\Omega }}_{\mathrm{m}}}^{}{\mathrm{\Lambda }}_{ijpq}(g_{pq}^{0\left(mn\right)}-g_{pq}^{\ast \left(mn\right)})\mathrm{d}{\mathrm{\Omega }}_{\mathrm{m}},$

where $g_{pq}^{0\left(mn\right)}$ denotes the unit test temperature gradients; in contrast to the unit test strain fields, the 2D case has two ($mn=11,22$) unit test temperature gradients (Fig.3(b)). $g_{pq}^{\ast \left(mn\right)}$ denotes the periodic characteristic gradient. Furthermore, in EBIGHM, the PUC is divided into $n$ elements, and the unit test temperature gradients are applied directly on the boundary of the basic cell, Eq. (16) is then rewritten as

(17)${\mathrm{\Lambda }}_{ijmn}^{\mathrm{h}\mathrm{o}\mathrm{m}}=\frac{1}{\left|{\mathrm{\Omega }}_{\mathrm{m}}\right|}\sum _{e=1}^n{\left({\mathit{t}}_e^{A\left(ij\right)}\right)}^{\mathrm{T}}{\mathit{k}}_e^{\mathrm{t}\mathrm{h}}{\mathit{t}}_e^{A\left(mn\right)},$

where ${\mathit{t}}_e^{A\left(ij\right)}$ is the control point temperature, and ${\mathit{k}}_e^{\mathrm{t}\mathrm{h}}$ is the element thermal conductivity matrix.

ln the 2D case, the homogenized thermal conductivity tensor ${\mathbf{\Lambda }}^{\mathrm{hom}}$ can be expanded as

(18)${\mathbf{\Lambda }}^{\mathrm{h}\mathrm{o}\mathrm{m}}=\left[\begin{array}{cc}{\mathrm{\Lambda }}_{1111}^{\mathrm{h}\mathrm{o}\mathrm{m}}& {\mathrm{\Lambda }}_{1122}^{\mathrm{h}\mathrm{o}\mathrm{m}}\\ {\mathrm{\Lambda }}_{2211}^{\mathrm{h}\mathrm{o}\mathrm{m}}& {\mathrm{\Lambda }}_{2222}^{\mathrm{h}\mathrm{o}\mathrm{m}}\end{array}\right].$

3 Multi-objective concurrent isogeometric topology optimization

In this section, a multi-objective concurrent ITO model is developed. As shown in Fig.5., macro and micro topologies are considered in the concurrent isogeometric topology design of multiscale structures through multi-objective optimization to combine structural performance and thermal conductivity.

Fig.5 Concurrent isogeometric topology optimization for multiscale structures.

Full size|PPT slide

3.1 Multi-objective optimization

The multi-objective optimization model is formulated in accordance with the weighting factor w using the linear weighted sum method to obtain the target multiscale structures with minimum mechanical compliance (or maximum mechanical stiffness) and minimum thermal compliance (or maximum thermal conductivity):

(19)$\{\begin{array}{c}\text{Find:}{\mathit{\rho }}_{\mathrm{M}}=\left[{\rho }_{\mathrm{M}}^1,{\rho }_{\mathrm{M}}^2,\dots ,{\rho }_{\mathrm{M}}^i,\dots \right],\\ {\mathit{\rho }}_{\mathrm{m}}=\left[{\rho }_{\mathrm{m}}^1,{\rho }_{\mathrm{m}}^2,\dots ,{\rho }_{\mathrm{m}}^j,\dots \right].\\ \text{Min:}C\left({\mathit{\rho }}_{\mathrm{M}},{\mathit{\rho }}_{\mathrm{m}}\right)=wC{}_{\mathrm{M}\mathrm{e}\mathrm{c}\mathrm{h}}+\left(1-w\right)C{}_{\mathrm{H}\mathrm{e}\mathrm{a}\mathrm{t}}.\\ \mathrm{s}.\mathrm{t}.\{\begin{array}{c}\mathit{K}\left({\mathit{\rho }}_{\mathrm{M}},{\mathit{\rho }}_{\mathrm{m}}\right){\mathit{U}}_{\mathrm{M}}=\mathit{F},\\ {\mathit{K}}^{\mathrm{t}\mathrm{h}}\left({\mathit{\rho }}_{\mathrm{M}},{\mathit{\rho }}_{\mathrm{m}}\right){\mathit{T}}_{\mathrm{M}}=\mathit{P},\end{array}\\ {\int }_{{\mathrm{\Omega }}_{\mathrm{M}}}^{}{\mathit{\rho }}_{\mathrm{M}}\mathrm{d}{\mathrm{\Omega }}_{\mathrm{M}}⩽V_{\mathrm{M}},0<{\rho }_{\mathrm{M}}^{min}<{\rho }_{\mathrm{M}}^i⩽1,\\ {\int }_{{\mathrm{\Omega }}_{\mathrm{m}}}^{}{\mathit{\rho }}_{\mathrm{m}}\mathrm{d}{\mathrm{\Omega }}_{\mathrm{m}}⩽V_{\mathrm{m}},0<{\rho }_{\mathrm{m}}^{min}<{\rho }_{\mathrm{m}}^j⩽1,\end{array}$

where $\mathit{K}$ and ${\mathit{K}}^{\mathrm{t}\mathrm{h}}$ are the stiffness matrix and thermal conductivity matrix of the global performance design domain. $\mathit{F}$ and ${\mathit{U}}_{\mathrm{M}}$ are the macrostructure control point force load and control point displacement, respectively, while $\mathit{P}$ and ${\mathit{T}}_{\mathrm{M}}$ are the macrostructure control point thermal load and control point temperature, respectively. The macrostructure ${\mathrm{\Omega }}_{\mathrm{M}}$ is discretized into $N$ isogeometric elements, while the microstructure $\mathrm{\Omega }{}_{\mathrm{m}}$ is discretized into $n$ isogeometric elements. The design variables (i.e., control point densities) at the macro and microscale are denoted by the ${\mathit{\rho }}_{\mathrm{M}}$ and ${\mathit{\rho }}_{\mathrm{m}}$, respectively; $V_{\mathrm{M}}$ and $V_{\mathrm{m}}$ are the prescribed structural volumes of the ${\mathit{\rho }}_{\mathrm{M}}$ and $\mathit{\rho }{}_{\mathrm{m}}$, respectively, i.e., the maximal volume fractions of macro and microstructure. ${\rho }_{\mathrm{M}}^{min}$ and ${\rho }_{\mathrm{m}}^{min}$ are the minimum values of macro and micro variables, respectively, to avoid numerical singularity in the optimization. $C$ denotes the multi-objective performance of multiscale structure. $C_{\mathrm{M}\mathrm{e}\mathrm{c}\mathrm{h}}$ and $C{}_{\mathrm{H}\mathrm{e}\mathrm{a}\mathrm{t}}$ are the single-objective performance of mechanical and heat, i.e., mechanical compliance and thermal compliance, respectively.

$C_{\mathrm{M}\mathrm{e}\mathrm{c}\mathrm{h}}$ can be expressed as

(20)$\begin{array}{rl}{C}_{\mathrm{M}\mathrm{e}\mathrm{c}\mathrm{h}}\left({\mathit{\rho }}_{\mathrm{M}},{\mathit{\rho }}_{\mathrm{m}}\right)& ={\mathit{U}}_{\mathrm{M}}^{\mathrm{T}}\mathit{K}\left({\mathit{\rho }}_{\mathrm{M}},{\mathit{\rho }}_{\mathrm{m}}\right)\mathit{U}{}_{\mathrm{M}}\\ & =\iota \left({\mathit{\rho }}_{\mathrm{M}}^{\mathrm{e}}\right)\sum _{e=1}^N{\mathit{U}}_e^{\mathrm{T}}{\mathit{K}}_e\left({\mathit{\rho }}_{\mathrm{m}}^{\mathrm{e}}\right){\mathit{U}}_e,\end{array}$

where $\iota \left({\mathit{\rho }}_{\mathrm{M}}^{\mathrm{e}}\right)$ donates the material interpolation function of macrostructure; ${\mathit{\rho }}_{\mathrm{M}}^{\mathrm{e}}$ is the isogeometric element densities of macrostructure, which can be derived by interpolating the control point densities ${\mathit{\rho }}_{\mathrm{M}}$ according to Eq. (4); ${\mathit{U}}_e$ is the eth macro element displacement vector; ${\mathit{K}}_e\left({\mathit{\rho }}_{\mathrm{m}}^e\right)={\int }_{{\mathrm{\Omega }}_{\mathrm{e}}}^{}{\mathit{B}}^{\mathrm{T}}{\mathit{E}}^{\mathrm{h}\mathrm{o}\mathrm{m}}\mathit{B}\mathrm{d}{\mathrm{\Omega }}_{\mathrm{e}}$ denotes the macro element stiffness matrix, and $\mathrm{\Omega }{}_{\mathrm{e}}$ is the macro element; $\mathit{B}$ is the strain–displacement matrix of macrostructure; ${\mathit{\rho }}_{\mathrm{m}}^{\mathrm{e}}$ is the isogeometric element densities of microstructure, which can be derived by interpolating the control point densities $\mathit{\rho }{}_{\mathrm{m}}$ according to Eq. (4).

In the same way, $C{}_{\mathrm{H}\mathrm{e}\mathrm{a}\mathrm{t}}$ can be expressed as

(21)$\begin{array}{rl}{C}_{\mathrm{H}\mathrm{e}\mathrm{a}\mathrm{t}}\left({\mathit{\rho }}_{\mathrm{M}},{\mathit{\rho }}_{\mathrm{m}}\right)& ={\mathit{T}}_{\mathrm{M}}^{\mathrm{T}}{\mathit{K}}^{\mathrm{t}\mathrm{h}}\left({\mathit{\rho }}_{\mathrm{M}},{\mathit{\rho }}_{\mathrm{m}}\right){\mathit{T}}_{\mathrm{M}}\\ & =\iota \left({\mathit{\rho }}_{\mathrm{M}}^{\mathrm{e}}\right)\sum _{e=1}^N{\mathit{T}}_e^{\mathrm{T}}{\mathit{K}}_e^{\mathrm{t}\mathrm{h}}\left({\mathit{\rho }}_{\mathrm{m}}^{\mathrm{e}}\right){\mathit{T}}_e,\end{array}$

where ${\mathit{T}}_e$ is the eth macro element temperature vector, ${\mathit{K}}_e^{\mathrm{t}\mathrm{h}}\left({\mathit{\rho }}_{\mathrm{m}}^{\mathrm{e}}\right)={\int }_{{\mathrm{\Omega }}_{\mathrm{e}}}^{}{\mathit{D}}^{\mathrm{T}}{\mathbf{\Lambda }}^{\mathrm{h}\mathrm{o}\mathrm{m}}\mathit{D}\mathrm{d}{\mathrm{\Omega }}_{\mathrm{e}}$ denotes the macro element thermal conductivity matrix, and $\mathit{D}$ is the conversion matrix of the temperature gradient to the control point temperature on the macrostructure.

3.2 Material interpolation scheme

The elastic tensors of the macro and microstructure are defined by the material interpolation scheme in the modified SIMP method:

(22)$\{\begin{array}{c}{\mathit{E}}_{\mathrm{M}}=\iota \left({\mathit{\rho }}_{\mathrm{M}}^{\mathrm{e}}\right){\mathit{E}}^{\mathrm{h}\mathrm{o}\mathrm{m}}=\left[\mathrm{\Delta }+{\left({\mathit{\rho }}_{\mathrm{M}}^{\mathrm{e}}\right)}^p\left(1-\mathrm{\Delta }\right)\right]{\mathit{E}}^{\mathrm{h}\mathrm{o}\mathrm{m}},\\ {\mathit{E}}_{\mathrm{m}}=\iota \left({\mathit{\rho }}_{\mathrm{m}}^{\mathrm{e}}\right){\mathit{E}}_0=\left[\mathrm{\Delta }+{\left({\mathit{\rho }}_{\mathrm{m}}^{\mathrm{e}}\right)}^p\left(1-\mathrm{\Delta }\right)\right]{\mathit{E}}_0,\end{array}$

where ${\mathit{E}}_{\mathrm{M}}$ and ${\mathit{E}}_{\mathrm{m}}$ are the elastic tensors of macro and microstructure, respectively; ${\mathit{E}}_0$ denotes the constitutive elastic tensor of the material; $\iota \left({\mathit{\rho }}_{\mathrm{m}}^e\right)$ denotes the material interpolation function of microstructure; Δ = 1E−9 is a constant to avoid singularities in the stiffness matrix; $p$ is the penalty parameter. As shown in Fig.6, the homogenized elastic tensor ${\mathit{E}}^{\mathrm{h}\mathrm{o}\mathrm{m}}$ is the bridge that connects the macro and microstructure, and it participates in the interpolation as the effective elastic tensor of the macro elements in Eq. (22), which in turn is calculated on the basis of the topology of the microstructure as follows:

Fig.6 Micro and macro coupling for concurrent topology optimization.

Full size|PPT slide

(23)$\begin{array}{rl}{\mathit{E}}^{\mathrm{h}\mathrm{o}\mathrm{m}}& ={\displaystyle \frac{1}{\left|{\mathrm{\Omega }}_{\mathrm{m}}\right|}}{\int }_{{\mathrm{\Omega }}_{\mathrm{m}}}^{}{\mathit{E}}_{\mathrm{m}}\left({\mathit{\epsilon }}^0-{\mathit{\epsilon }}^{\ast }\right)\mathrm{d}{\mathrm{\Omega }}_{\mathrm{m}}\\ & ={\displaystyle \frac{\iota \left({\mathit{\rho }}_{\mathrm{m}}^{\mathrm{e}}\right)}{\left|{\mathrm{\Omega }}_{\mathrm{m}}\right|}}{\int }_{{\mathrm{\Omega }}_{\mathrm{m}}}^{}{\mathit{E}}_0\left({\mathit{\epsilon }}^0-{\mathit{\epsilon }}^{\ast }\right)\mathrm{d}{\mathrm{\Omega }}_{\mathrm{m}}\\ & ={\displaystyle \frac{\iota \left({\mathit{\rho }}_{\mathrm{m}}^{\mathrm{e}}\right)}{\left|{\mathrm{\Omega }}_{\mathrm{m}}\right|}}\sum _{e=1}^n{\left({\mathit{u}}_e^{A\left(ij\right)}\right)}^{\mathrm{T}}{\mathit{k}}_0{\mathit{u}}_e^{A\left(mn\right)},\\ & \end{array}$

where ${\mathit{\epsilon }}^0-{\mathit{\epsilon }}^{\ast }$ denotes the superimposed strain field, the detailed realization of which has been specifically given in Section 2.2; ${\mathit{k}}_0={\int }_{{\mathrm{\Omega }}_0}^{}{\mathit{b}}^{\mathrm{T}}{\mathit{E}}_0\mathit{b}\mathrm{d}{\mathrm{\Omega }}_0$ denotes the micro element stiffness matrix, and ${\mathrm{\Omega }}_0$ is the micro element. $\mathit{b}$ is the strain–displacement matrix of microstructure.

Similarly, the thermal conductivity tensors of the macro and microstructure are defined by the material interpolation scheme in the modified SIMP method:

(24)$\{\begin{array}{c}{\mathbf{\Lambda }}_{\mathrm{M}}=\iota \left({\mathit{\rho }}_{\mathrm{M}}^{\mathrm{e}}\right){\mathbf{\Lambda }}^{\mathrm{h}\mathrm{o}\mathrm{m}},\\ {\mathbf{\Lambda }}_{\mathrm{m}}=\iota \left({\mathit{\rho }}_{\mathrm{m}}^{\mathrm{e}}\right){\mathbf{\Lambda }}_0,\end{array}$

where ${\mathbf{\Lambda }}_{\mathrm{M}}$ and ${\mathbf{\Lambda }}_{\mathrm{m}}$ are the thermal conductivity tensor of macro and microstructure, respectively; ${\mathbf{\Lambda }}_0$ denotes the intrinsic thermal conductivity tensor of the material. The homogenized thermal conductivity tensor ${\mathbf{\Lambda }}^{\mathrm{h}\mathrm{o}\mathrm{m}}$ is calculated on the basis of the topology of the microstructure as follows:

(25)$\begin{array}{rl}{\mathbf{\Lambda }}^{\mathrm{h}\mathrm{o}\mathrm{m}}& ={\displaystyle \frac{1}{\left|{\mathrm{\Omega }}_{\mathrm{m}}\right|}}{\int }_{{\mathrm{\Omega }}_{\mathrm{m}}}^{}{\mathbf{\Lambda }}_{\mathrm{m}}\left({\mathit{g}}^0-{\mathit{g}}^{\ast }\right)\mathrm{d}{\mathrm{\Omega }}_{\mathrm{m}}\\ & ={\displaystyle \frac{\iota \left({\mathit{\rho }}_{\mathrm{m}}^{\mathrm{e}}\right)}{\left|{\mathrm{\Omega }}_{\mathrm{m}}\right|}}\sum _{e=1}^n{\left({\mathit{t}}_e^{A\left(ij\right)}\right)}^{\mathrm{T}}{\mathit{k}}_0^{\mathrm{t}\mathrm{h}}{\mathit{t}}_e^{A\left(mn\right)},\\ & \end{array}$

where ${\mathit{g}}^0-{\mathit{g}}^{\ast }$ denotes the superimposed temperature gradient, ${\mathit{k}}_0^{\mathrm{t}\mathrm{h}}={\int }_{{\mathrm{\Omega }}_0}^{}{\mathit{d}}^{\mathrm{T}}{\mathbf{\Lambda }}_0\mathit{d}\mathrm{d}{\mathrm{\Omega }}_0$ denotes the micro element thermal conductivity matrix, and $\mathit{d}$ is the conversion matrix of the temperature gradient to the control point temperature on the microstructure.

3.3 Concurrent isogeometric topology optimization

The concurrent ITO formulation can be solved by many well-established gradient-based algorithms, and thus, the sensitivity information associated with the objective function must be obtained to perform the computation iteratively [48]. For the objective function shown in Eq. (19), the 1st-order derivative of the objective function is expressed as follows because the weighting factor w is constant with respect to the objective function:

(26)$\dot{C}\left({\mathit{\rho }}_{\mathrm{M}},{\mathit{\rho }}_{\mathrm{m}}\right)=w\dot{C}{}_{\mathrm{M}\mathrm{e}\mathrm{c}\mathrm{h}}+\left(1-w\right)\dot{C}{}_{\mathrm{H}\mathrm{e}\mathrm{a}\mathrm{t}}.$

The mechanical compliance sensitivity ${\dot{C}}_{\mathrm{M}\mathrm{e}\mathrm{c}\mathrm{h}}$ in Eq. (26) has two parts: the first one is about the sensitivity of the macro design variable ${\mathit{\rho }}_{\mathrm{M}}$, which is used to update the macro design domain:

(27)$\frac{\mathrm{\partial }{C}_{\mathrm{M}\mathrm{e}\mathrm{c}\mathrm{h}}}{\mathrm{\partial }{\mathit{\rho }}_{\mathrm{M}}}=\sum _{i\in {\mathrm{S}}_{\mathrm{M}}}^{}\frac{\mathrm{\partial }{C}_{\mathrm{M}\mathrm{e}\mathrm{c}\mathrm{h}}}{\mathrm{\partial }{\mathit{\rho }}_{\mathrm{M}i}^{\mathrm{e}}}\frac{\mathrm{\partial }{\mathit{\rho }}_{\mathrm{M}i}^{\mathrm{e}}}{\mathrm{\partial }{\mathit{\rho }}_{\mathrm{M}}},$

where ${\mathrm{S}}_{\mathrm{M}}$ denotes the set of elements affected by the control points corresponding to ${\mathit{\rho }}_{\mathrm{M}}$, ${\mathit{\rho }}_{\mathrm{M}i}^{\mathrm{e}}$ denotes the ith element density in ${\mathrm{S}}_{\mathrm{M}}$, ${\mathit{K}}_e\left({\mathit{\rho }}_{\mathrm{m}}^{\mathrm{e}}\right)$ is related to the micro element density but not to the macro element density; thus, the 1st-order derivative of mechanical compliance with respect to the element density is

(28)$\frac{\mathrm{\partial }{C}_{\mathrm{M}\mathrm{e}\mathrm{c}\mathrm{h}}}{\mathrm{\partial }{\mathit{\rho }}_{\mathrm{M}i}^{\mathrm{e}}}=-p{\left({\mathit{\rho }}_{\mathrm{M}i}^{\mathrm{e}}\right)}^{p-1}\left(1-\mathrm{\Delta }\right)\sum _{e=1}^N{\mathit{U}}_e^{\mathrm{T}}{\mathit{K}}_e\left({\mathit{\rho }}_{\mathrm{m}}^{\mathrm{e}}\right)U_e.$

Substituting Eq. (28) into Eq. (27) yields

(29)$\frac{\mathrm{\partial }{C}_{\mathrm{M}\mathrm{e}\mathrm{c}\mathrm{h}}}{\mathrm{\partial }{\mathit{\rho }}_{\mathrm{M}}}=-\sum _{i\in {\mathrm{S}}_{\mathrm{M}}}^{}\left[p{\left({\mathit{\rho }}_{\mathrm{M}i}^{\mathrm{e}}\right)}^{p-1}\left(1-\mathrm{\Delta }\right)\sum _{e=1}^N{\mathit{U}}_e^{\mathrm{T}}{\mathit{K}}_e\left({\mathit{\rho }}_{\mathrm{m}}^{\mathrm{e}}\right){\mathit{U}}_e\right]\frac{\mathrm{\partial }{\mathit{\rho }}_{\mathrm{M}i}^{\mathrm{e}}}{\mathrm{\partial }{\mathit{\rho }}_{\mathrm{M}}}.$

The second part deals with the sensitivity of the micro design variable ${\mathit{\rho }}_{\mathrm{m}}$:

(30)$\begin{array}{rl}\frac{\mathrm{\partial }{C}_{\mathrm{M}\mathrm{e}\mathrm{c}\mathrm{h}}}{\mathrm{\partial }\mathit{\rho }{}_{\mathrm{m}}}& =-\iota \left({\mathit{\rho }}_{\mathrm{M}}^{\mathrm{e}}\right)\sum _{e=1}^N{\mathit{U}}_e^{\mathrm{T}}\frac{\mathrm{\partial }{\mathit{K}}_e\left({\mathit{\rho }}_{\mathrm{m}}^{\mathrm{e}}\right)}{\mathrm{\partial }\mathit{\rho }{}_{\mathrm{m}}}{\mathit{U}}_e\\ & =-\iota \left({\mathit{\rho }}_{\mathrm{M}}^{\mathrm{e}}\right)\sum _{e=1}^N{\mathit{U}}_e^{\mathrm{T}}\left({\int }_{{\mathrm{\Omega }}_{\mathrm{e}}}^{}{\mathit{B}}^{\mathrm{T}}\frac{\mathrm{\partial }{\mathit{E}}^{\mathrm{h}\mathrm{o}\mathrm{m}}}{\mathrm{\partial }{\mathit{\rho }}_{\mathrm{m}}}\mathit{B}\mathrm{d}{\mathrm{\Omega }}_e\right){\mathit{U}}_e.\\ & \end{array}$

The 1st-order derivatives of the homogenized elastic tensor $E^{\mathrm{hom}}$ with respect to the micro design variables can likewise be transformed into 1st-order derivatives with respect to the micro element density:

(31)$\begin{array}{rl}\frac{\mathrm{\partial }{\mathit{E}}^{\mathrm{h}\mathrm{o}\mathrm{m}}}{\mathrm{\partial }{\mathit{\rho }}_{\mathrm{m}}}& =\sum _{i\in {\mathrm{S}}_{\mathrm{m}}}^{}\frac{\mathrm{\partial }{\mathit{E}}^{\mathrm{h}\mathrm{o}\mathrm{m}}}{\mathrm{\partial }{\mathit{\rho }}_{\mathrm{m}i}^{\mathrm{e}}}\frac{\mathrm{\partial }{\mathit{\rho }}_{\mathrm{m}i}^{\mathrm{e}}}{\mathrm{\partial }{\mathit{\rho }}_{\mathrm{m}}}\\ & =\sum _{i\in {\mathrm{S}}_{\mathrm{m}}}^{}\left[\frac{p{\left({\mathit{\rho }}_{\mathrm{m}}^{\mathrm{e}}\right)}^{p-1}(1-\mathrm{\Delta })}{\left|\mathrm{\Omega }{}_{\mathrm{m}}\right|}\sum _{e=1}^n{\left({\mathit{u}}_e^{A\left(ij\right)}\right)}^{\mathrm{T}}{\mathit{k}}_0{\mathit{u}}_e^{A\left(mn\right)}\right]\frac{\mathrm{\partial }{\mathit{\rho }}_{\mathrm{m}i}^{\mathrm{e}}}{\mathrm{\partial }{\mathit{\rho }}_{\mathrm{m}}},\end{array}$

where ${\mathrm{S}}_{\mathrm{m}}$ denotes the set of elements affected by the control points corresponding to ${\mathit{\rho }}_{\mathrm{m}}$, and ${\mathit{\rho }}_{\mathrm{m}i}^{\mathrm{e}}$ denotes the ith element density in ${\mathrm{S}}_{\mathrm{m}}$.

Similarly, the sensitivity for thermal compliance with respect to $\mathit{\rho }{}_{\mathrm{M}}$ is

(32)$\frac{\mathrm{\partial }{C}_{\mathrm{H}\mathrm{e}\mathrm{a}\mathrm{t}}}{\mathrm{\partial }{\mathit{\rho }}_{\mathrm{M}}}=-\sum _{i\in {\mathrm{S}}_{\mathrm{M}}}^{}\left[p{\left({\mathit{\rho }}_{\mathrm{M}i}^{\mathrm{e}}\right)}^{p-1}\left(1-\mathrm{\Delta }\right)\sum _{e=1}^N{\mathit{T}}_e^{\mathrm{T}}{\mathit{K}}_e^{\mathrm{t}\mathrm{h}}\left({\mathit{\rho }}_{\mathrm{m}}\right){\mathit{T}}_e\right]\frac{\mathrm{\partial }{\mathit{\rho }}_{\mathrm{M}i}^{\mathrm{e}}}{\mathrm{\partial }{\mathit{\rho }}_{\mathrm{M}}}.$

Sensitivity for thermal compliance with respect to ${\rho }_m$ is

(33)$\frac{\mathrm{\partial }{C}_{\mathrm{H}\mathrm{e}\mathrm{a}\mathrm{t}}}{\mathrm{\partial }{\mathit{\rho }}_{\mathrm{m}}}=-\iota \left({\mathit{\rho }}_{\mathrm{M}}^{\mathrm{e}}\right)\sum _{e=1}^N{\mathit{T}}_e^{\mathrm{T}}\left({\int }_{{\mathrm{\Omega }}_{\mathrm{e}}}^{}{\mathit{D}}^{\mathrm{T}}\frac{\mathrm{\partial }{\mathbf{\Lambda }}^{\mathrm{h}\mathrm{o}\mathrm{m}}}{\mathrm{\partial }{\mathit{\rho }}_{\mathrm{m}}}\mathit{D}\mathrm{d}{\mathrm{\Omega }}_{\mathrm{e}}\right){\mathit{T}}_e.$

The 1st-order derivative of the homogenized thermal conductivity tensor ${\mathbf{\Lambda }}^{\mathrm{hom}}$ with respect to ${\mathit{\rho }}_{\mathrm{m}}$ is

(34)$\frac{\mathrm{\partial }{\mathbf{\Lambda }}^{\mathrm{h}\mathrm{o}\mathrm{m}}}{\mathrm{\partial }{\mathit{\rho }}_{\mathrm{m}}}=-\sum _{i\in {\mathrm{S}}_{\mathrm{m}}}^{}\left[\frac{p{\left({\mathit{\rho }}_{\mathrm{m}}^{\mathrm{e}}\right)}^{p-1}\left(1-\mathrm{\Delta }\right)}{\left|{\mathrm{\Omega }}_{\mathrm{m}}\right|}\sum _{e=1}^n{\left({\mathit{t}}_e^{A\left(ij\right)}\right)}^{\mathrm{T}}{\mathit{k}}_0^{\mathrm{t}\mathrm{h}}{\mathit{t}}_e^{A\left(mn\right)}\right]\frac{\mathrm{\partial }{\mathit{\rho }}_{\mathrm{m}i}^{\mathrm{e}}}{\mathrm{\partial }{\mathit{\rho }}_{\mathrm{m}}}.$

The sensitivity for the material volume with respect to the design variables is

(35)$\{\begin{array}{rl}& \frac{\mathrm{\partial }{V}_{\mathrm{M}}}{\mathrm{\partial }{\mathit{\rho }}_{\mathrm{M}}}=\sum _{i\in {\mathrm{S}}_{\mathrm{M}}}^{}\frac{\mathrm{\partial }{V}_{\mathrm{M}}}{\mathrm{\partial }{\mathit{\rho }}_{\mathrm{M}i}^{\mathrm{e}}}\frac{\mathrm{\partial }{\mathit{\rho }}_{\mathrm{M}i}^{\mathrm{e}}}{\mathrm{\partial }{\mathit{\rho }}_{\mathrm{M}}}=\sum _{i\in {\mathrm{S}}_{\mathrm{M}}}^{}\frac{\mathrm{\partial }{\mathit{\rho }}_{\mathrm{M}i}^{\mathrm{e}}}{\mathrm{\partial }{\mathit{\rho }}_{\mathrm{M}}},\\ & \frac{\mathrm{\partial }{V}_{\mathrm{m}}}{\mathrm{\partial }{\mathit{\rho }}_{\mathrm{m}}}=\sum _{i\in {\mathrm{S}}_{\mathrm{m}}}^{}\frac{\mathrm{\partial }{V}_{\mathrm{m}}}{\mathrm{\partial }{\mathit{\rho }}_{\mathrm{m}i}^{\mathrm{e}}}\frac{\mathrm{\partial }{\mathit{\rho }}_{\mathrm{m}i}^{\mathrm{e}}}{\mathrm{\partial }{\mathit{\rho }}_{\mathrm{m}}}=\sum _{i\in {\mathrm{S}}_{\mathrm{m}}}^{}\frac{\mathrm{\partial }{\mathit{\rho }}_{\mathrm{m}i}^{\mathrm{e}}}{\mathrm{\partial }{\mathit{\rho }}_{\mathrm{m}}}.\end{array}$

${\displaystyle \frac{\mathrm{\partial }{\mathit{\rho }}_{\mathrm{M}i}^{\mathrm{e}}}{\mathrm{\partial }{\mathit{\rho }}_{\mathrm{M}}}}$ and ${\displaystyle \frac{\mathrm{\partial }{\mathit{\rho }}_{\mathrm{m}i}^{\mathrm{e}}}{\mathrm{\partial }{\mathit{\rho }}_{\mathrm{m}}}}$ in Eqs. (29), (31), (32), (34), and (35) can be obtained in accordance with Eq. (6). The heuristic updating scheme for the design variables can be formulated as follows:

(36)$\{\begin{array}{rl}{\mathit{\rho }}_{\mathrm{M}}^{\left(s+1\right)}& =\left({\mathrm{\Phi }}_{\mathrm{M}}^{\left(s\right)}\right){\mathit{\rho }}_{\mathrm{M}}^{\left(s\right)}\\ & =\left\{-{\displaystyle \frac{\mathrm{\partial }C\left({\mathit{\rho }}_{\mathrm{M}},{\mathit{\rho }}_{\mathrm{m}}\right)}{\mathrm{\partial }\mathit{\rho }{}_{\mathrm{M}}}}\cdot {\displaystyle \frac{1}{max\left(\mu ,{\lambda }_{\mathrm{M}}^{\left(s\right)}{\displaystyle \frac{\mathrm{\partial }{V}_{\mathrm{M}}}{\mathrm{\partial }\mathit{\rho }{}_{\mathrm{M}}}}\right)}}\right\}{\mathit{\rho }}_{\mathrm{M}}^{\left(s\right)},\\ {\mathit{\rho }}_{\mathrm{m}}^{\left(s+1\right)}& =\left({\mathrm{\Phi }}_{\mathrm{m}}^{\left(s\right)}\right){\mathit{\rho }}_{\mathrm{m}}^{\left(s\right)}\\ & =\left\{-{\displaystyle \frac{\mathrm{\partial }C\left({\mathit{\rho }}_{\mathrm{M}},{\mathit{\rho }}_{\mathrm{m}}\right)}{\mathrm{\partial }{\mathit{\rho }}_{\mathrm{m}}}}\cdot \frac{1}{max\left(\mu ,{\lambda }_{\mathrm{m}}^{\left(s\right)}{\displaystyle \frac{\mathrm{\partial }V{}_{\mathrm{m}}}{\mathrm{\partial }{\mathit{\rho }}_{\mathrm{m}}}}\right)}\right\}{\mathit{\rho }}_{\mathrm{m}}^{\left(s\right)},\end{array}$

where $s$ is the current iteration number. ${\mathrm{\Phi }}_{\mathrm{M}}^{\left(s\right)}$ and ${\mathrm{\Phi }}_{\mathrm{m}}^{\left(s\right)}$ are the updating factors for the two-scale design variables ${\rho }_{\mathrm{M}}^{\left(s\right)}$ and ${\rho }_{\mathrm{m}}^{\left(s\right)}$ at the sth iteration, respectively. μ = 1E−9 is a constant to avoid a denominator of 0. ${\lambda }_{\mathrm{M}}^{\left(s\right)}$ and ${\lambda }_{\mathrm{m}}^{\left(s\right)}$ are the Lagrange multipliers of the macro and micro design variables, respectively, which can be computed by the bisectioning method. The updating schemes for macro and micro design variable are the optimality criterion method, where macro and micro updating factors ${\mathrm{\Phi }}_{\mathrm{M}}^{\left(s\right)}$ and ${\mathrm{\Phi }}_{\mathrm{m}}^{\left(s\right)}$ contain macro and micro design variables, thus depending on ${\rho }_{\mathrm{M}}^{\left(s\right)}$ and ${\rho }_{\mathrm{m}}^{\left(s\right)}$. In other words, the micro design variable can influence the optimization of the macro topology, and the macro design variable also influence the updating of the micro topology, coupling in each iteration of the concurrent topology optimization of the multiscale structures to improve the associated performance [48].

In summary, the two-scale optimization problem in Eq. (19) can be separated into two, which is exclusively associated with macro and micro design variables, respectively. With the micro design variables fixed, the sensitivity of the macro design variable is obtained and macro design variable is updated accordingly. In the same optimization step, the microscopic optimization problem is solved similarly so that both scale design variables are updated. The decoupling multiscale analysis scheme can maintain the important hierarchical and mechanical relationships between the different scales, and the computational cost is substantially less than those relying on the coupling scheme [54]. The solution efficiency of the decoupling method can still be further improved. The study [55] skillfully transforms the sensitivity of the objective function to the micro design variables into the 1st-order derivative with respect to the equivalent elasticity matrix, thus effectively improving the computational efficiency. The efficient decoupled method is used instead in subsequent studies to improve the computational efficiency of the sensitivity. The flowchart of the algorithm for multi-objective concurrent ITO for multiscale structures is shown in Fig.7.

Fig.7 Flowchart of multi-objective concurrent isogeometric topology optimization.

Full size|PPT slide

4 Multi-objective concurrent isogeometric topology optimization

Three numerical examples are presented in this section to demonstrate the feasibility and advantages of the multi-objective concurrent ITO method for multiscale structures. Example 1 verifies the effectiveness of the proposed method and the advantages of concurrent multi-scale optimization by optimizing a cantilever beam. The classical Michell-type structure was optimized via Example 2 to show further the utility of the methodology proposed in this study, and the effects of volume fraction and different boundary conditions on the results are discussed. Example 3 extends the multi-objective concurrent ITO method to three dimensions. In all examples, the sizes of the material microstructure in all normal directions are defined as 0.1. The Young’s modulus E₀ and thermal conductivity Λ₀ of the solid micro elements are 2E7 and 1, respectively, in Examples 1 and 2, while in Example 3, they are 1E4 and 1, respectively, and the Young’s modulus and thermal conductivity of the hollow micro elements are both 1E−9. The Poisson’s ratio is 0.3, and the penalization parameter p is defined as 3. The basis degree of IGA for all examples in this section is equal to 2, and the NURBS degree of the filter is also set to 2. It will not be considered convergent until the maximum change is less than 0.02 in any density of macro and micro control points between two consecutive iterations. The following arithmetic examples are derived by running them in MATLAB R2022a.

4.1 Example 1

In this example, a multi-objective concurrent ITO of mechanical and thermal conductivity is performed for a cantilever beam, and Fig.8 shows the boundary conditions and initial design at two scales. The left edge of the beam is fixed, and the force F = 4000 is loaded on the lowermost right edge and directed downward (Fig.8(a)). Fig.8(b) describes the thermal boundary conditions, where a heat sink point is located at the middle of the left edge, where the temperature is set to zero, and the remaining boundaries of the cantilever beam are set to be adiabatic (i.e., zero heat transfer). The distributed heat P = 0.0005 is applied to the entire design domain. The initial design of the macrostructure consists of the same density of elements uniformly distributed to avoid local minimum designs, while the initial design of the microstructure was filled with some voids to avoid a uniform distribution of the sensitive field due to the initial imposed periodic boundary conditions [48,51,56]. In this example, the macroscopic initial structure is uniformly distributed with a volume fraction of 30% of the macro design domain, delineated by a grid of 150 × 75 isogeometric quadratic elements, and the final optimized volume is also 30% of the design domain. Fig.8(c) shows the microstructure with a circular hole in the center of the design domain with a radius of one-fourth of the size of the microstructure, allowing the material to be unevenly distributed. A 50 × 50 quadratic isogeometric element is used to discretize the microstructure, and the final optimized volume is 50% of the microscopic design domain.

Fig.8 Boundary conditions and initial designs at two scales. (a) Mechanical boundary conditions for the macrostructure, (b) thermal boundary conditions for the macrostructure, and (c) the microstructure.

Full size|PPT slide

Tab.1 presents the optimization results for the cantilever beam. For presentation and comparison, six different weighting factors w are defined as 0, 0.2, 0.4, 0.6, 0.8, and 1. Case 1 shows the results obtained by a multi-objective concurrent ITO approach, including topology of the macrostructure, the topology of material microstructure, and the value of the single-objective function, where the iterations are less than 300 for six weighting factors. The results obtained by a multi-objective isogeometric macrostructure topology optimization for a given microstructure are presented in Case 2, where the microstructure and its equivalent properties are detailed in Tab.2, and the volume fractions are consistent with Case 1. The single-scale optimization results in [57] are used as the standard structure.

Tab.1 Optimized results of the cantilever beam

w	Case 1	Case 2
0
0.2
0.4
0.6
0.8
1.0

Tab.2 Microstructure and its effective properties in Case 2

Microstructure	Volume fraction	Homogenized elastic tensor	Homogenized thermal conductivity tensor
	0.5	${10}^6\times \left[\begin{array}{ccc}5.1511& 1.5774& 0\\ 1.5774& 5.1511& 0\\ 0& 0& 1.4863\end{array}\right]$	$\left[\begin{array}{cc}0.3117& 0\\ 0& 0.3117\end{array}\right]$

The optimized topologies of the macro and microstructure change considerably with the weighting factor w, demonstrated in Case 1. When w is equal to 0 or 1, the proposed multi-objective optimization method degenerates into a single-objective optimization method for thermal conductivity and structure accordingly, so that only the corresponding objective function value is computed, and the value of the other objective function gradually tends to infinity, which is denoted by None. As the value of w increases, the tree-like branching structure that favors heat dissipation in the macro scale gradually tends to shrink and fade, evolving in the direction of increasing structural strength. Correspondingly, the micro scale is influenced by the macro boundary conditions, and when w is equal to 0, the material is mainly distributed in the corners of the design domain, which will be favorable for uniformly dispersing the energy [31]. When w equals 1, the stress distribution induced by macroscopic boundary conditions eventually leads to an asymmetrical arrangement of materials in the microstructure, characterized by horizontal, vertical, and approximately 45° to −45° orientations. With variations in the w value, the microstructure concurrently adjusts to accommodate the evolution of the macrostructure. The objective values in Tab.1 show that the multi-objective optimized structure exhibits less mechanical compliance and less thermal conductivity compared with the thermal conductivity and mechanical single-objective optimization, respectively. Therefore, in contrast to the existing multiscale topology optimization method [58,59] limited to single physical field problems, the outcomes by this method exhibits better multi-objective performance. The more detailed optimization of macro and microstructures with different weighting factors in Case 1 and their equivalent properties are shown in Appendix (Table A1), illustrating how two conflicting objectives compete through material redistribution in the design domain of two scales when the weight factor changes.

The Pareto curve obtained by varying the weighting factor from 0 to 1 in a span of 0.1 is shown in Fig.9. As the weighting factor increases, mechanical compliance decreases, while thermal compliance increases, indicating a strong shift between the two objective functions. Given that the two objectives are different physical quantities, the convexity of all functions within the multi-objective criteria might not exist when objective functions have different ranges and orders of magnitude [34,60]. In this case, the weighted sum of the multi-objective functions must be modified, and reasonable weighting factors, such as normative weighted sum [61], the combination of power and normalization forms [62], and adaptive adjustment of weight coefficients [5], must be set. Two conflicting objectives in this study form an approximate convex curve in the allowed design space, and a strong Pareto front is produced in Fig.9. This result indicates that the weighted sum form of the multi-objective function in Eq. (19) is appropriate for two objectives with different physical quantities in this study. To ensure the convexity of multi-objective function when extending the multi-objective in this study to other domains, the follow-up work will further improve the existing weighted sum form, and the corresponding research will be reported elsewhere. Fig.10 shows the iteration of the objective function and volume fraction at two scales, indicating that the concurrent ITO method can eventually converge. Some sharp jumps in the convergence process are attributed to the use of Heaviside projection filters to ensure an approximately black and white design.

Fig.9 Pareto front of Case 1.

Full size|PPT slide

Fig.10 Iterative history of Case 1 for w = 0.6.

Full size|PPT slide

The optimized structures in Case 1 have similar structural features as the standard structures in [57], demonstrating the effectiveness of the present method. Given that both have similar boundary conditions, the macrostructures have similar topology. However, the optimized microstructure affects the macro topology to some extent, resulting in subtle differences in the topology of the macrostructures between Case 1 and the standard structure. Moreover, the mechanical compliance and thermal compliance of the structure in Case 1 is less than that in Case 2. Therefore, the proposed multi-objective concurrent ITO has a greater advantage over the results obtained by multi-objective isogeometric macrostructure topology optimization with a given microstructure under the same conditions, considering mechanical and thermal compliance.

4.2 Example 2

In this section, the effect of macro and micro volume fractions on the structural performance is discussed on the basis of the optimization results of Michell-type structure obtained from the multi-objective concurrent ITO approach. Fig.11 shows the boundary conditions and the initial design for two scales. The lower left corner of the beam is fixed, while the lower right corner is supported by rollers. The force F = 3000 is loaded in the center of the lower edge. A heat sink also exists in the middle of the lower edge, where the temperature is set to zero and the rest of the boundary of the Michell-type beam is set to be adiabatic (i.e., zero heat transfer). The distributed heat P = 0.0005 is applied to the entire design domain. The macro initial structure is divided by a grid of 125 × 50 isogeometric quadratic elements. The microstructure is discretized by setting a circular hole in the center of the design domain with a radius of one-third of the size of the microstructure, slightly different from that in example 1, and 50 × 50 quadratic isogeometric elements are used to discretize the microstructure. Four different combinations of macro and micro volume fractions are defined, where the macro volume fractions are 0.3, 0.4, 0.5, and 0.6. The total volume fraction remains constant at 0.12, and the micro volume fractions are equal to 0.4, 0.3, 0.24, and 0.2, respectively.

Fig.11 Initial design and boundary conditions for two scales in Example 2. (a) Boundary conditions for the macrostructure and (b) the microstructure.

Full size|PPT slide

Tab.3 shows the macro and micro concurrent design of Michell-type structures and the corresponding two-objective function values for the four cases. With the same weighting factor, as the macro volume fraction increases, the micro material gradually transfers to the macrostructure. The basic configurations of the macro and microstructure remain unchanged, but the mechanical compliance and thermal compliance of the optimized structures gradually increase. To investigate the real performance of the single-scale structure with $V_{\text{m}}$ = 1, a set of single-scale Michell-type structures is further added, as shown in Fig.12. The results show that when w = 0, the thermal compliance of the single-scale design increases in comparison with the four cases in Tab.3. With the increase in the macro volume fraction, the thermal conductivity is first enhanced and then weakened, i.e., the thermal conductivity of the single-scale structure is inferior to that of the multiscale design. When w = 1, the mechanical compliance still maintains a decreasing trend, indicating that the stiffness performance of the single-scale structure is superior to that of the multiscale design, which is consistent with the conclusion in Ref. [48]. The multi-objective performance of most of the multiscale designs for w not equal to 0 or 1 outperforms that of the single-scale designs.

Tab.3 Optimized results of the Michell-type structure in three volume fraction cases

w	Case 3 V M = 0.3, V m = 0.4	Case 4 V M = 0.4, V m = 0.3	Case 5 V M = 0.5, V m = 0.24	Case 6 V M = 0.6, V m = 0.2
0
0.2
0.4
0.6
0.8
1.0

Fig.12 Single-scale optimized results of the Michell-type structure when VM = 0.12 and Vm = 1.

Full size|PPT slide

Fig.13 illustrates the Pareto curves for five sets of volume fractions at different weighting factors, where the single-objective optimization results for $w$ equal to 0 and 1 are removed. The more intuitive depiction of the observed pattern is shown: as the macro volume fraction decreases and the micro volume fraction increases, the mechanical and thermal compliance multi-objective values first decrease gradually and then increase inversely. Therefore, for the optimization of the two objectives, i.e., mechanical compliance and thermal compliance, the best design for multi-objective performance is the multiscale optimization in Case 3. In addition, multiscale structures exhibit better multi-objective performance than at single scales when generalized to a wide range of other topological design problems, such as two scales involving thermal insulation and stiffness [36] and natural frequency and thermal conductivity [8]. The iterative history of objective functions and volume fractions at both scales in Case 3 when $w$ = 0.2 is displayed in Fig.14.

Fig.13 Pareto front of Example 2.

Full size|PPT slide

Fig.14 Iterative history of Case 3 for w = 0.2.

Full size|PPT slide

With Case 3 taken as an example, when the weighting factor w is equal to 0, the optimized macrostructure presents the material extending from the middle of the lower side of the design domain to the ends of the upper side to form two mutually symmetric root-like structures as the main heat transfer path. Accordingly, the influence of the macro scale on the micro design compels the microscale materials to concentrate in regions parallel to the diagonal. The micro material distribution matches with the macro design structure, enhancing the efficient heat transfer from various corners of the macro design domain to the middle of the bottom edge along these two rooted structures, as illustrated in Fig.15(a). When the weighting factor w is equal to 1, the macro and microstructure similar to that in [48] are obtained. As the weight factor gradually increases, the mechanical compliance gradually decreases, while the thermal conductivity compliance gradually increases, and the macro and micro topological structures evolve toward increasing the structural strength. Balancing the decisions between these two conflicting objectives results in a series of compromise designs.

Fig.15 Influence of boundary conditions in Case 3 for w = 0. (a) Heat sink point at the middle of the bottom edge, (b) heat sink point at the left quarter bottom edge, (c) heat sink point at the lowermost left edge.

Full size|PPT slide

To explore further how the layout design of microstructure aligns with the macro scale under different boundary conditions to enhance target performance (manifested as thermal performance when w = 0), structural optimization was conducted for the volume fractions in Case 3 under three boundary conditions, as depicted in Fig.15. When the heat dissipation point is positioned at the middle of the bottom edge, Fig.15(a) reveals that micro materials concentrate in regions parallel to the two diagonals of the micro design domain, forming two primary thermal channels. This adaptation is to facilitate the uniformly distributed heat transfer from two approximately 45° and −45° rooted structures to the heat dissipation point in the macroscopic design. When micro elements are located in the left-rooted structure of the macroscopic design domain, thermal conduction primarily occurs through thermal channel 1, and correspondingly, micro elements in the right-rooted structure primarily rely on thermal channel 2. As the heat sink point is moved from the middle of the bottom edge to the left quarter, Fig.15(b) illustrates that two asymmetric root-like structures extend from the lower part of the macroscopic design domain to the upper ends. Compared with Fig.15(a), the left-rooted structure contracts, and the right-rooted structure expands. To match macrostructure, the width of the thermal channel 1 in the micro design is reduced accordingly, the width of the thermal channel 2 is increased, and the two channels are no longer symmetric. When the heat sink point is further moved to the far left, the width and direction of the two thermal channels in the micro design changed accordingly.

4.3 Example 3

In this example, the multi-objective concurrent ITO approach is extended to the 3D case, with the macro boundary conditions illustrated in Fig.16. Four corners at the bottom side of the 3D-supported structure are fixed. The downward force $F_y=4000$ and rightward force $F_{\mathit{\text{z}}}=4000$ are loaded at the center of the top surface along the y- and z-directions, respectively. Simultaneously, a heat sink is positioned at the center of top surface with a temperature of 0, while the remaining surfaces are adiabatic. A uniform thermal load of 0.0005 is distributed over the macro design domain. The macro initial structure is uniformly distributed with a 30% volume fraction of the macro design domain, divided into a 15 × 12 × 15 quadratic isogeometric element mesh, as shown in Fig.17. The final optimized volume is also 30% of the design domain. The micro initial structure has a spherical hole in the center of the design domain with a radius equal to one-sixth of the microstructure size. Discretization of the microstructure is achieved using a 15 × 15 × 15 quadratic isogeometric element mesh. The final optimized volume constitutes 40% of the microscopic design domain.

Fig.16 Initial design and boundary conditions for two scales in Example 3.

Full size|PPT slide

Fig.17 IGA mesh of macrostructure in Example 3.

Full size|PPT slide

Fig.18 shows the final optimized topology of macro and microstructures for weighting factors of 0, 0.6, and 1.0, respectively. Increasing the number of elements improves the clarity of the design but also incurs substantial computational cost. Therefore, a rough isogeometric mesh is chosen for the optimization analysis, and smoothing was applied to the final results. The density interpolation is performed for the coarse macroscopic and microscopic density elements, and then the isosurface technique is utilized to obtain smoother results. The final results show that the 3D-supported structure exhibits a root-like structure when w = 0 and divides into four main branches. The micro materials form a bottom-to-top heat transfer channel. When w = 1, the macrostructure of the 3D-supported structure is similar to [48]. The difference is that the material is more concentrated in two of the supports to resist the force along the z-direction. At w = 0.6, an intermediate shape between structural optimization and thermal optimization is presented. Tab.4 illustrates that the optimized microstructures and their equivalent properties, which are obtained by applying periodic boundary condition along three directions. Therefore, multi-objective concurrent ITO is equally effective for solving 3D optimization problems.

Tab.4 Microstructures and their equivalent properties in Example 3

Microstructure of concurrent design	Homogenized elastic tensor (107 ×)	Homogenized thermal conductivity tensor
	$\left[\begin{array}{cccccc}0.9339& 0.4977& 0.7282& 0& 0& 0\\ 0.4977& 2.6247& 0.4977& 0& 0& 0\\ 0.7282& 0.4977& 0.9339& 0& 0& 0\\ 0& 0& 0& 0.5345& 0& 0\\ 0& 0& 0& 0& 0.5345& 0\\ 0& 0& 0& 0& 0& 0.5139\end{array}\right]$	$\left[\begin{array}{ccc}0.1393& 0& 0\\ 0& 0.2327& 0\\ 0& 0& 0.1393\end{array}\right]$
	$\left[\begin{array}{cccccc}1.2016& 0.4754& 0.5269& 0.1515& 0& 0\\ 0.4754& 1.7558& 0.4582& 0.1938& 0& 0\\ 0.5269& 0.4582& 1.0222& 0.0045& 0& 0\\ 0.1515& 0.1938& 0.0045& 0.4769& 0& 0\\ 0& 0& 0& 0& 0.4412& 0.0077\\ 0& 0& 0& 0& 0.0077& 0.3889\end{array}\right]$	$\left[\begin{array}{ccc}0.1433& 0.019& 0\\ 0.019& 0.1801& 0\\ 0& 0& 0.1322\end{array}\right]$
	$\left[\begin{array}{cccccc}0.7768& 0.4235& 0.3545& 0.2285& 0& 0\\ 0.4235& 1.3937& 0.4192& 0.3465& 0& 0\\ 0.3545& 0.4192& 1.2488& 0.0499& 0& 0\\ 0.2285& 0.3465& 0.0499& 0.4208& 0& 0\\ 0& 0& 0& 0& 0.4127& 0.1056\\ 0& 0& 0& 0& 0.1056& 0.2841\end{array}\right]$	$\left[\begin{array}{ccc}0.1069& 0.0424& 0\\ 0.0424& 0.1503& 0\\ 0& 0& 0.1353\end{array}\right]$

Fig.18 Concurrent design of the 3D-supported structures. (a) w = 0; (b) w = 0.6; (c) w = 1.

Full size|PPT slide

5 Conclusions

In this study, a multi-objective concurrent ITO method based on the SIMP method is proposed for the design of multiscale structures with multiple physical properties, such as high thermal conductivity and stiffness. The multi-objective function is defined using a linear weighted sum method to assign weighting factors to each objective. Multiple effective properties of the material microstructure are computed by the energy-based homogenization method and sensitivity analyses are conducted at the macro and micro scales by the decoupled method to update the macro and microstructures simultaneously. Given the property of the NURBS mesh, traditional distance-based filters are unnecessary for ITO.

Numerical examples demonstrate that the proposed multi-objective concurrent ITO method is feasible and advantageous. By optimizing the structure at both scales simultaneously, the mechanical and thermal compliance of the multiscale structure is considerably reduced. The key transitions between conflicting objective functions lead to a strong Pareto front, indicating that the weighted sum form of the multi-objective function is appropriate in this study. The results at different volume fractions show that the transfer of macro materials to micro structures improves the mechanical and thermal conductivity. When the micro volume fraction is further increased, the multi-objective properties of mechanical and thermal conductivity are instead reduced, indicating that the multiscale design is superior to single-scale structure. In addition, the concurrent design not only adds the additional design freedom but also motivates the coordination of the micro layout with the macro structure to improve performance. The successfully realized 3D case further validates the effectiveness and practicality of the method.

Although the multi-objective concurrent ITO method was utilized for thermal conductivity and stiffness optimization in this study, it could be extended to wide range of other topological design problems. To ensure the convexity of multi-objective function, the existing weighted sum form should be further improved. The decoupled method is also not efficient enough for sensitivity analysis in concurrent topology optimization. In the future, we will improve this method to address the above issues, and the corresponding research will be reported elsewhere.

References

[1] ZhengX Y, Lee H, WeisgraberT H, ShusteffM, DeOtteJ, DuossE B, Kuntz J D, BienerM M, GeQ, Jackson J A, KucheyevS O, FangN X, Spadaccini C M. Ultralight, ultrastiff mechanical metamaterials. Science, 2014, 344(6190): 1373–1377

[2] ShaikeeaA J D, CuiH C, O’Masta M, ZhengX Y R, DeshpandeV S. The toughness of mechanical metamaterials. Nature Materials, 2022, 21(3): 297–304

[3] PejmanR, Najafi A R. Multiphysics topology optimization of a multifunctional structural battery composite. Structural and Multidisciplinary Optimization, 2023, 66(3): 46

[4] WangY H, Sha W, XiaoM, QiuC W, GaoL. Deep-learning-enabled intelligent design of thermal metamaterials. Advanced Materials, 2023, 35(33): 2302387

[5] ZhengY F, Fu Z J, WangY J, LuX, QuJ P, ZhangC Z. Hierarchical design of material microstructures with thermal insulation properties. International Journal of Heat and Mass Transfer, 2022, 186: 122514

[6] SuX N, Chen W J, LiuS T. Multi-scale topology optimization for minimizing structural compliance of cellular composites with connectable graded microstructures. Structural and Multidisciplinary Optimization, 2021, 64(4): 2609–2625

[7] Giraldo-LondoñoO, MirabellaL, Dalloro L, PaulinoG H. Multi-material thermomechanical topology optimization with applications to additive manufacturing: design of main composite part and its support structure. Computer Methods in Applied Mechanics and Engineering, 2020, 363: 112812

[8] ChenW J, Zheng Y F, WangY J. Multi-objective topology optimization filled with multiple microstructures. Composite Structures, 2023, 304: 116322

[9] BendsøeM P, KikuchiN. Generating optimal topologies in structural design using a homogenization method. Computer Methods in Applied Mechanics and Engineering, 1988, 71(2): 197–224

[10] SigmundO. A 99 line topology optimization code written in Matlab. Structural and Multidisciplinary Optimization, 2001, 21(2): 120–127

[11] AndreassenE, Clausen A, SchevenelsM, LazarovB S, Sigmund O. Efficient topology optimization in MATLAB using 88 lines of code. Structural and Multidisciplinary Optimization, 2011, 43(1): 1–16

[12] HuangX, Xie Y M. Bi-directional evolutionary topology optimization of continuum structures with one or multiple materials. Computational Mechanics, 2009, 43(3): 393–401

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

[14] GuoX, ZhangW S, ZhongW L. Doing topology optimization explicitly and geometrically—a new moving morphable components based framework. Journal of Applied Mechanics, 2014, 81(8): 081009

[15] KiranR, Nguyen-Thanh N, HuangJ Z, ZhouK. Buckling analysis of cracked orthotropic 3D plates and shells via an isogeometric-reproducing kernel particle method. Theoretical and Applied Fracture Mechanics, 2021, 114: 102993

[16] ZhangJ P, Luo T, ZhangD B, YinS H, HeH L, PengJ P. Multi-objective topology optimization of thermal-mechanical coupling anisotropic structures using the isogeometric analysis approach. Applied Mathematical Modelling, 2023, 117: 267–285

[17] HughesT J R, Cottrell J A, BazilevsY. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Computer Methods in Applied Mechanics and Engineering, 2005, 194(39–41): 4135–4195

[18] KatoJ, Yachi D, TeradaK, KyoyaT. Topology optimization of micro-structure for composites applying a decoupling multi-scale analysis. Structural and Multidisciplinary Optimization, 2014, 49(4): 595–608

[19] YanX, HuangX, ZhaY, XieY M. Concurrent topology optimization of structures and their composite microstructures. Computers & Structures, 2014, 133: 103–110

[20] ChenW J, Tong L Y, LiuS T. Concurrent topology design of structure and material using a two-scale topology optimization. Computers & Structures, 2017, 178: 119–128

[21] GroenJ P, Thomsen C R, SigmundO. Multi-scale topology optimization for stiffness and de-homogenization using implicit geometry modeling. Structural and Multidisciplinary Optimization, 2021, 63(6): 2919–2934

[22] YuC, WangQ F, XiaZ H, Wang Y J, MeiC, LiuY H. Multiscale topology optimization for graded cellular structures based on level set surface cutting. Structural and Multidisciplinary Optimization, 2022, 65(1): 32

[23] CoelhoP G, Cardoso J B, FernandesP R, RodriguesH C. Parallel computing techniques applied to the simultaneous design of structure and material. Advances in Engineering Software, 2011, 42(5): 219–227

[24] GaoJ, CaoX F, XiaoM, Yang Z Q, ZhouX Q, LiY, GaoL, YanW T, Rabczuk T, MaiY W. Rational designs of mechanical metamaterials: formulations, architectures, tessellations and prospects. Materials Science and Engineering R, 2023, 156: 100755

[25] LiS S, Zhu Y C, GuoX. Optimisation of spatially varying orthotropic porous structures based on conformal mapping. Computer Methods in Applied Mechanics and Engineering, 2022, 391: 114589

[26] LiuH, ZongH M, ShiT L, Xia Q. M-VCUT level set method for optimizing cellular structures. Computer Methods in Applied Mechanics and Engineering, 2020, 367: 113154

[27] GroenJ P, Wu J, SigmundO. Homogenization-based stiffness optimization and projection of 2D coated structures with orthotropic infill. Computer Methods in Applied Mechanics and Engineering, 2019, 349: 722–742

[28] LiuP, KangZ, LuoY J. Two-scale concurrent topology optimization of lattice structures with connectable microstructures. Additive Manufacturing, 2020, 36: 101427

[29] LiQ H, Xu R, WuQ B, LiuS T. Topology optimization design of quasi-periodic cellular structures based on erode–dilate operators. Computer Methods in Applied Mechanics and Engineering, 2021, 377: 113720

[30] Hu Y, LuoY F, LiuS T. Two-scale concurrent topology optimization method of hierarchical structures with self-connected multiple lattice-material domains. Composite Structures, 2021, 272: 114224

[31] AlAli M, Shimoda M. Investigation of concurrent multiscale topology optimization for designing lightweight macrostructure with high thermal conductivity. International Journal of Thermal Sciences, 2022, 179: 107653

[32] HuangX, Xie Y M, JiaB, LiQ, ZhouS W. Evolutionary topology optimization of periodic composites for extremal magnetic permeability and electrical permittivity. Structural and Multidisciplinary Optimization, 2012, 46(3): 385–398

[33] BayatM, Zinovieva O, FerrariF, AyasC, Langelaar M, SpangenbergJ, SalajegheR, Poulios K, MohantyS, SigmundO, HattelJ. Holistic computational design within additive manufacturing through topology optimization combined with multiphysics multi-scale materials and process modelling. Progress in Materials Science, 2023, 138: 101129

[34] AlAli M, Shimoda M. Toward multiphysics multiscale concurrent topology optimization for lightweight structures with high heat conductivity and high stiffness using MATLAB. Structural and Multidisciplinary Optimization, 2022, 65(7): 207

[35] TakezawaA, Yoon G H, JeongS H, KobashiM, Kitamura M. Structural topology optimization with strength and heat conduction constraints. Computer Methods in Applied Mechanics and Engineering, 2014, 276: 341–361

[36] YanX L, Huang X D, SunG Y, XieY M. Two-scale optimal design of structures with thermal insulation materials. Composite Structures, 2015, 120: 358–365

[37] XuB, HuangX, ZhouS W, Xie Y M. Concurrent topological design of composite thermoelastic macrostructure and microstructure with multi-phase material for maximum stiffness. Composite Structures, 2016, 150: 84–102

[38] WangY, Gao J, LuoZ, BrownT, ZhangN. Level-set topology optimization for multimaterial and multifunctional mechanical metamaterials. Engineering Optimization, 2017, 49(1): 22–42

[39] AlAli M, Shimoda M. On concurrent multiscale topology optimization for porous structures under hygro-thermo-elastic multiphysics with considering evaporation. International Journal for Numerical Methods in Engineering, 2023, 124(14): 3219–3249

[40] BaoY H, Wei Z S, JiaZ Y, WangD Z, ZhangX P, KangZ. Mechanical metamaterial design with the customized low-frequency bandgap and negative Poisson’s ratio via topology optimization. Extreme Mechanics Letters, 2024, 67: 102124

[41] GaoJ, LuoZ, XiaoM, Gao L, LiP G. A NURBS-based Multi-Material Interpolation (N-MMI) for isogeometric topology optimization of structures. Applied Mathematical Modelling, 2020, 81: 818–843

[42] KiyonoC Y, Vatanabe S L, SilvaE C N, ReddyJ N. A new multi-p-norm formulation approach for stress-based topology optimization design. Composite Structures, 2016, 156: 10–19

[43] WangY J, Xiao M, XiaZ H, LiP G, GaoL. From Computer-Aided Design (CAD) toward Human-Aided Design (HAD): an isogeometric topology optimization approach. Engineering, 2023, 22: 94–105

[44] SeoY D, Kim H J, YounS K. Isogeometric topology optimization using trimmed spline surfaces. Computer Methods in Applied Mechanics and Engineering, 2010, 199(49–52): 3270–3296

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

[46] GaiY D, Xing J, HuP. Efficient MATLAB implementation of NURBS-based IGA and material design using isogeometric topology optimization. Optimization and Engineering, 2023, 24(3): 1773–1808

[47] XiaL, Breitkopf P. Design of materials using topology optimization and energy-based homogenization approach in Matlab. Structural and Multidisciplinary Optimization, 2015, 52(6): 1229–1241

[48] GaoJ, LuoZ, XiaL, GaoL. Concurrent topology optimization of multiscale composite structures in Matlab. Structural and Multidisciplinary Optimization, 2019, 60(6): 2621–2651

[49] AndreassenE, Andreasen C S. How to determine composite material properties using numerical homogenization. Computational Materials Science, 2014, 83: 488–495

[50] WangY J, Liao Z Y, YeM, ZhangY, LiW H, XiaZ H. An efficient isogeometric topology optimization using multilevel mesh, MGCG and local-update strategy. Advances in Engineering Software, 2020, 139: 102733

[51] SigmundO. Materials with prescribed constitutive parameters: an inverse homogenization problem. International Journal of Solids and Structures, 1994, 31(17): 2313–2329

[52] WangD D, Xuan J C. An improved NURBS-based isogeometric analysis with enhanced treatment of essential boundary conditions. Computer Methods in Applied Mechanics and Engineering, 2010, 199(37–40): 2425–2436

[53] GaoJ, LiH, GaoL, XiaoM. Topological shape optimization of 3D micro-structured materials using energy-based homogenization method. Advances in Engineering Software, 2018, 116: 89–102

[54] KatoJ, Yachi D, KyoyaT, TeradaK. Micro-macro concurrent topology optimization for nonlinear solids with a decoupling multiscale analysis. International Journal for Numerical Methods in Engineering, 2018, 113(8): 1189–1213

[55] ZhaoJ P, Yoon H, YounB D. An efficient decoupled sensitivity analysis method for multiscale concurrent topology optimization problems. Structural and Multidisciplinary Optimization, 2018, 58(2): 445–457

[56] BendsøeM P, SigmundO. Topology Optimization: Theory, Methods, and Applications. 2nd ed. Berlin: Springer, 2003

[57] DenkM, Rother K, PaetzoldK. Multi-objective topology optimization of heat conduction and linear elastostatic using weighted global criteria method. In: Proceedings of the 31st Symposium Design for X (DFX2020). The Design Society, 2020, 91–100

[58] YuC, WangQ F, MeiC, XiaZ H. Multiscale isogeometric topology optimization with unified structural skeleton. Computer Modeling in Engineering & Sciences, 2020, 122(3): 779–803

[59] AminzadehM, Tavakkoli S M. Multiscale topology optimization of structures by using isogeometrical level set approach. Finite Elements in Analysis and Design, 2024, 235: 104167

[60] MarlerR T, Arora J S. The weighted sum method for multi-objective optimization: new insights. Structural and Multidisciplinary Optimization, 2010, 41(6): 853–862

[61] CarpinelliG, Caramia P, MottolaF, ProtoD. Exponential weighted method and a compromise programming method for multi-objective operation of plug-in vehicle aggregators in microgrids. International Journal of Electrical Power & Energy Systems, 2014, 56: 374–384

[62] WuF H, Wang Z H, HanJ X, PeiG B. Research on multiobjective topology optimization of diesel engine cylinder block based on analytic hierarchy process. Mathematical Problems in Engineering, 2019, 2019(1): 6194634

Acknowledgements

This work was supported by the National Key R&D Program of China (Grant No. 2022YFB3302900) and the National Natural Science Foundation of China (Grant No. 52475261).

Conflict of Interest

The authors declare no conflict of interest.

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 appropriate credit is given to the original author(s) and source, a link to the Creative Commons license is provided, and the changes made are indicated.

The images or other third-party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license 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.

Visit https://creativecommons.org/licenses/by/4.0/ to view a copy of this license.

RIGHTS & PERMISSIONS

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