Nicolas Chauvet
email: chauvet.nicolas@yahoo.fr

A steady-transient thermal coupling method for fluid-solid interfaces

Thermal analysis is of major concern in many fields where energy is present in various forms, like life science, geoscience, building and transports. Energy production and transfer is often associated with fluids and its mechanical exploitation is most of the time performed by mechanisms made of solid parts. This is especially the case in gas turbine engines, where the steady thermal states and the transient thermal evolutions in the solid parts play a crucial role not only on their lifetime, but also respectively on the maximum supplied power and on the acceleration capabilities. In order to predict the thermal resistance and deformation of the engine parts, engineers now make extensive use of numerical simulations. Under certain conditions, it is possible to separate the thermal computation in the solids from the thermo-dynamics computation of the fluid, which allows to reduce the overall computational cost. Though, if we take a look at the practices used to deal with the thermal exchanges at the interface between solids and fluids, they are still not satisfactory. Some "pragmatic methods", in the sense that they allow to get numerical results, ignore basic physical principles. "Mathematical methods", based on a mathematical analysis of equations, introduce numerical parameters which are not calculated by physical considerations and are set to arbitrary values by engineers.

The present analysis starts from simple mathematical and physical observations, then provides a real pragmatic method in the sense that the engineer has not to speculate on the proper numerical parameters to be used for its application. The usual condition employed in thermal solid computations to represent the influence of a fluid flowing on a portion of its boundary is the famous Newton's law of cooling. It is indeed a good starting point, but owing to its limitations, a more general model is needed. In the present coupling method, it is proposed to use a local linearization of the thermal fluid behavior on the fluid-solid interface and to perform an automatic identification of its coefficients. With the help of a slow solid -fast fluid assumption, the coupling method is finally extended from steady thermal analysis to transient thermal analysis.

Basic physical/modeling principles

In this section, two fundamental results of thermal coupling between a solid and a fluid, respectively pertaining to steady and transient issues, will be highlighted by considering a very simple problem. We consider a solid plate of constant thickness d. One face is maintained at a constant uniform temperature T bs while the opposite face is in contact with a flowing fluid, the bulk temperature T bf of which is controlled from the inlet. Conditions are such that gravity and radiation can be ignored. The distance d between the faces is short enough to consider locally the problem of temperature distribution as one-dimensional in the direction -→ n normal to the plate (variations in this direction are predominant over those in the streamwise direction; the heat flux and the temperature gradient are considered as collinear to -→ n ).

Steady solution of the basic thermal problem

In the fluid, the bulk fluid temperature T bf is imposed. At the interface with the solid, heat flux and temperature will be respectively referred to as q f and T if . Under particular assumptions, when the fluid reaches a steady state, its governing equations (mass conservation, momentum and energy, completed with thermodynamic state laws and constitutive laws like Newton's law on viscous stress and Fourier's law on conductive heat flux) lead to a proportional relation between the heat flux and the temperature difference:

q f = -h f .(T if -T bf ) (1)
Here q f is the heat flux going out of the fluid through the interface (i.e. entering the solid) and h f a positive constant called "heat transfer coefficient". This relation is known as "Newton's law of cooling" after Newton's early work on the temperature measurement of a metal cooled by a uniform wind [START_REF] Newton | Scala graduum caloris, calorum descriptiones & signa (scale of the degrees of heat)[END_REF]. But it can also be applied to the warming of a solid by a forced hot flow. It is worth noting that the term "law" is somewhat confusing because this relation is not universal. It is a consequence of general fundamental laws applied to particular fluids, but even more importantly, in particular circumstances. For instance, it applies when constitutive laws are linear, radiation is ignored and convection is only driven by an external pressure gradient, i.e. when natural convection is negligible.

The solid is an infinite plate on which one face is maintained at a given boundary temperature T bs and the opposite face is in contact with the fluid. On this interface, heat flux and temperature are respectively referred to as q s and T is . Fourier's law expresses the value of the heat flux in the solid as a proportional quantity of the temperature gradient. Since when the thermal field has reached a steady state, no energy is stored inside the volume, the integration of Fourier's law suffices to relate the heat flux going out of the solid through the interface to the difference between the bounding temperatures:

q s = - λ d (T is -T bs ) ( 2
)
where λ is the thermal conductivity and d the distance separating the two parallel faces of the solid plate.

The basic coupled problem consists in determining the temperature at the fluid-solid interface in the steady state, denoted T ∞ i . At any time, by continuity of the temperature field, the fluid temperature at the interface T if and the solid temperature at the interface T is are equal, and owing to the energy conservation, heat fluxes compensate each other:

T if = T is (3) q s = -q f ( 4
)
This also holds in the steady state where quantities are denoted by the upper-script ∞.

Matching conditions (3) and ( 4) combined with the steady solutions of each material domain ( 1) and ( 2) lead to the overall steady solution

T ∞ i = T ∞ if = T ∞ is : T ∞ i = λ d T bs + h f T bf λ d + h f = T bs + B i T bf 1 + B i (5)
where

B i = h f λ/d
is the Biot number. It shows the relative strength of convection with respect to conduction:

T ∞ i → T bs when B i → 0 and T ∞ i → T bf when B i → ∞.
In the next two subsections, we will consider that the steady state in the fluid and the steady state in the solid are calculated separately and look at the behavior of two coupling approaches. Both methods are iterative. We start with an initial estimation of the solid interface temperature T 0 is and, from an estimation T n-1 is at a iteration n -1, the coupling method provide a new estimation T n is at the next iteration n.

Thermal coupling by imposed T is and q f

In this first method of coupling, at a each iteration n, the solid temperature T n-1 is known from the previous iteration (n-1) is imposed to the fluid, then the thermal equilibrium in the fluid gives the heat flux at the interface q n f , this heat flux is imposed to the solid, and finally the thermal equilibrium in the solid gives a new estimate of the solid temperature T n is . This iterative process is initialized with a guess temperature T 0 s . In the simple problem described above, analytical steady solutions are known from (1) and ( 2):

q n f = -h f (T n-1 is -T bf ) (6) T n is = T bs + d λ q n f (7)
This allows to express the evolution of the difference of the temperature at iteration n with respect to the steady state value T ∞ i (5) as a geometric series:

(T n is -T ∞ i ) = - h f λ/d (T n-1 is -T ∞ i ) (8)
The recurrence relation (8) shows that this first method can converge toward the solution T ∞ i only if conduction dominates convection in the sense that B i < 1. Moreover, this convergence is oscillating, the temperature difference changing sign at each iteration. For conditions where B i > 1, the difference with the steady state temperature will grow indefinitely, this method will diverge.

To solve problems with B i > 1, one can apply a similar method with inverted choices of the quantities imposed as boundary conditions: at a each iteration n, the fluid temperature T n-1 if known from the previous iteration (n -1) is imposed to the solid, then the thermal equilibrium in the solid gives the heat flux at the interface q n s , this heat flux is imposed to the fluid, and finally the thermal equilibrium in the fluid gives a new estimate of the fluid temperature T n if . The iterative process is initialized with a guess temperature T 0 f . In the simple problem described above, analytical steady solutions are known from ( 1) and ( 2):

q n s = - λ d (T n-1 if -T bs ) (9)
T n if = T bf + 1 h f q n s ( 10
)
This allows to express the evolution of the difference of the temperature at iteration n with respect to the steady state value T ∞ i (5) as a geometric series:

(T n if -T ∞ i ) = - λ/d h f (T n-1 if -T ∞ i ) (11)
The recurrence relation (8) shows that this method can converge toward the solution T ∞ i only if convection dominates conduction in the sense that B i > 1. As previously, this convergence is oscillating, the temperature difference changing sign at each iteration. For conditions where B i < 1, the difference with the steady state temperature will grow indefinitely, this method will diverge.

Thermal coupling by imposed T is and (h, T ref )

In many cases, Newton's law of cooling ( 1) is a fair approximation of the thermal behavior of the fluid at the interface. More importantly, as shown later in this paper, its form can cover any flow in a limited range of wall temperature. Since T is = T if and q s = -q f , it is sensible to impose directly this relation as a boundary condition for the solid problem:

q s = h.(T is -T ref ) (12)
with h = h f and T ref = T bf . On the mathematical viewpoint, the boundary condition (12) is relaxed compared to an imposed temperature T is or an imposed heat flux q s , and the resolution of the solid problem provides both T is and q s . Now, let us consider the simple problem introduced previously and assume that we do not know that the fluid behaves like [START_REF] Newton | Scala graduum caloris, calorum descriptiones & signa (scale of the degrees of heat)[END_REF]. That means that we can use the relation (1) to determine q f from T if but we are not allowed to replace directly h by h f and T ref by T bf in (12), even if we already know that it is the only logical thing to do.

In this subsection, we will adopt the same coupling method as Verstraete et al. in [START_REF] Verstraete | Numerical study of the heat transfer in micro gasturbine[END_REF]: a relaxed boundary condition like (12) is imposed to the solid at the interface with the fluid, the value of the heat transfer coefficient h (mathematical, not physical) is chosen and remains constant during all the computation process, and the reference temperature T ref is determined by the following coupling procedure. At a each iteration n:

1) the solid temperature T n-1 is known from the previous iteration is imposed to the fluid; 2) then the thermal equilibrium in the fluid gives the heat flux at the interface q n f ; 3) this resulting heat flux q n f and the prescribed temperature T n-1 is are used to assess T ref that satisfies (12), h being arbitrarily chosen, with q s = q n f and T is = T n-1 is ; 4) h and T ref being known, the boundary condition ( 12) is imposed to the solid; 5) finally the thermal equilibrium in the solid gives a new estimate of the temperature T n is .

This iterative procedure is initialized with a guess temperature T 0 s .

In this simple problem, the analytical steady solution in the fluid is known from [START_REF] Newton | Scala graduum caloris, calorum descriptiones & signa (scale of the degrees of heat)[END_REF] and that in the solid is given by (5) where h f and T bf must be replaced by h and T ref :

q n f = -h f .(T n-1 is -T bf ) ( 13
)
T n is = λ d T bs + hT n ref λ d + h (14)
and T n ref is evaluated by equalizing (13) and (12):

T n ref = T n-1 is + h f h (T bf -T n-1 is ) (15)
Mixing (13), ( 14), (15) and (5) allows to express the evolution of the difference of the temperature at iteration n with respect to the steady state value T ∞ i (5) as a geometric series:

T n is -T ∞ i = h -h f h + λ/d (T n-1 is -T ∞ i ) (16)
The temperature difference T n is -T ∞ i is the error to the steady state solution at iteration n. The ratio of the temperature error from one iteration to the next one, or gain G, reads:

G = T n is -T ∞ i T n-1 is -T ∞ i = h h f -1 h h f + λ h f d = h h f -1 h h f + B -1 i with B i = h f λ/d ( 17
)
This gain is fully determined by the ratio of the numerical heat transfer coefficient h to the physical heat transfer coefficient h f . The coupling method converges toward the solution if and only if

|G| < 1, i.e. if h/h f > (1 -B -1 i )/2.
Then, if B i < 1, the coupling method converges whatever the finite value of h, since h > 0. The convergence is oscillating (the temperature error changing sign at each coupling iteration) if h/h f ∈](1 -B -1 i )/2, 1[ and it is monotonic if h/h f > 1. Inside the stable domain, the convergence is slower toward the boundaries h/h f → (1 -B -1 i )/2 + and h/h f → +∞. This is summed up in figure 1. The most important result concerning this method is that the convergence is the fastest and immediate (T 1 is = T ∞ i ) for the only sensible value h = h f . The second result is that, under this condition h = h f , the convergence of the iterative method in a single iteration, is independent of the Biot number and of the initial temperature. The conclusion of this simple analysis, even though limited to a simple problem, is that the coupling method must first of all determine the physical heat transfer coefficient of the fluid h f . In this simple model [START_REF] Newton | Scala graduum caloris, calorum descriptiones & signa (scale of the degrees of heat)[END_REF], two different wall temperatures (or two different heat fluxes) can be imposed to deduce the coefficients h = h f and T ref = T bf . On the whole, two steady fluid computations and one steady solid computation with the relaxed boundary condition (12) are enough to get the solution of the coupled thermal problem. This method can be easily extended to 3D geometries. In practice, fluid flows are often computed with adiabatic boundary conditions on the walls. This gives a temperature distribution T 1 (P ) on the walls (P point on the wall). A second fluid computation can be performed with either an imposed heat flux δq or an imposed temperature T 2 (P ) = T 1 (P )+δT on the walls. Due to the three dimensions, the heat flux can have a component tangent to the wall so that issues can arise in the computation of h(P ) = h f (P ) > 0 or T ref (P ) = T bf (P ) > 0. When it occurs, smaller values of increment δq or δT can be tried. Once the coefficients of (1) are known, it only remains to perform the thermal computation in the solid.

Slow solid -fast fluid assumption

Before addressing the question of transient thermal evolution in the solid, this section introduces the crucial assumption upon which the simplifications/accelerating procedure, exposed in next sections, can hold. Before all, if thermal analysis is used to assess the displacement of engine parts or their lifetime under low cyclic fatigue, conduction in the solid parts is the major topic of interest. Though, from the physical and mathematical point of view, it is the easiest to solve. The hardest issues to deal with are the thermal boundary conditions at the interfaces with the surrounding fluids. Radiation that will not be discussed at all in this paper is already complex when it involves only solid surfaces, but it is significantly made harder when the gases between facing surfaces are not perfectly transparent. Radiative fluxes depends on the relative positions and orientations of the radiative elements and of the receivers in the spatial domain but also on the wave frequencies and their 3D directions. That is why very simplified models are always used to deal with radiation. But hopefully, radiation is an instantaneous phenomenon, which means that it does not add more complexity to a transient problem. Also of great concern, contacts between solid parts will not be discussed herein. In this paper, only the heat exchanged by conduction through the contact between solids and fluids is treated.

When dealing with the thermal evolution of a system involving a solid and a fluid in contact, it is inevitable to compare the heat transfer by conduction to the heat transfer by convection. For thermal engineers focused on the temperature distribution in the solid, "heat transfer by conduction" means "heat transfer inside the solid", while "heat transfer by convection" means "heat transfer between the solid and a neighboring fluid", it is just a boundary condition for the solid. For a physicist focused on the fundamental mechanisms of heat transfer in the different material domains, heat transfers are limited to conduction inside the solids while they are a combination of conduction and convection inside the fluids. At the level of a material surface element, the heat flux by conduction is formalized by Fourier's law -λ∂T /∂n and the heat flux by convection is the internal energy e advected by the flow velocity -ρeV n1 . Approaching the wall, the fluid velocity vanishes, then the portion of heat convected into the volume of fluid in the vicinity of the solid vanishes too. By nature, the boundary layer is the place where velocity gradients and temperature gradient are the strongest in the fluid. Because of viscosity, velocity gradients induce dissipation, mainly µ(T )(∂u/∂y)2 (u longitudinal velocity, y coordinate normal to the wall). This is "thanks to" this dissipation (and to conduction) that kinetic energy is transformed into internal energy and is not negligible when the external velocity is high. When this effect is isolated from any external heat source, i.e. when the wall is insulated, this viscid effect rises the external static temperature T ∞ toward the warmer "recovery temperature" T r on the wall. In the particular case of fluids where P r = µc p /λ ≈ 1, all the mechanisms compensate in the boundary layers, which gives the illusion that we recover the external total temperature since

T if = T r ≈ T ∞ + V 2
∞ /2c p like if it was an inviscid process that brought velocity down to zero. In low Mach flows,

M 2 ∞ = V 2 ∞ /γRT ∞ 1, this effect is negligible so that T r ≈ T ∞ .
In the assumption of the Newton's law of cooling (1), the bulk fluid temperature T bf is equal to the recovery temperature T r . Overall, the steady heat transfer at the interface between a fluid and a solid then depends on the fields of density ρ and velocity V . In flows where density can be considered as constant, the velocity field is decoupled from the temperature field, it is only governed by pressure gradients through the momentum equation. In other words, the dynamic boundary layer is independent of the thermal boundary layer, whereas the thermal boundary layer always depends on the dynamic boundary layer because of the convection of internal energy. Strictly, in all other situations, (ρ, V , T ) are indissociable and coupled through the system of equations formed by the conservation of mass, the momentum balance, the internal energy balance, fluid state laws and constitutive laws. Moreover, gravity is always there. Temperature gradients in the fluid cause density gradients, hence buoyancy forces. The velocity field is then governed, mainly, through the momentum equation by both pressure gradients and temperature gradients. When the pressure gradients dominate the buoyancy forces, the flow regime is said in "forced convection". When the buoyancy forces caused by temperature gradients dominate the pressure gradients, the flow regime is said in "natural convection". Between these two regimes, we talk about "mixed convection". When the temporal evolution of the temperature field is analyzed, things are made more complex not only because heat is stored (term ∂e(T )/∂t in the energy balance equation) in any material domain, solid or fluid, but also because the heterogeneous dilatation in the fluid can disturb its organization, and takes its own characteristic time to follow the change in the boundary conditions. This rapid overview allows to realize that behind the expression "heat transfer by convection" seen as a simple boundary condition by the solid, is hidden a very complex physical machinery involving not only the conditions (T if , q f ) on the fluid-solid interface, but also the flow inlet and outlet conditions, the material properties of the fluid (viscosity, conductivity, heat capacity) and the external acceleration (gravity and non-inertial acceleration, e.g. on rotating walls). The purpose of this paper is not too go deeper into the mechanisms of flows but rather to consider their macroscopic thermal effect on the solid boundary.

Even if conduction in the solid is the easiest part of the problem to solve, it must be treated with care, especially in transient evolutions. When the temperature is changed at the surface of a solid body, the whole temperature field inside the body is affected instantaneously by conduction but with a magnitude that decays away from this surface, going deeper into the solid (see, for instance, Stokes first and second problems). Thus, the reasonable question to ask is, given a temperature increment, portion of the temperature variation at the surface, how much time τ does it take to feel this temperature increment at a given distance d from the surface ? A fundamental result of 1D conduction is that, as for any diffusion process, this time is proportional to the square of the distance: τ ∝ d 2 . This is just to remind that, even in a thick solid, thermal conduction can be very fast from an unsteady heat source to a probe, when their distance is short enough.

Finally, when considering the temporal evolution of the temperature distribution in a solid, it is important to know whether the boundary conditions depend on time or not. A boundary condition can explicitly depend on time t, e.g. T is (t) or q s (t). In this case, it is a forced unsteady boundary condition. A boundary condition can also depend on time through a temporal derivative, if for instance it is connected to a heat capacity, e.g. q s = CdT is /dt. And finally, a boundary condition can depend on time through its own state variables; it is the case with the relaxed condition (12), where q s (t) is related to T is (t). But in this third situation, the dependence on time is brought only by the instan-taneous solid state. Such a kind of boundary condition can be more generally written: f (q s , T is ) = 0, where f is a function which is not necessarily linear like (12). The physical interpretation of this form of boundary condition is that the external environment instantaneously adapt to the variations of the conditions on the interface with the solid, it does not add the influence of its own inertia. Said differently, when the boundary conditions (q s , T is ) change, the characteristic time for the flow to adapt is significantly lesser than the characteristic time needed to feel the temperature variation at a given distance d from the interface in the solid. And as underlined just above, this assumption cannot hold when d → 0. There is always a minimum non-zero distance under which this assumption cannot be satisfied. When the points of interest are located at a sufficient distance d from the fluid-solid interface, the boundary conditions on the interface can be approached by a relation of the form f (q s , T is ) = 0, i.e. neglecting the time for the flow to adapt to variation in (q s , T is ). In this framework, the flow computation can be separated from the transient solid computation and since the transient time to reach steady states in the flow are considered as negligible, steady flow computations can be performed.

Influence of the heat transfer coefficient in transient evolution

Still today, transient thermal computations are performed in industry, using the boundary condition ( 12) with arbitrary values of the heat transfer coefficient h as if it had no influence on the temperature evolution. That is why it is still necessary to show its influence. This is the purpose of this short subsection.

Typically, engineers have to produce numerical estimations quickly, so they want to perform the least computations as possible, and sometimes even less than what is necessary. If the fluid flow on a wall can be simply described by the Newton's law (1), only two computations are necessary but some thermal engineers want to perform only one computation (or just use the results of adiabatic computations performed by aerodynamicists). A practice, unfortunately quite common, is to probe a temperature value at an arbitrary distance from the wall (for instance the value in the first computed point above the wall, giving rise to a dependence on the mesh...). Once a temperature is chosen for the reference temperature, the results (T if and q f ) of the sole computation can be used in (12) to calculate a numerical value for h, which has no physical justification.

On the interface, (12) will hold for any time iteration, till the convergence toward the steady state (T ∞ is and q ∞ s ), then:

q s = q ∞ s -h.(T is -T ∞ is ) (18)
If we consider the simple 1D problem introduced previously, at a given time t, the variation in the internal energy of the solid will result from the sum of the fluxes on the boundaries. More particularly, on the fluid-solid interface, the heat flux will be given by ( 18). If we consider that the steady state (q ∞ s , T ∞ is ) is unchanged by the choice of the couple (h, T ref ), which is satisfied since we have considered that q ∞ s = h(T ∞ is -T ref ), the choice of h will change the heat flux q s during the transient phase and then the time required to reach the steady state (or more rigorously, to approach the steady state up to a chosen accuracy).

If the thermal behavior of the fluid is correctly reproduced by Newton's law of cooling (1), a value of the numerical parameter h in (12) different from the heat transfer coefficient h f in (1) will change the transient evolution of the temperature in the solid.

For instance, if the flow cools an initially warmer solid, T is > T ∞ is during the transient phase, and a higher value of h will increase the value of the heat transfered from the solid to the fluid, then reduce the time to approach the steady state. And vice versa, a lower value of h will increase the time to approach the steady state.

To sum up, when the physical conditions of the flow and the solid and the expected accuracy of the temperature field allow to model the thermal behavior of the fluid with Newton's law of cooling (1), the computation of the fluid and the computation of the temperature in the solid can be decoupled. In this case, the only logical choice is to directly impose Newton's law of cooling (1) as the boundary condition on the fluid interface of the solid, or said differently to use the relaxed boundary condition (12) with h = h f and T ref = T bf . This preliminary study showed that any other choice will rise numerical issues in the coupling process for steady computations (from slower convergence to divergence) and will lead to unphysical unsteady behavior in transient simulations. The advantage of decoupling the computations of fluid and solid when the slow solid -fast fluid assumption (1.4) is verified, is that long physical duration can be simulated, because the expensive unsteady fluid computations (for which a very short timestep is required compared to that in the solid) can be replaced by far cheaper steady fluid computations, which can benefit from acceleration techniques. The next two sections show how to generalize this method for the computation of steady or transient temperature fields.

Steady-state thermal coupling

In section 1.3, a coupling procedure was given to match a fluid computation with a thermal solid computation. It was relying upon the assumption that the thermal fluid behavior on the solid interface was correctly modeled by Newton's law of cooling [START_REF] Newton | Scala graduum caloris, calorum descriptiones & signa (scale of the degrees of heat)[END_REF]. In this relation, the heat transfer coefficient is independent of the wall temperature. In forced convection, its value is increasing with the flow velocity, because the flow velocity enhances the heat transfer by convection. In natural convection, the flow moves under the effect of buoyancy, itself provoked by the temperature gradients. In this case when the wall temperature changes, the heat transfer coefficient changes. That is why Newton's law of cooling with a constant heat transfer coefficient cannot correctly reproduce the heat exchange during the cooling or warming by natural convection. When the fluid flows in a cavity with a small mass flow rate, the heat exchange with the solid can change significantly the average temperature in the cavity, which must be reproduced by a change in the bulk fluid velocity T bf in [START_REF] Newton | Scala graduum caloris, calorum descriptiones & signa (scale of the degrees of heat)[END_REF]. Those two examples show how boundary conditions (T if and q f ) can affect the coefficients h f and T bf in [START_REF] Newton | Scala graduum caloris, calorum descriptiones & signa (scale of the degrees of heat)[END_REF]. In other words, the function f (q f , T if ) = 0 is not necessarily linear as claimed by Newton's law of cooling.

General steady method

In order to generalize the method presented in section 1.3, we will consider that the thermal behavior of the fluid at the interface can be written as:

q f = g(T if ) (19)
The function g is here not necessarily linear, but in the vicinity of a value of interface temperature T a if , the heat flux q b f at a temperature T b if close to T a if can be deduced from the Taylor series:

q b f = q a f + dq f dT if (T a if ).(T b if -T a if ) + o T b if -T a if ( 20
)
The iterative process described in section 1.3 can then be adapted to deal with the nonlinearity of g by assuming that its local linearization at each iteration is a good approximation. In the steady solid computation, this is represented by the linearized boundary condition :

q s = q n ref + h n .(T is -T n ref ) (21)
where T n ref , q n ref and h n are updated at each coupling iteration n such that q s = -q b f and

T is = T b if , and that T ref = T a if is close enough to T is = T b if .
During the iterative process, if at an iteration n, the temperature and the heat flux on the interface of the fluid are known, respectively referred to as T n if and q n f , they can be used as T a if and q a f in (20):

T a if = T n if (22) q a f = q n f (23)
With the help of an additional steady fluid computation with an imposed temperature distribution T n,bis if close to T n if , providing the heat flux q n,bis f , the first derivative

dq f
dT (T a if ) can be estimated by finite difference:

dq f dT if (T a if ) ≈ q n,bis f -q n f T n,bis if -T n if (24)
It is important to note that all the quantities in (24) are distributions on the interface I and that (24) can be defined only if:

∀P ∈ I : T n,bis if (P ) = T n if (P ) (25)
A simple choice to respect this condition is to set T n,bis if = T n if + δT with δT > 0. However, the additional fluid computation can be avoided, if this condition is already satisfied by the temperature field at the previous iteration n -1. In that case, its results can be used in (24), i.e. T n,bis if = T n-1 if and q n,bis f = q n-1 f . After choosing the accuracy δT > 0 expected on the temperature field (it is not only the maximum temperature difference between two successive iterations allowed to stop the loop, see later equations ( 27) and ( 28)), a tolerance δT 2 ≥ δT for the calculation of the heat transfer coefficient h n based on (24) and a small real parameter > 0 preventing zero division ( < δT ), the iterative process can then be described as follows.

|T n-

= T n if +δT .
This iterative process is not intended to converge in any situations because of the nonlinearity of the function q f = g(T if ). That is why the best pragmatic solution to test the convergence is to apply the method on each case of practical interest (as a reminder, Navier-Stokes solvers are tested this way everyday and everywhere in the world !) and change the initial conditions when it diverges. If q f = g(T if ) is linear, the procedure will converge in one step whatever the initial temperature distribution T 0 is . If q f = g(T if ) is not linear but continuous, there is always a neighborhood of the solution where T 0 is can be set such that the process converges toward the solution. Indeed, the analysis of the simple problem studied in section 1.3 suggests that in a close neighborhood of the solution T ∞ is , a slight deviation in the assessment of the heat transfer coefficient does not prevent convergence. The closer T 0 is from T ∞ is , the faster the convergence is. And the more non-linear the function q f = g(T if ), the more sensitive the process to an accurate numerical assessment of the first derivative ∂q f /∂T if , tuned by , δT and δT 2 .

When the iterative process is stopped at step 5 (assuming that it happens), as the boundary condition ( 21) is satisfied on the interface I, that its coefficients have been set by (26) at step 3 and that T n if = T n-1 is after step 1:

∀P ∈ I : |T n is (P ) -T n if (P ))| < δT (27)
and ∀P ∈ I :

|q n s (P ) -(-q n f (P ))| < h n .δT (28)
The last computed temperature fields in the fluid T n f and in the solid T n s satisfy the steady state equations in their respective domains. On the interface, the temperatures are continuous within the accuracy δT (27) and the heat fluxes respect the conservation of energy within the accuracy h n .δT (28). They all match when δT → 0.

It is worth noting that in the simple case where no heat source is acting in the volume and that conductivity is uniform, the extremal values of temperature in the solid are reached on its boundaries, so that the convergence check, step 5, can be performed only on the temperatures of the interface.

The procedure described above is general and can be simplified when the fluid has a linear thermal behavior on the interface. In this case, there is no need for additional fluid computations to update the coefficients T 2 ref q 2 ref and h 2 , since they would produce the same linearized boundary condition (21), and the second steady solid computation would then give the same result as the first one. This is exactly the situation described in section 1.3: the coupling procedure converges in only one iteration. Therefore, it can be stopped after the first coupling iteration.

Similarly, when the fluid has locally (in the vicinity of the current interface temperature T n if ) a linear thermal behavior, it is not necessary to update the coefficients (26).

In order to save the unnecessary second fluid computation made at step 2, the local linearity can be checked. This means checking if the new fluid data (T n if , q n if ), given by the fluid computation at step 1, is compliant with the last linearization of the fluid thermal behavior (21):

-q n f = q n-1 ref + h n-1 .(T n if -T n-1 ref ) (29)
However, because of numerical inaccuracy and to remain in line with the tolerance specified on the temperature field, this check must be performed within a tolerance, also on the temperature on the interface and with the same value δT :

∀P ∈ I : |q n f (P ) -q n-1 f (P ) + h n-1 (P ).(T n if (P ) -T n-1 if (P ))| ≤ h n-1 (P ).δT (30)
This linearity check is not mandatory. Its aim is only to reduce the overall computational time when the thermal fluid behavior approaches a linear behavior, avoiding potentially one fluid computation (step 2) and one solid computation (step 4) that would be performed by the baseline procedure in 5 steps before convergence is detected at step 5.

Restriction to usual solid solvers

In usual finite-element solvers devised for thermal analysis, the linearized boundary condition ( 21) is not available. Only classical conditions like an imposed temperature distribution, an imposed heat flux distribution or a relaxed condition (12), where h and T ref are distributions, can be applied on the solid boundaries. In this section, it is shown how the general method described in 2.1 can be adapted to be used with usual solid solvers.

A relaxed boundary condition like (12) can still be used in the steady solid computation:

q s = h n (T is -T n ref ) (31)
if its coefficients h n and T n ref are updated at each iteration n so that they reproduce the linearized approximation of the fluid behavior (20). This latter can be reformulated as:

q b f = dq f dT if (T a if ).   T b if -T a if + q a f dq f dT if (T a if )   ( 32
)
where q s = -q b f and T is = T b if , and under the condition that

dq f dT if (T a if ) = 0.
The identification of the remaining terms in (31) and (32) gives:

h n = - dq f dT if (T a if ) (33)
T n ref = T a if - q a f dq f dT if (T a if ) (34)
During the iterative process, h n and T n ref can then be updated in (31) in the same way as T n ref , q n ref and h n are updated in (21) of the general method. The procedure described in section 2.1 and all the comments still hold. It is just necessary to modify steps 3 and 4.

Then, with the same prerequisites, the coupling iteration n is computed in five steps:

1) a steady fluid computation is performed with the imposed temperature T n if = T n-1 is on the interface and gives a new estimation of the heat flux on the interface q n f .

2) if ∀P :

< |T n-1 if (P ) -T n if (P )| ≤ δT 2 , then T n,bis if = T n-1 if and q n,bis f = q n-1 f
, otherwise an additional steady fluid computation with the imposed temperature T n,bis if = T n if + δT on the interface is performed, providing the heat flux q n,bis f .

3) the coefficients of the relaxed boundary condition (31) are updated:

h n = - q n f -q n,bis f T n if -T n,bis if then T n ref = T n if + q n f h n (35)
4) a steady solid computation is performed with the relaxed boundary condition (31) and coefficients (35), it gives a new estimation of the temperature in the solid T n s ;

5) the new temperature distribution in the solid T n s is compared to the previous estimate T n-1 s : if their difference is less than the expected accuracy |T n s (P ) -T n-1 s (P )| < δT at any point P of the discretized solid S, the iterative process is stopped, otherwise it is continued.

This initialization is the same as that described in section 2.1.

It is important to keep in mind that T n ref , q n ref and h n in (21) are spatial distributions on the fluid-solid interface. The form of (21) could make think that the fluid behaves locally as a one-dimensional system in the direction normal to the wall, i.e. that (21) applied at each point of the surface can be independent of what happens on the sides. But this is not reality because the heat flux can have components tangent to the wall. More generally, it is not impossible that -q if = q n ref while T if = T n ref , i.e. h n = 0 in (21). This case does not raise any issue with the general method described in 2.1, because (21) allows to impose locally a heat flux (Neumann boundary condition) on the solid boundary. By contrast, this situation cannot be properly treated by the method adapted to the classical relaxed boundary condition (31), because that latter can only impose a heat flux when it is zero. Therefore, the restricted coupling method described in 2.2 cannot apply when there is a point on the interface where dq f dT if = 0 during the coupling iterations.

As a final remark before extending the method to transient thermal computations, it is worth reminding that the iterations in each separate computation, resp. solid and fluid, must be continued until an accuracy coherent with the accuracy δT chosen for the coupling procedure is reached. As for the fluid computations, at each coupling iteration, the initial field must be set equal to the last converged flow field in order to limit the number of iterations needed to reach convergence. the solid at the locations of interest. But the temperature evolution inside the solid remains dependent on the evolution in time of the temperature and heat flux at the interface with the fluid. That is why an update of the coefficients of the boundary condition (39) must be performed at each time step of the solid transient simulation, thanks to new steady fluid computations. It is worth noting that the time step of the transient solid computation is already significantly larger than the time step required to perform an unsteady fluid computation. Therefore, the approximate coupling method between steady fluid computations and unsteady solid computations already reduces the overall computational time, compared to a monolithic approach (fluid and solid both solved in the same solver with the same marching time step imposed by the fluid).

However, the transient thermal evolutions can last several minutes and then require many solid time steps to be simulated. That is why it can be interesting to reduce the number of fluid computations when possible, that is to say, compute several solid time steps before updating the boundary condition (39). This does not mean that the unsteady boundary condition (39) is reduced to a steady boundary condition like (21). This means that the functions q ref (t), T ref (t) and h(t) in (39) are not updated at each time step of the transient solid solver. They are updated at meeting dates, between which the solid solver marches in time with unchanged functions q ref (t), T ref (t) and h(t).

The aim of the game is then to predict the transient evolution of the temperature in the solid with an imposed (known) evolution of the operating conditions X(t). The time step in the solid computation is imposed by accuracy/stability considerations. The number of fluid computations used to update the unsteady boundary condition (39) on the fluid interface must be limited in order to minimize the overall computational time. To achieve this, the total time to be simulated is subdivided into several transient coupling intervals. A coupling interval is an elemental piece of the coupled simulation bounded by a starting time t 1 and an ending time t 2 . At the initial time t 1 , the instantaneous temperature field in the solid T s (t 1 ) is known, especially at the interface T is (t 1 ). Operating conditions are known at any time and especially in the interval [t 1 , t 2 ]. If the evolutions of the coefficients q ref (t), T ref (t) and h(t) of the boundary condition (39) are known from t 1 to t 2 , the solid computation can march in time till t 2 after several time steps. The new temperature field T s (t 2 ) is then known and especially at the interface T is (t 2 ). But the evolutions of the coefficients q ref (t), T ref (t) and h(t) are deduced from fluid computations performed only at the bounding times t 1 and t 2 . It is then only possible to estimate (q ref (t 1 ), T ref (t 1 ), h(t 1 )) and (q ref (t 2 ), T ref (t 2 ), h(t 2 )). The simplest way to deduce a continuous model of these coefficients during the coupling cycle is to perform a linear interpolation between the bounding values. Another issue is that T is (t 2 ) can be computed by the transient solid computation, once the evolution of (39) is known, and that latter can be known only after that fluid computations have been performed in the vicinity of T is (t 2 ). To sum up, T is (t 2 ) is the solution of an implicit equation of the form:

T is (t 2 ) = f (T is (t 2 )) ( 40
)
where f represents the process beginning with an evaluation of the boundary condition (39) at t 2 thanks to steady fluid computations around T is (t 2 ), and proceeding with a 5) a transient solid computation is performed with the linearized boundary condition (39) and coefficients (46), ( 47) and (48), it gives a new estimation of the evolution of the temperature in the solid T n s (t) and in particular a new estimation of the temperature at the final time T n s (t i+1 ); 6) the new temperature distribution in the solid at the final time T n s (t i+1 ) is compared to the previous estimate T n-1 s (t i+1 ): if their difference is less than the expected accuracy |[T s (t i+1 )] n (P ) -[T s (t i+1 )] n-1 (P )| < δT at any point P of the discretized solid S, the iterative process is stopped, otherwise it is continued.

To be clear on the vocabulary chosen here, a coupling iteration or coupling cycle refers to the sequence of the six operations described above. To sum up, the overall transient simulation, from t 0 to t N I , is carried out by a series of N I transient coupling intervals [t i , t i+1 ], i = 0, N I . The temperature field in the solid, T s (t), t ∈ [t i , t i+1 ], and the thermal boundary condition on its interface with the fluid, (39) with (q ref (t), T ref (t), h(t)), t ∈ [t i , t i+1 ], are determined by an iterative process in six steps. Since this iterative process cyclically travels in time on the interval [t i , t i+1 ], each one of its iterations is called a "coupling cycle". Several cycles are necessary to get an assessment of T s (t) and (q ref (t), T ref (t), h(t)) on [t i , t i+1 ] which satisfies the expected accuracy δT on the temperature field.

When the iterative process is stopped by the convergence check, step 6:

-The last computed fluid flow field satisfies the steady governing equations. On the interface, it gives (T n if , q n f ).

-The temperature field in the solid T n s (t i+1 ) is the result of the unsteady internal energy equation applied to the solid, the initial state T s (t i ), the boundary conditions depending on the known evolution of the operating conditions X(t), t ∈ [t i , t i+1 ] and the boundary condition (39) on the interface, the coefficients of which are considered to vary approximately linearly in time.

-On the interface I, the temperature is continuous within the accuracy δT :

∀P ∈ I : |[T is (t i+1 )] n (P ) -[T if (t i+1 )] n (P )| < δT (49) because ∀P ∈ I : |[T is (t i+1 )] n (P ) -[T is (t i+1 )] n-1 (P )| < δT (step 6) and [T if (t i+1 )] n = [T is (t i+1 )] n-1 (step 1)
. The heat flux is conserved within the accuracy :

∀P ∈ I : |[q s (t i+1 )] n (P ) -(-[q f (t i+1 )] n (P ))| < h n (t i+1 ).δT (50)
because (39) is satisfied by the solid computation step 5, especially at t i+1 , with the coefficients set at step 3, and with (49).

Of course, the larger the difference between t 1 and t 2 , the higher is the approximation due to the time linearization of the coefficients (step 4) of the linearized boundary condition (39).

5) a transient solid computation is performed with the linearized boundary condition (53) and coefficients (59), ( 60) and (61), it gives a new estimation of the evolution of the temperature in the solid T n s (t) and in particular a new estimation of the temperature at the final time T n s (t i+1 ).

The choice of a linear interpolation of the coefficients T n ref (t), q n ref (t) and h n (t) is the simplest one to get an approximate estimation of the instantaneous boundary condition (39), the error of this approximation vanishing when t i+1 → t i . But it is not the simplest choice to approximate the instantaneous boundary condition (53). For that latter, the simplest choice is to directly interpolate its own coefficients h n (t) and T n ref (t):

3) the values of the coefficients of the relaxed boundary condition (53) at t i+1 are updated:

[h(t i+1 )] n = - [q f (t i+1 )] n,bis -[q f (t i+1 )] n [T if (t i+1 )] n,bis -[T if (t i+1 )] n (63) and [ T ref (t i+1 )] n = [T if (t i+1 )] n - [q f (t i+1 )] n [h(t i+1 )] n (64)
4) the evolutions of the coefficients of the relaxed condition (53) are approximated by linear functions between the known values at t i and the estimated values at t i+1 :

h n (t) = t i+1 -t t i+1 -t i h(t i ) + t -t i t i+1 -t i [h(t i+1 )] n (65)
T n ref (t) = t i+1 -t t i+1 -t i T ref (t i ) + t -t i t i+1 -t i [ T ref (t i+1 )] n (

Comparison with other works

The strategy of coupling a transient solid computation to several steady fluid computations was already implemented in others works, and recently by Gimenez et al. [START_REF] Gimenez | A coupling numerical methodology for weakly transient conjugate heat transfer problems[END_REF] to simulate the temperature evolution of a flat plate, heated on one side and cooled by convection on the other side, during three hours.

The fundamental difference with the present approach is that the coupling method used in [START_REF] Gimenez | A coupling numerical methodology for weakly transient conjugate heat transfer problems[END_REF] and described by Errera et al. [START_REF] Errera | A quasi-dynamic procedure for coupled thermal simulations[END_REF] starts from a relaxed mathematical relation, between quantities at the interface of both solvers, which was devised to make the heat fluxes match q ∞ s = -q ∞ f when the temperatures match T ∞ is = T ∞ if at convergence. This relaxed condition, adapted to the notations used in the present text, the convention of positive outgoing fluxes and the indexation of the steps during a coupling cycle, reads:

q n s = -q n f + h n .(T n is -T n if ) (67)
where h n > 0 is a mathematical relaxation parameter.

The approach described herein starts from a local linearization of the heat flux with respect to the temperature at the interface of a steady flow. This allows to treat cases where the function q f (T if ) is globally not linear. Based on the slow solid -fast fluid assumption, and thanks to a second linearization in time, this relation is directly imposed to the solid in transient evolution. Because of the fixed point method chosen to solve this implicit problem, the present approach also leads to a relaxed relation between fluid and solid quantities at the interface, but without any indetermination of the coefficients.

Coefficients are related to physical quantities solved during the coupled simulation. In practice, this means that the spatial distribution of a relaxation parameter must be arbitrarily set by the user of the mathematical method described in [START_REF] Errera | A quasi-dynamic procedure for coupled thermal simulations[END_REF], while nothing has to be arbitrarily tuned in the present method.

In the academic configuration of the flat plate simulated by Gimenez et al. [START_REF] Gimenez | A coupling numerical methodology for weakly transient conjugate heat transfer problems[END_REF], it is trivial to calculate the heat transfer coefficient h f , because for the uniform low-Mach flow, the reference temperature T ref in (53) is fairly approached by the far-field static temperature T ∞ . Taking advantage of this simplicity, Gimenez et al. use a monolithic FLUENT computation (both solid and fluid are computed by only one conjugate heat transfer solver) to get the heat flux distribution, then the h f distribution. Their parametric study on the relaxation parameter h n (denoted α n f in [START_REF] Gimenez | A coupling numerical methodology for weakly transient conjugate heat transfer problems[END_REF]) allows to highlight that indeed the natural choice of h n = h f leads to the best/fastest convergence of the coupled computation of the segregated fluid and solid solvers compared to any other artificial choice. But, this simple case does not need a coupled method to be solved. In the real world, we want to avoid the accurate monolithic computation because of its cost on the target applications. That is exactly the raison d'être of the approximate coupled method. The other motivation is not to duplicate the thermal computation in solids which is already performed by solvers dedicated to thermal and mechanical analysis on complex structures/mechanisms. Practical applications have generally heterogeneous flows where it is not possible to define a location were a value of T ref could be probed. The value of the heat transfer coefficient of the fluid h f (in the general sense of h f = -∂q f /∂T if ) is the only sensible choice to set the relaxation parameter h n in (39) or (53). If an engineer had to use the mathematical method proposed by Errera et al. [START_REF] Errera | A quasi-dynamic procedure for coupled thermal simulations[END_REF], he/she would have to set the value of the relaxation parameter α n f at each point of the fluid-solid interface. In reality, the distribution of h f is not known a priori. In the best-case scenario, the engineer will impose crude piecewise constant estimations on the interfaces, inferred from academic cases which are far from the real configurations encountered in practical applications. Anyway, setting up the computation will be tedious.

In the present method, no arbitrary settings have to be done since the correct values of heat transfer coefficient are automatically set. Concerning the convergence of the coupling cycles, the analysis of the basic thermal problem shows that it will be the fastest in simple configurations. No mathematical proof of convergence is available for all possible configurations, then, the only pragmatic thing to do is to test, as it is commonplace in Computational Fluid Dynamics.

Conclusion

This paper proposes a new approach to model the thermal transfer between a solid and a fluid. It is first based on the slow solid -fast fluid assumption which allows to consider the adaptation of the fluid flow as instantaneous compared to the response of the solid at the points of interest (out of this assumption, the method does not apply). The fluid flow can then be solved by a steady solver. This permits to assess the thermal influence of the fluid only at the solid time steps and reduce the overall computational time. In order to reduce this time further, fluid computations can be conducted only at chosen solid time steps, but this requires to interpolate the fluid influence (through the boundary condition on the interface) between coupling dates. Of course, this approximation can hold only for a moderate elapsed time between coupling dates.

The method is based first on the linearization of the relation between the heat flux and the temperature on the wall. Since the flow is assumed to adapt instantaneously to boundary conditions, this relation is imposed by the steady fluid flow to the solid, all other things (i.e. operating conditions) being equal. This allows to propose a coupling method to solve steady states involving solids and fluids.

The approach developed for steady thermal computations can then be extended to transient simulations, by simply replacing the steady solid computation with a transient solid computation. The final approximative ingredient is the linearization in time of the coefficients of the boundary condition on the fluid-solid interface already linearized with respect to the wall temperature.

As far as I know, the transient coupling method proposed in this paper is today the only pragmatic solution to perform thermal analyses on systems involving solids and fluids over long periods of time in an industrial context.

Figure 1 :

 1 Figure 1: Stability domain of the coupling method by T is and (h, T ref ).

  66) 5) a transient solid computation is performed with the relaxed boundary condition (53) and coefficients (65) and (65), it gives a new estimation of the evolution of the temperature in the solid T n s (t) and in particular a new estimation of the temperature at the final time T n s (t i+1 ).

  1 if (P ) -T n if (P )| ≤ δT 2 with δT 2 ≥ δT ensures that δT ≤ |T n,bis if (P ) -T n if (P )| ≤ δT 2 . A high value of δT 2 > δT allows to avoid the second computation, to save some time. If δT 2 = δT , the second fluid computation is always performed, imposing T n,bis if

The discussion is here only on phenomenology, but in the numerical implementation of the equations of fluid, it is common to use the equation of the total energy E = e + V

/2 and gather the flux of total energy ρEV n and the flow work pV n to form the flux of total enthalpy ρHV n , with H = E + p/ρ.

At each iteration n ≥ 1 of the coupling procedure, the distributions on the interface T n-1 is and (T n-1 if , q n-1 f

) provided respectively by a steady thermal solid computation and a steady fluid computation are supposed to be known from the previous iteration n -1.

The new estimation at the coupling iteration n is computed in five steps:

1) a steady fluid computation is performed with the imposed temperature T n if = T n-1 is on the interface and gives a new estimation of the heat flux on the interface q n f .

2) if ∀P : < |T n-1 if (P ) -T n if (P )| ≤ δT 2 , then T n,bis if = T n-1 if and q n,bis f = q n-1 f , otherwise an additional steady fluid computation with the imposed temperature T n,bis if = T n if + δT on the interface is performed, providing the heat flux q n,bis f .

3) the coefficients of the linearized boundary condition (21) are updated:

4) a steady solid computation is performed with the linearized boundary condition (21) and coefficients (26), it gives a new estimation of the temperature in the solid T n s ;

5) the new temperature distribution in the solid T n s is compared to the previous estimate T n-1 s : if their difference is less than the expected accuracy |T n s (P ) -T n-1 s (P )| < δT at any point P of the discretized solid S, the iterative process is stopped, otherwise it is continued.

To be clear, a coupling iteration refers here to one sequence of these 5 operations.

In order to start the iterative process, the initial data T 0 is and (T 0 if , q 0 f ) must be set. Several choices are possible, but since it is crucial to respect the condition (25), two strategies will here be described. If a fair estimation of the temperature on the interface is known, it can be used as the initial distribution T 0 is . Then, choosing T 0 if = T 0 is + 2.δT 2 will force the second fluid computation, at step 2 of the first iteration, with a temperature distribution T 1,bis if = T 1 if + δT on the interface. The distribution q 0 f can be initialized with an arbitrary field, since it will not be used. If no estimation of the temperature on the interface is available, one can perform an adiabatic (q 0 f = 0) fluid computation which will provide a temperature distribution T 0 if on the interface. Then, choosing

and δT 2 > δT will avoid the second fluid computation at step 2 of the first iteration. With both strategies, (25) will be satisfied and the first derivative of q f = g(T if ) in h 1 (26) will be calculated with the prescribed precision δT .

At step 2, the condition |T

is (P )| ≥ δT (otherwise the iterative process would have been stopped at previous iteration) and the condition 3 Quasi-transient thermal coupling

introduction

The generalization of the coupling method presented for steady thermal fluid-solid problems (section 2.1) relies on the fast fluid -slow solid assumption (section 1.4). From the viewpoint of the fluid, the temperature in the solid varies so slowly that the boundary condition on the solid interface can be considered as stationary. From the viewpoint of the solid, the thermal field of the fluid, governed by its dynamics and thermodynamics, adapt so quickly to variations in its boundary conditions that, at the time-scale of the thermal response at the points monitored in the solid, the fluid seems to reach steady states instantaneously, and then to impose a boundary condition on the fluid interface of the solid that depends only on the instantaneous thermal quantities (T is , q s ). In the case of the internal gas flowing inside gas turbine engines, the boundary conditions in a given region of the engine are, for instance, the inlet total pressure(s) and total temperature(s), the outlet static pressure(s) and the rotation rate(s) of some walls. These parameters are globally called "operating conditions" in the remaining of this text and denoted X = (X p ) p=1,Np . The characteristic variation time of the operating conditions is indeed significantly larger than the acoustic and convective times of the gas in engines. That is why the thermal behavior of the fluid on the fluid-solid interface can be formally expressed without an explicit dependence on time, as:

where g is a function which is not necessarily linear with respect to any of these variables.

For given operating conditions X, the heat flux q f (T b if , X) at a temperature T b if close to a reference temperature T a if can be approached by Taylor series:

, the non-linear term can be neglected so that a boundary condition of the following form can reproduce the influence of the fluid on the solid:

In fact, during a transient regime, the operating conditions X are functions of time t. Nonetheless, in the present analysis, the function X(t) is considered as fixed input data.

In other words, the influence of the thermo-mechanical state of the solid on the operating conditions (like the influence of the thermal modification of blade-tip gaps on the pressure ratio and shaft work in a compressor or turbine) is not taken into account in this study.

Knowing the time evolution of the operating conditions X(t), (38) can be reformulated:

Rigorously, the fast fluid -slow solid allows to separate the fluid computation from the solid computation, because the fluid dynamics adapt far faster than the temperature in transient solid computation from t 1 to t 2 based on the linear interpolation of (39) between t 1 and t 2 , which allows to get another evaluation of T is (t 2 ). This problem can be solved iteratively, starting from an initial distribution T is (t 2 ) 0 and applying the recursive procedure:

Provided that a temperature accuracy δT has been chosen, the iterative process can be stopped when:

In fact, the solution to solve this transient coupling problem is very similar to that already described to solve the steady coupled problem in section 2.1. It is a fixed point iteration method. In the steady coupled problem, one has to find the solution of an implicit equation of the form:

where T ∞ is is the temperature on the fluid-solid interface in the steady state, i.e. after that the iterations of the steady solid solver has reached convergence, and f represents the sequence of operations detailed in section 2.1 (not the same f as f in (40)). Here, starting from the initial state T s (t 1 ) in the solid and applying a prescribed thermal loading (q n ref (t), T n ref (t), h n (t)) on the interface (conditions on other boundaries of the solid can also change in time, this influence will here be represented by the dependence on the operation conditions X(t)), the unsteady solid solver marches in time till t 2 , where it provides, in particular, a new estimate of

) is assessed thanks to a previous estimation [T is (t 2 )] n-1 , the new estimate [T is (t 2 )] n allows to update the estimation of the coefficients, that is why (40) is solved by an iterative process. For comparison, in the steady problem, a previous estimation

is allows to set the coefficients (q n ref , T n ref , h n ) of the linearized boundary condition (21). However, no proof is here provided that the iterative process in the case of an unsteady thermal evolution in the solid converges. In practice, there is no other method to solve a long lasting coupled thermal fluid-solid problem, therefore the only way is to test.

General quasi-transient method

After this introduction, a general procedure to solve a quasi-transient coupling between a solid and a fluid flow can be described as follows. The total period of time [t 0 , t N I ] to simulate must first be subdivided into N I time intervals [t i , t i+1 ], i = 0, N I . The overall simulation will then consist in a sequence of transient coupling intervals [t i , t i+1 ], i = 0, N I . For each transient coupling interval, one has to choose the accuracy δT > 0 expected on the temperature field, a tolerance δT 2 ≥ δT for the calculation of the heat transfer coefficient h n based on (24) and a small real parameter > 0 preventing zero division ( < δT ). The simplest choice is to set these parameters to constant values for all the intervals.

In each coupling interval [t i , t i+1 ], the evolution of the temperature in the solid is solved iteratively as presented in the introduction, by trial and error on the final temperature distribution on the interface T is (t i+1 ). This procedure can hold only if the slow solid-fast fluid assumption applies. At the beginning of a coupling interval, the initial temperature field in the solid T s (t i ) is known, especially at the interface ∀P ∈ I : [T is (t i )](P ) = [T s (t i )](P ), and the coefficients (q ref (t i ), T ref (t i ), h(t i )) of the linearized boundary condition (39) are known. At any time, the operating conditions X(t) are known and determine some of the fluid and solid boundary conditions. Since time is cyclically traveled from t i to t i+1 , each iteration of this trial and error process is called a "coupling cycle". At each coupling cycle n ≥ 1, the distributions on the interface [T is (t

provided respectively by a transient thermal solid computation and a steady fluid computation are supposed to be known from the previous coupling cycle n -1. They are an estimation of the unknown distributions T is (t i+1 ) and (T if (t i+1 ), q f (t i+1 )). At the coupling cycle n, the estimation of the temporal evolution of the temperature field in the solid T s (t) in the coupling interval [t i , t i+1 ], denoted T n s (t), associated to the linear estimation of the temporal evolution of the coefficients (q

), in turn associated to the final distribution T is (t i+1 ) estimated at the previous cycle [T is (t i+1 )] n-1 , is computed in 6 steps:

1) a steady fluid computation is performed with the boundary conditions adapted to the operating conditions X(t i+1 ) and with the imposed temperature

on the interface, it gives a new estimation of the heat flux on the interface [q f (t i+1 )] n at the final time t i+1 of the interval.

2) if ∀P

and [q f (t i+1 )] n,bis = [q f (t i+1 )] n-1 , otherwise an additional steady fluid computation with the imposed temperature [T if (t i+1 )] n,bis = [T if (t i+1 )] n + δT on the interface is performed, providing the heat flux [q f (t i+1 )] n,bis .

3) the values of the coefficients of the linearized boundary condition (39) at t i+1 are updated:

4) the evolutions of the coefficients of the linear boundary condition (39) are approximated by linear functions between the known values at t i and the estimated values at t i+1 :

Concerning the initialization, one has to distinguish the initialization of the overall transient coupled simulation (at t 0 ) and the initialization of the recursive procedure to solve each transient coupling interval [t i , t i+1 ] (at t i ).

Both require an initial temperature field T s (t i ) in the solid which is compliant with both the boundary conditions adapted to the operating conditions X(t i ) and the thermal condition imposed by the fluid flow on the interface (the boundary conditions of the fluid flow being also adapted to the operating conditions X(t i )). During the overall simulation, these data are provided to the current interval computation by the results of the previous interval computation (after convergence of several coupling cycles). At the first interval, these data can be either the result of a preliminary overall transient coupled simulation (in this case, the temperature field in the solid has not necessarily converged toward a steady state, but it is necessarily compliant on its interface with a steady fluid flow at the same operating conditions) or the result of a preliminary steady coupled simulation, whose procedure is described in section 2.1.

The recursive procedure to solve each transient coupling interval [t i , t i+1 ] needs a first estimation of the final distributions on each side of the interface, respectively [T is (t i+1 )] 0 and ([T if (t i+1 )] 0 , [q f (t i+1 )] 0 ). The simplest choice is to set:

In this way, the first estimation of the coefficients (44) and (45) will be calculated from two fluid computations with the respective wall temperatures T is (t i ) and T is (t i ) + δT . An alternative could be to extrapolate those values from two or more previous coupling times t i , t i-1 , ... Similarly, a higher order function-fitting could be used to interpolate the values of the coefficients of (39) in [t i , t i+1 ].

To finish, the convergence check, step 6, should be performed on the time and spatial distribution T n s (P, t), P ∈ S and t ∈ [t i , t i+1 ] but this would be too much time consuming.

Adaptation to usual solvers

In this section, it is shown how the general method described in 3.2 can be adapted to usual solid solvers. As said in section 2.2, the linearized boundary condition (39) is often not available in classical finite-element solvers used for thermal analysis. Only classical conditions like an imposed temperature distribution T is (t), an imposed heat flux distribution q is (t) or a relaxed boundary condition of the form:

where h(t) and T ref (t) are space distributions varying in time, can be applied on the solid boundaries. A relaxed boundary condition like (52) can still be used in the transient solid computation:

21 if its coefficients h n (t) and T n ref (t) are updated at each coupling iteration n so that they reproduce the linearized approximation of the fluid behavior (37). This latter can be reformulated as:

where q s = -q f (T b if , X) and T is = T b if , and under the condition that dq f

dT if (T a if , X) = 0. The identification of the remaining terms in (53) and (54) leads to:

During the iterative process, h n (t) and T n ref (t) can then be updated in (53) in the same way as T n ref (t), q n ref (t) and h n (t) are updated in (39) of the general method. The procedure described in section 3.2 and all the comments still hold, the initialization process included. It is just necessary to modify steps 3, 4 and 5. We will entirely rewrite these three steps in order to gather the information necessary to program the method, even if there are few changes. While keeping a linear interpolation of the coefficients in time, there are two possibilities in the choice of the coefficients to interpolate. If we want to reproduce exactly the same evolution as that described in section 3.2, step 3 can be unchanged, and the calculation of h n (t) and T n ref (t) must be inferred from T n ref (t), q n ref (t) and h n (t) after these latter have been interpolated at step 4:

3) the values of the coefficients of the linearized boundary condition (39) at t i+1 are updated:

and [h(t i+1 )] n = -[q f (t i+1 )]

then, the coefficients of (53) are deduced from that of (39), h n (t) is the same and: