# Topology optimization of contact-aided thermo-mechanical regulators

Anna Dalklint[][]footnotemark: [] Joe Alexandersen [][]footnotemark: []
Andreas Henrik Frederiksen[][]footnotemark: [] Konstantinos
Poulios[][]footnotemark: [] Ole Sigmund[][]footnotemark: [] Technical
University of Denmark, Department of Civil and Mechanical Engineering, Nils
Koppels Allé, Building 404, 2800 Kongens Lyngby, Denmark [][]footnotemark: []
University of Southern Denmark, Department of Mechanical and Electrical
Engineering, Odense 5230, Denmark [][]footnotemark: []

###### Abstract

Topology optimization is used to systematically design contact-aided thermo-
mechanical regulators, i.e. components whose effective thermal conductivity is
tunable by mechanical deformation and contact. The thermo-mechanical
interactions are modeled using a fully coupled non-linear thermo-mechanical
finite element framework. To obtain the intricate heat transfer response, the
components leverage self-contact, which is modeled using a third medium
contact method. The effective heat transfer properties of the regulators are
tuned by solving a topology optimization problem using a traditional gradient
based algorithm. Several designs of thermo-mechanical regulators in the form
of switches, diodes and triodes are presented.

###### keywords:

Thermo-mechanical regulators , Thermal components , Coupled non-linear thermo-
mechanical analysis , Third medium contact

##  1 Introduction

Managing heat transfer poses a significant challenge in modern technology.
While this challenge always has existed for conventional heating and cooling
devices like heat exchangers and refrigerators, it has become increasingly
important in modern applications to e.g. enhance the lifespan of electronics,
improve advanced manufacturing processes, minimize energy consumption in high-
performance computing or data storage clusters, and enable certain
functionalities in space applications. Traditionally, the thermal properties
of these systems are regulated using linear and passive thermal components,
whose thermal properties are mainly governed by their static thermal
conductivity and size. However, there is a growing interest in exploring
components that harness switchable and non-linear heat transfer mechanisms to
obtain tunable heat transfer control. This development has lead to the
realization of advanced components such as thermal switches, diode and
triodes, cf. Wehmeyer et al. [35] and Castelli et al. [7]. In this work, we
refer to components that exhibit non-linear heat transfer functions as
_thermal regulators_.

Thermal regulators are useful in various applications. For example, they can
be used to actively regulate the indoor climate of buildings (Miao et al.
[22]), or actively manage heat dissipation rates in e.g. electronics (Yang et
al. [38]) and spacecrafts (Milanez and Mantelli [24]). See review articles by
e.g. Wehmeyer et al. [35] as well as Wong et al. [36] and Ram et al. [25] for
an exhaustive exploration of application possibilities. Actuation mechanisms
and thereby design of the thermal regulators can vary. One approach to the
design of thermal regulators relies on active actuation, i.e. when the
effective conductivity of the regulator is actively controlled via e.g.
applied pressure (Cho et al. [8]), electric fields (Almanza et al. [2]) or
magnetic fields (Bosisio et al. [6]). Another approach relies on passive
actuation, such that the component self-regulates in response to the exterior
temperature changes via thermal expansion/shrinking (Heo et al. [18]) or
material phase changes (Gulfam et al. [17]). Existing designs of regulators
often leverage conduction via moving structural parts, i.e. mechanics, (Cho et
al. [8]), or gas/fluid convection (Yang et al. [38]) to change their effective
thermal conductivity during actuation. In this work, we focus on the former,
i.e. _thermo-mechanical_ regulators.

A thermal regulator should ultimately exhibit a non-linear or switchable heat
transfer function, but the desired shape of this transfer function depends on
the application. For example, a _thermal switch_ exhibits two distinct
functional states: an on state facilitating heat transfer and an off state
where heat transfer is impeded. A _thermal diode_ , on the other hand, should
rectify the heat flow, such that the heat transfer along a specific axis is
facilitated or impeded depending on the sign of the temperature gradient. An
even more complex component is the _thermal triode_ , which, analogously to an
electrical triode/transistor, should be able to amplify an incoming signal. A
natural and efficient way of controlling the effective thermal conductivity in
a thermo-mechanical regulator is to develop designs that utilize structural
contact. In this way, the effective heat conductivity is high when contact
occurs, and otherwise low, cf. the schematic of a contact-aided thermal switch
in Fig. 1.

Fig. 1.: A schematic of a thermal switch leveraging contact.

In the existing literature, thermal regulators are designed based on human
intuition and time-consuming design-and-experimental iterations. This process
can prove cumbersome when designing regulators with highly complex heat
transfer functions. Fortunately, this design process is an obvious candidate
for more efficient systematic design approaches such as topology optimization.
One hurdle to the topology optimization of contact-aided thermal regulators is
however their modeling complexity, caused by 1) the coupled thermo-mechanics
and 2) the intrinsically non-linear contact physics.

The topology optimization of thermo-mechanical problems appears in numerous
works. For example, Sigmund [29, 30] designed thermal actuators undergoing
finite deformations using topology optimization, which later were realized by
Jonsmann et al. [20]. Recently, Granlund et al. [15] considered the topology
optimization of coupled transient non-linear thermo-mechanical problems. A
topology optimization approach to inverse motion based form finding of coupled
static thermo-mechanical problems was investigated by Sui et al. [31]. Chung
et al. [9] used a level-set approach to the topology optimization of
structures undergoing large deformations due to thermal and mechanical loads.

Recently, Bluhm et al. [4] illustrated the power of the so-called third medium
contact method in topology optimization problems of contact-aided mechanisms.
This contact model, pioneered by Wriggers et al. [37], approximates contact
pressures via the reactive tractions that occur when a third medium, that must
be present in between the contacting surfaces, is compressed. Indeed, via
appropriate constitutive modeling, the third medium naturally exists in the
topology optimization formulation in the form of the ersatz material. Bluhm et
al. [4] employed a monolithic approach to the topology optimization problem,
i.e. solved the minimization of the objective and physical residuals
simultaneously using Newton’s method. In a recent work, Frederiksen et al.
[13] instead implemented the approach in the traditional staggered fashion,
and use MMA (Method of Moving Asymptotes, cf. Svanberg [32]) to update the
design. Other examples of the the third medium contact method in topology
optimization include Bluhm et al. [5], whereby non-linear springs are
designed, and Dalklint et al. [10], who design metamaterials that utilized
self-contact to obtain a prescribed macroscopic stress-strain behavior. Very
recently, Faltus et al. [12] extended a third medium contact model to
pneumatically actuated systems.

To summarize; in this work, we use topology optimization to design thermo-
mechanical regulators that harness self-contact to exhibit highly non-linear
heat transfer functions. To this end, we extend the prospect of the third
medium contact method to thermo-mechanical problems, such that we now seek a
third medium that, in addition to contact tractions, transfers heat when
contact occurs, but otherwise approximates an isolating material. The finite
element method is used to model the fully coupled non-linear thermo-mechanical
problem. The optimization problems are solved using the gradient-based
optimizer MMA, and the gradients of the cost and constraint functions are
obtained from an adjoint sensitivity analysis. The numerical results display
various designs of thermo-mechanical switches, diodes and triodes.

##  2 Preliminaries

In our steady-state thermo-mechanical formulation, the reference, undeformed
configuration $\Omega\in\mathbb{R}^{2}$ and the current, deformed
configuration $\Omega_{c}\in\mathbb{R}^{2}$ are distinguished. A material
point $\text{\boldmath$X$}\in\Omega$ is displaced to its deformed, spatial
position $\text{\boldmath$x$}\in\Omega_{c}$ via the displacement field
$\text{\boldmath$u$}(\text{\boldmath$X$}):\Omega\rightarrow\mathbb{R}^{2}$,
such that
$\text{\boldmath$x$}=\text{\boldmath$X$}+\text{\boldmath$u$}(\text{\boldmath$X$})$.
The local deformation is described by the deformation gradient
$\text{\boldmath$F$}=\text{\boldmath$\nabla$}\text{\boldmath$x$}=\text{\boldmath$1$}+\text{\boldmath$\nabla$}\text{\boldmath$u$}$,
where $\nabla$ is the gradient operator with respect to $X$, and
$J=\text{det}(\text{\boldmath$F$})$. The temperature field is denoted
$\theta(\text{\boldmath$X$}):\Omega\rightarrow\mathbb{R}$. The deformation
gradient is assumed to satisfy the multiplicative split
$\text{\boldmath$F$}=\text{\boldmath$F$}^{e}\text{\boldmath$F$}^{\theta}$,
where $\text{\boldmath$F$}^{e}$ and $\text{\boldmath$F$}^{\theta}$ are the
elastic and thermal parts of the deformation gradient, respectively. An
isotropic thermal deformation is assumed, i.e.
$\text{\boldmath$F$}^{\theta}=\left(J^{\theta}\right)^{1/3}\text{\boldmath$1$}$
with
$J^{\theta}=\text{det}(\text{\boldmath$F$}^{\theta})=e^{3\alpha(\theta-\theta_{o})}$,
where $\alpha$ is the thermal expansion coefficient and $\theta_{o}$ is the
reference temperature. For future reference, we denote
$J^{e}=\text{det}(\text{\boldmath$F$}^{e})$.

The undeformed boundary $\partial\Omega$ with normal $n$, is decomposed into
the complementary surfaces $\partial\Omega_{u}$ and $\partial\Omega_{t}$ for
the mechanical and $\partial\Omega_{\theta}$, $\partial\Omega_{q}$ and
$\partial\Omega_{q_{c}}$ for the thermal boundary value problems. For the
mechanical problem, we prescribe the displacement
$\text{\boldmath$u$}=\bar{\text{\boldmath$u$}}$ on $\partial\Omega_{u}$ and
the traction $\text{\boldmath$t$}=\text{\boldmath$0$}$ on
$\partial\Omega_{t}$, whereas for the thermal problem we prescribe the
temperature $\theta=\bar{\theta}$ on $\partial\Omega_{\theta}$, the heat flux
$\text{\boldmath$q$}\cdot\text{\boldmath$n$}=0$ on $\partial\Omega_{q}$ and
Newton convection
$\text{\boldmath$q$}\cdot\text{\boldmath$n$}=h\left(\theta-\theta_{\infty}\right)$
on $\partial\Omega_{q_{c}}$, where $h$ is the heat transfer coefficient and
$\theta_{\infty}$ is the temperature of the surrounding medium. For
simplicity, it is always assumed that
$\text{\boldmath$u$}=\text{\boldmath$0$}$ on $\partial\Omega_{q_{c}}$. The
heat flux in the undeformed configuration follows from Fourier’s law, i.e.

$\text{\boldmath$q$}=-J\text{\boldmath$F$}^{-1}\text{\boldmath$\kappa$}\text{\boldmath$F$}^{-T}\cdot\text{\boldmath$\nabla$}\theta,$
(1)

where $\text{\boldmath$\kappa$}=\kappa\text{\boldmath$1$}$ is the spatial
isotropic heat conductivity tensor and $1$ is the identity.

Following Miehe [23] and Holzapfel [19], we assume reversible, elastic,
mechanical response, such that the weak forms corresponding to the thermo-
mechanical problem are

$\begin{array}[]{ll}\displaystyle\mathcal{R}_{u}(\text{\boldmath$u$},\theta,\text{\boldmath$\lambda$}_{u};\delta\text{\boldmath$u$})&\displaystyle=\int_{\Omega}\Psi_{,u}(\text{\boldmath$u$},\theta;\delta\text{\boldmath$u$})\,dV+\int_{\partial\Omega_{u}}\text{\boldmath$\lambda$}_{u}\cdot\delta\text{\boldmath$u$}\,dS=0,\\\\[10.0pt]
\displaystyle\mathcal{R}_{\theta}(\text{\boldmath$u$},\theta,\lambda_{\theta};\delta\theta)&\displaystyle=\int_{\Omega}\text{\boldmath$\nabla$}\delta\theta\cdot\text{\boldmath$q$}(\text{\boldmath$u$},\theta)\,dV+\int_{\partial\Omega_{\theta}}\lambda_{\theta}\cdot\delta\theta\,dS+\int_{\partial\Omega_{q_{c}}}h\left(\theta-\theta_{\infty}\right)\delta\theta\,dS=0,\\\\[10.0pt]
\displaystyle\mathcal{R}_{\lambda_{u}}(\text{\boldmath$u$},\text{\boldmath$\lambda$}_{u};\delta\text{\boldmath$\lambda$}_{u})&\displaystyle=\int_{\partial\Omega_{u}}\left(\text{\boldmath$u$}-\bar{\text{\boldmath$u$}}\right)\cdot\delta\text{\boldmath$\lambda$}_{u}\,dS=0,\\\\[10.0pt]
\displaystyle\mathcal{R}_{\lambda_{\theta}}(\theta,\lambda_{\theta};\delta\lambda_{\theta})&\displaystyle=\int_{\partial\Omega_{\theta}}\left(\theta-(\theta_{o}+\bar{\theta})\right)\cdot\delta\lambda_{\theta}\,dS=0,\end{array}$
(2)

for all
$(\delta\text{\boldmath$u$},\delta\theta,\delta\text{\boldmath$\lambda$}_{u},\delta\lambda_{\theta})\in(\mathbb{R}^{2}\times\mathbb{R}\times\mathbb{R}^{2}\times\mathbb{R})$,
and where the boundary conditions are applied in a weak sense via the Lagrange
multiplier fields $\text{\boldmath$\lambda$}_{u}\in\partial\Omega_{u}$ and
$\lambda_{\theta}\in\partial\Omega_{\theta}$. It is noted that
$\text{\boldmath$\lambda$}_{u}$ and $\lambda_{\theta}$ represents unknown
reaction tractions and normal heat fluxes, i.e.
$\text{\boldmath$n$}\cdot\text{\boldmath$q$}$, respectively. In the above, the
Kelvin-Joule effect is neglected and vanishing heat sources are assumed within
$\Omega$. We further let $\Psi$ denote the strain energy density, which is
choosen according to the isotropic neo-Hookean type hyperelastic model

$\begin{array}[]{ll}\displaystyle\Psi(\text{\boldmath$u$},\theta)&\displaystyle=\frac{K}{2}\text{ln}(J^{e})^{2}+\frac{G}{2}\left((J^{e})^{-2/3}\left\lVert\text{\boldmath$F$}^{e}\right\rVert^{2}-3\right)\\\\[10.0pt]
&\displaystyle=\frac{K}{2}\left(\text{ln}(J)^{2}-6\alpha(\theta-\theta_{o})\text{ln}(J)+9\alpha^{2}(\theta-\theta_{o})^{2}\right)+\frac{G}{2}\left(J^{-2/3}\left\lVert\text{\boldmath$F$}\right\rVert^{2}-3\right),\end{array}$
(3)

where $K=E/3(1-2\nu)$, $G=E/2(1+\nu)$, $E$ and $\nu$ denote the initial bulk
modulus, shear modulus, Young’s modulus and Poisson’s ratio in the limit of
infinitesimal strains. From the strain energy density in Eq. (3), the first
Piola-Kirchhoff stress tensor is computed as

$\text{\boldmath$P$}(\text{\boldmath$\nabla$}\text{\boldmath$u$})=\frac{\partial\Psi}{\partial\text{\boldmath$F$}}=K\left(\text{ln}(J)\text{\boldmath$F$}^{-T}-3\alpha(\theta-\theta_{o})\text{\boldmath$F$}^{-T}\right)+GJ^{-2/3}\text{dev}(\text{\boldmath$F$}\text{\boldmath$F$}^{T})\text{\boldmath$F$}^{-T},$
(4)

where
$\text{dev}(\text{\boldmath$A$})=\text{\boldmath$A$}-\frac{1}{3}\text{tr}(\text{\boldmath$A$})\text{\boldmath$1$}$
for an arbitrary second order tensor $A$. In this way, we find
$\Psi_{,u}=\text{\boldmath$P$}:\text{\boldmath$\nabla$}\delta\text{\boldmath$u$}$,
cf. Eq. (2). For later convenience, we define

$\mathcal{R}(\text{\boldmath$a$};\delta\text{\boldmath$a$})=\mathcal{R}_{u}(\text{\boldmath$u$},\theta,\text{\boldmath$\lambda$}_{u};\delta\text{\boldmath$u$})+\mathcal{R}_{\theta}(\text{\boldmath$u$},\theta,\lambda_{\theta};\delta\theta)+\mathcal{R}_{\lambda_{u}}(\text{\boldmath$u$},\text{\boldmath$\lambda$}_{u};\delta\text{\boldmath$\lambda$}_{u})+\mathcal{R}_{\lambda_{\theta}}(\theta,\lambda_{\theta};\delta\lambda_{\theta})=0,$
(5)

where
$\text{\boldmath$a$}=(\text{\boldmath$u$},\theta,\text{\boldmath$\lambda$}_{u},\lambda_{\theta})$.

##  3 Model components

###  3.1 Topology optimization

The goal of the topology optimization is to distribute the material in the
design domain $\Omega$, such that an objective function is minimized. To this
end, let $z\in\Omega\rightarrow[0,1]$ denote the design variable volume
fraction field. To regularize the optimization problem, the design space is
restricted via the PDE filter (Lazarov and Sigmund [21]), such that fine scale
oscillations of $z$ are penalized. This provides a smoothed field
$\zeta\in\Omega\rightarrow[0,1]$, as a solution to the following boundary
value problem

$\mathcal{R}_{\zeta}(\zeta,z;\delta\zeta)=\int_{\Omega}l^{2}\text{\boldmath$\nabla$}\zeta\cdot\text{\boldmath$\nabla$}\delta\zeta\,dV+\int_{\Omega}\left(\zeta-z\right)\delta\zeta\,dV=0,$
(6)

for all $\delta\zeta\in\mathbb{R}$. In the above, $l$ is the filter radius.

This filtering invariably produces regions wherein
$\zeta(\text{\boldmath$X$})\in(0,1)$. In order to push the solution towards a
black and white design, thresholding (Guest et al. [16]) is employed via the
approximate Heaviside projection (Wang et al. [34])

$\bar{\zeta}=H_{\beta,\eta}(\zeta)=\frac{\text{tanh}(\beta\eta)+\text{tanh}(\beta(\zeta-\eta))}{\text{tanh}(\beta\eta)+\text{tanh}(\beta(1-\eta))},$
(7)

where $\beta$ and $\eta=0.5$ are numerical parameters that control the slope
and position of the Heaviside function. Additionally, intermediate
$\zeta(\text{\boldmath$X$})\in(0,1)$ regions are implicitly penalized via the
SIMP scheme (Bendsøe [3]), such that the material parameters ultimately are
interpolated as

$\begin{array}[]{cc}\displaystyle
K=\delta_{o}^{E}+(1-\delta_{o}^{E})\bar{\zeta}^{p}K_{o},\ \
G=\delta_{o}^{E}+(1-\delta_{o}^{E})\bar{\zeta}^{p}G_{o},\\\\[10.0pt]
\kappa=\delta_{o}^{\kappa}+(1-\delta_{o}^{\kappa})\bar{\zeta}^{p}\kappa_{o}\ \
\text{and}\ \ K\alpha=\bar{\zeta}^{p}K\alpha_{o}.\end{array}$ (8)

In the above, $K_{o},G_{o},\kappa_{o}$ and $\alpha_{o}$ are the reference
bulk-, shear modulus, thermal conductivity and thermal expansion coefficient,
respectively. Also, $p\geq 0$ is the penalization parameter and
$\delta_{o}^{E},\delta_{o}^{\kappa}\geq 0$ and $\delta_{o}^{\alpha}=0$ define
the ersatz material properties.

###  3.2 Internal contact

In topology optimization, it has recently been suggested to employ the so-
called third medium contact method (Wriggers et al. [37]) to model internal
contact, cf. Bluhm et al. [4], Frederiksen et al. [13], and Dalklint et al.
[10]. In this work, the prospect of the third medium contact method is further
extended to the problem of heat transfer between bodies in contact, such that
the third medium now transfers heat when contact occurs, but otherwise acts as
a thermally isolating material. Analogously to the purely mechanical problem,
it can be shown that this behavior can be obtained exclusively via appropriate
constitutive modeling of the third medium. This effect is readily seen by
examining the definition of the heat flux in Eq. (1). When the third medium is
compressed, $J\rightarrow 0$ and since $\text{\boldmath$F$}^{-1}$ is inversely
proportional to $J$, the heat flux in the normal direction of the surface in
contact $\text{\boldmath$n$}\cdot\text{\boldmath$q$}\rightarrow\infty$.

In all real world thermal contact problems, a so called _thermal contact
resistance_ exists due to irregularities in the contacting surfaces which
impede the heat transfer. This contact resistance will give rise to a
temperature jump across the contact interface at the macroscopic scale.
Experiments show that the thermal contact resistance decreases with increased
pressure, explained by the fact that the increased pressure causes more
surface area to come into contact, increasing the heat transfer.
Unfortunately, there exist no canonical form of the relationship between the
thermal contact resistance, or its inverse; the thermal contact conductance,
and the applied pressure. Indeed, this relationship is usually explored in an
experimental setting, cf. e.g. Rosochowska et al. [27] and Semiatin et al.
[28]. However, there exist numerical models which approximates the thermal
contact resistance. For example, Aalilija et al. [1] illustrated the use of a
diffuse contact interface model to capture the thermal contact resistance. As
we will show in the subsequent numerical examples, the third medium contact
method also naturally captures the effect of this thermal contact resistance
via the hyperelastic modeling of the third medium, i.e. the choice of the
strain energy density function and material parameters.

In this work, traditional finite elements are employed as basis for the
numerical solution of the boundary value problems from Eqs. (2) and (6).
Unfortunately, this numerical analysis is plagued by excessively distorted
finite elements in the third medium, which hinders the convergence. To improve
the deformed element quality, we follow Bluhm et al. [4] and augment the
strain energy density function Eq. (3) with a term that penalizes higher order
strain components, i.e.

$\Psi(\text{\boldmath$u$},\theta)\leftarrow\Psi(\text{\boldmath$u$},\theta)+\frac{k_{r}}{2}\mathbb{H}\text{\boldmath$u$}\,\leavevmode\hbox
to1.26pt{\vbox to6.43pt{\pgfpicture\makeatletter\hbox{\hskip
0.63055pt\lower-0.63055pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{
{}{{}}{}{{{}}{}{}{}{}{}{}{}{}}\pgfsys@beginscope\pgfsys@invoke{
}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@invoke{
}\pgfsys@color@gray@fill{0}\pgfsys@invoke{
}\definecolor[named]{pgffillcolor}{rgb}{0,0,0}\definecolor[named]{pgffillcolor}{rgb}{0,0,0}\pgfsys@color@gray@fill{0}\pgfsys@invoke{
}{}\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@moveto{0.43056pt}{0.0pt}\pgfsys@curveto{0.43056pt}{0.2378pt}{0.2378pt}{0.43056pt}{0.0pt}{0.43056pt}\pgfsys@curveto{-0.2378pt}{0.43056pt}{-0.43056pt}{0.2378pt}{-0.43056pt}{0.0pt}\pgfsys@curveto{-0.43056pt}{-0.2378pt}{-0.2378pt}{-0.43056pt}{0.0pt}{-0.43056pt}\pgfsys@curveto{0.2378pt}{-0.43056pt}{0.43056pt}{-0.2378pt}{0.43056pt}{0.0pt}\pgfsys@closepath\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@fillstroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{{}{}}{{}}{}{{{}}{}{}{}{}{}{}{}{}}\pgfsys@beginscope\pgfsys@invoke{
}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@invoke{
}\pgfsys@color@gray@fill{0}\pgfsys@invoke{
}\definecolor[named]{pgffillcolor}{rgb}{0,0,0}\definecolor[named]{pgffillcolor}{rgb}{0,0,0}\pgfsys@color@gray@fill{0}\pgfsys@invoke{
}{}\pgfsys@moveto{0.0pt}{2.58333pt}\pgfsys@moveto{0.43056pt}{2.58333pt}\pgfsys@curveto{0.43056pt}{2.82112pt}{0.2378pt}{3.01389pt}{0.0pt}{3.01389pt}\pgfsys@curveto{-0.2378pt}{3.01389pt}{-0.43056pt}{2.82112pt}{-0.43056pt}{2.58333pt}\pgfsys@curveto{-0.43056pt}{2.34554pt}{-0.2378pt}{2.15277pt}{0.0pt}{2.15277pt}\pgfsys@curveto{0.2378pt}{2.15277pt}{0.43056pt}{2.34554pt}{0.43056pt}{2.58333pt}\pgfsys@closepath\pgfsys@moveto{0.0pt}{2.58333pt}\pgfsys@fillstroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{{}{}}{{}}{}{{{}}{}{}{}{}{}{}{}{}}\pgfsys@beginscope\pgfsys@invoke{
}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@invoke{
}\pgfsys@color@gray@fill{0}\pgfsys@invoke{
}\definecolor[named]{pgffillcolor}{rgb}{0,0,0}\definecolor[named]{pgffillcolor}{rgb}{0,0,0}\pgfsys@color@gray@fill{0}\pgfsys@invoke{
}{}\pgfsys@moveto{0.0pt}{5.16666pt}\pgfsys@moveto{0.43056pt}{5.16666pt}\pgfsys@curveto{0.43056pt}{5.40445pt}{0.2378pt}{5.59721pt}{0.0pt}{5.59721pt}\pgfsys@curveto{-0.2378pt}{5.59721pt}{-0.43056pt}{5.40445pt}{-0.43056pt}{5.16666pt}\pgfsys@curveto{-0.43056pt}{4.92886pt}{-0.2378pt}{4.7361pt}{0.0pt}{4.7361pt}\pgfsys@curveto{0.2378pt}{4.7361pt}{0.43056pt}{4.92886pt}{0.43056pt}{5.16666pt}\pgfsys@closepath\pgfsys@moveto{0.0pt}{5.16666pt}\pgfsys@fillstroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\,\mathbb{H}\text{\boldmath$u$},$
(9)

where $\mathbb{H}$ is the Hessian operator. Following Bluhm et al. [4];
$k_{r}=\bar{k}_{r}H^{2}K_{o}$, where $\bar{k}_{r}$ is a small numerical
parameter and $H$ is the characteristic height of the domain $\Omega$.

###  3.3 Numerical solution procedure

The fully coupled thermo-mechanical system Eq. (5) is solved using Newton’s
method, i.e. in each step the full system Eq. (5) is linearized to obtain

$\displaystyle\mathcal{{R}}(\text{\boldmath$a$}+d\text{\boldmath$a$};\delta\text{\boldmath$a$})\approx\displaystyle\mathcal{{R}}(\text{\boldmath$a$};\delta\text{\boldmath$a$})+\mathcal{{R}}_{,a}(\text{\boldmath$a$};\delta\text{\boldmath$a$},d\text{\boldmath$a$})={0},$
(10)

which is iteratively solved for $d\text{\boldmath$a$}$, such that
$\text{\boldmath$a$}\leftarrow\text{\boldmath$a$}+d\text{\boldmath$a$}$.
Without loss of generality, we limit the numerical examples to plane strain
two dimensional problems. The presented thermo-mechanical topology
optimization framework is implemented in Python where all tangent matrices and
residual vectors are computed via the finite element library GetFEM (Renard
and Poulios [26]), which performs all underlying computations in C++. Design
updates are based on a Python implementation (Deetman [11]) of the Method of
Moving Asymptotes (MMA) (Svanberg [32, 33]).

###  3.4 Thermal contact

To illustrate how the third medium contact method facilitates thermal contact,
the simple rod problem depicted in Fig. 2 is considered. The rod is split into
three parts of lengths $L^{s}=L^{v}=L/3$, such that two solid parts
$\Omega_{s}$ with $E=1$ MPa, $\nu=0.4$, and $\kappa=10$ W/mK are separated by
a third medium region $\Omega_{v}$ with $E_{v}=\delta_{o}^{E}E=10^{-6}E$ and
$\kappa_{v}=\delta_{o}^{\kappa}\kappa$, cf. Fig. 2. We investigate different
values of $\delta_{o}^{\kappa}$ in the following analysis. For simplicity we
assign $\alpha=0$ and focus only on the mechanical effect on the thermal
problem, cf. Eq. (3). We set $h=1$ kW/m2K and $\theta_{\infty}=20$ $\degree$C.
Note that that $\partial\Omega_{u}$ is decomposed into the two complementary
sets $\partial\Omega_{u_{1}}$ and $\partial\Omega_{u_{2}}$.

Fig. 2.: The rod with dimensions $(L,H)$ = ($6$ cm $\times$ $1$ cm).

To obtain self-contact in the rod, the Dirichlet condition
$\text{\boldmath$u$}=(\bar{u},0)$ is applied on $\partial\Omega_{{u}_{2}}$,
where $\bar{u}=-0.5L$, cf. Eq. (2). On $\partial\Omega_{{u}_{1}}$, the
pointwise value of $u_{x}$ and the average of $u_{y}$ are set to zero,. Also,
we assign $\bar{\theta}=200$ $\degree$C. The rod is discretized by $120\times
20$ quadratic quadrilateral finite elements, and nine point Gauss-Lobatto
integration is used for improved stability in the third medium (Faltus et al.
[12]). The displacement $u$ and temperature $\theta$ fields are discretized by
incomplete Q2 elements, whereas the design variable $z$ and filtered $\zeta$
fields are discretized by Q0 and Q1 elements, respectively. The numerical
solution is compared to an analytical solution derived from the 1D
approximation of the thermal contact rod problem, which is obtained by solving
the 1D heat equation

$-\frac{d}{dX}\left(\kappa\frac{d\theta(X)}{dX}\right)=0,\quad X\in\Omega,$
(11)

over the deformed rod after contact has been initiated, i.e.
$\Omega=[0,L_{c}]$ where $L_{c}(\bar{u})=L+\bar{u}$ is the deformed length of
the rod. Equivalently to the 2D problem, Newton convection is present on the
left hand side, i.e. $q(0)=h\left(\theta_{\infty}-\theta(0)\right)$, and a
Dirichlet conditions is applied to the right hand side, i.e.
$\theta(L_{c})=\bar{\theta}$. Additionally, the thermal contact resistance is
accommodated via the temperature jump
$[[\theta(X=L_{c}/2)]]=-R_{th}q(L_{c}/2)$ at $X=L_{c}/2$, where $R_{th}$ is
the thermal contact resistance and $q=-\kappa\frac{d\theta}{dX}$ is the 1D
heat flux. We solve Eq. (11) over $\Omega$ by integrating twice over the two
separate regions $0\leq X\leq L_{c}/2$ and $L_{c}/2<X\leq L_{c}$, and utilize
the boundary conditions in combination with the aforementioned jump condition
to couple the domain solutions. The final 1D temperature distribution becomes

$\theta(X)=\begin{cases}-\frac{1}{\kappa}(AX+B_{1})\ \ \text{for}\ \ 0\leq
X\leq L_{c}/2,\\\\[10.0pt] -\frac{1}{\kappa}(AX+B_{2})\ \ \text{for}\ \
L_{c}/2<X\leq L_{c},\end{cases}$ (12)

where $A=\frac{\theta_{\infty}-\bar{\theta}}{L_{c}/\kappa+R_{th}+1/h}$,
$B_{1}=\kappa(\frac{A}{h}-\theta_{\infty})$ and
$B_{2}=-\kappa\bar{\theta}-AL_{c}$.

First, a sharp contact interface between $\Omega_{s}$ and $\Omega_{v}$ is
analyzed, i.e. the PDE filter is neglected. To investigate how the heat
conductivity in the middle void region, $\kappa_{v}$, affects the heat
transfer during contact, the average normal heat flux
$\frac{1}{|\partial\Omega_{\theta}|}\int_{\Omega_{\theta}}\lambda_{\theta}\,dS$
is plotted versus the applied displacement $\bar{u}\in[0,-0.5L]$ for varying
$\delta_{o}^{\kappa}$ in Fig. 3, where $|\partial\Omega_{\theta}|$ denotes the
area. The results are compared to the analytical solution in Eq. (12) with
perfect contact conditions, i.e. when $R_{th}=0$, and different levels of
deformation, i.e. deformed lengths $L_{c}$.

Fig. 3.: The analytical and numerical average normal reaction fluxes versus
the applied displacement for varying $\delta_{o}^{\kappa}$.

It is noted that that a large premature heat transfer occurs for large values
of $\delta_{o}^{\kappa}>10^{-3}$, cf. Fig. 3. To reduce this premature heat
transfer to acceptable levels, smaller values of $\delta_{o}^{\kappa}\leq
10^{-3}$ have to be considered. However, the simulations with
$\delta_{o}^{\kappa}\leq 10^{-3}$ depicts a drastically reduced final heat
flux magnitude in comparison to the theoretical maximum provided by the
analytical solution. Fortunately, this phenomena proves useful in our
analyzes, since it lets us simulate the aforementioned thermal contact
resistance. Similarly to Aalilija et al. [1], Fig. 3 shows that the magnitude
of the thermal contact resistance can be changed by tuning
$\delta_{o}^{\kappa}$.

Indeed, the 1D thermal contact resistance $R_{th}$ of our sharp interface
third medium can be analytically approximated as

$R_{th}=\frac{L_{c}^{v}}{\delta_{o}^{\kappa}\kappa},$ (13)

where $L_{c}^{v}$ is the deformed length of the third medium, i.e. void
region. In Fig. 4, the analytical and numerical temperature distributions in
the rod are compared for different thermal contact resistances, when
$L_{c}\approx 0.4$ m and $R_{th}$ is obtained from Eq. (13) with $L_{c}^{v}$
computed from the numerical solution, cf. Eq. (12). The results illustrate
good correspondence between the analytical and numerical solutions.

Fig. 4.: The analytical and numerical temperature distributions over the
$L_{c}\approx 0.04$ m rod for varying thermal contact resistances, i.e.
varying $\delta_{o}^{\kappa}$.

Next, the implicit relationship between the thermal contact conductance and
the applied pressure is investigated, cf. Eq. (13). The obtained relationship
is illustrated in Fig. 5 for varying $\delta_{o}^{\kappa}$.

Fig. 5.: The thermal contact conductance, $1/R_{th}$, versus the applied
pressure for varying $\delta_{o}^{\kappa}$.

The numerical investigation of the rod is concluded by considering a diffuse
contact interface via the introduction of the PDE filter, cf. Eq. (6). In Fig.
6 the average normal heat flux on $\partial\Omega_{\theta}$ is plotted versus
the applied displacement $\bar{u}$ for the sharp and diffuse interfaces with
varying filter radii and $\delta_{o}^{\kappa}$. As expected, it is confirmed
that the introduction of the diffuse interface increases the thermal contact
resistance.

Fig. 6.: The reaction flux versus the applied displacement for sharp and
diffuse interfaces with varying filter radii and $\delta_{o}^{\kappa}$.

Finally, the undeformed/deformed rods with sharp/diffuse interfaces for
$\delta_{o}^{\kappa}=10^{-4}$ are depicted in Fig. 7. A visual inspection of
the temperature distributions shown in the bottom row of Fig. 7 confirms the
increased contact resistance for the rods with diffuse interface. It is noted
that large transverse displacements occur at the contact interface when
introducing the filter, which are believed to be caused by the reduced
stiffness of the intermediate volume fraction regions (cf. the SIMP scheme in
Eq. (8)).

Fig. 7.: The undeformed/deformed rods for the sharp/diffuse interface
representations.

Based on the above numerical study we conclude that; 1) the third medium
contact method is capable of modeling contact with a non-vanishing thermal
contact resistance, and 2) the magnitude of the thermal contact resistance can
be controlled via the ersatz material parameters; more specifically
$\delta_{o}^{\kappa}$, and the gap size; i.e. $L_{c}^{v}$ in the 1D model. In
other words, $\delta_{o}^{\kappa}$ must not solely be viewed as a numerical
parameter, but must be calibrated based on experiments. In this work, we
nevertheless limit ourselves to numerical studies, wherefore
$\delta_{o}^{\kappa}\in\\{10^{-4},10^{-3}\\}$ is used, without loss of
generality. It is also concluded that the diffuse contact interfaces that
appear due to the filter increase the thermal contact resistance. However, in
our optimization, a continuation scheme for the thresholding (cf. Eq. (7)) is
used, wherefore we expect the final design to be approximately binary, i.e.
the above 2) conclusion still holds.

###  3.5 Optimization problem

The optimization problems solved can be written in the following format

$\mathbb{P}\ \begin{cases}\begin{array}[]{ll}\underset{z}{\text{min}}\
C({\zeta},\text{\boldmath$a$}),\\\\[8.0pt]
\text{s.t.}\begin{cases}{\mathcal{R}}(\text{\boldmath$a$};\delta\text{\boldmath$a$})={0},\\\\[5.0pt]
\mathcal{R}_{\zeta}(\zeta,z;\delta\zeta)={0},\\\\[5.0pt] \displaystyle
g=\int_{\Omega}H_{\beta,\eta_{d}}(\zeta)\,dV-V^{*}\leq 0,\\\\[8.0pt] 0\leq
z\leq 1,\end{cases}\end{array}\end{cases}$ (14)

where we define the objective function $C$ and an upper bound volume
constraint $g$, ensuring lightweight designs. In this work, thin features are
penalized in the designs via the inclusion of the volume constraint on the
dilated design field $H_{\beta,\eta_{d}}(\zeta)$, where $\eta_{d}=0.4$ (cf.
Frederiksen et al. [13]).

####  3.5.1 Sensitivity analsysis

The topology optimization problem is solved in a staggered fashion using the
adjoint method, wherefore in each design iteration, Eq. (2) is solved exactly.
In the adjoint method, the augmented objective function

$\begin{array}[]{ll}\displaystyle
C^{*}({\zeta},\text{\boldmath$a$})=&\displaystyle
C({\zeta},\text{\boldmath$a$})+\mathcal{{R}}_{a}({\zeta},\text{\boldmath$a$};\text{\boldmath$\mu$}_{a})+\mathcal{R}_{\zeta}(\zeta,z;\mu_{\zeta}),\end{array}$
(15)

is defined, where
$\text{\boldmath$\mu$}_{a}=(\text{\boldmath$\mu$}_{u},\mu_{\theta},\mu_{\lambda_{u}},\mu_{\lambda_{\theta}})$
is the adjoint vector. Taking the first variation with respect to the design
variable field $z$ renders

$\begin{array}[]{ll}\delta
C^{*}=&\mathcal{R}_{\zeta,z}(\zeta,z;\mu_{\zeta},\delta z)\\\\[10.0pt]
&+C^{*}_{,\zeta}+\mathcal{{R}}_{,\zeta}({\zeta},\text{\boldmath$a$};\text{\boldmath$\mu$}_{a},\delta\zeta)+\mathcal{R}_{\zeta,\zeta}(\zeta,z;\mu_{\zeta},\delta\zeta)\\\\[10.0pt]
&+C^{*}_{,a}+\mathcal{{R}}_{,a}({\zeta},\text{\boldmath$a$};\text{\boldmath$\mu$}_{a},\delta\text{\boldmath$a$}),\end{array}$
(16)

and to annihilate the implicit variation $\delta\text{\boldmath$a$}$,
$\text{\boldmath$\mu$}_{a}$ is assigned such that

$C^{*}_{,a}+\mathcal{{R}}_{,a}({\zeta},\text{\boldmath$a$};\text{\boldmath$\mu$}_{a},\delta\text{\boldmath$a$})=\text{\boldmath$0$}.$
(17)

Next, the implicit variation $\delta\zeta$ is annihilated by assign
$\mu_{\zeta}$ as the solution to

$C^{*}_{,\zeta}+\mathcal{{R}}_{,\zeta}({\zeta},\text{\boldmath$a$};\text{\boldmath$\mu$}_{a},\delta\zeta)+\mathcal{R}_{\zeta,\zeta}(\zeta,z;\mu_{\zeta},\delta\zeta)=0.$
(18)

Finally, the sensitivity expression reduces to

$\delta C^{*}=\mathcal{R}_{\zeta,z}(\zeta,z;\mu_{\zeta},\delta z).$ (19)

The sensitivity of the volume constraint $g$ is trivial and therefore omitted
in this presentation.

##  4 Numerical examples

In the following numerical examples, thermal regulators that actuate solely
via thermal expansion/shrinking are designed. Therefore, we assign a non-zero
thermal expansion coefficient $\alpha$ and only apply homogeneous Dirichlet
conditions on $\partial\Omega_{u}$. The same finite element spaces as
presented in the previous example is used to discretize the fields. The
material parameters used in the subsequent examples are stated in Tab. I.

Young’s Modulus $E$ [MPa] | $1$
---|---
Poisson’s ratio $\nu$ | $0.4$
Thermal conductivity $\kappa$ [W/mK] | $10$
Thermal expansion coefficient $\alpha$ [1/K] | $10^{-4}$
Reference temperature $\theta_{o}$ [°C] | $20$
Tab. I.: The material parameters.

The Heavside threshold parameter $\beta$ is initially set to
$\beta^{\text{ini}}$, and is after 300 design iterations increased
$\beta\leftarrow 2\cdot\beta$ every 50th design iteration until
$\beta^{\text{max}}$ is reached. Optimization and mesh parameters that are
common to all subsequent examples are presented in Tab. II.

Maximum volume $V^{*}$ | $0.4|\Omega|$
---|---
Mesh size | $160\times 80$
Element side length $l_{e}$ [m] | 0.0625
Filter radius $l$ [m] | $2l_{e}/2\sqrt{3}$
Ersatz material scaling $\delta_{o}^{E}$ | $10^{-6}$
Void regularization weight $\bar{k}_{r}$ | $10^{-6}$
Penalty parameter $p$ | 3
Tab. II.: The optimization and mesh parameters.

###  4.1 Thermal switch

In the first example, contact-aided thermal switches are designed. The design
domain contains an anode ($\partial\Omega_{\bar{\theta}_{1}}$) and a cathode
($\partial\Omega_{\bar{\theta}_{2}}$), at which the temperature is controlled,
cf. Fig. 8. The temperature at the anode is varied, whereas the temperature at
the cathode is kept constant. To mimic the behavior of a thermal switch, the
heat transfer between the anode and the cathode should vanish until a certain
critical temperature difference is reached, after which the heat transfer
should be facilitated. All boundaries without specific labels are part of
$\partial\Omega_{q}$ and $\partial\Omega_{t}$, i.e. are isolated and traction
free. Note that $\partial\Omega_{q_{c}}=\emptyset$ and that
$\partial\Omega_{\theta}$ is decomposed into the two complementary sets
$\partial\Omega_{{\theta}_{1}}$ and $\partial\Omega_{{\theta}_{2}}$ over which
$\bar{\theta}_{1}$ and $\bar{\theta}_{2}$ are prescribed, respectively. As
depicted in Fig. 8, a non-uniform volume fraction distribution with a centered
gap is utilized as initial guess. By moving the initial position of this gap,
similar albeit shifted designs are obtained, which indicates that we are
dealing with a highly non-convex optimization problem. It has been found that
a uniform initial design can also produce meaningful designs; however, the
convergence of the optimization process is significantly hindered.

Fig. 8.: The first design domain. The design variable field $z$ is prescribed
$z=1$ in the black regions throughout the optimization. The gray regions are
initially assigned $z=0.3$, whereas the white region is $z=10^{-4}$. All
dimensions are in millimeters.

To design thermal switches, the objective function

$C=\int_{\partial\Omega_{\bar{\theta}_{1}}}\left(-w_{1}\lambda_{\theta}^{(1)}+w_{2}\lambda_{\theta}^{(2)}\right)\,dS,$
(20)

is employed, where $(\,\cdot\,)^{(i)}$ denotes quantities evaluated at target
point $i=1,2$ and $w_{i}$ are weights. Since the direction of the temperature
gradient is known, the above objective function promotes vanishing heat flux
at the first target point ($i=1$) and a large (negative) heat flux at the
second target point ($i=2$). The parameters for the switch example are
summarized in Tab. IV.

Threshold parameters $(\beta^{\text{ini}},\beta^{\text{max}})$ | $(2,16)$
---|---
Target temperatures $(\bar{\theta}_{1}^{(1)},\bar{\theta}_{2}^{(1)})$ [$\degree$C] | $(200,0)$
Target temperatures $(\bar{\theta}_{1}^{(2)},\bar{\theta}_{2}^{(2)})$ [$\degree$C] | $(400,0)$
Weights $(w_{1},w_{2})$ | $(2\cdot 10^{3},10^{3})$
Tab. III.: The numerical parameters for the switch example, unless stated
otherwise. Note that only half the domain in Fig. 8 is discretized due to
symmetry.

In Fig. 9, two optimized thermal switch designs with different thermal contact
resistances, i.e. $\delta_{o}^{\kappa}$, are depicted. A glance at the designs
in their undeformed configurations reveals that the
$\delta_{i}^{\kappa}=10^{-4}$ design has a smaller initial gap between the
contact surfaces than the $\delta_{i}^{\kappa}=10^{-3}$ design. This is
expected since the thermal contact resistance of the
$\delta_{i}^{\kappa}=10^{-4}$ design is larger, i.e. it isolates better at the
first target point, and must exert a greater contact pressure than the
$\delta_{i}^{\kappa}=10^{-3}$ design to reach the same heat flux at the second
target point.

Fig. 9.: The thermal switch designs in their undeformed and deformed
configurations for $\delta_{o}^{\kappa}=10^{-3}$ (top row) and
$\delta_{o}^{\kappa}=10^{-4}$ (bottom row).

In Fig. 10 the temperature distributions in the deformed configurations are
shown for the Fig. 9 designs. Again, the effect of the thermal contact
resistances is noted. Indeed, the $\delta_{o}^{\kappa}=10^{-3}$ design
exhibits a smoother temperature gradient across the contact interface when
contact occurs, compared to the $\delta_{o}^{\kappa}=10^{-4}$ design.

Fig. 10.: The temperature distributions over the Fig. 9 designs for
$\delta_{o}^{\kappa}=10^{-3}$ (top row) and $\delta_{o}^{\kappa}=10^{-4}$
(bottom row).

The boundaries of the Fig. 9 designs are poorly resolved, which can adversely
affect the performance of the realized designs. Therefore, additional post-
processing analyzes are conducted, in which the performance of the structured
mesh is compared to body-fitted mesh representations. The body-fitted meshes
in Fig. 11 have been generated using GMSH (Geuzaine and Remacle [14]), based
on the $\zeta\approx 0.5$ contour line. The fitted meshes consist of Q2 and P2
elements which are integrated using nine and seven Gauss-points, respectively.
Again, for ease of implementation, the third medium contact method is used to
model the thermo-mechanical contact, wherefore both the solid and void regions
are meshed.

Fig. 11.: The body-fitted meshes of the Fig. 9 designs with
$\delta_{o}^{\kappa}=10^{-3}$ (left) and $\delta_{o}^{\kappa}=10^{-4}$
(right).

In Fig. 12, the average normal heat flux versus the applied temperature is
plotted for the Fig. 9 designs discretized on structured and body-fitted
meshes. As expected, the $\delta_{o}^{\kappa}=10^{-4}$ design exhibits less
heat flux at both target points compared to the $\delta_{o}^{\kappa}=10^{-3}$
design. Seemingly, there exists a trade-off between good isolating properties
(smaller $\delta_{o}^{\kappa}$), and good heat transfer properties (larger
$\delta_{o}^{\kappa}$). It is also noted that the designs exhibit premature
heat transfer, i.e. contact, on the structured mesh, in relation to the body-
fitted mesh. This is in line with the results from the rod problem when using
discrete or diffuse contact interfaces, cf. Fig. 6. Indeed, it is noted that
after a critical temperature of $\bar{\theta}\approx 500$ $\degree$C is
reached, the heat flux of the $\delta_{o}^{\kappa}=10^{-3}$ body-fitted design
exceeds its structured mesh counterpart, similar to the observed behavior of
the rod problem. The importance of the interface resolution was also
scrutinized in Bluhm et al. [4] and Dalklint et al. [10], who came to similar
conclusions. The $\delta_{o}^{\kappa}=10^{-4}$ design doesn’t exhibit the same
behavior at the applied temperature range due to the presence of a softening
instability mode, leading to bending of the right hand side top and bottom
structural members, effectively reducing the contact pressure and thereby heat
flux. To avoid instabilities at operating temperature conditions we could
follow Frederiksen et al. [13] and penalize small tangent stiffness values, or
introduce buckling constraints via a linearized buckling analysis, but this is
deemed outside the scope of this paper.

Fig. 12.: The average normal heat flux versus the applied temperature for the
Fig. 9 designs discretized on structured and body-fitted meshes. The dashed
black lines highlights the target values $\bar{\theta}_{1}\in\\{200,400\\}$
$\degree$C.

In Fig. 13 the temperature distributions over the deformed designs discretized
by body-fitted meshes are depicted. Qualitatively, the behavior of Fig. 10 is
confirmed, i.e. the $\delta_{o}^{\kappa}=10^{-3}$ design exhibits a smoother
temperature gradient across the contact interface than the
$\delta_{o}^{\kappa}=10^{-4}$ design, as expected.

Fig. 13.: The temperature distributions over the Fig. 11 designs for
$\delta_{o}^{\kappa}=10^{-3}$ (top row) and $\delta_{o}^{\kappa}=10^{-4}$
(bottom row) in their deformed configurations.

Based on this post-processing analysis, is is confirmed that there exist some
discrepancies between the analyzed and realized design. However, their extent
is limited in the sense that the general behavior of the design remains.
Therefore, for the following design examples, we refrain from conducting the
similar post-processing analysis.

###  4.2 Thermal diode

In the second example, a contact-aided thermal diode is designed. Again, the
design domain is that of Fig. 8, but now two load cycles are performed: 1)
$\partial\Omega_{\bar{\theta}_{1}}$ is heated whereas the temperature at
$\partial\Omega_{\bar{\theta}_{2}}$ constant, 2) the polarity is switched such
that $\partial\Omega_{\bar{\theta}_{1}}$ has constant temperature and
$\partial\Omega_{\bar{\theta}_{2}}$ is heated. Our diode should rectify the
heat current, such that the heat transfer between the terminals is facilitated
during load cycle 1), whereas it is impeded during load cycle 2). To this end,
the objective function

$C=w_{1}\int_{\partial\Omega_{\bar{\theta}_{1}}}\left(\lambda_{\theta}^{(1)}+\lambda_{\theta}^{(2)}\right)\,dS,$
(21)

is chosen, which promotes a large (negative) heat flux at the first target
point ($i=1$), i.e. during load cycle 1), and a vanishing heat flux at the
second target point ($i=2$), i.e. during load cycle 2). The parameters for the
diode example are summarized in Tab. IV.

Threshold parameters $(\beta^{\text{ini}},\beta^{\text{max}})$ | $(2,16)$
---|---
Ersatz material scaling $\delta_{o}^{\kappa}$ | $10^{-3}$
Target temperatures $(\bar{\theta}_{1}^{(1)},\bar{\theta}_{2}^{(1)})$ [$\degree$C] | $(400,0)$
Target temperatures $(\bar{\theta}_{1}^{(2)},\bar{\theta}_{2}^{(2)})$ [$\degree$C] | $(0,400)$
Weight $w_{1}$ | $10^{3}$
Tab. IV.: The numerical parameters for the diode example. Note that only half
the domain in Fig. 8 is discretized due to symmetry.

An optimized thermal diode design is depicted in Fig. 14. As envisioned, the
thermal expansion renders a contact connection between the terminals during
the load cycle 1), which facilitates the heat transfer. On the other hand, the
design avoids contact during load cycle 2) to minimize the heat transfer.

Fig. 14.: The thermal diode design in its undeformed and deformed
configurations.

The temperature distributions in Fig. 15 confirm a larger heat transfer
between the terminals in load cycle 1), compared to load cycle 2).

Fig. 15.: The temperature distributions over the Fig. 14 design in its
deformed configurations.

This behavior is again confirmed in Fig. 16, where the average normal heat
flux on $\partial\Omega_{\bar{\theta}_{1}}$ is plotted versus the temperature
difference between $\partial\Omega_{\bar{\theta}_{1}}$ and
$\partial\Omega_{\bar{\theta}_{2}}$. Clearly, our diode design exhibits larger
heat transfer when the temperature difference is positive, than when negative.

Fig. 16.: The average normal heat flux versus the difference in temperature on
$\partial\Omega_{\bar{\theta}_{1}}$ and $\partial\Omega_{\bar{\theta}_{2}}$
for the Fig. 14 design.

###  4.3 Thermal triode

In the last example, we design contact-aided thermal triodes. This design
domain comprises three terminals: one anode
($\partial\Omega_{\bar{\theta}_{1}}$), one cathode
($\partial\Omega_{\bar{\theta}_{2}}$) and one gate
($\partial\Omega_{\bar{\theta}_{3}}$), cf. Fig. 17. There exists a constant
large temperature difference between the anode and the cathode, but the
variable temperature at the gate dictates whether heat transfer between the
anode and cathode is facilitated or impeded. Again,
$\partial\Omega_{q_{c}}=\emptyset$ and all boundaries without specific labels
are part of $\partial\Omega_{q}$ and $\partial\Omega_{t}$, i.e. are isolated
and traction free. In this example, $\partial\Omega_{\theta}$ is decomposed
into the three complementary sets $\partial\Omega_{\bar{\theta}_{1}}$,
$\partial\Omega_{\bar{\theta}_{2}}$ and $\partial\Omega_{\bar{\theta}_{3}}$
over which $\bar{\theta}_{1}$, $\bar{\theta}_{2}$ and $\bar{\theta}_{3}$ are
prescribed, respectively. In this example, a uniform initial volume fraction
distribution is utilized.

Fig. 17.: The second design domain. The design variable field $z$ is
prescribed $z=1$ in the black regions throughout the optimization. The gray
region is initially assigned $z=0.3$. All dimensions are in millimeters.

To design the thermal triodes, the objective function

$\begin{array}[]{ll}\displaystyle
C=&\displaystyle\frac{w_{1}}{|\partial\Omega_{\bar{\theta}_{1}}|}\int_{\partial\Omega_{\bar{\theta}_{1}}}\left(-\lambda_{\theta}^{(1)}+\lambda_{\theta}^{(2)}\right)\,dS+\frac{w_{2}}{|\partial\Omega_{\bar{\theta}_{3}}|}\int_{\partial\Omega_{\bar{\theta}_{3}}}\left(\left(\lambda_{\theta}^{(1)}\right)^{2}+\left(\lambda_{\theta}^{(2)}\right)^{2}\right)\,dS,\end{array}$
(22)

is chosen. This objective function promotes vanishing heat transfer at the
first target point ($i=1$), and a large heat transfer between the anode and
the cathode at the second target point ($i=2$), whilst also promoting a
vanishing heat flux at the gate. The minus sign is introduced in the first
term of Eq. (22) since the heat flux in the normal direction will be negative
at the anode, whereas the third and fourth terms of Eq. (22) are squared since
the sign of the normal heat flux at the gate is unknown. Using similar
arguments, a positive sign is assigned to the second term of Eq. (22). The
influence of the penalty of the average normal heat flux at the gate is
investigated by designing two triodes. The first design (design 1) is obtained
using a weight $w_{2}>0$, whereas the second design (design 2) is obtained
when $w_{2}=0$. The parameters for the triode examples are summarized in Tab.
V.

Threshold parameters $(\beta^{\text{ini}},\beta^{\text{max}})$ | $(2,8)$
---|---
Ersatz material scaling $\delta_{o}^{\kappa}$ | $10^{-3}$
Target temperatures $(\bar{\theta}_{1}^{(1)},\bar{\theta}_{2}^{(1)},\bar{\theta}_{3}^{(1)})$ [$\degree$C] | $(380,50,0)$
Target temperatures $(\bar{\theta}_{1}^{(2)},\bar{\theta}_{2}^{(2)},\bar{\theta}_{3}^{(2)})$ [$\degree$C] | $(380,100,0)$
Weights design 1 $(w_{1},w_{2})$ | $(10^{3},10^{3})$
Weights design 2 $(w_{1},w_{2})$ | $(10^{3},0)$
Tab. V.: The numerical parameters for the triode examples.

The thermal triode designs are depicted in Fig. 18 in their undeformed and
deformed configurations. It is seen that as the temperature at the gate is
increased, the thermal expansion renders a connection between the anode and
the cathode via structural contact, such that the heat transfer is maximized.
It is noted that the gap to the anode is smaller in design 2, compared to
design 1. In Fig. 19, the temperature distributions are depicted in the
deformed configurations.

Fig. 18.: The thermal triode designs in their undeformed and deformed
configurations. Fig. 19.: The temperature distributions over the Fig. 18
designs in their deformed configurations.

Finally, the total heat transfer
$Q_{j}=\int_{\Omega_{\bar{\theta}_{j}}}\lambda_{\theta}\,dS$, at the anode
($j=1$), gate ($j=3$) and cathode ($j=2$) are plotted versus the applied
temperature in Fig. 20. It is confirmed that the optimization renders a design
that can amplify the control signal, i.e. the absolute value of the total heat
transfer at the gate is smaller than that of the cathode. We conclude that the
penalty to the average normal heat flux at the gate increases the efficiency
of the triode, in the sense that less heat is lost through the gate for design
1 than for design 2. It is also noted that the total heat transfer at the gate
changes sign during actuation, and that there are two points
($\bar{\theta}_{3}\approx 38$ $\degree$C and $\bar{\theta}_{3}\approx 70$
$\degree$C) at which the total heat transfer vanishes at the gate. It is
emphasized that this does not imply that the pointwise normal heat flux at the
gate vanishes.

Fig. 20.: The total heat transfer at the anode, gate and cathode versus the
applied temperature for the Fig. 18 designs. The dashed black line highlights
the target value $\bar{\theta}_{1}=50$ $\degree$C.

##  5 Conclusions

In this work, we have via several numerical examples exemplified how the third
medium contact method can be utilized to model complex thermo-mechanical
contact problems. Indeed, it has been shown that the method inherently
captures the thermal contact resistance that is present in real-world contact
problems via the third medium, i.e. ersatz material properties. It has
specifically been shown that the magnitude of the thermal contact resistance
can be tuned via the ersatz material thermal conductivity contrast,
$\delta_{o}^{\kappa}$. The implemented thermo-mechanical contact framework was
shown to correspond well to an analytical solution.

The fully coupled thermo-mechanical numerical framework was applied to solve
several topology optimization problems involving the design of thermal
regulators that leverage contact as a non-linear mechanism. More specifically,
we designed thermal switches, triodes and a diode. To validate the designs, a
post-processing analysis using body-fitted meshes was conducted.

## Acknowledgments

This work was sponsored by Independent Research Fond Denmark via the TopCon
Project (case number 1032- 00228B).

## References

  * Aalilija et al. [2021] Aalilija, A., Gandin, C.A., Hachem, E., 2021.  A simple and efficient numerical model for thermal contact resistance based on diffuse interface immersed boundary method.  International Journal of Thermal Sciences 166, 106817.
  * Almanza et al. [2018] Almanza, M., Depreux, L., Parrain, F., LoBue, M., 2018\.  Electrostatically actuated thermal switch device for caloric film.  Applied Physics Letters 112.
  * Bendsøe [1989] Bendsøe, M.P., 1989.  Optimal shape design as a material distribution problem.  Structural optimization 1, 193–202.
  * Bluhm et al. [2021] Bluhm, G.L., Sigmund, O., Poulios, K., 2021.  Internal contact modeling for finite strain topology optimization.  Computational Mechanics 67, 1099–1114.
  * Bluhm et al. [2023] Bluhm, G.L., Sigmund, O., Poulios, K., 2023.  Inverse design of mechanical springs with tailored nonlinear elastic response utilizing internal contact.  International Journal of Non-Linear Mechanics 157, 104552.
  * Bosisio et al. [2015] Bosisio, R., Valentini, S., Mazza, F., Benenti, G., Fazio, R., Giovannetti, V., Taddei, F., 2015.  Magnetic thermal switch for heat management at the nanoscale.  Physical Review B 91, 205420\.
  * Castelli et al. [2023] Castelli, L., Zhu, Q., Shimokusu, T.J., Wehmeyer, G., 2023\.  A three-terminal magnetic thermal transistor.  Nature Communications 14, 393\.
  * Cho et al. [2007] Cho, J., Wiser, T., Richards, C., Bahr, D., Richards, R., 2007.  Fabrication and characterization of a thermal switch.  Sensors and Actuators A: Physical 133, 55–63.
  * Chung et al. [2020] Chung, H., Amir, O., Kim, H.A., 2020.  Level-set topology optimization considering nonlinear thermoelasticity.  Computer Methods in Applied Mechanics and Engineering 361, 112735.
  * Dalklint et al. [2023] Dalklint, A., Sjövall, F., Wallin, M., Watts, S., Tortorelli, D., 2023.  Computational design of metamaterials with self contact.  Computer Methods in Applied Mechanics and Engineering 417, 116424.
  * Deetman [2021] Deetman, A., 2021.  GCMMA-MMA-Python: Python Code of the method of moving asymptotes. Retrieved 16/12/2022 from https://github.com/arjendeetman/GCMMA-MMA-Python.
  * Faltus et al. [2024] Faltus, O., Horák, M., Doškář, M., Rokoš, O., 2024.  Third medium finite element contact formulation for pneumatically actuated systems.  arXiv preprint arXiv:2405.01185 .
  * Frederiksen et al. [2023] Frederiksen, A.H., Sigmund, O., Poulios, K., 2023.  Topology optimization of self-contacting structures.  Computational Mechanics .
  * Geuzaine and Remacle [2009] Geuzaine, C., Remacle, J.F., 2009\.  Gmsh: A 3-d finite element mesh generator with built-in pre-and post-processing facilities.  International journal for numerical methods in engineering 79, 1309–1331.
  * Granlund et al. [2024] Granlund, G., Wallin, M., Günther-Hanssen, O., Tortorelli, D., Watts, S., 2024.  Topology optimization of compliant mechanisms under transient thermal conditions.  Computer Methods in Applied Mechanics and Engineering 418, 116478.
  * Guest et al. [2004] Guest, J.K., Prévost, J.H., Belytschko, T., 2004.  Achieving minimum length scale in topology optimization using nodal design variables and projection functions.  International journal for numerical methods in engineering 61, 238–254.
  * Gulfam et al. [2019] Gulfam, R., Zhang, P., Meng, Z., 2019.  Advanced thermal systems driven by paraffin-based phase change materials–a review.  Applied Energy 238, 582–611.
  * Heo et al. [2019] Heo, H., Li, S., Bao, H., Ju, J., 2019.  A passive thermal switch with kirigami-inspired mesostructures.  Advanced Engineering Materials 21, 1900225.
  * Holzapfel [2000] Holzapfel, A.G., 2000.  Nonlinear solid mechanics ii .
  * Jonsmann et al. [1999] Jonsmann, J., Sigmund, O., Bouwstra, S., 1999.  Compliant thermal microactuators.  Sensors and Actuators A: Physical 76, 463–469.
  * Lazarov and Sigmund [2011] Lazarov, B.S., Sigmund, O., 2011\.  Filters in topology optimization based on helmholtz-type differential equations.  International Journal for Numerical Methods in Engineering 86, 765–781.
  * Miao et al. [2022] Miao, R., Kishore, R., Kaur, S., Prasher, R., Dames, C., 2022.  A non-volatile thermal switch for building energy savings.  Cell Reports Physical Science 3.
  * Miehe [1995] Miehe, C., 1995.  Entropic thermoelasticity at finite strains. aspects of the formulation and numerical implementation.  Computer Methods in Applied Mechanics and Engineering 120, 243–269.
  * Milanez and Mantelli [2003] Milanez, F.H., Mantelli, M.B., 2003\.  Theoretical and experimental studies of a bi-metallic heat switch for space applications.  International Journal of Heat and Mass Transfer 46, 4573–4586.
  * Ram et al. [2023] Ram, B.R., Malik, V., Naik, B.K., Patel, K.S., 2023\.  A critical review on mechanical heat switches for engineering and space applications.  Heat Transfer Engineering 44, 1789–1802.
  * Renard and Poulios [2020] Renard, Y., Poulios, K., 2020\.  Getfem: Automated fe modeling of multiphysics problems based on a generic weak form language.  ACM Transactions on Mathematical Software (TOMS) 47, 1–31.
  * Rosochowska et al. [2003] Rosochowska, M., Balendra, R., Chodnikiewicz, K., 2003.  Measurements of thermal contact conductance.  Journal of Materials Processing Technology 135, 204–210.
  * Semiatin et al. [1987] Semiatin, S., Collings, E., Wood, V., Altan, T., 1987\.  Determination of the interface heat transfer coefficient for non-isothermal bulk-forming processes .
  * Sigmund [2001a] Sigmund, O., 2001a.  Design of multiphysics actuators using topology optimization–part i: One-material structures.  Computer methods in applied mechanics and engineering 190, 6577–6604.
  * Sigmund [2001b] Sigmund, O., 2001b.  Design of multiphysics actuators using topology optimization–part ii: Two-material structures.  Computer methods in applied mechanics and engineering 190, 6605–6627.
  * Sui et al. [2023] Sui, Q., Yan, J., Fan, Z., Wallin, M., Ristinmaa, M., Niu, B., 2023.  Topology optimization of thermo-hyperelastic structures utilizing inverse motion based form finding.  Engineering Optimization 55, 110–124.
  * Svanberg [1987] Svanberg, K., 1987.  The method of moving asymptotes—a new method for structural optimization.  International journal for numerical methods in engineering 24, 359–373.
  * Svanberg [2007] Svanberg, K., 2007.  Mma and gcmma-two methods for nonlinear optimization.  vol 1, 1–15.
  * Wang et al. [2011] Wang, F., Lazarov, B.S., Sigmund, O., 2011.  On projection methods, convergence and robust formulations in topology optimization.  Structural and multidisciplinary optimization 43, 767–784.
  * Wehmeyer et al. [2017] Wehmeyer, G., Yabuki, T., Monachon, C., Wu, J., Dames, C., 2017.  Thermal diodes, regulators, and switches: Physical mechanisms and potential applications.  Applied Physics Reviews 4.
  * Wong et al. [2021] Wong, M., Tso, C., Ho, T., Lee, H., 2021.  A review of state of the art thermal diodes and their potential applications.  International Journal of Heat and Mass Transfer 164, 120607.
  * Wriggers et al. [2013] Wriggers, P., Schröder, J., Schwarz, A., 2013.  A finite element method for contact using a third medium.  Computational Mechanics 52, 837–847.
  * Yang et al. [2019] Yang, T., Foulkes, T., Kwon, B., Kang, J.G., Braun, P.V., King, W.P., Miljkovic, N., 2019.  An integrated liquid metal thermal switch for active thermal management of electronics.  IEEE Transactions on Components, Packaging and Manufacturing Technology 9, 2341–2351.