|Computer Modeling in |
Engineering & Sciences
Topology and Shape Optimization of 2-D and 3-D Micro-Architectured Thermoelastic Metamaterials Using a Parametric Level Set Method
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: firstname.lastname@example.org
Received: 05 January 2021; Accepted: 24 February 2021
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
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.  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.  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.  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.  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.  developed a topology optimization method for porous composites with anisotropic negative CTEs and isotropic/anisotropic positive CTEs. Wang et al.  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.  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.  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 , while that at the macroscopic scale is . The two scales are linked through , 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
where (i = 1, 2, 3) denote the homogeneous parts, and (i = 1, 2, 3; k = 1, 2, 3,) represent local variations at the scale of heterogeneities. Note that is periodic in y (also known as Y-periodic), i.e., 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 with respect to the macroscopic coordinate x can be obtained as, with the help of the chain rule and the relation ,
In addition, the limit of the integral of a Y-periodic function as approaches zero is given by (e.g., )
where is the domain occupied by the material, Y represents a unit cell of the material, and 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, is the traction-prescribed part of the smooth closed boundary surface of , 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 , pi are the components of the traction on the interface , and 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., ):
where , and are, respectively, the effective (homogenized) elastic stiffness, thermal stress and coefficient of thermal expansion tensors, and are, respectively, the characteristic displacement fields induced by the mechanical and thermal loading, which can be obtained from (e.g., )
Note that Eqs. (5) and (6) can each be rewritten in a symmetric form as (e.g., [20,33])
where stands for three or six linearly independent unit test strain fields for the 2-D or 3-D case, and and represent 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 of Cijkl,
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. , which gives tight and sharp bounds compared to those presented in Schapery  and Rosen et al. , 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, with vfn being the volume fraction of each phase and N the number of phases, and are, respectively, the minimum and maximum shear moduli, and are the lower and upper Hashin-Shtrikman bounds on the 2-D bulk modulus given by [20,38]
and and are, respectively, the lower and upper Hashin-Shtrikman bounds on the 3-D bulk modulus given by [38,39]
where for the 2-D plane stress case and for the 3-D case, and for both the 2-D and 3-D cases (e.g., ).
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., ). In addition, for multi-material optimizations, Wang et al.  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 . In addition, evolutionary algorithms such as the genetic algorithm and particle swarm optimization have been widely used (e.g., ). 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.  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 , 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, with i representing the number of knot points in the design domain and k denoting two level set functions employed in the optimization, 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, and and and 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, and are, respectively, the lower and upper bounds of each design variable 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 , and a and l are the bilinear energy forms given by
in which is the locally varying elastic stiffness tensor that is related to the level set function , and 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 , 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 and represent two level set functions employed in the multi-material optimization, and C and (, 2) are, respectively, the elastic stiffness tensors and coefficient of thermal expansion tensors for the two phases, and H() (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., ): C, C, C, and C. To attain the isotropic symmetry, one additional constraint in the form of C will be imposed along with those for the cubic symmetry (e.g., ). 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 and constraint functions (including Poisson’s ratio and volume fractions Vf1 and Vf2) need to be computed with respect to the design variables . Hence, the derivatives of the effective elastic stiffness tensor and the effective thermal stress tensor with respect to are derived here.
From Eq. (10), the derivative of with respect to the pseudo-time t can 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-time t can 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 variables , Eq. (26) can be rewritten as, with the help of Eq. (21),
where the Dirac delta function is the derivative of the Heaviside function H(), i.e., ()/ (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., ). Through differentiating the zero-level set equation (x, t) = 0 with respect to the pseudo-time t, the Hamilton-Jacobi partial differential equation (HJ-PDE) can be readily obtained as
where k = 1 or 2 (with no sum on k), is the normal velocity of the kth level set function given by (e.g., )
in which (x) is a column matrix representing the n compactly supported radial basis functions (i = 1, 2,, n) defined over the fixed grid points, and mk(t) is a column matrix representing n generalized 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 of can also be expressed as, upon using the chain rule,
By comparing Eqs. (34) and (35), the derivative of with respect to can be obtained as
For the case with two level set functions (e.g., ), 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 variables can be obtained from Eq. (7) as
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., )
where Vf1 and Vf2 are, respectively, the volume fractions of the two solid materials 1 and 2.
Then, the sensitivity of each volume fraction with respect to the design variables can 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., ), and ANSYS  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., ) 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 of the eth element can be computed from Eqs. (21) and (22) through integration over the element domain as
where Ve is the volume (or area in a 2-D case) of the eth element.
When the unit cell is discretized into NE finite elements, the effective (homogenized) stiffness and thermal stress tensors given in Eqs. (10) and (11) can be written as
is the element stiffness matrix, is the strain-displacement matrix, denotes the element displacement under the applied strain , represents the global displacement vector associated with the element e under the mechanical loading obtainable from solving Eq. (14), and stands for the element displacement under and 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 function H() and the Dirac delta function involved in Eqs. (34)–(45) are approximated by the following piece-wise smooth functions :
where is a real positive constant that is usually chosen as 24 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 be and , in which c denotes 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 modulus KH is enforced to ensure the structural load-carrying capability. KH is 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., )
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 be E1 = 1 GPa, , 1/K for solid phase 1 (in blue) and E2 = 1 GPa, , 1/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 direction , with 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 fractions Vf1 and Vf2, respectively. The main results are summarized in Tab. 1, which contains the total number of iterations, minimum , effective 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 the x-direction 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., ). 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 . RBFs have been utilized in the studies on the topology optimization of 2-D mechanical metamaterials based on a parametric level set method  and on homogenization analyses of 3D printable interpenetrating phase composites by a meshfress method .
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.
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.
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., ).
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).
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.
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.
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.
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.
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 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 be E1 = 1 GPa, , 1/K for material 1 (in blue), and E2 = 1 GPa, , 1/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.
6.2.1 Minimum Anisotropic CTE with a Poisson’s Ratio Constraint
The first case is to obtain 3-D microstructures with the minimum under 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 unit cells are shown in Figs. 11b, 12b and 13b. In addition, Figs. 11c–11e, 12c–12e and 13c–13e display microstructural details from the x − y, x − z and y − z cross section views, respectively.
As shown in Figs. 11f–11h, the convergence is achieved after 116 iterations, giving a minimum of −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 minimum of −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 minimum of −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 under the constraint of Poisson’s ratio = 0, which is smaller than the magnitude of 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 .
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  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 mesh, which is not fine enough to produce smooth boundaries but remains a limitation of the current study for the 3-D three-phase 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., ).
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 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.
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., ). 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., ).
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.
|This work is licensed under a Creative Commons Attribution 4.0 International License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.|