# A coupling numerical methodology for weakly transient conjugate heat transfer problems

Guillaume Gimenez, Marc Errera, Dominique Baillis, Y. Smith, F. Pardo

► To cite this version:

Guillaume Gimenez, Marc Errera, Dominique Baillis, Y. Smith, F. Pardo. A coupling numerical methodology for weakly transient conjugate heat transfer problems. International Journal of Heat and Mass Transfer, 2016, 97, pp.975-989. 10.1016/j.ijheatmasstransfer.2016.02.037. hal-01413197

HAL Id: hal-01413197

https://hal.science/hal-01413197

Submitted on 17 Mar 2021

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
# A coupling numerical methodology for weakly transient conjugate heat transfer problems

G. Gimenez<sup>a,b,*</sup>, M. Errera<sup>c</sup>, D. Baillis<sup>a</sup>, Y. Smith<sup>b</sup>, F. Pardo<sup>b</sup>

This study deals with the development of a partitioned coupling strategy at the fluid-solid interface for weakly transient heat transfer problems. The thermal coupling is carried out by an iterative procedure (strong coupling) between a transient solid and a sequence of steady states in the fluid. Continuity of temperature and heat flux is ensured at each coupling time step.

Emphasis is put on the choice of interface conditions at the fluid-solid interface. Two fluid-solid transmission procedures are considered in this paper: Dirichlet-Robin and Neumann-Robin conditions. These conditions are theoretically examined and it is shown that the *Biot* number is a key parameter for determining relevant interface conditions. Stability diagrams are provided in each case and the most effective coupling coefficients are highlighted and expressed. Numerical thermal computations are then performed for two different *Biot* numbers. They confirm the efficiency of the interface conditions in terms of accuracy, stability and convergence. At the end of this paper a comparison between a partitioned and a monolithic approach is presented.

## 1. Introduction

The term conjugate heat transfer is used when the two modes of heat transfer – convection and conduction – are considered simultaneously. CHT procedures are today commonly found in many real-world environments in which accurate heat transfer predictions are needed to design efficient cooling or heating systems. The concept of CHT was first introduced by Perelman in the sixties [1].

Mathematically, a CHT problem is composed of a solid domain and a fluid domain, separated by an interface. Mass, momentum and energy conservation equations are solved in the fluid domain. Temperature and flux are continuous at the interface. Numerically, two main strategies can be employed to solve a CHT problem.

The first one is a monolithic approach. The equations are solved simultaneously, that is, they directly operate on the aggregated fluid and solid equations. In other words, the multi-physics interaction is accounted for in a single mathematical model. There are many monolithic solvers that treat coupled problems in this way in mechanical fluid-structure interactions [2,3] or in CHT [4,5]. The main advantage of the monolithic approach is that the mutual influence between the different domains is taken into account

directly. This approach has also a positive effect on stability, and no coupling iterations are required within a time step. In this paper, FLUENT capabilities will be used to implement that option.

As opposed to monolithic schemes, partitioned methods allow us to use efficient and specialized codes for each domain [6-9]. For partitioned methods, the physical domain is spatially decomposed into partitions and the solution is advanced in time over each partition. This strategy is very popular because it allows the direct use of specific solvers. Calculation codes communicate by exchanging interface conditions at coupling time steps. In this paper, a finite-volume fluid solver (FLUENT) and a finite-element solid solver (ANSYS) will be coupled to implement that option.

Moreover, strategies taking into account characteristic time discrepancies can be developed, in order to have reasonable computational costs. However, because of the sequential fluid/solid strategy, there is no continuity of flux and temperatures. Appropriate methods must be investigated to ensure flux and temperature continuity at the interface and the choice of interface conditions play a crucial role in stability and convergence speed.

In this study, both approaches will be exploited and compared. It is not our intention to discuss the pros and cons of these methods. Emphasis is clearly put on the definition of relevant conditions at a fluid-solid interface in a partitioned method. The monolithic procedure is just used as a means to evaluate and compare the results in terms of accuracy. A monolithic approach intrinsically

<footnote><sup>a</sup>Université de Lyon, CNRS, LaMCoS, INSA Lyon, UMR 5259, F-69621 Villeurbanne Cedex, France</footnote>

<footnote><sup>b</sup>Turbomeca (Safran Group), BP17, 64 511 Bordes Cedex, France</footnote>

<footnote><sup>c</sup>ONERA, The French Aerospace Lab, DMFN, 29 Avenue de la Division Leclerc, 92 322 Châtillon Cedex, France</footnote>

<footnote>* Corresponding author at: Université de Lyon, CNRS, LaMCoS, INSA Lyon, UMR 5259, F-69621 Villeurbanne Cedex, France.</footnote>
## Nomenclature

<table><tr><td>a</td><td>thermal diffusivity (m<sup>2</sup> s<sup>-1</sup>)</td></tr><tr><td>C<sub>p</sub></td><td>heat capacity (J kg<sup>-1</sup> K<sup>-1</sup>)</td></tr><tr><td>h</td><td>convection coefficient (W m<sup>-2</sup> K<sup>-1</sup>)</td></tr><tr><td>l</td><td>flat plate length (m)</td></tr><tr><td>L</td><td>flat plate thickness (m)</td></tr><tr><td>q</td><td>heat flux (W m<sup>-2</sup>)</td></tr><tr><td>T</td><td>temperature (K)</td></tr><tr><td>T<sub>ref</sub></td><td>reference temperature at bottom side (K)</td></tr><tr><td>U</td><td>fluid velocity (m s<sup>-1</sup>)</td></tr><tr><td>t<sub>c</sub></td><td>coupling time step (s)</td></tr><tr><td>Bi</td><td>Biot number (based on the convection coefficient h)</td></tr><tr><td>Bi</td><td>Biot number (based on the relaxation parameter α)</td></tr><tr><td>D</td><td>Fourier number</td></tr><tr><td>α</td><td>coupling relaxation parameter in Robin interface condition (W m<sup>-2</sup> K<sup>-1</sup>)</td></tr><tr><td>δt<sub>s</sub></td><td>solid time discretization (s)</td></tr><tr><td>Δt<sub>c</sub></td><td>coupling period (s)</td></tr></table>

<table><tr><td>Δx</td><td>space discretization (length) (m)</td></tr><tr><td>Δy</td><td>space discretization (width) (m)</td></tr><tr><td>λ</td><td>thermal conductivity (W m<sup>-1</sup> K<sup>-1</sup>)</td></tr><tr><td>ρ</td><td>mass density (kg m<sup>-3</sup>)</td></tr></table>

### Subscripts

<table><tr><td>f</td><td>fluid domain</td></tr><tr><td>s</td><td>solid domain</td></tr><tr><td>∞</td><td>free stream</td></tr></table>

### Superscripts

<table><tr><td>v</td><td>iteration step in coupling period</td></tr><tr><td>n</td><td>temporal index in the solid domain</td></tr><tr><td>( )</td><td>spatial mean quantity</td></tr><tr><td>( )</td><td>unknown quantity</td></tr></table>

ensures temperatures and flux continuity at the interface, and does not require the use of interface conditions and interpolations that may result in stability issues. But obviously, a monolithic approach based on the smallest time characteristic is impractical in complex industrial configurations. On the contrary, a partitioned approach based on an appropriate multiphysics strategy could be a viable method.

In recent years, the behavior (well-posedness, stability, convergence) of interface conditions in a CHT procedure in partitioned techniques has been studied in different ways. The most commonly used method is undoubtedly the normal mode analysis [10–14]. On this basis, a transition of the amplification factor was identified recently and as a result, optimal coefficients have been derived in steady CHT procedures [15]. Other models are available such as the energy method [16] or the matrix analysis [17]. This demonstrates that a great deal of effort has been dedicated to determine robust and efficient fluid-solid interface conditions. Therefore, many interesting stability studies are available and nowadays, steady CHT is applied to a great variety of problems.

On the contrary, the simulation of the transient heat load in solid structures via a fluid-solid coupling approach is much less common but starts to be increasingly used. It must be stressed that steady and unsteady CHT procedures have very little in common. These differences have been highlighted in a recent paper [18]. Unsteady CHT occurs for instance in the prediction of the dynamic thermal conditions in building simulations for modeling building heating, cooling and ventilating flows [19–21]. Accurate knowledge of the transient temperature field in the metallic structures plays also a major role, for example in gas turbine design. Recent fundamental studies of transient aerothermal analysis have already been performed [22–24]. These remarkable studies have been conducted through an entire flight cycle. This cycle is generally divided into ramps and in these ramps, linear distributions of the environment parameters are assumed. Each ramp is just a simplified scenario of “steady” or “unsteady” environment conditions. In steady conditions, internal air system conditions may change. The second case generally reflects severe conditions such as engine acceleration or deceleration.

Typically, the influence of unsteadiness in the fluid domain is negligible and the flow field is thus considered as a sequence of steady states. The solid simulation is treated as unsteady for the whole transient cycle. Many authors have already employed this quasi dynamic method [25,26]. Basically, a strong-coupling

algorithm is used, i.e., additional iterations are introduced to obtain a converged solution at each coupling time step [27,28].

The choice of relevant interface conditions in terms of stability and convergence speed is one of the main issues of partitioned methods. However few studies have been devoted to study interface conditions in transient CHT problems. The objective of this paper is to investigate numerically these conditions in the case of quasi-dynamic conditions. This work is based on Verstraete theory [29], initially developed for steady CHT problems. One of the goals of this paper is to extend the validity of this approach to weakly transient CHT problems.

## 2. Quasi-dynamic coupling strategy

### 2.1. Coupling algorithm

The convective time scale is approximately $\tau_{fluid} = L/\bar{U}$, and the solid diffusive time scale may be expressed as $\tau_{solid} = L^2/\bar{a}$. Hence the solid-fluid time scale ratio is $\tau_{solid}/\tau_{fluid} = LU/\bar{a}$. This ratio is in general very high, for instance in turbomachinery applications. Thus it is possible to assume that the influence of unsteadiness in the fluid domain is negligible and as a result, the flow field may be considered as a sequence of steady states. That is why it is legitimate to couple steady fluid calculations with transient solid calculations. This partitioned coupled method is called quasi dynamic [27], in which each subsystem is represented by an individual solution scheme. Note that if the solid-fluid time scale ratio decreases (for example in natural convection), the unsteady response of the fluid cannot be neglected anymore. In that case, transient calculations have to be performed in both the fluid and solid domains.

The quasi dynamic method is initialized by a steady fluid calculation, performed at the instant $t=0$, with the initial temperature of the solid imposed at the fluid interface. After convergence of this initial fluid calculation, interface conditions are given to the solid, for the beginning of the first coupling period.

Each coupling period, illustrated in Fig. 1 for a time period between $t_c$ and $t_c + \Delta t_c$, is composed of 4 steps.

The coupling period $t_{\Delta c}$ is divided into several solid time increments $\delta t_s(\Delta t_c = n\delta t_s)$.

Each coupling period is repeated until continuity of fluid and solid fluxes and temperatures at the interface, at every coupling
<figure><img src="image_2.png" /><figcaption>Fig. 1. Quasi-dynamic coupling algorithm.</figcaption></figure>

time step. At step ⑤ convergence criteria should be verified, defined at the vth iteration as:

$$ \frac{|\bar{q}_s^n - \bar{q}_f^n|}{|\bar{q}_s^v| + |\bar{q}_f^v|} < \varepsilon \qquad (1) $$

$$ \frac{|\bar{T}_s^v - \bar{T}_f^v|}{|\bar{T}_s^v| + |\bar{T}_f^v|} < \varepsilon' \qquad (2) $$

It is noticed that for quasi-adiabatic cases, the denominator of Eq. (1) may become or get close to zero. In any event, such case can only occur if both the fluid and solid heat flux are small, and thus continuity of heat flux is guaranteed.

At the first iteration of a coupling period between $t_c$ and $t_c + \Delta t_c$, interface conditions imposed on the solid side are constant, given by the fluid domain at the coupling time step $t_c$. Then at subsequent iterations within the same coupling period, interface conditions imposed on the solid side, at each time increment $\delta t_s$, are linearly interpolated between the converged fluid state at $t_c$ and the last known state at $t_c + \Delta t_c$.

## 2.2. Frequency of exchanges

The numerical strategy is driven by the need to analyze weakly transient conjugate heat transfer problems over a long period of time, at a reasonable computing cost. This period of time is in general described by a cycle. A typical cycle is composed of various conditions defined by a set of ramp points. Each point represents changes in the operating conditions defined either by unsteady ramps or by steady ramps. The ramp points are linked with operating conditions supplied for all the steady operating conditions and a linear variation of the boundary conditions are assumed between two ramp points. It is only at ramp points that a thermal coupling may take place and then these points can also be regarded as coupling times. A cycle example will be presented and illustrated in the description of the test case.

The duration of a coupling period $\Delta t_c$ is the time between exchanges of boundary conditions at the interface between the fluid and solid domains. In this study the frequency of exchanges is constant in the whole calculation.

## 2.3. Fluid-solid transmission conditions

### 2.3.1. General interface conditions

Several interface conditions can be imposed on the solid side.

Robin-type interface condition at time step $n(t_c < t \le t_c + \Delta t_c)$ and iteration $v + 1$ is defined as:

$$ \hat{q}_{s}^{n,v+1} = -q_{f}^{n,v} + \alpha_{f}^{n,v}(T_{f}^{n,v} - \hat{T}_{s}^{n,v+1}) \qquad (3) $$

where $\alpha_f^{n,v}$ is a relaxation parameter.

Convection boundary conditions can be imposed on the fluid side:

$$ \hat{q}_{f}^{n,v+1} = h^{n,v}(T_{ref}^{n,v} - \hat{T}_{s}^{n,v+1}) \qquad (4) $$

with

$$ h^{n,v} = \frac{q_{f}^{n,v}}{T_{ref}^{n,v} - T_{ref}^{n,v}} \qquad (5) $$

and $T_{ref}^{n,v}$ a reference temperature.

A convection condition is a particular case of Robin condition with:

$$ \left\{ \begin{array}{l} h^{n,v} = \alpha_f^{n,v} \\ T_{ref}^{n,v} = -\frac{q_f^{n,v}}{\alpha_f^{n,v}} + T_f^{n,v} \end{array} \right. \qquad (6) $$

Dirichlet boundary conditions (temperature imposed) and Neumann boundary condition (flux imposed) are particular cases of Robin boundary conditions:

$$ \alpha_f^{n,v} \to \infty : \quad \hat{T}_s^{n,v+1} = T_f^{n,v} \qquad (7) $$

$$ \alpha_f^{n,v} = 0 : \quad \hat{q}_s^{n,v+1} = -q_f^{n,v} \qquad (8) $$

Similar interface conditions can be imposed on the fluid side, except the convection condition because a convection coefficient cannot be defined in the solid side.

Robin condition on the fluid side is written:

$$ \hat{q}_{f}^{n,v+1} = -q_{s}^{n,v+1} + \alpha_{s}^{n,v+1}(T_{s}^{n,v+1} - \hat{T}_{f}^{n,v+1}) \qquad (9) $$

Dirichlet and Neumann interface conditions on the fluid side are respectively:

$$ \hat{T}_{f}^{n,v+1} = T_{s}^{n,v+1} \qquad (10) $$

$$ \hat{q}_{f}^{n,v+1} = q_{s}^{n,v+1} \qquad (11) $$

### 2.3.2. Stability analysis and appropriate interface quantities

Stability analysis of partitioned steady thermally coupled problems has been studied by Giles [10] by applying Godunov and Ryabenkii theory [30]. In most practical cases, for stability reasons, it is recommended to use Dirichlet-Neumann interface conditions, in other words the fluid domain is supplemented with a Dirichlet condition (temperature coming from the solid) and the Neumann condition (heat flux coming from the fluid) is imposed on the solid side. These conclusions are often used in CHT
literature. However Giles' pioneering work is based on hypothesis which are not always verified for every configuration.

Robin conditions have many attractive features because they can be formulated in such a way that the associated local problem is well posed. Moreover they locally provide an interface stiffness that considerably stabilize the coupled problem, if the coefficients are well chosen [15].

Another type of problem, which is quite different, is the analysis of the transient response of a solid to a change in the operating conditions over a long period of time. It is the major problem dealt here. Recently, a quasi dynamic process has been used to obtain consistent thermal solutions [27]. The fluid may be assumed to adjust instantaneously to changes and as result, a sequence of fluid steady states is coupled with a transient solid computation.

Robin conditions have already been used in transient problems [27], but theoretical analysis of relaxation parameters influence has not been carried out. The objective and contribution of this paper is to investigate numerically the influence of interface conditions in weakly transient CHT problems, and to find suitable choices of the relaxation parameters.

### 2.3.3. Transient and weakly transient in CHT analysis

A transient state is a state of non-equilibrium, when the solid temperatures are still changing with time. The Fourier number (D) and the Biot number (Bi) are dimensionless numbers that characterize transient conduction problems.

The Fourier number D is defined as:

$$D = \frac{a\Delta t}{\Delta^2} \qquad (12)$$

where Δ is a characteristic length of penetration of heat into the solid. The Fourier number is a dimensionless measure of time used in transient conduction problems. It is the dimensionless time for a temperature change to occur.

The Biot number Bi is given by:

$$Bi = \frac{\text{external convective conductance}}{\text{internal diffusivity conductance}} = \frac{h}{\frac{\lambda_s}{L}} \qquad (13)$$

This number determines whether or not the temperatures inside a body will vary significantly in space, while the body heats or cools over time, from a thermal gradient applied to its surface. This number plays a fundamental role in conduction problems involving surface convection effects.

At moderate or small Bi (smaller than 1) problems are generally thermally simple, due to relatively uniform temperature fields inside the body although this temperature may be changing, as heat passes into the solid from the surface. In other words, the "solid resistance" of the body is small relative to the "convective resistance" and as a result most of the temperature drop is in the fluid and the temperature gradients may be negligible in the solid.

In contrast, Biot numbers much larger than 1 characterize more difficult problems due to non-uniformity of temperature fields within the solid. If the material is thermally insulating (poorly conductive), such as PVC, the interior resistance to heat flow will exceed that of the fluid boundary. In this case, again, the Biot number will be greater than one and more complicated heat transfer equations for "transient heat conduction" will be required to describe the time-varying and non-spatially-uniform temperature field within the material body.

In what follows, small and high Biot numbers will be considered, in the context of weakly transient problems. In other words, the transients will weakly depend on the Fourier number. Another way to express this, is to consider that the slope of the solid heat flux can reasonably be assumed to be equal everywhere to $\frac{\lambda_s}{L}$. During a transient flight cycle for instance, this assumption is valid

during steady ramps, when the environment parameters of two ramp points are almost identical even if metal temperature may change or thermal gradients are not negligible.

This assumption allows us to introduce interface conditions and interface coupling coefficients that do not depend on the Fourier number but just on the Biot number. Verstraete [29] developed a physics-based approach for steady CHT and the goal of this paper is to assess whether, and to what extent, this theory can be applied to weakly transient problems.

## 3. Theoretical analysis: interface conditions

Several boundary conditions at the interface can be chosen. In the more general case, Robin boundary conditions (Eqs. (3)-(9)) are imposed on the fluid and solid sides [15]. If we take into account that the fluid and solid heat fluxes have opposite directions on their shared interface, we obtain:

$$\left\{ \begin{aligned} \hat{q}_s^{n,v+1} &= -q_f^{n,v} + \alpha_f^{n,v}(T_f^{n,v} - \hat{T}_s^{n,v+1}) \\ -\hat{q}_f^{n,v+1} &= q_s^{n,v+1} + \alpha_s^{n,v+1}(T_s^{n,v+1} - \hat{T}_f^{n,v+1}) \end{aligned} \right. \qquad (14)$$

where $\alpha_f^{n,v}$ and $\alpha_s^{n,v+1}$ are relaxation parameters.

It will be shown that if the fluid and solid domains converge, then Eqs. (1) and (2) are satisfied.

Let us assume that the solid and the fluid domains converge. This implies:

$$\left\{ \begin{aligned} T_s^v &= T_s^{v+1} = T_s^\infty \\ q_s^v &= q_s^{v+1} = q_s^\infty \\ \alpha_s^v &= \alpha_s^{v+1} = \alpha_s^\infty \end{aligned} \right. \qquad (15)$$

and

$$\left\{ \begin{aligned} T_f^v &= T_f^{v+1} = T_f^\infty \\ q_f^v &= q_f^{v+1} = q_f^\infty \\ \alpha_f^v &= \alpha_f^{v+1} = \alpha_f^\infty \end{aligned} \right. \qquad (16)$$

Thus Eq. (14) is rewritten:

$$\left\{ \begin{aligned} q_s^\infty &= -q_f^\infty + \alpha_f^\infty (T_f^\infty - T_s^\infty) \\ -q_f^\infty &= q_s^\infty + \alpha_s^\infty (T_s^\infty - T_f^\infty) \end{aligned} \right. \qquad (17)$$

By adding Eq. (17), we obtain:

$$(\alpha_f^\infty - \alpha_s^\infty)(T_f^\infty - T_s^\infty) = 0. \qquad (18)$$

If $\alpha_f^\infty \neq \alpha_s^\infty$, it can be deduced:

$$T_s^\infty = T_f^\infty. \qquad (19)$$

and then:

$$q_s^\infty = -q_f^\infty. \qquad (20)$$

As a result, in the general Robin-Robin transmission conditions (Eq. (14)), the solution obtained at convergence is independent of the coupling coefficients used at the interfaces and moreover, both temperature and flux continuity is satisfied, independently of the relaxation parameters.

Thus Eqs. (1) and (2) (flux and temperature continuity at the interface) are equivalent to the following criterion (flux convergence in solid), which is adopted in the numerical simulations:

$$\frac{|\overline{q_s^{v+1}} - \overline{q_s^v}|}{|q_s^{v+1}| + |q_s^v|} < \epsilon'' \qquad (21)$$

Note that Eq. (21) leads to the convergence of flux and temperatures in both fluid and solid domains.
Robin–Robin conditions contain a rich diversity of interfacial methods that include a large family of disparate schemes. We have therefore preferred to focus on the study of two specific and opposite transmission conditions: Dirichlet in the fluid domain, and then Neumann in this same medium. The first condition is widely used in the literature. In contrast, in the solid medium, a general Robin condition will be considered and studied numerically in terms of stability and convergence. We thus propose to begin by focusing on these two diametrically different conditions to provide a better understanding of these particular schemes before other numerical issues are dealt with.

Thus the objective is limited to find the optimal relaxation parameter $\alpha_f$ in Robin condition imposed on the solid side, in terms of stability and convergence speed. Verstraete [29] developed a stability theory for stationary CHT problems based on a 1D model, in which the Biot number determines the optimal choice of interface quantities.

### 3.1. Dirichlet–Robin conditions

Temperature imposed on the fluid side is a particular case of a Robin condition with $\alpha_s = \infty$.

The interface condition becomes:

$$ egin{cases} \dot{q}_s^{n,v+1} = -q_f^{n,v} + \alpha_f^{n,v} (T_f^{n,v} - \hat{T}_s^{n,v+1}) \\ \hat{T}_f^{n,v+1} = T_s^{n,v+1} \end{cases} \quad (22) $$

Note that in the literature, temperature is usually imposed on the fluid side.

Verstraete shows in [29] that:

$$ T_f^v = T_f^\infty + \left[ \frac{\tilde{B}_i - B_i}{\tilde{B}_i + 1} \right]^v \Delta T_0. \quad (23) $$

With:

$$ B_i = \frac{hL}{\lambda_s} \quad (24) $$

$$ \tilde{B}_i = \frac{\alpha_f^v L}{\lambda_s} \quad (25) $$

and $\Delta T_0$ the initialization error, defined by:

$$ \Delta T_0 = T_f^0 - T_f^\infty. \quad (26) $$

#### 3.1.1. Influence on stability

Thus the calculation is stable if:

$$ \left| \frac{\tilde{B}_i - B_i}{\tilde{B}_i + 1} \right| < 1 \quad (27) $$

Then:

$$ -1 < \frac{\tilde{B}_i - B_i}{\tilde{B}_i + 1} < 1. \quad (28) $$

Leading to:

$$ \text{stable for } \tilde{B}_i > \frac{B_i - 1}{2}. \quad (29) $$

It can be written as:

$$ \text{stable for } \alpha_f^v > \frac{h}{2} - \frac{\lambda_s}{2L} \quad (30) $$

Fig. 2 highlights two main zones separated by the line defined by the equation $\tilde{B}_i = \frac{B_i - 1}{2}$.

If $B_i < 1$, the method is stable regardless the choice of relaxation parameter $\alpha_f^v > 0$.

If $B_i > 1$, the method is stable provided an appropriate choice of $\alpha_f^v$, such as $\tilde{B}_i > \frac{B_i - 1}{2}$ or $\alpha_f^v > \frac{h}{2} - \frac{\lambda_s}{2L}$.

<figure><img src="image_3.png" /><figcaption>Fig. 2. Stability diagram (temperature imposed on fluid, Robin imposed on solid).</figcaption></figure>
### 3.1.2. Influence on convergence speed

From Eq. (27), we can notice that if $\tilde{Bi} = Bi$ (equivalent to $\alpha_f^* = h$), then $|\frac{\tilde{Bi}-Bi}{Bi+1}| = 0$. It means that no iteration is required to converge. However, in practice, the convection coefficient is not easy to determine in complex industrial configurations.

The smaller $|\frac{Bi-Bi}{Bi+1}|$, the faster the convergence speed. On the contrary, the higher $|\frac{Bi-Bi}{Bi+1}|$ the slower the convergence speed (as long as $|\frac{Bi-Bi}{Bi+1}| < 1$).

$|\frac{Bi-Bi}{Bi+1}|$ with respect to $\tilde{Bi}$ and $Bi$ is plotted on Fig. 3, representing stability and convergence speed. Blue areas are stable and red areas are unstable. It is noticed that stability increases for small Biot numbers.

Two examples of $|\frac{Bi-Bi}{Bi+1}|$ as a function of the relaxation parameter $\alpha_f$ are shown in Figs. 4 and 5, for respectively constant $Bi = 0.108$ and $Bi = 13.5$.

For a fixed Biot number, limits of $|\frac{Bi-Bi}{Bi+1}|$ with a zero or infinite relaxation parameter are:

$$ \left\{ \begin{array}{l} |\frac{\tilde{Bi}-Bi}{Bi+1}| \rightarrow |Bi| \quad \text{when } \tilde{Bi} \rightarrow 0 \\ |\frac{\tilde{Bi}-Bi}{Bi+1}| \rightarrow 1 \quad \text{when } \tilde{Bi} \rightarrow \infty \end{array} \right. \qquad (31) $$

If $Bi < 1$ (Fig. 4), a relatively low calculation time is obtained for $\tilde{Bi} < Bi$ (equivalent to $\alpha_f < h$). But the calculation time increases when $\alpha_f$ increases, for $\tilde{Bi} > Bi$. It confirms that the CHT calculation is stable for all relaxation parameters such that $\alpha_f > 0$. It is noticed that Dirichlet–Neumann interface conditions ($\alpha_f = 0$), often used in literature, offer a good convergence speed. However it is better to use a relaxation parameter close to the convection coefficient $h$.

If $Bi > 1$ (Fig. 5), the calculation may become unstable when $\tilde{Bi} < Bi$, and strongly increases when $\tilde{Bi} > Bi$. When $Bi > 1$, imposing temperature on the fluid side is not the most appropriate choice.

<figure><img src="image_4.png" /><figcaption>Fig. 4. $|\frac{Bi-Bi}{Bi+1}|$ with respect to $\alpha_f$ (for constant Biot number $Bi = 0.108$).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 5. $|\frac{Bi-Bi}{Bi+1}|$ with respect to $\alpha_f$ (for constant Biot number $Bi = 13.5$).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 3. $|\frac{Bi-Bi}{Bi+1}|$ with respect to $Bi$ and $Bi$.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Fig. 6. Stability diagram (flux imposed on fluid, Robin condition imposed on solid).</figcaption></figure>

### 3.2. Neumann–Robin conditions

The heat flux imposed on the fluid side is a particular case of Robin condition with $\alpha_s = 0$.

The interface condition becomes:

$$ \begin{cases} \bar{q}_s^{n,v+1} = -q_f^{n,v} + \alpha_f^{n,v}(T_f^{n,v} - \hat{T}_s^{n,v+1}) \\ \bar{q}_f^{n,v+1} = q_s^{n,v+1} \end{cases} \quad (32) $$

Verstraete shows in [29] that:

$$ q_f^v = q_f^\infty + \left[ \frac{\tilde{B}_i - B_i}{\tilde{B}_i(\tilde{B}_i + 1)} \right]^v q_0 \quad (33) $$

With $\Delta q_0$ the initialization error, defined by:

$$ \Delta q_0 = q_f^0 - q_f^\infty. \quad (34) $$

#### 3.2.1. Influence on stability

Thus the calculation is stable if:

$$ \left| \frac{\tilde{B}_i - B_i}{B_i(\tilde{B}_i + 1)} \right| < 1 \quad (35) $$

This criterion is similar to Dirichlet–Robin criterion, but multiplied by $\frac{1}{\tilde{B}_i}$.

Then:

$$ -1 < \frac{\tilde{B}_i - B_i}{B_i(\tilde{B}_i + 1)} < 1 \quad (36) $$

Leading to:

$$ \begin{cases} \text{if } B_i < 1 \rightarrow \text{stable for } \tilde{B}_i < \frac{2B_i}{1-\tilde{B}_i} \\ \text{if } B_i \ge 1 \rightarrow \text{always stable} \end{cases} \quad (37) $$

It can be written as:

$$ \begin{cases} \text{if } B_i < 1 \rightarrow \text{stable for } \alpha_f < \frac{2h}{1-\tilde{B}_i} \\ \text{if } B_i \ge 1 \rightarrow \text{always stable} \end{cases} \quad (38) $$

Fig. 6 highlights two main zones separated by the hyperbola defined by the equation $\tilde{B}_i = \frac{2B_i}{1-\tilde{B}_i}$.

If $B_i < 1$, the method is stable provided an appropriate choice of $\alpha_f^v$, such as $\tilde{B}_i < \frac{2B_i}{1-\tilde{B}_i}$.

If $B_i > 1$, the method is stable regardless the choice of the relaxation parameter $\alpha_f^v$.

#### 3.2.2. Influence on convergence speed

$|\frac{\tilde{B}_i - B_i}{B_i(\tilde{B}_i + 1)}|$ with respect to $\tilde{B}_i$ and $B_i$ is plotted on Fig. 7, representing stability and convergence speed. Blue areas are stable and red areas are unstable. It is noticed that stability increases for high Biot numbers.

<figure><img src="image_8.png" /><figcaption>Fig. 7. $|\frac{\tilde{B}_i - B_i}{B_i(\tilde{B}_i + 1)}|$ with respect to $\tilde{B}_i$ and $B_i$.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** $|\frac{\tilde{B}_i-B_i}{B_i(\tilde{B}_i+1)}|$ with respect to $\alpha_f$ (for constant Biot number $Bi = 0.108$).</figcaption></figure>

$|\frac{\tilde{B}_i-B_i}{B_i(\tilde{B}_i+1)}|$ is shown on Fig. 8 and Fig. 9 as a function of the relaxation parameter $\alpha_f$, for respectively constant $Bi = 0.108$ and $Bi = 13.5$.

For a fixed Biot number, limits of $|\frac{\tilde{B}_i-B_i}{B_i(\tilde{B}_i+1)}|$ with a zero or infinite relaxation parameter are:

$$ \left\{ \begin{array}{l} \left| \frac{\tilde{B}_i - B_i}{B_i(\tilde{B}_i + 1)} \right| \to 1 \quad \text{when } \tilde{B}_i \to 0 \\ \\ \left| \frac{\tilde{B}_i - B_i}{B_i(\tilde{B}_i + 1)} \right| \to \left| \frac{1}{\tilde{B}_i} \right| \quad \text{when } \tilde{B}_i \to \infty \end{array} \right. \qquad (39) $$

If $Bi < 1$ (Fig. 8), the calculation time is very high when $\tilde{B}_i < B_i$, and the method becomes rapidly unstable when $\tilde{B}_i^l > B_i$. When $Bi < 1$, imposing flux on the fluid side is not the best choice.

If $Bi > 1$ (Fig. 9), the calculation time increases significantly when $\alpha_f$ decreases, for $\tilde{B}_i < B_i$, and becomes infinitely high when $\tilde{B}_i = 0$. However it remains very low if $\tilde{B}_i < B_i$. So when $Bi > 1$, a

<figure><img src="image_10.png" /><figcaption>**Fig. 9.** $|\frac{\tilde{B}_i-B_i}{B_i(\tilde{B}_i+1)}|$ with respect to $\alpha_f$ (for constant Biot number $Bi = 13.5$).</figcaption></figure>

<table><caption>**Table 1**
Summary table of optimal interface conditions.</caption><thead><tr><th>Biot number <i>Bi</i></th><th>0</th><th>1</th><th>&infin;</th></tr></thead><tbody><tr><td>Fluid interface condition</td><td>Temperature</td><td>Flux</td><td></td></tr><tr><td>Solid interface condition</td><td>Robin (&alpha;<sub>f</sub><sup>l</sup> &le; <i>h</i>)</td><td>Robin (&alpha;<sub>f</sub><sup>l</sup> &ge; <i>h</i>)</td><td></td></tr></tbody></table>

relaxation parameter such as $\alpha_f^l > h$ (or $\tilde{B}_i > B_i$) is an excellent choice.

### 3.3. Summary

In steady conditions or in a weakly transient process, the Biot number determines the appropriate choice of interface conditions (see Table 1).

If $Bi < 1$, the best choice is to impose the solid temperature on the fluid side, and a Robin condition on the solid side with ideally $\alpha_f^l = h$, or $\alpha_f^l < h$. In practice $h$ may not be known, so a small $\alpha_f^l$ is a good choice. Physically, if $Bi \ll 1$, temperature is almost uniform in the solid. So it is logical to impose this temperature on the fluid side.

If $Bi > 1$, the best choice is to impose the solid heat flux on the fluid side, and a Robin condition on the solid side with ideally $\alpha_f^l = h$, or $\alpha_f^l > h$. In practice a high $\alpha_f^l$ is a good choice. Physically, if $Bi \gg 1$, temperature gradients in the solid are not negligible, so the heat flux is imposed on the fluid side.

### 4. Operating conditions

The test case is a flat plate cooled by convection on its upper face and heated on its lower face (see Fig. 10). Convection boundary conditions applied on the underside are time dependent. The fluid-solid coupling interface is the line $y = 0$. The temporal evolution of temperature at several points of coupling interface will be studied.
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Flat plate cooled by convection.</figcaption></figure>

<table><caption>**Table 2**
Fluid thermophysical properties (air).</caption><thead><tr><th>Notation</th><th>Name</th><th>Value</th></tr></thead><tbody><tr><td>ρ<sub>f</sub></td><td>Fluid density</td><td>1.225 kg m<sup>-3</sup></td></tr><tr><td>μ<sub>f</sub></td><td>Fluid dynamic viscosity</td><td>1.7894e-5 kg m<sup>-1</sup> s<sup>-1</sup></td></tr><tr><td>λ<sub>f</sub></td><td>Fluid thermal conductivity</td><td>0.0242 W m<sup>-1</sup> K<sup>-1</sup></td></tr><tr><td>C<sub>pf</sub></td><td>Fluid heat capacity</td><td>1006.43 J kg<sup>-1</sup> K<sup>-1</sup></td></tr></tbody></table>

<table><caption>**Table 3**
Solid thermophysical properties (solid 1).</caption><thead><tr><th>Notation</th><th>Name</th><th>Value</th></tr></thead><tbody><tr><td>ρ<sub>s</sub></td><td>Solid density</td><td>1200 kg m<sup>-3</sup></td></tr><tr><td>λ<sub>s</sub></td><td>Solid thermal conductivity</td><td>0.16 W m<sup>-1</sup> K<sup>-1</sup></td></tr><tr><td>C<sub>ps</sub></td><td>Solid heat capacity</td><td>1400 J kg<sup>-1</sup> K<sup>-1</sup></td></tr></tbody></table>

<table><caption>**Table 4**
Solid thermophysical properties (solid 2).</caption><thead><tr><th>Notation</th><th>Name</th><th>Value</th></tr></thead><tbody><tr><td>ρ<sub>s</sub></td><td>Solid density</td><td>1200 kg m<sup>-3</sup></td></tr><tr><td>λ<sub>s</sub></td><td>Solid thermal conductivity</td><td>40 W m<sup>-1</sup> K<sup>-1</sup></td></tr><tr><td>C<sub>ps</sub></td><td>Solid heat capacity</td><td>1400 J kg<sup>-1</sup> K<sup>-1</sup></td></tr></tbody></table>

<table><caption>**Table 5**
Characteristic time scales.</caption><thead><tr><th></th><th>Solid time step δt<sub>s</sub></th><th>Mesh Fourier <sup>λ<sub>f</sub>δt<sub>s</sub></sup>&frasl;<sub>ρ<sub>f</sub>C<sub>ps</sub>Δy<sub>i</sub><sup>2</sup></sub></th><th>Coupling period Δt<sub>c</sub></th></tr></thead><tbody><tr><td>Table 3</td><td>60 s</td><td>5.71</td><td>600 s</td></tr><tr><td>Table 4</td><td>60 s</td><td>1428.57</td><td>600 s</td></tr></tbody></table>

This simple test-case permits the use of the monolithic method with reasonable cost. Monolithic results are used for reference to test accuracy of partitioned methods.

Thermophysical properties of the fluid (air) are detailed in Table 2.

With the parameters of this application test, Reynolds number is $Re = \frac{\rho_f U_{\infty} l}{\mu_f} = 1.3410^6$.

Two different solids are used with different thermal conductivities (Tables 3 and 4), in order to have two different Biot numbers ($Bi < 1$ and $Bi > 1$). Other thermophysical properties (solid density and solid heat capacity) are identical for solid 1 and solid 2.

Convection conditions are applied at the bottom face ($y = -12$ mm). The convection coefficient is constant ($h = 500$ W m$^{-2}$ K$^{-1}$). However, the reference temperature ($T_{ref}$) is time dependent, as shown in Fig. 11. The temperature evolution is simulated over a long period of time ($t_{max} = 10,800$ s). The computation is divided into 18 coupling periods of 600 s (Fig. 11). Solid time step is $\delta t_s = 60$ s. Initial temperature in solid domain is $T_i = 300$ K.

The characteristic time scales of the numerical cases corresponding to Tables 3 and 4 are summarized in Table 5.

Let us recall that two different multiphysics solutions have been adopted in this paper:

- A reference solution provided by FLUENT in the framework of a monolithic approach

- A partitioned approach provided by FLUENT/ANSYS through WORKBENCH multiphysics platform and driven by PYTHON scripts.

<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Temporal evolution of reference temperature on bottom face ($y = -12$ mm).</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Fluid and solid mesh.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. 13.** Wall y$^+$ along the flat plate.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. 14.** Convection coefficient along the flat plate.</figcaption></figure>

A structured mesh is used for both the fluid and solid domains (Fig. 12). Fluid and solid meshes are not coincident at the interface (with 150 cells for the fluid along the interface, again 130 cells for the solid). Moreover the finite-volume flow solver is based on a cell-centered approach, while the finite-element solver defines the temperature at the vertices of the elements, requiring the use of spatial linear interpolation.

The width of the solid cell is $\Delta y_s = 10^{-3}$ m.

The k-ω SST turbulence model is employed, with a refined fluid mesh in the boundary layer. Criterion of $y^+$ in the first fluid cell is $y^+ \approx 1$ (Fig. 13).

A monolithic calculation with FLUENT gives the interface convection coefficient $h$, defined by Eq. (5) with a free stream reference temperature on the top face. The fluid mesh is refined near the leading edge, in order to capture high gradients of $h$ in this area (Fig. 14). The solid mesh is also refined near the leading edge, so that the solid can «see» the variations of the heat transfer coefficient near the leading edge, coming from the fluid side.

## 5. Numerical results

The influence of interface conditions on accuracy, stability and convergence speed is studied.

As mentioned previously, two types of boundary conditions are investigated at the coupling interface:

- Dirichlet conditions (temperature) imposed on the fluid side, Robin conditions imposed on the solid side;

- Neumann conditions (heat flux) imposed on the fluid side, Robin conditions imposed on the solid side.

As a result, coupling coefficients are used on the solid side only.

According to the theoretical analysis (see Section 3. Theoretical analysis: interface conditions), the choice of the interface conditions in terms of stability and convergence speed is dependent on the convection coefficient $h$. Thus the significant variations of this convection coefficient along the flat plate makes the choice of a local and time dependent relaxation parameter more relevant than a constant one. In this study, relaxation parameters are most often chosen proportional to the convection coefficient $h$. They also can be proportional to the instability critical relaxation parameters (Eqs. (30)–(38)), or equal to $10^6$ and $10^{-6}$ to respectively nearly impose the temperature and the heat flux on the solid side (Eqs. (7) and (8)).

### 5.1. Numerical influence on stability and convergence speed

The calculation time is nearly proportional to the number of internal iterations (see Section 2.1 Coupling algorithm).

In this study, convergence criterion $\epsilon'' = 10^{-4}$ is adopted (Eq. (21)).

#### 5.1.1. Analysis of Dirichlet–Robin conditions

5.1.1.1. Biot number $Bi > 1$. The number of iterations for each coupling period is plotted in Fig. 15. The calculation is performed with physical parameters of Table 3, leading to $Bi > 1$. Interface conditions are temperature imposed on the fluid side, and the
<figure><img src="image_16.png" /><figcaption>Fig. 15. Number of iterations for each coupling period for the relaxation parameters $\alpha_f = 1.7 * (\frac{h}{2} - \frac{\lambda_s}{2L})$, $\alpha_f = h$, $\alpha_f = 2h$ and $\alpha_f = 3h$ (Table 3, Robin condition imposed on solid, temperature imposed on fluid).</figcaption></figure>

Robin condition imposed on the solid side, for several relaxation parameters. It is noticed that the number of iterations is higher in transient phases than in steady phases.

The same tendencies are observed numerically and theoretically. Theoretically, the stability limit is obtained with $\alpha_f = \frac{h}{2} - \frac{\lambda_s}{2L}$ (Eq. (30)). Numerically, instability is observed below the critical relaxation parameter $\alpha_f = 1.7(\frac{h}{2} - \frac{\lambda_s}{2L})$. A discussion about the differences observed between Verstraete theory and numerical simulations is proposed in Section 5.1.3. However in practice, if $Bi > 1$, Neuman-Robin interface conditions should be preferably used to ensure unconditional stability (Eq. (38)).

Numerically, the calculation time reaches a minimum with $\alpha_f = h$, and then increases when $\alpha_f$ becomes smaller or larger than $h$, as expected theoretically.

Thus theoretical analysis can guide the appropriate choice of the relaxation parameter $\alpha_f$.

5.1.1.2. Biot number $Bi < 1$. The number of iterations for each coupling period is plotted on Fig. 16. The calculation is carried out with the physical parameters of Table 4, corresponding to $Bi < 1$. The interface conditions are temperature imposed on the fluid side, and a Robin condition imposed on the solid side.

<figure><img src="image_17.png" /><figcaption>Fig. 16. Number of iterations for each coupling period for the relaxation parameters $\alpha_f = 10^{-6}$, $\alpha_f = 0.5h$, $\alpha_f = h$ and $\alpha_f = 3h$ (Table 4, Robin condition imposed on solid, temperature imposed on fluid).</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Fig. 17.** Number of iterations for each coupling period for relaxation parameters $\alpha_f = 0.5h$, $\alpha_f = h$, $\alpha_f = 2h$ and $\alpha_f = 1e6$ (Table 3, Robin condition imposed on solid, flux imposed on fluid).</figcaption></figure>

It is noticed that the calculation is unconditionally stable, regardless the relaxation parameter $\alpha_f$. Results are consistent with theoretical analysis (Fig. 2). Best convergence speed is also obtained with $\alpha_f = h$.

### 5.1.2. Analysis of Neumann–Robin conditions

5.1.2.1. Biot number $Bi > 1$. The number of iterations for each coupling period is plotted on Fig. 17, with physical parameters of Table 3, corresponding to $Bi > 1$. Interface conditions are the heat

flux imposed on the fluid side, and a Robin condition imposed on the solid side.

The CHT calculation is always stable, for all positive relaxation parameters $\alpha_f$, as predicted by the theoretical analysis (Fig. 6). A fast convergence speed is again obtained with $\alpha_f = h$.

5.1.2.2. Biot number $Bi < 1$. The number of iterations for each coupling period is plotted on Fig. 18. Calculation is performed with physical parameters of Table 4, corresponding to $Bi < 1$. Interface

<figure><img src="image_19.png" /><figcaption>**Fig. 18.** Number of iterations for each coupling period for relaxation parameters $\alpha_f = 0.5h$, $\alpha_f = h$, $\alpha_f = 0.8 * (2h/(1-Bi))$ and $\alpha_f = 1.1 * (2h/(1-Bi))$ (Table 4, Robin condition imposed on solid, flux imposed on fluid).</figcaption></figure>
conditions are the heat flux imposed on the fluid side, and a Robin condition imposed on the solid side.

The calculation time reaches a minimum for $\alpha_f = h$. Theoretically, stability limit is obtained with $\alpha_f = \frac{2h}{1-Bi}$ (Eq. (38)). Numerically, the calculation is unstable above $\alpha_f = 1.1 * (\frac{2h}{1-Bi})$. Therefore, stability and convergence speed can be well estimated by the theoretical analysis.

### 5.1.3. Summary: relevance of a steady model

The numerical results presented in this paper in the case of weakly transient CHT problems show a good agreement with the results predicted by the steady model problem [29]. Qualitatively, an excellent agreement is found since the same behaviors are theoretically and numerically observed for stability and convergence speed with respect to the relaxation parameter $\alpha_f$. Quantitatively, a good numerical approximation of the theoretical stability limits is found. Nonetheless, slight discrepancies have been observed. Three sources of discrepancy will briefly be discussed:

(1) the relevance of a 1D model,

(2) the coupling numerical strategy adopted in this paper,

(3) the transient effects.

The first question is whether a 1D model can be relevant to predict the main numerical characteristics of a 2D/3D CHT computation? This study has shown that the general trends of theoretical stability bounds (stability limit, optimal coefficient) are found. Furthermore, the key role played by the Biot number as suggested by the theory has been confirmed. Thus we may note that a simplified model provides interesting results because one may reasonably assume that the modes that may be unstable are those whose variation is in the direction normal to the coupled interface and consequently a 1D model problem may capture the main unstable modes. Obviously, the 2D/3D effects are neglected in both the fluid and solid domains but a model problem remains a precious guide in multiphysics computation.

It is also fundamental to remember the second point regarding the coupling strategy. As mentioned in this paper (and shown in Fig. 1), a 2-way thermal coupling takes place only at coupling time

steps. The quantities are updated at intermediate temporal points through an interpolation and no coupling is performed at these instants on the fluid side (1-way coupling). This can be referred to as a “partial” coupling. This partial coupling may have different stability characteristics compared to a fully coupling approach. This could explain the discrepancy –a more restrictive stability bound– observed by the Dirichlet–Robin procedure.

The third point about the transient effects is essential. A “steady” theory has been applied in this paper for weakly transient problems with promising results. But is this theory still valid if transient effects are present? To reply to this question, it is first necessary to look at the time scales, and specifically at the mesh Fourier numbers respectively equal to 5.71 and 1428 (Table 5). As the solid domain contains 10 mesh points in the y-direction, this means that each increment gives the time for the heat wave to penetrate the wall thickness. On the contrary, as far as transient effects are concerned, a small Fourier number must be adopted (typically 0.5) and a Biot number cannot be defined anymore. There is a recent theory based on a normal mode stability analysis that gives optimal coefficients in the case in which a transient medium is considered [15]. The same tendencies are recovered, but the stability bounds are slightly different, directly dependent on the mesh Fourier number. This transient model problem is likely to be more suitable and justifiable for steep transient problems. But otherwise, and for the type of problems considered in this paper, it would be costly and not useful to employ a very small Fourier number. In sum, a steady model problem is easier to implement and works very well as long as transient effects are weak.

### 5.2. Partitioned method vs monolithic method

Interface conditions influence on accuracy is numerically investigated. Temperature evolution over time at the fluid–solid interface is studied. The results of the partitioned method are compared with those of the monolithic method, which can be regarded as the reference in terms of accuracy.

In Fig. 19, a comparison of the temperature evolution over time at the interface point $x = 17.5$ cm is shown, between a monolithic

<figure><img src="image_20.png" /><figcaption>**Fig. 19.** Evolution of the temperature over time at point $x = 17.5$ cm of upper face (Table 3). Comparison of monolithic method with partitioned method (Robin condition imposed on the solid with relaxation parameters $\alpha_f = h$, $\alpha_f = 2h$ and $\alpha_f = 10^6$).</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>**Fig. 20.** Zoom of Fig. 19, between time 5400 s and 7800 s.</figcaption></figure>

method and a partitioned method. With this latter, Robin conditions are imposed on the solid side, with several relaxation parameters. Calculations are performed with physical parameters specified in Table 3. Coupling time steps are also indicated with black points.

The error of the partitioned method compared to the monolithic method remains globally low. It is especially the case in steady phases (between 3000 s and 5400 s, and between 7200 s and 10800 s) where both methods give the same temperature, independently of the relaxation parameter chosen in the partitioned method.

Fig. 20 is a zoom of Fig. 19 during a transient phase, between times 5400 s and 7800 s. With the partitioned method, temperature at coupling time steps is independent of the relaxation parameter, in accordance with the theoretical analysis (Eq. (19)).

However, between the coupling time steps, at the beginning and the end of the transient phases (between 5400 s and 6000 s, and between 6600 s and 7200 s), the temperature is dependent on the relaxation parameters. It is caused by temporal linear interpolation of interface conditions imposed on the solid side (see Section 2.1 Coupling algorithm). Thus, with a relaxation parameter of $10^6$ (almost equivalent to a temperature condition imposed on the solid side), the interface temperature has a linear evolution over time, between coupling time steps. More generally, the temperature tends towards a linear evolution when the relaxation parameter imposed on the solid side increases. For instance, the temperature obtained with the relaxation parameter $2h$ is closer to a linear evolution than the temperature obtained with the relaxation parameter $h$ (Fig. 20). In reality, because of thermal inertia of the solid, temperature has an exponential type evolution at beginnings and ends of transient phases, as shown by the reference monolithic method. Similarly, the heat flux tends towards a linear evolution when the relaxation parameter decreases, and this may cause the same problems. That is why an « intermediate » relaxation parameter is usually necessary.

Finally, it is observed that the temperature field given by the quasi-dynamic method is slightly overestimated in transient phases, with respect to the temperature predicted by the monolithic method. The temperature differences are around $\Delta T = 0.17$ K between times 6000 s and 6600 s (Fig. 20). It can also be seen as a temporal advance of the quasi-dynamic method compared to the monolithic method, around $\Delta\tau = 28$ s. Indeed, with the quasi-dynamic method, interface conditions are linearly interpolated without taking into account the thermal inertia of the solid, and thus are in advance in transient phases.

## 6. Conclusion and perspectives

Two interface fluid-structure transmission procedures for weakly transient heat transfer problems have been considered and examined in the paper: Dirichlet–Robin and Neumann–Robin conditions. It has been shown that they give rise to disparate numerical behaviors and that the Biot number is a key parameter for determining efficient interface conditions, in terms of stability and convergence. At low or moderate Biot numbers, it is better to impose the temperature on the fluid side, and a Robin condition on the solid side with a low relaxation parameter. On the contrary, at high Biot numbers, it is more relevant to impose the heat flux on the fluid side, and a Robin condition on the solid side with a high relaxation parameter. Moreover, whatever the Biot number, the heat transfer coefficient is always the optimal relaxation parameter in terms of convergence speed. Thus the validity of Verstraete theory [29], developed for steady CHT problems, can be extended to weakly transient CHT problems. It is noticed that relevant interface conditions can be obtained even if the convection coefficient is not precisely known, which is often the case in complex industrial CHT computations.

Moreover the quasi-dynamic partitioned method has been validated in terms of accuracy by being compared with monolithic reference method. The agreement is excellent in steady phases. In transient phases, a small advance of the quasi-dynamic method is observed compared to the monolithic method. Moreover, with a quasi-dynamic method, the temperature at coupling time steps is independent of the interface conditions. However, between the coupling time steps, the temperature depends on relaxation parameters because of the linear interpolation of the interface conditions imposed on the solid side, and thus a relaxation parameter value close to the heat transfer coefficient is advised.

Future work must now include strongly transient conjugate heat transfer effects. This means considering physical phenomenon of thermal inertia via the Fourier number. The work presented in this paper constitutes a sound basis for this process.

## Acknowledgements

The authors wish to thank ANSYS for their help. The authors are also grateful to the ANRT for the financial support.

## References

[1] T.-L. Perelman, On conjugated problems of heat transfer, Int. J. Heat Mass Transf. 3 (1961) 293–303.

[2] M. Heil, An efficient solver for the fully coupled solution of large-displacement fluid-structure interaction problems, Comput. Methods Appl. Mech. Eng. 193 (2004) 1–23.

[3] B. Hubner, E. Walhorn, D. Dinkler, A monolithic approach to fluid-structure interaction using space–time finite elements, Comput. Methods Appl. Mech. Eng. 193 (2004) 2087–2104.

[4] F. Rahman, J.-A. Visser, R.-M. Morris, Capturing sudden increase in heat transfer on the suction side of a turbine blade using a Navier–Stokes solver, J. Turbomach. 127 (3) (2005) 552–556.

[5] J. Luo, E.-H. Razinsky, Conjugate heat transfer analysis of a cooled turbine vane using the V2F turbulence model, J. Turbomach. 129 (4) (2007) 773–781.

[6] C. Felippa, K. Park, Staggered transient analysis procedures for coupled dynamic systems: formulation, Comput. Methods Appl. Mech. Eng. 24 (1980) 61–112.

[7] S. Piperno, C. Farhat, B. Larrouturou, Partitionned procedures for the transient solution of coupled aeroelastic problems, Comput. Methods Appl. Mech. Eng. 124 (1980) 79–112.

[8] K.-C. Park, C.-A. Felippa, Partitioned analysis of coupled systems, Comput. Methods Transient Anal. (1983) 157–219.

[9] J. Degroote, Partionned simulation of fluid-structure interaction, Arch. Comput. Methods Eng. 3 (2013) 185–239.

[10] M.-B. Giles, Stability analysis of numerical interface conditions in fluid-structure thermal analysis, Int. J. Numer. Meth. Fluids 25 (1997) 421–436.

[11] B. Roe, R. Jaiman, A. Haselbacher, P.-H. Geubelle, Combined interface boundary method for coupled thermal simulations, Int. J. Numer. Meth. Fluids (2008) 329–354.
[12] W. Henshaw, K. Chand, A composite grid solver for conjugate heat transfer in fluid-structure systems, J. Comput. Phys. 228 (2009) 3708-3741.

[13] V. Kazemi-Kamyab, A. van Zuijlen, H. Bijl, Accuracy and stability analysis of a second-order time-accurate loosely coupled partitioned algorithm for transient conjugate heat transfer problems, Int. J. Numer. Meth. Fluids (2013).

[14] O. Joshi, P. Leyland, Stability analysis of a partitioned fluid-structure thermal coupling algorithm, J. Thermophys. Heat Transfer 28 (1) (2014).

[15] M.-P. Errera, S. Chemin, Optimal solutions of numerical interface conditions in fluid-structure thermal analysis, J. Comput. Phys. 245 (2013) 431-455.

[16] J. Lindstrom, J. Nordstrom, A stable and high-order conjugate heat transfer problem, J. Comput. Phys. 229 (2010) 5440-5456.

[17] F.-X. Roux, J.-D. Garaud, Domain decomposition methods methodology with Robin interface matching conditions for solving strongly coupled fluid-structure problems, Int. J. Multiscale Comput. Eng. 7 (2009) 9-16.

[18] M. Errera, G. Turpin, Temporal multiscale strategies for conjugate heat transfer problems, J. Coupled Syst. Multiscale Dyn. 1 (2013) 89-98.

[19] Z. Zhai, Q. Chen, Impact of determination of convective heat transfer coefficient on the coupled energy and CFD simulation for buildings, Proc. Build. Simul. Conf. 3 (2003) 1467-1474.

[20] I. Beausoleil-Morrison, The adaptive coupling of computational fluid dynamics with whole-building thermal simulation, in: 7th International IBPSA Conference, Rio de Janeiro, Brazil, 2001.

[21] Q. Chen, Z. Zhai, L. Wang, Computer modeling of multiscale fluid flow and heat and mass transfer in engineered spaces, Chem. Eng. Sci. 62 (2007) 3580-3588.

[22] Z. Sun, J. Chew, N. Hills, K. Volkov, C. Barnes, Efficient finite element analysis/computational fluid dynamics thermal coupling for engineering applications, J. Turbomach. 132 (2010).

[23] Z. Sun, J. Chew, N. Hills, L. Lewis, C. Mabilat, Coupled aerothermomechanical simulation for a turbine disk through a full transient cycle, J. Turbomach. 134 (2012).

[24] V. Ganine, U. Javiya, N. Hills, J. Chew, Coupled fluid-structure transient thermal analysis of a gas turbine internal air system with multiple cavities, J. Eng. Gas Turbines Power 134 (2012).

[25] D. Culler, J. McNamara, Studies on fluid-thermal-structural coupling for aerothermoelasticity in hypersonic flow, AIAA J. 48 (2010) 1721-1738.

[26] E. Thornton, P. Dechamphait, Coupled flow, thermal, and structural analysis of aerodynamically heated panels, J. Thermophys. Heat Transfer 11 (1997) 173-181.

[27] M.-P. Errera, B. Baqué, A quasi-dynamic procedure for coupled thermal simulations, Int. J. Numer. Meth. Fluids 72 (2013) 1183-1206.

[28] C. Shen, F.-X. Sun, X.-L. Xia, Analysis on transient conjugate heat transfer in gap-cavity-gap structure heated by high speed airflow, Int. J. Heat Mass Transf. 67 (2013) 1030-1038.

[29] T. Verstraete, Multidisciplinary turbomachinery component optimization considering performance, stress, and internal heat transfer (PhD thesis), Universiteit Gent, 2008.

[30] S. Godunov, V. Ryabenkii, The theory of difference schemes – an introduction, 1964.