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

Ellie Vineyard and Xin-Lin Gao

1PepsiCo R&D Center,Valhalla,NY 10595,USA

2Department of Mechanical Engineering,Southern Methodist University,Dallas,TX 75275-0337,USA

ABSTRACT 2-D and 3-D micro-architectured multiphase thermoelastic metamaterials are designed and analyzed using a parametric level set method for topology optimization and the finite element method.An asymptotic homogenization approach is employed to obtain the effective thermoelastic properties of the multiphase metamaterials.The ε-constraint multi-objective optimization method is adopted in the formulation.The coefficient of thermal expansion(CTE)and Poisson’s ratio(PR)are chosen as two objective functions,with the CTE optimized and the PR treated as a constraint.The optimization problems are solved by using the method of moving asymptotes.Effective isotropic and anisotropic CTEs and stiffness constants are obtained for the topologically optimized metamaterials with prescribed values of PR under the constraints of specified effective bulk modulus, volume fractions and material symmetry.Two solid materials along with one additional void phase are involved in each of the 2-D and 3-D optimal design examples.The numerical results reveal that the newly proposed approach can integrate shape and topology optimizations and lead to optimal microstructures with distinct topological boundaries.The current method can topologically optimize metamaterials with a positive,negative or zero CTE and a positive,negative or zero Poisson’s ratio.

KEYWORDS Topology optimization; thermoelastic metamaterial; level set method; sensitivity analysis; Poisson’s ratio;coefficient of thermal expansion; effective elastic properties

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

wherex,y)(i=1, 2, 3) denote the homogeneous parts, andx,y)(i=1, 2, 3;k=1, 2, 3,...)represent local variations at the scale of heterogeneities.Note thatx,y)is periodic in y (also known asY-periodic), i.e.,1+n1Y1,x2+n2Y2,x3+n3Y3)=1,x2,x3)for any integersn1,n2andn3and periodsY1,Y2andY3.

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

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

whereΩis the domain occupied by the material,Yrepresents 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])

whereΩεis the material domain with microstructures,Γtis the traction-prescribed part of the smooth closed boundary surface ofΩε,Sεis the interface (or hole surface),Cijklare the components of the elasticity (stiffness) tensor,viare the components of the virtual displacement,biare the components of the body force,tiare the components of the traction on the boundary partΓt,piare the components of the traction on the interfaceSε, andΔTis 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]):

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

wherestands for three or six linearly independent unit test strain fields for the 2-D or 3-D case, andandrepresent the locally varying strain and thermal strain fields induced,respectively, by the test strain modes and thermal load, which can be obtained from

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

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]

for 2-D cases, where

for 3-D cases, where

In Eqs.(16a), (16b), (17a) and (17b),represents the volume average of the property for the multi-phase composite, withvfnbeing the volume fraction of each phase andNthe number of phases,μminandμmaxare, respectively, the minimum and maximum shear moduli,andare the lower and upper Hashin-Shtrikman bounds on the 2-D bulk modulus given by [20,38]

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

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 ofmlevel 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

where the coefficients(k=1, 2;i=1, 2,...,n) serve as design variables in optimization, withirepresenting the number of knot points in the design domain andkdenoting two level set functions employed in the optimization,αjis the coefficient of thermal expansion serving as the objective function, withjbeing 1 and 2 for the 2-D case or 1, 2 and 3 for the 3-D case,andandandare, 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 CHdue 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,andare, respectively, the lower and upper bounds of each design variableto ensure a converged and stable solution, u and v are, respectively, the actual and virtual displacement fields,U is the kinematically admissible displacement space, u0is the prescribed displacement on the admissible Dirichlet boundaryΓD, andaandlare the bilinear energy forms given by

in which Cijkl(x,φ)is the locally varying elastic stiffness tensor that is related to the level set functionφ, andis 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

whereφ1andφ2represent 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, andH(φi) (i=1, 2) is the Heaviside function defined by

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αjand constraint functions (including Poisson’s ratioνand volume fractionsVf1andVf2) need to be computed with respect to the design variablesHence, the derivatives of the effective elastic stiffness tensorand the effective thermal stress tensorwith respect toare derived here.

From Eq.(10), the derivative ofwith respect to the pseudo-timetcan be obtained as

From Eq.(14), it follows that

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

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

From Eq.(15), it follows that

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

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

where the Dirac delta functionδφdis the derivative of the Heaviside functionH(φd), i.e.,δφd=∂H(φd)/∂φd(no sum ond).

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-timet, the Hamilton-Jacobi partial differential equation (HJ-PDE) can be readily obtained as

wherek= 1 or 2 (with no sum onk),is the normal velocity of thekth level set function given by (e.g., [51])

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

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

Alternatively, the derivative ofcan also be expressed as, upon using the chain rule,

By comparing Eqs.(34) and (35), the derivative ofwith respect tocan be obtained as

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

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

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

whered=1, 2 andb=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])

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

whereVf1andVf2are, respectively, the volume fractions of the two solid materials 1 and 2.

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

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 Ceand coefficient of thermal expansion tensorαeof theeth element can be computed from Eqs.(21) and (22) through integration over the element domain as

whereVeis the volume (or area in a 2-D case) of theeth element.

When the unit cell is discretized intoNEfinite elements, the effective (homogenized) stiffness and thermal stress tensors given in Eqs.(10) and (11) can be written as

where

is the element stiffness matrix, B(x)is the strain-displacement matrix,denotes the element displacement under the applied strain,represents the global displacement vector associated with the elementeunder the mechanical loading obtainable from solving Eq.(14), andstands for the element displacement underαpqand includes the components of the global displacement vector under thermal loading that can be computed using Eq.(15).

Further, to avoid the numerical singularity, the Heaviside functionH(φi) and the Dirac delta functionδφkinvolved in Eqs.(34)-(45) are approximated by the following piece-wise smooth functions [56]:

whereΔis a real positive constant that is usually chosen as 2～4 times of the background cell size, andξandγare small positive numbers to be selected to ensure that the overall stiffness of the structure is non-singular.In this study,Δis taken to be 3.5 times of the average grid knot spacing, andξandγare selected to be 0.0001 and 0.0005, respectively.

The method of moving asymptotes (MMA) is a good optimizer in tackling multiple constraints and normally can achieve a quick convergence [52,53].To stabilize the convergence of the optimization algorithm, the moving limit in this algorithm can be adjusted flexibly in the range of 0.05 to 0.1.In addition, the lower and upper bounds of the design variables are chosen to beandin whichcdenotes the iteration cycle,ζis a controlling parameter ranging from 1.5 to 2, accommodating different optimization problems.The results are considered convergent when the difference between two successive iterations for the objective function is less than 0.1%.

6.1 2-D Examples

To demonstrate the procedure of topology optimization of multiphase thermoelastic metamaterials, two sets of examples are presented in this sub-section under the plane stress conditions:one is to minimize anisotropic CTE, and the other is to minimize isotropic CTE, both with a Poisson’s ratio constraint.In all of the numerical examples, the constraint of the effective bulk modulusKHis enforced to ensure the structural load-carrying capability.KHis a measurement of resistance to the volumetric strain, which, for the 2-D case, can be expressed in terms of the effective elasticity tensor components as (e.g., [48])

In all the optimizations demonstrated below, the square design domain is discretized into 60 by 60 square four-node elements, in which four Gaussian quadrature points are distributed in each direction, resulting in a total of 16 integration points in each element.For illustration purposes,material properties used in the simulations here are taken to beE1=1 GPa,ν1=0.3,α1=1×10-51/K for solid phase 1 (in blue) andE2=1 GPa,ν2=0.3,α2=10×10-51/K for solid phase 2 (in pink).The volume fractions for the two solid materials are specified as 25% and 10%, respectively,with a tolerance of ±2%.Note that the elastic properties of the two constituent materials are chosen to be the same here to signify the effects of geometry and shape changes in the selected example problems.However, since the optimization method and the algorithm proposed here are mathematical in nature, they can be readily applied to other scenarios, where the two constituent materials have distinct (including high contrast) elastic properties, by simply changing values of the material parameters in the input file.This, of course, will lead to different optimal microstructures for the two-phase metamaterials, depending on the specified constituent properties.

Fig.2 shows the two initial guesses employed in the topology optimization of the 2-D composites with different material distributions.

6.1.1 Minimum Anisotropic CTE with a Poisson’s Ratio Constraint

The optimization herein is to minimize the CTE in the horizontal directionwith the constraint of a positive, a zero and a negative Poisson’s ratio, respectively.The optimization starts from the two initial material distributions illustrated in Fig.2.The constraint of 20% of the upper bound on the effective H-S bulk modulus predicted using Eq.(18b) is also enforced along with the horizontal and vertical geometrical symmetry.

Optimal microstructures with 3-by-3 unit cells (with the single unit cell denoted by the solid green line) from the initial distribution # 1 shown in Fig.2a are presented in Figs.3a, 3c and 3e with the constraint of Poisson’s ratio being 0.5, 0 and -0.5, respectively.Moreover, the final optimized microstructures from the initial distribution # 2 in Fig.2b are displayed in Figs.4a, 4c and 4e.The plots in the right column of Figs.3 and 4 show the iteration history curves of the objective function, Poisson’s ratioν, and two volume fractionsVf1andVf2, respectively.The main results are summarized in Tab.1, which contains the total number of iterations, minimumeffective elastic stiffness tensor and thermal strain tensor for all the optimization examples.It gives a side-by-side comparison of the optimization results based on the two initial distributions.It can be seen that for both of the initial distributions in Fig.2, with the decrease of the constrained Poisson’s ratio from 0.5 to -0.5, the absolute value of the attainable minimum CTE in thexdirection decreases under the same bulk modulus constraint.

Clearly, the final optimization results based on the two initial guesses differ from each other but do not display huge differences.This is consistent with the well-known fact that final optimal designs based on the conventional level set methods tend to strongly depend on initial guesses,since the Hamilton-Jacobi PDE (see Eq.(31)) satisfies a maximum principle which prohibits the hole nucleation inside the material domain (e.g., [60]).That is, without hole nucleation mechanism in the conventional level set methods, the final optimal design will largely depend on the initial guess which dictates the maximum number of holes allowed.In order to mitigate this inherent drawback, one common and simple approach is to choose an initial guess that contains sufficiently many small holes densely distributed in the design domain, which will be allowed to merge and evolve gradually in the optimization process (e.g., [60,61]).However, this simple approach has its limitations, and more advanced techniques are needed to reduce the dependence of the level set methods on initial guesses.One of such techniques is based on the use of multiquadric radial basis functions (RBFs) to construct implicit level set functions that do not require reinitialization [60].RBFs have been utilized in the studies on the topology optimization of 2-D mechanical metamaterials based on a parametric level set method [51] and on homogenization analyses of 3D printable interpenetrating phase composites by a meshfress method [62].

Moreover, it is observed from the optimal microstructures shown in the left column of Fig.3 that different Poisson’s ratio constraints result in different topological features, but similarities exist due to the same initial material distribution.A comparison of the microstructures displayed in Figs.3a, 3c and 3e shows that they all possess vertical struts with the higher CTE (in pink)that mainly control the CTE in the horizontal direction.However, the unit cell with the optimal minimum CTE under the constraint of PR = 0.5 shown in Fig.3a has struts with the lower CTE (in blue) along the vertical direction due to the positive Poisson’s ratio constraint.For the optimal microstructure under the constraint of PR = 0 displayed in Fig.3c, the solid material with the lower CTE (in blue) is mainly concentrated at the center of the unit cell, while most of the vertical struts comprise of the material with the higher CTE (in pink).In addition, for the optimal microstructure with the constraint of PR =-0.5 shown in Fig.3e, an anti-chiral structure, which is known to be a main mechanism responsible for the negative Poisson’s ratio, is clearly seen.Furthermore, through comparing all the optimized results displayed in Figs.3 and 4 based on the two initial distributions # 1 and # 2 shown in Fig.2, it is found that the optimized microstructures with the constraint of PR = 0.5 look almost the same.However, for the other two cases with PR = 0 and PR =-0.5, the final optimal configurations resulting from the two different initial distributions display distinct topologies.

Figure 2: Initial distributions of materials in a 2-D two-phase composite: (a) case # 1, (b) case # 2

Figure 3: Minimization of with a Poisson’s ratio constraint starting from the initial distribution# 1 shown in Fig.2a: (a) Three-by-three unit cells with PR = 0.5; (b) Iteration history with PR= 0.5 (totaling 108 iterations); (c) Three-by-three unit cells with PR = 0; (d) Iteration history with PR = 0 (totaling 190 iterations); (e) Three-by-three unit cells with PR = -0.5; (f) Iteration history with PR = -0.5 (totaling 149 iterations)

Figure 4: Minimization of with a Poisson’s ratio constraint starting from the initial distribution# 2 shown in Fig.2b: (a) Three-by-three unit cells with PR = 0.5; (b) Iteration history with PR= 0.5 (totaling 65 iterations); (c) Three-by-three unit cells with PR = 0; (d) Iteration history with PR = 0 (totaling 139 iterations); (e) Three-by-three unit cells with PR = -0.5; (f) Iteration history with PR = -0.5 (totaling 177 iterations)

Clearly, Tab.1 shows that the optimized 2-D composite in each case exhibits anisotropic CTEs and the cubic (for PR = 0.5 and -0.5) or isotropic (for PR = 0) elastic symmetry.

Table 1: Minimum α11 with the Poisson’s ratio constraint based on the two initial distributions in Fig.2

6.1.2 Minimum Isotropic CTE with a Poisson’s Ratio Constraint

The goal of this sub-section is to design microstructures with minimized isotropic CTEs and desirable Poisson’s ratios for the given horizontal and vertical symmetry.All the microstructures optimized herein possess cubic symmetry (e.g., [12]).

Figs.5 and 6 show the optimal microstructures with 3-by-3 unit cells (with the single unit cell denoted by the solid green line) under different bulk modulus constraints, i.e., 15%, 20%, 25%,30%, 35% and 40% of the H-R upper bound on the effective bulk modulus given by Eq.(18b),and the zero Poisson’s ratio constraint.The plots displayed in the right column of Figs.5 and 6 show the iteration history curves of the CTE, Poisson’s ratio and two volume fractions based on the initial material distribution # 1 displayed in Fig.2a.

It is clearly observed from Figs.5 and 6 that the optimized microstructures are different from each other, even though there is a detectable similarity due to the same initial distribution.Also,the microstructures with the 15%, 20% and 30% of the H-R upper bound of the effective bulk modulus displayed in Figs.5a, 5c and 6a exhibit almost the same topological features but with different unit cell sizes.In addition, for the microstructures with the 25%, 35% and 40% of the H-R upper bound of the effective bulk modulus constraints displayed in Figs.5e, 6c and 6e, they look slightly topologically similar but still show a significant difference.The main findings from the optimization results shown in Figs.5 and 6 are summarized in Tab.2, which gives a clear picture of how the effective thermoelastic properties of the optimized microstructure change with the bulk modulus constraint.It is observed that the obtained minimum isotropic CTE increases with the increase of the bulk modulus, which agrees with what is analytically predicted using Eq.(16a).

Figure 5: Minimization of isotropic CTE with the constraints of Poisson’s ratio = 0 and specified value of the bulk modulus based on the initial distribution # 1 in Fig.2: (a), (b) Three-by-three unit cells with the bulk modulus = 15% and the iteration history (118 iterations); (c), (d) Threeby-three unit cells with the bulk modulus = 20% and the iteration history (189 iterations); (e),(f) Three-by-three unit cells with the bulk modulus = 25% and the iteration history (176 iterations)

Figure 6: Minimization of isotropic CTE with the constraints of Poisson’s ratio = 0 and specified value of the bulk modulus based on the initial distribution # 1 in Fig.2: (a), (b) Three-by-three unit cells with the bulk modulus = 30% and the iteration history (117 iterations); (c), (d) Threeby-three unit cells with the bulk modulus = 35% and the iteration history (249 iterations); (e),(f) Three-by-three unit cells with the bulk modulus = 40% and the iteration history (302 iterations)

Clearly, Tab.2 shows that the topologically optimized 2-D three-phase composite in each case exhibits isotropic CTEs and the cubic or isotropic symmetry in the elastic stiffness constants.

Table 2: Minimum CTE with the constraints of different bulk modulus values and zero Poisson’s ratio from the initial distribution # 1 in Fig.2a

In addition to the topology optimization for minimizing CTE with the constraint of Poisson’s ratio = 0, simulations with the constraint of Poisson’s ratio = -0.5 are performed with the constraints of 15%, 20% and 25% of the H-R upper bounds on the effective bulk modulus given by Eq.(18b), respectively.

Figs.7 and 8 show the optimal results obtained from the two initial distributions displayed in Fig.2.In all six numerical examples, the constraint of Poisson’s ratio =-0.5 is well maintained.Also, it is seen that the two different initial distributions lead to different microstructures with distinct minimized isotropic CTEs.For the 15% bulk modulus constraint, a minimum of -1.442 is obtained from the initial distribution # 1 and of -1.959 from the initial distribution # 2.Further,the optimal microstructures from the initial distribution # 1 contain an anti-chiral unit that is mainly responsible for the negative Poisson’s ratio, while those from the initial configuration #2 possess a re-entrant structure, leading to a negative Poisson’s ratio of -0.5.The chirality and re-entrant structures are two well-known deformation mechanisms that can bring about negative Poisson’s ratios.Such observations are also true for the optimization with the 20% and 25% bulk modulus constraints, respectively.

Figure 7: Minimization of isotropic CTE with the constraints of Poisson’s ratio = -0.5 and specified value of the bulk modulus based on the initial distribution # 1 displayed in Fig.2:(a), (b) Three-by-three unit cells with the bulk modulus = 15% and the iteration history (179 iterations); (c), (d) Three-by-three unit cells with the bulk modulus = 20% and the iteration history(157 iterations); (e), (f) Three-by-three unit cells with the bulk modulus = 25% and the iteration history (241 iterations)

Figure 8: Minimization of isotropic CTE with the constraints of Poisson’s ratio = -0.5 and specified value of the bulk modulus based on the initial distribution # 2 shown in Fig.2: (a), (b)Three-by-three unit cells with the bulk modulus = 15% and the iteration history (168 iterations);(c), (d) Three-by-three unit cells with the bulk modulus = 20% and the iteration history (155 iterations); (e), (f) Three-by-three unit cells with the bulk modulus = 25% and the iteration history(182 iterations)

The optimization results displayed in Figs.7 and 8 are summarized in Tab.3, which includes the total number of iterations, the effective elastic stiffness tensor, and the coefficient of thermal expansion tensor.It is seen that the two different initial configurations lead to different optimization results, but no huge derivation exists.Moreover, the variation trend of the minimized CTE with respect to the bulk modulus constraint is the same: it increases with the increase of the bulk modulus in both cases.

Clearly, Tab.3 shows that the topologically optimized 2-D three-phase composite in each case exhibits isotropic CTEs and the cubic symmetry in the elastic stiffness constants.

Table 3: Minimum CTE with the constraints of different bulk modulus values and Poisson’s ratio= -0.5 based on the two initial distributions displayed in Fig.2

The effective properties of the 2-D composites with the optimal topologies shown in Figs.5-8 are displayed as three types of markers in Fig.9, in which the upper and lower bounds on the effective isotropic CTEs are plotted against the effective bulk modulus obtained from Eqs.(16a)and (16b).A total of 12 examples are included with the final optimal unit cells labeled from 1 to 12.As displayed at the right upper corner of the graph, all the numerical values of the effective CTE are close to the lower bound.However, through a close-up observation from Fig.9, it is seen that the effective CTE values for the optimized microstructures with the constraint of Poisson’s ratio = 0 are much closer to the lower bound than those with the constraint of Poisson’s ratio=-0.5.

Figure 9: Effective isotropic CTE obtained from the topology optimization compared with the upper and lower bounds of the effective CTE based on Eqs.(16a) and (16b) for the 12 sample cases with the indicated optimal unit cells initially shown in Figs.5-8

The results from Sections 6.1.1 and 6.1.2 reveal that Poisson’s ratio constraints can largely affect the extreme values of the effective CTE of a metamaterial with an optimized microstructure.The constraint of a positive Poisson’s ratio can lead to a negative CTE with a much larger magnitude than that of a negative Poisson’s ratio can.Therefore, optimal designs need to be performed to obtain metamaterials with desired negative CTEs and negative Poisson’s ratios.

6.2 3-D Examples

Similar to the 2-D multi-material and multi-objective optimizations described in Section 6.1,3-D optimizations are performed in this sub-section.The design domain consists of 24×24×24 8-node brick elements with three Gaussian quadrature points distributed in each direction, resulting in a total of 27 integration points in each element.The initial 3-D design with two solid materials and one void phase is displayed in Fig.10.For simplicity, the properties are taken to beE1=1 GPa,ν1=0.3,α1=1×10-51/K for material 1 (in blue), andE2=1 GPa,ν2=0.3,α2=10×10-51/K for material 2 (in pink).To demonstrate the design procedures for topology optimization of 3-D multi-phase thermoelastic metamaterials, two examples are presented here: one is to minimize anisotropic CTE, and the other is to minimize isotropic CTE.Both cases are constrained with a prescribed value of the effective Poisson’s ratio.In all the numerical examples shown below,constraints on the volume fractions are enforced along with material cubic symmetry conditions.Eqs.(46) and (47) are used to obtain the optimal effective elastic stiffness and CTE tensors for the 3-D examples presented herein.

Figure 10: Initial material distribution for the 3-D case

6.2.1 Minimum Anisotropic CTE with a Poisson’s Ratio Constraint

The first case is to obtain 3-D microstructures with the minimumunder the constraint of Poisson’s ratio being 0.1, 0 or -0.1.The volume fractions for the two solid materials are constrained to be 15% and 10%, respectively, with a tolerance of ±2%.

Figs.11a, 12a and 13a display the final optimal unit cells showing distributions of two solid phases.The 2×2×2 unit cells are shown in Figs.11b, 12b and 13b.In addition, Figs.11c-11e,12c-12e and 13c-13e display microstructural details from thex-y,x-zandy-zcross section views, respectively.

As shown in Figs.11f-11h, the convergence is achieved after 116 iterations, giving a minimumof -3.4584 and a Poisson’s ratio of 0.10177 that is very close to the constrained value of 0.1 in the topology optimization problem.Furthermore, the two volume fraction constraints are well maintained: one is 0.1499, and the other is 0.099, which are almost the same as the pre-defined volume fractions of 0.15 and 0.1.

Similarly, Figs.12f-12h show that the convergence is reached after 96 iterations, yielding a minimumof -2.7479 and a Poisson’s ratio of 0.000599, which is very close to the constrained value of zero Poisson’s ratio in the topology optimization.Moreover, the two volume fraction constraints are well maintained: one is 0.14997, and the other is 0.090172, which are almost the same as the pre-defined volume fractions of 0.15 and 0.1, with a tolerance of ±2%.

For the optimization with the constraint of Poisson’s ratio = -0.1, it can be observed from Figs.13f-13h that the convergence is achieved after 195 iterations, giving a minimumof-2.9405 and a Poisson’s ratio of -0.10213 which is quite close to the constrained value of Poisson’s ratio = -0.1 in the topology optimization.Also, the two volume fraction constraints are well maintained: one is 0.14795, and the other is 0.080033, which are very close to the pre-defined volume fractions of 0.15 and 0.1, with a tolerance of ±2%.

Note that the magnitude of= -3.4584 obtained with the constraint of Poisson’s ratio= 0.1 is larger than that of=-2.7479 under the constraint of Poisson’s ratio = 0, which is smaller than the magnitude of=-2.9405 obtained with Poisson’s ratio = -0.1.This indicates that the Poisson’s ratio and CTE are not directly related and can therefore be independently varied without affecting each another, which agrees with what was found in [19].

Figure 11: Minimization of under the constraint of Poisson’s ratio = 0.1: (a) Optimal unit cell; (b) Periodic microstructures with 2×2×2 unit cells; (c) Unit cell in the x-y cross section view; (d) Unit cell in the x-z cross section view; (e) Unit cell in the y-z cross section view; (f)Iteration history of CTE in the y-direction; (g) Iteration history of Poisson’s ratio; (h) Iteration history of two volume fraction constraints (total 116 iterations to achieve convergence)

Figure 12: Minimization of under the constraint of Poisson’s ratio = 0: (a) Optimal unit cell;(b) Periodic microstructures with 2×2×2 unit cells; (c) Unit cell in the x-y cross section view;(d) Unit cell in the x-z cross section view; (e) Unit cell in the y-z cross section view; (f) Iteration history of CTE in the y-direction; (g) Iteration history of Poisson’s ratio; (h) Iteration history of two volume fraction constraints (total 96 iterations to achieve convergence)

Figure 13: Minimization of under the constraint of Poisson’s ratio =-0.1: (a) Optimal unit cell; (b) Periodic microstructures with 2×2×2 unit cells; (c) Unit cell in the x-y cross section view; (d) Unit cell in the x-z cross section view; (e) Unit cell in the y-z cross section view; (f)Iteration history of CTE in the y-direction; (g) Iteration history of Poisson’s ratio; (h) Iteration history of two volume fraction constraints (total 195 iterations to achieve convergence)

Figs.11-13 show 3-D optimal microstructures with complex topologies achieved through the optimization evolution process originated from the initial design illustrated in Fig.10.This indicates that the approach employed in the current study is capable of merging and creating holes without implementing topological derivatives [63] or other techniques.Also, it is seen that the present approach can integrate shape and topology optimizations, in which the first few iterations are used to quickly achieve the optimal topology and the rest of iterations before the final convergence are mainly serving to complete the shape optimization.

The effective elastic stiffness and CTE tensors for the three optimal microstructures with different Poisson’s ratio constraints shown in Figs.11-13 are, respectively, obtained as

Clearly, Eqs.(52)-(54) show that the topologically optimized 3-D three-phase composite in each case exhibits anisotropic CTEs and elastic stiffness constants that satisfy the cubic symmetry.

It should be pointed out that the optimal geometrical boundaries shown in Figs.11-13 appear to have serrated shapes that are not sufficiently smooth.The reason for this is that the total number of elements used here is 13824 in a 24×24×24 mesh, which is not fine enough to produce smooth boundaries but remains a limitation of the current study for the 3-D threephase metamaterials.In addition, it can be observed that the same initial structures could lead to distinctive designs if they have been given different constraints even with the same objective function.Alternatively, the different initial distributions lead to different final optimal microstructures owing to the highly non-convex problems with multiple local regions existing in the current gradient-based optimization enabled by the PLSM (e.g., [60]).

From the unit cells displayed in Figs.11-13, it is not easy to tell underlying mechanisms that may lead to the desired negative CTE and Poisson’s ratio, unlike that for the 2-D cases described in Section 6.1.This is a feature inherent in topology optimization that could result in designs unattainable through conventional thinking.

6.2.2 Minimum Isotropic CTE with a Poisson’s Ratio Constraint

The goal here is to obtain the minimum isotropic CTE under the constraints of Poisson’s ratio being 0 and the volume fractions of the two solid materials at 15% and 6%, respectively,with a tolerance of ±2%.

Fig.14a displays the optimal unit cell, and Fig.14b shows 2×2×2 periodic unit cells.It can be seen that very clear and distinct material interfaces are formed between the two solid phases that can be reconstructed with ease from the mesh generated in topology optimization.In addition, Figs.14c-14e provide the 2-D cross-sectional views from three perspectives that further demonstrate fine geometrical features obtained from topology optimization.Furthermore,Figs.14f-14h give the convergence plots for the effective CTE, Poisson’s ratio, and two volume fraction constraints, in which 210 iterations are performed to achieve the final results with the CTE being -0.69308, Poisson’s ratio -0.0411 and two volume fractions 0.13005 and 0.04, respectively.For the solid material with the higher CTE displayed in pink, the volume fraction first jumps to about 90% and then gradually drops down to satisfy the volume fraction constraint in the topology optimization.On the other hand, for the solid phase with the lower CTE (marked in blue), it drops from the initial volume fraction of 48% down to near zero and then increases to achieve the required volume fraction.Moreover, the optimization approaches adopted here work pretty well concerning the hole nucleation and merging.

Figure 14: Minimization of isotropic CTE under the zero Poisson’s ratio constraint: (a) Optimal unit cell; (b) Periodic microstructures with 2×2×2 unit cells; (c) Unit cell in the x-y cross section view; (d) Unit cell in the x-z cross section view; (e) Unit cell in the y-z cross section view; (f) Iteration history of CTE; (g) Iteration history of Poisson’s ratio; (h) Iteration history of the two volume fraction constraints (totalling 210 iterations to achieve convergence)

The effective elastic stiffness and CTE tensors for the optimal microstructure generated in this case are obtained as

From Eq.(55), it is seen that the topologically optimized 3-D three-phase composite in this case shows an isotropic CTE and elastic stiffness constants satisfying the cubic symmetry.

Finally, it should be mentioned that fabricating the 3-D multi-phase metamaterial samples discussed here by using traditional subtractive manufacturing techniques is very challenging.However, additive manufacturing methods can be employed to print these optimally designed composite structures (e.g., [64]).Topologically optimized 2-D multiphase polymeric metamaterials with a tunable CTE have been produced by additive manufacturing through an Objet Connex 500 3D printer (e.g., [25]).

7 Conclusion

A parametric level set-based topology optimization method is used to design 2-D and 3-D multiphase thermoelastic metamaterials with minimized anisotropic and isotropic CTEs and prescribed values of Poisson’s ratio under the constraints of specified effective bulk modulus, volume fractions and material symmetry.The effective properties of the multiphase metamaterials are predicted using an asymptotic homogenization method.Theε-constraint multi-objective optimization method is utilized in the formulation, and the method of moving asymptotes is employed to solve the optimization problems.

In each design example considered, two solid materials and one void phase are included.It is found that optimal 2-D and 3-D microstructures can be generated by applying the newly proposed approach.Novel structures are obtained through the current topology optimization approach, which are difficult to achieve using conventional (heuristic) methods.In addition, the new approach can lead to topologically optimized metamaterials with a positive, negative or zero Poisson’s ratio and a positive, negative or zero CTE.

Acknowledgement:The authors gratefully thank Prof.Krister Svanberg for providing his MMA codes to us.We also would like to thank two anonymous reviewers for their encouragement and helpful comments on an earlier version of the paper.

Funding Statement:The authors received no specific funding for this study.

Conflicts of Interest:The authors declare that they have no conflicts of interest to report regarding the present study.