<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. 2023. hal-04108556v2

HAL Id: hal-04108556

https://hal.science/hal-04108556v2

Preprint submitted on 1 Sep 2023

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
Homogenization Based Topology
Optimization of a Coupled Thermal
Fluid-Structure Problem

Godfred Oheneba Agyekum¹*, Laurent Cangémi¹†
and François Jouve²†

¹*IFP Energies nouvelles, Rueil-Malmaison, 92852, France.

²Université Paris Cité, Laboratoire Jacques-Louis Lions (LJLL),
Paris, F-75006, France.

*Corresponding author(s). E-mail(s):

godfred.a.o.ezail@gmail.com;

Contributing authors: laurent.cangemi@ifpen.fr;

francois.jouve@u-paris.fr ;

†These authors contributed equally to this work.

Abstract

Homogenization method is applied to topology optimization of a weakly
coupled three physics problem, where structures are made of periodically
perforated material. The microscopic periodic cell is macroscopically
modulated, where the design is characterized by the material density
and its homogenized Hooke's law at each point of the domain. The cou-
pling is weak because the three physics involved are solved consecutively:
first, the coupled fluid flow is determined using Biot-Darcy's law, second,
the thermal model using the convection-diffusion equation and third, the
three-physic problem by solving the linear poro-thermo elasticity system;
our aim is to optimize the homogenized formulation of this system. This
approach permits a computationally low cost of evaluation of load sensi-
tivities using the adjoint-state method. Numerical two-dimensional and
moderately large-scale three-dimensional two or three-physic problems
are presented using the alternate directions algorithm. It is demonstrated
how the implementation can address a variety of design problems.

**Keywords:** Topology optimization, multi-scale, periodic homogenization, porous medium, adjoint methods, fluid-structure interaction, convective heat-transfer
# 1 Introduction

The ambition to develop simulation methods making it possible to predict the integrity or properties of use (e.g. mechanical, diffusive, thermal, electromagnetic, vibratory, etc.) of structures (e.g. industrial or natural), materials or processes involved in the development of new advanced technologies is growing consistently. Herein, homogenization-based method is proposed to investigate shape optimization problems for a weakly coupled model of heat propagation, fluid flow and structure strain; and making it possible to consider a weak coupling between the three physics at stake because one can assume that the fluid domain is fixed at first order. However, we should bear in mind that this weak coupling is a major simplification and therefore reduces the computational cost.

A comprehensive overview of shape optimization with the homogenization method is provided by [1] and, for a general summary of the homogenization method, we refer the reader to [2–8] and references therein. It should be noted that, this approach provides a consistent way for computing effective material with microstructures (i.e., composite materials) and that, once the optimal composite is obtained by homogenization-based topology optimization method, we might need to dehomogenize the solution; see [9] for periodically perforated materials. The design method described in this paper is strongly inspired by the works mentioned above as well as being related to modern production techniques such as additive manufacturing. It is worth noting that, there is yet another alternative approach, which amounts to couple (or incorporate) the homogenization method (inside a domain) with a geometric approach (for moving the domain’s boundary) in order to combine advantages of both methods [10].

A typical shape optimization problem arising in this context involves an objective function, depending on the geometries of the fluid and solid subdomains and where the whole domain is described by a density function (i.e., material density) that can take on values in the interval [0, 1], which has to be minimized under some constraints (e.g. volume or mass constraints). This allows to compute the sensitivities with respect to design variables using the adjoint-state method [11, 12], introducing adjoint states. For these adjoint states, which are to be solved, it turns out that the coupling is reversed for the adjoint system: the elasticity is solved first, followed by the convection-diffusion equation and the fluid model.

Shape optimization that involve pressure-loaded or thermal fluid-loaded boundaries has been conducted also by [13–25]. In [13–17], the authors deduce the topology or layout based on boundary identification schemes: in general, based on a *priori* chosen threshold density (i.e., iso-density curves/surfaces are identified). In [13, 14, 17], the authors employed the iso-density method to identify the pressure loading facets: Bézier spline curves were used to describe the pressure-loaded facets. This allows in [13, 14], to evaluate the sensitivities with respect to design variables using the finite difference formulation and, in [17], to provide an analytical method to calculate load sensitivities. Note
that in [13, 14, 17], the considered sensitivities were restricted to only the pressure-loaded boundaries.

In contrast, the works in [15, 16, 18] do not account for load sensitivities within their topology optimization setting: in [15], the pressure-loaded facets are pre-defined and an additional set of variables is used, which are optimized along with the design variables; whereas in [16], an element-based search method is employed to identify the pressure-loaded facets and in [18], an algorithm based on digital image processing and regional contour tracking is proposed to generate the pressure loading surface.

On the other hand, in [19], the authors deduce the topology optimization based on binary structures method to design structures that consider buckling constraints and loaded by design-dependent fluid pressure loads: it adopts binary design variables and handles multiple constraints solved by an integer linear programming scheme, where sensitivity filtering method is proposed.

In contrast to boundary identification schemes, in [20–22], the authors deduce the topology based on level-set methods: an implicit boundary description is available that can be used to define the pressure load. In [20], the Distance Regularized Level-Set Evolution is proposed to capture the structural boundary and using the zero level contour of a level-set function to represent the loaded-pressure boundary but did not account load sensitivities; whereas in [21], the Laplace’s equation approach is employed to compute hydrostatic fluid pressure fields, and also a flood fill procedure to capture the solid/fluid interface: shape sensitivities in conjunction with ersatz material interpolation method are used within their approach. Recently in [22], Hadamard’s method of shape differentiation is applied to shape and topology optimization of a coupled thermal fluid-structure problem in a level-set mesh evolution framework: sensitivity analysis is performed with respect to the geometry of the interface between the fluid and solid domain, using the Hadamard’s method of shape differentiation, introducing adjoint states.

Moreover, unlike boundary identification schemes or level-set methods, in [23, 24], the authors deduce the topology using density-based approach: shape optimization problems are transformed to material distribution problems using fictitious composite materials without identifying loading surfaces directly. In [23], a density-based topology optimization is proposed to design both structures and compliant mechanisms loaded by design-dependent pressure loads: Darcy’s law in conjunction with a drainage term is proposed to treat the pressure loads, which are transferred into a design dependent pressure field using a partial differential equation, which is solved using the finite element method; the load sensitivities are computed using the adjoint-variable method. Recently in [24], is developed according to the approach first reported in [23], a MATLAB implementation TOPress, using the method of moving asymptotes. In contrast, the density-based method presented in [25] is based on true composite materials: two material constituents (substance and void) are considered and the microscopic optimal void distribution is considered. An important feature of the procedure is that the homogenization method is applied to
determine macroscopic constitutive equations for the material with micro-
scopic material constituents. In [25], the porous material is described as the
Biot continuum derived by the homogenization of two decoupled problems:
deformation of a porous solid saturated by a slightly compressible static fluid,
first, and Stokes flow through the rigid porous structure, second. The effective
medium (or composite) properties are given by the drained skeleton elastic-
ity, the Biot stress coupling, the Biot compressibility coefficients and by the
hydraulic permeability of the Darcy flow model: these are computed using
characteristic responses of the representative unit cell constituted by an elas-
tic skeleton (solid) and by the fluid channel (void); the adjoint-state method
is proposed to evaluate sensitivities of objective functions constituted by the
Biot model coefficients with respect to the underlying pore shape described
by a B-spline box which embeds the whole representative cell and, where the
gradient-based method is employed to solve the optimization problems: the
shape derivatives of the homogenized coefficients are derived using the shape
sensitivity technique and the material derivative approach.

In this article, we present a new approach to design structures subjected to
thermal loads and cooled down by a fluid-pressure. This falls within the gen-
eral framework of density methods where phenomenological laws of equivalent
media are derived by the homogenization method. In this regard, the pre-
sented approach falls within the framework of recent work in [23, 25], where
Darcy's flow model is used to describe the fluid flow. However, the approach
in [23] is different since the continuous problem does not contain a model
that explicitly couples the fluid pressure to the solid skeleton and thus, it
induces difficulties in the modeling where a volumetric force is added intu-
itively in the elastic problem without being explicitly defined as the result of
a continuous physical law. Moreover, compared to [25], the porous material
is described as the Biot-Darcy continuum derived by the homogenization of a
weakly coupled three physics problem: deformation of a porous-thermoelastic
saturated by incompressible fluid pressure (satisfying Biot-Darcy's law), first,
the convection-diffusion equation for the temperature, second, and cool down
by a Biot-Darcy's flow through the rigid porous thermoelastic, third. The
effective medium properties are given by the undrained skeleton elasticity, the
Biot stress coupling, the Biot coefficient, the thermal coefficients and con-
ductivity for the convection-diffusion model and, by hydraulic permeability of
the Darcy flow model: these are computed using characteristic responses of
the representative unit cell, namely, a perforated hexagonal cell in 2-D or the
tetrakaidecahedron cell in 3-D, constituted by a solid phase and void; see Fig. 1
and Fig. 3. Consequently, the fluid pressure and the temperature together
determine the displacement of the solid part in the porous medium. This allows
to provide sensitivities of general "smooth enough" objective functions with
respect to design variables using the adjoint-state method. The homogenized
coefficients are differentiated using the Lagrangian method and the projected
gradient algorithm: in [25], the shape sensitivity technique and the material
derivative approach are employed. Of course, these shape derivatives are at the
basis of our gradient-based alternate directions algorithm [26], which is used
for our numerical simulations. In summary, we present the following aspects:

* homogenized Biot-Darcy's law is used to characterize the fluid flow through a true composite material,

* the convection-diffusion equation is used for the temperature,

* we weakly couple the fluid loads to the linearized thermoelasticity system for the solid displacement,

* the approach facilitates computationally inexpensive evaluation of the load sensitivities with respect to design variables using the adjoint-state method,

* the flow coefficient, Biot's coefficient and conductivity coefficients are derived using the homogenization method,

* the approach avoids explicit description of the loading boundary,

* the robustness and efficiency of the approach is demonstrated through several design problems, using the alternate directions algorithm.

The remainder of this paper is structured as follows: in Section 2, we briefly recall the necessary ingredients of the homogenization method and we explain our strategy. First, choose a parametrized periodicity cell. Second, we compute its effective properties for the entire range of its parameters. In Section 3, we give a precise account of our weakly coupled model of heat propagation, fluid flow and structure strain. In Section 4, we introduce the optimization problem formulation, which turns out to be a simple parametric optimization problem since our periodicity cell is parametrized. Section 5, is concerned with our topology optimization algorithm: it is an alternate directions algorithm [26], which successively computes the stress field through the solving of a weakly coupled three physics problem. Finally, our numerical results are presented in Section 6: 2-D and 3-D computations are displayed of various design problems involving two or three-physic for some objective functions. We summarize our findings and give an outlook in Section 7.

## 2 Periodic homogenization

The mathematical framework of the homogenization theory can be found in [1]. In this section, we briefly present the principles of this method and explain our strategy. We restrict our analysis to locally periodic hexagonal cell in 2-D and tetrakaidecahedron cell in 3-D. Note that, it should be possible to adapt the whole method (or at least part of it) to periodic square cells (in 2-D) or to other similar geometry as the tetrakaidecahedron cells (in 3-D). Our aim is to determine the homogenized physical properties of such materials when varying their parameter; it is a preprocessing stage, which can be performed off-line. It is independent of the objective function, computational domain, applied loads or boundary conditions. The Hooke's laws are computed by solving the so-called cell problems, that describe the deformation at the scale of the microstructure, which is a very classical task in homogenization theory.
## 2.1 Set of admissible microstructures

From now on, we restrain our analysis to a simple class of composites already introduced in [9], repeated periodically on the whole space.

### 2.1.1 Smooth honeycomb cell

Here, our interest for the smooth honeycomb relies on its smooth rounded corners, known to generate lower local concentration stress ([31, 32]), compared to the classical honeycomb; see Fig. 1 and Fig. 2.

<figure><img src="image_3.png" /><figcaption>Fig. 1 Isotropic design cells (images taken from [9])</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 2 Smooth honeycomb cell for different values of the density $\theta$</figcaption></figure>

Let $Y$ be the smooth honeycomb cell in 2-D: a unit regular hexagonal in 2-D with smooth hexagonal central hole, similar to the classical honeycomb cell, except that the interior corners of the central hexagonal hole are rounded; see Fig. 1. Note that, the parameter $e$ displayed by Fig. 1 (a), represents the relative width of bars of $Y$, defined by:

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

where the material density $\theta$, describes the proportion of solid phase in $Y$, which varies from 0 to 1. Note that, because of its rounded corners, the smooth
honeycomb cell can not reach complete void: i.e., the density $\theta$ going to zero is excluded. In this context, for practical reason, this kind of cell is parametrized by another parameter $h \in [0, 1]$, homogeneous to a distance and describes (similar to the parameter $e$) the relative pore diameter of the cell, that is:

$$h = 1 - \frac{2e}{\sqrt{3}}$$

Indeed, in order to design the smooth honeycomb cell $Y$, a parametric curve $\Gamma_h$ (which depends on $h$) is introduced and represents the boundary of the perforated smooth central hole: note that, this only applies in 2-D. We now introduce some notations before giving its polar equation. Let $v(t) = (\cos(t), \sin(t))^T$ and $n_i$ (for $i \in \{0, 1, 2\}$) represent the normal vectors of the three diagonals of $Y$, that are:

$$n_0 = \begin{pmatrix} 0 \\ 1 \end{pmatrix}, \quad n_1 = \begin{pmatrix} \frac{\sqrt{3}}{2} \\ \frac{1}{2} \end{pmatrix}, \quad n_2 = \begin{pmatrix} \frac{\sqrt{3}}{2} \\ -\frac{1}{2} \end{pmatrix}. \qquad (1)$$

The polar equation of the parametric curve $\Gamma_h$ (of the smooth hexagon hole) is defined by:

$$r(t) = h \frac{\sqrt{3}}{2} \left( \sum_{i=0}^{2} |v(t) \cdot n_i|^{k(h)} \right)^{\frac{-1}{k(h)}} \quad \text{with } t \in [0, 2\pi], \qquad (2)$$

where $k(h)$ is positive coefficient, which depends on $h$: for this work, we took $k(h) = 4+20h^2$. In the following, we give some remarks over the polar equation of the parametric curve $\Gamma_h$, which can be extended to other polygons. Let $\mathcal{H}$ be a regular unit hexagon, namely, the set of all points such that the maximal distance of a point in $\mathcal{H}$ from the three diagonals is equal to $\frac{\sqrt{3}}{2}$ and, let $M(r, t)$ be a point, with its polar coordinate denoted $(r, t)$. Thus, $M$ is a point in $\mathcal{H}$ if and only if, its polar coordinate $(r, t)$ satisfies

$$r \max_i |v(t) \cdot n_i| = \frac{\sqrt{3}}{2}. \qquad (3)$$

Hence, the polar equation of $\mathcal{H}$ verifies

$$r(t) = \frac{\sqrt{3}}{2} \left( \max_i |v(t) \cdot n_i| \right)^{-1}. \qquad (4)$$

And, we recall that,

$$\left( \sum_{i=0}^{2} |v(t) \cdot n_i|^k \right)^{\frac{1}{k}} \rightarrow_{k \rightarrow \infty} \max_i |v(t) \cdot n_i|. \qquad (5)$$
Note that, the polar equation of $Γ_h$ comes from combining the polar equation of $H$ and the above limit; the parameter $h$ is added in order to adjust the diameter of its inner hole.

### 2.1.2 Tetrakaidecahedron cell

In 3-D, isotropic cells are not so easy to design. Here, we design a regular tetrakaidecahedron, known as the Kelvin foam, similar to the one introduced in [9, 34]. 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 3 and Fig 4. Several approaches are possible in order to design isotropic cells from this one, parametrized by one parameter, the local density. A naive approach consists in adding a central tetrakaidecahedron inclusion characterized by its size in [0, 1].

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

<figure><img src="image_6.png" /><figcaption>Fig. 4 Tetrakaidecahedron (or Kelvin foam) for different values of the material density $θ$.</figcaption></figure>

Based on the construction method described above, the material density $θ$ has been used as a control parameter to generate foams of chosen porosities. Upon increasing $θ$, one would reach a point where the structure is full of material (i.e., the central tetrakaidecahedron inclusion will be closed).
## 2.2 Cell problem and homogenized elasticity tensor

Here, we only give a few important results on the theory of homogenization. The interested reader will find more details in [1].

Assume that, in a given macroscopic domain $\Omega$, there is a periodic distribution of holes inside an isotropic elastic material, with constant elastic tensor $A$. Let $\epsilon > 0$ be the periodic size and, let $Y$ be the unit periodic pattern: either the regular smooth hexagonal in 2-D or the regular tetrakaidecahedron in 3-D. The periodicity of $Y$ is defined by the same displacement on two opposite and parallel faces; hence, for the hexagonal cell in 2-D, there are three directions of periodicity and for the tetrakaidecahedron cell in 3-D, there are seven directions of periodicity: see Fig. 1 and Fig. 3. Let $Y_0$ be the solid part in $Y$ and, we denote by $|Y|$, the volume of the periodic cell $Y$; let $\Gamma$ be the boundary of the holes (i.e., the solid/void interface) and $n$ be the normal vector to the boundary $\Gamma$.

In addition, assume that whenever $\epsilon$ tends to zero, the porous medium can be considered homogeneous, with effective tensor $A^*(x)$ (defined at each point $x \in \Omega$ of the domain). To compute the homogenized tensor $A^*$, one needs the so-called correctors $w_{ij}$, corresponding to the local displacements in the periodic cell $Y$, defined for each pair $(i, j) \in \{1, \dots, N\}$ ($N \in \{2, 3\}$) as the solutions to the following set of equations:

$$ \left\{ \begin{array}{ll} -\mathrm{div}(A(e_{ij} + e(w_{ij}))) = 0 & \text{in } Y \\ A(e_{ij} + e(w_{ij})) \cdot n = 0 & \text{on } \Gamma \\ y \mapsto w_{ij}(y) & \text{Y-periodic,} \end{array} \right. \qquad (6) $$

where $e_{ij} = \frac{1}{2}(e_i \otimes e_j + e_j \otimes e_i)$ represents the basis of the symmetric tensors of order 2. As a consequence, the variational formulation associated to (6) is: 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 (7) $$

which admits a unique solution (up to a rigid displacement field). The entries of $A^*$ are then defined in terms of the correctors $w_{ij}$ (solutions of (6)):

$$ 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 (8) $$

We emphasize that, if the homogenized tensor $A^*$ is supposed to be isotropic, only two of its coefficients are needed (e.g., $A_{1122}^*$ and $A_{1212}^*$) to fully characterize $A^*$. However, all the coefficients were numerically computed to verify that the homogenized material is isotropic (or quasi-isotropic). Herein, we use a linear material model with Young's modulus $E = 12 \times 10^9 Nm^{-2}$ (i.e., 12GPa) and Poisson's ratio $\nu = 0.35$. It is worth noting that, if the effective tensor $A^*$
is isotropic, it can be written as:

$$A^* = 2\mu^* I_{2N} + (\kappa^* - \frac{2\mu^*}{N}) I_N \otimes I_N, \quad (9)$$

where $\kappa^*$ and $\mu^*$ are the bulk and shear moduli of the homogenized Hooke's law $A^*$, with its first Lamé coefficient defined by: $\lambda^* = \kappa^* - \frac{2\mu^*}{N}$. $I_{2N}$ and $I_N$ in (9), represent the fourth order symmetric identity and the identity tensor of order 2. The entries are defined by:

$$\left\{
\begin{array}{l}
\mu^* = A_{ijij}^* \\
\lambda^* = A_{iijj}^* \\
\kappa^* = A_{iijj}^* + \frac{2}{N} A_{ijij}^*
\end{array}
\right.
\quad \forall i, j \in \{1, \dots, N\} \text{ and } i \neq j \quad (10)$$

where 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_{iijk}^* = 0; \quad \text{if } j \neq k \\
A_{iiii}^* = A_{jjjj}^* \\
A_{iijj}^* = A_{kkll}^*; \quad \text{if } i \neq j, k \neq l \\
A_{iiii}^* = A_{ijij}^* + A_{iijj}^*; \quad \text{if } i \neq j
\end{array}
\right.
\quad (11)$$

## Numerical results

The homogenized tensor $A^*(\theta)$ has been computed for both periodic cells, for discrete values $(\theta_i)_{i=1,...,n_s}$ of $\theta$: a preprocessing stage, which is performed offline. Here, a table of the size $n_s = 1000$ is built, which is then used to compute the local composites during the optimization process: a linear interpolation approach is used to update the homogenized tensor.

The relative errors with respect to the equalities in (11) are depicted by Fig. 5. We note that the residuals with respect to the equalities in (11), for the smooth honeycomb and the Kelvin foam, are sufficiently small, which validates the isotropy assumption: the relative errors are less than $10^{-4}$.

## 2.3 Sensitivity of the homogenized elasticity tensor

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

**Definition 1** Let $\delta \in W_{\#}^{1,\infty}(Y; \mathbb{R}^N)$. The shape derivative of a function $F(Y)$ is defined as the Fréchet derivative in $W_{\#}^{1,\infty}$ at 0 of the application $\delta \to F((Id + \delta)Y)$,
<figure><img src="image_7.png" /><figcaption>Fig. 5 The relative errors with respect to the equalities in (11), for the smooth honeycomb cell in 2-D and the Kelvin foam cell in 3-D.</figcaption></figure>

that is to say:

$$F((Id + \delta)Y) = F(Y) + \langle F'(Y), \delta \rangle + o(\delta) \quad \text{with} \quad \lim_{\delta \to 0} \frac{|o(\delta)|}{\|\delta\|_{W^{1,\infty}}}$$

where $F'(Y)$ is a continuous linear form on $W^{1,\infty}_#(Y; \mathbb{R}^N)$. Now, let $M_N$ be the set of squared $N \times N$ matrices and $M_N^s$ the subset of symmetric ones. Then
**Lemma 1** The shape derivative of $A_{ijkl}^*$ is given by:

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

where $\operatorname{de}(w)$ is a linear operator from $W_{\#}^{1,\infty}(Y; \mathbb{R}^N)$ to $L^2(Y; \mathcal{M}_N^s)$ defined for every $w \in H_{\#}^1(Y; \mathbb{R}^N)$ by

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

See [1, 35], for a proof of Lemma 1: it relies on the Lagrangian method, which amounts to introduce a Lagrangian, defined as the sum of formula (8) for $A_{ijkl}^*$ and of the variational formulation (7). 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 (12). Differentiating with respect to the shape leads to the final result.

## 2.4 Cell problem and homogenized model in conductivity

We now consider a general model problem of thermal or pressure conductivity, similar to the elasticity system defined in Sect. 2.2. Starting from a microscopic description of a problem, one seeks a macroscopic or effective model problem in conductivity $K^*$ which can either be the permeability of a cold fluid flowing through a porous material or the conductivity of a hot fluid flowing through a two-phase material composed of solid and fluid phases. We introduce the so-called cell problems, similar to the elasticity system. Howbeit, since the considered cell $Y$ is isotropic, only one of its coefficient (e.g. $(K^*)_{11}$) could be computed in order to fully characterize $K^*$.

Let $(e_i)_{i=1,...,N}$ be the canonical basis of $\mathbb{R}^N$ and, for each unit vector $e_i$, we consider the following conductivity problem in the periodic cell $Y$:

$$ \left\{ \begin{array}{ll} -\operatorname{div}(K(e_i + \nabla w_i))) = 0 & \text{in } Y \\ y \mapsto w_i(y) & \text{Y periodic,} \end{array} \right. \quad (13) $$

where $w_i(y)$ is the local variation of pressure created by an averaged (or macroscopic) gradient $e_i$. The homogenized conductivity $K^*$ is then given in terms of the correctors $w_i$ (solutions of (13)), defined by:

$$ (K^*)_{ij} = \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\} \quad (14) $$

where the tensor $K^*$ describes the effective or homogenized properties of the heterogeneous microstructure. We recall that, $K^*$ does not depend on the
choice of domain $\Omega$, source term, or boundary conditions and, $Y = Y_0 \cup (Y \setminus Y_0)$ is a disjoint reunion of the solid and void phases, where in thermal conductivity, $K$ can be 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 (15)$$

where, $k_s$ and $k_f$ are some fixed constant thermal conductivities inside the solid and void phases. Though, in pressure conductivity, $K$ is:

$$K(y) = \begin{cases} \epsilon_0 & \text{in } y \in Y_0 \\ 1 & \text{in } y \in Y \setminus Y_0 \end{cases} \qquad (16)$$

where $\epsilon_0$ is a given small positive value to avoid degeneracy. Note that, in all generality, $K^*$ is a tensor of order 2; however, because the chosen cell is isotropic, we have: $\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}^*$.

## Numerical results

The homogenized conductivity $K^*$ has been computed for the smooth honeycomb in 2-D and tetrakaidecahedron in 3-D, for discrete values of $\theta$. On Fig. 6, we plot (with respect to the smooth honeycomb and Kelvin foam) the homogenized conductivity $\frac{\alpha^*}{||\alpha^*||_\infty}$ for discrete values $(\theta_i)_{i=1,...,n_s}$, in comparison to a given smooth function:

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

where $\epsilon_0$ and $K_\infty$ are some fixed constants: $\epsilon_0 = 10^{-4}$, $K_\infty = 10^3$. Note that, $||\alpha^*||_\infty = \sup_{\theta_i} \alpha^*(\theta_i)$. As expected, $\frac{\alpha^*}{||\alpha^*||_\infty}$ is a decreasing function with respect to the density $\theta$. We emphasize that, the residuals $|K_{11}^* - K_{22}^*|/K_{11}^* \le 10^{-6}$ and $K_{12}^*/K_{11}^* \le 10^{-3}$ are sufficiently small, which validates the isotropy assumption. In addition, we note that the homogenized conductivity coefficient $\frac{\alpha^*}{||\alpha^*||_\infty}$ can be approximated by the smooth function $K_D^*$.

## 3 Setting of the three-physic problem

Here, we investigate the weakly coupled model of heat propagation, fluid flow and structure deformation. First, the fluid flow is described using the Biot-Darcy's law, second, the heat propagation is characterized using the convection-diffusion equation and third, the three-physic problem is described by the linearized poro-thermo elasticity system, for the mechanical displacement.
<figure><img src="image_8.png" /><figcaption>Fig. 6 The homogenized conductivity $\frac{\alpha^*}{||\alpha^*||_\infty}$ as a function of $\theta$ for the smooth honeycomb (left) and Kelvin foam (right), compared to $K_D^*$.</figcaption></figure>

Let $\Omega$ be 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 one parameter $\theta$ (the material density), which is to be optimized. Let $n$ be the normal vector to the boundary $\partial\Omega$, pointing outward to the domain $\Omega$. The domain $\Omega$ is described by the three physics system which are governed by three coupled models, that are:

* the motion of the fluid inside the domain $\Omega$ described by the pressure field $p$, satisfying the Biot-Darcy's law

* the diffusion of heat inside the whole domain $\Omega$, and its transport by convection in the fluid domain, resulting in a temperature field $T$;

* the deformation of the solid phase, as a result of the stress exerted by the fluid part and of the dilation induced by thermoelastic effects, characterized by a mechanical displacement $u$

The physical equations chosen for the modeling of the state variables $p$, $T$ and $u$ with their relevant set of boundary conditions are now described in strong form in Section 3.1, Section 3.2 and Section 3.3.

## 3.1 Hydraulic law of Biot-Darcy-type for the pressure variable

We give a precise account of the weakly coupled fluid model that is based on the Biot-Darcy's law; the effect of which is to establish the pressure field as a function of material density $\theta$.

From a fluid point of view, Darcy's law describes the fluid ability to flow through a porous media such as soil, sandstone or rock; it states that, the fluid flowing through a unit area is directly proportional to the pressure drop
per unit length $\nabla p$ and, inversely, that the resistance of the porous medium is proportional to the flow $\mu_f$ ([37]), which is defined by:

$$ q := -\frac{\kappa_f}{\mu_f} \nabla p = -K_D^* \nabla p, \qquad (18) $$

where $q$, $\kappa_f$, $\mu_f$, $\nabla p$, and $K_D^*$ represent the flux ($kgm^2s^{-1}$), permeability ($m^2$), fluid viscosity ($Nm^{-2}s$), pressure gradient ($Nm^{-3}$), and the conductivity coefficient ($m^4N^{-1}s^{-1}$) (which defines the fluid ability to flow through a porous medium). In order to smoothly and continuously distribute the pressure drop in the fluid domain and differentiate between solid phase and void in the whole domain, the homogenized conductivity coefficient $K^*$ (which we numerically compute, see Section 2.4) is approximated by a smooth function:

$$ K_D^*(\theta) := \min \left( \frac{\epsilon_0 + (1 - \epsilon_0)(1 - \theta)}{\theta(x)}, K_\infty \right), \qquad (19) $$

where $\epsilon_0$ and $K_\infty$ are given positive values: $\epsilon_0 = 10^{-4}$, $K_\infty = 10^3$. We recall that an homogenization method was performed on the conductivity coefficient and we notice that, the homogenized flow coefficient $K^*$ can be approximated by the above function; see Fig. 6.

We now assume that the pressure field $p$ is satisfying (a more accurate description of) the Biot's law:

$$ p := Mm - Mbe_{vol}, \qquad (20) $$

where $m, M, b$ and $e_{vol}$ represent the mass of the fluid flow, the so called Biot modulus, Biot coefficient and the volume variation of the solid phase, defined as follows:

$$ m(\theta) := (1 - \theta)\rho, \quad M(\theta) := \frac{1 - \theta}{\kappa_v} - \frac{b(\theta) - (1 - \theta)}{\kappa_s}, \quad e_{vol} := \operatorname{div}(u), \qquad (21) $$

where $\rho$, $\kappa_v$ and $\kappa_s$ represent the density of the flux and the compressibility of the void and solid phases. Moreover, in case of an isotropic porous medium, the Biot's coefficient $b(\theta)$ is explicitly defined as follows:

$$ b(\theta) := 1 - \frac{\kappa_s(\theta)}{\kappa}, \qquad (22) $$

where, $\kappa$ and $\kappa_s(\theta)$ represent the bulk moduli of the elastic tensor $A$ and the homogenized tensor $A^*(x)$ (at each point $x \in \Omega$ of the design domain): note that, $A^*$ tends to $A$, when $\theta$ tends to 1 and thus, $\kappa_s(\theta)$ tends to $\kappa$. In addition, the Biot's law (20) satisfied by the pressure $p$, is assumed to be related to Darcy's law (18):

$$ q := mv_f = -K_D^* \nabla p, \qquad (23) $$
where $v_f$ represents the velocity of the flux. Equation (23), allows to render gradually the pressure drop from the inner pressure boundary to the outer pressure boundary. We emphasize that, this penetrating pressure of Biot-Darcy's law is similar to that introduced in [23], which makes this pressure loading boundary a smeared-out version of an applied pressure load on a sharp boundary. Now, for sake of simplicity, we opted here for a simplified version, that is, in addition to the Biot-Darcy equation (23), we assume that the state equation satisfies the law of conservation of mass (in view of incompressible fluid):

$$ \frac{\partial \mathbf{m}}{\partial t} := -\mathrm{div}(\mathbf{q}) = \mathrm{div}(K_D^* \nabla p) = 0 \qquad (24) $$

And by using the Biot-Darcy's law (24), the weakly coupled fluid model can be defined as follows:

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

where, $f_f$ is an applied Neumann isoflux condition (or mass flow rate) for the pressure variable $p$. The boundary of the fluid phase is the disjoint reunion

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

of a Dirichlet (or inlet) part $\Gamma_D^f$ where the flow enters with a given pressure $p = p_0$, a Neumann (or outlet) part $\Gamma_N^f$ where normal stress (or mass flow rate) $\mathbf{q}_\Gamma \cdot \mathbf{n} = f_f$ is imposed, and free interface $\Gamma^f$ of $\partial\Omega$. At this stage, it is assumed that the deformation of the solid domain is sufficiently small so that no slip boundary conditions hold on: $\mathbf{q}_\Gamma \cdot \mathbf{n} = 0$. Therefore, the variable $p$ depends solely on the material density $\theta(x)$, for all $x \in \Omega$.

## 3.2 Convection-diffusion for the temperature variable

We give a precise account of the weakly coupled diffusion of the heat model, that 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. [38]). Here, we chose to use an empirical form of the heat equation for the whole domain featuring fluid and solid phases (see e.g. [39]). Some authors employ two-equation modeling for each phase of the medium (see e.g [40]). 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.

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 void phases; let $K^*$ denote the homogenized thermal conductivity inside the porous media (which we numerically computed, see Section 2.4), 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:

$$ egin{cases} -div(K^*\nabla T) + \phi(\rho c_p)_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{cases} \quad (26) $$

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 flux $h_0$ is applied and a free interface $\Gamma^T$. Note that, $\phi = 1 - \theta$ is the volume fraction of fluid: i.e., the proportion of the void phase at each point $x \in \Omega$. The convection-diffusion equation (26) is now recast as thermal Biot-Darcy equation:

$$ egin{cases} -div(K^*\nabla T) + (-c_p K_D^* \nabla p)_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{cases} \quad (27) $$

where $\phi \rho v_f = mv_f = -K_D^* \nabla p$ is the fluid flow, satisfying the Biot-Darcy model (25). As a consequence, the convection-diffusion model is defined by (27).

Similar to the Biot-Darcy model (25), 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 void inclusions. As seen in Section 2, to compute the homogenized conductivity $K^*$, we introduce the cell problems. 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 Elasticity with fluid-structure interaction

Finally, the pressure variable $p$ and the temperature $T$ together determine the displacement $u$ of the solid phase in $\Omega$, which we assume to be isotropic
thermoelastic composite material, with homogenized Lamé coefficients denoted $\lambda^*$, $\mu^*$. Let $\alpha$ be the thermal expansion parameter and $T_{ref}$ be the temperature at rest; then, the weakly coupled thermal fluid-elastic model is defined by the linear thermo-elasticity (of Biot-Coussy type):

$$ \left\{ \begin{array}{ll} -\mathrm{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^s, \end{array} \right. \qquad (28) $$

where, the homogenized stress tensor is defined as follows:

$$ \left\{ \begin{array}{ll} \sigma(u, T) = A^*(e(u) - \alpha(T_s - T_{ref})I_N) & \text{and where,} \\ A^*e(u) = 2\mu^*e(u) + \lambda^*\mathrm{Tr}(e(u))I_N & \end{array} \right. \qquad (29) $$

$I_N$ is the identity tensor. Note that, the source term in (28) is the body force (i.e., $-b\nabla p$) exerted by the fluid part. The boundary $\partial\Omega$ is split into 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$. We recall that the Biot coefficient, defined by Eq. (22) is only valid in the case of a linear isotropic elastic tensor microscopically homogeneous. Moreover, we emphasize that the above model (28) is a simplified version of a genuine thermal fluid-solid coupling:

$$ \text{(Biot-Coussy)} \quad \left\{ \begin{array}{ll} -\mathrm{div}(\tilde{\sigma}) = 0 & \text{in } \Omega, \\ \tilde{\sigma} = \sigma(u, T) - bpI, \\ \sigma(u, T) = A^*(e(u) - \alpha(T - T_{ref})I), \\ p = M\mathbf{m} - Mbe_{vol} \end{array} \right. $$

where the fluid pressure $p$ is satisfying a more accurate description of fluid-structure interaction (which features a transition regime and an inertia regime):

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

where $C_{for}$ is an inertia parameter of the fluid flow, called Forchheimer coefficient. Here, for sake of simplicity, we opted for a simplified version, which is justifiable insofar as we wish to obtain a first qualitative result of microstructure without however sizing as accurately as possible the system. Thanks to this simplification, the systems (25), (27) and (28) are only weakly coupled: the resolution is achieved by solving the fluid system (25), then using the fluid stress resulting from the pressure $p$ in the heat transfer equation (27), and finally using the fluid stress and the temperature $T$ to solve the linear poro-thermo elasticity system (28).
*Remark 1:* Regarding the domain of validity, several characteristics length scales have been proposed in the literature, i.e., various definitions in the form of mor- phological 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 [38], 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 simulation, 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

with respect to circular strut cross-section having variable pore diameters and poros- ity for a given cell size: constant fluid properties were used. As aforementioned, we opted here for simplified regime, the Darcy regime (which is justifiable insofar as we wish to obtain a first qualitative result of microstructure without however sizing as accurately as possible the system). In addition, we opted for isotropic compos- ite microstructures (i.e., composite periodically perforated by hexagonal cells in 2-D and tetrakaidecahedron cells in 3-D) and as such, the thermal or hydraulic parame- ter (e.g. thermal or pressure conductivity) can be identified to scalar values, which depend solely on the material density $\theta$. Furthermore, regarding the regimes of strong couplings with respect to those of weak ones, it is worth to recall that, from a mod- eling 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, indepen- dent 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. 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).

## 4 The optimization problem formulation and its sensitivity analysis

We present the optimization problem formulation associated to the weakly cou- pled three-physic problem and discuss the sensitivity analysis for such design problems. The announced goal is the resolution of the relaxed constrained
version of the original optimization problem:

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

where $\theta$ denotes the material density of the structure to be optimized. $J^*$ refers to a given objective function which quantifies the system performance and which is to be minimized. Classical objective functions commonly encountered in shape optimization include the compliance of a mechanical structure, the drag force induced by an airfoil, or the heat stored into a thermal system. Equality and inequality constraints $g_i$ and $h_j$ model physical load specifications; they refer to target values some physical quantities 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). Both objective function $J^*$ and constraints $g_i, h_j$ depend on the optimization variable $\theta$ and on the physical variables $p(\theta), T(\theta), u(\theta)$, which depend themselves on $\theta$ through physical state equations.

For its resolution, we shall rely on the alternate directions algorithm [26] detailed in Section 5.1, which requires to compute the sensitivity of the above functionals: here, the adjoint-state method [1] is used to determine the sensitivities of the objective functions and constraints with respect to the design variable $\theta$. It is worth mentioning that, likewise in [25], the authors employed the adjoint-state method to evaluate the sensitivities of objective functions (constituted by the Biot model coefficients with respect to the underlying pore shape described by a B-spline box which embeds the whole representative cell) and the gradient-based method to solve the optimization problems. Howbeit, the considered shape derivatives of the homogenized coefficients are different: they rely on the shape sensitivity technique and the material derivative approach. Here, the "shape" is described by the density $\theta$ and, the shape derivatives are performed using the Lagrangian method and the projected gradient algorithm. In this context, 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 $\mathcal{U}_{ad}$, defined as follows:

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

In contrast, in [25], a general optimization variable $\alpha$ is introduced which is related to the effective medium parameters: it determines the homogenized coefficients for any position $x \in \Omega$.
## 4.1 A fully Lagrangian setting for the sensitivity analysis

We use a fully Lagrangian setting to compute the derivative of general objective functionals in the simplified setting of Section 3.1 to Section 3.3. The sensitivities of the state variables $p(\theta)$, $T(\theta)$ and $u(\theta)$ are calculated first, in order to obtain the derivative of an arbitrary objective functional in volume form. Then, under sufficient regularity assumptions, the adjoint method together with suitable augmented Lagrangian functional yield derivative formulas.

## 4.2 A modified objective functional and Lagrangian derivative

In a gradient-based topology optimization, it is essential to determine sensitivities of the objective functional and the constraints with respect to the design variable(s). The starting remark is that the relaxed functional $J^*$, although appearing naturally in the formulation of the optimization problem (30) is not so convenient for the mathematical analysis. Indeed, the domain of definition of $J^*(\theta, \dots, \theta)$ is a functional space which depends on the first argument $\theta$. In order to discuss the precise mathematical settings of this three-physic problem, we introduce the functional spaces which are required, that are:

$$ 
\begin{aligned}
V(\Gamma_D^f) &= \{q \in H^1(\Omega) \mid q = 0 \text{ on } \Gamma_D^f\}, \text{ for the pressure variable } p \\
V(\Gamma_D^T) &= \{S \in H^1(\Omega) \mid S = 0 \text{ on } \Gamma_D^T\}, \text{ for temperature variable } 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 (31) $$

In addition, note that the non-homogeneous Dirichlet boundary data $u_0, p_0$ and $T_0$ featured in (25) to (28) 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}{ll}
\displaystyle \int_{\Omega} K_D^* \nabla p \cdot \nabla q \, dx + \int_{\Gamma_N^f} f_f q \, ds = 0, & \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, & \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, & \forall v \in V(\Gamma_D^s)
\end{array}
\right.
\quad (32) $$

In order to address the sensitivity of an arbitrary objective function, the classical idea is to work within a Lagrangian framework. Therefore, we consider
the corresponding Lagrangian:

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

where again, we point out that $J^*$ is assumed to be a "smooth" enough function, otherwise we cannot apply the adjoint-state method [1, 11].

## 4.3 Sensitivity analysis for the three-physic problem

The sensitivities are evaluated by differentiating the Lagrangian (33) with respect to state variables $u$, $p$ and $T$ in directions $\phi_u \in H^1(\Omega)^N$, $\phi_p \in H^1(\Omega)$ and $\phi_T \in H^1(\Omega)$. Let $(u, \underline{u}) \in H_0^1(\Omega; \mathbb{R}^N)^2$ be a stationary point of $\mathcal{L}$. Then, the derivative of the Lagrangian (33) with respect to $u$, in direction $\phi_u \in H^1(\Omega)^N$ is given 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 \frac{\partial \sigma(\hat{u}, \hat{T})}{\partial u}, \phi_u \right\rangle : e(\hat{\underline{u}}) \right) dx, \quad (34) $$

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

$$
\begin{equation}
\begin{split}
\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 \hat{\underline{u}} \, dx + \int_{\Omega} (-c_p \underline{T} K_D^* \nabla \phi_p) \cdot \nabla \hat{T} \\
& + \int_{\Omega} K_D^* \nabla \phi_p \cdot \nabla \underline{p} \, dx
\end{split}
\tag{35}
\end{equation}
$$

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

$$
\begin{equation}
\begin{split}
\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(\hat{\underline{u}}) \, dx \\
& + \int_{\Omega} K^* \nabla(\phi_T) \cdot \nabla \underline{T} \, dx + \int_{\Omega} (-c_p \underline{T} K_D^* \nabla \hat{p}) \cdot \nabla (\phi_T) \, dx,
\end{split}
\tag{36}
\end{equation}
$$

which, when equations (34), (35) and (36) vanish, are the variational formulation of the adjoint-states. Moreover, the derivatives with respect to $\underline{u}$, $\underline{p}$ and $\underline{T}$, in directions $\phi_u \in H^1(\Omega)^2$, $\phi_p \in H^1(\Omega)$ and $\phi_T \in H^1(\Omega)$ are simply 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. 
\quad (37) $$

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

$$ \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( \begin{aligned} & \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. + \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} \end{aligned} \right) \, dx \quad (38) $$

Note that Eq. (38) is defined using the adjoint-state method. The term $\frac{\partial J^*}{\partial \bar{\theta}}$, $\bar{\theta} > 0$ is the partial derivative of the objective function $J^*$ with respect to $\theta$, in direction $\bar{\theta}$, while the term

$$ \int_{\Omega} \left( \begin{aligned} & \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} \end{aligned} \right) \, dx $$

is used to determine the load sensitivities of the constraints with respect to $\theta$. Here, $\underline{u}$, $\underline{p}$ and $\underline{T}$ are the adjoint-state variables, solutions to the adjoint equations (34)-(36) (which we have to solve first). $\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_s}$ of $\theta$, using formula (12) (a preprocessing stage, which is performed offline) and, a linear interpolation approach is used to update the derivative with respect to $\theta$ of $\sigma$: $(\sigma(\theta))' = (A^*(\theta))'$: $e(u)$. As for the derivative with respect to $\theta$ of $K^*(\theta)$, $K_D^*(\theta)$ and the Biot coefficient $b(\theta)$, they are approximated by the derivative of smooth "enough" functions; see equations (17) and (22).
# 5 Topology optimization of modulated periodic composite materials

Here, our numerical algorithm is proposed; we describe how the methodology applies to the weakly coupled three physics system of (25), (27) and (28).

## 5.1 Alternate directions method

The problem (30) is solved using the alternate directions algorithm [26].

### 5.1.1 Minimizing over the stress field

For a given design field $\theta$, the minimization with respect to the stress field $\sigma$ amounts to solve the linear poro-thermo elasticity problem (28), with a material of elasticity tensor equal to $A^*$. The composite design $(\theta, A^*(\theta))$ is computed by linear interpolation in the table obtained after the preprocessing stage, performed offline; see Section 2.2.

### 5.1.2 Minimizing over the density field

For a given stress field $\sigma$, the minimization with respect to the density $\theta$ is performed using the projected gradient algorithm. Since the problem (30) is not self-adjoint, one needs to define the associated adjoint problem. As a consequence, the descent direction is given by the derivative of $\mathcal{L}$ with respect to $\theta$:

$$ \left\langle \frac{\partial \mathcal{L}}{\partial \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}) + \left( \frac{\partial b}{\partial \theta} \nabla p \cdot \underline{u} \right) + \frac{\partial K_D^*}{\partial \theta} \nabla p \cdot \nabla \underline{p} \right) dx \\ + \frac{\partial K^*}{\partial \theta} \nabla T \cdot \nabla \underline{T} + \left( -c_p \underline{T} \frac{\partial K_D^*}{\partial \theta} \nabla p \right) \cdot \nabla T + \ell \right) \bar{\theta} \, dx, \quad (39) $$

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. \quad (40) $$

It 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} \right) dx + \frac{\partial K^*}{\partial \theta} \nabla T \cdot \nabla \underline{T} \\ + \left( -c_p \underline{T} \frac{\partial K_D^*}{\partial \theta} \nabla p \right) \cdot \nabla T + \ell \quad (41) $$
with $< \frac{\partial J^*}{\partial \bar{\theta}}, \bar{\theta} > = \int_{\Omega} \frac{\partial J^*}{\partial \bar{\theta}} \cdot \bar{\theta} \, dx$. At iteration $n$, the optimal density $\theta$ is then updated by performing the projected gradient:

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

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

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

where $\eta$ is a small coefficient, which typically depends on the size of the elements of the mesh: thanks to this coefficient, 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, structured as follows:

1. Initialization of the design variable $\theta$ 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$:

(a) Computation of the state variable $p^n$ through the Biot-Darcy model (25), with composite design $\{\theta^n, A^*(\theta^n)\}$

(b) Computation of the state variable $T^n$ through the convection-diffusion equations (27), with composite design $\{\theta^n, A^*(\theta^n)\}$

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

(d) Updating the design variable $\theta^{n+1}$ using formulas (42), for the descent direction $d\theta^n$ and then updating the composite design $\{\theta^{n+1}, A^*(\theta^{n+1})\}$, by linear interpolation.

## 6 Numerical results and discussion

We introduce our numerical results in the three physics context detailed in Section 3. A variety of 2-D and 3-D test cases are presented to demonstrate
that our alternate directions algorithm produces physically correct results. The algorithm (5.1.3) has been implemented in FreeFem++[41], where all the unknowns are discretized using $P_1$ finite elements. For all our computations, a linear material model with Young's modulus $E = 12 \times 10^9 \text{Nm}^{-2}$ (12GPa) and Poisson's ratio $\nu = 0.35$ are considered. The void (i.e., $\theta = 0$) is replaced with a very compliant material: 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.

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

## 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 fluid is entering with a given mass flow rate $f_f$ on its Neumann (or outlet) part $\Gamma_N^f$ (of length $L/5$), while submitted to a pressure load $p_0$ on its Dirichlet (or inlet) part $\Gamma_D^f$ (of length $L/5$). On the free interface $\Gamma^f$, zero normal stress is applied. See Fig. 7 for a sketch of this test case.

<figure><img src="image_9.png" /><figcaption>**Fig. 7** Setting for the fluid-structure interaction problem of Section 6.1. The boundary parts $\Gamma_D^f$, $\Gamma_N^f$ and $\Gamma_N^s$ are centered.</figcaption></figure>

We neglect the thermal effects (namely, Eq. (27) is ignored), so that Eq. (28), 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
<table><caption>**Table 1** Numerical values of the physical parameters for test case of Section 6.1</caption><thead><tr><th>L (m)</th><th>f<sub>f</sub> (kg/s)</th><th>p<sub>0</sub> (bar)</th></tr></thead><tbody><tr><td>0.1</td><td>1.5</td><td>0</td></tr></tbody></table>

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 \underbrace{\left( \int_{\Omega} A^* e(u) : e(u) \, dx \right)}_{\text{Elastic strain energy}} + (1 - \omega) \underbrace{\left( - \int_{\Omega} K^* \nabla p \cdot \nabla p \, dx \right)}_{\text{Hydraulic power}} \quad (44)$$

for some fixed coefficient $\omega \in [0, 1]$: it measures the relative weight given to each term in (44). The objective function $J^*$ is minimized using the alternate directions algorithm of Section 5.1 and subject to a volume constraint $\Theta = 44\%$ of volume $|\Omega|$:

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

In this example, the numerical values of the considered physical parameters are supplied by Table 1. Here, we consider two cases in (44): either $\omega = 1$ or $\omega = 2/5$. Note that, unfortunately the fixed coefficient $\omega$ is not optimal and so, the terms composing the objective (44) 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). On Fig. 8, the optimal density is displayed, while on Fig. 9, the corresponding pressure field and solid displacement are displayed.

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

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 (with microstructures periodically perforated by smooth honeycomb cells). 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 zone at the center. 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. 9. On Fig. 10, we plot the convergence history for this calculation: with $\omega = 1$, 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.

<figure><img src="image_11.png" /><figcaption>**Fig. 9** (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>

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

This second example is similar to that introduced and solved in [22], 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 (or outlet) part $\Gamma_N^f$ (of length $L/5$) and, while submitted to a pressure load $p_0$ on its Dirichlet (or inlet) part $\Gamma_D^f$ (of length $L/5$). In addition, a temperature $T_1$ is entering through 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. 11 for a sketch of this test case.

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

For this example, we neglect the elastic contribution (namely, Eq. (28) is ignored), so that the example involves only a coupling of the flow (25) and the heat equation (27). 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 \underbrace{\left( - \int_{\Omega} K^* \nabla p \cdot \nabla p \, dx \right)}_{\text{Hydraulic power}} + (1 - \omega) \left( \int_{\Omega} c_p K_D^* \nabla p \cdot \nabla T \, dx \right), \quad (45)$$

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\%$ of volume $|\Omega|$. The numerical values of the considered physical parameters are supplied by Table 2. The balance coefficient is set to $\omega = 1/2$. The optimal density, the corresponding pressure $p$ and temperature $T$ are displayed on
<table><caption>**Table 2** Numerical values of the physical parameters in the convective heat transfer problem of Section 6.2.</caption><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>

Fig. 12 and Fig. 13. The topology of the result is quite similar to that obtained in [22] (where the solution is autopenalized).

<figure><img src="image_14.png" /><figcaption>(a) Density</figcaption></figure>

**Fig. 12** The optimal density for the two-physic topology optimization problem of Section 6.2.

<figure><img src="image_15.png" /><figcaption>**Fig. 13** Pressure $p$ (bar) and temperature $T$ (C) for the optimized configuration of the two-physic problem of Section 6.2.</figcaption></figure>

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. 13. On Fig. 14, we plot the convergence history for this calculation: smooth and relatively rapid convergence is observed.

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

## 6.3 A thermal fluid-structure interaction problem

We finally turn to a topology optimization example in the full three-physic setting presented in Section 3; a test case which is new to the best of our knowledge. The computational domain is a square of dimension $L$, fixed on its boundary $\Gamma_D^s$ (of length $L/10$). A fluid is entering with a given mass flow rate $f_f$ on its Neumann part $\Gamma_N^f$ (of length $L/5$). It is submitted to a pressure load $p_0$ on its Dirichlet part $\Gamma_D^f$ (of length $L/5$). In addition, a temperature $T_1$ is entering through the upper and lower walls (i.e., on its Dirichlet part $\Gamma_D^T$ where $T = T_1$), while a temperature $T_{in}$ is imposed on the middle left part of its boundary $\Gamma_D^T$ (of length $L/5$). On the free interface $\Gamma = \Gamma^f = \Gamma^T$, zero normal stress is imposed for both the temperature and the pressure. See Fig. 15 for a sketch of this test case.
Here, the objective

$$ J^*(\theta, p(\theta), T(\theta), u(\theta)) := \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), \quad (46) $$
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, for some given coefficients $\omega_1, \omega_2, \omega_3 \in [0, 1]$, such that $\omega_1 + \omega_2 + \omega_3 = 1$.

<figure><img src="image_17.png" /><figcaption>**Fig. 15** Setting of the three-physic problem of Section 6.3. Here, the boundary parts $\Gamma_D^f, \Gamma_N^f$ and $\Gamma_D^T$ are centered.</figcaption></figure>

<table><caption>**Table 3** Numerical values of the physical parameters in the three-physics problem of Section 6.3.</caption><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>ref</sub></th><th>T<sub>1</sub></th><th>α</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>0</td><td>10</td><td>1</td><td>10</td><td>1</td></tr></tbody></table>

The objective function $J^*$ is minimized with a volume constraint $\Theta = 23\%$ of volume $|\Omega|$. The thermal capacity of the fluid is set to $c_p = 0.5$ ($J/kg.K$) and the remaining numerical values of the considered physical parameters are supplied by Table 3. We consider two cases: $\omega_1 = \omega_2 = \omega_3 = 1/3$ and $\omega_1 = 1/5$, $\omega_2 = \omega_3 = 2/5$. The optimal density is displayed on Fig. 16. Very interestingly, we retrieve the fact that in the case where $\omega_1 = \omega_2 = \omega_3 = 1/3$, the topology of the solution contains a large composite zone at the center and a large contact surface with the fluid at entrance, namely, the middle right wall, so to reduce the effect of the pressure source; the algorithm tends to distribute more material in regions where the pressure gradient is high and seems to be driven by the elastic strain energy. In the case $\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. On Fig. 17 and Fig. 18, the corresponding pressure $p$, temperature $T$ and displacement $u$ are displayed, for the case $\omega_1 = \omega_2 = \omega_3 = 1/3$. We notice
a cool down of the heat at the center of the optimal shape and becomes hot elsewhere, with high pressure gradient located mainly at large solid regions.

<figure><img src="image_18.png" /><figcaption>Fig. 16 The optimal density for the three-physic problem of Section 6.3, in both cases.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 17 Pressure $p$ and temperature $T$ for the optimized configuration of the three-physic problem of Section 6.3, with $\omega_1 = \omega_2 = \omega_3 = 1/3$.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>**Fig. 18** The solid displacement for the three-physic topology optimization problem of Section 6.3, with $\omega_1 = \omega_2 = \omega_3 = 1/3$ (displacement amplified by a factor 3).</figcaption></figure>

On Fig. 19, we plot the convergence history for the elastic strain energy, hydraulic power and thermal strain energy (for $\omega_1 = \omega_2 = \omega_3 = 1/3$) and, the corresponding objective history (for both cases considered). We note that, the objective $J^*$ (as well as the elastic strain energy) increases sometimes substantially (in the first part of the optimization) due to the fact that the volume constraint is not yet satisfied.

<figure><img src="image_21.png" /><figcaption>**Fig. 19** The convergence history for the three-physic topology optimization problem of Section 6.3: to the right, the elastic strain energy (1st term in $J^*$ (46)), the hydraulic power (2nd term in $J^*$) and thermal strain energy (3rd term in $J^*$), in the case $\omega_1 = \omega_2 = \omega_3 = 1/3$.</figcaption></figure>

## 6.4 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 temperature $T = T_{in}$ is imposed on the middle back and left walls $\Gamma_D^f$ (of size $\frac{L}{5} \times \frac{L}{5}$). On $\Gamma = \Gamma^f = \Gamma^T$, zero normal stress is imposed for both the temperature and the pressure. See Fig. 20 for a sketch of this test case.

<figure><img src="image_22.png" /><figcaption>**Fig. 20** The setting of the 3-D two-physical problem of Section 6.4, subjected to fluid-thermal loads.</figcaption></figure>

We recall that, for this example, the elastic model (28) is ignored. Thus, its
involves only a coupling of the flow (25) and heat equation (27). 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 (45) for the optimization problem. The functional $J^*(45)$, is minimized
under the volume constraint $\Theta = 44\%$ of the volume $|\Omega|$. 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. 21 represents the iso-surface $\theta \ge 0.5$ of the density. In this example, the
iso-surfaces are smooth and embedded into each other as $\theta$ increases. The
topology of solution is quite similar to that obtained in 2-D, namely, 3-D pipe-
like "shape" is observed, although it contains a large composite zone at its
center.
<figure><img src="image_23.png" /><figcaption>**Fig. 21** Two views of the solution of the 3-D thermal fluid interaction problem of Section 6.4. Composite solution represented as the part of the domain where $\theta \ge 0.5$.</figcaption></figure>

## 6.5 3-D thermal fluid-structure interaction problem

This last example is the 3-D analog of third example of Section 6.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. 22 for a sketch of this test case.

<figure><img src="image_24.png" /><figcaption>**Fig. 22** The setting of the 3-D three-physic problem of Section 6.5, subjected to fluid-thermal loads. On the red rectangles, the solid is clamped.</figcaption></figure>
We recall that the objective is to achieve a trade-off between the maximization of the hydraulic power dissipated by the fluid, the maximization of the heat convected by the fluid and minimization of the compliance of the solid body. The objective function $J^*$ (46) is submitted to volume constraint $\Theta = 30\%$ of the volume $|\Omega|$. Like in the 2-D test case (6.3), we consider two different sets of coefficients in Eq. (46): $\omega_1 = \omega_2 = \omega_3 = 1/3$ and $\omega_1 = 1/5$, $\omega_2 = \omega_3 = 2/5$. On Fig. 23, we display the iso-surface $\theta \ge 0.3$ of the density for both cases. On Fig. 24, the corresponding solid displacement is displayed.

<figure><img src="image_25.png" /><figcaption>**Fig. 23** 3-D three-physic problem of Section 6.5: composite solution displayed as part of the domain filled with values of $\theta \ge 0.3$, in both cases.</figcaption></figure>

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

Very interestingly, we retrieve the fact that the topology of the result is quite similar to that obtained in the 2-D case, for both cases; howbeit, in the case
where $\omega_1 = \omega_2 = \omega_3 = 1/3$, the topology of the solution contains a large composite zone at the center 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.

# 7 Conclusions

In this study, we have demonstrated the relevance of shape and topology optimization for generating unconventional design involving two or three-physic interactions using the homogenization method. In our proposed method, Biot-Darcy's law and the convection-diffusion equation are employed to characterize the pressure (of the fluid flow) and the temperature (of the heat flux), the effect of which is weakly coupled to the solid phase by solving the associated PDEs using the standard finite element method. The porosity of each finite element is related to the material density through a smooth enough function to ensure a smooth transition between void and solid phases. The physical parameters (for fluid and heat flow) are numerically computed in the case of isotropic porous medium. The method facilitates calculation of the load sensitivities with respect to the design variables, using the adjoint-state method. It is noticed that consideration of load sensitivities within the approach does alter the composite designs and are particularly important when designing multi-physic systems. In contrast to methods that use explicit boundary tracking, the Biot-Darcy's model, together with our simplified heat transfer equation offer the potential for relatively straightforward extension to 3-D problems. The effectiveness and robustness of the proposed homogenization method is verified by minimizing several objective functionals.
**Acknowledgments.** This research project was fully sponsored by IFPEN: IFP energies nouvelles; and was supported by Pierre Viot: design and simulation engineer at IFPEN.

## Declarations

On behalf of all authors, the corresponding author states that there is no conflict of interest.

## 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.

## References

[1] Allaire, G., Shape Optimization by the Homogenization Method, Springer, Applied Mathematical Sciences, 146 (2002).

[2] Murat, F., Contre-exemples pour divers problèmes où le controle intervient dans les coefficients, Ann. Mat. Pura Appl., 112 (1997) 49-68.

[3] Kohn, R.V. and Strang, G., Optimal design and relaxation of variational problems, I. Comm. Pure Appl. Math, 39(1) (1986) 113-137.

[4] Lurie, K.A. and Cherkaev, A.V. and Fedorov, A.V., Regularization of optimal design problems for bars and plates, I, II. J. Optim. Theory Appl., 37(4) (1982) 499-522, 523-543.

[5] Bendsøe, M.P. and Kikuchi, N., Generating optimal topologies in structural design using a homogenization method, Comput. Methods Appl. Mech. Engrg, 71(2) (1988) 197-224.

[6] Bendsøe, M.P., Optimal shape design as a material distribution problem, Structural and multidisciplinary optimization, 1(4) (1989) 193-202.

[7] Bendsøe, M.P. and Sigmund, O., Topology optimization: theory, methods, and applications, Springer Science & Business Media, (2003).

[8] Groen, J. and Sigmund, O., Homogenization-based topology optimization for high-resolution manufacturable microstructures, Internat. J. Numer. Methods Engrg., 113(8) (2018) 1148-1163.

[9] Geoffroy-Donders, P., Homogenization method for topology optimization of structures built with lattice materials, PhD thesis, Université Paris Saclay (COmUE), 2018.
[10] Geoffroy-Donders, P. and Allaire, G. and Michailidis, O. and Pantz, O.,
Coupled optimization of macroscopic structures and lattice infill, IJNME,
123(13), (2022) pp.2963-2985.

[11] Allaire, G., Conception optimale de structures, Springer Berlin, Heidelberg, (2007).

[12] Céa, J., Conception optimale ou identification de formes, calcul rapide de la dérivée directionnelle de la fonction coût, ESAIM: Mathematical Modelling and Numerical Analysis, 30(6) (1986) 371-402.

[13] Hammer, V.B. and Olhoff, N., Topology optimization of continuum structures subjected to pressure loading, J. Structural and Multidisciplinary Optimization, 19(2) (2000) 85-92.

[14] Du, J. and Olhoff, N., Topological optimization of continuum structures with design-dependent surface loading - Part I: New computational approach for 2D problems, J. Structural and Multidisciplinary Optimization 27(3) (2004) 151-165.

[15] Fuchs, M.B. and Shemesh, NNY., Density-based topological design of structures subjected to water pressure using a parametric loading surface, J. Structural and Multidisciplinary Optimization, 28(1) (2004) 11-19.

[16] Zheng, B. and Chang, C.J. and Gea H.C., Topology optimization with design-dependent pressure loading, J. Structural and Multidisciplinary Optimization, 38(6) (2009) 535-543.

[17] Lee, E. and Martins, JRRA., Structural topology optimization with design-dependent pressure loads., J. Computer Methods in Applied Mechanics and Engineering, 40-48 (2012) 233-236.

[18] Li, Z-m and Yu, J. and Yu, Y. and Xu, L., Topology optimization of pressure structures based on regional contour tracking technology, J. Structural and Multidisciplinary Optimization, 58(2) (2018) 687-700.

[19] Mendes, E. and Sivapuram, R. and Rodriguez, R. and Sampaio, M. and Picelli, R., Topology optimization for stability problems of submerged structures using the TOBS method, Computers & Structures, 259 (2022) 106-685.

[20] Wang, C. and Zhao, M. and Ge, T., Structural topology optimization with design-dependent pressure loads, Structural and Multidisciplinary Optimization, 53(5):1005-1018 (2016).

[21] Picelli, R. and Neofytou, A. and Kim, HA., Topology optimization for design-dependent hydrostatic pressure loading via the level-set method,
Structural and Multidisciplinary Optimization, 60(4):1313-1326 (2019).

[22] Feppon, F. and Allaire, G. and Bordeu, F. and Cortial, J. and Dapogny, C., Shape Optimization of a Coupled Thermal Fluid-Structure Problem in a Level Set Mesh Evolution Framework, J. Boletín de la Sociedad Española de Matemática Aplicada, Springer 76(3) (2019) 413-458.

[23] Kumar, P. and Frouws, J.S. and Langelaar, M., Topology Optimization of Fluidic Pressure Loaded Structures and Compliant Mechanisms using the Darcy Method, Springer: Structural and Multidisciplinary Optimization, 61 (2020) 1637-1655.

[24] Kumar, P., TOPress: a MATLAB implementation for topology optimization of structures subjected to design-dependent pressure loads, Springer: Structural and Multidisciplinary Optimization 66, 97 (2023).

[25] Hübner, D. and Rohan, E. and Lukeš, V. and Stingl, M., Optimization of the porous material described by the Biot model, Elsevier, International Journal of Solids and Structures, 156-157 (2019) 216-233.

[26] Allaire, G. and Bonnetier, E. and Francfort, G. and Jouve, F., Shape optimization by the homogenization method, Numerische Mathematik, 76 (1997) 27-68.

[27] Murat, F., and Simon, J., Sur le contrôle par un domaine géométrique, rapport 189, univ, Paris VI, (1976).

[28] Henrot, A. and Pierre, M., Variation et optimisation de formes: une analyse géométrique, Springer Science & Business Media, (2006).

[29] Murat, F. and Tartar, L., Calcul des variations et homogénéisation, Les lectures de l'homogénéisation : théorie et applications en physique, (Eyrolles, 1985) 319-369.

[30] Pantz, O. and Trabelsi, K., A post-treatment of the homogenization method for shape optimization, SIAM Journal on Control and Optimization, 47(3) (2008) 1380-1398.

[31] 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) 4547.

[32] Abad, K.M.E. and Khanoki, A.S. and Pasini, D.F., Fatigue design of lattice materials via computational mechanics: Application to lattices with smooth transitions in cell geometry, International Journal of Fatigue, 47 (2013) 126-136.
[33] Vigdergauz, S., Energy-minimizing inclusions in a planar elastic structure with macroisotropy, Structural optimization, 17(2-3) (2013) 104-112.

[34] Lin, C.Y. and Kikuchi, N. and Hollister, J., A novel method for biomaterial scaffold internal architecture design to match bone elastic properties with desired porosity, Journal of biomechanics, 43(5) (2004) 623-636.

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

[36] Hashin, Z. and Shtrikman, S., A variational approach to the theory of the elastic behaviour of multiphase materials, Journal of the Mechanics and Physics of Solids, 11(2) (1963) 127-140.

[37] Batchelor, G., An introduction to fluid dynamics, Cambridge university press, 2000.

[38] Jobic, Y. and Kumar, P. and Topin, F. and Occelli, R., Transport properties of solid foams having circular strut cross section using pore scale numerical simulations, Heat Mass Transfer, 54 (2018) 2351-2370.

[39] Vafai, K. and Tien, C.L., Boundary and inertia effects on flow and heat transfer in porous media, International Journal of Heat and Mass Transfer, 24(6) (1981) 195-203.

[40] Wang, P. and Vafai, K. and Liu, D.Y., Analysis of the volumetric phenomenon in porous beds subject to irradiation, Numerical Heat Transfer, Part A: Applications, 70(2) (2016) 567-580.

[41] Hecht, F., J. Math. Anal. Appl., New development in freefem++, J. Numerical Mathematics, (2012).