Multi-Physics Inverse Homogenization for the Design of Innovative Cellular Materials: Application to Thermo-Elastic Problems

Abstract

We present a new algorithm to design lightweight cellular materials with required properties in a multi-physics context. In particular, we focus on a thermo-elastic setting by promoting the design of unit cells characterized both by an isotropic and an anisotropic behavior with respect to mechanical and thermal requirements. The proposed procedure generalizes the microSIMPATY algorithm to a thermo-elastic framework by preserving all the good properties of the reference design methodology. The resulting layouts exhibit non-standard topologies and are characterized by very sharp contours, thus limiting the post-processing before manufacturing. The new cellular materials are compared with the state-of-art in engineering practice in terms of thermo-elastic properties, thus highlighting the good performance of the new layouts which, in some cases, outperform the consolidated choices.

1. Introduction

Cellular materials represent an effective solution for structural applications where conventional monolithic materials fail to satisfy the design constraints [1]. The fast advancements in additive manufacturing technologies experienced in the last few years have further amplified the interest toward metamaterials. In addition, the possibility to employ a large variety of bulk materials in manufacturing processes (e.g., metals, polymers, ceramics [2,3,4]) has enabled the design of new metamaterials, featuring innovative combinations of physical effective properties. The possibility to blend different materials in order to reach diverse objectives proved to have a great impact in all the contexts where multi-functionality is required. For example, in [5,6,7,8], biocompatible 3D-printed metal bone implants promoting bone ingrowth are proposed by properly tailoring the material microstructure in order to reproduce the elastic modulus and the permeability of the human bone. Other applications range from thermal-cloaking systems fitly combining microstructure geometry and orientation [9,10] to lattice-based heat exchangers, where good thermal conductivity and convection properties are exploited to enhance the devices’ performance [11,12].

From a modeling viewpoint, the proposal of innovative multifunctional cellular materials can benefit from the most recent advancements in topology optimization [13], properly combined with direct and inverse homogenization processes [14,15,16]. Several optimization approaches can be exploited in the context of metamaterial design. The layout of the employed microstructures can be selected a priori, starting from consolidated dictionaries of unit cells [16,17,18,19,20,21], or designed from scratch to match the expected effective properties [22,23,24,25,26,27,28]. In this context, a single- or a multi-objective topology optimization at the microscale can drive the design of new unit cells matching target properties at the macroscale, potentially in a multi-physics framework. For instance, the optimization of homogenized elastic properties is tackled in [29,30,31] with the aim of maximizing the bulk (or shear) modulus. To this aim, the authors control specific components of the homogenized elastic tensor or resort to the minimization of the compliance of a given structural part. Other works focus on a multi-physics optimization (for instance, by considering elastic, thermal, and electrical properties) by providing microstructures optimized with respect to diverse objectives and physics [32,33,34,35].

Nevertheless, it is well-known that standard topology optimization techniques suffer from typical issues that may compromise the effective performance and manufacturability of the new layouts. Among the most recurrent, we mention the possible presence of intermediate densities, the non-smooth contours of the final design and the generation of unit cells which turn out to be unprintable since presenting too thin struts. All these drawbacks are strictly related to the selected computational grid: a coarse mesh promotes jagged boundaries and a diffused void/material interface; vice versa, an extremely fine mesh leads to a non-affordable computational effort and fosters the generation of too complex structures. Filtering offers a possible remedy to address all these concerns by alternating smoothing with sharpening phases to be properly tuned. Such a tuning is not a trivial task and may often lead to non-optimal design solutions [13,33,36].

The selection of a computational mesh customized to the design problem has been proved to be instrumental in order to limit the main issues of topology optimization. For instance, in [37], the combination of a standard density-based method for topology optimization with an anisotropic mesh adaptation procedure has been used to get rid of intermediate densities, irregular boundaries, and thin struts in the design of structures at the macroscale. The proposed algorithm, named SIMPATY (SIMP with mesh AdaptiviTY), is based on a robust mathematical tool, namely an a posteriori estimator for the discretization error, and it leads to final designs characterized by reliable mechanical properties as well as by free-form features. The same procedure has been successfully exploited at the microscale, with the proposal of the microSIMPATY algorithm [26]. So far, this procedure has been used for the design of unit cells with optimized mechanical properties in a linear elasticity setting [27,38].

In this work, we propose a new pipeline for the design of new cellular materials by extending the microSIMPATY algorithm to a multi-physics context. The objective is to obtain lightweight metamaterials with prescribed requirements on the elastic and thermal conductivity properties, which are characterized by a ready-to-print topology. The design strategy here developed is confined to a 2D setting and has to be meant as a proof-of-concept, preliminary to a 3D implementation. However, to corroborate the effectiveness of the proposed methodology, we perform a cross-comparison between the new cells and the standard ones in thermo-elastic applications.

The paper is organized as follows. Section 2 represents the core of the paper. It provides the physical problem constraining the optimization process, outlines the main theoretical tools to perform the optimization, and formalizes the thermo-elastic design procedure in the MultiP-microSIMPATY algorithm. Three design cases are considered in Section 3 to apply the MultiP-microSIMPATY algorithm to diverse scenarios. Section 4 further analyzes the results in the previous section by comparing the new designs with the state-of-the-art. Finally, Section 5 outlines the most remarkable contributions of the work together with some future perspectives.

2. Methods

In this paper, we refer to a multi-physics framework, in order to provide new layouts for the design of cellular materials. A standard issue consists in optimizing the microscale to ensure desired properties at the macroscale. To deal with this two-scale setting, it is crucial to properly transfer the physical characterization of the micro- to the macroscale. Direct and inverse homogenization represent widespread solutions in such a direction [14,39,40,41]. In particular, the direct approach incorporates the microscopic effects into a homogenized macroscopic model. As a consequence, the microscopic behavior is known, whereas we have to identify the (homogenized) macroscopic characterization. Vice versa, inverse homogenization starts from desired macroscopic physical properties and designs the microscale in order to match such features, thus swapping the role played by known and unknown scales with respect to direct homogenization.

In this paper, we focus on a two-dimensional setting and on linear thermo-elastic properties, so that, at the macroscale, the reference models are the linear elasticity equation [42] and the linear thermal conduction problem, as identified by the standard stress–strain ($\sigma$-$\epsilon$) and heat flux-temperature ($\mathbf{q}$-$\theta$) relations, given by

$$\sigma \left(\mathbf{u}\right)=\left[\begin{array}{c}{\sigma }_{11}\left(\mathbf{u}\right)\\ {\sigma }_{22}\left(\mathbf{u}\right)\\ {\sigma }_{12}\left(\mathbf{u}\right)\end{array}\right]=\left[\begin{array}{c}{E}_{1111}{E}_{1122}{E}_{1112}\\ E_{2211}{E}_{2222}{E}_{2212}\\ E_{1211}{E}_{1222}{E}_{1212}\end{array}\right]\left[\begin{array}{c}{\epsilon }_{11}\left(\mathbf{u}\right)\\ {\epsilon }_{22}\left(\mathbf{u}\right)\\ 2{\epsilon }_{12}\left(\mathbf{u}\right)\end{array}\right]=\mathbf{E}\epsilon \left(\mathbf{u}\right),$$

(1)

and

$$\begin{array}{c}\mathbf{q}(\theta )=\left[\begin{array}{c}{q}_1(\theta )\\ q_2(\theta )\end{array}\right]=\left[\begin{array}{c}{k}_{11}{k}_{12}\\ k_{21}{k}_{22}\end{array}\right]\left[\begin{array}{c}{\displaystyle \frac{\partial \theta }{\partial x_1}}\\ {\displaystyle \frac{\partial \theta }{\partial x_2}}\end{array}\right]=\mathbf{k}\nabla \theta ,\end{array}$$

(2)

respectively, where $\mathbf{E}$ and $\mathbf{k}$ are the stiffness and the conductivity tensors characterizing the considered solid material.

We observe that in view of the homogenization procedures, the constitutive laws (1) and (2) have to be properly modified to include the effects of the microscale, into

$${\sigma }^H\left(\mathbf{u}\right)={\mathbf{E}}^H\epsilon \left(\mathbf{u}\right),{\mathbf{q}}^H(\theta )={\mathbf{k}}^H\nabla \theta ,$$

where ${\mathbf{E}}^H$ and ${\mathbf{k}}^H$ denote the homogenized stiffness and thermal conductivity tensors, as detailed in the next section.

2.1. Inverse Homogenization

Inverse homogenization is the procedure that allows us to design microstructures with prescribed properties at the macroscale. The required features are mathematically commuted into a goal functional $\mathcal{J}$ and into suitable constraints driving a topology optimization process to be solved in the unit cell $Y\subset {\mathbb{R}}^2$ whose periodic repetition yields the cellular material  [26,27,29,38,43]. According to a density-based approach, a standard way to perform such an optimization leads us to define an auxiliary scalar field, $\rho$, that models the relative material density at the microscale. A priori, it is assumed that $\rho =1$ labels the material, while $\rho =0$ identifies the void. However, since density $\rho \in L^{\infty }\left(Y,[0,1]\right)$ can take all the values in $[0,1]$, it is standard to penalize the intermediate values (i.e., intermediate material densities) that are not physically consistent. To this aim, we resort to the SIMP method, which modifies the reference state equations by weighting the constitutive laws with a suitable power of the density [13].

In general, the optimization problem we are interested to solve is

$${\displaystyle \underset{\rho \in L^{\infty }\left(Y,[0,1]\right)}{min}\mathcal{J}(z(\rho ),\rho ):\left\{\begin{array}{c}{a}_{\rho }\left(z(\rho ),w\right)=F_{\rho }\left(w\right)\forall w\in W\hfill \\ {\mathbf{L}}_B\le \mathbf{C}(z(\rho ),\rho )\le {\mathbf{U}}_B,\hfill \end{array}\right.}$$

(3)

where $z=z(\rho )$ denotes the state variable depending on the density field, $a_{\rho }(·,·)$ together with $F_{\rho }(·)$ defines the state equation constraining the topology optimization, W is a suitable function space [44], and the box inequality includes specific design and physical requirements, with $\mathbf{C}(·,·)$ being the vector gathering the quantities to be controlled through the corresponding lower and upper bounds, ${\mathbf{L}}_B$ a ${\mathbf{U}}_B$.

In the analysis below, we pick the objective functional $\mathcal{J}$ as

$$\mathcal{M}(\rho )={\int }_Y\rho dY$$

(4)

since we are interested in minimizing the total mass, $\mathcal{M}$, of the cellular structure, i.e., to design lightweight materials.

According to a standard homogenization procedure, as the state equation, we select the linear elasticity model at the microscale weighed by the density function, which describes the Y-periodic displacement field fluctuations, ${\mathbf{u}}^{*,ij}$, induced by the test displacement fields ${\mathbf{u}}^{0,ij}$, with ${\mathbf{u}}^{0,11}={[x,0]}^{\mathrm{T}}$, ${\mathbf{u}}^{0,22}={[0,y]}^{\mathrm{T}}$ and ${\mathbf{u}}^{0,12}={[y,0]}^{\mathrm{T}}$. This leads us to identify the forms $a_{\rho }(·,·)$ and $F_{\rho }(·)$ in (3) with

$$\begin{array}{ccc}\hfill a_{\rho }^{E,ij}\left({\mathbf{u}}^{*,ij}(\rho ),\mathbf{v}\right)& =& {\displaystyle {\displaystyle \frac{1}{\left|Y\right|}}{\int }_Y{\rho }^p\sigma \left({\mathbf{u}}^{*,ij}\right):\epsilon \left(\mathbf{v}\right)dY,}\hfill \\ \hfill F_{\rho }^{E,ij}\left(\mathbf{v}\right)& =& {\displaystyle {\displaystyle \frac{1}{\left|Y\right|}}{\int }_Y{\rho }^p\sigma \left({\mathbf{u}}^{0,ij}\right):\epsilon \left(\mathbf{v}\right)dY,}\hfill \end{array}$$

(5)

respectively, with $p\in {\mathbb{R}}^{+}$, $ij\in I=\{11,22,12\}$, and where the superscript E refers to the elasticity setting. The state equation associated with (5) is completed with fully periodic conditions on the cell boundary $\partial Y$, according to the asymptotic homogenization theory. Thus, W in (3) coincides with the space ${\mathcal{U}}_{\#}^2={\left[H_{⥁}^1\left(Y\right)\right]}^2$ of the $H^1\left(Y\right)$-vector functions satisfying periodic boundary conditions along $\partial Y$.

To include also the thermal component in the topology optimization, we further constrain the process with the $\rho$-weighed thermal conductivity model at the microscale, which are characterized by the forms

$$\begin{array}{ccc}\hfill a_{\rho }^{k,m}({\theta }^{*,m}(\rho ),v)& =& {\displaystyle \frac{1}{\left|Y\right|}{\int }_Y{\rho }^s\mathbf{q}\left({\theta }^{*,m}\right):\nabla vdY,}\hfill \\ \hfill F_{\rho }^{k,m}\left(v\right)& =& {\displaystyle \frac{1}{\left|Y\right|}{\int }_Y{\rho }^s\mathbf{q}\left({\theta }^{0,m}\right):\nabla vdY,}\hfill \end{array}$$

(6)

with $s\in {\mathbb{R}}^{+}$, $m\in J=\{1,2\}$, and where the superscript k refers to the thermal framework. Analogously to (5), we complete the thermal state equation identified by (6) with periodic boundary conditions along $\partial Y$, so that ${\theta }^{*,m}$ and $v\in {\mathcal{U}}_{\#}^1=H_{⥁}^1\left(Y\right)$, where ${\theta }^{*,m}$ denotes the temperature fluctuations associated with the test temperature fields ${\theta }^{0,m}$ (namely, ${\theta }^{0,1}=x$ and ${\theta }^{0,2}=y$).

The two problems at the microscale, (5) and (6), are instrumental to define the homogenized elastic tensor, ${\mathbf{E}}^H$, and the homogenized thermal conductivity tensor, ${\mathbf{k}}^H$, which will be involved in the box constraints in (3). The component-wise definition of ${\mathbf{E}}^H$ and ${\mathbf{k}}^H$ is

$${\displaystyle E_{ijkl}^H=\frac{1}{\left|Y\right|}{\int }_Y{\rho }^p\left[\sigma \left({\mathbf{u}}^{0,ij}\right)-\sigma \left({\mathbf{u}}^{*,ij}(\rho )\right)\right]:\left[\epsilon \left({\mathbf{u}}^{0,kl}\right)-\epsilon \left({\mathbf{u}}^{*,kl}(\rho )\right)\right]dY,}$$

(7)

$${\displaystyle k_{mn}^H=\frac{1}{\left|Y\right|}{\int }_Y{\rho }^s\left[\mathbf{q}\left({\theta }^{0,m}\right)-\mathbf{q}\left({\theta }^{*,m}(\rho )\right)\right]:\left[\nabla {\theta }^{0,n}-\nabla {\theta }^{*,n}(\rho )\right]dY,}$$

(8)

respectively, with $ij,kl\in I$ and $m,n\in J$. In particular, the two-sided inequality in (3) will be exploited to promote diverse mechanical and thermal behaviors along the different spatial directions. To this aim, we constrain the two ratios $E_{2222}^H/E_{1111}^H$ and $k_{22}^H/k_{11}^H$ so that they vary in suitable ranges. This choice allows us to penalize the mechanical and the thermal contributions in a different way along the two directions, as shown in the numerical assessment. An additional two-sided control is enforced on the first and the last diagonal terms, $E_{1111}^H$ and $E_{1212}^H$, of the homogenized elastic tensor, as well as on the first diagonal term, $k_{11}^H$, of the homogenized thermal conductivity tensor.

To sum up, the optimization setting we are led to deal with coincides with the following problem:

$${\displaystyle \underset{\rho \in L^{\infty }\left(Y,[0,1]\right)}{min}\mathcal{M}(\rho ):\left\{\begin{array}{c}{a}_{\rho }^{E,ij}\left({\mathbf{u}}^{*,ij}(\rho ),\mathbf{v}\right)=F_{\rho }^{E,ij}\left(\mathbf{v}\right)\forall \mathbf{v}\in {\mathcal{U}}_{\#}^2,ij\in I\hfill \\ a_{\rho }^{k,m}\left({\theta }^{*,m}(\rho ),v\right)=F_{\rho }^{k,m}\left(v\right)\forall v\in {\mathcal{U}}_{\#}^1,m\in J\hfill \\ E_{1111}^{\mathrm{low}}\le E_{1111}^H\le E_{1111}^{\mathrm{up}}\hfill \\ E_{1212}^{\mathrm{low}}\le E_{1212}^H\le E_{1212}^{\mathrm{up}}\hfill \\ {\left({\displaystyle \frac{E_{2222}}{E_{1111}}}\right)}^{\mathrm{low}}\le {\displaystyle \frac{E_{2222}^H}{E_{1111}^H}}\le {\left({\displaystyle \frac{E_{2222}}{E_{1111}}}\right)}^{\mathrm{up}}\hfill \\ k_{11}^{\mathrm{low}}\le k_{11}^H\le k_{11}^{\mathrm{up}}\hfill \\ {\left({\displaystyle \frac{k_{22}}{k_{11}}}\right)}^{\mathrm{low}}\le {\displaystyle \frac{k_{22}^H}{k_{11}^H}}\le {\left({\displaystyle \frac{k_{22}}{k_{11}}}\right)}^{\mathrm{up}}\hfill \\ {\rho }_{\min }\le \rho \le 1\hfill \end{array}\right.}$$

(9)

where all the bound values, ${(·)}^{\mathrm{low}}$ and ${(·)}^{\mathrm{up}}$, will be set according to the application at hand and in order to avoid an unfeasible solution (inappropriate constraints might lead to an empty solution space). Finally, the last inequality in (9) is meant to ensure the well-posedness of both the elasticity and the thermal problems (5) and (6), ${\rho }_{min}$ being a suitable value in $(0,1)$ (see Section 3 for more details).

2.2. Discretization on Anisotropic Adapted Meshes

With a view to the solution of problem (9), all the quantities involved in the state equations, as well as in the constraints, have to be discretized on a suitable tessellation of the unit cell Y. For this purpose, we resort to a computational mesh ${\mathcal{T}}_h=\left\{K\right\}$ customized to the problem at hand and characterized by stretched elements (i.e., a so-called anisotropic adapted mesh). Mesh ${\mathcal{T}}_h$ is employed to discretize both the test and the trial functions in the state equations, as well as the density function $\rho$, by means of a finite element scheme [44].

The anisotropic reference setting is the one proposed in [45]. In particular, the anisotropic features of each element K coincide with the lengths, ${\lambda }_{1,K}$, ${\lambda }_{2,K}$, and the directions, ${\mathbf{r}}_{1,K}$, ${\mathbf{r}}_{2,K}$, of the semi-axes of the ellipse circumscribed to K through the standard affine map, $T_K:\widehat{K}\to K$, between the reference element $\widehat{K}$ and the triangle K.

Concerning the adaptation procedure, we resort to a metric-based approach driven by an a posteriori estimator for the discretization error associated with the density function $\rho$. Among the error estimators available in the literature [46,47], we refer to an a posteriori recovery-based error analysis. Following the seminal work by O.C. Zienkewicz and J.Z. Zhu [48], we control the $H^1$-seminorm of the discretization error on the density, $e_{\rho }=\rho -{\rho }_h$. The selection of such an estimator is motivated by the fact that the density $\rho$ exhibits strong gradients (i.e., large values for the $H^1$-seminorm) across the material–void interface. This feature will yield meshes whose elements are crowded along the boundaries of the structure, thus promoting the design of very smooth layouts. To this aim, we exactly integrate the so-called recovered error, ${\mathbf{E}}_{\nabla }=P\left(\nabla {\rho }_h\right)-\nabla {\rho }_h$, namely,

$$\begin{array}{cc}\hfill |e_{\rho }{|}_{H^1\left(Y\right)}^2& {\displaystyle =\parallel \nabla e_{\rho }{\parallel }_{L^2\left(Y\right)}^2={\int }_Y{\left|\nabla \rho -\nabla {\rho }_h\right|}^{2}dY}\hfill \\ \hfill & {\displaystyle \simeq \parallel {\mathbf{E}}_{\nabla }{\parallel }_{L^2\left(Y\right)}^2={\int }_Y{|P\left(\nabla {\rho }_h\right)-\nabla {\rho }_h|}^{2}dY,}\hfill \end{array}$$

(10)

where ${\rho }_h$ denotes the finite element discretization of $\rho$ in the space $V_h^r$ of the piecewise polynomials of degree $r\in \mathbb{N}$ associated with ${\mathcal{T}}_h$. The operator $P:{\left[V_h^{r-1}\right]}^2\to {\left[V_h^s\right]}^2$ in (10), with $s\in \mathbb{N}$, denotes the recovered gradient, which, in general, provides a more accurate estimate of the exact gradient $\nabla \rho$ with respect to the discrete gradient $\nabla {\rho }_h$. Several recipes are available in the literature to define P [49,50,51,52]. In particular, we select operator $P:{\left[V_h^0\right]}^2\to {\left[V_h^0\right]}^2$ as the area-weighted average of $\nabla {\rho }_h$ over the patch of the elements, ${\Delta }_K=\{T\in {\mathcal{T}}_h:T\cap K\ne \varnothing \}$, associated with K; i.e., we opt for

$${\displaystyle P\left(\nabla {\rho }_h\right)\left(\mathbf{x}\right)={\displaystyle \frac{1}{|{\Delta }_K|}}\sum _{T\in {\Delta }_K}\left|T\right|\nabla {\rho }_h{|}_T\forall \mathbf{x}\in K,}$$

(11)

with $|\omega |$ the area of the generic domain $\omega \subset {\mathbb{R}}^2$, where we have set the degree of the finite element space for ${\rho }_h$ to $r=1$. Space $V_h^1$ is also adopted to discretize the components of the displacement vectors ${\mathbf{u}}^{*,ij}$ as well as the temperature fields ${\theta }^{*,m}$ in (9), with $ij\in I$ and $m\in J$.

According to [53,54], we here adopt the anisotropic generalization of (10). This estimator essentially exploits the anisotropic counterpart of the definition of the $H^1$-seminorm [45], based on the symmetric semidefinite positive matrix $G_{{\Delta }_K}$, with entries

$${\displaystyle {\left[G_{{\Delta }_K}(\nabla g)\right]}_{i,j}=\sum _{T\in {\Delta }_K}{\int }_T{\displaystyle \frac{\partial g}{\partial x_i}}{\displaystyle \frac{\partial g}{\partial x_j}}dTi,j=1,2,}$$

(12)

with $g\in H^1\left(Y\right)$, and where it is understood $x_1=x$ and $x_2=y$. Thus, the squared $H^1$-seminorm $|e_{\rho }{|}_{H^1\left(Y\right)}^2$ is evaluated by the (global) error estimator ${\eta }^2={\sum }_{K\in {\mathcal{T}}_h}{\eta }_K^2$, where

$${\displaystyle {\eta }_K^2=\frac{1}{{\lambda }_{1,K}{\lambda }_{2,K}}\sum _{i=1}^2{\lambda }_{i,K}^2\left({\mathbf{r}}_{i,K}^{\mathrm{T}}{G}_{{\Delta }_K}\left({\mathbf{E}}_{\nabla }\right){\mathbf{r}}_{i,K}\right),}$$

(13)

defines the local error estimator. The contribution between brackets coincides with the projection of the squared $L^2$-norm of the recovered error along the anisotropic directions, while the scaling factor ${\left({\lambda }_{1,K}{\lambda }_{2,K}\right)}^{-1}$ guarantees the consistency with the isotropic case (for more details, see [53]).

The new adapted mesh is generated after commuting the error estimator ${\eta }_K$ into a new mesh spacing (the metric), $\mathcal{M}$, consisting of the triplet $\{{\lambda }_{1,K}^{adapt},{\lambda }_{2,K}^{adapt},{\mathbf{r}}_{1,K}^{adapt}\}$, where the direction ${\mathbf{r}}_{2,K}^{adapt}$ is automatically defined being ${\mathbf{r}}_{1,K}^{adapt}·{\mathbf{r}}_{2,K}^{adapt}=0$, for each element $K\in {\mathcal{T}}_h$. This operation is performed by taking into account three different criteria, namely, (i) the minimization of the mesh cardinality $\#{\mathcal{T}}_h$; (ii) an accuracy requirement on the discretization error $|e_{\rho }{|}_{H^1\left(Y\right)}$ (i.e, on the error estimator $\eta$), controlled up to a user-defined tolerance TOL; (iii) the equidistribution of the error throughout the mesh elements (i.e., ${\eta }_K^2={\mathtt{TOL}}^2/\#{\mathcal{T}}_h$). These three criteria lead us to solve a constrained minimization problem on each triangle $K\in {\mathcal{T}}_h$. The solution to this local optimization problem can be analytically derived, as proved in [55], being

$$\begin{array}{c}{\lambda }_{1,K}^{adapt}=g_2^{-1/2}{\left(\frac{{\mathtt{TOL}}^2}{2\#{\mathcal{T}}_h\left|{\widehat{\Delta }}_K\right|}\right)}^{1/2},{\mathbf{r}}_{1,K}^{adapt}={\mathbf{g}}_2,\\ {\lambda }_{2,K}^{adapt}=g_1^{-1/2}{\left(\frac{{\mathtt{TOL}}^2}{2\#{\mathcal{T}}_h\left|{\widehat{\Delta }}_K\right|}\right)}^{1/2},{\mathbf{r}}_{2,K}^{adapt}={\mathbf{g}}_1\end{array}$$

(14)

where $g_1$, $g_2$ and ${\mathbf{g}}_1$, ${\mathbf{g}}_2$ are the eigenvalues and the eigenvectors of the scaled matrix ${\widehat{G}}_{{\Delta }_K}\left({\mathbf{E}}_{\nabla }\right)=G_{{\Delta }_K}\left({\mathbf{E}}_{\nabla }\right)/\left|{\Delta }_K\right|$, with $g_1\ge g_2>0$.

Finally, the metric $\mathcal{M}={\{{\lambda }_{1,K}^{adapt},{\lambda }_{2,K}^{adapt},{\mathbf{r}}_{1,K}^{adapt}\}}_{K\in {\mathcal{T}}_h}$ has to be changed into a quantity associated with the vertices of ${\mathcal{T}}_h$ received as an input by the selected mesh generator. A standard choice consists in an arithmetic mean formula applied to the patch of elements associated with each vertex in ${\mathcal{T}}_h$ [54].

The anisotropic mesh adaptation based on the metric (14) is customized to a topology optimization problem in the algorithm SIMPATY, proposed in [37]. This procedure has been successfully employed for the design of structures at the macroscale [37,56,57] as well as for the design of new metamaterials with the proposal of the algorithm microSIMPATY [26,27]. Moreover, a combination of topology optimization at the macroscale and at the microscale is carried out in [38]. In particular, a multiscale topology optimization process is used for the design of orthotic devices for 3D printing manufacturing with the proposal of patient-specific innovative solutions.

It has been verified that the adoption of an adapted anisotropic mesh leads to free-form layouts characterized by very smooth boundaries both at the macroscale and at the microscale, mitigating some of the well-known drawbacks of standard topology optimization, such as the massive employment of filtering, the staircase effect, and the generation of too complex structures [13,33,36]. However, in [57], it has been observed that the presence of deformed elements inside the structures makes the finite element analysis less reliable. To overcome this issue, the authors suggest a hybrid approach. Thus, the mesh is kept isotropic, with a uniform diameter $h^{\mathrm{iso}}$ in the full-material regions, $\{\mathbf{x}\in Y:{\rho }_h\left(\mathbf{x}\right)>{\rho }^{\mathrm{th}}\}$ with ${\rho }^{\mathrm{th}}$ as a user-defined threshold, whereas the stretched triangles are preserved along the material–void interface. Actually, these hybrid meshes ensure an effective balance between the smoothness of the structure and robust engineering performances (the interested reader can find a quantitative investigation of the benefits of the hybrid approach in terms of accuracy in (Section 5, [57])). For this reason, we resort to hybrid meshes in the sequel.

2.3. Multi-Physics Optimization Algorithm

In this section we propose the multi-physics adaptive inverse homogenization procedure, which generalizes the algorithm proposed in [26]. The discretization of the state equations associated with (5) and (6) is performed with the open-source finite element solver FreeFEM [58], which provides the ideal environment to implement an anisotropic mesh adaptation procedure in Section 2.2 through the built-in mesh generator BAMG (Bidimensional Anisotropic Mesh Generator).

The developed multi-physics optimization procedure is listed in the pseudocode in Algorithm 1. The main loop (lines 3–12) includes an optimization step, a filtering phase, and the mesh adaptation. At each global iteration k, the optimization problem is solved (line 4, function $\mathtt{optimize}$) by taking into account all the constraints on the components of the elastic and of the thermal conductivity tensors in (9). To this aim, we use the interior point algorithm IPOPT [59], although any other optimization tool can be selected [60]. IPOPT requires as input the functional $\mathcal{J}$ to be minimized; the vector $\mathcal{C}$ gathering the constrained quantities in the optimization procedure; the two vectors ${\mathbf{c}}^l$ and ${\mathbf{c}}^u$ of the lower and upper bounds for the components in $\mathcal{C}$; the array $\mathcal{G}$ collecting the derivative of the functional $\mathcal{J}$ and of the constraints $\mathcal{C}$ with respect to $\rho$, computed by the adjoint Lagrangian approach (for more details, we refer to [27]); the initial guess ${\rho }_h^{\mathtt{k}}$ to start the optimization process; the accuracy $\mathtt{TOPT}$ for the minimization problem; the maximum number of iterations $\mathtt{IT}$ to stop the optimization. In particular, in the numerical assessment of Section 3, we set $\mathtt{TOPT}={10}^{-5}$, and $\mathtt{IT}=100$ for $\mathtt{k}=0$ and $\mathtt{IT}=10$ for all the successive iterations. The higher value for $\mathtt{IT}$ for $\mathtt{k}=0$ takes into account that the initial guess ${\rho }_h^{\mathtt{0}}$ can be completely arbitrary with respect to the minimum to be reached. On the contrary, a smaller value for $\mathtt{IT}$ is sufficient for $\mathtt{k}>0$, since the initial guess, ${\rho }_h^{\mathtt{k}}$, coincides with the output of a previous optimization step.

Algorithm 1 MultiP-microSIMPATY

1: Input: CTOL, kmax, ${\mathbf{c}}^l$, ${\mathbf{c}}^u$, ${\rho }_h^0$, TOPT, IT, kfmax, $\tau$, $\beta$, ${\mathcal{T}}_h^0$, TOL, $\mathtt{HYB}$ 2: Set: k = 0, errC = 1+CTOL; 3: while errC > CTOL & k< kmax do 4:     ${\rho }_h^{\mathtt{k}}+\mathtt{1}$ = optimize($\mathcal{J}$, $\mathcal{C}$, ${\mathbf{c}}^l$, ${\mathbf{c}}^u$, $\mathcal{G}$, ${\rho }_h^{\mathtt{k}}$, TOPT, IT); 5:     if $\mathtt{k}<\mathtt{kfmax}$ then 6:         ${\rho }_h^{\mathtt{k}}+\mathtt{1}$ = helmholtz(${\rho }_h^{\mathtt{k}}+\mathtt{1}$, $\tau$); 7:         ${\rho }_h^{\mathtt{k}}+\mathtt{1}$ = heaviside(${\rho }_h^{\mathtt{k}}+\mathtt{1}$, $\beta$); 8:     end if 9:     ${\mathcal{T}}_h^{\mathtt{k}+\mathtt{1}}$ = adapt(${\mathcal{T}}_h^{\mathtt{k}}$, ${\rho }_h^{\mathtt{k}+\mathtt{1}}$,TOL, $\mathtt{HYB}$); 10:     errC = $|\#{\mathcal{T}}_h^{\mathtt{k}+\mathtt{1}}-\#{\mathcal{T}}_h^{\mathtt{k}}|/\#{\mathcal{T}}_h^{\mathtt{k}}$; 11:     k = k+1; 12: end while 13: ${\mathcal{T}}_h$ = ${\mathcal{T}}_h^{\mathtt{k}}$; 14: ${\rho }_h={\rho }_h^{\mathtt{k}}$; 15: $\left[{\mathbf{E}}^H,{\mathbf{k}}^H\right]$ = homogenize(${\rho }_h$); 16: return ${\mathcal{T}}_h$, ${\rho }_h$, ${\mathbf{E}}^H$, ${\mathbf{k}}^H$

Function $\mathtt{optimize}$ returns the density ${\rho }_h^{\mathtt{k}+\mathtt{1}}$, which is successively processed by means of Helmholtz and Heaviside filters (lines 6–7, functions $\mathtt{helmholtz}$ and $\mathtt{heaviside}$) [61,62]. The two filtering operations work in a complementary way. The Helmholtz partial differential equation is instrumental to remove too thin features, although promoting intermediate densities along the layout contour. In more detail, it consists of a low-pass filter based on a diffusion kernel with radius $\tau \in {\mathbb{R}}^{+}$. On the contrary, the Heaviside filter, coinciding with a $\beta$-dependent regularization of the Heaviside function with $\beta \in {\mathbb{R}}^{+}$, penalizes the intermediate material densities, also due to the Helmholtz filter, thus increasing the sharpness of the material/void interface. The combined filtering takes place for the first $\mathtt{kfmax}$ global iterations only. This choice leads to start the mesh adaptation procedure with a density field, which is free from too complex features, while exhibiting a clear alternation between void and material. The filtering phase becomes redundant when the optimization loop approaches the minimum, so that mesh adaptation alone suffices to ensure well-defined structures. In the next section, filtering parameters $\tau$ and $\beta$ are set equal to $0.02$ and 5 respectively, while $\mathtt{kfmax}=25$.

The next step coincides with the mesh adaptation procedure detailed in Section 2.2 and here represented by function $\mathtt{adapt}$ (line 9). The input parameter $\mathtt{TOL}$ establishes the accuracy of the error estimator $\eta$ through the predicted metric in (14). Parameter $\mathtt{HYB}$ is a boolean flag that, in correspondence with the full material, switches the employment of an isotropic mesh on or off.

The main loop is controlled by a check on the stagnation of the relative difference between the cardinality of two consecutive meshes (line 10), up to a maximum number of global iterations $\mathtt{kmax}$ (line 3). The choices TOL $={10}^{-5}$ and kmax$=100$ are preserved throughout all the numerical assessment below.

Algorithm MultiP-microSIMPATY returns the final adapted mesh ${\mathcal{T}}_h$, the optimized density ${\rho }_h$, and the homogenized elastic and conductivity tensors, ${\mathbf{E}}^H$ and ${\mathbf{k}}^H$, which are computed by function $\mathtt{homogenize}$ (line 15), based on (7) and (8).

We remark that the procedure itemized in Algorithm 1 is fully general and it can be applied in a straightforward way to different multi-physics contexts after properly modifying the formulation in (9).

3. Results

We analyze three different cases of microstructure design according to (9). In order to highlight the interplay between the different (thermal and mechanical) physics involved, we consider configurations where the thermal conductivity and the elastic stiffness requirements act along different directions. For instance, a high shear stiffness combined with a high thermal conductivity along the x-direction orient the material along two opposite directions, with the prescription of a conflict configuration.

The whole verification below shares common choices for some physical quantities and discretization parameters. In particular, the unit cell $Y\subset {\mathbb{R}}^2$ is identified with the unitary square, $Y={(0,1)}^2$. Moreover, we set the Young modulus, $E_Y$, and the Poisson ratio, $\nu$, to 1 and $0.3$, respectively, and we consider an isotropic solid material with unitary thermal conductivity by setting $k_{11}=k_{22}=1$ and $k_{12}=k_{21}=0$. These choices allow us to obtain normalized homogenized mechanical and thermal properties for the cellular structures. Following [26], both the SIMP-powers, p and s, in (7) and (8) are chosen equal to 4 to penalize intermediate densities.

Concerning the discretization frame, we choose a random density field, ${\rho }_h^0$, as the initial guess for the optimization process, in order to avoid any bias. In particular, ${\rho }_h^0$ is defined on an initial structured mesh characterized by 30 subdivisions per side, and with values ranging from ${\rho }_{min}={10}^{-4}$ to 1 (see Figure 1 for an example).

Finally, to ensure a reliable finite element analysis, we resort to the hybrid mesh adaptation procedure ($\mathtt{HYB}$ = 1 in function $\mathtt{adapt}$). In particular, we choose the threshold value ${\rho }^{\mathrm{th}}=0.9$ to manage the alternation between isotropic and anisotropic elements, and the isotropic tessellation is characterized by the uniform diameter $h^{\mathrm{iso}}=0.03$ (approximately $1/30$ of the design domain dimension).

After the optimization, we perform a verification step to check the actual mechanical and thermal properties of the material yielded by a periodic repetition of the optimized unit cell. To this aim, we use the Abaqus software (Abaqus, Dassault Systèmes Simulia Corp, Johnston, RI, USA). The layouts provided by Algorithm 1 are imported in Abaqus after a thresholding, which neglects the density smaller than 0.75. The obtained geometry is remeshed on a uniform isotropic triangular mesh with an average size equal to $0.01$, while the displacement and temperature fields are discretized with quadratic finite elements, completed with periodic boundary conditions. The verification here performed can be considered as a preliminary step toward the integration of the MultiP-microSIMPATY algorithm into a common workflow for structural analysis.

3.1. Design Case 1

The main goal of this first optimization process is to design a lightweight unit cell characterized by isotropic mechanical homogenized properties and, vice versa, anisotropic thermal homogenized features. This problem can be cast in setting (9), after making the following choices for the constraints:

$$\left\{\begin{array}{c}0.05\le E_{1111}^H\le 0.08\hfill \\ 0.055\le E_{1212}^H\le 0.080\hfill \\ 1\le {\displaystyle \frac{E_{2222}^H}{E_{1111}^H}}\le 2\hfill \\ 0.01\le k_{11}^H\le 1.00\hfill \\ 0.00\le {\displaystyle \frac{k_{22}^H}{k_{11}^H}}\le 0.58.\hfill \end{array}\right.$$

(15)

The isotropic mechanical behavior and the anisotropic thermal properties are enforced by the constraints in (15)${}_3$ and (15)${}_5$. In particular, we expect ratios $E_{2222}^H/E_{1111}^H$ and $k_{22}^H/k_{11}^H$ to coincide with the corresponding lower and upper bounds, respectively. Moreover, since a control on the ratios does not ensure $E_{1111}^H$, $E_{2222}^H$, $k_{11}^H$, and $k_{22}^H$ to be in a physically admissible range of values, we further constrain the optimization through the box inequalities (15)${}_1$ and (15)${}_4$. Finally, a control on the component $E_{1212}^H$ of the homogenized stiffness tensor closes the minimization problem, thus further restricting the solution space.

For the values set for the input parameters, the MultiP-microSIMPATY algorithm converges in 51 global iterations. Figure 2 shows the layout and the associated anisotropic adapted mesh at three different iterations.

We remark that the final topology of the layout is already detected at the first iteration, although the quality of the solution is improved throughout the optimization process. In particular, at the first iteration ($\mathtt{k}=1$), we observe a significant staircase effect together with the presence of intermediate densities along the microstructure interface. At the end of the filtering phase ($\mathtt{k}=24$), the jagged boundaries are fully smoothed, despite the intermediate densities still blurring the design. The spreading effect along the material/void interface is gradually reduced when switching off the filtering, i.e., for $\mathtt{k}>24$, as shown by the last column in Figure 2. Thus, the final optimized solution ($\mathtt{k}=51$) shows an extremely sharp transition from material to void and smooth boundaries, which make the structure ready for printing or manufacturing, with a limited need for post-processing.

Concerning the adapted mesh, we recognize the effect of the hybrid approach, which combines stretched elements to discretize the strong gradients of the density field, coarse anisotropic triangles outside the structure, and isotropic elements in correspondence with the material.

Additional quantitative information on the MultiP-microSIMPATY algorithm is provided by

Table 1. Design cases 1, 2, and 3: values of the constraints and of the objective functional computed with Abaqus software, together with the lower and the upper bounds, $c^l$ and $c^u$, involved in the optimization.

	${\mathit{E}}_{1111}^{\mathit{H}}$	${\mathit{E}}_{1212}^{\mathit{H}}$	$\frac{{\mathit{E}}_{2222}^{\mathit{H}}}{{\mathit{E}}_{1111}^{\mathit{H}}}$	${\mathit{k}}_{11}^{\mathit{H}}$	$\frac{{\mathit{k}}_{22}^{\mathit{H}}}{{\mathit{k}}_{11}^{\mathit{H}}}$	$\mathcal{M}$
			Design Case 1
$c^u$	0.080	0.080	2.000	1.000	0.580
c	0.038	0.056	1.299	0.199	0.566
$c^l$	0.050	0.055	1.000	0.010	0.000	0.292
			Design Case 2
$c^u$	0.350	0.150	2.000	1.000	2.000
c	0.250	0.086	0.299	0.317	0.597
$c^l$	0.230	0.080	0.300	0.300	0.000	0.412
			Design Case 3
$c^u$	0.150	0.100	1.100	0.400	1.100
c	0.151	0.083	1.074	0.260	1.002
$c^l$	0.100	0.080	1.000	0.250	1.000	0.415

and by the diagrams in Figure 3, which show the evolution of the objective function and of the constrained quantities (top), together with the trend of the mesh cardinality (bottom), over the global iterations. Notice that the values of the constraints have been normalized between 0 and 1 (see the highlighted area in the top panel of Figure 3). It is evident that the mass exhibits a completely different trend when compared with the constrained quantities. The value of the objective function oscillates with values between 0.325 and 0.475 over the first 35 iterations, and it eventually converges toward a stable phase. On the contrary, all the constrained quantities are characterized by mild oscillations. In particular, $k_{11}^H$ remains essentially constant over the whole optimization process. The plot of the ratios $E_{2222}^H/E_{1111}^H$ and $k_{22}^H/k_{11}^H$ confirms that the two inequalities are in conflict so that the active constraints are the lower and upper bound, respectively. Moreover, from the values in Table 1, it can be observed that the stiffness component along the x-direction, $E_{1111}^H$, reaches a value that is about 25% lower than the corresponding $c^l$. This can be ascribed to the presence of very thin struts generated by the severe thresholding (${\rho }_h<0.75$) applied before performing the analyses in Abaqus.

The evolution of the topology in Figure 2 is consistent with the trend in Figure 3 (top panel). The topology does not essentially vary during the optimization process, according to the almost constant trend of the constraints. On the other hand, the highly oscillatory trend of $\mathcal{M}$ in the first optimization stage is related to the effect of the smoothing and of the sharpening operations, which are confined to the first 24 iterations. From $\mathtt{k}=25$, only the minimization process and the mesh adaptation contribute to a mass variation, with less striking changes.

Finally, in Figure 4 (left), we show the $3\times 3$-cell material generated by a periodic repetition of the optimized unit cell.

3.2. Design Case 2

The second MultiP-microSIMPATY run aims at designing a microstructure that provides high stiffness and thermal conductivity along the x-direction and a high shear stiffness. As for the Design Case 1, these requirements might originate a set of conflicting constraints. In fact, the two former demands are expected to orient the material along the x-direction, while the latter requirement prescribes also the presence of material along the diagonal of the cell Y, which could react by tension to shear loading. This design setting is formalized by problem (9) when completed by the following set of constraints:

$$\left\{\begin{array}{c}0.23\le E_{1111}^H\le 0.35\hfill \\ 0.08\le E_{1212}^H\le 0.15\hfill \\ 0.3\le {\displaystyle \frac{E_{2222}^H}{E_{1111}^H}}\le 2.0\hfill \\ 0.3\le k_{11}^H\le 1.0\hfill \\ 0\le {\displaystyle \frac{k_{22}^H}{k_{11}^H}}\le 2.\hfill \end{array}\right.$$

(16)

We highlight that the bounds for the stiffness tensor components to be promoted, $E_{1111}^H$ and $E_{1212}^H$, are set by taking into account the mass minimization goal, i.e., by keeping them considerably lower than 1.

Algorithm 1 stops in 56 iterations due to mesh stagnation. Figure 5 gathers the density field distribution together with the associated anisotropically adapted computational mesh at iterations $\mathtt{k}=$ 5, 20, and 56. At the fifth iteration, the cell presents very thin struts that are progressively erased by the combined action of the Helmholtz and the Heaviside filters. For $\mathtt{k}=20$, the topology essentially coincides with the final optimized one, although the layout still exhibits intermediate density values along the boundaries. The structure contours become sharper and sharper throughout the next iterations when filtering is switched off and thanks to the mesh adaptation procedure.

Concerning the final topology, we observe that most of the material is aligned along the two main diagonals of Y. This guarantees high shear stiffness, while ensuring a low stiffness along the y-direction, so that the lower bound for $E_{2222}^H/E_{1111}^H$ is reached. On the other side, the requirements on $E_{1111}^H$ and $k_{11}^H$ are taken into account by the two thinner struts along the x-direction, which improves the corresponding stiffness and the thermal conductivity. Figure 4 (center) provides a sketch of the metamaterial associated with the optimized cell in a $3\times 3$ cellular pattern.

For a more quantitative characterization of the optimized structure in terms of mass and reached constraints, we refer to Table 1. We notice that to address the conflict among the several requirements, the optimization process pushes all the constrained quantities toward the lower bound of the corresponding range, while increasing the mass of the structure if compared, for instance, with the previous design case.

3.3. Design Case 3

As a third design, we carry out the optimization of a microcell characterized by similar stiffness and thermal conductivity along the x- and y-directions and by a high shear stiffness. This leads to solve problem (9) when the following constraints are enforced:

$$\left\{\begin{array}{c}0.10\le E_{1111}^H\le 0.15\hfill \\ 0.08\le E_{1212}^H\le 0.10\hfill \\ 1.0\le {\displaystyle \frac{E_{2222}^H}{E_{1111}^H}}\le 1.1\hfill \\ 0.25\le k_{11}^H\le 0.40\hfill \\ 1.0\le {\displaystyle \frac{k_{22}^H}{k_{11}^H}}\le 1.1.\hfill \end{array}\right.$$

(17)

The limited range for the two ratios $E_{2222}^H/E_{1111}^H$ and $k_{22}^H/k_{11}^H$ is consistent with the request for comparable stiffness and thermal conductivities along the two directions, whereas the mass minimization goal justifies the tight variation for the other tensors components.

The MultiP-microSIMPATY algorithm resorts to 35 loops before satisfying the stopping criterion. Figure 6 shows the density field and the mesh for three different global iterations of the algorithm. As for the previous design cases, thin features are removed by filtering during the first 24 iterations, while intermediate densities are erased in the second part of the process by the mesh adaptation procedure. As a consequence, the final microstructure exhibits very sharp density gradients, so that little post-processing has to be applied. In the final layout, most of the material is allocated along the two main diagonals of the domain, which ensures the required high shear stiffness as well as the balance between stiffness and thermal conductivity with respect to the horizontal and vertical directions.

Table 1 offers some additional quantitative information regarding the optimized structure. All the box constraints are satisfied (with a slight violation for the component $E_{1111}^H$) in the presence of a structure mass comparable with the one obtained for Design Case 2 (about $40\%$ with respect to the full material configuration). We refer to Figure 4 (right) for an example of the microcellular material associated with the optimized cell.

4. Discussion of Results

This section is meant to highlight the benefits of the MultiP-microSIMPATY algorithm. To this aim, we compare the layouts provided by the proposed methodology with unit cells available in engineering practice and with cellular materials designed by a standard inverse homogenization procedure, which does not exploit mesh adaptation.

4.1. Comparison with Off-The-Shelf Designs

This first investigation is carried out by comparing each of the three designs in the previous section with state-of-the-art unit cells in terms of mechanical and thermal performance, after setting a reference value for the overall mass. The quantities involved in such a comparison are the homogenized elastic modulus, $E_x^H$ and $E_y^H$, associated with the direction x and y, which coincide with the inverse of the diagonal entries, $C_{11}^H$ and $C_{22}^H$, of the compliance matrix ${\mathbf{C}}^H={\left({\mathbf{E}}^H\right)}^{-1}$; the homogenized shear modulus, $G^H,$ equal to the inverse of the third diagonal entry of matrix ${\mathbf{C}}^H$; the homogenized thermal conductivities, $k_{11}^H$ and $k_{22}^H$, along the x- and y-direction. The results of this analysis are summarized in

Table 2. Comparison between the MultiP-microSIMPATY optimized structures and off-the-shelf designs in terms of homogenized elastic and thermal properties, for comparable volume fraction values.

		${\mathit{E}}_{\mathit{x}}^{\mathit{H}}$	${\mathit{E}}_{\mathit{y}}^{\mathit{H}}$	${\mathit{G}}^{\mathit{H}}$	${\mathit{k}}_{11}^{\mathit{H}}$	${\mathit{k}}_{22}^{\mathit{H}}$
			Design Case 1
D1		0.012	0.015	0.056	0.200	0.113
A		0.009	0.009	0.075	0.163	0.163
B		0.095	0.042	0.059	0.198	0.131
			Design Case 2
D2		0.126	0.039	0.082	0.317	0.126
C		0.341	0.116	0.002	0.432	0.125
			Design Case 3
D3		0.070	0.070	0.082	0.260	0.261
L		0.188	0.188	0.072	0.255	0.255

.

Concerning Design Case 1, we perform two comparisons. Since the geometry provided by MultiP-microSIMPATY is similar to a square cell rotated by 45${}^{\circ }$, we choose simple squares (A and B) characterized by the same rotation as state-of-the-art unit cells. The basic squares in layout A fully couple mechanical and thermal features, thus excluding this cell for the purpose addressed in the first design case. This justifies the selection of cell B where the reinforcing horizontal strut mimics the very thin diagonal member connecting the adjacent sides in the proposed layout (D1). From a structural perspective, the horizontal strut in B increases the nodal connectivity and reacts with tension/compression to a load applied along the x-axis. This fact is confirmed by the non-isotropic elastic behavior of the material (compare the values $E_x^H$ and $E_y^H$). Regarding thermal conduction, the strut promotes heat transfer along the horizontal direction, as highlighted by the discrepancy between $k_{11}^H$ and $k_{22}^H$. In the optimized layout D1, the thin member is instead slightly inclined and does not connect two opposite nodes. Thus, the elastic modulus along the two directions is similar, since the strut reacts by bending to a load applied along the x-axis. Moreover, the thin member promotes the heat transfer along the x-direction, thus decoupling the ratios $E_y^H/E_x^H$ and $k_{22}^H/k_{11}^H$.

The unit cell D2 has been designed to ensure high stiffness and conductivity along the x-direction as well as a high shear modulus. As reference layout, we consider a square cell characterized by a rectangular cavity. This choice offers us a trivial solution to optimize stiffness and conductivity along direction x. The optimization performed by MultiP-microSIMPATY is corroborated by the values of $G^H$. In fact, cell D2 is characterized by a shear modulus, which is approximately 40 times higher when compared with the reference layout, although the values of $E_x^H$ and $k_{11}^H$ for cell D2 are, on average, 30% lower with respect to cell C.

Finally, the Design Case 3 aims at ensuring equal elastic modulus and conductivity along the x- and y-directions, as well as a high shear modulus. The paradigm for an isotropic stretch-based lattice, namely the standard triangular cell (L), is assumed as the off-the-shelf layout. A comparison between the corresponding values in Table 2 shows a 15% increment in the shear modulus of cell D3. In addition, both cells D3 and L exhibit the requested isotropic behavior in terms of the selected mechanical and thermal properties.

4.2. Comparison with Standard Inverse Homogenization

This section is meant to verify the benefits led by mesh adaptation in the context of thermo-elastic inverse homogenization, which is in accordance with the preliminary remarks in Section 3.

To this aim, we carry out a comparison between the MultiP-microSIMPATY algorithm and a standard inverse homogenization procedure. This comparison is performed in terms of mass. We expect that the employment of mesh adaptation leads to efficiently allocate the available material, thus promoting the mass minimization. As a reference standard approach, we implement a non-adaptive version of Algorithm 1, where the adaptation loop (lines 3–12) is replaced by the single call

$${\rho }_h=\mathtt{o}\mathtt{p}\mathtt{t}\mathtt{i}\mathtt{m}\mathtt{i}\mathtt{z}\mathtt{e}(\tilde{\mathcal{J}},\tilde{\mathcal{C}},{\mathbf{c}}^l,{\mathbf{c}}^u,\tilde{\mathcal{G}},{\rho }_h^{\mathtt{0}},\mathtt{T}\mathtt{O}\mathtt{P}\mathtt{T},\mathtt{I}\mathtt{T}).$$

We refer to this variant of Algorithm 1 as MultiP-microSIMP. In this case, the optimization is performed on the filtered density, so that the goal functional, the constraints, and the associated derivatives are modified accordingly (this justifies the new notation $\mathcal{Q}\to \tilde{\mathcal{Q}}$, with $\mathcal{Q}=\mathcal{J},\mathcal{C},\mathcal{G}$, where $\tilde{\mathcal{Q}}$ refers to quantities dependent on the filtered density). This choice is recurrent in topology optimization [61,62]. As far as all the parameters required by the optimization are concerned, we preserve the same values as in Section 2.3, while the computational mesh coincides with a $50\times 50$ structured mesh.

Figure 7 compares the optimized layouts delivered by MultiP-microSIMP (top) and MultiP-microSIMPATY (bottom) for the three design cases in Section 3. The topologies characterizing the three cells vary when resorting to mesh adaptation. In general, MultiP-microSIMPATY provides more complex layouts, which however are still manufacturable. The presence of intermediate densities in the cells yielded by MultiP-microSIMP is highligthed by the blurred structure contours, promoted by the massive employment of filtering.

Table 3. Comparison between the optimized cells delivered by MultiP-microSIMPATY and a standard inverse homogenization algorithm in terms of mass.

	D1	D2	D3
MultiP-microSIMP	0.330	0.443	0.486
MultiP-microSIMPATY	0.292	0.412	0.415
Mass reduction [%]	11.5%	7.0%	14.6%

quantitatively assesses the optimization performance of the two algorithms by collecting the mass of the corresponding unit cells, together with the percentage mass reduction ensured by MultiP-microSIMPATY. On average, a mass saving of approximately 10% is guaranteed by the sharp detection of the material/void interface, i.e., by the removal of intermediate densities.

The use of filtering deserves further discussion. In particular, we prove the redundancy of the filtering phase after a sufficiently large number of global optimization iterations. To this aim, we run Algorithm 1 for $\mathtt{kfmax}=25$ and $\mathtt{kfmax}=\mathtt{kmax}$ (i.e., smoothing and sharpening filters in lines 6–7 are applied at each global iteration). Figure 8 compares the output associated with these two choices. The final topology provided by both the procedures is the same. This confirms that filtering is instrumental only in the identification of the final layout, and this takes place during the first iterations. From the top-left panel, the slightly diffusive action of the selected filtering is also evident, giving rise to intermediate densities along the layout boundaries. On the other hand, the removal of filtering allows mesh adaptation to sharply detect gradients from material to void, thus increasing the quality of the final output (compare the two panels on the left panel). The improvement in terms of boundary detection is confirmed also by the final adapted mesh, which captures the steep gradients of the density with thinner refined areas (compare the two panels on the right).

Finally, we highlight that the presence of blurred interfaces may raise issues in the extraction of the final geometry after the optimization procedure. In fact, the extracted geometry strongly depends on the cut-off threshold, with a possible significant alteration of the overall mass and the expected thermo-elastic properties.

5. Conclusions and Perspectives

In this paper, we provide a new methodology for the design of cellular materials optimized by means of multi-physics inverse homogenization, which was discretized on customized computational meshes. The inverse homogenization problem is modeled by a standard density-based topology optimization at the microscale; the grid is generated by exploiting an anisotropic a posteriori error estimator that drives a mesh adaptation procedure. These two phases are iteratively coupled in the MultiP-microSIMPATY algorithm in order to deliver layouts characterized by clear-cut contours. In particular, the goal of the analyzed test cases is the design of lightweight structures with prescribed elastic and thermal properties, according to a multi-physics framework.

The main results of this work can be outlined as follows:

(i)

The MultiP-microSIMPATY algorithm provides original design solutions, complying also with conflicting requirements;

(ii)

The good performance of microSIMPATY has been confirmed also in a thermo-elastic context. Standard issues typical of topology optimization, such as the presence of intermediate densities, of jagged boundaries, and of too complex structures, is mitigated by the employment of a mesh customized to the design process (see Figure 7 and Table 3);

(iii)

The new cellular materials have been successfully compared with consolidated solutions in terms of mechanical and thermal properties (see Table 2);

(iv)

Filtering can be considerably limited thanks to the use of mesh adaptation. This turns into an improvement in terms of accuracy of the optimization process (see Figure 8);

(v)

The employment of an anisotropic mesh adaptation provides advantages with a view to a manufacturing phase. Indeed, the unit cells designed by MultiP-microSIMPATY exhibit very smooth geometries which demand for a very limited post-processing;

(vi)

The procedure here settled turns out to be fully general with respect to the selected multi-physics context.

Possible future developments include the extension of the MultiP-microSIMPATY design procedure to a 3D setting. The proposed methodology could also be exploited in a multiscale topology optimization framework [38], inspired by the many possible applications in engineering practice (including medicine, aerospace, automotive, and architecture). In such a context, with a view to the manufacturing step, another issue that deserves further investigation is represented by the handling of the transition area between different cellular materials. Finally, innovative techniques, such as model reduction or machine learning, still represent topics of high relevance in topology optimization for a future examination [63,64,65].