Computer Modeling in Engineering & Sciences

DOI: 10.32604/cmes.2021.015688

ARTICLE

Topology and Shape Optimization of 2-D and 3-D Micro-Architectured Thermoelastic Metamaterials Using a Parametric Level Set Method

Ellie Vineyard1 and Xin-Lin Gao2,*

1PepsiCo R&D Center, Valhalla, NY 10595, USA
2Department of Mechanical Engineering, Southern Methodist University, Dallas, TX 75275-0337, USA
*Corresponding Author: Xin-Lin Gao. Email: xlgao@smu.edu
Received: 05 January 2021; Accepted: 24 February 2021

1  Introduction

Micro-architectured thermoelastic metamaterials are a new class of materials with unusual thermal and elastic properties such as a negative Poisson’s ratio and a non-positive coefficient of thermal expansion (CTE) (e.g., [1–7]). For such metamaterials, micro-architectures play a crucial role in attaining targeted or extremal properties beyond those of their constituents. It offers additional degrees of freedom to achieve exotic properties that are not exhibited by naturally occurring or conventionally designed materials.

In developing these metamaterials, heuristic approaches have typically been used (e.g., [1,8–14]). However, such approaches are limited to simple geometrical designs or loading conditions. For complex configurations and deformation mechanisms, topology optimization has emerged as a promising method (e.g., [15–18]). Sigmund et al. [19,20] optimally designed materials with a zero or negative CTE based on the solid isotropic material with penalization (SIMP) method. They found that materials with a negative CTE can be obtained by using two solid materials (each having a positive CTE) and one void phase. Schwerdtfeger et al. [21] optimized a 3-D structure with a negative Poisson’s ratio and increased the negativity of Poisson’s ratio by a factor 2 using a SIMP based inverse homogenization approach. Andreassen et al. [22] designed a manufacturable 3-D elastic microstructure with a Poisson’s ratio of −0.5 through topology optimization with a manufacturing constraint. Wang et al. [23] proposed multi-phase metamaterials with unusual thermoelastic properties by using a parametric level set-based topology optimization method combined with a finite element approach. Vogiatzis et al. [24] employed a reconciled level set-based topology optimization method to design single- and multiple-phase metamaterials with negative Poisson’s ratios in both 2-D and 3-D configurations. Takezawa et al. [25] developed a topology optimization method for porous composites with anisotropic negative CTEs and isotropic/anisotropic positive CTEs. Wang et al. [26] designed multi-phase and multi-functional metamaterials with targeted effective elastic moduli and CTEs and proposed some periodic microstructures that can produce negative Poisson’s ratios and negative CTEs. Ye et al. [27] developed an optimization framework for gradually stiffer mechanical metamaterials with a negative Poisson’s ratio using a parametric level set method and a numerical homogenization approach. Li et al. [28] proposed a robust topology optimization model for thermoelastic properties of multiphase metamaterials by considering hybrid interval-random uncertainties in properties of constituent materials. However, these authors did not consider 3-D cases or anisotropic CTEs under constraints of Poisson’s ratio varying from negative to positive values. This motivated the current study.

In the present paper, topology optimization of micro-architectured multiphase thermoelastic metamaterials is conducted using a parametric level set method combined with a finite element analysis. The ε-constraint multi-objective optimization method is adopted. An asymptotic homogenization approach is employed to obtain the effective thermoelastic properties of the metamaterials, and the method of moving asymptotes is applied to solve the optimization problem. Both 2-D and 3-D design examples are included to illustrate the new approach, which differs from what was done in [20,26] and other existing studies on topology optimization of multiphase metamaterials with extreme thermoelastic properties, where only 2D configurations were considered in their design examples. In Section 2, the asymptotic homogenization method for predicting effective thermoelastic properties of composites is briefly introduced. Existing analytical formulas for bounds on the effective CTE of heterogeneous composites are reviewed in Section 3. In Section 4, the topology optimization based on the parametric level set method (PLSM) is formulated, which is followed by a sensitivity analysis in Section 5. Numerical results for 2-D and 3-D design examples are presented in Section 6. The paper concludes in Section 7 with a summary.

2  Asymptotic Homogenization of Thermoelastic Properties

Asymptotic homogenization is a widely used approach in which two spatial scales (i.e., microscopic and macroscopic) are considered (e.g., [29–31]). The coordinate system at the microscopic scale is y=(y1,y2,y3), while that at the macroscopic scale is x=(x1,x2,x3). The two scales are linked through y=x/ε, where ε ( > 0) is a small dimensionless parameter.

Using the double-scale asymptotic expansion, the displacement field in a periodic heterogeneous (composite) material can be expressed as

uiε(x,y)=ui0(x,y)+εui1(x,y)+ε2ui2(x,y)+⋅⋅⋅,(1)

where ui0(x,y) (i = 1, 2, 3) denote the homogeneous parts, and uik(x,y) (i = 1, 2, 3; k = 1, 2, 3,…) represent local variations at the scale of heterogeneities. Note that uik(x,y) is periodic in y (also known as Y-periodic), i.e., uik(x1+n1Y1,x2+n2Y2,x3+n3Y3)=uik(x1,x2,x3) for any integers n1, n2 and n3 and periods Y1, Y2 and Y3.

It then follows from Eq. (1) that the derivative of the displacement field uiε(x,y)with respect to the macroscopic coordinate x can be obtained as, with the help of the chain rule and the relation y=x/ε,

∂uiε(x,y)∂xj=∂ui0(x,y)∂xj+1ε∂ui0(x,y)∂yj+ε∂ui1(x,y)∂xj+∂ui1(x,y)∂yj+ε2∂ui2(x,y)∂xj+ε∂ui2(x,y)∂yj+⋅⋅⋅.(2)

In addition, the limit of the integral of a Y-periodic function Ψ(y) as ε approaches zero is given by (e.g., [29])

limε→0∫ΩΨ(xε)dΩ=∫Ω[1|Y|∫YΨ(y)dY]dΩ,(3)

where Ω is the domain occupied by the material, Y represents a unit cell of the material, and |Y| is the unit cell area (in a 2-D case) or volume (in a 3-D case).

According to the principle of virtual work, the equilibrium equation of a material undergoing thermoelastic deformations can be written as (e.g., [29,32])

∫ΩεCijkl∂ukε∂xl∂vi∂xjdΩ=∫ΩεbividΩ+∫ΓttividΓ+∫SεpividS+∫ΩεCijklαklΔT∂vi∂xjdΩ,(4)

where Ωε is the material domain with microstructures, Γt is the traction-prescribed part of the smooth closed boundary surface of Ωε, Sε is the interface (or hole surface), Cijkl are the components of the elasticity (stiffness) tensor, vi are the components of the virtual displacement, bi are the components of the body force, ti are the components of the traction on the boundary part Γt, pi are the components of the traction on the interface Sε, and ΔT is the temperature change.

Using Eqs. (2) and (3) in Eq. (4) gives the following three hierarchical equations based on the order of ε (e.g., [32]):

CijklH=1|Y|∫Y(Cijkl-Cijmn∂χmkl∂yn)dY,(5)

βijH=1|Y|∫YCijmn(αmnΔT-∂ϒm∂yn)dY,(6)

αijH=(CijklH)-1βklH,(7)

where CijklH, βijH and αijH are, respectively, the effective (homogenized) elastic stiffness, thermal stress and coefficient of thermal expansion tensors, χmkl and ϒm are, respectively, the characteristic displacement fields induced by the mechanical and thermal loading, which can be obtained from (e.g., [20])

∫YCijmn∂χmkl∂yn∂vi∂yjdY=∫YCijmnεmn0(kl)∂vi∂yjdY,(8)

∫YCijkl∂ϒk∂yl∂vi∂yjdY=∫YCijklαkl∂vi∂yjdY.(9)

Note that Eqs. (5) and (6) can each be rewritten in a symmetric form as (e.g., [20,33])

CijklH=1|Y|∫Y(εpq0(kl)-εpq*(kl))Cpqrs(εrs0(ij)-εrs*(ij))dY,(10)

βklH=1|Y|∫Y(αpqΔT-εpqα)Cpqrs(εrs0(kl)-εrs*(kl))dY,(11)

where εrs0(ij) stands for three or six linearly independent unit test strain fields for the 2-D or 3-D case, and εpq*(kl) and εpqα represent the locally varying strain and thermal strain fields induced, respectively, by the test strain modes and thermal load, which can be obtained from

εpq(kl)=εpq(χkl)=12(∂χpkl∂yq+∂χqkl∂yp),(12)

εpqα=εpq(ϒα)=12(∂ϒpα∂yq+∂ϒqα∂yp).(13)

Also, Eqs. (8) and (9) can be rewritten as, with the help of Eqs. (12) and (13) and the symmetry of Cijkl,

∫YCpqrsεpq*(χkl)∂vr∂ysdY=∫YCpqrsεpq0(χkl)∂vr∂ysdY,(14)

∫YCpqrsεpq(ϒα)∂vr∂ysdY=∫YCpqrsαpq∂vr∂ysdY.(15)

3  Bounds on the Effective Coefficient of Thermal Expansion

Several analytical and semi-analytical formulas have been provided to evaluate the effective CTE of heterogeneous composites, which include those reported in [34–37]. The model by Gibiansky et al. [37], which gives tight and sharp bounds compared to those presented in Schapery [35] and Rosen et al. [36], is adopted in this study. The upper and lower bounds on the effective CTE read [20,37]

α2DUα2DL}=1κ*(κ2DU-κ2DL)[(κ2DU-κ*)(κ2DL+μmin)〈κακ+μmin〉+(κ*-κ2DL)(κ2DU+μmax)〈κακ+μmax〉±Ψ2D1/2(κ2DU-κ*)1/2(κ*-κ2DL)1/2](16a)

for 2-D cases, where

Ψ2D=-(κ2DU+μmax)(κ2DL+μmin){〈κακ+μmin〉-〈κακ+μmax〉}2+(μmax-μmin){(κ2DL+μmin)〈κακ+μmin〉2-(κ2DU+μmax)〈κακ+μmax〉2}+(κ2DU-κ2DL){μmax〈κα2κ+μmax〉-μmin〈κα2κ+μmin〉},(16b)

and

α3DUα3DL}=1κ*(κ3DU-κ3DL)[(κ3DU-κ*)(κ3DL+4μmin3)〈3κα3κ+4μmin〉+(κ*-κ3DL)(κ3DU+4μmax3)〈3κα3κ+4μmax〉±Ψ3D1/2(κ3DU-κ*)1/2(κ*-κ3DL)1/2](17a)

for 3-D cases, where

Ψ3D=-(κ3DU+4μmax3)(κ3DL+4μmin3)[〈3κα3κ+4μmin〉-〈3κα3κ+4μmax〉]2+4(μmax-μmin)3[(κ3DL+4μmin3)〈3κα3κ+4μmin〉2-(κ3DU+4μmax3)〈3κα3κ+4μmax〉2]+(κ3DU-κ3DL)[4μmax〈κα23κ+4μmax〉-4μmin〈κα23κ+4μmin〉].(17b)

In Eqs. (16a), (16b), (17a) and (17b), 〈X〉=∑n=1NvfnXn represents the volume average of the property for the multi-phase composite, with vfn being the volume fraction of each phase and N the number of phases, μmin and μmax are, respectively, the minimum and maximum shear moduli, κ2DL and κ2DU are the lower and upper Hashin-Shtrikman bounds on the 2-D bulk modulus given by [20,38]

κ2DL=〈1κ+μmin〉-1-μmin,(18a)

κ2DU=〈1κ+μmax〉-1-μmax,(18b)

and κ3DL and κ3DU are, respectively, the lower and upper Hashin-Shtrikman bounds on the 3-D bulk modulus given by [38,39]

κ3DL=〈33κ+4μmin〉-1-4μmin3,(19a)

κ3DU=〈33κ+4μmax〉-1-4μmax3,(19b)

where κ=E/[2(1-ν)] for the 2-D plane stress case and κ=E/[3(1-2ν)] for the 3-D case, and μ=E/[2(1+ν)] for both the 2-D and 3-D cases (e.g., [40]).

4  Topology Optimization Using the PLSM

The parametric level set method (PLSM) developed in [41,42] is a powerful shape and topology optimization approach that can overcome some shortcomings of the conventional level set method, such as no regularization, velocity extension, and numerical time step limitations (e.g., [43]). In addition, for multi-material optimizations, Wang et al. [44] proposed the use of m level set functions to conduct multi-phase structural optimizations on the basis of the PLSM, which can prevent overlapping between different phases and suppress redundant regions in the design domain when compared to other multi-material approaches (e.g., [45,46]).

In order to obtain multi-functional composites with optimal properties, multi-objective topology optimization methods need to be used. A few algorithms have been developed to achieve multi-objective optimization. The weighted sum optimization method is a classical approach that combines all the objectives into a single objective function, where the designer prescribes the weights a priori. Different weights can be assigned to different terms in the objective function based on their relative importance. These weights influence the final design. A major disadvantage of the weighted sum optimization is its strong dependence on the weights that needs to be carefully tailored according to specific applications (e.g., [47,48]). The ε-constraint method is another popular method, in which one of the objective functions is optimized, while the other objective functions are regarded as constraints that can limit the influence from the assigned weights. A comparison between these two methods can be found in [49]. In addition, evolutionary algorithms such as the genetic algorithm and particle swarm optimization have been widely used (e.g., [50]). However, with a large number of design variables involved in a multi-objective topology optimization problem, it is quite challenging to use evolutionary algorithms due to the lack of speedy convergence and consistency.

Sigmund et al. [19] showed that there is no direct relationship between a negative coefficient of thermal expansion (CTE) and a negative Poisson’s ratio (PR). CTEs and PRs are not competing properties. Therefore, each of them can change without affecting the other.

In this study, the ε-constraint multi-objective optimization approach is adopted. The CTE and PR are chosen as two objective functions, but only the CTE is optimized and the PR is treated as a constraint. This differs from that in [26], where the weighted sum optimization method was employed by combining the thermal strain and elasticity tensors into one objective function using two sets of weighted factors.

The current topology optimization algorithm can be described as

{Findmik(i=1,2,…,n;k=1,2)Minimize:αjSubject to:V1L≤Vf1≤V1U,V2L≤Vf2≤V2U,g(CH)=0,ν≤ν*,a(u,v,ϕ)=l(v,ϕ),u|ΓD=u0,∀v∈U,mikL≤mik≤mikU,(20a)

where the coefficients mik (k = 1, 2; i = 1, 2,…, n) serve as design variables in optimization, with i representing the number of knot points in the design domain and k denoting two level set functions employed in the optimization, αj is the coefficient of thermal expansion serving as the objective function, with j being 1 and 2 for the 2-D case or 1, 2 and 3 for the 3-D case, V1L and V1U and V2L and V2U are, respectively, the prescribed lower and upper bounds on the volume fraction for solid materials 1 and 2, g(CH) denotes the constraint on the stiffness tensor CH due to the material symmetry, ν and ν* are, respectively, Poisson’s ratio and its prescribed upper bound to be used in the constraint to obtain a positive, negative or a near zero Poisson’s ratio, mikL and mikU are, respectively, the lower and upper bounds of each design variable mik to ensure a converged and stable solution, u and v are, respectively, the actual and virtual displacement fields, U is the kinematically admissible displacement space, u0 is the prescribed displacement on the admissible Dirichlet boundary ΓD, and a and l are the bilinear energy forms given by

a(u,v,ϕ)=∫Yεij(u)Cijkl(x,ϕ)εkl(v)dY,(20b)

l(v,ϕ)=∫Yεij0(u)Cijkl(x,ϕ)εkl(v)dY,(20c)

in which Cijkl(x,ϕ) is the locally varying elastic stiffness tensor that is related to the level set function ϕ, and εij0 is the applied strain.

The parametric level set-based topology optimization method (e.g., [41,51]) is adopted in this study. The method of moving asymptotes (MMA) developed by Svanberg [52,53] is used to solve the optimization problem defined in Eq. (20a). The MMA is a well-known optimizer that is very efficient for structural optimization by minimizing sequential convex approximations of the original function. It is also widely employed to handle multiple constraints. These are different from what was done in [20], where the sequential linear programming method was employed as the optimization method.

According to the multi-material parametric level set method (e.g., [44,54]), the elastic stiffness tensor C and the coefficient of thermal expansion tensor α at an arbitrary point x can be determined as

C(x,ϕ)=H(ϕ1)[(1-H(ϕ2))C1+H(ϕ2)C2],(21)

α(x,ϕ)=(1-H(ϕ2))α1+H(ϕ2)α2,(22)

where ϕ1 and ϕ2 represent two level set functions employed in the multi-material optimization, and Cδ and αδ (δ=1, 2) are, respectively, the elastic stiffness tensors and coefficient of thermal expansion tensors for the two phases, and H(ϕi) (i = 1, 2) is the Heaviside function defined by

H(ϕi(x))={0,ϕi(x)<0,1,ϕi(x)≥0.(23)

Note that the cubic and isotropic elastic symmetries are considered in the microstructure design in the current study. To ensure the cubic symmetry, the constraints on the elastic constants given by g(CH) = 0 in Eq. (20a) have the explicit expressions (e.g., [55]): C1111=C2222=C3333, C1122=C1133=C2233, C1123=C1112=C1113=C2223=C2212=C2213=C3323=C3312=C3313=C2312=C2313=C1213=0, and C2323=C1212=C1313. To attain the isotropic symmetry, one additional constraint in the form of C2323=(C1111-C1122)/2 will be imposed along with those for the cubic symmetry (e.g., [12]). The above-mentioned constraints can be readily implemented in the MMA optimization algorithm, which is capable of handling multiple-constraint problems.

5  Sensitivity Analysis

In order to employ the MMA, the first derivatives of the objective function αj and constraint functions (including Poisson’s ratio ν and volume fractions Vf1 and Vf2) need to be computed with respect to the design variables mik. Hence, the derivatives of the effective elastic stiffness tensor CijklH and the effective thermal stress tensor βklH with respect to mik are derived here.

From Eq. (10), the derivative of CijklH with respect to the pseudo-time t can be obtained as

∂CijklH∂t=1|Y|∫Y(εrs0(ij)-εrs*(ij))∂Cpqrs(x,ϕ)∂t(εpq0(kl)-εpq*(kl))dY+21|Y|∫Y(εpq0(kl)-εpq*(kl))Cpqrs(x,ϕ)∂(εrs0(ij)-εrs*(ij))∂tdY.(24)

From Eq. (14), it follows that

∫Y[∂Cpqrs∂t(εpq0(kl)-εpq*(kl))+Cpqrs∂(εpq0(kl)-εpq*(kl))∂t]∂vr∂ysdY=0.(25)

Using Eq. (25) in Eq. (24) leads to

∂CijklH∂t=-1|Y|∫Y(εrs0(ij)-εrs*(ij))∂Cpqrs(x,ϕ)∂t(εpq0(kl)-εpq*(kl))dY.(26)

Similarly, the derivative of the effective thermal stress tensor in Eq. (11) with respect to the pseudo-time t can be determined as

∂βklH∂t=1|Y|∫Y(αpq-εpqα)∂Cpqrs(x,ϕ)∂t(εrs0(kl)-εrs*(kl))dY+1|Y|∫Y∂(αpq-εpqα)∂tCpqrs(x,ϕ)(εrs0(kl)-εrs*(kl))dY+1|Y|∫Y(αpq-εpqα)Cpqrs(x,ϕ)∂(εrs0(kl)-εrs*(kl))∂tdY.(27)

From Eq. (15), it follows that

∫Y[∂Cpqrs∂t(αpq-εpqα)+Cpqrs∂(αpq-εpqα)∂t]∂vr∂ysdY=0.(28)

Substituting Eq. (28) into Eq. (27) yields, with the help of Eq. (26),

∂βklH∂t=-1|Y|∫Y(αpqΔT-εpqα)∂Cpqrs(x,ϕ)∂t(εrs0(kl)-εrs*(kl))dY.(29)

To obtain the derivative of a function with respect to the design variables mik, Eq. (26) can be rewritten as, with the help of Eq. (21),

∂CijklH∂t=-1|Y|∫Y(εrs0(ij)-εrs*(ij))∑d=12[∂Cpqrs(x,ϕ)∂H(ϕd)δϕd∂ϕd∂t](εpq0(kl)-εpq*(kl))dY,(30)

where the Dirac delta function δϕd is the derivative of the Heaviside function H(ϕd), i.e., δϕd=∂H(ϕd)/∂ϕd (no sum on d).

In a level set-based topology optimization method, the structural design boundary is implicitly represented by the zero-level set of a one-dimensional-higher level set function with the Lipschitz continuity (e.g., [56]). Through differentiating the zero-level set equation ϕk(x, t) = 0 with respect to the pseudo-time t, the Hamilton-Jacobi partial differential equation (HJ-PDE) can be readily obtained as

∂ϕk(x,t)∂t-vnk|∇ϕk|=0,(31)

where k = 1 or 2 (with no sum on k), vnk is the normal velocity of the kth level set function given by (e.g., [51])

vnk=[φk(x)]Tṁk(t)|(∇φk)Tmk(t)|,(32)

in which φk(x) is a n×1 column matrix representing the n compactly supported radial basis functions φik (i = 1, 2,…, n) defined over the fixed grid points, and mk(t) is a n×1 column matrix representing n generalized expansion coefficients mik (i = 1, 2,…, n) serving as design variables. Note that in reaching Eq. (32) use has been made of the relation

ϕk(x,t)=φk(x)Tmk(t).(33)

Substituting Eqs. (31) and (32) into Eq. (30) gives, with the help of Eq. (33),

∂CijklH∂t=-1|Y|∫Y(εrs0(ij)-εrs*(ij))[∑d=12∂Cpqrs(x,ϕ)∂H(ϕd)δϕdφd(x)Tṁd(t)](εpq0(kl)-εpq*(kl))dY=-1|Y|∫Y(εrs0(ij)-εrs*(ij))[∑d=12∑b=1n∂Cpqrs(x,ϕ)∂H(ϕd)δϕdφbd(x)∂mbd(t)∂t](εpq0(kl)-εpq*(kl))dY.(34)

Alternatively, the derivative of CijklH can also be expressed as, upon using the chain rule,

∂CijklH∂t=∑d=12∑b=1n∂CijklH∂mbd∂mbd(t)∂t.(35)

By comparing Eqs. (34) and (35), the derivative of CijklH with respect to mbd can be obtained as

∂CijklH∂mbd=-1|Y|∫Y(εrs0(ij)-εrs*(ij))∂Cpqrs(x,ϕ)∂H(ϕd)δϕdφbd(x)(εpq0(kl)-εpq*(kl))dY.(36)

For the case with two level set functions (e.g., [26]), the two derivatives can be obtained from Eqs. (36) and (21) as

∂CijklH∂mb1=-1|Y|∫Y(εpq0(ij)-εpq*(ij))[(1-H(ϕ2))Cpqrs1+H(ϕ2)Cpqrs2](εrs0(kl)-εrs*(kl))φb1(x)δϕ1dY,(37)

∂CijklH∂mb2=-1|Y|∫Y(εpq0(ij)-εpq*(ij))H(ϕ1)(Cpqrs2-Cpqrs1)(εrs0(kl)-εrs*(kl))φb2(x)δϕ2dY.(38)

Similarly, the derivatives of the effective thermal stress tensor can be obtained from Eqs. (29) and (21) as

∂βklH∂mb1=-1|Y|∫Y(αpqΔT-εpqα)[(1-H(ϕ2))Cpqrs1+H(ϕ2)Cpqrs2](εrs0(kl)-εrs*(kl))φb1(x)δϕ1dY,(39)

∂βklH∂mb2=-1|Y|∫Y(αpqΔT-εpqα)H(ϕ1)(Cpqrs2-Cpqrs1)(εrs0(kl)-εrs*(kl))φb2(x)δϕ2dY.(40)

For the effective CTE tensor, its derivative with respect to the design variables mbd can be obtained from Eq. (7) as

∂αijH∂mbd=(CijklH)-1∂βklH∂mbd+∂(CijklH)-1∂mbdβklH,(41)

where d = 1, 2 and b = 1, 2,…, n.

For a multi-material system with two solid phases and one void phase, as shown in Fig. 1, the volume fractions of the two solid materials can be written as (e.g., [51])

Vf1=1|Y|∫YH(ϕ1)[1-H(ϕ2)]dY,(42a)

Vf2=1|Y|∫YH(ϕ1)H(ϕ2)dY,(42b)

where Vf1 and Vf2 are, respectively, the volume fractions of the two solid materials 1 and 2.

Figure 1: Multi-material system with two solid phases and one void phase via two level set functions

Then, the sensitivity of each volume fraction with respect to the design variables mbd can be obtained from Eqs. (42a) and (42b) as, with the help of Eq. (33),

∂Vf1∂mi1=1|Y|∫Yφi1(x)(1-H(ϕ2))δϕ1dY,(43a)

∂Vf1∂mi2=-1|Y|∫Yφi2(x)H(ϕ1)δϕ2dY,(43b)

∂Vf2∂mi1=1|Y|∫YH(ϕ2)φi1(x)δϕ1dY,(43c)

∂Vf2∂mi2=1|Y|∫YH(ϕ1)φi2(x)δϕ2dY.(43d)

6  Numerical Results and Discussion

The linear thermoelasticity equations given in Eqs. (14) and (15) are solved using the finite element method (FEM). A fixed and rectilinear mesh is used in the homogenization-based topology optimization (e.g., [57]), and ANSYS [58] is employed as the computational tool. In addition, the level set knots are positioned at the finite element nodes for simplicity. In the examples included here, an ‘ersatz’ material model (e.g., [59]) is used to approximate material properties for those elements crossed by the moving level set boundaries, i.e., the zero level set. The volume-averaged elastic stiffness tensor Ce and coefficient of thermal expansion tensor αe of the eth element can be computed from Eqs. (21) and (22) through integration over the element domain as

Ce=1Ve∫Ωe[H(ϕ1)(1-H(ϕ2))C1+H(ϕ2)C2]dΩ,(44)