# Time-averaged algorithm for solving the topology optimization problem for unsteady laminar, turbulent and anisothermal flows

Delphine Ramalingom, Alain Bastide, Pierre-Henri Cocquet, Michaël Rakotobé, David Marti

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

## ► To cite this version:

Delphine Ramalingom, Alain Bastide, Pierre-Henri Cocquet, Michaël Rakotobé, David Marti. Time-averaged algorithm for solving the topology optimization problem for unsteady laminar, turbulent and anisothermal flows. Journal of Computational Physics, 2025, 538, pp.114175. 10.1016/j.jcp.2025.114175. hal-05113168v2

HAL Id: hal-05113168

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

Submitted on 16 Jul 2025

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

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

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

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

# Time-averaged algorithm for solving the topology optimization problem for unsteady laminar, turbulent and anisothermal flows

Delphine Ramalingom ID a,*, Alain Bastide ID a, Pierre-Henri Cocquet ID b,
Michaël Rakotobé a, David Marti a

## ARTICLE INFO

**Keywords:**

Unsteady flow

Turbulent flow

Anisothermal flow

Large eddy simulation

Topology optimization

Reynolds-averaged method

Estimated turbulent kinetic viscosity

Estimated turbulent thermal diffusivity

## ABSTRACT

This paper proposes a new algorithm to solve topology optimization problems for laminar unsteady or turbulent flows. Instead of computing the gradient of the cost function after solving the direct and adjoint (both unsteady) PDE on the full time interval, our algorithm uses averaged physical quantities on a smaller unspecified time interval to define a (steady) Reynolds-Averaged Method (RAM) model which is then used as constraint in an optimization problem to update the design variable. Another feature of the proposed method is that the RAM model can be defined whatever the initial model and CFD turbulence models initially chosen to compute the instantaneous physical quantities. The RAM model involves turbulent quantities such as turbulent kinetic viscosity and turbulent thermal diffusivity are estimated instead of using the concept of "frozen turbulence". In contrast with the classical methods built to solve unsteady topology optimization problems, the main advantage of the proposed algorithm is that it updates the design variable by solving an auxiliary steady topology optimization problem. Three configuration cases are studied to illustrate the ability of our algorithm to optimize pressure losses and heat transfer by adding material to smooth the laminar unsteady or turbulent flows. We also calculate the number of required design parameter updates to obtain an optimized design. Thus, our algorithm overcomes three major scientific challenges in solving optimization problems in turbulence, namely leveraging efficient temporal turbulence models or a Direct Numerical Simulation (DNS) model, computational cost and data storage requirements.

## 1. Introduction

Topology optimization for fluid flow problems aims at finding the location and physical parameters of a solid inside a fluid to minimize a given physical criterion (called the cost or objective function). Thus, it involves solving an optimization problem, while satisfying a set of design constraints and governing equations that define the physical behaviour of the system. This process can be defined as an iterative, gradient-based design method. Material distribution within the domain is represented by a density field, which assumes a value zero in the fluid regions and value one in the solid regions. To apply gradient-based optimization techniques, the discrete zero-one problem is relaxed, allowing the density field to take intermediate values. Brinkman penalization is typically utilized to define a variable porous material corresponding to these intermediate density values [1] (see also e.g. [2–10] for other

<footnote>ª Physique et Ingénierie Mathématique pour l'Énergie et l'Environnement (PIMENT), Université de la Réunion, 2 rue Joseph Wetzellr, 97490, Sainte-Clotilde country,</footnote>

<footnote>ᵇ Laboratoire des Sciences de l'Ingénieur Appliquées à la Mécanique et au Génie Electrique (SIAME), E2S-UPPA, Université de Pau et des Pays de l'Adour, 64000, Pau, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: delphine.ramalingom@univ-reunion.fr (D. Ramalingom).</footnote>
works using the Brinkmann penalization). While not all fluid flows reach a steady state, the majority of studies have focused on steady-state flow problems (see the review [11, Section 2.1.2]). However, there are situations where the unsteady motion can be exploited by the design, such as transients in onset flows or vortex shedding obstacles. Turbulent flow can be simulated through various numerical approaches [12], including Direct Numerical Simulations (DNS), Large Eddy Simulations (LES), and the Reynolds-averaged Navier-Stokes (RANS) equations, which assume the statistical time-averaged equations of motion for fluid flow. DNS fully resolves the turbulence spectrum without any modelling, requiring a very fine mesh to capture even the smallest scale of turbulence. LES necessitates fine discretizations near the body due to small turbulence scales present in the thin boundary layer along the walls [13].

The RANS method provides good turbulence resolution and handles complex geometries without needing excessively fine wall level discretizations. These approximation theories [13] are based on the physical interpretation of the eddy viscosity models with time-averaged velocity gradients for specific scientific and engineering applications. Although topology optimization of fluid flows has predominantly focused on steady flows, unsteady flows are common in reality, and research on optimizing unsteady flows using topology optimization method has been conducted.

Based on the standard approach proposed by Petersson and Borrvall [1], Kreissl et al. [14] proposed a topology optimization method for unsteady flow channel design problems, formulating the optimization problem as one of minimizing mean pressure drop minimization problem. Similarly, Deng et al. [15] examined how optimized channel configurations depend on inlet velocity profiles. This optimization procedure employs the continuous adjoint method and the finite element method. Numerical examples illustrate the feasibility of this topology optimization method for unsteady Navier-Stokes flows. Furthermore, Deng et al. [16] extended their approach to unsteady flow channel design problems, considering body forces driving fluid flows within the channels. The topology optimization problem is analyzed by the continuous adjoint method and solved by the finite element method combined with gradient-based methods. The numerical results confirm that the method effectively achieve topology optimization of flows driven by body forces. For the optimization of unsteady flows, a straightforward global implementation requires the entire flow solution history to be available during the reversed time calculation of the adjoint equations. Consequently, three-dimensional optimization problems are extremely challenging due to their massive computational cost. This has been a significant research topic.

Yamaleev et al. [17] presented a new local-in-time (LT) discrete adjoint-based methodology for solving design optimization problems in unsteady aerodynamic applications. This new methodology avoids the storage requirements associated with the global adjoint-based optimization, which stores the entire flow solution history for all time levels. They demonstrated that the method could converge to a local minimum of the original unsteady optimization problem. Yaji et al. [18] showed that incorporating the LT method, the lattice Boltzmann method, and parallel computing could solve large-scale topology optimization problems. Notably, some discretization methods, such as lattice Boltzmann methods, are inherently transient. Chen et al. [19] presented a local-in-time discrete adjoint-based topology optimization method for unsteady incompressible viscous flows incorporating the lattice Boltzmann method. They approximate the true time-dependent adjoint equations by splitting the time series into subintervals rather than using the full-time series forwards and then backwards. This significantly reduces the storage requirements but introduces an approximation error due to the discontinuity between subintervals, where adjoint variables may not fully capture the time dynamics of the flow. Some researchers, such as Yamaleev et al. [17], Theulings et al. [20], explored similar local-in-time adjoint methods and highlighted that while storage is reduced, the accuracy of the optimization process is compromised. Recently, Nobis et al. [21] investigate the applicability of a high-order Spectral Element Method (SEM) to density-based topology optimization of unsteady flows in two dimensions. Direct Numerical Simulations (DNS) were conducted using Brinkman penalization to describe the presence of solid within the domain. The optimization procedure utilized the adjoint-variable method to compute gradients and a checkpointing strategy to reduce storage requirements.

Optimization of turbulent flows has only been treated in a few publications [11]. Although turbulence is inherently time-dependent, many studies on topology optimization of turbulent flow restrict themselves to the steady-state time-averaged approximation of turbulence, specifically using the RANS equations. In the context of topology optimization for turbulent flows, Other [22], Pietropaoli et al. [23] employed a simplified sensitivity analysis through the so-called "frozen turbulence" assumption, which neglects variations in the turbulence field with respect to design variables. However, following the derivation of a continuous adjoint formulation for the Spalart-Allmaras turbulence model [24], Kontoleontos et al. [25] demonstrated that the typical "frozen turbulence" assumption can result in sensitivities with incorrect sign in some cases. Besides Dilgen et al. [26] differentiate the turbulence models and demonstrate that frozen turbulence can lead to wrong sensitivities. Yoon [27] presented a discrete adjoint approach to density-based optimization of turbulent flow problems using the Spalart-Allmaras turbulence model and a modified wall equation. However, this study used meshes that were too coarse to capture turbulence accurately. Subsequently, Yoon [28] investigated the impact of turbulent flow in topology optimization by employing a $k-\epsilon$ turbulent finite element model, one of the RANS equations. They designed 2D flow components minimizing turbulent energy, and thus reducing noise.

In addition to the computational and data storage costs associated with optimizing unsteady flows, turbulent flows exhibit significant dynamics, necessitating careful selection from various models in the literature. Extending topology optimization to unsteady Navier-Stokes flows is therefore highly desirable. In this paper, we present a novel algorithm based on a Reynolds-Averaged Method (RAM) to solve the unsteady topology optimization problem for buoyancy-driven turbulent flows. We solve the incompressible Navier-Stokes equations and instantaneous quantities are produced. These are used to build mean velocities, mean temperatures and mean tensors over iterations. Once these mean quantities stabilize, we estimate turbulent viscosity and turbulent thermal diffusivity. We obtain the model based on Reynolds-Averaged Method with stationary equations whose resolution closely approximates the averaged quantities. This RAM model is then used as a constraint to solve the optimization problem. The process is repeated until the optimization converges. Our algorithm requires only the means of the physical variables, and the adjoint problem is solved in only one iteration.
tion per algorithm loop. We apply this algorithm to the optimization of unsteady laminar and 2D/3D turbulent flows. The method is illustrated on three configuration cases: an unsteady two dimensional laminar flow around a circular cylinder with a Reynolds number of $Re = 140$, the classical Pitz-Daily turbulence case where $Re = 12700$ and a turbulent anisothermal flow by adapting the previous Pitz-Daily case. We demonstrate the algorithm capability to reduce pressure losses and enhance heat transfer by adding material within the computational domain and show that the unsteady laminar flow becomes steady due to the optimized design. We report the number of loops required by the algorithm to solve the optimization problem. That corresponds to the number of updates for the design parameter. The value of the cost function obtained after optimization is also specified.

The paper is organized as follows: Section 2 formulates the general topology optimization problem with incompressible Navier-Stokes equations under the Boussinesq assumption. Section 3 presents the algorithm and its main steps. Section 4 details the Reynolds-Averaged Method for solving topology optimization of unsteady flows. In Section 5, we illustrate our algorithm with the three case studies.

## 2. Topology optimization formulation

Although the methodology presented in our paper can be adapted to any transport model with a constitutive equation, we focus on the incompressible Navier-Stokes equations under the Boussinesq assumption.

$$ egin{cases} 
\nabla \cdot \vec{u} = 0, \\
\partial_t \vec{u} + \nabla \cdot (\vec{u} \otimes \vec{u}) = -\nabla p + A \Delta \vec{u} - h_r(\alpha)\vec{u} + f(\alpha) B \theta \vec{e}_y, \\
\partial_t \theta + \nabla \cdot (\theta \vec{u}) = C \nabla \cdot (k_c(\alpha)\theta),
\end{cases} 
\quad (1) $$

where $\vec{u}, p$ and $\theta$ are respectively the dimensionless velocity, pressure and temperature, $t$ is the time. The parameter $\alpha$ represents the effective inverse permeability of the solid. As stated in the introduction, $\alpha$ takes the value one for cells filled with impermeable solid material and zero for fluid. Intermediate values can be interpreted as porous material with different porosities. The goal of the optimization is to find the porosity distribution that minimizes the objective function and satisfies the set of constraints. The functions $x \mapsto \{h_r(\alpha(x)), k_c(\alpha(x))\}$ are both smooth regularization of Heaviside step functions, whose expressions can be found, for instance, in [29]. These functions facilitate the interpolation between the fluid and the solid zones in the computational domain and allow to work with an optimization parameter taking continuous values $\alpha: x \in \Omega \mapsto \alpha(x) \in [0, \alpha_{\max}]$ instead of $\vec{\alpha}: x \in \Omega \mapsto \vec{\alpha}(x) \in \{0, \alpha_{\max}\}$. However, we do not provide their specific definitions here since the methods introduced in the work do not require particular interpolation functions. The function $\alpha \mapsto f(\alpha) = 1 - \alpha/\alpha_{\max}$ cancels the effects of buoyancy in the solid zones, where $\alpha_{\max}$ represents the maximum of the inverse permeability. The physical constants A, B, C for mixed natural convection are given as

$$ A = \mathrm{Re}^{-1}, \quad B = \mathrm{Ri}, \quad C = \mathrm{Re}^{-1} \mathrm{Pr}^{-1}, $$

with Pr is the Prandlt number and Ri is the Richardson number.

Assuming the boundary $\partial\Omega = \Gamma$ can be decomposed as $\Gamma = \Gamma_{in} \cup \Gamma_{out} \cup \Gamma_{w} \cup \Gamma_{d}$, where $\Gamma_{in}$, $\Gamma_{out}$, $\Gamma_{w}$, and $\Gamma_{d}$ are respectively the inlet, the outlet, the heated walls, and the adiabatic walls of the domain, we consider the following classical boundary conditions:

$$ egin{cases} 
\vec{u} = \vec{u}_{\mathrm{in}}, & \theta = 0, & \text{on } \Gamma_{\mathrm{in}}, \\
\vec{u} = 0, & \partial_n \theta = 0, & \text{on } \Gamma_{\mathrm{d}}, \\
\vec{u} = 0, & \theta = \Phi, & \text{on } \Gamma_{\mathrm{w}}, \\
\partial_n \vec{u} = 0, & \partial_n \theta = 0, & p=0 \text{ on } \Gamma_{\mathrm{out}}.
\end{cases} 
\quad (2) $$

where $\Phi$ is a given heat flux and $\vec{u}_{\mathrm{in}}$ is a given inlet velocity. The initial conditions are denoted by

$$ \vec{u}(0, x) = \vec{u}_0(x), \quad \theta(0, x) = \theta_0(x), \quad x \in \Omega, $$

whose expressions will be specified when needed. In this paper, we present an algorithm to solve the next general time-dependent topology optimization problem

$$ egin{gathered} \min \mathcal{J}(\alpha, \vec{u}, \theta, p) \\ \text{s.t. } (\vec{u}, \theta, p) \text{ solution of (1) – (2) parametrized by } \alpha, \end{gathered} 
\quad (3) $$

where $\mathcal{J}$ is a given cost function such as pressure losses or heat transfers, which can be written as:

$$ \mathcal{J}(\alpha; \vec{u}, \theta, p) = \int_{0}^{T} \left( \int_{\Omega} J_{\Omega}(\alpha, \vec{u}, \theta, p) d\Omega + \int_{\Gamma} J_{\Gamma}(\alpha, \vec{u}, \theta, p) d\Gamma \right) dt. $$

with $T$ being the final time considered in our numerical simulations.

## 3. The time-averaged algorithm

We consider an unsteady flow that, on average, becomes steady over time. Our algorithm addresses the general time-dependent topology optimization problem (3) by iteratively computing averaged physical quantities to define a model (referred below as the RAM model) which is then used as constraint in the optimization problem. Topological optimization involves strategically adding solid materials to achieve specific objectives. The proposed algorithm is presented in three stages: firstly, a schematic diagram (Fig. 1)
<figure><img src="image_4.png" /><figcaption>Fig. 1. Diagram illustrating the RAM method with the steps indicated in our Algorithm 3.1.</figcaption></figure>

offers a visual overview; secondly, a step-by-step text provides explanation; thirdly, a main flowchart of the proposed algorithm is introduced (Algorithm 3.1); and finally, each step is examined in detail in the following sections.

We now explain each step of the Algorithm 3.1.

**Step 1:** This step relies on averaged physical quantities noted ($\bar{\bar{u}}, \bar{\bar{\theta}}, \bar{\bar{p}}$). Instantaneous quantities are calculated using either a temporal turbulence model or a DNS model. At each time step $\Delta t$, we monitor the change in mean velocity from its previous value over preceding time step. This condition is typically satisfied after a certain amount of time steps, indicating that the average flow variation with time is minimal. Denoting by $\bar{\bar{\psi}}_{[t_1, t_2]} = \frac{1}{t_2-t_1} \int_{t_1}^{t_2} \psi(s) ds$ the mean value of $\psi$ over $[t_1, t_2]$, it is worth noting that we have

$$ \begin{aligned} \bar{\bar{\psi}}_{[t_0, t+\Delta t]} - \bar{\bar{\psi}}_{[0, t]} &= \frac{1}{t+\Delta t} \int_0^{t+\Delta t} \psi(s) ds - \frac{1}{t} \int_0^t \psi(s) ds \\ &= \left( \frac{1}{t+\Delta t} - \frac{1}{t} \right) \int_0^t \psi(s) ds + \frac{1}{t+\Delta t} \int_t^{t+\Delta t} \psi(s) ds \\ &= -\frac{\Delta t}{t+\Delta t} \bar{\bar{\psi}}_{[0, t]} + \frac{\Delta t}{t+\Delta t} \bar{\bar{\psi}}_{[t, t+\Delta t]} \end{aligned} $$

**Algorithm 3.1:** Main flowchart of the proposed algorithm.

**Data:** Set $\alpha_0 = 0$, initial condition $\bar{\bar{u}}_0, \bar{\bar{p}}_0$ and $\theta_0$, time step $\Delta t$ and tolerances $\text{tol}_u, \text{tol}_p, \text{tol}_\theta, \text{tol}_J$.

**Step 1:** Find $t$ such that

$$ \| \bar{\bar{u}}_{[0,t+\Delta t]} - \bar{\bar{u}}_{[0,t]} \|_{L^2(\Omega)} \le \text{tol}_u, $$

$$ \| \bar{\bar{p}}_{[0,t+\Delta t]} - \bar{\bar{p}}_{[0,t]} \|_{L^2(\Omega)} \le \text{tol}_p, $$

$$ \| \bar{\bar{\theta}}_{[0,t+\Delta t]} - \bar{\bar{\theta}}_{[0,t]} \|_{L^2(\Omega)} \le \text{tol}_\theta, $$

where, by definition

$$ \bar{\bar{\psi}}_{[t_1, t_2]}(x) = \frac{1}{t-t_1} \int_{t_1}^{t_2} \psi(t, x) dt $$

and $(\bar{\bar{u}}, \bar{\bar{p}}, \bar{\bar{\theta}})$ satisfies the transport model with the associated constitutive equations, parameterized by $\alpha_0$ on the time interval $[t_0, t+\Delta t]$ with $\bar{\bar{u}}|_{t=0} = \bar{\bar{u}}_0$.

Compute $\bar{\bar{J}}_0 = \mathcal{J}(\bar{\bar{u}}, \bar{\bar{p}}, \bar{\bar{\theta}})$.

**Step 2:** Define a RAM model where the turbulent viscosity and turbulent thermal diffusivity are replaced by estimated turbulent quantities which only involve means of physical variables.

**Step 3:** Compute $\alpha$ that (at least approximately) minimizes $\mathcal{J}_{\text{ram}}$ with the RAM model from **Step 2** as a constraint.

**Step 4:**

if $\| \bar{\bar{J}}_{[0,t+\Delta t]} - \bar{\bar{J}}_{[0,t]} \|_{L^2(\Omega)} \le \text{tol}_J$ then

else

Return to **Step 1**;

end
and thus, as expected, the condition involved in Step 1 could be verified for small enough $\Delta t$ even for flows whose average is highly varying in time.

**Step 2:** In this step, we define a steady model where turbulent quantities, i.e. turbulent viscosity and turbulent thermal diffusivity, are replaced by estimated turbulent parameters bases solely on the mean physical quantities variables. Since we rely here on a Reynolds Averaged Method (RAM) to define this steady model, we called this steady PDE the RAM model.

The RAM model aims to streamline the topology optimization process by reducing the complexity associated with solving the full unsteady turbulent flow equations. The averaging process leads to additional unknown quantities and requires the definition of a closure problem Our main strategy at this step is the estimation of turbulence parameters, such as turbulent viscosity and turbulent thermal diffusivity. It is important to note that the average quantities are computed from the instantaneous resolved quantities, rather than modelled quantities. Additionally, regardless of the system of partial differential equations resulting from a transport model combined with a constitutive law, instantaneous quantities produced are intended to represent the fields of velocity, pressure, and temperature at a certain spatial and temporal scales, and could be used to define a RAM model (which aims to best approximate the average quantities of the Navier-Stokes solutions) and then proceed with the optimization.

**Step 3:** Here, we update the design variable $\alpha$ by minimizing the cost function $J_{\text{ram}}$ which is defined from the cost function $J$ through the identification

$$ 
\begin{aligned}
J(\alpha, \bar{u}, \theta, p) &= \int_{0}^{T} \left( \int_{\Omega} J_{\Omega}(\alpha, \bar{u}, \theta, p) d\Omega + \int_{\Gamma} J_{\Gamma}(\alpha, \bar{u}, \theta, p) d\Gamma \right) dt \\
&= \int_{0}^{T} (J_{\text{ram}}(\alpha, \bar{u}, \theta, p)) dt
\end{aligned}
$$

hence giving

$$ J_{\text{ram}}(\alpha, \bar{u}, \bar{\theta}, \bar{p}) = \int_{\Omega} J_{\Omega}(\alpha, \bar{u}, \bar{\theta}, \bar{p}) dx + \int_{\Gamma} J_{\Gamma}(\alpha, \bar{u}, \bar{\theta}, \bar{p}) d\Gamma. $$

This optimization problem is now constrained by the RAM model established in Step 2. We solve this PDE-constrained optimization problem with gradient-based algorithm where the gradient of the cost function is computed using the adjoint method. The adjoint equations of the RAM model are given in A and solving both the direct and adjoint RAM models allows the design variable $\alpha$ to be updated.

**Steps 4:** Finally, we iterate the process from Step 1 until the cost function $J$ between two iterations does not vary more than a specified tolerance $\text{tol}_J$. This iterative refinement updates the design variable $\alpha$ and adjusts the flow variables $(\bar{u}, p, \theta)$ accordingly to achieve the desired optimization objectives.

Our method introduces several distinctive features compared to existing approaches [15]:

* Since we update the design variable $\alpha$ though the minimization of the RAM cost function constrained by the RAM model, our method makes it possible to build sensitivities and perform topology optimization, independently of the turbulence model and the model class (DNS, LES, RANS, etc.)

* We maintain a classical scheme, namely a finite volume method, for solving both direct and adjoint problems. This ensure robustness and reliability in the optimization process and significantly reduces computational costs.

* In our RAM algorithm, each update of the design variable is based on solving a topology optimization problem with a steady PDE as constraint. This contrasts with classical methods used to solve unsteady topology optimization problems which require solving unsteady direct and adjoint problems to update the design variable.

## 4. Reynolds-averaged method for topology optimization of unsteady flows

Turbulent flow, characterized by its multiscale nature, presents challenges for computational methods due to its complex dynamics. Among the various theories and numerical approaches proposed in the literature, the Reynolds-Averaged Navier-Stokes equations (RANS) and Large Eddy Simulation (LES) are prominent in engineering applications [26,27]. These methods utilize time-averaged equations of motion to model fluid flow, suitable for stationary turbulence where the flow statistics remain steady over time. The key concept involves Reynolds decomposition, which decomposes the flow into a mean component and fluctuations.

### 4.1. Reynolds decomposition

Reynolds decomposition, introduced by Reynolds [30] considers a smooth function $\varphi \in [0, T] \times \Omega \mapsto \varphi(t, x) \in \mathbb{R}$ and the time-averaged or Reynolds-averaged quantity $\bar{\varphi}$ is defined as:

$$ \bar{\varphi}_{[t_1, t_2]}(x) = \frac{1}{t_2 - t_1} \int_{t_1}^{t_2} \varphi(t, x) dt. $$

Since, in this section, we do not need to specify the time interval on which the mean value is computed, we thus drop the subscript $[t_1, t_2]$ and only rely on the notation $\bar{\varphi}_{[t_1, t_2]} = \bar{\varphi}$. The Reynolds decomposition of $\varphi$ writes it as the sum of a steady mean component $\bar{\varphi}(x)$ and a time varying fluctuating component $\varphi'(t, x)$ with zero mean value. It is given as

$$ \varphi(t, x) = \varphi'(t, x) + \bar{\varphi}(x). \quad (4) $$
This decomposition allows us to separate the flow into its time-averaged mean and fluctuating parts. This time-averaging satisfies the next properties [13] that can be checked by direct computations

$$ 
\begin{aligned}
\overline{\overline{\varphi}} &= \overline{\varphi}, \quad \overline{\varphi}' = 0, \quad \overline{\varphi' \overline{\varphi}} = 0, \\
\overline{\partial_j \varphi} &= \partial_j \overline{\varphi}, \\
\partial_t \varphi' &= (\partial_t \varphi)' , \quad \overline{\partial_t \varphi} = \overline{\partial_t \varphi'} = 0,
\end{aligned} 
\qquad (5) $$

where $\varphi' = \varphi - \overline{\varphi}$.

We then form the mean of the continuity and Navier-Stokes equations term by term. The RAM model associated to equations (1) is now defined as the problem satisfied by $(\overline{\overline{u}}, \overline{\overline{\theta}}, \overline{\overline{p}})$. The latter can be obtained by incorporating the Reynolds decomposition and next by time-averaging Navier-Stokes equations (1). This gives the next system of partial differential equations:

$$ 
\left\{
\begin{aligned}
\nabla \cdot \overline{\overline{u}} &= 0, \\
\nabla \cdot (\overline{\overline{u}}' \otimes \overline{\overline{u}}') + \nabla \cdot (\overline{\overline{u}} \otimes \overline{\overline{u}}) &= -\nabla \overline{p} + A \Delta \overline{u} - h_{\tau}(\alpha) \overline{u} + f(\alpha) B \overline{\theta} \overline{e}_{y}, \\
\nabla \cdot (\overline{\theta}' \overline{\overline{u}}') + \nabla \cdot (\overline{\theta} \overline{\overline{u}}) &= C \nabla \cdot (k_{\tau}(\alpha) \nabla \overline{\theta}).
\end{aligned}
\right.
\qquad (6) $$

This averaging process leads to the appearance of additional unknown quantities. Indeed, the Reynolds stress $\overline{\overline{u}}' \otimes \overline{\overline{u}}'$ and the turbulent heat flux $\overline{\theta}' \overline{\overline{u}}'$ in equations (6), are unclosed terms that require a closure.

## 4.2. Closure problem

To close the previous equations (6), we employ the Boussinesq approximation [31] for the Reynolds stresses and an analogous approach for the turbulent heat flux. Specifically, we approximate the unclosed terms as follows:

$$ 
\begin{aligned}
\overline{\overline{u}}' \otimes \overline{\overline{u}}' &\approx -2v_{\text{ram}}S(\overline{\overline{u}})^d, \\
\overline{\theta}' \overline{\overline{u}}' &\approx -k_{\text{ram}}\nabla\overline{\theta},
\end{aligned}
 $$

where $S(\overline{\overline{u}}) = (\nabla \overline{u} + \nabla \overline{u}')/2$ is the rate-of-strain tensor, and $v_{\text{ram}}$ and $k_{\text{ram}}$ are the turbulent viscosity and turbulent thermal diffusivity, respectively. The superscript $d$ represents the deviating part of the considering tensors. The first equation represents the Boussinesq hypothesis, which relates the Reynolds stress to the mean strain rate through a turbulent viscosity $v_{\text{ram}}$. This hypothesis is a common approach in modelling turbulent flows across various fields [13]. The second equation is an analogous expression for the turbulent heat flux, where $k_{\text{ram}}$ is the turbulent thermal diffusivity or eddy diffusivity.

Using these approximations, the equations (6) can be closed by determining the turbulent parameters $v_{\text{ram}}$ and $k_{\text{ram}}$. We achieve this by minimizing the residuals $E_{\text{v}}$ and $E_{\text{k}}$ which quantify the errors in the approximations:

$$ 
\begin{aligned}
E_{\text{v}} &= \overline{\overline{u}}' \otimes \overline{\overline{u}}' + 2v_{\text{ram}}S(\overline{\overline{u}})^d, \\
E_{\text{k}} &= \overline{\overline{\theta}}' \overline{\overline{u}}' + k_{\text{ram}}\nabla\overline{\theta}.
\end{aligned}
\qquad (7) $$

We compute $v_{\text{ram}}$ and $k_{\text{ram}}$ by minimizing respectively $E_{\text{v}} : E_{\text{v}}$ and $E_{\text{k}} : E_{\text{k}}$. This actually amounts to compute

$$ 
\begin{aligned}
\frac{\partial(E_{\text{v}} : E_{\text{v}})}{\partial v_{\text{ram}}} &= 0, \\
\frac{\partial(E_{\text{k}} : E_{\text{k}})}{\partial k_{\text{ram}}} &= 0,
\end{aligned}
\qquad (8) $$

from which we obtain

$$ 
\begin{aligned}
v_{\text{ram}} &= - \frac{(\overline{\overline{u}}' \otimes \overline{\overline{u}}')^d : S(\overline{\overline{u}})^d}{2S(\overline{\overline{u}})^d : S(\overline{\overline{u}})^d} \quad \text{and} \quad k_{\text{ram}} = - \frac{\overline{\overline{\theta}}' \overline{\overline{u}}' \cdot \nabla \overline{\theta}}{\overline{\nabla \overline{\theta}} \cdot \overline{\nabla \overline{\theta}}},
\end{aligned}
\qquad (9) $$

However, these expressions require storing all quantities at all time steps of a simulation to determine the fluctuating components $\overline{\overline{u}}' = \overline{\overline{u}} - \overline{\overline{u}}$ and $\overline{\overline{\theta}}' = \overline{\overline{\theta}} - \overline{\overline{\theta}}$. This can be impractical for simulations with a large number of discretization cells, as noted by several researchers [17–19,21]. Thus, while Equation (9) provides a way to determine $v_{\text{ram}}$ and $k_{\text{ram}}$, it also highlights the challenges associated with storing and processing large datasets. These challenges necessitate efficient computational strategies to manage and utilize the fluctuating quantities effectively $\overline{\overline{u}}'$ and $\overline{\overline{\theta}}'$.

By using a result of the Reynolds decomposition, one has $(\overline{\overline{u}}' \otimes \overline{\overline{u}}') = (\overline{\overline{u}} \otimes \overline{\overline{u}}) - (\overline{\overline{u}} \otimes \overline{\overline{u}})$ which is easily computable during a numerical simulation without excessive storage. For the energy equation, the expression for the term including fluctuation is expressed as $\overline{\overline{\theta}}' \overline{\overline{u}}' = \overline{\overline{\theta}} \overline{\overline{u}} - \overline{\overline{\theta}} \overline{\overline{u}}$. Equations (9) thus assume the following expressions:

$$ 
\begin{aligned}
v_{\text{ram}} &= - \frac{[\overline{\overline{u}} \otimes \overline{\overline{u}}] - (\overline{\overline{u}} \otimes \overline{\overline{u}})^d}{2S(\overline{\overline{u}})^d : S(\overline{\overline{u}})^d},
\end{aligned}
 $$
$$k_{\text{ram}} = - \frac{[\vec{\theta}\vec{u} - \vec{\theta}\vec{u}] \cdot \nabla \vec{\theta}}{\nabla \vec{\theta} \cdot \nabla \vec{\theta}}. \qquad (10)$$

The expressions $k_{\text{ram}}$ and $v_{\text{ram}}$ given in equations (10) provide a practical way to compute these turbulent parameters without explicitly storing all quantities at each time step.

It is worth noting that equations (10) depend on the physical quantities and their average quantities. The fact that they do not depend solely on mean quantities is an additional obstacle to overcome if we wish to consider in our development a dependence of this viscosity on mean quantities. This dependence can then be estimated by exploiting the immense literature in the field of turbulent viscosity modelling initiated previously by [32–35]. In this study, we propose to exploit an algebraic model as a function of $\vec{u}$ and $\vec{u}\vec{\theta}$.

### 4.3. Least-square determination of estimated turbulent parameters

As an alternative to the frozen turbulence concept [22], our strategy consists to estimate the expressions of (10) from the quantities $(\vec{u}, \vec{\theta})$ which are more reachable. In the following, the estimated quantities will be noted with a superscript E. We now wish to determine constants $c_j$ such that $v_{\text{ram}}$ and $k_{\text{ram}}$ can be approximated as follow

$$v_{\text{ram}} \approx v_{\text{ram}}^{E} = c_{0} + \vec{H}, \quad k_{\text{ram}} \approx k_{\text{ram}}^{E} = c_{0} + \vec{G}, \qquad (11)$$

where $\vec{H} = \sum_j c_j v_j$, $\vec{G} = \sum_j c_j T_j$ and where $v_j, T_j$ are fixed functions that only involve averaged quantities of $(\vec{u}, \theta, p)$. We emphasize that $\vec{H}$ belongs to $\text{span}(S)$ which is the set of linear combinations of elements of $S := \{v_1, \dots, v_N\}$. The same holds for $\vec{G}$. One may choose for instance

$$S := \{ |\vec{u}|^j, |S(\vec{u})|^k, |\nabla \vec{\theta}|^l, j = 1, \dots, N_j, l = 1, \dots, N_l, k = 1, \dots, N_k \},$$

for chosen $N_j, N_k, N_l$. The set $\mathcal{T}$ can be defined similarly. We emphasize that many other examples of sets $\mathcal{T}$ and $S$ can be used. Indeed, we are going to define our estimated turbulent quantities $v_{\text{ram}}^E, k_{\text{ram}}^E$ as 11 where the constants $\vec{C} = (c_0, c_1, \dots, c_N)$ are determined so that the next functions

$$f(\vec{C}) = \int_{\Omega} (v_{\text{ram}} - v_{\text{ram}}^{E}(\vec{C}))^2 d\Omega,$$

$$g(\vec{C}) = \int_{\Omega} (k_{\text{ram}} - k_{\text{ram}}^{E}(\vec{C}))^2 d\Omega,$$

are minimized. Alternatively, this reduces to chose $\vec{H} = \sum_j c_j^* v_j, \vec{G} = \sum_j d_j^* T_j$ where $\vec{C}^*, \vec{D}^*$ are given as

$$\vec{C}^* = \arg \min_{\vec{C} \in \mathbb{R}^N} (f(\vec{C})), \quad \vec{D}^* = \arg \min_{\vec{C} \in \mathbb{R}^N} (g(\vec{C})). \qquad (12)$$

Note that $\vec{C} \mapsto \psi_l^E(\vec{C})$ with $\psi \in \{v_{\text{ram}}, k_{\text{ram}}\}$ are affine functions. Since the application $N : s \mapsto \int_{\Omega} (v_{\text{ram}} - s)^2 d\Omega$ is convex, we obtain that $\vec{C} \mapsto f(\vec{C})$ (idem for $g$) are convex functions. As a result, $\vec{C}^*$ is given by the solution to

$$\nabla \vec{C} f(\vec{C}^*) = \vec{0}. \qquad (13)$$

A computation gives

$$\begin{aligned} \nabla_{\vec{C}} f(\vec{C}^*) &= N'(v_{\text{ram}}^E(\vec{C}^*)) \nabla_{\vec{C}} v_{\text{ram}}^E \\ &= -2 \left( \int_{\Omega} (v_{\text{ram}} - v_{\text{ram}}^E) \nabla_{\vec{C}} v_{\text{ram}}^E d\Omega \right) \\ &= -2 \left( \int_{\Omega} v_{\text{ram}} \nabla_{\vec{C}} v_{\text{ram}}^E d\Omega - \int_{\Omega} v_{\text{ram}}^E \nabla_{\vec{C}} v_{\text{ram}}^E d\Omega \right). \end{aligned}$$

The condition (13) reduces then to solve

$$\int_{\Omega} v_{\text{ram}}^{E} \nabla_{\vec{C}} v_{\text{ram}}^{E} d\Omega = \int_{\Omega} v_{\text{ram}} \nabla_{\vec{C}} v_{\text{ram}}^{E} d\Omega. \qquad (14)$$

We emphasize that $\nabla_{\vec{C}} v_{\text{ram}}^E$ is a vector that only depends on averaged quantities but not on $\vec{C}$. Therefore, (14) is a linear system. Remark that, at this point, $v_{\text{ram}}$ is already determined (see Step 3 of Algorithm 3.1) since we have computed the instantaneous physical quantities and thus the estimated turbulent quantities are fully determined and also only depend on time-averaged variables.

**Remark 1.** Using that $v_{\text{ram}}^E = \sum_{j=0}^N c_j v_j$, we have $\nabla_{\vec{C}} v_{\text{ram}}^E = (v_i)_{i=1,...,N}$ and the linear system (14) can be recast as

$$\sum_{j=0}^{N} c_j \int_{\Omega} v_i v_j d\Omega = \int_{\Omega} v_{\text{ram}} v_i d\Omega.$$

In matrix form, this reads as

$$K \vec{C} = \vec{F},$$
with

$$K = \left( \int_{\Omega} \nu_1 \nu_j d\Omega \right)_{i,j=0,...,N}, \quad \vec{F} = \int_{\Omega} \nu_{ram} \nu_j d\Omega.$$

We emphasize that K is then the Gram matrix associated to {$1, \nu_1, ..., \nu_N$} and is thus invertible as long as the vectors are linearly independent in the Hilbert $L^2(\Omega)$ equipped with the usual inner product.

**Remark 2.** We can take any non-negative, differentiable convex function F to define $\vec{C}^*$ that minimize $F(v_{ram} - v_{ram}^E)$. In that case, the optimal parameters also satisfy

$$\nabla \vec{C} F(\vec{C}^*) = 0.$$

Note however, that this may not lead to a linear system.

**Remark 3 (Examples).** We give here some examples of the previous procedure in the cases $v_{ram}^E = C f(\bar{u}, \bar{\theta}, \bar{p})$ for a given function f and $v_{ram}^E = C_1 + C_2 |\bar{u}|^2$.

In case we consider $v_{ram}^E = C f(\bar{u}, \bar{\theta}, \bar{p})$ then (14) gives that the optimal constant is

$$C_{opt} = \frac{1}{\int_{\Omega} f(\bar{u}, \bar{\theta}, \bar{p}) d\Omega} \int_{\Omega} v_{ram} d\Omega. \qquad (15)$$

We now look for $v_{ram}^E$ given by $v_{ram}^E = C_1 + C_2 |\bar{u}|^2$. Thanks to (14), these constants satisfy the next linear system

$$\begin{pmatrix} \gamma_1 & \gamma_2 \\ \gamma_2 & \gamma_3 \end{pmatrix} \begin{pmatrix} C_1^* \\ C_2^* \end{pmatrix} = \begin{pmatrix} \beta_1 \\ \beta_2 \end{pmatrix}$$

where

$$\begin{align*}
\gamma_1 &= |\Omega|, \quad \gamma_2 = \int_{\Omega} |\bar{u}|^2 d\Omega, \quad \gamma_3 = \int_{\Omega} |\bar{u}|^4 d\Omega, \\
\beta_1 &= \int_{\Omega} v_{ram} d\Omega, \quad \beta_2 = \int_{\Omega} v_{ram} |\bar{u}|^2 d\Omega.
\end{align*}$$

We have

$$C_1^* = \frac{\gamma_3\beta_1 - \gamma_2\beta_2}{\gamma_1\gamma_3 - \gamma_2^2}, \quad C_2^* = \frac{\gamma_1\beta_2 - \gamma_2\beta_1}{\gamma_1\gamma_3 - \gamma_2^2}. \qquad (16)$$

In the subsequent part of our article, we are going to use $v_{ram}^E$ and $k_{ram}^E$ defined as :

$$\begin{align}
v_{ram}^E &= C_2 |\bar{u}|^2, \nonumber \\
k_{ram}^E &= C_{2,\rho} |\bar{u}| |\bar{\theta}|, \tag{17}
\end{align}$$

where the optimal constants have been computed in Equation (15).

We emphasize that there are a lot of different models that could be used to approach $v_{ram}^E$ and $k_{ram}^E$. We chose here these two estimated quantities because we find them to be mathematically simple and an effective way to illustrate the method's potential.

### 4.4. Optimization problem with the closed RAM model and estimated turbulent parameters

We now introduce the estimated turbulent quantities $v_{ram}^E, k_{ram}^E$ that only depend on $(\bar{u}, \bar{p}, \bar{\theta})$ in equations (6). We neglected here dependence on $\alpha$ as these estimations have been done for a given $\alpha$ which will then be updated after solving an optimization problem with the steady-RAM PDE as constraint. As a result, the closed RAM model considered in Step 2 from Algorithm 3.1 can be written as

$$\left\{
\begin{array}{l}
\displaystyle \nabla \cdot \bar{u} = 0, \\
\\
\displaystyle -\nabla \cdot (\bar{u} \otimes \bar{u}) - \nabla \bar{p} + 2\nabla \cdot (\{A + v_{ram}^E\} S(\bar{u})) - h_{\tau}(\alpha)\bar{u} + f(\alpha)\bar{B}\bar{\theta}\bar{e}_y = 0, \\
\\
\displaystyle -\nabla \cdot (\bar{\theta}\bar{u}) + \nabla \cdot (\{C k_{\tau}(\alpha) + k_{ram}^E\} \nabla \bar{\theta}) = 0,
\end{array}
\right.
\qquad (18)$$

where we used the incompressibility condition to get $A\Delta\bar{u} = A\nabla \cdot S(\bar{u})$. We remind that the idea here is to produce a suitable RAM model whose solution aims at approximating at best the time-averaged physical variable. In addition, such RAM model will be used as constraint in the minimization of the cost function $J_{ram}$ whose optimized solution will serve as update for the design variable $\alpha$.
Regarding the boundary conditions, we note that the time-averaging operator commutes with the boundary operator (the Dirichlet and Neumann traces). Therefore, the boundary conditions (2) translates to

$$ 
\begin{aligned}
\bar{\bar{u}} &= \bar{u}_{\text{in}}, & \bar{\bar{\theta}} &= 0, & \text{on } \Gamma_{\text{in}}, \\
\bar{\bar{u}} &= 0, & \bar{\bar{\theta}} &= \Phi, & \text{on } \Gamma_{\text{w}}, \\
\bar{\bar{u}} &= 0, & \partial_n \bar{\bar{\theta}} &= 0, & \text{on } \Gamma_{\text{d}}, \\
\partial_n \bar{\bar{u}} &= 0, & \partial_n \bar{\bar{\theta}} &= 0, & \bar{p} = 0, & \text{on } \Gamma_{\text{out}},
\end{aligned} 
\quad (19) $$

where $\bar{u}_{\text{in}}$ and $\Phi$ are the same source terms as for the initially considered time-dependent problem (2).

**Remark 4.** It is worth noting that the solution to (18)–(19), although denoted by $(\bar{\bar{u}}, \bar{\bar{\theta}}, \bar{\bar{p}})$, are not the mean values of the solution to the unsteady Navier-Stokes equations (1)-(2). Nevertheless, we keep the notation with the overbar since the RAM model has been built using the knowledge of the time-averaged instantaneous fields.

In Step 3 of our time-averaged algorithm, we need to solve the following optimization problem

$$ 
\begin{array}{l}
\min \mathcal{J}_{\text{ram}}(\alpha, \bar{\bar{u}}, \bar{\bar{p}}, \bar{\bar{\theta}}) \\
\text{s.t. } (\bar{\bar{u}}, \bar{\bar{p}}, \bar{\bar{\theta}}) \text{ solution of (18) – (19) parametrized by } \alpha,
\end{array} 
\quad (20) $$

where the cost function $\mathcal{J}_{\text{ram}}$ is obtained from the cost function $\mathcal{J}$ involved in the original topology optimization problem for unsteady flows (3). We recall we assumed it can be written as

$$ \mathcal{J}_{\text{ram}}(\bar{\bar{u}}, \bar{\bar{p}}, \bar{\bar{\theta}}) = \int_{\Omega} \bar{\bar{J}}_{\Omega}(\bar{\bar{u}}, \bar{\bar{p}}, \bar{\bar{\theta}}) d\Omega + \int_{\Gamma} \bar{\bar{J}}_{\Gamma}(\bar{\bar{u}}, \bar{\bar{p}}, \bar{\bar{\theta}}) d\Gamma. $$

We chose to solve (20) with an *optimize-then-discretize* approach [36] which allows to compute the gradient of the cost function with respect to $\alpha$ thanks to the continuous adjoint model. To compute the latter, we first introduce the Lagrangian associated to (20) which reads

$$ 
\begin{aligned}
\mathcal{L}(\vec{X}, \vec{X}^*) = & \mathcal{J}_{\text{ram}}(\alpha, \bar{\bar{u}}, \bar{\bar{\theta}}, \bar{\bar{p}}) - \int_{\Omega} \bar{\bar{p}}^* \nabla \cdot \bar{\bar{u}} d\Omega \\
& - \int_{\Omega} (-\nabla \bar{\bar{p}} + 2\nabla \cdot (\{A + v_{\text{ram}}^E(\vec{X})\} S(\bar{\bar{u}})) - h_{\tau}(\alpha)\bar{\bar{u}} + f(\alpha)B\bar{\bar{\theta}}\bar{\bar{e}}_y) \cdot \bar{\bar{u}}^* d\Omega \\
& - \int_{\Omega} (\nabla \cdot (\{C k_{\tau}(\alpha) + k_{\text{ram}}^E(\vec{X})\} \nabla \bar{\bar{\theta}})) \bar{\bar{\theta}}^* d\Omega \\
& - \int_{\Gamma_{\text{in}}} (\bar{\bar{u}} - \bar{\bar{u}}_{\text{in}}) \cdot \bar{\bar{\Phi}}_{\text{in}} + \bar{\bar{\theta}} \bar{\bar{\psi}}_{\text{in}} d\Gamma - \int_{\Gamma_{\text{w}}} (\bar{\bar{u}} \cdot \bar{\bar{\Phi}}_{\text{w}} + \partial_n \bar{\bar{\theta}} \bar{\bar{\psi}}_{\text{w}}) d\Gamma \\
& - \int_{\Gamma_{\text{d}}} (\bar{\bar{u}} \cdot \bar{\bar{\Phi}}_{\text{d}} + \partial_n \bar{\bar{\theta}} \bar{\bar{\psi}}_{\text{d}}) d\Gamma \\
& - \int_{\Gamma_{\text{out}}} (\partial_n \bar{\bar{u}} \cdot \bar{\bar{\Phi}}_{\text{out}} + \partial_n \bar{\bar{\theta}} \bar{\bar{\psi}}_{\text{out}} + \bar{\bar{p}} \bar{\bar{q}}_{\text{out}}) d\Gamma \\
& + \int_{\Omega} (\nabla \cdot (\bar{\bar{u}} \otimes \bar{\bar{u}}^*)) \cdot \bar{\bar{u}}^* d\Omega \\
& + \int_{\Omega} (\nabla \cdot (\bar{\bar{\theta}} \bar{\bar{u}})) \bar{\bar{\theta}}^* d\Omega
\end{aligned}
$$

where $\vec{X} = (\bar{\bar{u}}, \bar{\bar{\theta}}, \bar{\bar{p}})$ and $\vec{X}^* = (\bar{\bar{u}}^*, \bar{\bar{\theta}}^*, \bar{\bar{p}}^*, \mathcal{A}^*)$ are respectively the primal and adjoint variables. The additional adjoint variables

$$ \mathcal{A}^* = (\bar{\bar{\Phi}}_{\text{in}}, \bar{\bar{\psi}}_{\text{in}}, \bar{\bar{\Phi}}_{\text{w}}, \bar{\bar{\psi}}_{\text{w}}, \bar{\bar{\Phi}}_{\text{d}}, \bar{\bar{\psi}}_{\text{d}}, \bar{\bar{\Phi}}_{\text{out}}, \bar{\bar{\psi}}_{\text{out}}, \bar{\bar{q}}_{\text{out}}), $$

are only introduced to enforce the boundary conditions and are not going to be involved in the definition of the adjoint model (see e.g. [29, Remark 1]). The adjoint model on $\Omega$ is computed as detailed in A and reads as

$$ 
\left\{
\begin{array}{l}
-\nabla \cdot \bar{\bar{u}}^* + 2\partial_{\bar{\bar{p}}} v_{\text{ram}}^E (S(\bar{\bar{u}}): \nabla \bar{\bar{u}}^*) + \partial_{\bar{\bar{p}}} k_{\text{ram}}^E (\nabla \bar{\bar{\theta}} \cdot \nabla \bar{\bar{\theta}}^*) + \partial_{\bar{\bar{p}}} \bar{\bar{J}}_{\Omega} = 0, \\
\nabla \bar{\bar{p}}^* + h_{\tau}(\alpha) \bar{\bar{u}}^* - 2\nabla \cdot (\{A + v_{\text{ram}}^E(\vec{X})\} S(\bar{\bar{u}}^*)) - \bar{\bar{\theta}} \nabla \bar{\bar{\theta}}^* - 2S(\bar{\bar{u}}^*) \bar{\bar{u}} \\
+ 2\partial_{\bar{\bar{u}}} v_{\text{ram}}^E (S(\bar{\bar{u}}): \nabla \bar{\bar{u}}^*) + \partial_{\bar{\bar{u}}} k_{\text{ram}}^E (\nabla \bar{\bar{\theta}} \cdot \nabla \bar{\bar{\theta}}^*) + \partial_{\bar{\bar{u}}} \bar{\bar{J}}_{\Omega} = 0, \\
-(\nabla \cdot (\{C k_{\tau}(\alpha) + k_{\text{ram}}^E\} \nabla \bar{\bar{\theta}}^*)) - \bar{\bar{u}} \cdot \nabla \bar{\bar{\theta}}^* - f(\alpha)B\bar{\bar{e}}_y \cdot \bar{\bar{u}}^* \\
+ 2\partial_{\bar{\bar{\theta}}} v_{\text{ram}}^E (S(\bar{\bar{u}}): \nabla \bar{\bar{u}}^*) + \partial_{\bar{\bar{\theta}}} k_{\text{ram}}^E (\nabla \bar{\bar{\theta}} \cdot \nabla \bar{\bar{\theta}}^*) + \partial_{\bar{\bar{\theta}}} \bar{\bar{J}}_{\Omega} = 0.
\end{array}
\right. 
\quad (21) $$
For the boundary conditions on the adjoint velocity $\bar{u}$, we have

$$ 
\left\{ 
\begin{array}{l}
\bar{n} \cdot \bar{u}^* = -\partial_{\bar{p}} \bar{J}_{\Gamma}, \quad \bar{t} \cdot \bar{u}^* = 0 \qquad \text{on } \Gamma_{\text{in}} \cup \Gamma_{\text{w}} \cup \Gamma_{\text{d}}, \\
-\bar{n} \bar{p}^* + 2 \left\{ A + \nu_{\text{ram}}^{\text{E}}(\bar{X}) \right\} S(\bar{u}^*) \bar{n} \\
\qquad + \bar{n} \bar{\theta} \bar{\theta}^* + (\bar{u}^* \otimes \bar{n}) \bar{u} + (\bar{n} \otimes \bar{u}) \bar{u}^* \\
\qquad - \partial_{\bar{u}} k_{\text{ram}}^{\text{E}} \nabla \bar{\theta} \cdot \bar{n} \bar{\theta}^* - 2 \partial_{\bar{u}} \nu_{\text{ram}}^{\text{E}} S(\bar{u}) \bar{n} \cdot \bar{u}^* = -\partial_{\bar{u}} \bar{J}_{\Gamma} \quad \text{on } \Gamma_{\text{out}},
\end{array}
\right. 
\quad (22) 
$$

where $\bar{t}$ is some unit tangent vector. For the adjoint time-averaged temperature $\bar{\theta}^*$, we have the next boundary conditions:

$$ 
\left\{ 
\begin{array}{l}
\bar{\theta}^* = 0 \qquad \text{on } \Gamma_{\text{in}}, \\
\left\{ C k_{\tau}(\alpha) + k_{\text{ram}}^{\text{E}}(\bar{X}) \right\} \nabla \bar{\theta}^* \cdot \bar{n} + \bar{\theta}^* \bar{u} \cdot \bar{n} \\
- \partial_{\bar{\theta}} k_{\text{ram}}^{\text{E}} \nabla \bar{\theta} \cdot \bar{n} \bar{\theta}^* - 2 \partial_{\bar{\theta}} \nu_{\text{ram}}^{\text{E}} S(\bar{u}) \bar{n} \cdot \bar{u}^* = - \partial_{\bar{\theta}} \bar{J}_{\Gamma} \quad \text{on } \Gamma_{\text{w}} \cup \Gamma_{\text{d}} \cup \Gamma_{\text{out}}.
\end{array}
\right. 
\quad (23) 
$$

Finally, the gradient of the cost function is

$$ 
\left\{ 
\begin{array}{l}
\nabla_a \mathcal{J}_{\text{ram}} = \partial_a \bar{J}_{\Omega} + \partial_a h_{\tau}(\alpha) \bar{u} \cdot \bar{u}^* + C \partial_a k_{\tau}(\alpha) \nabla \bar{\theta} \cdot \nabla \bar{\theta}^* - \partial_a f(\alpha) B \bar{\theta} \bar{e}_y \cdot \bar{u}^* \quad \text{on } \Omega, \\
\nabla_a \mathcal{J}_{\text{ram}} = \partial_a \bar{J}_{\Gamma} - C \partial_a k_{\tau}(\alpha) \phi \bar{\theta}^* \quad \text{on } \Gamma_{\text{w}}.
\end{array}
\right. 
\quad (24) 
$$

where $(\bar{u}^*, \bar{\theta}^*, \bar{p}^*)$ satisfies the continuous adjoint problem (21)–(23). To solve the RAM optimization problem so defined, we employ a gradient descent method to iteratively update $\alpha$. So, applying the proposed approach, our algorithm gives us an optimized design for which the functional is well minimized.

**Remark 5.** Since our aim is to minimize pressure losses, our algorithm will tend to promote a steady flow regime, which is known to stabilize the flow and reduce large temporal fluctuations. This behavior will be clearly seen in our simulations, particularly through the evolution of the time step $\Delta T$, which adjusts to flow stability and progressively decreases. This observation supports our decision to retain only the final value of the design parameter $\alpha$ because this final state yields an optimal configuration consistent with the selected criterion.

## 5. Validation cases

In this section, we evaluate the accuracy and performance of the proposed time-averaged method for optimizing unsteady laminar and turbulent flows. We examine three cases:

* **Case 1:** Unsteady two dimensional laminar flow around a circular cylinder with $Re = 140$ (see Fig. 2a)

* **Case 2:** The classic Pitz-Daily turbulence case [37] with $Re = 12700$ (cf. Fig. 2b)

* **Case 3:** The Pitz-Daily turbulence case with an imposed heat flux on the lower wall.

For each of these cases, we have represented the horizontal velocity component in Figs. 3(a),5(a), 8(a) before optimization process. Our observations confirm the presence of unstable flow in the two-dimensional case and turbulent flow in the three-dimensional cases. For Cases 1 & 2, we employ the classic pressure losses objective function [15, Eq. (22)] which is equivalent to [38, Section 3.2]:

$$ J = J_1 = - \int_0^T \int_{\Gamma} \left( p + \frac{1}{2} |\bar{u}|^2 \right) \bar{u} \cdot \bar{n} \, dt \, d\Gamma, \quad (25) $$

from which we obtain

$$ J_{\text{ram}} = - \int_{\Gamma} \left( \bar{p} + \frac{1}{2} |\bar{u}|^2 \right) \bar{u} \cdot \bar{n} \, d\Gamma. $$

For the Case 3, we combine two objective functions: the previously pressure drop minimization $J_1$ and the heat transfer maximization [25], equivalent to [38, Section 3.3] and defined as $J_2$:

$$ J_2 = \int_0^T \int_{\Gamma} \bar{n} \cdot \bar{u} \, \theta \, dt \, d\Gamma, \quad (26) $$

from which we obtain

$$ J_{\text{ram}} = - \beta_1 \int_{\Gamma} \left( \bar{p} + \frac{1}{2} |\bar{u}|^2 \right) \bar{u} \cdot \bar{n} \, d\Gamma - \beta_2 \int_{\Gamma} \bar{n} \cdot \bar{u} \, \bar{\theta} \, d\Gamma, $$

where $\beta_1 = \beta_2 = 5$ in all our numerical experiments.

For the three studied cases, $\nu_{\text{ram}}^{\text{E}}$ and $k_{\text{ram}}^{\text{E}}$ are those defined in equations (17)).
<table><caption>**Table 1**
Summary of parameters used in our numerical experiments.</caption><thead><tr><th></th><th>Case 1</th><th>Case 2</th><th>Case 3</th></tr></thead><tbody><tr><td>Reynolds number Re</td><td>140</td><td colspan="2">12700</td></tr><tr><td>Prandtl number Pr</td><td colspan="3">0.71</td></tr><tr><td>Richardson number Ri</td><td>0.338</td><td>1e-8</td><td>0.338</td></tr></tbody></table>

<table><thead><tr><th>α<sub>max</sub></th><th>α<sub>0</sub></th><th>τ</th><th>k<sub>s</sub><br>(Aluminum)</th><th>k<sub>f</sub><br>(Air)</th><th>tol<sub>u</sub></th><th>tol<sub>θ</sub></th><th>tol<sub>J</sub></th></tr></thead><tbody><tr><td>60</td><td>20</td><td>0.6</td><td>6.4e-5</td><td>2.125e-5</td><td>1e-4</td><td>1e-4</td><td>1e-2</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>**Fig. 2.** Three configuration cases.</figcaption></figure>

In addition, the next fluid-solid interpolation functions [29] are used

$$h_{\tau}(\alpha) = \alpha_{\max} \left( \frac{1}{1 + e^{-\tau(\alpha - \alpha_0)}} - \frac{1}{1 + e^{\tau\alpha_0}} \right),$$

$$k_{\tau}(\alpha) = \left( \frac{k_s}{k_f} - 1 \right) \left( \frac{1}{1 + e^{-\tau(\alpha - \alpha_0)}} - \frac{1}{1 + e^{\tau\alpha_0}} \right) + 1,$$

where $α_0$ represents the midpoint of the sigmoid function, $τ$ is the slope parameter, and $α_{max}$ is the maximum value that the design parameter $α$ can take. $k_s$ and $k_f$ are the thermal diffusivities of the solid and fluid domains, respectively.

Parameters used in the simulations are listed in Table 1.

Finally, Table 2 summarizes the results obtained using 3.1. We report the maximum number of internal iterations performed in Step 1. This corresponds to the number of loops to obtain $Δt$ for which the time variation of the averaged quantities ($\bar{u}$, $\bar{\theta}$, $\bar{p}$) do not vary over a given tolerance. This row actually indicates the number of time we solved the time-dependent PDE (Equation (1) supplemented with boundary conditions 2). We indicate the total number of updates of the design parameter $α$, corresponding to the number of optimization loops. This corresponds to the number of steady-state problems that are solved during the whole optimization process. Finally, we show the values of the cost function $J$ before and after optimization and the relative decrease of cost function.
<figure><img src="image_6.png" /><figcaption>(a) Case 1: Horizontal velocities without optimization at $t = 150$ in the transient simulation</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(b) Case 1: Horizontal velocities after optimization at time-step number 3 corresponding to $t = 16$ in the transient simulation</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(c) Case 1: Horizontal velocities after optimization at time-step number 21 corresponding to $t = 52$ in the transient simulation</figcaption></figure>

<caption>Fig. 3. Horizontal velocities obtained for Case 1 without and during optimization process.</caption>

## 5.1. Numerical methods

In this study, the governing and adjoint equations are discretized using the Finite Volume Method (FVM). Simulations are carried out with OpenFOAM v2206. The discretized equations (1)) are resolved through the use of the Pressure-Implicit with Splitting of Operators (PISO) algorithm. This corresponds to Step 1 of Algorithm 3.1. The Semi-Implicit Method for Pressure-Linked Equations (SIMPLE) algorithm and a segregated approach for the pressure-velocity coupling are implemented to resolve the discretized RAM equations, corresponding to Step 3 of Algorithm 3.1. We employed the Hestenes-Stiefel conjugate gradient method in our calculations when solving the adjoint RAM problem during Step 3, and used the Fletcher-Reeves method for the descent algorithm. The convective, diffusive and temporal terms are handled using a specific second-order scheme. To ensure oscillation-free solutions while maintaining
<figure><img src="image_9.png" /><figcaption>(a) Optimized design at time-step number 3 corresponding to t = 16 in the transient simulation</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(b) Optimized design at time-step number 21 corresponding to t = 51 in the transient simulation</figcaption></figure>

<caption>Fig. 4. Optimized design obtained for Case 1 during optimization process.</caption>

<caption>Table 2</caption>

Summary of the results obtained with the Algorithm 3.1.

<table><thead><tr><th></th><th>Case 1</th><th>Case 2</th><th>Case 3</th></tr></thead><tbody><tr><td>Max iterations of Step 1</td><td>21</td><td>5</td><td>30</td></tr><tr><td>Updates of α in Step 3</td><td>21</td><td>211</td><td>71</td></tr><tr><td>J value before optimization</td><td>8.57e-1</td><td>1.5e-3</td><td>5.96e-3</td></tr><tr><td>J value after optimization</td><td>9.61e-3</td><td>1.33e-3</td><td>1.18e-3</td></tr><tr><td>Reduction rate of J</td><td>98.87%</td><td>11.33%</td><td>80.20%</td></tr></tbody></table>

second-order accuracy, Total Variation Diminishing (TVD) schemes are employed for computing the convective term in the time-dependent equations.

For the turbulent cases (Case 2 & 3), we used, instead of (1), a Large Eddy Simulation model described by Yoshizawa [39], considering buoyancy effects. In addition, we employed a subgrid scale Prandtl number of 0.6 [40]. This value is chosen to model the turbulent heat flux. The RAM model obtained from the implementation of the Yoshizawa model gave equations (18) and the Algorithm 3.1 can be applied to these turbulent cases without any modifications. Furthermore, following the approach of several authors such as Lazarov [41], a topology filter is applied to ensure a homogeneous transition between fluid and solid regions. This approach consisted in applying a Helmholtz-type partial differential equation as an alternative to standard density filtering in topology optimization problems. In all our numerical experiments, the filter radius is set to 0.05. These numerical methods ensure accurate and stable solutions for complex turbulent flows while keeping computational costs reasonable.

## 5.2. Unsteady 2D laminar flow

The flow field around the cylinder is modelled in two dimensions with the axes of the cylinder perpendicular to the direction of flow. The cylinder is modelled as a circle and a rectangular flow domain is created surrounding the cylinder. We consider the left-to-right flow of a cylinder of diameter D = 2R, immersed in an incompressible fluid. The computational domain consists of an upstream 13 times the radius to downstream 100 times the radius and the width of the domain is 13 times the radius of the cylinder. The problem configuration and important dimensions are shown in Fig. 2(a). A typical hexahedra mesh was used for the numerical simulations of this study. This particular mesh has 300,862 nodes, 597,707 faces and 148,536 cells with considerable mesh concentration around the cylinder. A rectangle of 3 times the radius of the cylinder to 6 times the radius of the cylinder is created around the cylinder with smaller meshing in regions where there is more movement, including boundary layers, to ensure more accurate results. The wall boundary conditions used in this case are those of non-slip condition namely $\vec{u} = (u, v) = (0, 0)$. At the inlet, we set $\vec{u} = U\vec{e}_x$ with $U = 0.105$ corresponding to $Re = 140$. Also the flow exit is treated as a pressure outlet (see equation (2)). The simulation was performed with a time duration of $T = 400$.

Fig. 4 shows optimized design obtained after minimizing the pressure losses cost function (equation (25)) using Algorithm 3.1. The optimized design is computed when the unsteady flow is fully developed as shown in Fig. 3(a). The algorithm adds material
<figure><img src="image_11.png" /><figcaption>(a) Case 2: Horizontal velocities without optimization at time t = 524 in the transient simulation</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(b) Case 2: Horizontal velocities during optimization at time-step number 151 corresponding to t = 200 in the transient simulation</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(c) Case 2: Horizontal velocities during optimization at time-step number 271 corresponding to t = 320 in the transient simulation</figcaption></figure>

<caption>**Fig. 5.** Horizontal velocities for Case 2 without optimization and during optimization process.</caption>

just behind the cylinder in order to reduce pressure losses, resulting in a shape similar to a well-known NACA wing profile. Fig. 10 shows the convergence history of the cost function $J$ over the iterations of the Algorithm 3.1. Thus, the optimized shape is obtained after 21 iterations. By the way, oscillations of horizontal velocity observed in Fig. 3(a) were suppressed by Algorithm 3.1 which added material behind the cylinder (see Fig. 3(b)-(c)). As expected, to minimize pressure losses, our algorithm tries to favor a steady laminar flow, known to stabilize flow and reduce significant temporal variations. This behavior justifies our choice of keeping the final value of the design parameter $\alpha$, assuming that the final state provides an optimal configuration representative of the criterion under consideration.

We can conclude first that the algorithm succeeds in reducing the pressure losses of 98.9% (Table 2). The flow becomes steady as shown in Fig. 3 thanks to the optimized design obtained in Fig. 4. Calculations related to our Algorithm 3.1 are faster than those of other methods used in the literature [15,17,19,21,26]. Indeed, in our method, the adjoint RAM problem is solved only once after the convergence of each average quantity. Particularly for this studied case, convergence of average quantities, such as mean velocities, is reached also at most in 21 iterations of the PISO algorithm (see Step 1 of Algorithm 3.1). It should be noted that during the optimization process, Step 1 of the algorithm takes less and less time, since as the optimization progresses, the flow stabilizes and,
<figure><img src="image_14.png" /><figcaption>(a) Case 2: Optimized design at time-step number 151 corresponding to $t = 200$ in the transient simulation</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(b) Case 2: Optimized design at time-step number 271 corresponding to $t = 320$ in the transient simulation</figcaption></figure>

<caption>Fig. 6. Optimized design obtained for Case 2 during optimization process.</caption>

consequently, so do the average quantities. In this Case 1, we also update the design parameter $\alpha$ (cf. Step 3 in Algorithm 3.1) 21 times.

### 5.3. Turbulent flow

This studied case is inspired from the experiment performed in 1983 by Robert W Pitz and John W Daily [42] in which they aimed to assess the effect of combustion on the mean flow field properties. Le et al. [43] also studied turbulent flow over a backward-facing step by using direct numerical solution of the Navier-Stokes equations. This is a simple test case to evaluate and verify the proposed method RAM when applied to unsteady turbulent flow. The simulation was done with a three-dimensional mesh consisting of a short inlet, a backward facing step, and a converging nozzle as the outlet as shown in Fig. 2b. The domain consists of an entry region of length 2D before the step with a width of 3.5 and an exit region of length 11.5D after the step, the entry region has a height of D as the step height. A typical hexahedral mesh was used for simulation of this study. This particular mesh has 799,068 nodes, 2,327,135 faces and 764,250 cells. The wall boundary conditions used in this case are those of non-slip condition i.e., $u = 0, v = 0$. A fixed value constraint with velocity $U = 1.0$ corresponding to $Re = 12700$ is applied at the inlet boundary. Also the flow exit is treated as a pressure outlet.

Fig. 5(a) illustrates a turbulent flow with a recirculation zone below the corner of the step, identifiable by the negative values of the horizontal component of velocity, at the time $t = 524$ in the transient simulation. The optimized design obtained by our algorithm can be seen in Fig. 6 with the corresponding optimization histories, i.e. optimized design are displayed after a number of updates of $\alpha$ parameter (see Step 3 of our Algorithm 3.1). Once again, we see that, to limit pressure losses, our algorithm tends to add material to stabilize the flow and attenuate fluctuations over time. This behavior justifies our choice to retain the final value of the design parameter $\alpha$ since this final state represents the most satisfactory solution according to the chosen optimization criterion. As observed by Hammond et al. [44], we obtained an elongated ‘ramp’ length that suppresses recirculation regions. Finally, the cost function is reasonably converged after 211 iterations of our Algorithm 3.1. The convergence of mean physical quantities, corresponding to Step 1 in our Algorithm 3.1, i.e the PISO algorithm, is achieved in maximum 5 repetitions. The algorithm succeeds in reducing pressure losses by 11.3% in the case of the Pitz-Daily turbulent flow.

### 5.4. Turbulent anisothermal flow

For this specific studied case, we apply the wall boundary conditions detailed in the previous section, with the exception of the bottom wall where a constant heat flux of $\phi = 1$ is imposed and we prefer a domain width of 1 to lighthen the overall numerical computations. The mesh configuration remains unchanged, and we define a steady RAM model where turbulent quantities are substituted with estimated values $v_{ram}^E$ and $k_{ram}^E$ defined in Equation (17). First of all, we emphasize that this simulation is particularly longer than those of previous cases regarding the amount of time steps required to solve steady-state for the adjoint RAM model. Fig. 9(a) shows the temperature field before optimization. The optimized shape is achieved after 71 updates of the design parameter
<figure><img src="image_16.png" /><figcaption>(a) Case 3: Horizontal velocities before optimization at time t = 526 in the transient simulation</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(b) Case 3: Horizontal velocities after optimization at time-step number 21 corresponding to t = 100 in the transient simulation</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>(c) Case 3: Optimized design at time-step number 72 corresponding to t = 200 in the transient simulation</figcaption></figure>

<caption>**Fig. 7.** Horizontal velocities for Case 3 without optimization and during optimization process.</caption>

$\alpha$, corresponding to 71 iterations of Algorithm 3.1. The convergence of mean physical quantities, corresponding to Step 1 in our Algorithm 3.1, is achieved in approximately 30 repetitions. This optimized design successfully reduces pressure losses and enhances thermal transfer efficiency by 80.2%. As in the other two cases, the material added by the algorithm stabilizes the flow (Fig. 7) and reduces significant temporal variations. This allows us to retain the last $\alpha$ obtained as the final solution to our optimization problem. Thus, Fig. 8(c) illustrates the effectiveness of our algorithm in creating a design that gives a steady laminar flow. We can therefore conclude that the RAM method reduces pressure losses and increases heat transfer by 80.2% in the case of turbulent anisothermal flow.
<figure><img src="image_19.png" /><figcaption>(a) Case 3: Optimized design at time-step number 21 corresponding to t = 100 in the transient simulation</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>(b) Case 3: Optimized design at time-step number 72 corresponding to t = 200 in the transient simulation</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>(a) Case 3: Temperature field without optimization at time-step number corresponding to t = 526 in the transient simulation</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>(b) Case 3: Temperature field after optimization at time-step number corresponding to t = 200 in the transient simulation</figcaption></figure>

<caption>Fig. 9. Case 3: Temperature field for Case 3 (a) before optimization - (b) after optimization.</caption>
<figure><img src="image_23.png" /><figcaption>Fig. 10. Evolution of cost function *J* during iterations of Algorithm 3.1 beyond convergence.</figcaption></figure>

## 6. Conclusion

In this paper, we have proposed a new algorithm based on a Reynolds-Averaged Method to solve topology optimization problem for buoyancy-driven turbulent flows and for laminar unsteady flows. The incompressible Navier-Stokes equations or Large Eddy Simulation and averaged physical variables, such as mean velocities, are computed over a series of iterations, counted between 5 and 30, using the PISO (Pressure-Implicit with Splitting of Operators) algorithm. Once these averaged quantities stabilize according to given criteria on their mean value over a specified time interval, a RAM model is implemented, in which the turbulent kinetic viscosity and turbulent thermal diffusivity are replaced by some estimated values derived from the averages of the physical variables. This RAM model is then used as a constraint in an optimization problem whose solution allow us to update the design variable $\alpha$. The process is repeated until the cost function reaches a local minima. Our algorithm has been illustrated with three configuration cases: an unsteady two dimensional laminar flow around a circular cylinder with a Reynolds number of $Re = 140$, the classical Pitz-Daily turbulence case where $Re = 12700$ and an anisothermal turbulent flow by adapting the previous Pitz-Daily case. We demonstrate the algorithm capability to reduce pressure losses between 11.3% and 98.9%. The design parameter $\alpha$ is updated between 21 times to 211 times in order to achieve the optimized design. Another remarkable feature of our algorithm is its ability to design an optimal shape so that the flow to become stable. Thus, our algorithm overcomes three major scientific challenges in solving optimization problems in turbulence: leveraging efficient temporal turbulence models or a DNS model, computational cost and data storage requirements.
We emphasize that the algorithm remains adapted to a specific objective function, namely minimizing pressure losses. In future work, it should be extended and tested on other objective functions in order to assess the flexibility and robustness of the algorithm over a wider range of cost functions.

## CRediT authorship contribution statement

**Delphine Ramalingom:** Writing – review & editing, Writing – original draft, Investigation, Conceptualization; **Alain Bastide:** Writing – review & editing, Validation, Supervision, Conceptualization; **Pierre-Henri Cocquet:** Writing – review & editing, Validation, Supervision; **Michaël Rakotobé:** Writing – original draft, Software, Methodology; **David Marti:** Writing – original draft, Software, Methodology.

## Data availability

Data will be made available on request.

## Declaration of competing interest

The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: Reports a relationship with that includes: Has patent pending to. If there are other authors, they declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Appendix A. Gradient computation with adjoint equations for the RAM model

In our time-averaged algorithm, we need to solve the following optimization problem

$$ \begin{array}{l} \min \mathcal{J}_{\text{ram}}(\alpha, \bar{u}, \bar{p}, \bar{\theta}) \\ \text{s.t. } (\bar{u}, \bar{p}, \bar{\theta}) \text{ solution of (18) – (19) parametrized by } \alpha, \end{array} \qquad (A.1) $$

where $\mathcal{J}_{\text{ram}}$ is a given cost function that may correspond to the cost function constrained by the RAM model described in this paper. We chose to solve (A.1) with a *differentiate then discretize* approach [36] which allows to compute the gradient of the cost function with respect to $\alpha$ thanks to the adjoint model. To compute the latter, we first introduce the Lagrangian associated to (A.1) which reads

$$ \begin{aligned} \mathcal{L}(\vec{X}, \vec{X}^*) = & \mathcal{J}_{\text{ram}}(\alpha, \bar{u}, \bar{\theta}, \bar{p}) - \int_{\Omega} \bar{p}^* \nabla \cdot \bar{u} d\Omega \\ & - \int_{\Omega} (-\nabla \bar{p} + 2\nabla \cdot (\{A + v_{\text{ram}}^E(\vec{X})\} S(\bar{u})) - h_{\tau}(\alpha)\bar{u} + f(\alpha)B\bar{\theta}\bar{c}_y) \cdot \bar{u}^* d\Omega \\ & - \int_{\Omega} (\nabla \cdot (\{C k_{\tau}(\alpha) + k_{\text{ram}}^E(\vec{X})\} \nabla \bar{\theta})) \bar{\theta}^* d\Omega \\ & - \int_{\Gamma_{\text{in}}} (\bar{u} - \bar{u}_{\text{in}}) \cdot \bar{\Phi}_{\text{in}} + \bar{\theta} \bar{\psi}_{\text{in}} d\Gamma - \int_{\Gamma_{\text{w}}} (\bar{u} \cdot \bar{\Phi}_{\text{w}} + \partial_n \bar{\theta} \bar{\psi}_{\text{w}}) d\Gamma \\ & - \int_{\Gamma_{\text{d}}} (\bar{u} \cdot \bar{\Phi}_{\text{d}} + \partial_n \bar{\theta} \bar{\psi}_{\text{d}}) d\Gamma \\ & - \int_{\Gamma_{\text{out}}} (\partial_n \bar{u} \cdot \bar{\Phi}_{\text{out}} + \partial_n \bar{\theta} \bar{\psi}_{\text{out}} + \bar{p} \bar{q}_{\text{out}}) d\Gamma \\ & + \int_{\Omega} (\nabla \cdot (\bar{u} \otimes \bar{u})) \cdot \bar{u}^* d\Omega \\ & + \int_{\Omega} (\nabla \cdot (\bar{\theta} \bar{u})) \bar{\theta}^* d\Omega \end{aligned} $$

where $\vec{X} = (\bar{u}, \bar{\theta}, \bar{p})$ and $\vec{X}^* = (\bar{u}^*, \bar{\theta}^*, \bar{p}^*, \mathcal{A}^*)$ are respectively the primal and adjoint variables. We remind that $A = \text{Re}^{-1}$, $B = \text{Ri}$ and $C = \text{Re}^{-1} \text{Pr}^{-1}$. The additional adjoint variables

$$ \mathcal{A}^* = (\bar{\Phi}_{\text{in}}, \bar{\psi}_{\text{in}}, \bar{\Phi}_{\text{w}}, \bar{\psi}_{\text{w}}, \bar{\Phi}_{\text{d}}, \bar{\psi}_{\text{d}}, \bar{\Phi}_{\text{out}}, \bar{\psi}_{\text{out}}, \bar{q}_{\text{out}}), $$

are only introduced to enforce the boundary conditions and are not going to be involved in the definition of the adjoint model (see e.g. [29, Remark 1]). We set $\partial_X F(X)[\delta X]$ as the Fréchet derivative of $F$ which is defined as

$$ \partial_X F(X)[\delta X] = \lim_{h \to 0} \frac{F(X + h\delta X) - F(X)}{h}. $$

One can then see that the constraints equations (18)-(19) are recovered with

$$ \partial_{\vec{X}^*} \mathcal{L}(\vec{X}, \vec{X}^*)[\delta \vec{X}^*] = 0, \forall \delta \vec{X}^*, $$
while the so-called adjoint model is given as

$$ \partial_{\vec{X}} \mathcal{L}(\vec{X}, \vec{X}^*)[\delta \vec{X}] = 0, \forall \delta \vec{X}. $$

Noting $\overline{\vec{X}, \vec{X}^*}$ the solution to the primal and adjoint equations respectively, the gradient of the cost function with respect to $\alpha$ can be finally computed with

$$ \partial_a \mathcal{L}(\overline{\vec{X}, \vec{X}^*})[\delta \alpha] = 0, \forall \delta \alpha. $$

We refer to [36] for additional details and now move to the computation of the adjoint model. It is worth noting that the adjoint model have been computed in [29,38] in the case where the estimated turbulent parameters $v_{ram}^E, k_{ram}^E$ vanishes. A computation gives:

$$
\begin{aligned}
\partial_{\vec{X}} \mathcal{L}(\vec{X}, \vec{X}^*)[\delta \vec{X}] &= \partial_{\vec{X}} \mathcal{J}_{ram}(\alpha; \vec{X})[\delta \vec{X}] - \int_{\Omega} \bar{p}^* \nabla \cdot \delta \bar{u} \, d\Omega \\
&\quad - \int_{\Omega} (-\nabla \delta \bar{p} + 2\nabla \cdot (\{A + v_{ram}^E(\vec{X})\} S(\delta \bar{u})) - h_{\tau}(\alpha) \delta \bar{u} + f(\alpha) B \delta \bar{\theta} \bar{e}_y) \cdot \bar{u}^* \, d\Omega \\
&\quad - \int_{\Omega} (\nabla \cdot (\{C k_{\tau}(\alpha) + k_{ram}^E(\vec{X})\} \nabla \delta \bar{\theta})) \bar{\theta}^* \, d\Omega \\
&\quad - \int_{\Omega} (\nabla \cdot (\partial_{\vec{X}} k_{ram}^E(\vec{X})[\delta \vec{X}] \nabla \bar{\theta})) \bar{\theta}^* + (\nabla \cdot \partial_{\vec{X}} v_{ram}^E(\vec{X})[\delta \vec{X}] S(\bar{u})) \cdot \bar{u}^* \, d\Omega \\
&\quad - \int_{\Gamma_{in}} \delta \bar{u} \cdot \bar{\Phi}_{in} + \delta \bar{\theta} \bar{\psi}_{in} \, d\Gamma - \int_{\Gamma_{w}} (\delta \bar{u} \cdot \bar{\Phi}_{w} + \delta \bar{\theta} \bar{\psi}_{w}) \, d\Gamma \\
&\quad - \int_{\Gamma_{d}} (\delta \bar{u} \cdot \bar{\Phi}_{d} + \partial_{n} \delta \bar{\theta} \bar{\psi}_{d}) \, d\Gamma \\
&\quad - \int_{\Gamma_{out}} (\partial_{n} \delta \bar{u} \cdot \bar{\Phi}_{out} + \partial_{n} \delta \bar{\theta} \bar{\psi}_{out} + \delta \bar{p} \bar{q}_{out}) \, d\Gamma, \\
&\quad + \int_{\Omega} (\nabla \cdot (\delta \bar{u} \otimes \bar{u}) + \nabla \cdot (\bar{u} \otimes \delta \bar{u})) \cdot \bar{u}^* \, d\Omega \\
&\quad + \int_{\Omega} (\nabla \cdot (\delta \bar{\theta} \bar{u}) + \nabla \cdot (\bar{\theta} \delta \bar{u})) \bar{\theta}^* \, d\Omega
\end{aligned}
$$

and it now remains to integrate by parts to remove the derivatives on $\delta \vec{X}$ and next get the adjoint model. We recall first the next Green's formulas

$$
\begin{aligned}
\int_{\Omega} \nabla \cdot \mathcal{A} \cdot \bar{v} \, d\Omega &= \int_{\Gamma} \mathcal{A} \bar{n} \cdot \bar{v} \, d\sigma - \int_{\Omega} \mathcal{A} : \nabla \bar{v} \, d\Omega, \\
\int_{\Omega} \nabla \cdot (b S(\bar{u})) \cdot \bar{v} \, d\Omega &= \int_{\Gamma} b S(\bar{u}) \bar{n} \cdot \bar{v} - b \bar{u} \cdot S(\bar{v}) \bar{n} \, d\sigma \\
&\quad + \int_{\Omega} \bar{u} \cdot \nabla \cdot (b S(\bar{v})) \, d\Omega
\end{aligned}
$$

where $\mathcal{A}: \Omega \mapsto \mathbb{R}^{d \times d}$ is a smooth tensor field and $b$ a smooth scalar. Using these integration by parts formulas, we obtain

$$
\begin{aligned}
\partial_{\vec{X}} \mathcal{L}(\vec{X}, \vec{X}^*)[\delta \vec{X}] &= \partial_{\vec{X}} \mathcal{J}_{ram}(\alpha; \vec{X})[\delta \vec{X}] - \int_{\Gamma} \delta \bar{u} \cdot \bar{n} \bar{p}^* \, d\Omega + \int_{\Omega} \delta \bar{u} \cdot \nabla \bar{p}^* \, d\Omega \\
&\quad + \int_{\Gamma} \delta \bar{p} \bar{n} \cdot \bar{u}^* \, d\Gamma - \int_{\Omega} \delta \bar{p} \nabla \cdot \bar{u}^* \, d\Omega - \int_{\Omega} (-h_{\tau}(\alpha) \delta \bar{u} + f(\alpha) B \delta \bar{\theta} \bar{e}_y) \cdot \bar{u}^* \, d\Omega \\
&\quad - \int_{\Gamma} 2 \{A + v_{ram}^E(\vec{X})\} (S(\delta \bar{u})) \bar{n} \cdot \bar{u}^* - \delta \bar{u} \cdot S(\bar{u}^*) \bar{n} \, d\Gamma \\
&\quad - \int_{\Omega} \delta \bar{u} \cdot 2 \nabla \cdot (\{A + v_{ram}^E(\vec{X})\} S(\bar{u}^*)) \, d\Omega \\
&\quad - \int_{\Gamma} \{C k_{\tau}(\alpha) + k_{ram}^E(\vec{X})\} (\nabla \delta \bar{\theta} \cdot \bar{n} \bar{\theta}^* - \nabla \bar{\theta}^* \cdot \bar{n} \delta \bar{\theta}) \, d\Gamma \\
&\quad - \int_{\Omega} \delta \bar{\theta} (\nabla \cdot (\{C k_{\tau}(\alpha) + k_{ram}^E(\vec{X})\} \nabla \bar{\theta}^*)) \, d\Omega \\
&\quad - \int_{\Gamma} \partial_{\vec{X}} k_{ram}^E(\vec{X})[\delta \vec{X}] \nabla \bar{\theta} \cdot \bar{n} \bar{\theta}^* \, d\Gamma + \int_{\Omega} \partial_{\vec{X}} k_{ram}^E(\vec{X})[\delta \vec{X}] (\nabla \bar{\theta} \cdot \nabla \bar{\theta}^*) \, d\Omega \\
&\quad - 2 \int_{\Gamma} \partial_{\vec{X}} v_{ram}^E(\vec{X})[\delta \vec{X}] S(\bar{u}) \bar{n} \cdot \bar{u}^* \, d\Gamma + 2 \int_{\Omega} \partial_{\vec{X}} v_{ram}^E(\vec{X})[\delta \vec{X}] (S(\bar{u}) : \nabla \bar{u}^*) \, d\Omega \\
&\quad - \int_{\Gamma_{in}} \delta \bar{u} \cdot \bar{\Phi}_{in} + \delta \bar{\theta} \bar{\psi}_{in} \, d\Gamma - \int_{\Gamma_{w}} (\delta \bar{u} \cdot \bar{\Phi}_{w} + \delta \bar{\theta} \bar{\psi}_{w}) \, d\Gamma
\end{aligned}
$$
$$
\begin{aligned}
& - \int_{\Gamma_d} (\delta \bar{u} \cdot \bar{\Phi}_d + \partial_n \delta \bar{\theta} \bar{\psi}_d) d\Gamma \\
& - \int_{\Gamma_{\text{out}}} (\partial_n \delta \bar{u} \cdot \bar{\Phi}_{\text{out}} + \partial_n \delta \bar{\theta} \bar{\psi}_{\text{out}} + \delta \bar{p} \bar{q}_{\text{out}}) d\Gamma \\
& + \int_{\Gamma} (\delta \bar{u} \otimes \bar{u}) + (\bar{u} \otimes \delta \bar{u}) \bar{n} \cdot \bar{u}^* d\Gamma - \int_{\Omega} (\delta \bar{u} \otimes \bar{u}) + (\bar{u} \otimes \delta \bar{u}) : \nabla \bar{u}^* d\Omega \\
& + \int_{\Gamma} (\delta \bar{\theta} \bar{u} \cdot \bar{n}) + (\bar{\theta} \delta \bar{u} \cdot \bar{n}) \bar{\theta}^* d\Gamma - \int_{\Omega} (\delta \bar{\theta} \bar{u}) + (\bar{\theta} \delta \bar{u}) \cdot \nabla \bar{\theta}^* d\Omega.
\end{aligned}
$$

Assuming now that the cost function can be written as

$$
\mathcal{J}_{\text{ram}} = \int_{\Omega} \bar{J}_{\Omega}(\alpha; \bar{u}, \bar{\theta}, \bar{p}) d\Omega + \int_{\Gamma} \bar{J}_{\Gamma}(\alpha; \bar{u}, \bar{\theta}, \bar{p}) d\Gamma,
$$

the adjoint model on the computation domain Ω then reads

$$
\left\{
\begin{aligned}
& -\nabla \cdot \bar{u}^* + 2\partial_{\bar{p}} v_{\text{ram}}^{\text{E}} (S(\bar{u}): \nabla \bar{u}^*) + \partial_{\bar{p}} k_{\text{ram}}^{\text{E}} (\nabla \bar{\theta} \cdot \nabla \bar{\theta}^*) + \partial_{\bar{p}} \bar{J}_{\Omega} = 0, \\
& \nabla \bar{p}^* + h_{\tau}(\alpha) \bar{u}^* - 2\nabla \cdot \left( \left\{ A + v_{\text{ram}}^{\text{E}}(\bar{X}) \right\} S(\bar{u}^*) \right) - \bar{\theta} \nabla \bar{\theta}^* - 2S(\bar{u}^*) \bar{u} \\
& \qquad + 2\partial_{\bar{u}} v_{\text{ram}}^{\text{E}} (S(\bar{u}): \nabla \bar{u}^*) + \partial_{\bar{u}} k_{\text{ram}}^{\text{E}} (\nabla \bar{\theta} \cdot \nabla \bar{\theta}^*) + \partial_{\bar{u}} \bar{J}_{\Omega} = 0, \\
& -\left( \nabla \cdot \left\{ C k_{\tau}(\alpha) + k_{\text{ram}}^{\text{E}} \nabla \bar{\theta}^* \right\} \right) - \bar{u} \cdot \nabla \bar{\theta}^* - f(\alpha) B \bar{e}_y \cdot \bar{u}^* \\
& \qquad + 2\partial_{\bar{\theta}} v_{\text{ram}}^{\text{E}} (S(\bar{u}): \nabla \bar{u}^*) + \partial_{\bar{\theta}} k_{\text{ram}}^{\text{E}} (\nabla \bar{\theta} \cdot \nabla \bar{\theta}^*) + \partial_{\bar{\theta}} \bar{J}_{\Omega} = 0,
\end{aligned}
\right.
\quad (21)
$$

where we emphasize that $\partial_{\bar{u}} v_{\text{ram}}^{\text{E}}(\bar{X})$ and $\partial_{\bar{u}} k_{\text{ram}}^{\text{E}}(\bar{X})$ are vectors. In addition, it is worth noting that we recover the incompressibility condition $\nabla \cdot \bar{u}^* = 0$ when $\partial_{\bar{p}} v_{\text{ram}}^{\text{E}} = 0$, $\partial_{\bar{p}} k_{\text{ram}}^{\text{E}} = 0$ and $\partial_{\bar{p}} J_{\text{ram}} = 0$. For the boundary conditions on the adjoint velocity $\bar{u}$, we have

$$
\left\{
\begin{array}{ll}
\bar{n} \cdot \bar{u}^* = -\partial_{\bar{p}} \bar{J}_{\Gamma}, & \bar{t} \cdot \bar{u}^* = 0 \quad \text{on } \Gamma_{\text{in}} \cup \Gamma_{\text{w}} \cup \Gamma_{\text{d}}, \\
-\bar{n} \bar{p}^* + 2 \left\{ A + v_{\text{ram}}^{\text{E}}(\bar{X}) \right\} S(\bar{u}^*) \bar{n} \\
\qquad + \bar{n} \bar{\theta} \bar{\theta}^* + (\bar{u}^* \otimes \bar{n}) \bar{u} + (\bar{n} \otimes \bar{u}) \bar{u}^* \\
\qquad -\partial_{\bar{u}} k_{\text{ram}}^{\text{E}} \nabla \bar{\theta} \cdot \bar{n} \bar{\theta}^* - 2 \partial_{\bar{u}} v_{\text{ram}}^{\text{E}} S(\bar{u}) \bar{n} \cdot \bar{u}^* = -\partial_{\bar{u}} \bar{J}_{\Gamma} & \text{on } \Gamma_{\text{out}},
\end{array}
\right.
\quad (22)
$$

where $\bar{t}$ is some unit tangent vector. For the adjoint time-averaged temperature $\bar{\theta}^*$, we have the next boundary conditions, in particular with a constant heat flux imposed on $\Gamma_w$

$$
\left\{
\begin{array}{ll}
\bar{\theta}^* = 0 & \text{on } \Gamma_{\text{in}}, \\
\left\{ \begin{array}{l} C k_{\tau}(\alpha) + k_{\text{ram}}^{\text{E}}(\bar{X}) \nabla \bar{\theta}^* \cdot \bar{n} \\ + \bar{\theta}^* \bar{u} \cdot \bar{n} \end{array} \right\} & \text{on } \Gamma_{\text{in}} \cup \Gamma_{\text{w}} \cup \Gamma_{\text{d}}, \\
-\partial_{\bar{\theta}} k_{\text{ram}}^{\text{E}} \nabla \bar{\theta} \cdot \bar{n} \bar{\theta}^* - 2 \partial_{\bar{\theta}} v_{\text{ram}}^{\text{E}} S(\bar{u}) \bar{n} \cdot \bar{u}^* = -\partial_{\bar{\theta}} \bar{J}_{\Gamma} & \text{on } \Gamma_{\text{w}} \cup \Gamma_{\text{d}} \cup \Gamma_{\text{out}}.
\end{array}
\right.
\quad (23)
$$

Finally, the gradient of $J_{\text{ram}}$ is defined thanks to $\partial_{\alpha} \mathcal{L}(\bar{X}, \bar{X}^*)[\delta\alpha] = 0, \forall \delta\alpha$. The gradient of the cost functional then reads

$$
\left\{
\begin{array}{ll}
\nabla_{\alpha} \mathcal{J}_{\text{ram}} = \partial_{\alpha} \bar{J}_{\Omega} + \partial_{\alpha} h_{\tau}(\alpha) \bar{u} \cdot \bar{u}^* + C \partial_{\alpha} k_{\tau}(\alpha) \nabla \bar{\theta} \cdot \nabla \bar{\theta}^* & \\
-\partial_{\alpha} f(\alpha) \bar{B} \bar{\theta} \bar{e}_y \cdot \bar{u}^* & \text{on } \Omega, \\
\nabla_{\alpha} \mathcal{J}_{\text{ram}} = \partial_{\alpha} \bar{J}_{\Gamma} - C \partial_{\alpha} k_{\tau}(\alpha) \phi \bar{\theta}^* & \text{on } \Gamma_{\text{w}},
\end{array}
\right.
\quad (24)
$$

where $(\bar{u}^*, \bar{\theta}^*, \bar{p}^*)$ is the solution to the continuous adjoint problem (21)–(23).

**Remark 6.** The adjoint model (21)-(22) is only defined if $v_{\text{ram}}^{\text{E}}$, $k_{\text{ram}}^{\text{E}}$ do not depend on $\partial_{\bar{u}} \bar{u}$, $\partial_{\bar{u}} \bar{\theta}$. In such cases, one have to integrate by part further the terms involving $\partial_{\bar{X}} v_{\text{t}}^{\text{E}}(\bar{X})[\delta \bar{X}]$, $\partial_{\bar{X}} \lambda_{\text{t}}^{\text{E}}(\bar{X})[\delta \bar{X}]$ to remove all derivatives on $\delta \bar{X}$. In our numerical results, we will not consider estimated turbulent quantities of this kind and thus we are not going to compute the adjoint model in this extended setting enven if this does not yield additional difficulties.
**Remark 7.** In this paper, we deal with the case where we impose a constant heat flux on the heated wall. In the case where the boundary condition used is a constant temperature, we obtain this term in the Lagrangian expression:

$$ \int_{\Gamma_w} (\bar{\bar{u}} \cdot \bar{\bar{\Phi}}_w + (\bar{\bar{\theta}} - \theta_h) \bar{\bar{\psi}}_w) d\Gamma $$

where $\theta_h$ is a given constant temperature.

As the same, for the adjoint time-averaged temperature $\bar{\bar{\theta}}^*$, we would have the next boundary conditions instead of (23)

$$ \left\{ \begin{array}{ll} \bar{\bar{\theta}}^* = 0 & \text{on } \Gamma_{in} \cup \Gamma_w, \\ \left\{ C k_{\varepsilon}(\alpha) + k_{\text{ram}}^E(\vec{X}) \right\} \nabla \bar{\bar{\theta}}^* \cdot \vec{n} \\ + \bar{\bar{\theta}}^* \bar{\bar{u}} \cdot \vec{n} & \\ - \partial_{\bar{\bar{\theta}}} k_{\text{ram}}^E \nabla \bar{\bar{\theta}} \cdot \vec{n} \bar{\bar{\theta}}^* - 2 \partial_{\bar{\bar{\theta}}} v_{\text{ram}}^E S(\bar{\bar{u}}) \vec{n} \cdot \bar{\bar{u}}^* = - \partial_{\bar{\bar{\theta}}} \bar{\bar{J}}_{\Gamma} & \text{on } \Gamma_d \cup \Gamma_{out}. \end{array} \right. $$

Then, the gradient of the cost functional on $\Gamma_w$ in (24) would be

$$ \nabla_{\alpha} J_{\text{ram}} = \partial_{\alpha} \bar{\bar{J}}_{\Gamma} \quad \text{on } \Gamma_w $$

## References

[1] J. Petersson, T. Borrvall, Topology optimization of fluids in Stokes flow, Int J Numer Meth Fl 41 (2003) 77–107.

[2] J. Alexandersen, N. Aage, C.S. Andreasen, O. Sigmund, Topology optimisation for natural convection problems, Int. J. Numer. Methods Fluids 76 (10) (2014) 699–721.

[3] J. Alexandersen, O. Sigmund, N. Aage, Large scale three-dimensional topology optimisation of heat sinks cooled by natural convection, Int. J. Heat Mass Transf. 100 (2016) 876–891.

[4] T.E. Bruns, Topology optimization of convection-dominated, steady-state heat transfer problems, Int. J. Heat Mass Transf. 50 (15–16) (2007) 2859–2873.

[5] A. Vieira, A. Bastide, P.-H. Cocquet, Topology optimization for steady-state anisothermal flow targeting solids with piecewise constant thermal diffusivity, Appl. Math. Optim. 85 (3) (2022) 41.

[6] A. Vieira, P.-H. Cocquet, The bousinesq system with nonsmooth boundary conditions: existence, relaxation, and topology optimization, SIAM J. Control Optim. 60 (6) (2022) 3457–3484.

[7] Y.H. Choi, G.H. Yoon, A new density filter for pipes for fluid topology optimization, J. Fluid Mech. 986 (2024) A9.

[8] F. Navah, M.-É. Lamarche-Gagnon, F. Ilinca, Thermofluid topology optimization for cooling channel design, Appl. Therm. Eng. 236 (2024) 121317.

[9] A. Evgrafov, The limits of porous materials in the topology optimization of Stokes flows, Appl. Math. Optim. 52 (3) (2005) 263–277.

[10] H. Garcke, M. Hinze, C. Kahle, K.F. Lam, A phase field approach to shape optimization in Navier–Stokes flow with integral state constraints, Adv Comput Math 44 (2018) 1345–1383.

[11] J. Alexandersen, C.S. Andreasen, A review of topology optimization for fluid-based problems, Fluids 5 (1) (2020) 29. https://doi.org/10.3390/fluids5010029

[12] C.D. Argyropoulos, N.C. Markatos, Recent advances on the numerical modelling of turbulent flows, Appl. Math. Model. 39 (2) (2015) 693–732.

[13] D.C. Wilcox, et al., Turbulence Modeling for CFD, 2, DCW Industries La Canada, CA, 1998.

[14] S. Kreissl, G. Pingen, K. Maute, Topology optimization for unsteady flow, Int. J. Numer. Methods Eng. 87 (13) (2011) 1229–1253.

[15] Y. Deng, Z. Liu, P. Zhang, Y. Liu, Y. Wu, Topology optimization of unsteady incompressible Navier–Stokes flows, J. Comput. Phys. 230 (17) (2011) 6688–6708. https://doi.org/10.1016/j.jcp.2011.05.004

[16] Y. Deng, Z. Liu, Y. Wu, Topology optimization of steady and unsteady incompressible Navier–Stokes flows driven by body forces, Struct. Multidiscip. Optim. 47 (4) (2012) 555–570. https://doi.org/10.1007/s00158-012-0847-8

[17] N.K. Yamaleev, B. Diskin, E.J. Nielsen, Local-in-time adjoint-based method for design optimization of unsteady flows, J. Comput. Phys. 229 (14) (2010) 5394–5407. https://doi.org/10.1016/j.jcp.2010.03.045

[18] K. Yaji, M. Ogino, C. Chen, K. Fujita, Large-scale topology optimization incorporating local-in-time adjoint-based method for unsteady thermal-fluid problem, Struct. Multidiscip. Optim. 58 (2) (2018) 817–822. https://doi.org/10.1007/s00158-018-1922-6

[19] C. Chen, K. Yaji, T. Yamada, K. Izui, S. Nishiwaki, Local-in-time adjoint-based topology optimization of unsteady fluid flows using the lattice Boltzmann method, Mech. Eng. J. 4 (3) (2017) 17–00120.

[20] M. Theulings, R. Maas, L. Noël, F. van Keulen, M. Langelaar, Reducing time and memory requirements in topology optimization of transient problems, Int. J. Numer. Methods Eng. 125 (14) (2024) e7461.

[21] H. Nobis, P. Schlatter, E. Wadbro, M. Berggren, D.S. Henningson, Topology optimization of unsteady flows using the spectral element method, Comput. Fluids 239 (2022) 105387.

[22] C. Othmer, A continuous adjoint formulation for the computation of topological and surface sensitivities of ducted flows, Int. J. Numer. Methods Fluids 58 (8) (2008) 861–877.

[23] M. Pietropaoli, R. Ahlfeld, F. Montomoli, A. Ciani, M. D’Ercole, Design for additive manufacturing: internal channel optimization, J. Eng. Gas Turbine Power 139 (10) (2017).

[24] A.S. Zymaris, D.I. Papadimitriou, K.C. Giannakoglou, C. Othmer, Continuous adjoint approach to the Spalart–Allmaras turbulence model for incompressible flows, Comput. Fluids 38 (8) (2009) 1528–1538.

[25] E.A. Kontoleontos, E.M. Papoutsis-Kiachagias, A.S. Zymaris, D.I. Papadimitriou, K. Giannakoglou, Adjoint-based constrained topology optimization for viscous flows, including heat transfer, Eng. Optim. 45 (8) (2013) 941–961.

[26] C.B. Dilgen, S.B. Dilgen, D.R. Fuhrman, O. Sigmund, B.S. Lazarov, Topology optimization of turbulent flows, Comput. Methods Appl. Mech. Eng. 331 (2018) 363–393.

[27] G.H. Yoon, Topology optimization for turbulent flow with Spalart–Allmaras model, Comput. Methods Appl. Mech. Eng. 303 (2016) 288–311.

[28] G.H. Yoon, Topology optimization method with finite elements based on the k-ε turbulence model, Comput. Methods Appl. Mech. Eng. 361 (2020) 112784.

[29] D. Ramalingom, P.-H. Cocquet, A. Bastide, A new interpolation technique to deal with fluid-porous media interfaces for topology optimization of heat transfer, Comput. Fluids 168 (2018) 144–158. https://doi.org/10.1016/j.complfluid.2018.04.005

[30] O. Reynolds, On the Dynamical Theory of Incompressible Viscous Fluids and the Determination of the Criterion (Reprinted from Papers on Mechanical and Physical Subjects, Vol 2, Pg 535-577, 1901), in: Proceedings of the Royal Society-Mathematical and Physical Sciences, 451, Royal Soc 6–9 Carlton House Terrace, London SW1Y 5AG, England, 1995, pp. 5–47.

[31] J. Bousinesq, Essai sur la théorie des eaux courantes, Imprimerie Nationale, Paris, France, 1877.

[32] A.N. Kolmogorov, Equations of turbulent motion in an incompressible fluid, in: Dokl. Akad. Nauk SSSR, 30, 1941, pp. 299–303.
[33] P.Y. Chou, On velocity correlations and the solutions of the equations of turbulent fluctuation, Q. Appl. Math. 3 (1) (1945) 38–54. https://doi.org/10.1090/qam/11999

[34] J. Rotta, Statistische Theorie nichthomogener Turbulenz, Zeitschrift für Physik 129 (6) (1951) 547–572. https://doi.org/10.1007/bf01330059

[35] S. Zeierman, M. Wolfshtein, Turbulent time scale for turbulent-flow calculations, AIAA J. 24 (10) (1986) 1606–1610. https://doi.org/10.2514/3.9490

[36] M.D. Gunzburger, Perspectives in Flow Control and Optimization, SIAM, 2002.

[37] H.K. Versteeg, An Introduction to Computational Fluid Dynamics the Finite Volume Method, 2/E, Pearson Education India, 2007.

[38] D. Ramalingom, P.-H. Cocquet, R. Maleck, A. Bastide, A multi-objective optimization problem in mixed and natural convection for a vertical channel asymmetrically heated, Struct. Multidiscip. Optim. 60 (5) (2019) 2001–2020. https://doi.org/10.1007/s00158-019-02306-7

[39] A. Yoshizawa, Statistical theory for compressible turbulent shear flows, with the application to subgrid modeling, Phys. Fluids 29 (7) (1986) 2152-2164.
https://doi.org/10.1063/1.865552

[40] D. Li, Revisiting the subgrid-scale Prandtl number for large-eddy simulation, J. Fluid Mech. 802 (2016) R2.

[41] B.S. Lazarov, O. Sigmund, Filters in topology optimization based on Helmholtz-type differential equations, Int. J. Numer. Methods Eng. 86 (6) (2010) 765-781.
https://doi.org/10.1002/nme.3072

[42] R.W. Pitz, J.W. Daily, Combustion in a turbulent mixing layer formed at a rearward-facing step, AIAA J. 21 (11) (1983) 1565–1570.

[43] H. Le, P. Moin, J. Kim, Direct numerical simulation of turbulent flow over a backward-facing step, J. Fluid Mech. 330 (1997) 349–374.

[44] J. Hammond, M. Pietropaoli, F. Montomoli, Topology optimisation of turbulent flow using data-driven modelling, Struct. Multidiscip. Optim. 65 (2) (2022) 49.