<figure><img src="image_1.png" /></figure>

HAL
open science

# Homogenization based topology optimization of a coupled thermal fluid-structure problem

Godfred Oheneba Agyekum, Laurent Cangémi, François Jouve

► To cite this version:

Godfred Oheneba Agyekum, Laurent Cangémi, François Jouve. Homogenization based topology optimization of a coupled thermal fluid-structure problem. Computers & Structures, 2025, 308, pp.107636. 10.1016/j.compstruc.2024.107636. hal-04919691

HAL Id: hal-04919691

https://ifp.hal.science/hal-04919691v1

Submitted on 29 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Homogenization based topology optimization of a coupled thermal fluid-structure problem

Godfred Oheneba Agyekum<sup>a, ID,*</sup>, Laurent Cangémi<sup>a, ID,†</sup>, François Jouve<sup>b,†</sup>

## ARTICLE INFO

### Keywords:

Topology optimization

Multi-scale

Periodic homogenization

Porous medium

Adjoint methods

Fluid-structure interaction

Convective heat-transfer

## ABSTRACT

This article focuses on the topology optimization of a weakly coupled three physics problem. The structures are made of periodically perforated material, where the microscopic periodic cell is macroscopically modulated. The objective is to optimize the homogenized formulation of this system, where the coupling is weak because the three physics involved are solved consecutively: first, a coupled fluid flow is determined using the Biot-Darcy's law for the fluid domain, second, a thermal model using the convection-diffusion equation for the whole domain, and third, a three-physics problem by solving the linear poro-thermo elasticity problem in the solid domain. This approach allows low computational cost of evaluation of load sensitivities using the adjoint-state method. Two-dimensional and three-dimensional numerical problems are presented using the alternate directions algorithm. It is demonstrated how the implementation makes it possible to treat a variety of design problems.

## 1. Introduction

Shape optimization that involve pressure-loaded or thermal fluid-loaded boundaries is a very classical problem, and it is acknowledged that keeping a clear definition of the boundary is a key ingredient for efficiency in many cases. Indeed, using boundary identification schemes, Hammer and Olhoff [1] (2000) proposed an iso-density approach to identify the pressure loading boundaries. However, the proposed technique presented separate isolines-islands, and to resolve this issue, Du and Olhoff [2] (2004) introduced a modified isolines technique that generates the appropriate loading surface. In analogue, using the density-based technique, Fuchs and Shemesh [3] (2004) introduced an additional set of variables to treat the evolving pressure loading boundary. Subsequent research in [4] (2007), extended this to develop a topology optimization of regions of Darcy and Stokes flow, using a regularization method. Zheng et al. [5] (2009) introduced a pseudo equal-potential function to treat the evolving position and direction of pressure loading. In analogue, using the SIMP model, Lee and Martins [6] (2012) introduced a material boundary identification scheme by iteratively connecting points of equal density, with the use of a predefined void region. Subsequent research in [7] (2016), presented a new scheme based on image segmentation technique to recover the loading surface: the distance

regularized level set evolution [8] (2010) (DRLSE) method is utilized to detect the image edge. Similarly, Li et al. [9] proposed a new approach based on digital image processing and regional contour tracking technology to generate the appropriate loading surface.

For thermal fluid-loaded problems, many studies have been conducted with using an explicit boundary between the fluid and structure. Yaji et al. [10] (2015) presented a topology optimization method for a coupled thermal-fluid problem, using a level set boundary expressions and a Tikhonov-based regularization scheme. In analogue, Feppon et al. [11] (2019) proposed a level set mesh evolution framework: the Hadamard's method of shape differentiation is used to maintain distinct fluid and solid domains. Similarly, Picelli et al. [12] (2019) introduced a level-set strategy for design-dependent hydrostatic pressure loading: the Ersatz material approach with fixed grid is applied together with work equivalent load integration.

Subsequent research extended this to develop a topology optimization method using homogenization based approach that avoid detecting the loading surface. Andreasen and Sigmund [13,14] (2011 and 2013) proposed a multiple-scale method applied to periodic media: the resulting model couples the Stokes flow in the pores of the structure with the deformation of the elastic skeleton through a macroscopic Darcy-type flow law. In analogue, using the SIMP method, Chen et al. [15] (2001)

<footnote><sup>a</sup> IFP Energies nouvelles, Rueil-Malmaison, 92852, France</footnote>

<footnote><sup>b</sup> Université Paris Cité, Laboratoire Jacques-Louis Lions (LJLL), Paris, F-75006, France</footnote>

<footnote>* Corresponding author.
E-mail addresses: godfred.a.o.ezail@gmail.com (G. Oheneba Agyekum), laurent.cangemi@ifpen.fr (L. Cangémi), francois.jouve@u-paris.fr (F. Jouve).</footnote>

<footnote>† These authors contributed equally to this work.</footnote>
proposed a new approach based on applying a fictitious thermal load-
ing to simulate a three-phase material distribution problem within the
design domain in which the solid, void, and hydrostatic fluid phases are
optimally distributed. Similarly, Yoon [16] (2012) developed a new and
rigorous monolithic analysis and optimization framework for an electro-
fluid-thermal-compliant (EFTC) microactuator. Furthermore, Marck et
al. [17] (2013) proposed a finite-volume method (FVM), for the direct
solver, with the discrete adjoint approach, to tackle both fluid dynamic
and heat transfer optimization in the frame of laminar flows. Similarly,
Marck and Privat [18] (2014), presented a new homogenization strat-
egy to two configurations subject to conductive and convective heat
transfers. Subsequent research based on a review about the engineer-
ing design of optimal heat transfer systems, using topology optimiza-
tion techniques is presented in [19] (2017). Furthermore, Mendes et
al. [20] (2022) introduced a topology optimization of binary structures
(TOBS) method to solve problems that consider buckling constraints
and design-dependent loads. By virtue of the homogenization, Hübner
et al. [21] (2019) introduced a shape optimization of microstructures
generating porous locally periodic materials saturated by viscous flu-
ids: the porous material is described by the Biot model, and by the
hydraulic permeability of the Darcy flow model. Subsequent research
extended this to develop a new density-based topology optimization ap-
proach using Darcy's law with a drainage term to provide a treatment of
design-dependent fluidic pressure loads [22] (2020). In analogue, Ku-
mar [23] (2023) presented a compact 100-line MATLAB code, TOPress,
developed according to the approach first reported in Kumar et al. [22]
(2020). Our approach follows a similar strategy developed in [21,22],
which has not been reported before.

As seen above, several optimization methods exist to optimize the
topology, and therefore the holes in a structure. The SIMP method and
the level-set method are based on finite elements and so require to
mesh the working domain. Consequently, if we consider lattice mate-
rials (i.e., microscopic holes in the structures), the mesh should be fine
enough to capture their shapes, and in practice, such approaches are
prohibitive to compute numerically. A solution is given by the homoge-
nization methods [24–28]: allow for fine mixtures of void and material
on a scale which is much smaller than the mesh used for the working
domain. However, the set of effective tensors resulting from the mixture
in fixed volume fraction of two elastic materials is unknown, albeit its
conductivity analogue is known [27]. For a general summary of the ho-
mogenization method, we refer the reader to [29–38], and references
therein. In this context, we adopt an approach developed in [29], that
was first numerically implemented and popularized. Thereby, we define
a new equivalent full media, with a behavior similar to that of lattice
materials. Then a relatively coarse mesh can be used to run the numer-
ical optimization algorithm.

This paper proposes an homogenization based topology optimiza-
tion, where the effective medium properties are computed using char-
acteristic responses of the representative unit cell, namely a perforated
hexagonal cell in 2-D and a tetrakaidecahedron cell in 3-D, constituted
by void and solid phases. We propose an extension of the Biot-Darcy's
model introduced in our preliminary work [39] (see also our work
in [40]), in two-dimensional. In this regard, we propose to coupled
the Biot-Darcy's law to a convective heat propagation, and to optimize
the homogenized formulation of a weakly coupled thermal-fluid sys-
tems, fluid-loaded structures and thermal fluid-loaded structures, where
it is possible to consider a weak coupling between the multi-physics at
stake because one can assume that the fluid domain is fixed at first or-
der; howbeit, this weak coupling is a major simplification and therefore
reduces the computational cost. The sensitivities of general “smooth”
objective functions with respect to design variables, are provided us-
ing the adjoint-state method (see [41–43] and the references therein for
details), and the homogenized coefficients are differentiated using the
Lagrangian method and the design (or optimization) variable is opti-
mized using the projected gradient algorithm [32].

The remainder of this paper is organized as follows: in Section 2, the
necessary ingredients of the homogenization method are briefly recalled
and we describe our methodology. First, we choose a representative unit
cell. Second, we compute (offline and once and for all) its homoge-
nized properties for the entire range of densities (that will be stored
once and for all to be used in the later optimizations). In Section 3,
our weakly coupled model of heat propagation, fluid flow and structure
strain is presented. In Section 4, the optimization problem formulation
and the associated sensitivity analysis are introduced, which turns out
to be a simple parametric optimization problem since our representa-
tive unit cell is parametrized. Section 5, is concerned with our topology
optimization algorithm: an alternate directions algorithm [32], which
successively computes the stress field through the solving of a weakly
coupled three physics problem. Finally in Section 6, we present our
numerical results: 2-D and 3-D computations are displayed of various
design problems involving two or three-physic for some objective func-
tions. We summarize our findings and give an outlook in Section 7.

## 2. Periodic homogenization

We briefly present the principles of the homogenization method, and
explain our methodology; the mathematical framework of this theory
can be found in [35], and for a general summary of the homogeniza-
tion method, we refer the reader to [24–35] and references therein.
We restrict our analysis to locally periodic hexagonal cell in 2-D and
tetrakaidecahedron cell in 3-D. The objective is to determine the ho-
mogenized physical properties of such materials when varying their
parameter; this approach is independent of the objective function, com-
putational domain, applied loads or boundary conditions.

### 2.1. Set of admissible microstructures

From now on, we restrain our analysis to a simple class of isotropic
microstructures: a smooth honeycomb cell in 2-D (Fig. 1(b)) and
tetrakaidecahedron cell in 3-D (Fig. 3). This class of isotropic cells is
of great interest for shape optimization: first, there is no need to ori-
entate them locally in the domain (which simplifies the optimization
process but also the deshomogenization one), and second, the geome-
try is parametrized by only one parameter (the local density); therefore,
only one parameter has to be optimized, which makes the optimization
process very close to the SIMP method.

#### 2.1.1. Smooth honeycomb cell

Our interest for the smooth honeycomb relies on its smooth rounded
corners, known to generate lower local concentration stress [24,45],
compared to the classical honeycomb (see Fig. 1).

Let $Y$ be a unit regular hexagonal in 2-D, similar to the so called
classical honeycomb (see Fig. 1(a)), except that the interior corners of
the central hexagonal hole are rounded. As a consequence, the cen-
tral smooth hexagonal can not reach the outer regular hexagonal: i.e.,
the density $\theta$ going to zero is excluded. Thus, for practical reason, an-
other parameter $h \in [0, 1]$ (homogeneous to a distance) is introduced to
parametrize the smooth honeycomb; its describes (similar to the param-
eter $e$) the relative pore diameter of the cell, defined by

$$h = 1 - \frac{2e}{\sqrt{3}}, \qquad (1)$$

where the relative width $e$ is defined by

$$e = \frac{1}{2\sqrt{3}}(1 - \sqrt{1 - \theta}), \qquad (2)$$

with $e$ ranging from 0 to $\frac{1}{2\sqrt{3}}$. In this context, a parametric curve $\Gamma_h$ is introduced, and represents the inner smooth central hole: note that, this only applies in 2-D. For its geometric definition, we refer the reader to our preliminary work [39], Section 2.1.
<figure><img src="image_4.png" /><figcaption>Fig. 1. 2D isotropic design cells: (a) Classical honeycomb compared to (b) Smooth honeycomb (images taken from [44]).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 2. Tetrakaidecahedron: an isotropic design cell in 3-D (image taken from [44]).</figcaption></figure>

### 2.1.2. Tetrakaidecahedron cell

In 3-D, isotropic cells are not so easy to design; here, we design a regular tetrakaidecahedron, also known as the Kelvin foam, similar to the one introduced in [37,44]. It is shown to yield isotropic (or quasi isotropic) homogenized Hooke's laws. It features fourteen faces: six unit squares and eight regular unit hexagons; see Fig. 2 and Fig. 3.

Several approaches are possible in order to design isotropic cells from this one, parametrized by one parameter (i.e., the local density, which is defined by the ratio of the homothety between the hole in a face and the external border of this face). A naive approach consists in adding a central tetrakaidecahedron inclusion characterized by its size in [0, 1].

Based on the construction method described above, the material density $\theta$ has been used as a control parameter to generate foams of chosen porosities. Upon increasing the local density $\theta$, one would reach a point where the structure is full of material: the central tetrakaidecahedron inclusion will be closed. For more details, we refer the reader to Appendix A2, Sect. A.1.

### 2.1.3. Deshomogenization process

Concerning the mapping process between the actual microstructure geometry and the design variable densities, we refer the reader to Appendix A2, Sect. A.2. Fortunately, when the orientation of the lattice cells is uniform, this is not a big challenge: the periodic size of the final microstructure is chosen, then on each cell we use the average parameters to design the local microstructure. However, the connectivity between adjacent cells has to be carefully ensured.

### 2.2. Cell problem and homogenized elasticity tensor

We briefly recall some important results on the theory of homogenization; the interested reader can refer to [35], for more details.

Let Y be the unit periodic hexagonal in 2-D or the tetrakaidecahedron in 3-D, with its volume denoted $|Y|$. Then, the homogenized (or effective) tensor $A^*$ is defined by first solving the variational formulation, that is: for each pair $(i, j) \in \{1, ..., N\} (N = 2, 3)$, find $w_{ij} \in H^1_\#(Y, \mathbb{R}^N) = \{w \in H^1(Y, \mathbb{R}^N) | w \text{ is Y-periodic}\}$ such that

$$ \forall \phi \in H^1_\#(Y, \mathbb{R}^N), \quad \int_Y Ae(w_{ij}) : e(\phi) \, dy + \int_Y Ae_{ij} : e(\phi) \, dy = 0, \quad (3) $$

where $w_{ij}$ is the so-called correctors (corresponding to the local displacements in Y), and $e_{ij} = \frac{1}{2}(e_i \otimes e_j + e_j \otimes e_i)$ is the basis of the symmetric tensors of order 2. The entries of the effective tensor $A^*$ are then computed in terms of the correctors $w_{ij}$:

$$ A_{ijkl}^* = \frac{1}{|Y|} \int_Y A(e_{ij} + e(w_{ij})) : (e_{kl} + e(w_{kl})) \, dy \quad \forall i, j, k, l \in \{1, \dots, N\} \quad (4) $$

Note that, Y is termed periodic because the periodicity of the cell is defined by the same displacement on two opposite and parallel faces; thus, there are three directions of periodicity in 2-D (for the hexagonal cell) and seven directions of periodicity in 3-D (for the tetrakaidecahedron cell).

### 2.2.1. Numerical results

Here, we used for all finite element operations a typical mesh consisting of triangular elements (in 2-D) and tetrahedral elements (in 3-D); the range of the density is regularly discretized with 60 elements, in 2D and with 20 elements in 3D. For both periodic cells, we have computed the homogenized tensor $A^*(\theta)$ for discrete values $(\theta_i)_{i=1,...,n_s}$ of $\theta$, with constant elastic tensor $A$ of Young's modulus $E = 12 \times 10^9 \text{ N m}^{-2}$ (12 GPa) and Poisson's ratio $\nu = 0.35$; a preprocessing stage, which is performed offline: a table of surrogate model of the size $n_s = 1000$ is constructed, and it is used to compute (via a linear interpolation approach) the local composites during the optimization process. The isotropy of $A^*$ is satisfied if we have the following relations:

$$ \forall i, j, k, l \in \{1, \dots, N\} \quad \left\{ \begin{array}{l} A_{ijkl}^* = A_{klij}^* \\ A_{ijkl}^* = 0; \text{ if } j \neq k \\ A_{ijjj}^* = A_{kkll}^*; \text{ if } i \neq j, k \neq l \\ A_{iiii}^* = A_{ijij}^* + A_{iijj}^*; \text{ if } i \neq j \end{array} \right. \quad (5) $$

The relative errors with respect to the equalities in Eq. (5), for the Kelvin foam (i.e., the tetrakaidecahedron cell in 3-D) are depicted by Fig. 4: we emphasize that the residuals are sufficiently small, namely the relative errors are less than $10^{-4}$, which validates the isotropy assumption. Note that, the above results still hold for the smooth honeycomb cell in 2-D: see our preliminary work [39], in section 2.2 for justification.

### 2.3. Sensitivity of the homogenized elasticity tensor

The computation of the sensitivity of the homogenized elasticity tensor with respect to the design parameter of the periodic cell $Y(\theta)$ is based on the notion of shape derivative. We define $W^1_\#(Y(\theta); \mathbb{R}^N)$ as the set of $Y$-periodic Lipschitz maps from $Y(\theta)$ with values in $\mathbb{R}^N$, for a given $\theta \in [0, 1]$.

**Lemma 1.** Let $\theta \in W^1_\#(Y; \mathbb{R}^N)$. The shape derivative of $A_{ijkl}^*$ is given by:

$$ \langle \langle (A_{ijkl}^*)' (Y), \theta \rangle \rangle = \int_Y A(e_{ij} + e(w_{ij})) : (e_{kl} + e(w_{kl})) \, div(\theta) \, dy \\ - \int_Y A(e_{ij} + e(w_{ij})) : \langle de(w_{kl}), \theta \rangle \, dy \quad (6) \\ - \int_Y A(e_{kl} + e(w_{kl})) : \langle de(w_{ij}), \theta \rangle \, dy, $$

where $de(w)$ is a linear operator from $W^1_\#(Y; \mathbb{R}^N)$ to $L^2(Y; M_N^s)$ defined for every $w \in H^1_\#(Y; \mathbb{R}^N)$ by
<figure><img src="image_6.png" /><figcaption>Fig. 3. Tetrakaidecahedron (or Kelvin foam) for different values of the material density $\theta$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 4. The relative errors with respect to the equalities in (5), for the Kelvin foam cell.</figcaption></figure>

$$ \langle de(w), \vartheta \rangle = \frac{1}{2}(\nabla w \nabla \vartheta + \nabla \vartheta^T \nabla w^T), $$

and the function $A_{ijkl}^*(Y)$ is defined as the Fréchet derivative in $W^{1,\infty}$ at 0 of the application $\vartheta \to A_{ijkl}^*((Id + \vartheta)Y)$, with $A_{ijkl}^* \in M_N^s$ in the set of squared $N \times N$ and symmetric ones.

See [35,46], for a proof of Lemma 1: it relies on the Lagrangian method, which amounts to introduce a Lagrangian, defined as the sum of formula (4) for $A_{ijkl}^*$ and of the variational formulation (3). Differentiating with respect to the state variable gives the adjoint system. It turns out that the problem is self-adjoint, so no adjoint appears in Eq. (6). Differentiating with respect to the shape leads to the final result.

## 2.4. Derivatives of the homogenized elasticity tensor

In order to compute the derivative of the homogenized tensor, with respect to the design parameter $\theta$ (the local density), we choose a specific deformation field $\vartheta$ in the shape derivative formula (6). Introducing the smooth $(0, 1)^2$-periodic field $\vartheta$ defined by:

$$ \vartheta = c \sin(2\pi y), \qquad (7) $$

where $c = -\sin(\pi\theta)^{-1}$ is the rescaling coefficient, and such that:

$$ Y(\theta + \delta\theta) = (Id + \delta\theta\vartheta)Y(\theta) $$

where $\delta\theta$ is a small increment. It follows that:

$$ \frac{\partial A^*}{\partial \theta}(\theta) = \langle (A^*)'(Y), \vartheta \rangle $$

Therefore, the shape derivative formula (6) leads to the sensitivities of $A^*(\theta)$ with respect to the design parameter $\theta$.

## 2.5. Cell problem and effective model in conductivity

Similar to the elasticity case in Section 2.2, we consider a general model problem of thermal or pressure conductivity, with piecewise constant conductivity $K$. The homogenized conductivity $K^*$ is defined by first solving the variational formulation, that is: for each $i \in \{1, ..., N\}$ ($N = 2, 3$), find $w_i \in H_\#^1(Y, \mathbb{R}) = \{w \in H^1(Y, \mathbb{R}) | w \text{ is Y-periodic}\}$ such that

$$ \forall \phi \in H_\#^1(Y, \mathbb{R}), \quad \int_Y K \nabla w_i : \nabla \phi \, dy + \int_Y K e_i : \nabla \phi \, dy = 0, \qquad (8) $$

where $w_i$ is the so-called cell problems (corresponding to the local variation of pressure created by an averaged gradient), and $(e_i)_{i=1,...,N}$ is the canonical basis of $\mathbb{R}^N$. The entries of the homogenized conductivity $K^*$ are then computed in terms of the correctors $w_i$:

$$ (K^*)_i_j = \frac{1}{|Y|} \int_Y K(e_i + \nabla w_i) \cdot (e_j + \nabla w_j) \, dy \quad \forall i, j \in \{1, 2\} \qquad (9) $$

### 2.5.1. Numerical results

For both periodic cells, the homogenized conductivity $K^*(\theta)$ has been computed for discrete values $(\theta_i)_{i=1,...,n_s}$ of $\theta$, with piecewise constant $K$ defined by

$$ K(y) = \begin{cases} k_s & \text{in } y \in Y_0 \\ k_f & \text{in } y \in Y \setminus Y_0 \end{cases} \qquad (10) $$

where $k_s$ and $k_f$ are some fixed thermal or pressure conductivities inside the solid and fluid phases: here, the chosen solid to fluid thermal conductivity ratio is $k_s/k_f = 10$ (with $k_f = 10$ W/mK, similar to work in [47], for isotropic metal foams), and the solid to fluid pressure per-
<figure><img src="image_8.png" /><figcaption>Fig. 5. The homogenized conductivity $\frac{\alpha^*}{||\alpha^*||_\infty}$ as a function of $\theta$ for the Kelvin foam, compared to $K_D^*$.</figcaption></figure>

meability ratio is $k_s/k_f = 0.001$ (with $k_f = 1 \text{ m}^2$). We emphasize that, if $Y$ is isotropic (e.g. the honeycomb in 2-D or Kelvin foam in 3-D), then we get: $\forall i, j \in \{1, ..., N\}, K_{ii}^* = K_{jj}^*$ and $K_{ij}^* = 0$, if $i \neq j$. Thus, $K^*$ is proportional to the identity tensor $I_N$, and as such, it can be identified to the scalar $\alpha^* = K_{11}^*$.

On Fig. 5, we display for the Kelvin foam the effective conductivity $\frac{\alpha^*}{||\alpha^*||_\infty}$ (with $||\alpha^*||_\infty = \max_{\theta_i}|\alpha^*(\theta_i)|$) for discrete values $(\theta_i)_{i=1,...,n_s}$, compared to a given smooth function:

$$K_D^* = \min \left( \frac{\epsilon_0 + (1 - \epsilon_0)(1 - \theta)}{\theta}, K_\infty \right), \quad (11)$$

where $\epsilon_0$ and $K_\infty$ are some fixed constants: i.e., $\epsilon_0 = 10^{-3}$, $K_\infty = 10^3$. We note that, the homogenized conductivity is a decreasing function with respect to the density $\theta$, and the residual errors $|K_{11}^* - K_{22}^*|/K_{11}^* \le 10^{-6}$ and $K_{12}^*/K_{11}^* \le 10^{-3}$ are small enough. Note that, the above results still hold for the smooth honeycomb cell in 2-D: see the work in [39], section 2.3 for justification.

### 3. Setting of the three-physic problem

Throughout this paper, we assume that $\Omega$ is a fixed domain in $\mathbb{R}^N$ ($N = 2$, or 3), filled with composite material periodically perforated by the hexagonal cell in 2-D or tetrakaidecahedron cell in 3-D, and characterized by $\theta$, with the normal vector to the boundary $\partial\Omega$ denoted $n$, pointing outward to $\Omega$.

We now describe in a formal way in Section 3.1, Section 3.2 and Section 3.3, the coupled three models, with their relevant set of boundary conditions.

#### 3.1. Hydraulic model

The motion of the fluid inside the porous medium $\Omega$ is based on the Darcy's law [48], characterized by the pressure field $p$, satisfying a Biot's law [39] and a law of conservation of mass (in view of incompressible fluid):

$$\left\{
\begin{array}{ll}
-div(K_D^* \nabla p) = 0 & \text{in } \Omega, \\
p = p_0 & \text{on } \Gamma_D^f, \\
q_{\Gamma} \cdot n = f_f & \text{on } \Gamma_N^f, \\
q_{\Gamma} \cdot n = 0 & \text{on } \Gamma^f = \partial \Omega \setminus (\Gamma_D^f \cup \Gamma_N^f)
\end{array}
\right.
\quad (12)$$

where $f_f$ is a given mass flow rate (satisfying $-K_D^* \frac{\partial p}{\partial n} = f_f$, a Neumann isoflux condition). The boundary of the fluid phase is the disjoint re-union

$$\partial\Omega = \Gamma_D^f \cup \Gamma_N^f \cup \Gamma^f$$

of a Dirichlet (or inlet) part $\Gamma_D^f$ where the fluid enters with a given pressure $p = p_0$, a Neumann (or outlet) part $\Gamma_N^f$ where normal stress (i.e., mass flow rate) $q_{\Gamma} \cdot n = f_f$ is imposed, and free interface $\Gamma^f$. Note that, a more accurate description of the hydraulic model (12) and the physical equation chosen for the modeling of the state variable $p$ can be found in [39], Section 3.1. It is worth mentioning that in 2-D, the hexagonal cell geometry cannot yield the Darcy equation (12) because the fluid is not connected. Therefore model (12) is just a crude approximation and not a true homogenized model.

#### 3.2. Convection-diffusion for the temperature variable

The coupled diffusion of heat inside the whole domain $\Omega$ and its transport in the fluid domain is based on the convection-diffusion equation: the fluid velocity $v_f$ determines the physical behavior of the temperature $T$ in the whole domain, as a result of convection and diffusion effects inside the fluid domain and of pure diffusion inside the solid domain (see e.g. [47]). In this work, we have chosen to use an empirical form of the heat equation for the whole domain featuring fluid and solid phases (see e.g. [49]). Some authors employ two-equation modeling for each phase of the medium (see e.g. [50]). This modeling difference lies in the fact that the empirical form at one temperature does not take into account certain local turbulence effects, which degrades the heat transfer across the interface between the fluid and solid phases. However, even in this case, the thermal capacity is indeed that of the fluid alone.

We now assume that, when $\epsilon$ tends to 0, the periodic microstructure $\Omega_\epsilon$ tends to an homogeneous domain $\Omega$, filled with fine mixtures of solid and fluid phases; let $K^*$ denote the homogenized thermal conductivity inside the porous media (which we numerically computed; see Section 2.5.1), and let $c_p$ be the thermal capacity of the fluid. Then, the weakly coupled diffusion of the heat model is defined by the convection-diffusion equation:

$$\left\{
\begin{array}{ll}
-div(K^* \nabla T) + \phi(\rho c_p)_{f_f} v_f \nabla T = 0 & \text{in } \Omega, \\
T = T_0 & \text{on } \Gamma_D^T, \\
-K^* \frac{\partial T}{\partial n} = h_0 & \text{on } \Gamma_N^T, \\
-K^* \frac{\partial T}{\partial n} = 0 & \text{on } \Gamma^T = \partial \Omega \setminus (\Gamma_D^T \cup \Gamma_N^T)
\end{array}
\right.
\quad (13)$$

where the subscript $f$ denotes the restriction to the fluid phase in $\Omega$. The boundary $\partial\Omega = \Gamma_D^T \cup \Gamma_N^T \cup \Gamma^T$ is split into a Dirichlet part $\Gamma_D^T$, where a temperature $T_0$ is imposed, a Neumann part $\Gamma_N^T$, where a given incoming heat flow $h_0$ is applied and a free interface $\Gamma^T$. Note that, $\phi = 1 - \theta$ defines the volume fraction of the fluid phase: i.e., the proportion of the void phase at each point $x \in \Omega$. The convection-diffusion equation (13) is now recast as thermal-Biot-Darcy equation:

$$\left\{
\begin{array}{ll}
-div(K^* \nabla T) + (-c_p K_D^* \nabla p)_{f_f} \nabla T = 0 & \text{in } \Omega, \\
T = T_0 & \text{on } \Gamma_D^T, \\
-K^* \frac{\partial T}{\partial n} = h_0 & \text{on } \Gamma_N^T, \\
-K^* \frac{\partial T}{\partial n} = 0 & \text{on } \Gamma^T
\end{array}
\right.
\quad (14)$$

where $\phi \rho v_f = m v_f (= -K_D^* \nabla p)$ is the fluid flow, satisfying the Biot-Darcy model (12). Similar to the Biot-Darcy model (12), starting from a microscopic description of a problem, one seeks a macroscopic (or effective) model problem in conductivity; we consider a model problem of thermal flow in a periodic medium: an heterogeneous domain obtained by mixing periodically two different phases, one being the solid phase and the other being the fluid inclusions. As seen in Section 2, to
compute the effective conductivity $K^*$, we introduce the cell problems, and since the considered cell $Y$ is specifically chosen in order to design isotropic composites, only one of its coefficient (e.g. $(K^*)_{11}$) could be computed in order to fully characterize $K^*$.

### 3.3. Fluid-thermal structure interaction

The solid deformation in $\Omega$, is determined as a result of the stress exerted by the fluid part, the dilation induced by thermoelastic effects and of external boundary forces. We now consider that $\Omega$ is filled with an isotropic composite material, with homogenized Lamé coefficients denoted $\lambda^*, \mu^*$; let $\alpha$ be the linear thermal expansion parameter (in an isotropic solid) and $T_{ref}$ be the temperature at rest. Then, the weakly coupled fluid-thermal elastic model (of Biot-Coussy type) is defined by

$$ \begin{cases} -div(\sigma(u, T)) = -b \nabla p & \text{in } \Omega, \\ u = u_0 & \text{on } \Gamma_D^s, \\ \sigma(u, T) \cdot \mathbf{n} = f_s & \text{on } \Gamma_N^s, \\ \sigma(u, T) \cdot \mathbf{n} = 0 & \text{on } \Gamma_D^s, \end{cases} \quad (15) $$

where

$$ \sigma(u, T) = A^*(e(u) - \alpha(T_s - T_{ref})I_N) $$

represents the effective thermoelastic tensor, and with the effective elastic tensor defined by

$$ A^* e(u) = 2\mu^* e(u) + \lambda^* \mathrm{Tr}(e(u)) I_N, $$

where $I_N$ is the identity tensor. In Eq. (15), the boundary $\partial\Omega$ is split into respectively a Dirichlet part $\Gamma_D^s$ (where a displacement $u = u_0$ is prescribed), a Neumann part $\Gamma_N^s$ (where a surface force $f_s$ is imposed), and a free part $\Gamma^s$; the coefficient $b$ appearing in the source term $-b\nabla p$ (equivalent to the body force exerted by the fluid part) is the so-called Biot coefficient, which in case of linear isotropic elastic tensor (microscopically homogeneous) is explicitly defined by

$$ b(\theta) = 1 - \frac{\kappa_s(\theta)}{\kappa}, \quad (16) $$

where $\kappa$ and $\kappa_s(\theta)$ represent the bulk moduli of the elastic tensor $A$ and the homogenized tensor $A^*$. In this context, note that the above system (15) is a simplified version of a genuine three-physics model: we refer the reader to Appendix B.1 for more details on the complex version, and its domain of validity. Thanks to this simplification, the systems (12), (14) and (15) are only weakly coupled: the resolution is achieved by solving the fluid system (12), then using the fluid stress resulting from the pressure $p$ in the heat transfer equation (14), and finally using the fluid stress and the temperature $T$ to solve the linear poro-thermo elasticity system (15).

## 4. The optimization problem formulation and its sensitivity analysis

The optimization problem formulation for the weakly coupled three-physics system (15) is introduced and the associated sensitivity analysis is discussed. Our objective is the resolution of the homogenized constrained version of the original optimization problem:

$$ \begin{array}{l} \min_{\theta} J^*(\theta, p(\theta), T(\theta), u(\theta)) \\ \text{s.t} \quad \begin{cases} g_i(\theta, p(\theta), T(\theta), u(\theta)) = 0, & 1 \le i \le p, \\ h_j(\theta, p(\theta), T(\theta), u(\theta)) \le 0, & 1 \le j \le q, \end{cases} \end{array} \quad (17) $$

where $J^*$ is a given objective (or cost) function, namely, a trade-off between the compliance of solid, the hydraulic power dissipated by the fluid and the heat convected by the fluid:

$$ \begin{aligned} J^*(\theta, p, T, u) ={}& \omega_1 \underbrace{\left( \int_{\Omega} A^* e(u) : e(u) \, dx \right)}_{\text{Elastic strain energy}} \\ & + \omega_2 \underbrace{\left( - \int_{\Omega} K_D^* \nabla p \cdot \nabla p \, dx \right)}_{\text{Hydraulic power}} \\ & + \omega_3 \left( \int_{\Omega} c_p K_D^* \nabla p \cdot \nabla T \, dx \right), \end{aligned} \quad (18) $$

upon which we shall impose adequate equality and inequality constraints $g_i$ and $h_j$: the physical load specifications needed to reach (e.g. a desired volume or mass for the structure) or that should not be exceeded (e.g. an upper bound limit for the overall temperature). Here, $\theta$ denotes the material density of the structure to be optimized, and $\omega_j \in [0, 1]$ (for $j \in \{1, 2, 3\}$) are some given coefficients, such that $\omega_1 + \omega_2 + \omega_3 = 1$. Usually, $J^*$ refers to a given objective function which quantifies the system performance: classical objective functions commonly encountered include the compliance of a mechanical structure, the drag force induced by an airfoil, or the heat stored into a thermal system. Both objective function $J^*$ and constraints $g_i, h_j$ depend on the optimization variable $\theta$ and on the physical variables $p, T, u$, which depend themselves on $\theta$ through physical state equations.

For its resolution, we shall rely on the alternate directions algorithm [32] detailed in Section 5.1, which requires to compute the sensitivity of the above objective function: the adjoint-state method [35] is used to determine the sensitivities of the objective function and constraints with respect to the design variable $\theta$. The “shape” is characterized by the material density $\theta$, and the shape derivatives are performed using the Lagrangian method and the design is optimized using the projected gradient algorithm. In this framework, this means computing the derivative of the mapping:

$$ \theta \mapsto J^*(\theta, p(\theta), T(\theta), u(\theta)), $$

where $\theta$ belongs in general to the set of admissible design variables $U_{ad}$, defined as follows:

$$ U_{ad} := \left\{ \theta \in L^{\infty}(\Omega; \mathbb{R}_+) \mid \theta(x) \in [0, 1], \forall x \in \Omega \right\} $$

It worth mentioning that, in a recent work in [21], the optimization variable was related to its effective medium parameters.

### 4.1. A fully Lagrangian setting for the sensitivity analysis

We emphasize that the homogenized functional $J^*$, appearing naturally in (17) is not so convenient for the mathematical analysis. The domain of definition of the objective $J^*(\theta, \dots, \theta)$ is a functional space, which depends on the material density $\theta$. In order to discuss the precise mathematical settings, we introduce the following functional spaces

$$ \begin{aligned} V(\Gamma_D^f) &= \{q \in H^1(\Omega) \mid q = 0 \text{ on } \Gamma_D^f\}, & \text{for the pressure } p \\ V(\Gamma_D^T) &= \{S \in H^1(\Omega) \mid S = 0 \text{ on } \Gamma_D^T\}, & \text{for the temperature } T \\ V(\Gamma_D^s) &= \{v \in H^1(\Omega)^N \mid v = 0 \text{ on } \Gamma_D^s\}, & \text{for the displacement } u \end{aligned} \quad (19) $$

We emphasize that, the non-homogeneous Dirichlet boundary data $u_0, p_0$ and $T_0$ featured in Eq. (12) to (15) are seen as the traces on the boundary of $H^1(\Omega)$ functions; the state variables $p, T$ and $u$ are then solutions to the following variational problems: find $(p, T, u) \in \{p_0 + V(\Gamma_D^f), T_0 + V(\Gamma_D^T), u_0 + V(\Gamma_D^s)\}$ such that,
$$
\left\{
\begin{array}{l}
\displaystyle \int_{\Omega} K_D^* \nabla p \cdot \nabla q \, dx + \int_{\Gamma_N^f} f_f q \, ds = 0, \quad \forall q \in V(\Gamma_D^f), \\[2ex]
\displaystyle \int_{\Omega} K^* \nabla T \cdot \nabla S \, dx + \int_{\Gamma_N^T} h_0 S \, ds + \int_{\Omega} (-c_p S K_D^* \nabla p) \cdot \nabla T \, dx = 0, \\[2ex]
\forall S \in V(\Gamma_D^T), \\[2ex]
\displaystyle \int_{\Omega} \sigma(u, T_s) : e(v) \, dx - \int_{\Gamma_N^s} f_s \cdot v \, ds + \int_{\Omega} b \nabla p \cdot v \, dx = 0, \\[2ex]
\forall v \in V(\Gamma_D^s),
\end{array}
\right.
\tag{20}
$$

We now address the sensitivity analysis of $J^*$, which is performed within a Lagrangian framework. The classical idea is to define the corresponding Lagrangian to the optimization problem (17), which is defined by:

$$
\begin{align}
\mathcal{L}(\hat{\theta}, \hat{u}, \hat{u}, \hat{p}, \underline{\hat{p}}, \underline{\hat{T}}, \underline{\hat{T}}, \ell) &= J^*(\hat{\theta}) + \int_{\Omega} (\sigma(\hat{u}, \hat{T}) : e(\underline{\hat{u}}) + b(\hat{\theta})\nabla\hat{p} \cdot \underline{\hat{u}}) \, dx \nonumber \\
&\quad + \int_{\Omega} K^* \nabla \hat{T} \cdot \nabla \underline{\hat{T}} \, dx - \int_{\Omega} c_p \underline{\hat{T}} K_D^* \nabla \hat{p} \cdot \nabla \hat{T} \, dx \tag{21} \\
&\quad + \int_{\Omega} K_D^* \nabla \hat{p} \cdot \nabla \underline{\hat{p}} \, dx + \ell \left( \int_{\Omega} \hat{\theta} \, dx - \Theta \right), \nonumber
\end{align}
$$

where $J^*$ is assumed to be a sufficiently “smooth” function (otherwise we cannot apply the adjoint-state method [35,43]).

4.2. Sensitivity analysis for the three-physics problem

In order to evaluate the sensitivity, the corresponding Lagrangian $\mathcal{L}$ (21) is differentiated with respect to the state variables $u, p$ and $T$ in the directions $\phi_u \in H^1(\Omega)^N$, $\phi_p \in H^1(\Omega)$ and $\phi_T \in H^1(\Omega)$; let $(u, u) \in H_0^1(\Omega; \mathbb{R}^N)^2$ be a stationary point of $\mathcal{L}$. Then, the derivative of $\mathcal{L}$ (21) with respect to $u$, in the direction $\phi_u \in H^1(\Omega)^N$ is defined by

$$
\left\langle \frac{\partial \mathcal{L}}{\partial u}, \phi_u \right\rangle = \left\langle \frac{\partial J^*}{\partial u}, \phi_u \right\rangle + \int_{\Omega} \left( \left\langle \left( \frac{\partial \sigma(\hat{u}, \hat{T})}{\partial u} \right), \phi_u \right\rangle : e(\hat{u}) \right) dx, \quad (22)
$$

while the derivative with respect to $p$, in the direction $\phi_p \in H^1(\Omega)$ is
defined by,

$$
\begin{equation}
\left\langle \frac{\partial \mathcal{L}}{\partial p}, \phi_p \right\rangle = \left\langle \frac{\partial J^*}{\partial p}, \phi_p \right\rangle + \int_{\Omega} b \nabla \phi_p \cdot \underline{\hat{u}} \, dx + \int_{\Omega} (-c_p \underline{\hat{T}} K_D^* \nabla \phi_p) \cdot \nabla \hat{T} \tag{23} \\
\qquad + \int_{\Omega} K_D^* \nabla \phi_p \cdot \nabla \underline{\hat{p}} \, dx
\end{equation}
$$

and the derivative with respect to $T$, in the direction $\phi_T \in H^1(\Omega)$ is
defined by:

$$
\begin{align}
\left\langle \frac{\partial \mathcal{L}}{\partial T}, \phi_T \right\rangle &= \left\langle \frac{\partial J^*}{\partial T}, \phi_T \right\rangle + \int_{\Omega} \left\langle \frac{\partial \sigma}{\partial T}, \phi_T \right\rangle : e(\underline{\hat{u}}) \, dx \tag{24} \\
&\quad + \int_{\Omega} K^* \nabla \phi_T \cdot \nabla \underline{\hat{T}} \, dx \nonumber \\
&\quad + \int_{\Omega} (-c_p \underline{\hat{T}} K_D^* \nabla \hat{p}) \cdot \nabla \phi_T \, dx, \nonumber
\end{align}
$$

which, when equations (22), (23) and (24) vanish, are the variational formulation of the adjoint-states; and the derivatives with respect to $\underline{u}$, $\underline{p}$ and $\underline{T}$, in the directions $\phi_u \in H^1(\Omega)^2$, $\phi_p \in H^1(\Omega)$ and $\phi_T \in H^1(\Omega)$ are the variational form of the state equations:

$$
\left\{
\begin{aligned}
\left\langle \frac{\partial \mathcal{L}}{\partial \underline{u}}, \phi_u \right\rangle &= \int_{\Omega} (\sigma(\hat{u}, \hat{T}) : e(\phi_u) + b(\hat{\theta}) \nabla \hat{p} \cdot \phi_u) \, dx, \\
\left\langle \frac{\partial \mathcal{L}}{\partial \underline{p}}, \phi_p \right\rangle &= \int_{\Omega} K_D^* \nabla \hat{p} \cdot \nabla \phi_p \, dx, \\
\left\langle \frac{\partial \mathcal{L}}{\partial \underline{T}}, \phi_T \right\rangle &= \int_{\Omega} K^* \nabla \hat{T} \cdot \nabla \phi_T + \int_{\Omega} (-c_p \phi_T K_D^* \nabla \hat{p}) \cdot \nabla \hat{T} \, dx,
\end{aligned}
\right.
\tag{25}
$$

Note that, in Eq. (21), the integration by parts has removed boundary terms (involving the source terms $f_s$, $h_0$, $f_f$), so that equations in (25), when vanishing, are consistent with homogeneous problems. Finally, the partial derivative of $\mathcal{L}$ with respect to $\theta$, in the direction $\bar{\theta} \in L^\infty(\Omega; \mathbb{R})$ at the stationary point $(u, \underline{u}, p, \underline{p}, T, \underline{T})$ is defined by

$$
\begin{equation}
\left\langle \frac{d\mathcal{L}}{d\theta}, \bar{\theta} \right\rangle = \left\langle \frac{\partial J^*}{\partial \theta}, \bar{\theta} \right\rangle + \int_{\Omega} \left( \frac{\partial \sigma}{\partial \theta} : e(\underline{u}) + \frac{\partial b}{\partial \theta} \nabla p \cdot \underline{u} + \frac{\partial K_D^*}{\partial \theta} \nabla p \cdot \nabla \underline{p} \right. \\
\left. \qquad + \frac{\partial K^*}{\partial \theta} \nabla T \cdot \nabla \underline{T} + (-c_p \underline{T} \frac{\partial K_D^*}{\partial \theta} \nabla p) \cdot \nabla T + \ell \right) \bar{\theta} \, dx \tag{26}
\end{equation}
$$

Note that in Eq. (26), the term $\langle \frac{\partial J^*}{\partial \bar{\theta}}, \bar{\theta} \rangle$ is the partial derivative of the objective function $J^*$ with respect to $\theta$, in the direction $\bar{\theta}$, while the term

$$
\int_{\Omega} \left( \frac{\partial \sigma}{\partial \theta} : e(\underline{u}) + \frac{\partial b}{\partial \theta} \nabla p \cdot \underline{u} + \frac{\partial K_D^*}{\partial \theta} \nabla p \cdot \nabla \underline{p} + \frac{\partial K^*}{\partial \theta} \nabla T \cdot \nabla \underline{T} \right. \\
\left. + (-c_p \underline{T} \frac{\partial K_D^*}{\partial \theta} \nabla p) \cdot \nabla T + \ell \right) \bar{\theta} \, dx
$$

is used to determine the load sensitivities of the constraints with re-
spect to $\theta$; here, $\underline{u}$, $\underline{p}$ and $\underline{T}$ are the adjoint-state variables, solutions
to the adjoint equations (22)-(24) (which we have to solve first), and
$\ell$ is the Lagrange multiplier associated to the volume constraint: it is
numerically adjusted at each iteration to fulfill the volume constraint.
Note that, the homogenized tensor $A^{*}(\theta)$ has been computed for both
periodic cells, for discrete values $(\theta_i)_{i=1,...,n}$ of $\theta$, using formula (6) (a
preprocessing stage, which is performed offline), and a linear interpola-
tion approach is proposed to update the derivative with respect to $\theta$ of
$\sigma$: $\sigma'(\theta) = A^{*'}(\theta) : e(u)$; the derivative with respect to $\theta$ of $K^*(\theta)$, $K_D^*(\theta)$
and the Biot coefficient $b(\theta)$, are approximated by some smooth func-
tions: see Equations (11) and (16).

5. Topology optimization of modulated periodic composite materials

We now propose our numerical algorithm and describe how the
methodology applies to the weakly coupled three physics system of (12),
(14) and (15).

5.1. Alternate directions method

The optimization problem (17) is solved based on the alternate di-
rections algorithm [32] (see Appendix C, for more details).

5.1.1. Minimizing over the stress field

Assume the design field $\theta$ is given, then the minimization over the stress field $\sigma$ comes down to solving the linear poro-thermo elasticity problem (15), with a material of elasticity tensor $A^*(\theta)$. The composite design $(\theta, A^*(\theta))$ is then computed by linear interpolation over the table (i.e., the surrogate model of $A^*(\theta)$) obtained after the preprocessing stage (performed offline); see Section 2.2.1.

5.1.2. Minimizing over the density field

Assume the stress field $\sigma$ is given, then the minimization over the
density field $\theta$ is done using the projected gradient algorithm. And
since, the considered optimization problem (17) is not self-adjoint, it is then essential to define and solve the associated adjoint problem; this is performed by computing the descent direction, which is given by the derivative of $\mathcal{L}$ with respect to $\theta$ (26), where the descent direction $\bar{\theta} = d\theta$ has to be selected such that:

$$ \left\langle \frac{\partial \mathcal{L}}{\partial \theta} (\theta, u, \underline{u}, p, \underline{p}, T, \underline{T}, \ell), d\theta \right\rangle < 0, \qquad (27) $$

which is achieved by choosing

$$ d\theta = - \left( \frac{\partial J_1^*}{\partial \theta} (\theta) + \frac{\partial \sigma}{\partial \theta} : e(\underline{u}) + \frac{\partial b}{\partial \theta} \nabla p \cdot \underline{u} + \frac{\partial K_D^*}{\partial \theta} \nabla p \cdot \nabla \underline{p} \, dx \right. \\ \left. + \frac{\partial K^*}{\partial \theta} \nabla T \cdot \nabla \underline{T} \right. \\ \left. + \left( -c_p T \frac{\partial K_D^*}{\partial \theta} \nabla p \right) \cdot \nabla T + \ell \right) \qquad (28) $$

where $\langle \frac{\partial J_1^*}{\partial \theta}, \bar{\theta} \rangle = \int_{\Omega} \frac{\partial J_1^*}{\partial \theta} \cdot \bar{\theta} \, dx$. At iteration $n$, we then update the optimal design $\theta$ by computing the projected gradient:

$$ \theta^{n+1} = P_{[0,1]}(\theta^n - \delta d\theta), \qquad (29) $$

where $P_{[0,1]}$ is the projection operator in $[0, 1]$ and $\delta > 0$ is the step size; here, $\delta$ is an adaptive step: at each iteration, if the newly computed design is accepted (i.e., if the current objective function $J^*(\theta^n)$ is lower than previous one $J^*(\theta^{n-1})$), $\delta$ is increased by 20%, else, if it is rejected, $\delta$ is divided by 2. In practice, the partial derivative of the Lagrangian $\frac{\partial \mathcal{L}}{\partial \theta}$ is regularized using an equivalent $H^1$-norm:

$$ \int_{\Omega} \left( \frac{\partial \mathcal{L}}{\partial \theta} \bar{\theta} + \eta^2 \nabla \frac{\partial \mathcal{L}}{\partial \theta} \cdot \nabla \bar{\theta} \right) dx = \left\langle \frac{\partial J^*}{\partial \theta}, \bar{\theta} \right\rangle \qquad (30) \\ + \int_{\Omega} \left( \left( \frac{\partial \sigma}{\partial \theta} : e(\underline{u}) + \frac{\partial b}{\partial \theta} \nabla p \cdot \underline{u} \right) \right. \\ \left. + \frac{\partial K^*}{\partial \theta} \nabla T \cdot \nabla \underline{T} + \left( -c_p \frac{\partial K^*}{\partial \theta} \nabla p \right) \cdot \nabla T \underline{T} + \ell \right) \bar{\theta} \, dx, $$

where $\eta$ is a small coefficient, which depends on the size of the elements of the mesh: thanks to $\eta$, we are able to numerically regularize the partial derivative on a length scale of order $\eta$ and to limit local instabilities on the density $\theta$ (e.g. checkerboard effect).

### 5.1.3. Complete optimization algorithm

The alternate directions algorithm is an iterative method, which is structured as follows in Algorithm 1.

**Algorithm 1** Alternate directions algorithm.

1: For $n = 0$, initialize $\theta^0$ such that, it satisfies the volume constraint $\Theta$:

$$ \forall x \in \Omega, \quad \theta^0(x) = \frac{\Theta}{\int_{\Omega} 1 \, dx} $$

2: Iteration until convergence, for $n \ge 0$:

1: Computation of the pressure $p^n$ through Biot-Darcy model (12), using the composite design $\{\theta^n, A^*(\theta^n)\}$

2: Computation of the temperature $T^n$ through the convection-diffusion equations (14), using the composite design $\{\theta^n, A^*(\theta^n)\}$

3: Computation of the stress tensor $\sigma^n$ through the linear poro-thermo elasticity (15), using the composite design $\{\theta^n, A^*(\theta^n)\}$ and the descent direction $d\theta^n$, for a given stress tensor $\sigma^n$, using formulas (28)

4: Update the design $\theta^{n+1}$, using formulas (29) for the descent direction $d\theta^n$, and then update the composite design $\{\theta^{n+1}, A^*(\theta^{n+1})\}$, by linear interpolation.

<table><caption>**Table 1**
Numerical values of the physical parameters for test case of Section 6.1.</caption><thead><tr><td>L (m)</td><td>f<sub>f</sub> (kg/s)</td><td>p<sub>0</sub> (bar)</td></tr></thead><tbody><tr><td>0.1</td><td>1.5</td><td>0</td></tr></tbody></table>

## 6. Numerical results and discussion

Our numerical results in the three physics context detailed in Section 3 is presented. A variety of 2-D and 3-D test cases are introduced to demonstrate that the alternate directions Algorithm 1 produces physically correct results; FreeFem++ software [51] is used here to implement the algorithm, where all the unknowns are discretized using $P_1$ finite elements. For all finite element operations a typical mesh consisting of 15228 triangular elements (in 2-D) and 119172 tetrahedral elements (in 3-D) were used. The void (i.e., $\theta = 0$) is replaced with a very compliant material: here, $\theta$ is fixed at $10^{-4}$, in order to avoid singularities of the effective tensor when the elasticity problem is solved.

We propose four test cases, where two out of four are new to the best of our knowledge: the first two are 2-D examples and the last two are in 3-D. First, we propose a fluid-mechanical problem (where the convection-diffusion model (14) is not taken into account), second, a convective heat transfer (where the linear elasticity model (15) is not taken into account) and, for the remaining test cases, a 3-D analogue of the second example is provided and 3-D example, where all the three physics are involved.

### 6.1. A fluid-structure interaction problem

In this example, the bounding box of the structure is a square of dimension $L$, fixed on its boundary $\Gamma_D^s$ (of length $L/10$) and subject to surface loading $f_s$ on its boundary $\Gamma_N^s$ (of length $L/15$). In addition, a mass flow rate $f_f$ is imposed on its Neumann part $\Gamma_N^f$ (of length $L/5$), while submitted to a pressure load $p_0$ on its Dirichlet part $\Gamma_D^f$ (of length $L/5$). On the free interface $\Gamma^f$, zero normal stress is applied. See Fig. 6 for a sketch of this test case.

Here, we neglect the thermal effects (i.e., Eq. (14) is ignored), so that Eq. (15), becomes a poro-linear elasticity system with the forcing induced by the mechanical load and the fluid. Our aim is to achieve a trade-off between the minimization of the compliance of the solid body and the maximization of the hydraulic power transferred by the fluid

$$ J^*(\theta, p(\theta), u(\theta)) = \omega \left\langle \underbrace{\int_{\Omega} A^* e(u) : e(u) \, dx}_{\text{Elastic strain energy}} + (1-\omega) \left( \underbrace{-\int_{\Omega} K^* \nabla p \cdot \nabla p \, dx}_{\text{Hydraulic power}} \right) \right\rangle \qquad (31) $$

for some fixed coefficient $\omega \in [0, 1]$; it measures the relative weight given to each term in Eq. (31). The objective function $J^*$ is subjected to a volume constraint $\Theta = 44\%$:

$$ \frac{1}{|\Omega|} \int_{\Omega} \theta \, dx = \Theta $$

In this example, the numerical values of the considered physical parameters are supplied by Table 1; we consider two cases in Eq. (31): either $\omega = 1$ or $\omega = 2/5$. Note that, unfortunately the fixed coefficient $\omega$ is not optimal and so, the terms composing the objective (31) are very heterogeneous: several values of $\omega$ were considered (e.g. $\omega \in [0, 1]$) but we did not achieve a better results (or trade-off).
<figure><img src="image_9.png" /><figcaption>Fig. 6. Setting for the fluid-structure interaction problem of Section 6.1. The boundary parts $Γ_D^f, Γ_N^f$ and $Γ_N^u$ are centered.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 7. The optimal density for the two-physic problem of Section 6.1, in both cases.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 8. (a) State variable $p$ (bar) and (b) the solid displacement (amplified by a factor 3) for the two-physic problem of Section 6.1, with $\omega = 2/5$.</figcaption></figure>

On Fig. 7, the optimal density is displayed, while on Fig. 8, the corresponding pressure field and solid displacement are displayed.

The optimal density $\theta$ is represented in a gray scale: regions where $\theta = 1$ are black (pure material), whereas white regions correspond to voids. The gray regions correspond to the composite design. The topology of the solution is quite similar to that expected prior to the optimization. We note that, the algorithm tends to distribute more material in regions where the pressure gradient is high, of which we can clearly see a silhouette, although contains a large composite design. However, in the case where $\omega = 2/5$, the topology is more diffuse: it seems to be driven by the hydraulic power term. For the pressure field $p$, we emphasize that the boundary conditions are respected: regions with high pressure gradient are located mainly on the silhouette “shape” and its becomes diffuse elsewhere. For this latter, the forcing induced by the fluid and the mechanical load are displayed by the displacement of solid $u$; see Fig. 8.

<figure><img src="image_12.png" /><figcaption>Fig. 9. The convergence history of the objective function $J^*$, for the two-physic topology optimization problem of Section 6.1, in both cases.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 10. Setting of the convective heat transfer test case of Section 6.2.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 11. The optimal density for the two-physic topology optimization problem of Section 6.2.</figcaption></figure>

On Fig. 9, we plot the convergence history for this calculation: with $\omega = 2/5$, the objective $J^*$ turns out to be very sensitive with respect to very small variations of the topology and increases (in the first part of the optimization) due to the fact that the volume constraint is not yet satisfied.

## 6.2. A convective heat transfer problem

This second example is similar to that introduced and solved in [11], where the level-set method and a Stokes model for the fluid flow is employed. The computational domain is a square of dimension $L$, where a fluid is entering with a given mass flow rate $f_f$ on its Neumann part $Γ_N^f$ (of length $L/5$), and while submitted to a pressure load $p_0$ on its Dirichlet part $Γ_D^f$ (of length $L/5$). In addition, a temperature $T_1$ is entering through
<caption>**Table 2**</caption>

Numerical values of the physical parameters in the convective heat transfer problem of Section 6.2.

<table><thead><tr><th>L (m)</th><th>f<sub>f</sub> (kg/s)</th><th>p<sub>0</sub> (bar)</th><th>T<sub>in</sub> (C)</th><th>T<sub>1</sub></th><th>c<sub>p</sub> (J/kg.K)</th><th>k<sub>s</sub></th><th>k<sub>f</sub> (W/m.K)</th></tr></thead><tbody><tr><td>0.1</td><td>1</td><td>0</td><td>0</td><td>10</td><td>100</td><td>10</td><td>1</td></tr></tbody></table>

<figure><img src="image_15.png" /><figcaption>(a) Pressure field p</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(b) Temperature T</figcaption></figure>

<caption>**Fig. 12.** Pressure *p* (bar) and temperature *T* (C) for the optimized configuration of the two-physic problem of Section 6.2.</caption>

the upper and lower walls (i.e., on its Dirichlet part $\Gamma_D^T : T = T_1$), while an inlet temperature $T_{in}$ is flowing through the middle left part of its boundary (i.e., on its Dirichlet part $\Gamma_D^T : T = T_{in}$ of length $L/5$). On the free interface $\Gamma = \Gamma^f = \Gamma^T$, zero normal stress is applied for both the temperature and the pressure; see Fig. 10 for a sketch of this test case.

For this example, we neglect the elastic contribution (i.e., Eq. (15) is ignored), so that the example involves only a coupling of the fluid (12) and the heat equation (14); here, our aim is to achieve a trade-off between the maximization of the hydraulic power dissipated by the fluid and the maximization of the heat convected by the fluid:

$$J^*(\theta, p(\theta), T(\theta)) = \omega \left( \underbrace{- \int_{\Omega} K^* \nabla p \cdot \nabla p \, dx}_{\text{Hydraulic power}} \right) + (1-\omega) \left( \int_{\Omega} c_p K_D^* \nabla p \cdot \nabla T \, dx \right), \quad (32)$$

for some fixed coefficient $\omega \in [0, 1]$ and a thermal capacity $c_p$ (Jkg$^{-1}$K$^{-1}$). The objective function $J^*$ is minimized with a volume constraint $\Theta = 44\%$, and the balance coefficient is set to $\omega = 1/2$. The numerical values of the considered physical parameters are supplied by Table 2.

The optimal density, the corresponding pressure $p$ and temperature $T$ are displayed on Fig. 11 and Fig. 12. The topology of the result is quite similar to that obtained in [11] (where the obtained solution is autopedalized). We note that, the algorithm tends to distribute more material in regions where the pressure gradient is low and temperature gradient is high. We can clearly see a silhouette of a pipe. For the pressure $p$, we emphasize that the boundary conditions are respected: regions with high pressure gradient are located mainly at the center of the pipe and becomes diffuse elsewhere. For the temperature $T$, we notice a cool down of the heat at the center of the pipe and becomes hot elsewhere, with high temperature gradient; see Fig. 12.

On Fig. 13, we plot the convergence history for this calculation: smooth and relatively rapid convergence is observed.

### 6.3. 3-D thermal fluid interaction problem

This example is an attempt to compute a more realistic structure; it is the 3-D analog of the second example of Section 6.2. The computational domain is a 3-D cubic box of dimension $L$, where a fluid is entering with a given mass rate $f_f$ on its Neumann part $\Gamma_N^f$ (of size $L/5 \times L/5$). A temperature $T_1$ is imposed on the upper and lower walls $\Gamma_D^T$ and a

<figure><img src="image_17.png" /><figcaption>**Fig. 13.** Convergence history of the objective function $J^*$ for $\omega = 1/2$, for the two-physic topology optimization problem of Section 6.2.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Fig. 14.** The setting of the 3-D two-physic problem of Section 6.3, subjected to fluid-thermal loads.</figcaption></figure>

temperature $T = T_{in}$ is imposed on the middle back and left walls $\Gamma_D^T$ (of size $L/5 \times L/5$). On $\Gamma = \Gamma^f = \Gamma^T$, zero normal stress is imposed for both the temperature and the pressure; see Fig. 14 for a sketch of this test case.

We recall that, for this example, the elastic model (15) is ignored. Thus, it involves only a coupling of the fluid (12) and heat equation (14).

The objective is to achieve a trade-off between the maximization of the hydraulic power dissipated by the fluid and the maximization of the heat convected by the fluid. See Eq. (32), for the optimization problem.
<figure><img src="image_19.png" /><figcaption>Fig. 15. Cross-sectional views of the solution of the 3-D thermal fluid interaction problem of Section 6.3. Composite solution represented as the part of the domain filled with values of $\theta \ge 0.5$.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Fig. 16. The setting of the 3-D three-physics problem of Section 6.4, subjected to fluid-thermal loads. On the red rectangles, the solid is clamped.</figcaption></figure>

The functional $J^*$ (32), is minimized under the volume constraint $\Theta = 44\%$. The numerical values of the considered physical parameters are supplied by Table 2. We used for all finite element operations a mesh consisting of 119172 tetrahedral elements and started from an initial configuration $\theta^0 = \Theta$, throughout. The algorithm produces a symmetric layout. The total FreeFEM running CPU time on Intel Xeon 2.60 GHz for this calculation is 5331 seconds.

Fig. 15 represents the iso-surface (or threshold) of the domain according to the optimized density $\theta$. In this example, the iso-surfaces are smooth and embedded into each other as $\theta$ increases. The topology of the solution is quite similar to that obtained in 2-D, namely, a 3-D pipe-like “shape” is observed, although it contains a large composite zone.

### 6.4. 3-D thermal fluid-structure interaction problem

We finally turn to a topology optimization example in the full three-physics setting 3; the computational domain is a 3-D cubic box of dimension $L$, fixed on its boundary $\Gamma_D^s$ (of size $\frac{L}{10} \times \frac{L}{10}$). A fluid is entering with a given mass rate $f_f$ on its Neumann part $\Gamma_N^f$ (of size $\frac{L}{5} \times \frac{L}{5}$). A temperature $T_1$ is imposed on the upper and lower walls $\Gamma_D^T$ and a temperature $T = T_{in}$ is imposed on the middle back and left walls $\Gamma_D^T$ (of size $\frac{L}{5} \times \frac{L}{5}$). On $\Gamma = \Gamma^f = \Gamma^T$, zero normal stress is applied for both the temperature and the pressure; see Fig. 16 for a sketch of this test case.

Here, the objective is to achieve a trade-off between minimization of the compliance of solid, the maximization of the hydraulic power dissipated by the fluid and the maximization of the heat convected by the fluid:

$$ J^*(\theta, p, T, u) = \omega_1 \left( \underbrace{\int_{\Omega} A^* e(u) : e(u) \, dx}_{\text{Elastic strain energy}} + \omega_2 \left( -\int_{\Omega} K_D^* \nabla p \cdot \nabla p \, dx \right) \right) + \omega_3 \left( \int_{\Omega} c_p K_D^* \nabla p \cdot \nabla T \, dx \right), \quad (33) $$

for are some given coefficients $\omega_j \in [0, 1]$ (with $j \in \{1, 2, 3\}$), such that $\omega_1 + \omega_2 + \omega_3 = 1$: we consider two different sets of coefficients $\omega_1 = \omega_2 = \omega_3 = 1/3$ and $\omega_1 = 1/5$, $\omega_2 = \omega_3 = 2/5$. The volume constraint is set to $\Theta = 30\%$.

On Fig. 17, we display the iso-surface $\theta \ge 0.3$ of the density for both cases, and on Fig. 18, the corresponding solid displacement is displayed. Very interestingly, we retrieve the fact that the topology of the solution contains a large composite zone and a large contact surface with the fluid at entrance, so to reduce the effect of the pressure source. In the case where $\omega_1 = 1/5$, $\omega_2 = \omega_3 = 2/5$, the topology of the solution is homogeneous and seems to achieve a trade-off between the minimization of the compliance (of the solid) induced by the fluid, the maximization of the hydraulic power and of the heat convected by the fluid.

### 6.5. Discussion

We emphasize that the obtained optimized structures for real isotropic microstructures are not autopenalized (except for the convective heat transfer Problem 6.2); this suggests that the microstructures are optimal for the objective function $J^*$. In practice, the local microstructure has to be designed in order to support the stress in privileged directions, namely the principal directions of the stress. However, we note that in our numerical results, the local microstructure is able to support stresses in several directions: this means that the proposed isotropic cells could then constitute a great advantage. In 3-D, the optimized structure is clearly not a black and white design: high densities are reached only around where Dirichlet and Neumann boundary conditions are applied, and the remaining part of the domain is filled with intermediate or grey densities. Hence, the considered isotropic microstructures are adapted to our three physics problems. However, other kind of microstructures, not necessarily isotropic, should be investigated.

## 7. Conclusions

In this work, we demonstrated that the alternate directions method applied to the weakly coupled three physics system produces physically correct results for unconventional design problems. In our methodology,
<figure><img src="image_21.png" /><figcaption>(a) iso-surface $\theta \ge 0.3$</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>(b) $\omega_1 = \omega_2 = \omega_3 = 1/3$</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>(c) $\omega_1 = 1/5, \omega_2 = \omega_3 = 2/5$</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>**Fig. 18.** The solid displacement for the three-physic topology optimization problem of Section 6.4, in the case $\omega_1 = \omega_2 = \omega_3 = 1/3$ (displacement amplified by a factor 3).</figcaption></figure>

Biot-Darcy's law and the convection-diffusion equation are proposed to describe the pressure (of the fluid flow) and the temperature (of the heat flux) in the fluid domain, which we weakly coupled to the solid phase by solving the associated PDEs using the standard finite element method. Furthermore, smooth enough functions have been introduced in order to treat the transition between the solid and void phases, which are related to the material density. Here, the adjoint-state method is proposed to compute the load sensitivities (with respect to the design variables) and, it is noticed that the consideration of load sensitivities within the methodology is essential when designing multi-physic systems. We emphasize that, the Biot-Darcy's model, together with our simplified heat

transfer equation allow for relatively straightforward extension to 3-D problems.

### CRediT authorship contribution statement

Godfred Oheneba Agyekum: Writing – original draft. Laurent Cangémi: Supervision. François Jouve: Supervision.

### Replication of results

On behalf of all authors, the corresponding author states that there is no codes available as supplementary material because this work is carried out as part of a thesis in partnership with IFPEN. The latter holds the intellectual rights.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgements

This research project was fully sponsored by IFPEN: IFP energies nouvelles; and was supported by Pierre Viot: design and simulation engineer at IFPEN.

## Appendix A. Composite designs to graded microstructures

As aforementioned in the introduction 1, topology optimization involves generalized shape optimization problems in which the degree of
<figure><img src="image_25.png" /><figcaption>Fig. A1. Presentation of an octahedron cell, composed of six identical square and eight hexagonal faces.</figcaption></figure>

topological connectedness (or equivalently, the number of holes in 2D) is not prescribed a priori and may vary during the optimization process. For this reason, we adopt here an approach developed by Bendsoe and Kikuchi [29]. By this mean, we define a new equivalent full media, with a behavior similar to that of lattice materials, where the geometry of the chosen lattice is characterized by only one parameter (the local density). In the following, we present the geometry of the chosen representative unit cell in 3-D, and then, we recall some limitations of topology optimization and how to possibly resolve them.

## A.1. Kevin foam cell

A Kevin foam cell, also denoted as tetrakaidecahedron cell is a truncated octahedron as displayed in Fig. A2(c). Here, a regular tetrakaidecahedron unit cell is designed; it is shown to yield isotropic (or quasi isotropic) homogenized Hooke's laws. It features fourteen faces: six identical square and eight regular hexagonal faces formed by a network of 36 ligaments. Each of the 24 vertices rigidly connects the ends of coinciding ligaments. The cell presents an isometry and cubic material symmetry, where the dihedral angle is approximately 109.471° at edges shared by two hexagons or 125.263° at edges shared by a hexagon and a square faces; see Fig. A2, for an example of construction of Kelvin foam cell, which is based on a truncated octahedron (Fig. A1).

Several approaches are possible in order to design isotropic cells from this one, parametrized by one parameter (i.e., the local density, which is characterized by the ratio of the homothety between the hole in a face and the external border of this face). Here, a naive approach is proposed, which consists in adding a central tetrakaidecahedron inclusion characterized by its size (see Fig. A2(a)). Based on the construction method described above, the material density $\theta$ has been used as a control parameter to generate foams of chosen porosities. Upon increasing $\theta$, one would reach a point where the structure is full of material.

## A.2. Deshomogenization process

Unfortunately, the approach presented in [29] is not straightforward: two limitations put it aside progressively in favor of the SIMP method. First, the lattice structures are not easy to manufacture using the traditional fabrication processes, such as foundry, because of their microscopic structure and their graded (or complex) microstructure. However, thanks to the recent improvements in additive manufacturing processes, this problem is partly solved. Graded microstructures (or lattice materials) are indeed now manufacturable. But many researches are currently

being carried out in order to improve the process and model it, particularly from a thermal point of view in order to characterize the residual stresses or possible defects in the structure. Nevertheless, this led to a resumption of the homogenization based method.

The second limitation is about the deshomogenization process, because the homogenization method yields a generalized shape (or composite), which is not a genuine shape, as in the SIMP or level-set method (straightly manufacturable). Therefore, we have to construct from composite (or effective) design to graded microstructure (or genuine shape). This deshomogenization process replace the traditional penalization step. Fortunately, when the orientation of the lattice cells is uniform, this is not a big challenge: the periodic size of the deshomogenized cells is chosen, then on each cell we use the average parameters to design the local microstructure. However, the connectivity between adjacent cells has to be carefully ensured. But when there is need to optimally orientate the cells locally in the domain, then the orientation has to be taken into account during the deshomogenization process. Geoffroy-Donders [44,52], proposed two main approaches to construct a sequence of genuine shapes that actually converges toward the homogenized structure:

* First method is a local approach and leads to regular cells, in the sense that each cell is exactly that described in section 2.1. However, this local approach is very time consuming since each cell has to be constructed separately.

* Second method is a global one, the final structure is described through a level-set function and is really fast to be computed.

The mathematical framework of the deshomogenization method can be found in [44,53].

# Appendix B. The regime analysis

## B.1. Domain of validity

It is worth mentioning that the presented system in (15), is a simplified version of a genuine three-physics model:

$$ egin{cases} -\mathrm{div}(\bar{\sigma}) = 0 & \text{in } \Omega, \\ \bar{\sigma} = \sigma(u, T) - bpI, \\ \sigma(u, T) = A^*(e(u) - \alpha(T - T_{ref})I), \\ p = Mm - Mbe_{vol} \end{cases} $$

where $p$ satisfies a hydraulic model, which features a transition regime and an inertia regime:

$$ q = -K_D^* \nabla p = mv_f + \rho C_{for} v_f^2 \quad \text{in } \Omega, $$

where $C_{for}$ is the so called Forchheimer coefficient: an inertia parameter of the fluid flow. Regarding the domain of validity, several characteristic length scales have been proposed in the literature, namely, various definitions in the form of morphological parameter (e.g. pore diameter) or hydraulic parameter (e.g. permeability) have been used as a characteristic length choice. No general consensus has been ever achieved on this matter and varies from one author to another. However in [47], the authors numerically demonstrate that, it is insightful to obtain the threshold range of Reynolds number $Red_p ( = \frac{\rho v_f d_p}{\mu_f})$ based on pore diameter $d_p$, which can be used to distinguish the flow regimes and choose the flow law accordingly. In their numerical simulations, the authors chose a Kelvin-like structure (i.e., tetrakaidecahedron cell in 3-D) as a computational domain and they numerically distinguished three regimes:

* $Red_p \le 0.3$: Darcy regime,

* $0.3 < Red_p \le 30$: Cubic regime,

* $Red_p > 30$: Weak inertia regime
<figure><img src="image_26.png" /><figcaption>Fig. A2. Construction of tetrakaidecahedron model of an isotopic Kelvin foam cell.</figcaption></figure>

with respect to circular strut cross-section having variable pore diameters and porosity for a given cell size: constant fluid properties were used. As aforementioned, we opted here for simplified regime, namely, the Darcy regime in an isotropic composite domain $\Omega$ (periodically perforated by hexagonal cells in 2-D and tetrakaidecahedron cells in 3-D). As such, the thermal or hydraulic parameter (e.g. thermal or pressure conductivity) can be identified to scalar values, which depend solely on the material density $\theta$.

## B.2. Strong couplings

Regarding the regimes of strong couplings with respect to those of weak ones, it is worth to recall that, from a modeling point of view, very often, the mechanical deformations and displacements are small. This implies that the fluid domain is fixed at first order, that is, independent of the deformation of the structure. It is therefore natural to consider only a weak coupling between the three physics at stake, and of course this weak coupling is a major simplification and it dramatically reduces the computational cost since no monolithic coupled system has to be solved. In contrast, a strong coupling between the three physics at stake is difficult to solve from a numerical point of view and in particular when solving the optimization problem using the alternate directions algorithm, because the resulting computations are very costly and there is no guaranty that our algorithm converge.

We emphasize that at $\theta = 0$, the stress field $\sigma = 0$ (for instance) vanishes, and which calls for special care in the gradient method. This lack of continuity at $\theta = 0$ is the mathematical manifestation of the presence of holes in the computational domain. In practice, the smallest admissible value of $\theta$ is fixed at $10^{-4}$, in order to avoid singularities of the effective tensor when the elasticity problem is solved (i.e., the ersatz method). We emphasize that at $\theta = 0$, the stress field $\sigma = 0$ (for instance) vanishes, and which calls for special care in the gradient method. This lack of continuity at $\theta = 0$ is the mathematical manifestation of the presence of holes in the computational domain. In practice, the smallest admissible value of $\theta$ is fixed at $10^{-4}$, in order to avoid singularities of the effective tensor when the elasticity problem is solved (i.e., the ersatz method).

## Appendix C. Alternate directions algorithm

In practice the alternate directions algorithm [32] is inspired by the homogenization method for the optimization of the compliance, for which we can formulate the relaxed problem as a double minimization problem: the minimization of the effective complementary energy (i.e., over the stress field $\sigma$, which amounts to solve the linear elasticity problem) for a given composite design ($\theta, A^*(\theta)$) and the minimization of the

compliance for a given stress field $\sigma$. We note that in all computations involving complementary energies, high degree finite elements have to be used for stress accuracy and the resulting computations are very costly and limited. In addition, the relaxed (or effective) complementary energy is not smooth at $\sigma = 0$, which calls for special care in the gradient method: the convergence to the minimum is usually fairly slow.

Hence, we prefer the alternate directions algorithm, which is based on two key ideas: first one is to consider the relaxed problem $J^*$ as a minimization problem not only for the stress field $\sigma$, but also for the structural parameters (i.e., $\theta$ and $A^*(\theta)$), and the second is not to try to minimize directly in the triplet of variables ($\sigma, \theta, A^*(\theta)$), but rather to adopt an iterative approach and minimize separately and successively in the design variables ($\theta, A^*(\theta)$) and in the field variable $\sigma$; here, we can not formulate the direct problem (which amounts to solve the weakly coupled thermal fluid-elastic problem) under this form (i.e., as in the case of pure-elastic problem) and as a consequence, there is no guaranty that the energy decreases at each step, but it is still very natural to proceed in this way.

## Data availability

The data that has been used is confidential.

## References

[1] Hammer VB, Olhoff N. Topology optimization of continuum structures subjected to pressure loading. Struct Multidiscip Optim 2000;19(2):85-92.

[2] Du J, Olhoff N. Topological optimization of continuum structures with design-dependent surface loading - part I: new computational approach for 2D problems. Struct Multidiscip Optim 2004;27(3):151-65.

[3] Fuchs MB, Shemesh NNY. Density-based topological design of structures subjected to water pressure using a parametric loading surface. Struct Multidiscip Optim 2004;28(1):11-9.

[4] Wiker N, Klarbring A, Borrvall T. Topology optimization of regions of Darcy and Stokes flow. Int J Numer Methods Eng 2007;69:1374-404. https://doi.org/10.1002/nme.1811.

[5] Zheng B, Chang CJ, Gea HC. Topology optimization with design-dependent pressure loading. Struct Multidiscip Optim 2009;38(6):535-43.

[6] Lee E, Martins JRRA. Structural topology optimization with design-dependent pressure loads. Comput Methods Appl Mech Eng 2012;40(48):233-6.

[7] Wang C, Zhao M, Ge T. Structural topology optimization with design-dependent pressure loads. Struct Multidiscip Optim 2016;53(5):1005-18.

[8] Li C, Xu C, Fox MD. Distance regularized level set evolution and its application to image segmentation. IEEE Trans Image Process 2010;19(12):3243-54.

[9] Li Z-m, Yu J, Yu Y, Xu L. Topology optimization of pressure structures based on regional contour tracking technology. Struct Multidiscip Optim 2018;58(2):687-700.

[10] Yaji K, Yamada T, Kubo S, Izui K, Nishiwaki S. A topology optimization method for a coupled thermal-fluid problem using level set boundary expressions. Int J Heat Mass Transf 2015;81:878-88. https://doi.org/10.1016/j.ijheatmasstransfer.2014.11.005.

[11] Feppon F, Allaire G, Bordeu F, Cortial J, Dapogny C. Shape optimization of a coupled thermal fluid-structure problem in a level set mesh evolution framework. Bol Soc Esp Mat Apl 2019;76(3):413-58.
[12] Picelli R, Neofytou A, Kim HA. Topology optimization for design-dependent hydrostatic pressure loading via the level-set method. *Struct Multidiscip Optim* 2019;60(4):1313–26.

[13] Andreasen CS, Sigmund O. Saturated poroelastic actuators generated by topology optimization. *Struct Multidiscip Optim* 2011;43:693–706.

[14] Andreasen CS, Sigmund O. Topology optimization of fluid-structure-interaction problems in poroelasticity. *Comput Methods Appl Mech Eng* 2013;258:55–62.

[15] Chen BC, Kikuchi N. Topology optimization with design dependent loads. *Finite Elem Anal Des* 2001;37(1):57–70.

[16] Yoon GH. Topological layout design of electro-fluid-thermal-compliant actuator. *Comput Methods Appl Mech Eng* 2012;209212(28):44.

[17] Marck G, Nemer M, Harion J-L. Topology optimization of heat and mass transfer problems: laminar flow. *Numer Heat Transf, Part B, Fundam* 2013;63:508539.

[18] Marck G, Privat Y. On some shape and topology optimization problems in conductive and convective heat transfers. In: Papadrakakis M, Karlaftis M, Lagaros N, editors. *OPTI 2014, an international conference on engineering and applied sciences optimization; June 2014*. p. 16401657.

[19] Dbouk T. A review about the engineering design of optimal heat transfer systems using topology optimization. *Appl Therm Eng* 2017;112:841854.

[20] Mendes E, Sivapuram R, Rodriguez R, Sampaio M, Picelli R. Topology optimization for stability problems of submerged structures using the TOBS method. *Comput Struct* 2022;259:106–685.

[21] Hübner D, Rohan E, Lukeš V, Stingl M. Optimization of the porous material described by the Biot model. *Int J Solids Struct* 2019;156:157–216–33.

[22] Kumar P, Frouws JS, Langelaar M. Topology optimization of fluidic pressure loaded structures and compliant mechanisms using the Darcy method. *Struct Multidiscip Optim* 2020;61:1637–55.

[23] Kumar P. TOPress: a Matlab implementation for topology optimization of structures subjected to design-dependent pressure loads. *Struct Multidiscip Optim* 2023;66:97.

[24] Neuber H. Theory of notch stresses: principles for exact calculation of strength with reference to structural form and material. *USAEC Office of Technical Information*; 1961. p. 4547.

[25] Hashin Z, Shtrikman S. A variational approach to the theory of the elastic behavior of multiphase materials. *J Mech Phys Solids* 1963;11:127–40.

[26] Lurie KA, Cherkvaev AV, Fedorov AV. Regularization of optimal design problems for bars and plates, I, II. *J Optim Theory Appl* 1982;37(4):499–522. pp. 523–543.

[27] Murat F, Tartar L. Calcul des variations et homogénéisation, les méthodes de l'homogénéisation théorie et applications en physique. In: Dir. études et recherches EDF, Eyrolles; 1985. p. 319–69. English translation. In: Kohn RV, editor. *Topics in the mathematical modeling of composite materials*. In: *Progress in Nonlinear Differential Equations and their Applications*, vol. 3. Boston: Birkhauser; in press.

[28] Kohn RV, Strang G. Optimal design and relaxation of variational problems, I. *Commun Pure Appl Math* 1986;39(1):113–37.

[29] Bendsøe MP, Kikuchi N. Generating optimal topologies in structural design using a homogenization method. *Comput Methods Appl Mech Eng* 1988;71(2):197–224.

[30] Suzuki K, Kikuchi N. A homogenization method for shape and topology optimization. *Comput Methods Appl Mech Eng* 1991;93(3):291–318.

[31] Murat F. Contre-exemples pour divers problèmes où le contrôle intervient dans les coefficients. *Ann Mat Pura Appl* 1997;112:49–68.

[32] Allaire G, Bonnetier E, Francfort G, Jouve F. Shape optimization by the homogenization method. *Numer Math* 1997;76:27–68.

[33] Bendsøe MP. Optimal shape design as a material distribution problem. *Struct Multidiscip Optim* 1989;1(4):193–202.

[34] Chen B-C, Kikuchi N. Topology optimization with design dependent loads. *Finite Elem Anal Des* 2001;37(1):57–70. https://doi.org/10.1016/S0168874X(00)00021-4.

[35] Allaire G. Shape optimization by the homogenization method. Springer: Applied Mathematical Sciences; 2002. p. 146.

[36] Bendsøe MP, Sigmund O. Topology optimization: theory, methods, and applications. Springer Science & Business Media; 2003.

[37] Lin CY, Kikuchi N, Hollister J. A novel method for biomaterial scaffold internal architecture design to match bone elastic properties with desired porosity. *J Biomech* 2004;43(5):623–36.

[38] Behrou R, Ghanem MA, Macnider BC, Verma V, Alvey J, Hong J, et al. Topology optimization of nonlinear periodically microstructured materials for tailored homogenized constitutive properties. *Compos Struct* 2021;266:113729. https://doi.org/10.1016/j.compstruct.2021.113729.

[39] Agyekum-Oheneba G, Cangémi L, Jouve F. Homogenization based topology optimization of fluid-pressure loaded structures using the Biot-Darcy model. *Optim Eng* 2023;1573–2924.

[40] Agyekum-Oheneba G. Mise en oeuvre d'une approche d'optimisation topologique "composite" poro-thermo-mécanique. PhD thesis. Université Paris Cité; 2022. Available from: https://theses.fr/2022UNIP7023.

[41] Demš K, Mróž Z. Variational approach by means of adjoint systems to structural optimization and sensitivity analysis. I. *Int J Solids Struct* 1983;19(8):677–92.

[42] Céa J. Conception optimale ou identification de formes, calcul rapide de la dérivée directionnelle de la fonction coût. *Math Model Numer Anal* 1986;30(6):371–402.

[43] Allaire G. Conception optimale de structures. Berlin, Heidelberg: Springer; 2007.

[44] Geoffroy-Donders P. Homogenization method for topology optimization of structures built with lattice materials. PhD thesis. Université Paris Saclay (COMUE); 2018.

[45] Abad KME, Khanoki AS, Pasini DF. Fatigue design of lattice materials via computational mechanics: application to lattices with smooth transitions in cell geometry. *Int J Fatigue* 2013;47:126–36.

[46] Michailidis G. Manufacturing constraints and multi-phase shape and topology optimization via a level-set method. PhD thesis. Ecole Polytechnique X; 2014.

[47] Jobic Y, Kumar P, Topin F, Occelli R. Transport properties of solid foams having circular strut cross section using pore scale numerical simulations. *Heat Mass Transf* 2018;54:2351–70.

[48] Batchelor G. An introduction to fluid dynamics. Cambridge University Press; 2000.

[49] Vafai K, Tien CL. Boundary and inertia effects on flow and heat transfer in porous media. *Int J Heat Mass Transf* 1981;24(6):195–203.

[50] Wang P, Vafai K, Liu DY. Analysis of the volumetric phenomenon in porous beds subject to irradiation. *Numer Heat Transf, Part A, Appl* 2016;70(2):567–80.

[51] Hecht F. New development in freefem++. *J Numer Math* 2012.

[52] Geoffroy-Donders P, Allaire G, Michailidis O, Pantz O. Coupled optimization of macroscopic structures and lattice infill. *Int J Numer Methods Eng* 2022;123(13):2963–85.

[53] Pantz O, Trabelsi K. A post-treatment of the homogenization method for shape optimization. *SIAM J Control Optim* 2008;47(3):1380–98.