# Simulation of turbulent flow subjected to conjugate heat transfer via a dual immersed boundary method

G.F. Narváez, Eric Lamballais, E.B. Schettini

► To cite this version:

G.F. Narváez, Eric Lamballais, E.B. Schettini. Simulation of turbulent flow subjected to conjugate heat transfer via a dual immersed boundary method. Computers and Fluids, 2021, 229, pp.105101. 10.1016/j.compfluid.2021.105101. hal-04217800

HAL Id: hal-04217800

https://hal.science/hal-04217800v1

Submitted on 17 Mar 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Simulation of turbulent flow subjected to conjugate
heat transfer *via* a dual immersed boundary method

G. F. Narváezª, E. Lamballaisᵇ, E. B. Schettiniª

ªInstituto de Pesquisas Hidráulicas (IPH), Universidade Federal do Rio Grande do Sul,
Porto Alegre, RS, Brazil

ᵇPPRIME Institute, Incompressible Turbulence and Control Group - Université de Poitiers,
CNRS, ISAE-ENSMA Futuroscope Chasseneuil Cedex France

Abstract

In the context of turbulent flows, this study presents a new methodology for
the high-fidelity simulation of conjugate heat transfer between solid and fluid.
The numerical strategy is mainly based on the immersed boundary method to
ensure the expected thermal boundary condition at the fluid/solid interface.
The approach is dual since the solid domain is the immersed region for the fluid
and, conversely, the fluid for the solid. The resulting method is an extension
of an already existing immersed boundary method based on the reconstruc-
tion of the solution inside the immersed region. The duality of the technique
requires to define one temperature field for the fluid and another one for the
solid. The interaction between the two temperature fields is ensured, at the
fluid/solid interface, through an efficient weak coupling preserving numerical
stability and accuracy. Thanks to the reconstruction in its own immersed zone,
each temperature field solution is perceived as smooth by the numerical dif-
ferentiation method, for any fluid-to-solid ratio of thermal conductivities and
diffusivities, while representing a desirable sharp interface. This feature is cru-
cial for the high-order finite difference schemes used in this work. The method
is validated for the plane channel and pipe flow configurations and a demanding
pure conduction case with variable conductivity in a composite wall. Its ability
to preserve second-order accuracy in space is checked in the laminar and pure
conduction cases by comparison with analytical solutions. In the turbulent case, the assessment is based on comparisons of basic turbulent statistics with reference data obtained by direct numerical simulation. The step-by-step analysis of thermal statistics shows that the present immersed boundary method can ensure accurately imposed temperature, imposed heat flux, and conjugate heat transfer conditions at the fluid/solid interface.

*Keywords:* Turbulent channel flow · Turbulent pipe flow · Conjugate heat transfer · Direct numerical simulation · Immersed boundary method ·

High-order numerical schemes

# 1. INTRODUCTION

A moving fluid is an efficient medium to transport heat from one source to another. In the industry, especially in heating/cooling devices, flows are common in the turbulent regime. Thus, the efficient design of these devices requires a better understanding of turbulent heat transfer processes for modelling purposes. In many applications, the influence of temperature fluctuations on the fluid dynamics can be neglected, leading to the general framework of forced convection in which the temperature can be considered as a passive scalar. The present study is based on this framework while focusing on thermal fluid/solid interaction occurring in wall-bounded turbulent flows subjected to temperature gradients. This particular interaction, usually referred to as *Conjugate Heat Transfer (CHT)*, corresponds to an unsteady problem where turbulent heat transfers in the fluid and heat conduction in the solid are fully coupled.

Accurate representation of fluid/solid *CHT* is fundamental for the prediction of heat transfer in heating/cooling systems, especially in extreme conditions, as in the context of the nuclear industry, for instance. Despite the fully recognised importance of *CHT* mechanisms in industrial applications, their physical
understanding suffers from a certain deficit which is against their relevant mod-
elling. The main lack of knowledge is related to phenomena occurring in the
neighbourhood of the fluid/solid interface, namely inside the thermal boundary
layer. Former papers have studied analytically and experimentally the CHT
problem (e.g. POLIAKOV, 1974; KHABAKHPASHEVA, 1984; SINAI, 1987; KASAGI
et al., 1989; SOMMER et al., 1995). However, the simplifications imposed by the
analytical approaches and the difficulty to measure experimentally turbulent
quantities inside very thin boundary layers (SAHA et al., 2014) can only give
a partial picture of the physical processes associated with CHT. This situation
makes highly attractive numerical simulation as a tool to accurately investi-
gate the role of turbulent structures near the fluid/solid interface (KASAGI and
SHIKAZONO, 1995; DUCHAINE et al., 2009, 2013) in CHT. As a prototype flow
involving turbulent heat transfer, the turbulent heated channel flow has been
investigated by Direct Numerical Simulations (DNS) (e.g. KIM (1989); KASAGI
et al. (1992); KAWAMURA et al. (1998); TISELJ et al. (2001a,b); PIROZZOLI et al.
(2016); ORLANDI et al. (2016); FLAGEUL et al. (2015, 2017)), as well as, tur-
bulent heated pipe flow (e.g. SATAKE et al., 2000; REDJEM-SAAD et al., 2007;
SAHA et al., 2011). In most of these studies, CHT is idealized using simplified
boundary conditions at the fluid/solid interface consisting in (i) locally imposed
temperature (**IT**), (ii) locally imposed heat flux (**IF**) or (iii) combination of **IT**
and **IF**, as a Robin condition. These conditions require the *a priori* knowledge
of the temperature and/or heat flux at the wall. In practical situations involv-
ing CHT, these conditions are rarely known, especially when thermal fluid/solid
diffusivities of the solid are of the same order (FLAGEUL et al., 2015), leading
to unsatisfactory predictions of heat transfer.

Fluid and solid energy equations have to be coupled at the fluid/solid in-
terface to numerically solve the *CHT* problem. A first strategy is to solve
each equation separately in its domain while matching their common boundary
conditions to ensure their physical coupling. Technically, this coupling can be
strong, using, for instance, an iterative process (WANG and GEORGIADI, 1996),
or can be weak, by tolerating a small temperature discontinuity at the interface
(YOUNG and VAFAI, 1998; TISELJ et al., 2001a; FLAGEUL et al., 2015, 2017).
In this approach, distinct computational domains have to be employed, provid-
ing flexibility for the choice of numerical methods (which can be different in
fluid/solid regions) but also complexity for the treatment of complex geometry.
A second strategy is to include all the fluid/solid regions in a single compu-
tational domain while ensuring the switching from fluid to solid temperature
equation. This switching is incorporated through a change of the diffusion coef-
ficients (viscosity or thermal conductivity) and the introduction of a volumetric
source term, depending on the solid/fluid regions considered. This alternative
essentially corresponds to an *Immersed boundary method (IBM)* initially intro-
duced by PESKIN (1972), for the mass and momentum balance equations, but
this time applied to the thermal energy equation (KIM and CHOI (2004)). For
a detailed review on *IBM*, the reader is referred to MITTAL and IACCARINO
(2005). A well-known advantage of *IBM* is the flexibility of the computational
mesh generation, which is not required to be body-fitted, allowing, for instance,
to consider moving complex solid geometries using a rudimentary stationary
Cartesian grid for the spatial discretization. Another convenient feature is the
possibility to represent both no-slip and thermal boundary conditions (KIM
and CHOI, 2004; PACHECO et al., 2005; IACCARINO and MOREAU, 2006; ZHANG
et al., 2008; WANG et al., 2009; KANG et al., 2009; REN et al., 2012, 2013;
LUO et al., 2017). Most of the previous *IBM* studies only consider *IT* on the
immersed boundary because this situation corresponds to a Dirichlet boundary
condition, similarly to the no-slip condition for the velocity. The implementa-
tion of *IF* using *IBM* represents a bigger challenge because it corresponds to a Neumann boundary condition for which formally, the forcing is not applied to the solution but to its spatial derivatives (REN et al., 2013). This point is addressed in the first part of this study, as a preliminary step towards the *CHT* imposition.

In the context of *IBM*, some strategies have been proposed to deal with the *CHT*. SONG et al. (2017) have developed and evaluated for laminar cases an *IBM* based on the effective thermal conductivity concept with a modification of the convection terms. For turbulent *CHT* in the finite volume framework, some *RANS-IBM* (e.g. IACCARINO and MOREAU (2006), DE MARINIS et al. (2016)), *LES-IBM* (e.g. HE (2019) and LEE et al. (2019)) and *DNS-IBM* (DAS et al. (2019), FAVRE et al. (2019)) solvers have been developed. On the other hand, LUO et al. (2017) developed a *ghost-cell IBM* in explicit finite differences to simulate compressible flows under temperature boundary conditions of Dirichlet, Neumann, and Robin type. This methodology is applied to 2D external flows without including *CHT* (at low Reynolds numbers or Mach number close to one). To the best of our knowledge, in the compact finite difference framework, neither *DNS* nor *LES* study considering turbulent *CHT* by means of *IBM* has been documented yet.

This paper presents the development of an original, straightforward, and accurate *IBM* combined with high-order compact finite differences schemes to study *CHT* in wall-bounded turbulent flows by *DNS/LES*. The highly accurate and massively parallel code *Incompact3d* (LAIZET and LAMBALLAIS, 2009) is used while implementing the new *IBM* for *CHT*, in order to preserve its massively parallel features. For that purpose, the customised *IBM* of GAUTIER et al. (2014) is extended to the temperature treatment by the implementation of *IT* and *IF* conditions up to the realistic case of *CHT*. Since the basic mechanisms
of the fluid/solid heat transfer have already been investigated through DNS in turbulent channel flow (TISELJ et al., 2001a; FLAGEUL et al., 2015, 2017), the focus is this canonical configuration. The pipe flow is also considered under IT condition, as a required step for CHT modelling while having data of reference for this more demanding geometry with respect to the Cartesian mesh used in this study. In the next section, the governing equations, boundary conditions, flow configuration and parameters are presented as well as the connections between input and global output data such as the friction Reynolds and Nusselt numbers. Section 3 is dedicated to the algorithm presentation, including the IBM to impose Dirichlet or Neumann temperature conditions at the wall and the fluid/solid coupling strategy. The developments for laminar and turbulent channel and pipe flows are assessed in Section 4. Main concluding remarks are summarised in Section 5.

## 2. GOVERNING EQUATIONS AND FLOW CONFIGURATION

In this section, the governing equations corresponding to fluid/solid thermal interaction are presented, in the framework of incompressible flows representing the temperature as a passive scalar. Attention is paid to the boundary conditions at the walls and in the streamwise direction for which the periodicity requires to impose a specific forcing term. We also discuss the heat transfer equations to represent the fluid/solid thermal coupling.

### 2.1. *Fluid dynamics equations*

In this study, $x$ is the streamwise direction, whereas $yz$ is a cross-section of the channel/pipe (Figure 1). In channel configuration $y$ and $z$ are the wall-normal and spanwise directions, respectively. The fluid flow in the channel develops in $-H/2 \le y \le H/2$, where $H$ is the channel height. In pipe configuration, the fluid flow develops in $r \le D/2$, where $D$ is the internal pipe diameter
and $r = \sqrt{y^2 + z^2}$ the radial coordinate. To express governing equations in di-

<figure><img src="image_2.png" /><figcaption>Figure 1: Channel and pipe flow configurations.</figcaption></figure>

mensionless form, the reference quantities employed are the bulk velocity $U_b$, the channel height $H$ or the pipe diameter $D$, and the density of the fluid $\rho$. Incompressible Navier-Stokes dimensionless equations read

$$ \frac{\partial \mathbf{u}}{\partial t} = -\nabla p - \frac{1}{2} [\nabla (\mathbf{u} \cdot \mathbf{u}) + (\mathbf{u} \cdot \nabla)\mathbf{u}] + \frac{1}{Re} \nabla^2 \mathbf{u} + \mathbf{f}_{\tau}, \quad (1) $$

$$ \nabla \cdot \mathbf{u} = 0, \quad (2) $$

where $\mathbf{u} = [u_x, u_y, u_z]^T$ is the velocity field and $p$ the pressure field. The fluid is assumed to have constant density $\rho$ and kinematic viscosity $\nu$. The streamwise forcing term $\mathbf{f}_{\tau} = [f_{\tau}, 0, 0]^T$, constant in space but variable in time, is imposed to compensate the net wall viscous friction, maintaining constant in time the bulk velocity $U_b$, defined as

$$ U_b = \frac{1}{H} \int_{-H/2}^{+H/2} \langle u_x \rangle dy, \quad \text{in channel,} \quad (3) $$

$$ U_b = \frac{8}{D^2} \int_{0}^{D/2} \langle u_x \rangle r dr, \quad \text{in pipe,} \quad (4) $$

where $\langle \cdot \rangle$ represents time/longitudinal and spanwise (channel) or azimuthal (pipe) average, then $\langle u_x \rangle$ depends only on $y$ (channel) or $r$ (pipe). This approach corresponds to the classical treatment of streamwise periodicity in the
simulation of channels and ducts. The global Reynolds number $Re$ is defined using the bulk velocity and channel height or pipe diameter like $Re = U_bH/\nu$ or $Re = U_bD/\nu$, respectively. Note that in our IBM, no forcing term has to be introduced in (1) using a technique where boundary conditions are implicitly considered for the computation of the spatial derivatives in the governing equations.

For the full computational domain, periodic boundary conditions are imposed in every direction, whereas solid walls are modelled by the IBM. A specific forcing is applied to ensure both the velocity and temperature boundary conditions at the fluid/solid interface (Section 3.2).

## 2.2. Heat transfer equations

The fluid temperature is assumed to be a passive scalar subjected to the convection-diffusion equation

$$ \partial_t T = -\mathbf{u} \cdot \nabla T + \alpha \nabla^2 T, \qquad (5) $$

where $\alpha$ is the constant thermal diffusivity of the fluid related to the kinematic viscosity by the Prandtl number $Pr = \nu/\alpha$. The thermal conductivity $\lambda$ is related to the thermal diffusivity $\alpha$ by $\alpha = \frac{\lambda}{\rho c_p}$ where $c_p$ is the specific heat capacity assumed constant in the fluid.

For the present flow configurations, heat transfer at the wall makes the temperature field statistically non-homogeneous in the streamwise direction with a non-zero streamwise gradient of the bulk $T_b(x)$ and mean-wall $T_w(x)$ temperatures of the fluid defined, in channel, as

$$ T_b(x) = \frac{1}{HU_b} \int_{-H/2}^{+H/2} \langle u_x T \rangle_{t,z} dy, \qquad (6) $$

$$ T_w(x) = \langle T \rangle_{y=\pm H/2} \qquad (7) $$
and, in pipe, as

$$T_b(x) = \frac{8}{D^2 U_b} \int_0^{D/2} \langle u_x T \rangle_{t,\theta} r dr, \quad (8)$$

$$T_w(x) = \langle T \rangle_{r=D/2} \quad (9)$$

where, $\langle \cdot \rangle_{t,z}$ represents time/spanwise average, then $\langle u_x T \rangle_{t,z}$ depends only on $x$ and $y$; $\langle \cdot \rangle_{t,\theta}$ represents time/azimuthal average, then $\langle u_x T \rangle_{t,\theta}$ depends only on $x$ and $r$. This means that, in a fully developed heated channel/pipe flow, the mean temperature profile is expected to evolve along the streamwise direction $x$.

To make compatible the streamwise dependency of temperature with the use of periodic boundary conditions in $x$, it is common to introduce a *scaled* dimensionless temperature $\Theta$, defined to ensure statistical homogeneity in the streamwise direction with a constant bulk temperature $\Theta_b = \text{cst}$, defined as in (6 and 8). This dimensionless temperature has a different expression depending on the type of wall boundary conditions for temperature, as explained in the two next sections. For simplicity, in the following dimensionless definitions, we refer to $H$ as a generic characteristic length, which can be seen as the channel height (as already defined) or as the pipe diameter ($H \to D$).

### 2.3. Imposed temperature boundary condition

Using the dimensionless temperature

$$\Theta = \frac{T_w - T}{T_w - T_b}, \quad (10)$$

where $T_w$ is the imposed temperature at the wall, it is easy to show that the bulk dimensionless temperature is constant (with $\Theta_b = 1$), as a condition to make statistically homogeneous $\Theta$ in the streamwise direction while enabling the
imposition of periodic boundary condition in this direction as for the velocity.
Then, imposed temperature boundary conditions (IT) simply consist in

$$ \Theta_w = 0 \qquad (11) $$

while leaving free the wall-normal gradient $\partial_n \Theta|_w$. This corresponds to the mixed-type condition described by PILLER (2005).

Using the wall-normal heat flux $\varphi_w = \lambda \partial_n T|_w$, a Nusselt number can be defined as $Nu = \frac{H \varphi_w}{\lambda(T_w - T_b)}$ and simply related to $\Theta$ by $Nu = -H \partial_n \Theta|_w$ with an associated mean value $\langle Nu \rangle = -H \partial_n \langle \Theta \rangle|_w$. The resulting boundary condition is mixed in the sense that it is neither isothermal (because of the streamwise evolution of $T_w$) nor isoflux (because of the wall fluctuations of $\varphi_w$) as described by PILLER (2005). However, it has in common with isothermal boundary conditions that no turbulent fluctuations are possible at the wall.

The governing equation for $\Theta$ can be easily obtained by substitution of (10) into (5) leading to

$$ \partial_t \Theta = -\mathbf{u} \cdot \nabla \Theta + \alpha \nabla^2 \Theta + f_{\Theta}, \qquad (12) $$

where the source term $f_{\Theta}$ is proportional to $u_x$, such that $\Theta_b = 1$ at any time.

## 2.4. Imposed heat flux boundary condition

For the case of an isoflux boundary condition (IF), the convenient scaled dimensionless temperature becomes

$$ \Theta = -\frac{T - T_b}{Nu \frac{T_w - T_b}{T_w}} \qquad (13) $$

leading again to a constant bulk dimensionless temperature which is, in this case, $\Theta_b = 0$. With this definition, the imposition of the wall-normal gradient
temperature $\partial_n T|_w$ corresponds to

$$ \partial_n \Theta|_w = \frac{1}{H}, \qquad (14) $$

while leaving free the wall values $\Theta_w$. The mean Nusselt number is connected to its mean value at the wall by means of $\langle Nu \rangle = 1/\langle \Theta_w \rangle$.

For the IF thermal condition, the governing equation for $\Theta$ is formally identical to (12) with a source term $f_\Theta$ remaining proportional to $u_x$, but with a different adjusted magnitude in order to ensure $\Theta_b = 0$ at any time.

For IF and IT thermal conditions, the computational domain is assumed periodic in all directions for the computation of spatial derivatives of the corresponding scaled dimensionless temperature $\Theta$ and for the velocity dimensionless $u$. The immersed plates representing the channel walls have a thickness of $W = H/20$, enabling a smooth reconstruction of the solution, as explained in section 3.3. A buffer zone is also used to consider the plate as fully immersed with a fluid at rest in the two subdomains $-L_y/2 \le y \le -H/2 - W$ and $H/2 + W \le y \le L_y/2$ (see the schematic view in Figure 2-top). Analogously, the pipe wall has a thickness $W = D/20$ with a buffer zone in the subdomain $r \ge D/2 + W$ (see cross-flow projection in Figure 3). No physical information is expected in these immersed and buffer zones, which are only used for convenience in the present IBM technique.

## 2.5. Conjugate heat transfer boundary condition

Conjugate heat transfer (CHT) represents the fluid/solid thermal coupling. Thus, the temperature equation in the fluid domain $\Omega$ and the conduction equation in the solid domain $\Omega_s$ are coupled by the thermal boundary conditions at the fluid/solid interface $\Omega \cap \Omega_s$ and the outer solid boundary. The dimensionless
temperature equations read

$$ \partial_t \Theta = -\mathbf{u} \cdot \nabla \Theta + \alpha \nabla^2 \Theta + f_\Theta \quad \text{in } \Omega, \quad (15) $$

$$ \partial_t \Theta_s = \frac{\alpha}{G} \nabla^2 \Theta_s \quad \text{in } \Omega_s, \quad (16) $$

$$ \Theta = \Theta_s \quad \text{at } \Omega \cap \Omega_s, \quad (17) $$

$$ \partial_{\mathbf{n}} \Theta_s = \frac{1}{G_2} \partial_{\mathbf{n}} \Theta \quad \text{at } \Omega \cap \Omega_s, \quad (18) $$

$$ \partial_{\mathbf{n}} \Theta_s = \frac{\varphi_w}{G_2} \quad \text{at outer } \Omega_s \text{ boundary,} \quad (19) $$

where $\Theta_s$ is the solid temperature, $G = \alpha/\alpha_s$ is the fluid-to-solid ratio of thermal diffusivities and $G_2 = \lambda_s/\lambda$ is the solid-to-fluid ratio of thermal conductivities.

For the CHT case, the solid regions correspond to the two subdomains $-L_y/2 \le y \le -H/2$ and $H/2 \le y \le L_y/2$, inside which the solid temperature is physically meaningful. The corresponding set of thermal boundary conditions is schematically summarised in Figure 2-bottom. The way to impose thermal boundary conditions at the fluid/solid interfaces $y = \pm H/2$ is explained in section 3.5, where it will be shown that technically, it requires to deal with IT conditions for the fluid and IF conditions for the solid. At the outer faces of the thick solid plates, corresponding to the upper and lower boundaries of the computational domain $y = \pm L_y/2$, a thermal isoflux boundary condition is applied to ensure consistency with the definition of the scaled dimensionless temperature $\Theta$, as defined in (13) and (14).

In pipe flow, only IT condition is assessed. This is an intermediate step towards the full CHT implementation (which requires an IT condition in the fluid) for a geometry disconnected from the computational mesh organisation. The limitation to IT condition enables us to make comparison with DNS reference data of turbulent pipe, a flow configuration for which CHT conditions have not been addressed yet, as far as we know.
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /><figcaption>Figure 2: Scheme of the domain configuration and the thermal boundary conditions at the domain and immersed boundaries. Top: imposed temperature (IT) or heat flux (IF); bottom: Conjugate heat transfer (CHT) conditions.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 3: *zy* projection of immersed pipe.</figcaption></figure>

# 3. NUMERICAL METHODOLOGY

## 3.1. *Numerical treatment of the governing equations*

The code *Incompact3d* is used, in order to solve numerically the fluid dynamic governing equations (1, 2). In this work, the code was adapted to solve the heat transfer governing equations (15, 16), employing the *IBM*. This open-source code (License GNU GPL v3, https://www.incompact3d.com/) was developed at the Université de Poitiers and Imperial College of London (LARDEAU et al., 2002; LAIZET and LAMBALLAIS, 2009; LAIZET and LI, 2011). The code is based on 6<sup>th</sup>-order *compact schemes of finite differences* for the spatial differentiation in a Cartesian mesh. A 3<sup>rd</sup>-order *Adams-Bashforth scheme* is adopted for the time advancement.

To treat the incompressibility condition, a fractional step method is employed, requiring to solve a Poisson equation for the pressure. Solving this equation directly in the spectral space using 3D Fast Fourier Transforms compatible with ghost boundary conditions, and thanks to the concept of modified wavenumber, the divergence-free condition is ensured up to the machine accuracy (LAIZET and LAMBALLAIS, 2009). Pressure is solved on a staggered grid with respect to the collocated mesh nodes of velocity components and temper-
ature.

In this code, a complex wall geometry is represented through a customised
IBM, which is an adaptation of a direct forcing method (GAUTIER et al., 2014).
The adaptation consists in reconstructing the solution inside the immersed mesh
to improve its smoothness (in terms of derivative continuity) at the fluid/solid
interface, allowing better performance of the high-order finite difference schemes
used in *Incompact3d*.

For more details on the code *Incompact3d* for the fluid dynamic resolution,
the reader is referred to LAIZET and LAMBALLAIS (2009). In the rest of this
section, the heat transfer treatment is detailed as the main novelty of this study
in terms of numerical development.

3.2. *Immersed boundary treatment*

As already mentioned, the presence of solid walls, in this study, is exclusively
modelled using an *IBM*. For the plane geometry, it is easy to ensure body-fitted
conditions using a simple Cartesian mesh as in the code *Incompact3d*. However,
in order to assess the potential of *IBM* to deal with *CHT*, we consider the turbu-
lent channel flow in this framework with the advantage that reliable data have
been documented by FLAGEUL et al. (2015, 2017) for this specific configuration.
This type of assessment is a crucial first step before going further toward com-
plex geometry for which, in the context of *CHT*, the availability of *DNS* data is
very recent and limited. It also enables straightforward access to the wall-normal
information so that the reconstruction technique of the fluid/solid solution can
be presented and validated with accuracy. The same strategy of *IBM* is also
validated for the pipe flow but only for the fluid through the application of an
*IT* condition.

For the velocity field, the main purpose of the IBM is to ensure no-slip condi-
tion (as a Dirichlet boundary condition) at the fluid/solid interface. In a direct
forcing framework, GAUTIER et al. (2014) developed a reconstruction method by
alternating direction, which can reduce spurious oscillations in the neighbour-
hood of the body while enabling to deal with complex geometries through the
preservation of the massively parallel features of the code *Incompact3d*. Initially,
the principle of this technique was to use information from the fluid domain to
define the solution inside the body region using Lagrange interpolating poly-
nomials to improve the accuracy of spatial differentiation. This interpolation
process will be referred to as *reconstruction*. In the present context of *CHT*,
this principle will be extended to the thermal field by using information from
the fluid domain to define the temperature in the solid domain, and conversely.
The reciprocity of the reconstruction will require to introduce two temperature
fields defined everywhere in the computational domain, one for the fluid and
another one for the solid.

3.3. Imposition of Dirichlet boundary condition at the fluid/solid interface

In *CHT*, the temperature continuity (17) corresponds to a Dirichlet Boundary
ary Condition (*DBC*) at the fluid/solid interface, similarly to the simpler *IT*
condition (11). Then, thermal *DBC* can be imposed using the same technique
as for no-slip condition of velocity, as explained in GAUTIER et al. (2014). The
principle of the temperature reconstruction is illustrated in Figure 4, for the
general case in which the solid boundary does not intercept a mesh node. In
the illustration of Figure 4, reconstruction is performed inside the solid domain
as an extension of the temperature in the fluid domain. The reconstruction is
essentially 1D and based on 5<sup>th</sup>-order Lagrange interpolating polynomials re-
quiring the use of temperature values at 4 mesh nodes inside the fluid domain
and also 2 boundary points at the fluid/solid interface where the *DBC* is im-
posed. At this polynomial order, it can be shown that a third-order accuracy
can be reached for the computation of first and second spatial derivatives. It
is well known that high-order Lagrange interpolation, as for the 5<sup>th</sup>-order here,
can lead to Runge's phenomenon with the occurrence of small-scale oscillations.
It would be against the smoothness expectation of the present reconstruction
technique. In practice, we have never observed small-scale (i.e. at scale not
accurately described by the mesh) oscillations in the IB region. One possible
reason for this good behaviour is that the reference points are not equispaced,
as a favourable condition to avoid Runge's phenomenon.

The location of the fluid/solid interface and the wall temperature Θ<sub>w</sub> are
used by the reconstruction, as a *soft* procedure to impose the *DBC*; where the
term soft means that no forcing is applied in the fluid region, but at the interface
location just as the sharp physical transition from fluid to solid. The closest fluid
mesh node to the fluid/solid interface is never used by the reconstruction to avoid
redundancy or singularity when it is almost located on the fluid/solid interface.
The temperature field is then defined everywhere in the computational domain
(i.e., in the solid as well as in the fluid) so that the spatial differentiation can
be performed in the reconstruction direction without any particular treatment
near the fluid/solid interface. The smoothness of the reconstruction improves
the accuracy of the differentiation scheme, especially for the high-order compact
finite difference schemes used in the code *Incompact3d*. In a 3D case, this
reconstruction scheme must be performed three times, namely, once before the
computation of all derivatives, in every direction. Even if only a simple geometry
is considered here, it is worth mentioning that this *IBM* strategy has been
shown to be efficient and robust for complex geometry without any additional
adaptation thanks to its 1D nature (GAUTIER et al., 2014). This conclusion
will be confirmed here by the *DNS* of a pipe flow with heat transfer, as a
prototype of complex geometry for which the wall geometry is disconnected
from the computational mesh organisation.
<figure><img src="image_6.png" /><figcaption>Figure 4: Illustration of 1D reconstruction scheme of the temperature $\Theta$ inside the solid domain using information from the fluid domain. Example for a 5$^{th}$-order Lagrange interpolating polynomial using 4 mesh nodes in the fluid domain and 2 points at the fluid/solid interface ($\otimes$).</figcaption></figure>

### 3.4. Imposition of Neumann boundary condition at the fluid/solid interface

In *CHT*, the heat flux continuity, as expressed in (18), corresponds to a Neumann Boundary Condition (NBC) at the fluid/solid interface, similarly to the simpler *IF* condition (14). As it will be explained in the next section, this *NBC* is used to solve the temperature field in the solid region, through the present coupling technique to treat *CHT*. For simplicity, and in coherence with the convention used in the previous section, the reconstruction example inside the solid region is also used for *IF* condition. In *CHT* case, the reconstruction inside the fluid region from the solid region is strictly identical, and it can be performed through a simple switch between fluid and solid, for the used and reconstructed information.

The imposition of *NBC* by an *IBM* is not an easy task. In this work, a straightforward strategy is proposed for a case where the wall normal corresponds to a main direction of the mesh, so that the reconstruction remains a 1D procedure as in GAUTIER et al. (2014). This wall-normal direction is referred as *y* by convention. Instead of trying to force the wall-normal derivative, it is suggested to only adjust the wall temperature $\Theta_w$ required to obtain the expected gradient. This approach, which consists of a soft conversion of *NBC* into
DBC through an approximation, is quite common, especially in the context of finite difference schemes outside the IBM framework. In particular, it has been used successfully by DAIRAY (2013) to perform DNS of an impinging jet on a heated plate subjected to IF condition using the code Incompact3d. Here, it is suggested to adapt this technique to the present IBM for which $\Theta_w$ is weakly imposed through the reconstruction scheme explained in the previous section.

Then, for a given heat flux $\partial_y \Theta|_w$, the temperature at the fluid/solid interface $\Theta_w$ can be estimated by a 3rd-order non-centred scheme

$$ \partial_y \Theta|_w = a \Theta|_w + b \Theta_{j+1} + c \Theta_{j+2} + d \Theta_{j+3}, \quad (20) $$

leading to the 2nd-order approximation

$$ \Theta|_w = \frac{-b\Theta_{j+1} - c\Theta_{j+2} - d\Theta_{j+3} + \partial_y \Theta|_w}{a}, \quad (21) $$

where,

$$ a = -b - c - d, $$

$$ b = \frac{1}{\Delta y_{int}} + \frac{\Delta y_{int} + \Delta y_1 + \Delta y_2}{\Delta y_1 \Delta y_2} - \frac{\Delta y_{int} + \Delta y_1}{\Delta y_2 (\Delta y_1 + \Delta y_2)}, $$

$$ c = -\Delta y_{int} \frac{\Delta y_{int} + \Delta y_1 + \Delta y_2}{\Delta y_1 \Delta y_2 (\Delta y_{int} + \Delta y_2)}, $$

$$ d = \Delta y_{int} \frac{\Delta y_{int} + \Delta y_1}{\Delta y_2 (\Delta y_{int} + \Delta y_1 + \Delta y_2) (\Delta y_1 + \Delta y_2)}. $$

Figure 5 illustrates the resulting estimation of $\Theta_w$ from the known fluid temperature and the expected value $\partial_y \Theta|_w$. Once this value of $\Theta_w$ is estimated, the same reconstruction scheme, as for a DBC, can be used as detailed in the previous section and illustrated in Figure 4. Note that in coherence with the
reconstruction scheme and to avoid singular cases, the closest mesh node from
the solid/fluid interface (corresponding to the node $j$) has no contribution in
the approximation scheme (21) as also illustrated in Figure 5.

<figure><img src="image_7.png" /><figcaption>Figure 5: Approximation of $\Theta_w$ from the knowledge of $\partial_y\Theta|_w$ and from the temperature values at three mesh nodes inside the fluid close to the fluid/solid interface. Scheme of conversion of a Neumann boundary condition into a Dirichlet boundary condition in order to impose the heat flux (i. e. wall-normal first derivative) to the fluid at the fluid/solid interface.</figcaption></figure>

The present imposition technique of NBC is particularly soft in the sense
that it is based on a conversion into a DBC, which is itself only weakly prescribed
in the present IBM. That is why its behaviour must be carefully assessed in
the laminar case as well as in the more demanding turbulent case. These two
flow regimes are considered here for the channel flow with the application of
Neumann BC at $y = \pm L_y/2$ for both IF and CHT conditions.

3.5. *Fluid/solid coupling*

In a general *IBM* framework, a single computational domain is used to cover
both solid and fluid regions. The solid is said to be immersed in the fluid, and
for a non-moving body, the main unknowns are related to the fluid. That is why
the fluid solution stored at the mesh nodes inside the solid, without physical
meaning, can be considered as available for various techniques of interpola-
tion/reconstruction to improve the accuracy of the boundary condition at the
fluid/solid interface.
In the *CHT* problem, the temperature field is defined in both fluid and solid regions, namely everywhere in the computational domain. This can be seen as an advantage of *IBM* in the sense that a single field can represent the temperature solution in the whole computational domain. However, this optimal storage of the solution has two drawbacks. First, the values of the temperature at the mesh nodes inside the solid, which are, in this case, physically meaningful, cannot be modified through a reconstruction procedure to improve the imposition of the boundary condition at the fluid/solid interface. This constraint makes the present reconstruction technique less convenient. The second drawback is more fundamental because it lies in the regularity of the temperature field at the fluid/solid interface. In many practical situations, conductivities are different in the fluid and the solid. As a consequence, because of the heat flux continuity (18), the wall-normal derivative of the temperature becomes discontinuous at the solid/fluid interface making only $C^0$ a temperature field defined everywhere in the computational domain. This lack of smoothness is clearly a disadvantage of an *IBM* approach based on a single field because of the difficulty of numerically dealing with a loss of differentiability.

Here, to preserve the efficiency of our IBM while enabling a straightforward treatment of the wall-normal derivative discontinuity of the temperature at the fluid/solid interface, two temperature fields are defined: one for the fluid $\Theta$ and another one for the solid $\Theta_s$. Each field has its own physical domain (i.e., where its governing equation is solved) and its own reconstruction domain. As a consequence, spatial derivatives are computed independently for each field.

Naturally, this dual strategy requires an extra computational effort. In terms of storage, by comparison to the requirement for one temperature field in the fluid with a time advancement based on four levels, as for the present 3$^{rd}$-order Adams-Bashforth scheme (*AB3*), the cost may be seen as doubled. However, in
a very recent study, Lamballais et al. (2021) has proposed a new approach called the viscous filtering which has the potential to limit the storage to only one level for a purely conductive equation. Thanks to this new feature, the extra-cost becomes only +25% by comparison to an AB3 scheme for the temperature equation alone. Relatively to the whole storage requirement of the fluid solver, the actual extra-cost is less than +5%.

In terms of computational effort, thanks to the simplicity of the heat conduction equation, the extra-cost is also moderate. By comparison to the cost of the temperature equation in the fluid region, we only observe an increase of +50%. With respect to the total cost of the fluid/solid solver, the actual extra-cost related to the resolution of a second temperature equation is less than +10%.

These additional computational efforts can be considered as moderate in view of the flexibility and accuracy offered by the present dual approach, especially when the fluid and solid thermal conductivities are different.

In order to preserve computational efficiency, governing equations are solved using the weak coupling technique for the time advancement. This technique allows the sequential solving of the fluid dynamic equations (1-2) and of the thermal equations (15-16) at every time step, whereas a strong coupling would require to solve all these equations simultaneously through, for instance, an iterative procedure. A strong coupling represents a drastic increase of the computational effort, highly undesirable in the context of DNS/LES. On the other hand, because of its explicit nature in time, a weak coupling can introduce severe numerical stability limitations. Fortunately, GILES (1997) have shown that the restriction on the time step can be alleviated if Dirichlet/Neumann boundary conditions are applied for the fluid/solid, respectively. In particular, to produce the DNS/LES reference data used in this study for validation, FLAGEUL et al. (2015) have observed that the restriction on the time step from the CFL con-
dition was stronger in comparison with the stability condition from the weak thermal coupling. Besides, FLAGEUL et al. (2015) have found that the temperature discontinuity at the solid/fluid interface resulting from the weak coupling procedure was negligible.

To summarise, the weak coupling strategy proposed by GILES (1997), and used successfully by FLAGEUL et al. (2015) for CHT in a turbulent channel, is also employed in this study, through a sequential time advancement for the fluid and the solid with the following two main steps:

1. The fluid temperature at the fluid/solid interface is updated by the half-sum of the fluid and solid values at the previous time step with

$$ \Theta|_{w}^{n+1} = \frac{1}{2} (\Theta|_{w}^{n} + \Theta|_{w}^{n}). \qquad (22) $$

This update is used as a DBC to solve the fluid temperature equation (15) leading to the prediction of $\Theta^{n+1}$ everywhere in the fluid. The DBC is imposed in the framework of IBM with the reconstruction technique explained in section 3.3.

2. Using the knowledge of $\Theta^{n+1}$, the temperature derivative at the fluid/solid interface $\partial_n \Theta|_w^{n+1}$ is computed at the wall using scheme (20). This value is used to provide a NBC for the solid using condition (18) updated at the new time step by

$$ \partial_n \Theta_s|_w^{n+1} = \frac{1}{G_2} \partial_n \Theta|_w^{n+1}. \qquad (23) $$

This NBC is imposed, in the framework of IBM, with the reconstruction technique explained in section 3.4. By solving the conduction equation (16), $\Theta_s^{n+1}$ is predicted everywhere in the solid. Because of the lack of any iteration, a small temperature difference is introduced at the wall
between the fluid and solid predictions, with a first-order accuracy in
time for the corresponding error $\Theta|_{w}^{n+1} - \Theta_{s}|_{w}^{n+1} = O(\Delta t)$. This error is
controlled throughout the simulation thanks to the scheme (22) used at
the beginning of every time step.

Although the methodology has been developed for constant thermal prop-
erties, the present dual *IBM* approach can be adapted to problems involving
variable thermal physical properties. In this more general context, the tem-
perature continuity and energy conservation at the interface (Eq. 17 and 18)
are preserved, remaining the previous steps 1. and 2. unaltered. The mod-
ification lies in the computation of the diffusion terms (in eq. 15 and 16),
which are no longer only Laplacian operators when $\alpha(\mathbf{x}, t)$ and $\alpha/G = \alpha_s(\mathbf{x}, t)$
are variable in space. In this case, the diffusion terms can be expanded as
$$ \nabla \cdot (\alpha \nabla \Theta) = \alpha \nabla^2 \Theta + \nabla \alpha \cdot \nabla \Theta \text{ and } \nabla \cdot \left(\frac{\alpha}{G} \nabla \Theta_s\right) = \frac{\alpha}{G} \nabla^2 \Theta_s + \nabla \frac{\alpha}{G} \cdot \nabla \Theta_s, $$
requiring the extra estimation of the product of diffusivity and temperature gra-
dients by numerical differentiation (here, 6$^{th}$-order first-derivative schemes after
solution reconstruction in the immersed region). For generality, the fluid/solid
diffusivities are expressed in terms of spatial and time dependency, but most of
the time, this dependency comes from a known functional relation with the tem-
perature as expressed by $\alpha(\Theta)$ and $\alpha_s(\Theta_s)$. In section 4.1, we have addressed
this more demanding situation of variable conductivity through a formal test-
case based on an analytical solution for a purely conductive problem with a set
of numerical *BC* equivalent to those for *CHT*.

Typical temperature profiles obtained for *CHT* in a laminar channel are
presented in Figure 6. Their reconstructed part is shown to illustrate the reg-
ularity condition separately satisfied by $\Theta$ and $\Theta_s$ at the fluid/solid interface,
for the benefit of the spatial differentiation accuracy. In this example, where
$G_2 = 2$, the wall-normal derivative discontinuity of the physically meaningful
temperature (i.e., Θ in the fluid and Θ_s in the solid) can be clearly observed
through the check of condition (23). These laminar results are discussed more
quantitatively in the next section.

<figure><img src="image_8.png" /><figcaption>Figure 6: Fluid streamwise velocity and temperature fields reconstruction. This example represents an instantaneous field, for $G_2 = 2$, before complete convergence to the laminar solution.</figcaption></figure>

4. VALIDATION RESULTS

To assess the present dual *IBM* to model fluid/solid geometry in the context
of *DNS* including *CHT* phenomena, a validation decomposed into three stages
is adopted, considering pure conduction through a composite wall, *IT*/*IF*/*CHT*
in channel and *IT* in pipe flow configuration. First, in section 4.1, the pure con-
duction problem through a composite wall is studied to evaluate the capability
and accuracy of the dual *IBM* to deal with variable conductivity. Secondly, the
laminar case is addressed in sections 4.2 and 4.3, in order to assess the actual ac-
curacy of the schemes while checking if they can ensure numerical stability. For
this preliminary step, the ability to impose *IT*, *IF* conditions (which are both
analytically equivalent in the laminar case) and *CHT* will be assessed with a rig-
orous evaluation of the error introduced by the various approximations. Then,
the turbulent case is investigated in sections 4.4 and 4.5 imposing *IT*, *IF* and
*CHT* conditions. This last case is more demanding in the sense that it requires
the local (in time and space) imposition of various thermal boundary conditions at the fluid/solid interface, as explained in section 3.5.

### 4.1. Pure conduction through a composite wall

Here, the temperature inside a composite wall is considered. The wall is assumed to be infinite in the xz-plane with a temperature field fully established both temporally and spatially making one-dimensional the problem with formally $T(y)$. Between $-L_y/2 \le y \le -H/2$ and $H/2 \le y \le L_y/2$, the conductivity $\lambda_s$ is constant. In the central region the conductivity is variable following the law

$$ \lambda_f = a_f + b_f T \quad \forall y \in [H/2, -H/2] \qquad (24) $$

where $a_f, c_f$ are constants. In this framework, the heat equation reads

$$ \frac{\partial}{\partial y} \left( \lambda_f \frac{\partial T}{\partial y} \right) = 0 \qquad (25) $$

By analytical integration of this non-linear equation, the exact temperature solution can be obtained with

$$ T(y) = \frac{-a_f + \sqrt{a_f^2 + 2b_f(c_1y + c_2)}}{b_f} \quad \forall y \in [-H/2, H/2] \qquad (26) $$

where the constants $c_1$ and $c_2$ are related to the Dirichlet boundary conditions $T(-H/2) = T_1$ and $T(H/2) = T_2$ as expressed by

$$ c_1 = \frac{a_f}{H} (T_2 - T_1) + \frac{b_f}{2H} (T_2^2 - T_1^2) \qquad (27) $$

$$ c_2 = \frac{a_f}{2} (T_1 + T_2) + \frac{b_f}{4} (T_1^2 + T_2^2) \qquad (28) $$
In the lower and upper walls, the Dirichlet boundary conditions

$$T(-L_y/2) = T_0 \qquad (29)$$

$$T(L_y/2) = T_3 \qquad (30)$$

leads to the simple expressions

$$T(y) = \frac{1}{L_y - H} [2(T_1 - T_0)y + T_1 L_y - T_0 H] \quad \forall y \in [-L_y/2, -H/2] \quad (31)$$

$$T(y) = \frac{1}{L_y - H} [2(T_3 - T_2)y + T_2 L_y - T_3 H] \quad \forall y \in [H/2, L_y/2] \quad (32)$$

where $T_0$ and $T_3$ can be deduced from the continuity of the normal heat flux

$$\vec{\phi}_w = -c_1 \vec{y} \qquad (33)$$

leading to

$$T_0 = T_1 - c_1 \frac{L_y - H}{2\lambda_s} \qquad (34)$$

$$T_3 = T_2 + c_1 \frac{L_y - H}{2\lambda_s} \qquad (35)$$

The parameters were set to be highly demanding for the dual numerical approach in terms of temperature gradient discontinuity and thermal interaction at the interface. With $a_f = 0.5 W m^{-1} K^{-1}$, $b_f = 0.5 W m^{-1} K^{-1}/^{\circ}C$, $(\rho c_p)_f = (\rho c_p)_s = 10^{-7} J m^{-3} K^{-1}$, $T_1 = 10^{\circ}C$ and $T_2 = 20^{\circ}C$, at the interface, the solid-to-fluid conductivity ratio is $G_2 = \lambda_s/\lambda_f = 20$; the fluid-to-solid diffusivity ratio is $G = \alpha_f/\alpha_s = 0.05$ and a thermal activity ratio $K = \lambda_f/\lambda_s \sqrt{\alpha_s/\alpha_f}$, which regulates the intensity of the thermal coupling, is $K = 1/\sqrt{20}$. Figure 7 illustrates the temperature solution, revealing the dual IBM capability required to deal with variable conductivity and sharp temper-
ature gradients at the interface.

<figure><img src="image_9.png" /><figcaption>Figure 7: Temperature solution in the composite wall.</figcaption></figure>

The validation tests have been performed in a computational domain of $L_y = 2H$, using five spatial resolutions $n_y = (103, 203, 403, 803)$. The numerical solution is initialised with the exact solution to reduce the convergence time. The convergence criterion is based on the difference between two distant time steps ensuring that the difference remains at least seven orders of magnitude lower than the maximum error. The numerical convergence is evaluated for 5$^{th}$- and 7$^{th}$-order polynomial reconstruction in the immersed region. As observed in figure 8, a 5$^{th}$-order polynomial reconstruction leads to a 2$^{nd}$-order convergence of the norm $L_\infty$ with a slightly higher convergence of $L_2$-norm. The convergence order is found to increase with the polynomial order through a 3$^{rd}$-order accuracy for the 7$^{th}$-order reconstruction. In the present one-dimensional and purely conductive case based on 6$^{th}$-order finite-difference schemes, this behaviour is expected because the main error contribution comes from the solution reconstruction. However, this possibility to increase the accuracy order in this way can be seen as useless for a fluid solver subjected to the 2$^{nd}$-order limitation of the pressure treatment in presence of IBM.
<figure><img src="image_10.png" /><figcaption>Figure 8: Numerical convergence of the temperature inside a composite wall of variable con ductivity. Comparison between 5$^{th}$- and 7$^{th}$-order polynomial reconstruction in the immersed regions. Left: L$_{2}$-norm. Right: L$_{∞}$-norm.</figcaption></figure>

4.2. Laminar Poiseuille channel flow

In laminar Poiseuille channel flow, with − H/ 2 ≤ y ≤ H/ 2, the velocity analytical solution is

$$
u_x(y) = \frac{3}{2} \left(1 - 4 \frac{y^2}{H^2}\right) U_b, \quad u_y = u_z = 0, \qquad (36)
$$

and, for IT, the temperature solution is

$$
\Theta(y) = \frac{Nu}{16} \left( 5 - 24 \frac{y^2}{H^2} + 16 \frac{y^4}{H^4} \right), \quad (37)
$$

while, for IF/CHT conditions

$$
\Theta(y) = \frac{1}{16} \left( 5 - 24 \frac{y^2}{H^2} + 16 \frac{y^4}{H^4} \right) - \frac{1}{Nu}, \quad (38)
$$

with Nu = 140 / 34 ≈ 4 . 118 at any Reynolds and Prandtl numbers. When CHT is considered, the temperature profile inside the solid is linear, given by

$$
\Theta_s(y) = \frac{1}{G_2} \left( \frac{1}{2} - \frac{|y|}{H} \right) - \frac{1}{Nu} \quad \text{for } H/2 \le |y| \le L_y/2. \quad (39)
$$

These solutions are used as initial conditions and as an exact reference to estimate the numerical error corresponding to a deviation from them. This
deviation can be observed at every mesh node (with a maximum absolute value equal to the $L_\infty$-norm) or estimated globally through the computation of the $L_2$-norm. Even if any initial condition can lead to the same established flow, the use of (36, 37) as initial conditions is a mean to minimise the transient stage during which the numerical solution has to adapt. In practice, for every test, it has been checked that the solution was fully established and converged in time, enabling to consider the numerical errors as only due to spatial discretization. The convergence criterion is based on the difference between two distant time steps ensuring that the difference remains at least seven orders of magnitude lower than the maximum error.

Despite the 1D nature of the laminar solution, the validation tests have been performed using the 3D version of the code by choosing arbitrary the streamwise and spanwise lengths $L_x = L_z = H/2$ of the computational domain and the associated numbers of mesh nodes $n_x = n_z = 10$. Naturally, it has been observed that in practice, the results are not dependent on these parameters. To exhibit the convergence order, every case among IT, IF and CHT has been computed using four spatial resolutions in the wall-normal direction as indicated in Table 1.

<table><caption>Table 1: Computational parameters for the laminar test cases under different Thermal Boundary Conditions on the wall (TBC). For all CHT simulation, the diffusivity ratios is set to $G=1$.</caption><thead><tr><th>TBC</th><th>n<sub>y</sub></th><th>L<sub>y</sub>/H</th><th>&Delta;y/H &asymp;</th><th>G<sub>2</sub></th></tr></thead><tbody><tr><td>IT/IF</td><td>( 66, 130, 258, 514)</td><td>1.28</td><td>(20, 10, 5, 2.5) &times; 10<sup>-3</sup></td><td>&mdash;</td></tr><tr><td>CHT</td><td>(103, 203, 403, 803)</td><td>2</td><td>(20, 10, 5, 2.5) &times; 10<sup>-3</sup></td><td>(1, 2, 4)</td></tr></tbody></table>

The resulting convergence curves are presented in Figure 9, for two error norms, namely the maximum ($L_\infty$-norm) and the mean square ($L_2$-norm) error of the numerical solution relative to the analytical one. For the velocity, as explained in GAUTIER et al. (2013, 2014), because of the lack of any specific
treatment of the pressure at the fluid/solid interface, the expected convergence
order is only two for the velocity. This is a common feature of most of IBM
when applied on a complex geometry in the context of incompressible flow. For
the present channel geometry, a higher convergence order could be obtained for
the temperature field alone using the exact velocity field. However, because
in a practical calculation, the velocity error also governs the temperature er-
ror, the choice in this section is to use a reconstruction scheme that is only
third-order accurate for the *IT* case and second-order accurate for the *IF* and
*CHT* cases (see sections 3.3 and 3.4). The convergence curves presented in
Figure 9 are consistent with these expectations. Both error norms show that
the velocity solution is slightly better than second-order accurate. This is an
artefact of the minor role of the second-order accurate pressure correction for
the present channel geometry. In this particular situation, the convergence re-
flects the third-order accuracy of the reconstruction scheme. For the *IT* case,
the same behaviour is logically recovered for the temperature. On the con-
trary, for the *IF* case, a remarkable second-order convergence can be observed
for the temperature. This behaviour is consistent with the formal accuracy of
the reconstruction scheme, which is only second-order accurate for this type of
*BC* because of the use of scheme (21). This convergence test shows that the
present *IBM* for an *IF* condition is numerically stable while validating its actual
order of accuracy. The same can be concluded for the *CHT* case for which a
second-order accuracy is also recovered.

This first success in the *CHT* case deserves further comments. First, it worth
noting that despite its fully explicit nature, the present implementation of the
thermal coupling is numerically stable while enabling to represent smoothly
situations in which the normal derivative of the temperature is discontinuous
at the fluid/solid interface. This discontinuity must occur when the fluid and
<figure><img src="image_11.png" /><figcaption>Figure 9: Convergence of $\Theta$ and $u_x$. From top to bottom: IT case, IF case, CHT case. Left: $L_2$-norm. Right: $L_\infty$-norm.</figcaption></figure>
solid conductivities are different, a pervasive situation in industrial applications. For the present laminar Poiseuille flow, the discontinuous wall-normal derivative obtained in present CHT simulations is illustrated in Figure 10, for $G = 1$ and conductivity ratios $G_2 = 1, 2, 4$. The temperature continuity condition ($\Theta = \Theta_s$) and the heat-flux continuity condition ($\partial_n \Theta_s = \frac{1}{G_2} \partial_n \Theta$) can be clearly observed at the interface while ensuring the regularity of each $\Theta$ and $\Theta_s$ profile which are only meaningful in their own physical domain, namely the fluid and solid domain, respectively.

<figure><img src="image_12.png" /><figcaption>Figure 10: Temperature profiles under conjugate heat transfer (CHT), for diffusivity ratio $G = 1$ and different conductivity ratios ($G_2 = 1, 2, 4$). Dashed lines indicate non-physical reconstructed field, as a part of the dual IBM present technique.</figcaption></figure>

### 4.3. Laminar Poiseuille pipe flow

In laminar Poiseuille pipe flow, with $r \le D/2$, the velocity analytical solution is

$$u_x(r) = 2 \left( 1 - 4 \frac{r^2}{D^2} \right) U_b, \quad u_\theta = u_r = 0. \qquad (40)$$
where $r$ is the radial coordinate and $u_\theta/u_r$ the azimuthal/radial velocity. The dimensionless temperature solution, for IF/IT, is

$$ \Theta(r) = 2Nu \left[ \frac{3}{16} + \frac{r^4}{D^4} - \frac{r^2}{D^2} \right] \quad \text{with} \quad Nu = -D \frac{\partial \Theta}{\partial r} \Big|_{r=D/2} = 4.36 \quad (41) $$

The validation tests have been performed choosing the streamwise length $L_x = D/2$ and the associated numbers of mesh nodes $n_x = 10$. The solutions have been computed using four spatial resolutions in the cross-section plan $n_y = n_z = (64, 128, 256, 512)$ in a domain size $L_y/D = L_z/D = 1.28$.

Despite the disconnection of the mesh organisation from the wall geometry for the present pipe flow configuration, the convergence behaviour is similar to what has been observed for the channel flow. Figure 11 shows the convergence curves for the two types of norm considered. The $L_2$-norm shows higher than second-order of accuracy for both velocity and temperature. When the maximum error $L_\infty$-norm is examined, a second-order accuracy is recovered as expected. It can be easily checked that this maximum error is located near the fluid/solid interface. More generally, the main contribution to the $L_2$-norm comes from at the nearest nodes from the wall consistently with the observations reported in GAUTIER et al. (2013) for a cylindrical but quite different flow configuration (stationary flow over a cylinder at low Reynolds number).

It should be pointed out that, conversely to channel flow, in pipe flow, the distance of the solid boundary to the nearest fluid mesh node varies azimuthally along the pipe wall. It means that the magnitude of the error changes along the fluid/solid interface. Despite this important change by comparison to the plane channel geometry, it is worth noting that the order of convergence remains unaffected, confirming that the current 1D reconstruction strategy is consistent irrespective of the wall orientation in relation to the computational mesh. This is a very important feature in the context of IBM that is checked here in the
context of heat transfer.

<figure><img src="image_13.png" /><figcaption>Figure 11: Convergence of $\Theta$ and $u_x$, for IT condition. Left: $L_2$-norm. Right: $L_\infty$-norm.</figcaption></figure>

To close these assessments of the convergence order, and as a transition to the next section, some comments must me made. First, it has to be recognised that various approximation orders are used in our IBM from the 5$^{th}$-order of the polynomial reconstruction up the 2$^{nd}$-order of the pressure treatment. In practice, the latter is the limiting factor so that the use of any higher-order approximation may be seen as useless. Even if this view is correct for the formal convergence order, our experience is that an increase of the reconstruction order and its related approximations can increase the smoothness of the solution with an improvement of the accuracy at marginal resolution when turbulent flows are addressed. That is why we have chosen the present combination of 5$^{th}$- and 7$^{th}$-order approximations. Another legitimate issue is to question the interest of the use of a 6$^{th}$-order solver with an essentially 2$^{nd}$-order IBM. For laminar flow configurations, there are neither advantages nor drawbacks in the use of a high-order fluid solver. But for the computation of turbulence, a clear increase of accuracy can be observed for the prediction of basic turbulent statistics. For this second perspective, our experience is that the use of high-order schemes leads to a better ratio in terms of accuracy/cost. Some previous investigations about these points, namely 2$^{nd}$-order vs. high-order schemes with and without IBM, can be found in Parnaureau et al. (2004); LAIZET and LAMBALLAIS (2009);
GAUTIER et al. (2013); James et al. (2013).

### 4.4. Turbulent channel flow

In this section, three DNS of turbulent channel flow are presented depending on the type of thermal boundary conditions considered. The main physical and numerical parameters are summarised in Tables 2 and 3. These parameters were chosen to match those of previous studies to enable comparisons for validation purposes. All the presented turbulent statistics correspond to fully established flows.

The results are validated using the DNS database of FLAGEUL et al. (2015, 2017) as reference. In these previous studies, the same code Incompact3d was used in the fluid domain while coupling it with a conduction solver in the solid domain for the CHT cases. The spatial differentiation of this companion solver is based on the same finite difference schemes as in Incompact3d in the streamwise and spanwise directions but uses a Chebyshev method in the wall-normal direction. As a consequence, the methods used by FLAGEUL et al. (2015, 2017) are very close to those used here, the main difference being the present treatment of wall boundary conditions by IBM.

Another difference is the absence of mesh refinement near the fluid/solid interface in this study. This choice was made for two reasons. First, technically, the code Incompact3d enables an easy mesh refinement at the boundaries of the computational domain $y = \pm L_y/2$ or at the center $y = 0$ but not in an intermediate region like $y = \pm H/2$. This constraint is connected to the choice to solve the Poisson equation in the Fourier space (LAIZET and LAMBALLAIS, 2009). A specific numerical development could make more flexible the mesh node distribution, but it would be against the philosophy of the present IBM in which the mesh organisation is not directly connected to the body geometry, enabling the use of a simple Cartesian mesh to handle complex geometry. This
is the second reason why a uniform mesh is used even in the wall-normal direction. This makes the following validation more representative of a typical IBM application. By comparison to the reference resolution, the present wall-normal resolution is finer far from the walls and coarser close to them, as reported in Table 3. In particular, the cell size nearby the wall $\Delta y^+ = 1.5$ does not allow an accurate description of the near-wall dynamics in the viscous sublayer. This marginal near-wall resolution leads to a slight underestimation of the friction velocity of about 0.7%. This level of accuracy was considered acceptable for the present validation.

<table><caption>Table 2: Main simulation parameters, for turbulent cases. Thermal Boundary Conditions (TBC): imposed temperature (IT), isoflux of heat (IF), conjugate heat transfer (CHT).</caption><thead><tr><th>TBC</th><th>Re</th><th>Pr</th><th>n<sub>x</sub> × n<sub>y</sub> × n<sub>z</sub></th><th>(L<sub>x</sub> × L<sub>y</sub> × L<sub>z</sub>) / H</th></tr></thead><tbody><tr><td>IT</td><td>4560</td><td>0.71</td><td>256 × 256 × 256</td><td>12.8 × 1.28 × 4.26</td></tr><tr><td>IF</td><td>4560</td><td>0.71</td><td>256 × 256 × 256</td><td>12.8 × 1.28 × 4.26</td></tr><tr><td>CHT</td><td>4560</td><td>0.71</td><td>256 × 401 × 256</td><td>12.8 × 2 × 4.26</td></tr></tbody></table>

<table><caption>Table 3: Simulation parameters and the reference ones.</caption><thead><tr><th>Work</th><th>Current</th><th>FLAGEUL et al. (2015)</th></tr></thead><tbody><tr><td>⟨Re<sub>τ</sub>⟩</td><td>148</td><td>149</td></tr><tr><td>Re'<sub>τ rms</sub></td><td>0.77</td><td>—</td></tr><tr><td>Pr</td><td>0.71</td><td>0.71</td></tr><tr><td>⟨N<sub>uD</sub>⟩</td><td>15.1</td><td>—</td></tr><tr><td>N'<sub>uD rms</sub></td><td>0.27</td><td>—</td></tr><tr><td>[L<sup>+</sup><sub>x</sub>, L<sup>+</sup><sub>y</sub>, L<sup>+</sup><sub>z</sub>]</td><td>[3789, 296, 1259]</td><td>[3814, 298, 1269]</td></tr><tr><td>Δy<sup>+</sup></td><td>1.5</td><td>[0.49; 4.8]</td></tr><tr><td>[Δx<sup>+</sup>, Δz<sup>+</sup>]</td><td>[14.8, 4.9]</td><td>[14.8, 4.9]</td></tr><tr><td>Δt<sup>+</sup></td><td>0.03</td><td>0.02</td></tr><tr><td>Duration</td><td>2900</td><td>29000</td></tr></tbody></table>

Mean velocity and Reynolds stress profiles are presented in Figure 12. A good agreement with the reference data can be observed, the slight discrepancy corresponding to the underestimation of the friction velocity as a consequence of the use of a coarse resolution in the near-wall region.

Interestingly, temperature statistics are less affected by this choice of near-
<figure><img src="image_14.png" /><figcaption>Figure 12: Velocity statistics in wall-units (continuous line). Comparison between present (continuous line) and FLAGEUL et al. (2015) data (dashed line). Left: Mean streamwise velocity profile $\langle u_x^+ \rangle$. Right: Reynolds stress profiles.</figcaption></figure>

wall resolution with a remarkable collapse between the present and reference data of FLAGEUL et al. (2015) for the *IT* case (see the profiles of mean temperature, root mean square fluctuating temperature and turbulent heat flux presented in Figure 13). A similar agreement can be noticed in Figure 14 for the *IF* case except for $\Theta_{\text{rms}}^{+}$ where a small deviation from the reference data can be observed in the viscous sublayer. However, this effect is found to be limited with a relative error on the statistics less than 6%. As for the velocity statistics, this slight discrepancy can be related to the lack of any mesh refinement near the wall. By comparison to the *IT* case, it is coherent to observe a slightly stronger sensitivity to the use of a coarse near-wall resolution, since the imposition of the flux requires an extra approximation. However, in view of the mesh resolution, the present results can be considered as fully satisfactory for both the *IT* and *IF* cases as the first crucial successful step for dealing with CHT, which can be viewed, technically, as a combination of these two types of boundary conditions.
<figure><img src="image_15.png" /><figcaption>Figure 13: Temperature statistics in wall-units for the *IT* case. Comparison between present (continuous line) and FLAGEUL et al. (2015) data (dashed line). Top row: Mean temperature profile. Center row: Temperature root mean square profile. Bottom row: Turbulent heat flux profile. Left/right column: wall distance in lin/log scale.</figcaption></figure>

Figure 15 shows the temperature statistics, for the *CHT* case. For this validation case, unitary ratios of thermal conductivities and diffusivities ($G = G_2 = 1$) are considered, identical to the ratios defined by FLAGEUL et al. (2015). The use of equal diffusivities is a condition to remain distant from both *IT* (i.e. $G_2 \to \infty$) and *IF* (i.e. $G_2 \to 0$) conditions, allowing the *CHT* processes to be fully interactive. In this demanding context, the agreement between present
and FLAGEUL et al. (2015) data (Figure 15) can be evaluated as satisfactory. The discrepancy of $\Theta_{rms}^{+}$ remains less than 6% as for the *IF* case whereas the turbulent heat flux is better predicted by comparison to the data of FLAGEUL et al. (2015). In view of the coarse near-wall resolution, this error, by comparison to a case free from any *IBM* approximation, as in FLAGEUL et al. (2015), can be considered as acceptable. Furthermore, the strong damping of the temperature fluctuations in the solid region with the fluid/solid interface distance can be observed in Figure 15. Exponential damping is well recovered here, as suggested by a simplified spectral analysis (FLAGEUL et al., 2015).

Figure 16 presents an instantaneous visualisation of the temperature field in both the fluid and solid domains. This visualisation well illustrates the thermal fluid/solid interaction. In particular, the turbulent structures are clearly visible on the temperature distribution in the core of the flow as well as in the near-wall region. Note that, in the solid, these temperature fluctuations correspond to deviation from a linear mean temperature profile in *y* as correctly predicted in present results. This global agreement for the *CHT* turbulent case shows that the present dual *IBM* remains numerically stable for an unsteady flow while preserving accuracy for the prediction of basic turbulent statistics.
<figure><img src="image_16.png" /><figcaption>Figure 14: Temperature statistics in wall-units for the *IF* case. Comparison between present (continuous line) and FLAGEUL et al. (2015) data (dashed line). Top row: Mean temperature profile. Centre row: Temperature root mean square profile. Bottom row: Turbulent heat flux profile. Left/right column: wall distance in lin/log scale.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Figure 15: Temperature statistics in wall-units for the CHT case. Comparison between present (continuous line) and FLAGEUL et al. (2015) data (dashed line). Top row: Mean temperature profile. Centre row: Temperature root mean square profile. Bottom row: Turbulent heat flux profile. Left/right column: wall distance in lin/log scale.</figcaption></figure>

4.5. *Turbulent pipe flow*

In this section, velocity/temperature statistics in turbulent pipe flow for
IT condition are validated, with their reference DNS counterparts, at various
Reynolds and Prandtl numbers using the same spatial resolution (n_x×n_y×n_z) =
<figure><img src="image_18.png" /><figcaption>Figure 16: From top to bottom, instantaneous temperature field $\Theta$, temperature fluctuation $\Theta'$ and absolute value of the fluctuations $|\Theta'|$. This last visualisation is in logarithmic scale. The visualisation corresponds to CHT simulation (see Table 2), with unitary ratios of thermal diffusivities and conductivities ($G = G_2 = 1$).</figcaption></figure>

(640 × 256 × 256) for the same computational domain ($L_x \times L_y \times L_z$)/D = (12.5 × 1.28 × 1.28) and with the same time step $\Delta t = 8 \times 10^{-4} D/U_b$. A 2D view of the cross-flow uniform mesh is shown in Figure 3. This computational configuration where a pipe flow geometry is addressed with a Cartesian mesh
thanks to the same IBM and fluid solver is documented in DAIRAY et al. (2017).

The low Reynolds number case $Re = 5500$ is addressed by DNS whereas for $Re = 17000, 19000$, the implicit LES strategy of LAMBALLAIS et al. (2011); DAIRAY et al. (2017) is adopted. The latter is based on the introduction of a controlled numerical dissipation corresponding to an implicit spectral vanishing viscosity which preserves the 6^th^-order accuracy of the differentiation. The corresponding nominal values of friction Reynolds number $Re_τ = u_τD/(2ν)$ are $Re_τ = 186, 500, 550$ respectively.

Tables 4 and 5 show numerical and physical parameters for current and reference simulations. The experimental work of ISSHIKI et al. (1993) is included to assess the Nusselt number estimation. It can be observed that in the core of the pipe, the mesh refinement associated to the current use of a regular Cartesian mesh is higher than EL KHOURY et al. (2013) and SAHA et al. (2014) in the radial direction but coarser than any reference studies in the wall-normal direction close to the wall. Despite this difference, an excellent agreement with (REDJEM-SAAD et al., 2007; SAHA et al., 2014) is found for the prediction of the friction velocity at low Reynolds number. The conclusion is the same for the present LES at $Re = 17000$ and $Re = 19000$, consistently with the trends reported in DAIRAY et al. (2017).

Concerning the Nusselt number, an over-prediction of about 8% is obtained at low Reynolds number by comparison to the experimental value documented in ISSHIKI et al. (1993). This discrepancy can be accepted considering that the experiment was carried out at slightly lower Reynolds and Prandtl numbers. At $Re = 17000$, Nu is under-predicted by about 11% from the value reported by SATAKE et al. (2000) while being within a deviation of 2% from the empirical correlation provided by GNIELINSKI (1976) as illustrated in Figure 17. This level of agreement with respect to experimental results of reference can be considered
as a satisfactory validation of the present IBM for heat transfer prediction.

<table><caption>Table 4: DNS references. Friction Reynolds number, Prandtl number and Nusselt number.</caption><thead><tr><th>Work</th><th>Re</th><th>⟨Re<sub>τ</sub>⟩</th><th>Re'<sub>τ rms</sub></th><th>Pr</th><th>⟨N<sub>u</sub>⟩</th><th>N'<sub>u rms</sub></th></tr></thead><tbody><tr><td>Current DNS</td><td>5500</td><td>186.46</td><td>1.223</td><td>1.0</td><td>21.73</td><td>0.631</td></tr><tr><td>REDJEM-SAAD et al. (2007)</td><td>5500</td><td>186</td><td>—</td><td>1.0</td><td>—</td><td>—</td></tr><tr><td>SAHA et al. (2014)</td><td>5500</td><td>186</td><td>—</td><td>1.0</td><td>—</td><td>—</td></tr><tr><td>ISSHIKI et al. (1993)</td><td>5300</td><td>180</td><td>—</td><td>0.71</td><td>20.06</td><td>—</td></tr><tr><td>Current LES1</td><td>17000</td><td>495.79</td><td>1.675</td><td>0.71</td><td>46.05</td><td>0.637</td></tr><tr><td>Current LES2</td><td>19000</td><td>550.51</td><td>1.694</td><td>0.71</td><td>50.83</td><td>0.577</td></tr><tr><td>SATAKE et al. (2000)</td><td>17000</td><td>500</td><td>—</td><td>0.71</td><td>≈ 52</td><td>—</td></tr><tr><td>EL KHOURY et al. (2013)*</td><td>19000</td><td>550</td><td>—</td><td>—</td><td>—</td><td>—</td></tr></tbody></table>

<table><caption>Table 5: DNS references and the present simulations. Spatial resolution in wall-units and method. FDⁿ: Finite differences method of nᵗʰ-order; FV: finite volumes method; SE: spectral elements method; Exp.: Experimental.</caption><thead><tr><th>Work</th><th>Δz<sup>+</sup></th><th>“Radial”</th><th>“Azimuthal”</th><th>Method</th></tr></thead><tbody><tr><td>Current DNS</td><td>7.5</td><td>Δx<sup>+</sup> = Δy<sup>+</sup> = 1.9</td><td>Δx<sup>+</sup> = Δy<sup>+</sup> = 1.9</td><td>FD<sup>6</sup></td></tr><tr><td>REDJEM-SAAD et al.</td><td>10.5</td><td>Δr<sup>+</sup> ∈ [0.29; 1.04]</td><td>Δ(Rθ)<sup>+</sup> = 8.84</td><td>FD<sup>2</sup></td></tr><tr><td>SAHA et al.</td><td>9.3</td><td>Δr<sup>+</sup> ∈ [0.061; 3.85]</td><td>Δ(Rθ)<sup>+</sup> = 9.13</td><td>SE</td></tr><tr><td>ISSHIKI et al.</td><td>—</td><td>—</td><td>—</td><td>Exp.</td></tr><tr><td>Current LES1</td><td>19.4</td><td>Δx<sup>+</sup> = Δy<sup>+</sup> = 5.0</td><td>Δx<sup>+</sup> = Δy<sup>+</sup> = 5.0</td><td>FD<sup>6</sup></td></tr><tr><td>Current LES2</td><td>21.5</td><td>Δx<sup>+</sup> = Δy<sup>+</sup> = 5.5</td><td>Δx<sup>+</sup> = Δy<sup>+</sup> = 5.5</td><td>FD<sup>6</sup></td></tr><tr><td>SATAKE et al.</td><td>14.6</td><td>Δr<sup>+</sup> ∈ [0.1; 2.6]</td><td>Δ(Rθ)<sup>+</sup> = 8.18</td><td>FV</td></tr><tr><td>EL KHOURY et al.</td><td>[3.06; 9.99]</td><td>Δr<sup>+</sup> ∈ [0.15; 4.49]</td><td>Δ(Rθ)<sup>+</sup> = [1.45; 4.75]</td><td>SE</td></tr></tbody></table>

<figure><img src="image_19.png" /><figcaption>Figure 17: Nusselt number for turbulent pipe flow.</figcaption></figure>

To project the velocity/temperature information from the Cartesian mesh

<footnote>* In EL KHOURY et al. (2013) work it is only solved the mometum equation. It was included to evaluate the velocity statistics because of the spectral code employed.</footnote>
toward the radial coordinate $r$, $2^{nd}$-order procedure proposed by DAIRAY et al. (2017) is applied. This procedure yields azimuthally-averaged profiles as a function of $r$ with a radial resolution $\Delta r$ that can be chosen freely with respect to the Cartesian mesh resolution where $\Delta y = \Delta z$. In particular, profiles can be generated at a finer resolution than the computational mesh, namely $\Delta r < \Delta y$. Here, for the velocity (Figure 18) and temperature (Figure 19 and 20) profiles, a radial resolution such as $\Delta r = \Delta y/8$ is used, enabling to capture the near-wall behaviour at a scale significantly smaller than the cell size of the computational mesh.

At $Re = 5500$, an excellent agreement of basic velocity statistics with their reference counterpart can be observed in figure 18. This is the confirmation that the present numerical methodology, based on a regular Cartesian mesh while using an *IBM* with a coarse near-wall resolution, can provide the *DNS* accuracy. This agreement is slightly less favourable at the higher Reynolds numbers $Re = 17000, 19000$ for which the near-wall resolution is even coarser in wall units. The main discrepancy is related to the streamwise Reynolds stress component $\langle u_x u_x \rangle$ which is underestimated by about 7.5%. This behaviour is expected in the LES framework where only a part of the kinetic energy can be captured with also an unavoidable loss of accuracy due to the errors introduced by the subgrid-scale modelling.
<figure><img src="image_20.png" /><figcaption>Figure 18: Current velocity statistics (continuous line) compared to DNS references REDJEM-SAAD et al. (2007), SATAKE et al. (2000) and EL KHOURY et al. (2013) (dashed lines). Left: Mean velocity profiles $\langle u_x \rangle^+$. Right: Reynolds stresses profiles. $\langle u'_i u'_j \rangle^+$ Top: $Re = 5500$ ($Re_τ = 186$); center: $Re = 17000$ ($Re_τ = 500$); bottom: $Re = 19000$ ($Re_τ = 550$). Azimuthal average with $\Delta r = \Delta y/8$.</figcaption></figure>

Similarly, the temperature statistics presented in Figure 19 exhibits a remarkable agreement with the reference DNS of REDJEM-SAAD et al. (2007); SAHA et al. (2014) at $Re = 5500$. The peaks of the temperature fluctuations $\Theta_{rms}^+$ and the turbulent heat flux $\langle u_x' \Theta' \rangle^+$ are found to be slightly higher than the reference ones with an overestimation of about 2% and 3% respectively. As expected, the implicit LES performed at higher Reynolds numbers $Re = 17000, 19000$ lead to a less favourable agreement with reference data. The
main discrepancy with the *DNS* data of SATAKE et al. (2000) is located in the core of the pipe as shown in figure 20. It is worth noting that the corresponding overestimation of temperature fluctuations can be related to the typical protuberance of turbulent statistics in the pipe core region ($y^+ > 100$) as a signature of large-scale structures observed in wall turbulence at high Reynolds number.

<figure><img src="image_21.png" /><figcaption>Figure 19: Temperature statistics, at $Re_τ = 186$ ($Re = 5500$) and $Pr = 1$ compared to REDJEM-SAAD et al. (2007) and SAHA et al. (2014). Top row: Mean temperature profile $\langle\Theta\rangle^+$. Center row: Temperature root mean square profile $\Theta_{rms}^+$. Bottom row: Turbulent heat flux profile $\langle u_x' \Theta' \rangle^+$. Azimuthal average with $\Delta r = \Delta y/8$.</figcaption></figure>

More qualitatively, instantaneous velocity and temperature fields can be observed in a cross-section (see figure 21) or in the whole computational domain (see figure 22). The similar pattern obtained for the temperature and longitudinal velocity is as expected through the dominance of convective effects on the temperature fluctuations. The most energetic near-wall scales, namely the high- and low-speed streaks, can be clearly observed with an averaged spacing of about 100 wall-units for the low and high Reynolds number considered. This value, widely recovered in various wall-turbulence configurations, suggests that the dynamics of the main turbulence structures is correctly captured in the
near-wall region where the influence of the numerical errors related to the *IBM* can be assumed to be the strongest. Another remarkable feature is the lack of spurious oscillations, even for the implicit LES at high Reynolds numbers. In particular, the quality of the imposition of the thermal boundary conditions by the present *IBM* does not seem to be affected by the use of a coarse mesh resolution or to the application of an artificial dissipation as a procedure of implicit subgrid-scale modelling. This favourable feature open the way to high-fidelity simulation of turbulent heat-transfer at low computational cost as an attractive compromise in the framework of *DNS/LES*.

<figure><img src="image_22.png" /><figcaption>Figure 20: Temperature statistics (at $Re_τ$ = 500, $Re_τ$ = 550 and $Pr$ = 0.71) compared to SATAKE et al. (2000) ($Re_τ$ = 500). Mean temperature profile (Top) and Temperature root mean square profile (Bottom). Azimuthal average with $\Delta r = \Delta y/8$.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>Figure 21: Instantaneous temperature $\Theta$ (first row) and streamwise velocity $u_x$ (second row) fields at cross section in $x = L_x/2$. Left: $Re_\tau = 186$ ($Re = 5500$) and $Pr = 1$; Right: $Re_\tau = 550$ ($Re = 19000$) and $Pr = 0.71$.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Figure 22: Iso-surfaces of temperature. Top: $Re_\tau = 186$ and $Pr = 1$; Bottom: $Re_\tau = 550$ and $Pr = 0.71$.</figcaption></figure>
5. CONCLUSIONS

This paper is an investigation of the potential of immersed boundary method
for the treatment of conjugate heat transfer. For that purpose, a stepwise ap-
proach is adopted by considering successively: ideal locally imposed temper-
ature, ideal locally imposed heat flux, and real conjugate heat transfer as a
fluid/solid thermal coupling problem. As a generic solver for direct numerical
simulation, the highly accurate code *Incompact3d* is used with a customised
immersed boundary method enabling to treat accurately and efficiently any
complex geometry thanks to one-dimensional reconstruction schemes inside the
immersed region. For the temperature field, the solid is the immersed region
for the fluid and conversely. This dual definition of the immersed zone leads to
the incorporation of two temperature fields, one for the fluid and another one
for the solid.

The validation of the full set of thermal boundary conditions is based on a
classical plane channel flow configuration. This academic geometry enables a
rigorous estimation of the formal error by comparison with analytical solutions
for the laminar case. Additionally, for the turbulent case, accurate thermal
statistics from direct numerical simulation can be used as a reference for vali-
dation of the three types of thermal boundary conditions examined here. The
dual approach capability to deal with the demanding case of variable thermal
conductivity has also been addressed. The formal test-case employed is based
on an analytical solution for a purely conductive problem of a heated composite
wall. In this way, the local error introduced by the coupling through the double
solution reconstruction at the interface could be accurately estimated. As a
prototype of computational configuration where the mesh organisation is dis-
connected from the wall geometry, as enables the immersed boundary method,
pipe flow calculations are presented for the validation of imposed temperature
condition in laminar and turbulent regimes.

The analysis of numerical convergence shows that at least second-order ac-
curacy can be reached for the three types of thermal boundary conditions ev-
erywhere in the fluid and solid domain. For the turbulent case, an excellent
overall agreement with the reference data is obtained for the basic thermal
statistics. The most accurate results are obtained for the imposed temperature
case, closely followed by the imposed heat flux case, which is less straightforward
to impose in the context of the immersed boundary method. For the conjugate
heat transfer case, for which extra approximations are required when imposed
with a dual immersed boundary method, a slight loss of accuracy is observed
in the near-wall region when a coarse wall-normal resolution is used. However,
even in this unfavourable situation, the reached accuracy is considered to be
satisfactory and compatible with the standards of direct numerical simulation.

The very small errors introduced by the present immersed boundary method,
even for the demanding case of conjugate heat transfer, opens new perspectives
for the simulation of turbulent flows involving fluid/solid interfaces of complex
geometry. This is the next step of this work which is focused here on the
new concept of dual immersed boundary method with an assessment based on
academic flow configurations (i.e. with simple wall geometry) in order to enable
comparisons with well documented results of reference for a rigorous estimation
of the accuracy. As the conjugate heat transfer problem is posed, present results
suggest that the adaptation to complex geometry is straightforward for the fluid
solver in the sense that it only requires the imposition of the temperature at
the fluid/solid interface. For the conduction solver, the imposition of Neumann
condition is more delicate. The technique proposed here, as a conversion from
a Neumann condition to a soft Dirichlet condition, has led to successful results
for a plane wall for which the wall-normal keeps the same direction. However, it
has to be recognised that this principle has to be generalised and validated for any
wall geometry with a specific reconstruction based on a decomposition between
the wall-normal and wall-parallel heat flow in the solid side. This work is still
in progress with difficulties to find fully comparable documented results in the
context of conjugate heat transfer involving a geometrically complex fluid/solid
interface. Once this milestone successfully reached, the present versatile strategy
may become an efficient way to handle internal flows subjected to fluid/solid
thermal interactions where both the fluid and the solid can be described in
the same computational domain. Depending on the compactness of the device
under study, the use of a Cartesian mesh could even be a very favourable option
and an excellent compromise in terms of simplicity, efficiency and accuracy.

ACKNOWLEDGEMENTS

Gabriel Narváez Campo, recipient of a doctoral scholarship, gratefully acknowledges the support of the National Council of Technological and Scientific Development (Conselho Nacional de Desenvolvimento Científico e Tecnológico, CNPq), Brazil. This work was partially funded by the “Investissements d’Avenir” Program of the France Government (LABEX INTERACTIFS, reference ANR-11-Labex-0017-01). This work was granted access to the HPC resources of the “Núcleo avançado de computação de alto desempenho” (NACAD/COPPE - UFRJ), under the platform LoboC, and the “Très Grand Centre de Calcul” (TGCC) under the allocation A0032A07624 made by the “Le Grand Équipement National de Calcul Intensif” (GENCI). The authors acknowledge the National Laboratory for Scientific Computing (LNCC/MCTI, Brazil) for providing HPC resources of the SDumont supercomputer, which have contributed to the research results reported within this paper. URL: http://sdumont.lncc.br.
References

A. F. POLIAKOV, Wall effect on temperature fluctuations in the viscous sublayer, Teplofizika Vysokikh Temperatur 12 (1974) 328–337.

Y. M. KHABAKHPASHEVA, Experimental investigation of turbulent momentum and heat transfer in the proximity of the wall, in: Heat transfer 1986, Proceedings of The 8th International Heat Transfer Conference, San Francisco, CA - United States, 1984.

Y. L. SINAI, A wall function for the temperature variance in turbulent flow adjacent to a diabatic wall, J. Heat Transfer 109 (1987) 861–865.

N. KASAGI, A. KURODA, M. HIRATA, Numerical investigation of near-wall turbulent heat transfer taking into account the unsteady heat conduction in the solid wall, J. Heat Transfer 111 (1989) 385–392.

T. P. SOMMER, R. M. C. SO, H. S. ZHANG, Heat transfer modeling and the assumption of zero wall temperature fluctuations., J. Heat Transfer 116 (1995) 855–863.

S. SAHA, A. S. H. OOI, H. M. BLACKBURN, Validation criteria for DNS of turbulent heat transfer in pipe flow, in: Procedia Engineering, 10th International Conference on Mechanical Engineering, Dhaka - Bangladesh, 2014, pp. 599–604.

N. KASAGI, N. SHIKAZONO, Contribution of direct numerical simulation to understanding and modelling turbulent transport, in: Proceedings of the Royal Society of London, Oct. 9. Series A: Mathematical and Physical Sciences, volume 451, London, 1995, pp. 257–292.
F. DUCHAINE, A. CORPRON, L. PONS, V. MOUREAU, F. NICOUD, T. POINSOT, Development and assessment of a coupled strategy for conjugate heat transfer with large eddy simulation: Application to a cooled turbine blade, International Journal of Heat and Fluid Flow 30 (2009) 1129 – 1141. doi:https://doi.org/10.1016/j.ijheatfluidflow.2009.07.004.

F. DUCHAINE, N. MAHEU, V. MOUREAU, G. BALARAC, S. MOREAU, Large-Eddy Simulation and Conjugate Heat Transfer Around a Low-Mach Turbine Blade, Journal of Turbomachinery 136 (2013). doi:10.1115/1.4025165, 051015.

J. KIM, Classic ercoftac database at university of surrey, Available in www.ercoftac.org (1989).

N. KASAGI, K. HORIUTI, Y. MIYAKE, T. MIYAUCHI, Y. NAGANO, Test case: fully developed 2-d channel flow, Available in www.thtlab.t.u-tokyo.ac.jp/ (1992).

H. KAWAMURA, K. OHSAKA, H. ABE, K. YAMAMOTO, DNS of turbulent heat transfer in channel flow with low to medium-high Prandtl number fluid, Int. J. Heat and Fluid Flow 19(5) (1998) 482–491.

I. TISELJ, R. BERGANT, B. MAVKO, I. BAJSIĆ, G. HETSRONI, DNS of turbulent heat transfer in channel flow with heat conduction in solid wall, J. Heat Transfer 123 (2001a) 849–857.

I. TISELJ, E. POGREBNYAK, C. LI, A. MOSYAK, G. HETSRONI, Effect of wall boundary condition on scalar transfer in a fully developed turbulent flume, Phys. Fluids 13 (2001b) 1028–1039.

S. PIROZZOLI, M. BERNARDINI, P. ORLANDI, Passive scalars in turbulent channel flow at high reynolds number, J. Fluid Mech. 788 (2016) 614 – 639.
P. ORLANDI, D. SASSUN, S. LEONARDI, Dns of conjugate heat transfer in presence of rough surfaces, International Journal of Heat and Mass Transfer 100 (2016) 250 - 266. URL: http://www.sciencedirect.com/science/article/pii/S0017931015317683. doi:https://doi.org/10.1016/j.ijheatmasstransfer.2016.04.035.

C. FLAGEUL, S. BENHAMANDOUCHE, E. LAMBALLAIS, D. LAURENCE, DNS of turbulent channel flow with conjugate heat transfer: Effect of thermal boundary conditions on the second moments and budgets, Int. J. Comp. Fluid Dyn. 28 (2015) 393-410.

C. FLAGEUL, S. BENHAMANDOUCHE, E. LAMBALLAIS, D. LAURENCE, On the discontinuity of the dissipation rate associated with the temperature variance at the fluid-solid interface for cases with conjugate heat transfer, Int. J. Comp. Fluid Dyn. 111 (2017) 321-328.

S. SATAKE, T. KUNUGI, R. HIMENO, High reynolds number computation for turbulent heat transfer in a pipe flow., in: High Performance Computing, Third International Symposium., Tokyom Japan, 2000, pp. 514-523.

L. REDJEM-SAAD, M. OULD-ROUISS, G. LAURIAT, Direct numerical simulation of turbulent heat transfer in pipe flow., Int. J. Heat and Fluid Flow 28 (2007) 847-861.

S. SAHA, H. M. BLACKBURN, A. S. H. OOI, The influence of pipe length on thermal statistics computed from DNS of turbulent heat transfer., Int. J. Heat and Fluid Flow 32 (2011) 1083-1097.

M. WANG, J. G. GEORGIADI, Conjugate forced convection in crossflow over a cylinder array with volumetric heating, Int. J. Heat and Mass Transfer 39 (1996) 1351-1361.
T. J. YOUNG, K. VAFAI, Convective cooling of a heated obstacle in a channel, Int. J. Heat and Mass Transfer 41 (1998) 3131–3148.

C. S. PESKIN, Flow patterns around heart valves: a digital computer method for solving the equations of motion, Ph.D. thesis, Physiol., Albert Einstein Coll. Med., Univ. Microfilms. Bronx, NY, 1972.

J. KIM, H. CHOI, An immersed-boundary finite-volume method for simulation of heat transfer in complex geometries, KSME international journal 18 (2004) 1026–1035.

R. MITTAL, G. IACCARINO, Immersed boundary methods, Annu. Rev. Fluid Mech. 37 (2005) 239–261.

J. R. PACHECO, A. PACHECO-VEGA, T. RODI, R. E. PECK, Numerical simulations of heat transfer and fluid flow problems using an immersed-boundary finite-volume method on nonstaggered grids, Numerical Heat Transfer, Part B: Fundamentals 48 (2005) 1–24.

G. IACCARINO, S. MOREAU, Natural and forced conjugate heat transfer in complex geometries on cartesian adapted grids, Journal of Fluids Engineering 128 (2006) 838–846.

N. ZHANG, Z. C. ZHENG, S. ECKELS, Study of heat-transfer on the surface of a circular cylinder in flow using an immersed-boundary method, Int. J. Heat and Fluid Flow 29 (2008) 1558–1566.

Z. WANG, J. FAN, K. LUO, K. CEN, Immersed boundary method for the simulation of flows with heat transfer, Int. J. Heat and Mass Transfer 52 (2009) 4510–4518.

S. KANG, G. IACCARINO, F. HAM, Dns of buoyancy-dominated turbulent flows
on a bluff body using the immersed boundary method, J. Comp. Phys. **228**
(2009) 3189–3208.

W. W. REN, C. SHU, J. WU, W. M. YANG, Boundary condition-enforced
immersed boundary method for thermal flow problems with dirichlet temper-
ature condition and its applications, Comp. Fluids **57** (2012) 40–51.

W. REN, C. SHU, W. YANG, An efficient immersed boundary method for ther-
mal flow problems with heat flux boundary conditions, Int. J. Heat and Mass
Transfer **64** (2013) 694–705.

K. LUO, C. MAO, Z. ZHUANG, J. FAN, N. E. L. HAUGEN, A ghost-cell immersed boundary method for the simulations of heat transfer in compressible flows under different boundary conditions part-ii: Complex geometries, Int. J. Heat and Mass Transfer **104** (2017) 98–111.

J. C. SONG, J. AHN, J. S. LEE, An immersed-boundary method for conjugate
heat transfer analysis, Journal of Mechanical Science and Technology **31**
(2017) 2287–2294.

D. DE MARINIS, M. D. DE TULLIO, M. NAPOLITANO, G. PASCAZIO, Improv-
ing a conjugate-heat-transfer immersed-boundary method, Int. J. Numer.
Methods for Heat and Fluid Flow **26** (2016) 1272 – 1288.

L. HE, Closely coupled fluid-solid interface method with moving-average for les based conjugate heat transfer solution, Int. J. Heat and Fluid Flow **79** (2019) 108440.

S. LEE, , W. HWANG, Development of an efficient immersed-boundary method with subgrid-scale models for conjugate heat transfer analysis using large eddy simulation, Int. J. Heat and Mass Transfer **134** (2019) 198–208.
S. DAS, A. PANDA, N. G. DEEN, J. A. M. KUIPERS, A sharp-interface immersed boundary method to simulate convective and conjugate heat transfer through highly complex periodic porous structures, Chemical Eng. Sci. **191** (2019) 1–18.

F. FAVRE, O. ANTEPARA, C. OLIET, O. LEHMKUHL, C. D. PEREZ-SEGARRA, An immersed boundary method to conjugate heat transfer problems in complex geometries. application to an automotive antenna, Appl. Therm. Eng. **148** (2019) 907–928.

S. LAIZET, E. LAMBALLAIS, High-order compact schemes for incompressible flows: A simple and efficient method with quasi-spectral accuracy, J. Comp. Phys. **228** (2009) 5989–6015.

R. GAUTIER, S. LAIZET, E. LAMBALLAIS, A DNS study of jet control with microjets using an immersed boundary method, Int. J. Comp. Fluid Dyn. **28** (2014) 393–410.

M. PILLER, Direct numerical simulation of turbulent forced convection in a pipe, Int. J. Numer. Methods Fluids **49** (2005) 583–602.

S. LARDEAU, E. LAMBALLAIS, J. P. BONNET, Direct numerical simulation of a jet controlled by fluid injection, J. Turbulence **3** (2002).

S. LAIZET, N. LI, Incompact3d: A powerful tool to tackle turbulence problems with up to $o(10^5)$ computational cores, Int. J. Numer. Methods Fluids **67** (2011) 1735–1757.

T. DAIRAY, High fidelity simulation of the aerothermal of an impinging jet, Theses, Université de Poitiers, 2013. URL: https://hal.archives-ouvertes.fr/tel-01101235.
E. Lamballais, R. Vicente Cruz, R. Perrin, Viscous and hyperviscous filtering for direct and large-eddy simulation, Journal of Computational Physics 431 (2021) 110115.

M. B. GILES, Stability analysis of numerical interface conditions in fluid-structure thermal analysis, Int. J. Numer. Methods Fluids 25 (1997) 421–436.

R. GAUTIER, D. BIAU, E. LAMBALLAIS, A reference solution of the flow over a circular cylinder at Re=40., Journal Computers and Fluids 75 (2013) 103–111.

P. Parnaudeau, D. Heitz, E. Lamballais, J. H. Silvestrini, Combination of the immersed boundary method with compact schemes for dns of flows in complex geometry, in: Direct and Large-Eddy Simulation V, Springer, 2004, pp. 581–590.

N. James, D. Biau, M. Pierre, E. Lamballais, Assessment of immersed boundary method toward high accuracy, in: EUROMECH COLLOQUIUM 549 Immersed Boundary Methods : Current Status and Future Research Directions, Leiden, The Netherlands, 2013.

T. DAIRAY, E. LAMBALLAIS, S. BENHAMADOUCHE, Mesh node distribution in terms of wall distance for large-eddy simulation of wall-bounded flows, Flow, Turbulence and Combustion (2017) 1–10.

E. LAMBALLAIS, V. FORTUNÉ, S. LAIZET, Straightforward high-order numerical dissipation via the viscous term for direct and large eddy simulation, J. Comp. Phys. 230 (2011) 3270–3275.

T. DAIRAY, E. LAMBALLAIS, S. LAIZET, J. C. VASSILISCOS, Numerical dissipation vs. subgrid-scale modelling for large eddy simulation, J. Comp. Phys. (2017) 1–10.
S. ISSHIKI, T. OBATA, N. KASAGI, N. HIRATA, An experimental study on heat transfer in a pulsating pipe flow (1st report, time-averaged turbulent characteristics), Bulletin of the Japan Society of Mechanicals Engineers **59** (1993) 2245-2251.

G. K. EL KHOURY, P. SCHALATTER, A. NOORANI, P. F. FISCHER, G. BRETHOUWER, A. V. JOHANSSON, Direct numerical simulation of turbulent pipe flow at moderately high reynolds numbers., Flow, Turbulence and Combustion **91** (2013) 475-495.

V. GNIELINSKI, New equations for heat and mass transfer in the turbulent pipe and channel flow, International Chemical Engineering **16** (1976) 359-368.