# Parametric Level Set-Based Multimaterial Topology Optimization of Heat Conduction Structures

Authors: Yadong Shen; Jianhu Feng
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

This paper presents a parametric level set-based method (PLSM) for multimaterial topology optimization of heat conduction structures with volume constraints. A parametric level set-based optimization model of heat conduction structures is built with multimaterial level set (MM-LS) model, which describes the boundaries of different materials by the combination of all level set functions. The heat dissipation efficiency which means the quadratic temperature gradient is conducted as the objective function. The adjoint method is utilized to calculate the sensitivities of the objective function with respect to expansion coefficients of the compactly supported radial basis functions (CSRBFs). The optimal configuration is achieved by updating the expansion coefficients gradually with the method of moving asymptotes (MMA). Several numerical examples are discussed to demonstrate effectiveness of the proposed method for multimaterial topology optimization of heat conduction structures.


---

## Body

## 1. Introduction

In a high temperature environment, the performance, and efficiency of the structure tend to decrease. Therefore, the layout design of heat conduction structures is an important issue in the field of structural design. In conventional methods, size optimization and shape optimization are used in design of heat conduction structures. Unlike these passive methods, topology optimization is an active method to find the optimal layout of materials for achieving the minimized objective within a given design domain. Based on the advantages of convenience and high efficiency, topology optimization becomes an effective way for the design of heat conduction structures [1].In recent years, different topology optimization methods have been applied to the layout design of heat conduction structures. Variable density method with material interpolation model of solid isotropic material with penalization (SIMP) was first introduced by Sigmund [1] to solve the heat conduction problem by using optimality criteria (OC) in two-dimensional space. And it is also extended into three-dimensional space by Liu et al. [2]. Gersborg-Hansen et al. [3] applied the variable density method into heat conduction problem, and the optimal structure was obtained by the finite volume method. Li et al. [4, 5] employed the evolutionary structural optimization (ESO) to solve the steady heat conduction problem. In this method, material elements which had less contribution to the objective function were allowed to be removed gradually. Bidirectional evolutionary structural optimization (BESO) also had been utilized to solve heat conduction problem by He et al. [6]. The optimal structure was achieved by gradually adding material elements as well as removing according to a certain criterion. However, these methods involve the problems of checkerboards pattern and mesh-dependent phenomena. Furthermore, the optimal structures which have zigzag boundaries are hard to manufacture.Comparing to SIMP, ESO, and other material distribution methods, level set method (LSM) has the following advantages:(1) the structural boundary is described by geometrical parameters explicitly which is easier to convert to CAD model for manufacturing; (2) it is without checkerboards pattern and mesh-dependent phenomena; and (3) the optimal configuration does not depend on the finite element meshing. Based on these advantages of LSM, it has been widely studied in different problems, such as elastic structure problem [7, 8], vibration problem [9], and compliant mechanisms problem [10]. For heat conduction problem [11], Ha and Cho [12] developed the weak form of equilibrium equation for heat conduction problem. The adjoint method and the Lagrange multiplier method were employed to derive the shape design sensitivity. Zhuang et al. [13] investigated the topology optimization of heat conduction problem under multiple load cases and developed an effective numerical procedure to implement the topological structure. Kim et al. [14] developed the topology optimization method which used the combination of topological derivatives to improve convergence of optimization process as well as avoid local minimum. However, there are also some shortcomings in classic level set method. Firstly, as the optimal structure is obtained by solving the Hamilton-Jacobi partial differential equation (HJ-PDE) with explicit time-marching algorithm, and the equations only can be implemented under some strict conditions [15, 16], it has a low computational efficiency. Secondly, the Courant-Friedrichs-Lewy (CFL) condition and reinitialization are also needed to keep the numerical stability. In addition, the optimal structure generally depends on the initial topology as the new holes cannot be created in the design domain.Parametric level set-based method (PLSM) solves the topology optimization problems effectively without directly solving the HJ-PDE. In this method, the implicit level set function is interpolated by radial basis functions (RBFs); the time dependence of the implicit function is transformed into the time dependence of interpolation coefficients. Therefore, the HJ-PDE is converted into a series of ordinary differential equations which are easier to handle, and the propagation of boundaries is driven by updating interpolation coefficients using mathematical programming methods. The PLSM not only inherits the favorable features of the conventional LSM but also avoids the difficulties in solving the HJ-PDE. Wang et al. [17] applied globally supported radial basis functions (GSRBFs) into topology optimization of minimum compliance design. Their work illustrated the efficiency and accuracy of RBF-based LSM. Luo et al. [18] introduced the compactly supported radial basis functions (CSRBFs) [19, 20] into topology optimization of compliant mechanisms and established the optimization model of a nonconvex objective function with specified constraints. The evolution of the boundary was the process of updating the interpolation coefficients using method of moving asymptotes (MMA) [21]. Luo et al. [22, 23] also proposed effective methods for structural stiffness design which combined the CSRBFs with well-established optimization algorithms, such as MMA and OC. With the initial level set functions being interpolated by CSRBFs, the optimization algorithms were applied to iteratively updating the interpolation coefficients to achieve the optimal structure.Up to date, multimaterial topology optimization design has been widely investigated [24–27]. LSM based method is an effective way that the boundaries of different materials can be described explicitly using the combination of all the level set functions. Wang et al. [28] developed the level set model of multimaterial with the minimization of the mean compliance under the volume constraints. They also showed the effectiveness of this approach in compliant mechanisms design of multimaterial [10]. Zhuang et al. [29] presented the multimaterial level set model of the heat conduction structures. The optimization problem was formulated as minimizing the heat compliance under the volume constraints for different materials. PLSM also had been applied into multimaterial topology optimization. Wang et al. [30] proposed a new multimaterial level set (MM-LS) model, in which M materials and void material were represented by M level set functions. The optimal configuration was achieved by iteratively updating the interpolation coefficients of different level set functions until convergence. Wang et al. [31] investigated the optimization design of mechanical metamaterials for multimaterial by PLSM. In their work, they utilized the homogenization method to evaluate the properties of microstructure, while PLSM was used to implement the evolution of microstructure. The combination of all level set functions represented the boundaries of multimaterial microstructure. Chu et al. [32] applied the stress interpolation scheme into design of compliant mechanisms for multimaterial. The final results with the proper distribution of multimaterial decreased the output displacement and the compliance of compliant mechanisms and satisfied the stress constraints of different materials simultaneously.In this paper, the multimaterial topology optimization model of the heat conduction structures using PLSM is investigated. The optimization is conducted to minimize the quadratic temperature gradient in the whole design domain and subject to the volume constraints of different materials. The remainder of this paper is organized as follows. Section2 reviews the MM-LS model for topology optimization of multimaterial. Section 3 discusses the multimaterial topology optimization model for heat conduction structures based on PLSM. Section 4 gives the optimization algorithm. Section 5 presents some numerical examples to demonstrate the validity of the proposed method. Conclusions are given in Section 6.


## 2. MM-LS Model

In this section, the MM-LS model for multimaterial topology optimization is reviewed, in which level set functions are used to describe materials and the void phase.


### 2.1. Level Set Representation for Multimaterial

In level set-based topology optimization, zero level sets embedded in a higher dimensional scalar function are used to describe the structural boundaries of different materials. The level set functions are used to denote the following:(1)Φkx>0,∀x∈Ωk\ΓkΦkx=0,∀x∈Γkk=1,2,…,MΦkx<0,∀x∈D\Ωk∪Γkwhere Φk(x) is the kth level set function, D is the design domain, Ωk is the domain with the positive value of the kth level set function, and Γk denotes the zero level set of the kth level set function.A pseudo-timet is introduced to manage Φk(x). The HJ-PDE can be solved by differentiating both sides of Φkxt,t=0 with respect to t:(2)∂Φkx,t∂t-vnk∇Φkx,t=0where vnk is the normal velocity of the kth level set function, as the normal component of velocity contributes to the shape evolution of the boundary. The finite difference method like up-wind scheme [33] is used to solve the HJ-PDE.


### 2.2. MM-LS Model

To solve the multimaterial design problem, the MM-LS model is applied to represent the different materials. In this method, a combination ofM level set functions is used to describe the M+1 phases, including M materials and the void area. For an M-material structure, the characteristic function ρiΦ for the ith material can be defined as follows:(3)ρiΦ=HΦ1,M=1∏k=1iHΦk1-HΦi+1,i=1,2,…,M-1∏k=1iHΦk,i=M,M≥2where HΦk is the Heaviside function of kth level set function. The MM-LS model which involves two materials and the void area by the combination of two level set functions is illustrated in Figure 1.Figure 1
MM-LS model for two materials and the void area.Based on (3), the heat conductivity CΦ at computational point x for design problem can be written as(4)CΦ=∑i=1MρiΦciwhere ci is the heat conductivity of ith material.


## 2.1. Level Set Representation for Multimaterial

In level set-based topology optimization, zero level sets embedded in a higher dimensional scalar function are used to describe the structural boundaries of different materials. The level set functions are used to denote the following:(1)Φkx>0,∀x∈Ωk\ΓkΦkx=0,∀x∈Γkk=1,2,…,MΦkx<0,∀x∈D\Ωk∪Γkwhere Φk(x) is the kth level set function, D is the design domain, Ωk is the domain with the positive value of the kth level set function, and Γk denotes the zero level set of the kth level set function.A pseudo-timet is introduced to manage Φk(x). The HJ-PDE can be solved by differentiating both sides of Φkxt,t=0 with respect to t:(2)∂Φkx,t∂t-vnk∇Φkx,t=0where vnk is the normal velocity of the kth level set function, as the normal component of velocity contributes to the shape evolution of the boundary. The finite difference method like up-wind scheme [33] is used to solve the HJ-PDE.


## 2.2. MM-LS Model

To solve the multimaterial design problem, the MM-LS model is applied to represent the different materials. In this method, a combination ofM level set functions is used to describe the M+1 phases, including M materials and the void area. For an M-material structure, the characteristic function ρiΦ for the ith material can be defined as follows:(3)ρiΦ=HΦ1,M=1∏k=1iHΦk1-HΦi+1,i=1,2,…,M-1∏k=1iHΦk,i=M,M≥2where HΦk is the Heaviside function of kth level set function. The MM-LS model which involves two materials and the void area by the combination of two level set functions is illustrated in Figure 1.Figure 1
MM-LS model for two materials and the void area.Based on (3), the heat conductivity CΦ at computational point x for design problem can be written as(4)CΦ=∑i=1MρiΦciwhere ci is the heat conductivity of ith material.


## 3. PLSM for Multimaterial Topology Optimization of Heat Conduction Structures

### 3.1. Parametric Level Set

The PLSM with CSRBF which is a kind of radial basis function centered at a particular point with compact support is applied into multimaterial topology optimization of heat conduction structures. The level set surface can be approximated by CSRBF over the design domain. The level set function can be expressed as(5)Φkx,t=∑j=1Nwjxαjkt,k=1,2,…,Mwhere N is the total number of computational points. wj denotes the CSRBF at jth point. αjk is the expansion coefficient at jth point in kth level set function.In this work, CSRBF with C2 continuity [19, 20] is employed to approximate the level set function, which can be formulated as follows:(6)wjx=max⁡0,1-rjx44rjx+1j=1,2,…,Nwhere r is the radius of support, and it is usually defined in a 2D Euclidean space as(7)rjx=x-xjR,j=1,2,…Nwhere xj is the coordinate of the jth computational point and R is the radius of support domain of the computational point x, which is defined as 2.0-3.0.Substituting (5) into (2), the HJ-PDE in (2) can be rewritten as follows:(8)∑j=1Nwjxα˙jkt-vnk∑j=1N∇wjxαjkt=0,k=1,2,…,MThe normal velocity of the level set function can be calculated by(9)vnk=∑j=1Nwjxα˙jkt∑j=1N∇wjxαjktwhere α˙jk(t)=dαjk(t)/dt.


### 3.2. Formulations of Heat Conduction Problem

For a steady heat conduction optimization problem in two dimensions, the governing equation is given as(10)∇·C∇u+Q=0inΩu=0onΓDC∇u·n=0onΓNwhere C is the heat conduction coefficient, u is the temperature field, Q is the heat source, and n is the unit normal vector of boundary. Ω is a closed domain with boundaries; the boundaries are decomposed into three disjoint parts ΓD, ΓN, and ΓH; i.e., ∂Ω=ΓD∪ΓN∪ΓH. ΓD is the Dirichlet boundary, ΓN is the nonhomogenous Neumann boundary, and ΓH is the homogenous Neumann boundary which is the only part to be optimized in present study.Introducing the virtual temperature fieldv and applying integration to (10), the weak form of state equation is given as(11)au,v,Φ=lv,Φ∀v∈UThe symmetric bilinear formau,v,Φ and the linear thermal load form lv,Φ are written as(12)au,v,Φ=∑i=1M∫Ωci∇uT∇vρiΦdΩ(13)lv,Φ=∑i=1M∫ΩpvρiΦdΩ+∫ΓhvdΓwhere p is the heat source per unit volume, ρiΦ is the characteristic function of the ith material, and h is the ambient temperature imposed on ΓN.


### 3.3. The Multimaterial Topology Optimization Model Using PLSM

In this paper, in order to minimize the quadratic temperature gradient subject to volume constraints of different materials, the multimaterial topology optimization model of heat conduction structures using PLSM is proposed. The optimization problem is mathematically expressed as follows:(14)find:αjk,k=1,2,…,M;j=1,2,…,Nmin:Ju,Φ=∫ΩCΦ∇u2dΩs.t.au,v,Φ=lv,ΦGkΦ≤fkV0,k=1,2,…,Mα_k≤αjk≤α¯k,j=1,2,…,Nwhere Ju,Φ denotes the objective function, u is the temperature field, and v is the virtual temperature field. V0 is the volume of design domain; GkΦ is the kth volume constraints; i.e., GkΦ=∫Ω∏i=1kHidΩ,k=1,2,…,M, and fk denotes the prescribed volume fraction according to kth volume constraints. α_k and α¯k are the lower and the upper bounds of the design variables in the kth level set function.In this topology optimization model, the volume constraints impose a restriction to each individual material, and the number of volume constraints is equal to the number of materials. For the volume constraints of three materials, volume constraints can be described by(15)G1Φ=∫ΩH1dΩ,G2Φ=∫ΩH1H2dΩ,G3Φ=∫ΩH1H2H3dΩ,where the first volume constraint restricts volume ratio of total materials, the second one refers to the sum of volumes for material two and three, and the third one denotes the usage of the material three.A smoothed approximation ofHΦk is used [7] to avoid the regenerating the element mesh during the updating of Φkx,t. The regularized Heaviside function and its derivatives are shown in Figure 2; it can be expressed as follows:(16)HΦk=χ,Φk≤-Δ31-χ4ΦkΔ-13ΦkΔ3+1+χ2,-Δ≤Φk≤Δ1,Φk≥Δwhere χ and γ are small positive numbers to avoid the singularity and Δ is the width of the numerical approximation.Figure 2
Regularized Heaviside function and its derivatives. (a) Heaviside function and its regularized Heaviside function; (b) the derivatives of regularized Heaviside function.
(a) (b)


### 3.4. Sensitivity Analysis

In order to obtain the structural optimal configuration, gradient based optimization algorithms are employed, such as method of moving asymptotes (MMA) and optimality criteria (OC). Thus, sensitivity analysis for the objective function and volume constraints with respect to the design variables are required. In this section, the adjoint method is used to obtain the sensitivities for multimaterial topology optimization problem of heat conduction structures using PLSM.The Lagrangian function is defined as(17)L=Ju,Φ+au,v,Φ-lv,Φ+∑k=1MλkGkΦ-fkV0where λk is Lagrangian multiplier for the kth volume constraint. The variation of the boundaries with respect to the pseudo-time t is considered. For the mth level set function, the shape derivative of the Lagrangian is given as follows:(18)∂L∂tΦm=∂Ju,Φ∂tΦm+∂au,v,Φ∂tΦm-∂lv,Φ∂tΦm+∑k=1MλkGkΦ∂tΦmwhere the derivatives (∂Ju,Φ/∂t)Φm, ∂au,v,Φ/∂tΦm, (∂lv,Φ/∂t)Φm, and (∂GkΦ/∂t)Φm are expressed by the following:(19)∂Ju,Φ∂tΦm=∑i=1M∫Ω2ciρiΦ∇u∇u˙dΩ+∑i=1M∫Ωci∇u2∂ρiΦ∂Φm∂Φm∂tdΩ(20)∂au,v,Φ∂tΦm=∑i=1M∫ΩciρiΦ∇uT∇v˙dΩ+∑i=1M∫ΩciρiΦ∇u˙T∇vdΩ+∑i=1M∫Ωci∇uT∇v∂ρiΦ∂Φm∂Φm∂tdΩ(21)∂lv,Φ∂tΦm=∑i=1M∫Ωpv˙ρiΦdΩ+∑i=1M∫Ωpv∂ρiΦ∂Φm∂Φm∂tdΩ+∫Γhv˙dΓ(22)GkΦ∂tΦm=∫Ω∏z=1,z≠mkHzδΦm∂Φm∂tdΩCollecting all the terms that include∇u˙ on the right side of (18) and making their sum to be zero, we can obtain the weak form of the state equation as follows:(23)∑i=1M∫Ω2ciρiΦ∇u∇u˙dΩ+∑i=1M∫ΩciρiΦ∇u˙T∇vdΩ=0Collecting all the terms that include∇v˙ on the right side of (18) and making their sum to be zero, we can obtain the weak form of the adjoint equation as follows:(24)∑i=1M∫ΩciρiΦ∇uT∇v˙dΩ-∑i=1M∫Ωpv˙ρiΦdΩ-∫Γhv˙dΓ=0The shape derivative of the Lagrangian can be obtained by substituting (19) to (24) into (18); that is,(25)∂L∂tΦm=∑i=1M∫Ωci∇u2∂ρiΦ∂Φm∂Φm∂tdΩ+∑i=1M∫Ωci∇uT∇v∂ρiΦ∂Φm∂Φm∂tdΩ+∑k=1Mλk∫Ω∏z=1,z≠mkHzδΦm∂Φm∂tdΩIt also can be written by substituting (5) into (25); that is,(26)∂L∂tΦm=∑j=1N∑i=1M∫Ωci∇u2+ci∇uT∇v∂ρiΦ∂Φmwjm∂αjmt∂tdΩ+∑j=1N∑k=1Mλk∫Ω∏z=1,z≠mkHzδΦmwjm∂αjmt∂tdΩThe derivative of the Lagrangian function with respect tot also can be obtained by using the chain rule in terms of the expansion coefficients as follows:(27)∂L∂tΦm=∑j=1N∂Ju,Φ∂αjmt+∑j=1N∑k=1Mλk∂GkΦ∂αjmt∂αjmt∂tComparing the corresponding terms in (26) and (28), the derivatives of the objective function and the volume constraints with respect to the expansion coefficients can be obtained by(28)∂Ju,Φ∂αjmt=∑i=1M∫Ωci∇u2+ci∇uT∇v∂ρiΦ∂ΦmwjmdΩ(29)∂GkΦ∂αjmt=∫Ω∏z=1,z≠mkHzδΦmwjmdΩ


## 3.1. Parametric Level Set

The PLSM with CSRBF which is a kind of radial basis function centered at a particular point with compact support is applied into multimaterial topology optimization of heat conduction structures. The level set surface can be approximated by CSRBF over the design domain. The level set function can be expressed as(5)Φkx,t=∑j=1Nwjxαjkt,k=1,2,…,Mwhere N is the total number of computational points. wj denotes the CSRBF at jth point. αjk is the expansion coefficient at jth point in kth level set function.In this work, CSRBF with C2 continuity [19, 20] is employed to approximate the level set function, which can be formulated as follows:(6)wjx=max⁡0,1-rjx44rjx+1j=1,2,…,Nwhere r is the radius of support, and it is usually defined in a 2D Euclidean space as(7)rjx=x-xjR,j=1,2,…Nwhere xj is the coordinate of the jth computational point and R is the radius of support domain of the computational point x, which is defined as 2.0-3.0.Substituting (5) into (2), the HJ-PDE in (2) can be rewritten as follows:(8)∑j=1Nwjxα˙jkt-vnk∑j=1N∇wjxαjkt=0,k=1,2,…,MThe normal velocity of the level set function can be calculated by(9)vnk=∑j=1Nwjxα˙jkt∑j=1N∇wjxαjktwhere α˙jk(t)=dαjk(t)/dt.


## 3.2. Formulations of Heat Conduction Problem

For a steady heat conduction optimization problem in two dimensions, the governing equation is given as(10)∇·C∇u+Q=0inΩu=0onΓDC∇u·n=0onΓNwhere C is the heat conduction coefficient, u is the temperature field, Q is the heat source, and n is the unit normal vector of boundary. Ω is a closed domain with boundaries; the boundaries are decomposed into three disjoint parts ΓD, ΓN, and ΓH; i.e., ∂Ω=ΓD∪ΓN∪ΓH. ΓD is the Dirichlet boundary, ΓN is the nonhomogenous Neumann boundary, and ΓH is the homogenous Neumann boundary which is the only part to be optimized in present study.Introducing the virtual temperature fieldv and applying integration to (10), the weak form of state equation is given as(11)au,v,Φ=lv,Φ∀v∈UThe symmetric bilinear formau,v,Φ and the linear thermal load form lv,Φ are written as(12)au,v,Φ=∑i=1M∫Ωci∇uT∇vρiΦdΩ(13)lv,Φ=∑i=1M∫ΩpvρiΦdΩ+∫ΓhvdΓwhere p is the heat source per unit volume, ρiΦ is the characteristic function of the ith material, and h is the ambient temperature imposed on ΓN.


## 3.3. The Multimaterial Topology Optimization Model Using PLSM

In this paper, in order to minimize the quadratic temperature gradient subject to volume constraints of different materials, the multimaterial topology optimization model of heat conduction structures using PLSM is proposed. The optimization problem is mathematically expressed as follows:(14)find:αjk,k=1,2,…,M;j=1,2,…,Nmin:Ju,Φ=∫ΩCΦ∇u2dΩs.t.au,v,Φ=lv,ΦGkΦ≤fkV0,k=1,2,…,Mα_k≤αjk≤α¯k,j=1,2,…,Nwhere Ju,Φ denotes the objective function, u is the temperature field, and v is the virtual temperature field. V0 is the volume of design domain; GkΦ is the kth volume constraints; i.e., GkΦ=∫Ω∏i=1kHidΩ,k=1,2,…,M, and fk denotes the prescribed volume fraction according to kth volume constraints. α_k and α¯k are the lower and the upper bounds of the design variables in the kth level set function.In this topology optimization model, the volume constraints impose a restriction to each individual material, and the number of volume constraints is equal to the number of materials. For the volume constraints of three materials, volume constraints can be described by(15)G1Φ=∫ΩH1dΩ,G2Φ=∫ΩH1H2dΩ,G3Φ=∫ΩH1H2H3dΩ,where the first volume constraint restricts volume ratio of total materials, the second one refers to the sum of volumes for material two and three, and the third one denotes the usage of the material three.A smoothed approximation ofHΦk is used [7] to avoid the regenerating the element mesh during the updating of Φkx,t. The regularized Heaviside function and its derivatives are shown in Figure 2; it can be expressed as follows:(16)HΦk=χ,Φk≤-Δ31-χ4ΦkΔ-13ΦkΔ3+1+χ2,-Δ≤Φk≤Δ1,Φk≥Δwhere χ and γ are small positive numbers to avoid the singularity and Δ is the width of the numerical approximation.Figure 2
Regularized Heaviside function and its derivatives. (a) Heaviside function and its regularized Heaviside function; (b) the derivatives of regularized Heaviside function.
(a) (b)


## 3.4. Sensitivity Analysis

In order to obtain the structural optimal configuration, gradient based optimization algorithms are employed, such as method of moving asymptotes (MMA) and optimality criteria (OC). Thus, sensitivity analysis for the objective function and volume constraints with respect to the design variables are required. In this section, the adjoint method is used to obtain the sensitivities for multimaterial topology optimization problem of heat conduction structures using PLSM.The Lagrangian function is defined as(17)L=Ju,Φ+au,v,Φ-lv,Φ+∑k=1MλkGkΦ-fkV0where λk is Lagrangian multiplier for the kth volume constraint. The variation of the boundaries with respect to the pseudo-time t is considered. For the mth level set function, the shape derivative of the Lagrangian is given as follows:(18)∂L∂tΦm=∂Ju,Φ∂tΦm+∂au,v,Φ∂tΦm-∂lv,Φ∂tΦm+∑k=1MλkGkΦ∂tΦmwhere the derivatives (∂Ju,Φ/∂t)Φm, ∂au,v,Φ/∂tΦm, (∂lv,Φ/∂t)Φm, and (∂GkΦ/∂t)Φm are expressed by the following:(19)∂Ju,Φ∂tΦm=∑i=1M∫Ω2ciρiΦ∇u∇u˙dΩ+∑i=1M∫Ωci∇u2∂ρiΦ∂Φm∂Φm∂tdΩ(20)∂au,v,Φ∂tΦm=∑i=1M∫ΩciρiΦ∇uT∇v˙dΩ+∑i=1M∫ΩciρiΦ∇u˙T∇vdΩ+∑i=1M∫Ωci∇uT∇v∂ρiΦ∂Φm∂Φm∂tdΩ(21)∂lv,Φ∂tΦm=∑i=1M∫Ωpv˙ρiΦdΩ+∑i=1M∫Ωpv∂ρiΦ∂Φm∂Φm∂tdΩ+∫Γhv˙dΓ(22)GkΦ∂tΦm=∫Ω∏z=1,z≠mkHzδΦm∂Φm∂tdΩCollecting all the terms that include∇u˙ on the right side of (18) and making their sum to be zero, we can obtain the weak form of the state equation as follows:(23)∑i=1M∫Ω2ciρiΦ∇u∇u˙dΩ+∑i=1M∫ΩciρiΦ∇u˙T∇vdΩ=0Collecting all the terms that include∇v˙ on the right side of (18) and making their sum to be zero, we can obtain the weak form of the adjoint equation as follows:(24)∑i=1M∫ΩciρiΦ∇uT∇v˙dΩ-∑i=1M∫Ωpv˙ρiΦdΩ-∫Γhv˙dΓ=0The shape derivative of the Lagrangian can be obtained by substituting (19) to (24) into (18); that is,(25)∂L∂tΦm=∑i=1M∫Ωci∇u2∂ρiΦ∂Φm∂Φm∂tdΩ+∑i=1M∫Ωci∇uT∇v∂ρiΦ∂Φm∂Φm∂tdΩ+∑k=1Mλk∫Ω∏z=1,z≠mkHzδΦm∂Φm∂tdΩIt also can be written by substituting (5) into (25); that is,(26)∂L∂tΦm=∑j=1N∑i=1M∫Ωci∇u2+ci∇uT∇v∂ρiΦ∂Φmwjm∂αjmt∂tdΩ+∑j=1N∑k=1Mλk∫Ω∏z=1,z≠mkHzδΦmwjm∂αjmt∂tdΩThe derivative of the Lagrangian function with respect tot also can be obtained by using the chain rule in terms of the expansion coefficients as follows:(27)∂L∂tΦm=∑j=1N∂Ju,Φ∂αjmt+∑j=1N∑k=1Mλk∂GkΦ∂αjmt∂αjmt∂tComparing the corresponding terms in (26) and (28), the derivatives of the objective function and the volume constraints with respect to the expansion coefficients can be obtained by(28)∂Ju,Φ∂αjmt=∑i=1M∫Ωci∇u2+ci∇uT∇v∂ρiΦ∂ΦmwjmdΩ(29)∂GkΦ∂αjmt=∫Ω∏z=1,z≠mkHzδΦmwjmdΩ


## 4. Optimization Algorithm

In this paper, the MMA is employed to solve the optimization problems. In this method, a disturbance for asymptote parameter is introduced to the optimization model, and then implicit optimization problem is converted into a series of explicit and strictly convex approximation optimization subproblems. The explicit subproblems can be solved by dual method or main dual method. As a well-founded mathematical programming algorithm, it has been widely used in structural topology optimization.The process of optimization is shown in Figure3; the corresponding introduction is presented as follows:Figure 3
Flowchart of multimaterial topology optimization.Step 1.
Initialize the design domain of multimaterial with PLSM. Define the boundary conditions, and set the volume constraints.Step 2.
Carry out the finite element method and calculate the temperature field by solving the heat conduction problems.Step 3.
Calculate the derivatives of objective function and the volume constraints with respect to the expansion coefficients in terms of (28) and (29).Step 4.
Calculate the values of objective function and the volume constraints. Update the expansion coefficients with MMA optimization algorithm, and then the updated level set function is obtained.Step 5.
Check whether the convergence condition shown in (30) is satisfied. If the terminal condition is met, the optimal structure is found. Otherwise, go back to Step 2. The terminal condition is defined as(30)Js-Js-1Js-1≤εwhere J(s) denotes the value of objective function at sth iteration, and ε=0.001.


## 5. Numerical Examples

In this section, numerical examples of two materials and three materials are given to demonstrate the effectiveness of the proposed method in this paper for multimaterial topology optimization of heat conduction structures. In all the examples, the design domain is discretized40×40 with linear quadrilateral elements for finite element analysis. The computational points in level set function are corresponding to the finite element nodes. The parameters in (14) are α_k=6×mini(αik)0 and α¯k=6×maxi(αik)0, where (αik)0 denotes the initial values for the kth level set function of expansion coefficients; the support radius for computational points is set to be R=2.5. Δ=1,χ=0.001 and γ=0.0005 are parameters given in (16) and (17).


### 5.1. Optimization for a Two Materials Problem

Example 1.
The design domain and boundary conditions are shown in Figure4. The heat supply per unit volume is set to be p=10; the values of the Dirichlet boundary conditions are given as u=0. The heat conductivities of material 1 and material 2 are c1=1 and c2=3, respectively. The volume fractions are given as f1=0.5 and f2=0.3, which means the maximum material volume of material 1 and material 2 is 0.2 and 0.3, respectively.
The optimization process of Example1 is shown in Figure 5. The blue regions represent the first material with heat conductivity c1, while the red regions represent the second material with heat conductivity c2. The initial configuration of multimaterial with a combination of two level set functions is shown in Figure 5(a). During the optimization process, each point in the design domain denotes one material, without overlaps between different materials. The configurations with smooth boundaries for different materials are achieved by solving two level set functions. The optimal configuration in Figure 5(f) shows that the second material with higher heat conductivity is distributed in the heat conduction path, and the first material with lower heat conductivity is distributed around the second material as an auxiliary phase. The iteration histories of objective function and volume ratios are given in Figure 6. It is found that the objective function decreases gradually until converging at the value of 42.58 after 85 iterations, which shows that the proposed method has a higher computational efficiency than the conventional level set method. And the volume constraints can be well handled.Figure 4
The design domain and boundary conditions for Example1.Figure 5
The optimization process of Example1. (a) Initial configuration; (b-e) intermediate configurations; (f) optimal configuration.
(a) (b) (c) (d) (e) (f)Figure 6
Iteration histories of objective function and volume ratios of Example1.Example 2.
The design domain and boundary conditions are shown in Figure7. The heat supply per unit volume is set to be p=2; the values of the Dirichlet boundary conditions are given as u=0. The other parameters are the same as those of Example 1.
The optimization process of Example2 is shown in Figure 8. The blue regions represent the first material with heat conductivity c1, while the red regions represent the second material with heat conductivity c2. The initial configuration of multimaterial with a combination of two level set functions is shown in Figure 8(a). During the optimization process, each point in the design domain denotes one material, without overlaps between different materials. The configurations with smooth boundaries for different materials are achieved by solving two level set functions. The optimal configuration in Figure 8(f) shows that the second material with higher heat conductivity is distributed in the heat conduction path, and the first material with lower heat conductivity is distributed around the second material as an auxiliary phase. The iteration histories of objective function and volume ratios are given in Figure 9. It is found that the objective function decreases gradually until convergence at the value of 27.35 after 94 iterations, which shows that the proposed method has a higher computational efficiency than the conventional level set method. And the volume constraints can be well handled.Figure 7
The design domain and boundary conditions for Example2.Figure 8
The optimization process of Example2. (a) Initial configuration; (b-e) intermediate configurations; (f) optimal configuration.
(a) (b) (c) (d) (e) (f)Figure 9
Iteration histories of objective function and volume ratios of Example2.


### 5.2. Optimization for a Three Materials Problem

Example 3.
The design domain and boundary conditions are shown in Figure10. The heat supply per unit volume is set to be p=2; the values of the Dirichlet boundary conditions are given as u=0. The heat conductivities of three materials are c1=1, c2=3, and c3=5, respectively. The volume fractions are given as f1=0.5, f2=0.3, and f3=0.15, which means the maximum material volume of three materials is 0.2, 0.15, and 0.15, respectively.
The optimization process of Example3 is shown in Figure 11. The yellow regions, blue regions, and red regions denote the first material with heat conductivity c1, the second material with heat conductivity c2, and the third material with heat conductivity c3, respectively. Figure 11(a) shows the initialization of multimaterial with a combination of three level set functions. During the optimization process, the configurations with clear boundaries for different materials are achieved by updating the expansion coefficients, and the smooth boundaries can be well retained. The optimal configuration in Figure 11(f) shows that the third material with the highest heat conductivity is distributed in regions of the heat supply and the heat conduction path, the second material with the higher heat conductivity is distributed around the third material, and the first material with the lower heat conductivity is distributed in the heat conduction path as an auxiliary phase. The iteration histories of objective function and volume ratios are given in Figure 12. It is found that the objective function increases gradually until convergence at the value of 30.45 after 185 iterations. The volumes of materials decrease to the volume constrains. And the volume constraints can be well handled.
From these examples we can observe that the optimal configurations obtained by the proposed method are essentially identical to the optimal configurations obtained by the conventional method [29]; the optimal results show the effectiveness of the proposed method. Comparing with the conventional method, the proposed method possesses the following advantages: (1) it naturally has a smoother boundary as LSFs are interpolated by CSRBFs; (2) the sensitivities are obtained easily as the explicit mathematical expression of MM-LS model and CSRBFs; (3) without solving the HJ-PDE equation, the computational efficiency is improved.Figure 10
The design domain and boundary conditions for Example3.Figure 11
The optimization process of Example3. (a) Initial configuration; (b-e) intermediate configurations; (f) optimal configuration.
(a) (b) (c) (d) (e) (f)Figure 12
Iteration histories of objective function and volume ratios of Example3.


## 5.1. Optimization for a Two Materials Problem

Example 1.
The design domain and boundary conditions are shown in Figure4. The heat supply per unit volume is set to be p=10; the values of the Dirichlet boundary conditions are given as u=0. The heat conductivities of material 1 and material 2 are c1=1 and c2=3, respectively. The volume fractions are given as f1=0.5 and f2=0.3, which means the maximum material volume of material 1 and material 2 is 0.2 and 0.3, respectively.
The optimization process of Example1 is shown in Figure 5. The blue regions represent the first material with heat conductivity c1, while the red regions represent the second material with heat conductivity c2. The initial configuration of multimaterial with a combination of two level set functions is shown in Figure 5(a). During the optimization process, each point in the design domain denotes one material, without overlaps between different materials. The configurations with smooth boundaries for different materials are achieved by solving two level set functions. The optimal configuration in Figure 5(f) shows that the second material with higher heat conductivity is distributed in the heat conduction path, and the first material with lower heat conductivity is distributed around the second material as an auxiliary phase. The iteration histories of objective function and volume ratios are given in Figure 6. It is found that the objective function decreases gradually until converging at the value of 42.58 after 85 iterations, which shows that the proposed method has a higher computational efficiency than the conventional level set method. And the volume constraints can be well handled.Figure 4
The design domain and boundary conditions for Example1.Figure 5
The optimization process of Example1. (a) Initial configuration; (b-e) intermediate configurations; (f) optimal configuration.
(a) (b) (c) (d) (e) (f)Figure 6
Iteration histories of objective function and volume ratios of Example1.Example 2.
The design domain and boundary conditions are shown in Figure7. The heat supply per unit volume is set to be p=2; the values of the Dirichlet boundary conditions are given as u=0. The other parameters are the same as those of Example 1.
The optimization process of Example2 is shown in Figure 8. The blue regions represent the first material with heat conductivity c1, while the red regions represent the second material with heat conductivity c2. The initial configuration of multimaterial with a combination of two level set functions is shown in Figure 8(a). During the optimization process, each point in the design domain denotes one material, without overlaps between different materials. The configurations with smooth boundaries for different materials are achieved by solving two level set functions. The optimal configuration in Figure 8(f) shows that the second material with higher heat conductivity is distributed in the heat conduction path, and the first material with lower heat conductivity is distributed around the second material as an auxiliary phase. The iteration histories of objective function and volume ratios are given in Figure 9. It is found that the objective function decreases gradually until convergence at the value of 27.35 after 94 iterations, which shows that the proposed method has a higher computational efficiency than the conventional level set method. And the volume constraints can be well handled.Figure 7
The design domain and boundary conditions for Example2.Figure 8
The optimization process of Example2. (a) Initial configuration; (b-e) intermediate configurations; (f) optimal configuration.
(a) (b) (c) (d) (e) (f)Figure 9
Iteration histories of objective function and volume ratios of Example2.


## 5.2. Optimization for a Three Materials Problem

Example 3.
The design domain and boundary conditions are shown in Figure10. The heat supply per unit volume is set to be p=2; the values of the Dirichlet boundary conditions are given as u=0. The heat conductivities of three materials are c1=1, c2=3, and c3=5, respectively. The volume fractions are given as f1=0.5, f2=0.3, and f3=0.15, which means the maximum material volume of three materials is 0.2, 0.15, and 0.15, respectively.
The optimization process of Example3 is shown in Figure 11. The yellow regions, blue regions, and red regions denote the first material with heat conductivity c1, the second material with heat conductivity c2, and the third material with heat conductivity c3, respectively. Figure 11(a) shows the initialization of multimaterial with a combination of three level set functions. During the optimization process, the configurations with clear boundaries for different materials are achieved by updating the expansion coefficients, and the smooth boundaries can be well retained. The optimal configuration in Figure 11(f) shows that the third material with the highest heat conductivity is distributed in regions of the heat supply and the heat conduction path, the second material with the higher heat conductivity is distributed around the third material, and the first material with the lower heat conductivity is distributed in the heat conduction path as an auxiliary phase. The iteration histories of objective function and volume ratios are given in Figure 12. It is found that the objective function increases gradually until convergence at the value of 30.45 after 185 iterations. The volumes of materials decrease to the volume constrains. And the volume constraints can be well handled.
From these examples we can observe that the optimal configurations obtained by the proposed method are essentially identical to the optimal configurations obtained by the conventional method [29]; the optimal results show the effectiveness of the proposed method. Comparing with the conventional method, the proposed method possesses the following advantages: (1) it naturally has a smoother boundary as LSFs are interpolated by CSRBFs; (2) the sensitivities are obtained easily as the explicit mathematical expression of MM-LS model and CSRBFs; (3) without solving the HJ-PDE equation, the computational efficiency is improved.Figure 10
The design domain and boundary conditions for Example3.Figure 11
The optimization process of Example3. (a) Initial configuration; (b-e) intermediate configurations; (f) optimal configuration.
(a) (b) (c) (d) (e) (f)Figure 12
Iteration histories of objective function and volume ratios of Example3.


## 6. Conclusions

In this paper, we have presented PLSM to solve the multimaterial topology optimization of heat conduction structures with volume constraints. The combination of different level set functions is used to implicitly represent the boundaries for different materials. The sensitivities can be obtained easily by MM-LS model, which has an explicit mathematical formation, and there is no overlap among different materials during the optimization process. The boundaries evolve by updating the expansion coefficients of PLSM, in which the CSRBFs with C2 smoothness are utilized. In this paper, the optimization model is developed by PLSM, the derivatives of objective function and volume constraints with expansion coefficients are derived, and the expansion coefficients are updated by MMA.Numerical examples illustrate the effectiveness of the proposed method for multimaterial topology optimization of heat conduction structures. Smooth geometric boundaries for different materials are obtained by combining level set functions. The problem of transient heat topology optimization for minimizing quadratic temperature gradient with multimaterial will be considered in further research.



---

*Source: 9804123-2018-09-16.xml*