M-P Errera
email: marc.errera@onera.fr

R El Khoury
email: roch.khoury@gmail.com

NUMERICAL STUDY OF TWO OPTIMIZED COUPLING INTERFACE TREATMENTS FOR STEADY CONJUGATE HEAT TRANSFER PROBLEMS

Keywords: Conjugate heat transfer, stability, optimal, fluid-solid, Dirichlet, Neumann, Robin

This paper presents two transmission interface treatments, Dirichlet-Robin and

Neumann-Robin procedures, that may be employed for conjugate heat transfer problems. These conditions are analyzed on the basis of a 1D simplified model problem. In the first part of the paper, the Dirichlet-Robin procedure is presented. This interface treatment is the most widely employed in the literature. The same analysis is then performed with a Neumann-Robin procedure. On the basis of the model problem, the general expression of the amplification factor, the stability bounds and the optimal coefficients are provided. It is shown that the two interface treatments are opposite and complementary. Moreover, the so-called optimal coefficient provides the best results in terms of stability and convergence in the Dirichlet-Robin procedure. A criterion is expressed to choose the most appropriate transmission procedure and its importance is underlined by a test case.

INTRODUCTION

Conjugate heat transfer (CHT) analysis refers to the ability to address the thermal interaction between a body and a fluid flowing over or through it. The objective of CHT modeling is the accurate prediction of temperature and heat flux distribution in space and time in a body and on its boundaries. Conjugate heat transfer problems occur whenever the fluid convection and the solid material conduction are taken into account simultaneously. This mutual interaction is becoming increasingly important because in many numerical simulations it is no longer acceptable to consider heat transport in the fluid only, with some ideal conditions at the fluidsolid interface.

Even if coupled procedures are the most correct and reasonable way to use when accurate heat transfer predictions are needed, in most cases, arbitrary relaxation parameters are used to stabilize the coupling procedure. This may have a significant negative impact on the numerical properties (low convergence rate, oscillations, and instabilities). Our goal in this paper is to present two interface conditions derived from a stability analysis and to present their main numerical characteristics.

Recently, we have shown [START_REF] Errera | Optimal solutions of numerical interface conditions in fluid-structure thermal analysis[END_REF] using a 1D thermal model problem that in a coupled system, a numerical transition can be identified. This fundamental result has been derived from a normal mode stability analysis based on the theory of Godunov-Ryabenkii [2][3]. This transition can be regarded as an optimal choice in terms of stability and convergence. Ideally, a Robin condition on each side of the interface should be considered. But in this case, a general Robin-Robin interface condition leads to a too large family of schemes. Thus, in this paper we will confine ourselves to two commonly used conditions :

-Dirichlet-Robin procedure : the temperature coming from the solid is applied on the fluid side and a "relaxed heat flux" is in turn used as a boundary condition for the solid.

-Neumann-Robin procedure : the heat flux coming from the solid is applied on the fluid side and a Robin condition is in turn used as a boundary condition for the solid.

These conditions have been chosen in order to deal with most situations likely to arise. They will be analysed on the basis of a 1D coupled aerothermal model problem. In the first part of the paper, Dirichlet-Robin interface treatments are presented. This transmission interface treatment is the most widely employed in the literature. The same analysis is then performed in the second part with a Neumann-Robin procedure. These two treatments will then be compared.

MODEL EQUATIONS

Model problem

It has been often been stressed that the nature of instabilities derived from a 1D model can give insight into the potential instabilities in 2D/3D computations. As a result, the behavior of interface conditions in CHT is often studied using a normal mode analysis. This is 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. A 1D model is composed of two partitions with a shared interface. Partitioned techniques [START_REF] Felippa | Staggered transient analysis procedures for coupled mechanical systems: Formulation[END_REF] [START_REF] Felippa | Partitioned analysis of coupled mechanical systems[END_REF] are very popular, because they allow the direct use of specifically designed solvers for different fields and may offer significant benefits in terms of efficiency over the monolithic techniques. In this study, the fluid-solid system will be decomposed into two partitions: the (+) subdomain

(Ω + , index 0 ≥ j ; 0 ≥ x ) and the (-) subdomain (Ω -, index 0 ≤ j ; 0 ≤ x
) as illustrated in Figure 1.

Interface conditions on the solid side

Interface conditions are needed on either side of the shared interface (x = 0), where coupling conditions are applied. Our goal is to ensure a stable CHT process and to avoid destabilizing effects. It is well known that Robin conditions have many attractive features since such a condition introduces an interface stiffness forcing the boundary to behave in the same way as the boundary of the other domain resulting in much better stability properties. A Robin condition on the solid side is written simply as

[ ] [ ] f f f s f s T Q T Q α α + = + ˆ (1)
The subscripts f and s denote the fluid and solid domain respectively and the (^) notation indicates the sought values. Q is the interfacial heat flux (W.m -2 ) and T is the interface temperature. The general Robin condition (1) introduces the numerical coupling parameter f α (W.m -2 .K -1 ) the choice of which directly influences the stability of the coupling process.

Interface conditions on the fluid side

A similar expression could be used on the fluid side [START_REF] Godunov | The theory of difference schemes. An introduction[END_REF]. We readily obtain

[ ] [ ] s s s f s f T Q T Q α α + = + ˆ (2)
s f Q Q = ˆ (4) External BC Fluid (+) Solid (-) ∆x f -J -2 -1
Interface "0"

Equations in the fluid and solid domains

The CHT strategy used to obtain rapidly a fluid-structure steady solution relies on the significant discrepancy of the characteristic physical times of the two domains, namely a fast transient process in the fluid, a very slow one in the structure. In the fluid subsystem, the Navier-Stokes (NS) equations are solved generally to steady-state by a time marching scheme. As a consequence, a time marching scheme will be employed here in the fluid domain in the CHT model.

On the contrary, in the solid domain, if we are only interested in computing the steadystate solution, then a second order ordinary differential equation can be solved directly and coupled to the unsteady fluid domain. Thus, only a Laplace equation is considered in the solid.

Characteristic factors

Three dimensionless numbers play a major role in CHT problems and appear in the stability analysis performed in this paper.

The first one is the mesh Fourier number defined as follows

2 f f f f x t a D ∆ ∆ = (5)
where f a is the fluid diffusivity (m 2 .s -1 ). f D is the mesh Fourier number that characterizes transient heat diffusion in the fluid domain.

The second one is a conventional Biot number defined by

s K h Bi = (6)
where h (W.m -2 .K -1 ) is the heat transfer coefficient and s K is the thermal conductance of the solid. Bi determines whether or not the temperatures inside a solid body vary significantly in space, while the body heats or cools over time, from a thermal gradient applied to its surface.

The third number, a local Biot number, naturally appears in the stability analysis performed hereafter and takes the following form

s f K K Bi = ∆) ( (7) ) (∆
Bi is a "new" parameter that can be regarded as the ratio of the thermal conductance of the 1 st fluid cell over the thermal conductance of the whole solid domain. This number represents the balance at each coupling time between the solid resistance of the entire solid body (only a steady-state is considered in this domain) and the resistance of the transient fluid domain represented by the diffusion in the first near wall interface.

It is also interesting to introduce a normalized Fourier number f D (see [START_REF] Errera | Comparative study of Coupling Coefficients in Dirichlet-Robin Procedure for Fluid-Structure AeroThermal Simulations[END_REF] for more details) defined by

f f f f D D D D 2 1 1 + + + = (8) The domain of f D is 0 ≥ f D and the range is 1 0 < ≤ f D .

Normal mode analysis

We have performed a stability analysis by considering normal mode solutions in the form

     ≤ > = - 0 , 0 , 1 j z j z T j s n j f n n j κ κ (9)
where z and κ are two complex functions. z is the "temporal amplification factor" and κ is the "spatial amplification factor". The stability analysis is very similar to the standard Fourier stability method except that the Fourier analysis ignores boundary conditions and as these may affect stability, the theory of Godunov & Ryabenkii [START_REF] Godunov | The theory of difference schemes. An introduction[END_REF][3] is preferable. After elementary transformations not reported in this paper (see details in [START_REF] Errera | Optimal solutions of numerical interface conditions in fluid-structure thermal analysis[END_REF]), we find that the temporal amplification factor z (i.e. each mode increases in amplitude by a ratio z ) takes the following forms in the two interface conditions considered in this paper :

(1) Dirichlet-Robin interface condition : ( ) ( )

∞ ≥ = , 0 , f s f α α α
The temporal amplification factor can be written as

[ ] ) ( ) , ( 1 ) , ( f f f f f s f K z K z g z α α κ α α - - + = = (10) (2) Neumann-Robin interface condition : ( ) ( ) 0 , 0 , > = f s f α α α
The temporal amplification factor can be written as

      - + + = = ∆) ( ) , ( . 1 ) , ( Bi K z K z g z f f f f f f s f α α κ α α α (11)

A fundamental transition

The complex function g is a rather complicated nonlinear equation for z. After the change of variable z z 1 →

, this function becomes holomorphic on the open set 1 < z

. As a result, on the basis of the maximum modulus principle in complex analysis, the maximum value of g is achieved on the boundary, i.e., at some point on the unit circle 1 = z .

Actually, the maximum is obtained either at 1

+ = z or 1 - = z
. But, under certain conditions, there is a sudden transition from one point to another resulting in an amplification factor composed of two half-lines with a singular point at the intersection of these two lines. At this intersection the amplification factor turns back and attains its absolute minimum always located in the stable zone 1 < z .

The point where the maximum is transferred from

1 - = z to 1 + = z
is a fundamental transition in the aerothermal coupling.

DIRICHLET-ROBIN INTERFACE CONDITIONS

Stability bounds

The stability condition 1 ) , ( < f z g α applied to Eq.(10) leads, after some basic calculus manipulations, to a lower stability bound min

f α > f α 2 ) 1 ( 2 min s f f f K D K - - = α (12)

A dynamical Biot number

Substituting the definition of ) (∆ Bi (Eq.( 7)), min

f α becomes min f α = [ ] 1 ) 1 ( 2 ) ( - - ∆ f s D Bi K (13)
As the coupling coefficient is always positive, two zones are clearly identified :

• 1 ) 1 ( ) ( ≤ - ∆ f D Bi : the coupling process is stable 0 ≥ ∀ f α • 1 ) 1 ( ) ( ≥ - ∆ f D Bi
: the coupling procedure exhibits a lower stability bound min f α This demonstrates how stability depends mainly on the ratio of thermal resistances, but also on the dynamics of the transient fluid system. The higher the local Biot number, the more difficult it will be to stabilize the coupling. Of course, the transient effects of the fluid system also play a major role and can stabilize any process via the term ) 1 ( f D -, but at a significant cost.

Optimal procedure

It is very noteworthy that the modulus of the amplification factor does not have a monotonic variation (in terms of f α ) but goes through an absolute minimum, denoted opt f α . In other words, the existence of a transition value for f α can be identified. At this transition value, the shape of the curve of the amplification factor switches and turns back as can be seen in Figure 1. In this Figure, the two curves have been plotted using the following characteristics :

• Local Biot number 2 ) ( 10 25 . 1 = ∆ Bi • Mesh Fourier numbers : (1) ⇒ = 40 f D 8 . 0 = f D (2) ⇒ = 000 , 40 f D 992 . 0 = f D Thus, the stabilizing effect of the term ) 1 ( f D - is noticeable. The most striking point is that at the remarkable coefficient opt f f α α =
, the modulus of the amplification attains its absolute minimum value, always less than one.

Low values of

f α ) ( opt f f α α <
will result in a rapid convergence, but in this case however f D must be sufficiently large to allow heat diffusion on the fluid side. Otherwise, a low Fourier number will soak up a lot of heat. It will be then necessary to enhance stability by increasing the coupling parameter. is the optimal choice in the case of the model equation adopted herein. .000

Dirichlet-Robin procedure

Summary

The general behavior of the Dirichlet-Robin coupling procedure for steady CHT, in terms of f α is illustrated in Table 1. --------------opt f α -----------------

f α 0 min f α -
- ∞ c o n d i t i o n g U N S T A B L E 1 opt g 1 1 ) 1 ( ) ( ≥ - ∆ f D Bi g 0 g opt g 1 1 ) 1 ( ) ( ≤ - ∆ f D Bi Table 1: Numerical behavior of Dirichlet-Robin procedure vs f α ( ) 0 ; ≥ ∞ = f s α α

NEUMANN-ROBIN INTERFACE CONDITIONS

Stability bounds

The stability condition 1 ) , ( < f z g α applied to Eq.( 10) leads, after some basic calculus manipulations, to an upper stability bound max

f α < f α ) 1 ( 2 max f f s f s f D K K K K + - = α (14)
Substituting the definition of ) (∆ Bi (Eq.( 7)), max

f α becomes max f α = ) 1 ( 1 2 ) ( f f D Bi K + -∆ (15)
As the coupling coefficient is always positive, two stability regions are highlighted :

• 1 ) 1 ( ) ( ≥ + ∆ f D Bi : the coupling process is stable 0 ≥ ∀ f α • 1 ) 1 ( ) ( ≤ + ∆ f D Bi
: the coupling procedure exhibits an upper stability bound max f α

Optimal procedure

This time, in the framework of the Neumann-Robin procedure, three zones can be identified. The first zone exhibits an upper stability limit as just mentioned. There is also in this zone an optimal coefficient for which the modulus of the amplification factor attains an absolute minimum.

The second zone is very narrow. It also presents an optimal coefficient, but in contrast to the previous case, there is no stability bound.

The third zone is also unconditionally stable but no optimal coefficient can be defined. That is to say the minimum of the amplification factor is obtained for ∞ = f α , i.e. for the Neumann-Dirichlet procedure (heat flux imposed on the fluid side & temperature imposed on the solid side) with no relaxation. Another way of explaining it, is that the amplification factor is a monotonic function.

Summary

The general behavior of the Neumann-Robin coupling procedure for steady CHT, in terms of f α is illustrated in Table 2. ---------------max f α ------------------∞ c o n d i t i o

f α 0 opt f α -
n g 1 opt g 1 U N S T A B L E 1 ) 1 ( ) ( ≤ + ∆ f D Bi g 1 opt g ) (∞ g 2 ) 1 ( 1 ) ( ≤ + ≤ ∆ f D Bi g 1 ) (∞ g 2 ) 1 ( ) ( ≥ + ∆ f D Bi

COMPARISON OF THE TWO INTERFACE CONDITIONS

Stabilizing effect of the Fourier number

As just stated, the two interface conditions considered in this paper may be regarded as two complementary conditions, and it seems that we can go as far as to say that they are almost opposite. The first condition is unconditionally stable at small Biot number while the second one is unconditionally stable at large Biot numbers. But there is one major difference. In the case of the Dirichlet-Robin procedure, the term ) 1 ( f D may become as small as necessary and as a result the unconditional stability condition 1 ) 1 (

) ( ≤ - ∆ f D Bi
(2 nd row of Table 1) can always be satisfied. In other words, the lower stability bound can always be removed by an appropriate choice of f D . Conversely, in the condition 1 ) 1 (

) ( ≤ + ∆ f D Bi (1 st row of Ta- ble 2), the term ) 1 ( f D +
has no stabilizing effect.

Overlapping zone

It is not quite accurate to say that the two interface conditions are opposite. Actually they exhibit an overlapping zone where both may be considered as stable. On the basis of the stability limits presented in the preceding sections, it is seen that there is an overlapping zone where both procedures are unconditionally stable and defined by

f f D Bi D - ≤ ≤ + ∆ 1 1 1 1 ) ( (16)
This overlapping zone gets narrower as the normalized Fourier number f D gets smaller. On the contrary, this zone is significantly extended for large Fourier numbers.

Existence and non-existence of optimal procedures

The two interface conditions presented in this paper exhibit 5 different zones (2 in Table 1 & 3 in Table 2). In one of them (3 rd row of Table 2), there is no optimal coefficient. This does however not mean that the coupling procedure will converge slowly but simply that the condi-

α ( ) 0 ; 0 ≥ = f s α α tion ∞ = f α
is a relevant choice. In other words the Neumann-Dirichlet transmission is recommended. In the first procedure, the optimal choice requires a certain amount of relaxation and this time a Dirichlet-Robin procedure is recommended in conjunction with the use of the optimal coefficient. Actually, there exist two types of optimal procedures. The first one is not subject to any condition and always exhibits an optimal coefficient. The second one presents a conditional optimal coefficient. Clearly, the transient effects of the fluid system play a major role in the stability analysis.

Criterion for a relevant interface treatment

We have just seen how to optimize a specific transmission condition procedure. But this study suggests that it can be interesting to select adequately the most relevant procedure and not necessarily to optimize a particular one, more or less arbitrarily chosen. Two options might be considered :

(1) A "dynamic" option The dynamic option is supported by the fact that the optimal coefficient of the Dirichlet-Robin condition is able to stabilize any procedure. Bi is very high, for instance in CFD computations in which a high wall resolution is adopted. In this case, it is always possible to exploit to advantage the term ) 1 (

f D - such that 1 ) 1 ( ) ( ≤ - ∆ f D Bi
but this can have a substantial cost.

Link between the local and the conventional Biot numbers

This local Biot number introduced by Eq. ( 7) takes directly into account the thermal and dynamic response of the boundary layer and directly participates in the stability of the coupled process. If a transient fluid flow is employed during the CHT computation it is this number that will drive and control the stability of the CHT problem.

On the other hand, the conventional Biot defined by Eq.( 6) is a measure of the resistance to heat flow within the solid relative to the resistance presented by the convection processes at the surface. This parameter is also a key parameter, but only at steady state since it determines the stability of the fluid-solid equilibrium. However, it is clear that this conventional Biot number, not easily defined during the transients, cannot be used to set up a numerical CHT procedure, as long as a transient fluid state is involved in this procedure.

At steady state, the thermal fluid resistance is equal to 1/h and thus the various stability criteria and stability bounds presented in this paper may easily be extended to steady state by simply imposing 0 = f D and h K f = .

Numerical example

A complex CHT computation has been performed in a 3D test case representing an industrial reheat furnace with heteroclite fluid-structure interface (details can be found in [START_REF] El Khoury | Development and validation of coupling methods for the modelling of thermal systems[END_REF]). As a first step, the Dirichlet-Robin condition has been used and a systematic and comparative study of various coupling coefficients has been undertaken. This first study showed that the optimal coupling coefficient presented in this paper was by far, the best-performing procedure.

But by looking more closely at the fluid-solid interfaces of the furnace, high local Biot number can be found, for instance This ratio strongly suggests the use of a Neumann-Dirichlet procedure as mentioned in § 5.4. Figure 2 presents the convergence history of these CHT tests. All the tests performed are not indicated in this figure. It is sufficient to mention that the so-called, optimal Dirichlet-Robin procedure is always stable and oscillation-free and clearly the most efficient coupling method. But, as suggested by the criterion proposed in this study, it is instructive to see that a relevant procedure, even with no relaxation parameter can be even more performing. This means that when the local Biot number ) (∆ Bi is high, stability can be obtained with the Dirichlet-Robin condition but at the expense of computational efficiency (CPU time). In this case, a natural physics-based approach is to impose a Neumann condition.

CONCLUSIONS

On the basis of a simplified model problem, the general expression of the amplification factor, the stability bounds, the optimal coefficient and the general numerical behavior have been presented for two different and complementary interface treatments. It was shown that the numerical properties depend on the ratio of fluid and thermal resistances and that the mesh Fourier number plays a crucial role. Furthermore, a local Biot number has been introduced and it was shown that this number "drives" the overall coupling process.

The two interface treatments have been compared and it has been highlighted that these interface treatments are almost opposite and complementary. When the first method is unconditionally stable, the second one exhibits an upper stability bound. When the second method is unconditionally stable, the first one exhibits a lower stability bound. These two interface schemes present an overlapping area where both of them are stable. It is also shown that the so-called optimal coefficient provides the best results in terms of stability and convergence in the Dirichlet-Robin procedure. The numerical criteria establishing the nature and character of the most relevant interface treatment have been identified and expressed. A comparison in a CHT test case has emphasized the importance of a physics-based numerical approach.

Condition ( 2 )

 2 introduces another coefficient s α computed in the solid side and implemented in the fluid side. The simultaneous consideration of two coefficients could be complicated and thus it has been decided in this paper to focus on two complementary treatments :A Dirichlet condition on the fluid side, obtained easily by setting ∞ =

Fig. 1 :

 1 Fig.1 : Schematic of the fluid & solid domains

But

  to an extremely slow convergence since this corresponds to a slow diffusion of heat through the fluid subsystem. But it should be pointed out that relatively large Fourier numbers ( this case will be unnecessarily frozen by f α .It should be noted that these situations might occur in the same coupled computation or even in the same fluid-structure interface. It is the reason why it is crucial to use local coupling coefficients and it has been shown that

Figure 1 :

 1 Figure 1: Amplification factor for

  the overlapping zone, the procedure with the term ) A "static" option The other option does not take into consideration the dynamics of the fluid system ( option takes advantage of the stabilizing effect of the mesh Fourier number in the transient fluid domain. The second option is more secure especially when the ratio ) (∆

Figure 2 :

 2 Figure 2 : Convergence history and comparison between Dirichlet and Neumann procedures

Table 2 :

 2 Numerical behavior of Neumann-Robin procedure vs f


Antoine Tordeux
email: a.tordeux@fz-juelich.de

Guillaume Costeseque
email: guillaume.costeseque@inria.fr

Michael Herty
email: herty@igpm.rwth-aachen.de

Armin Seyfried
email: a.seyfried@fz-juelich.de

From traffic and pedestrian follow-the-leader models with reaction time to first order convection-diffusion flow models

Keywords: First order traffic flow models, micro/macro connection, hyperbolic conservation laws, Godunov scheme, numerical simulation. AMS: 35F20, 70F45, 90B20, 65M12

In this work, we derive first order continuum traffic flow models from a microscopic delayed follow-the-leader model. Those are applicable in the context of vehicular traffic flow as well as pedestrian traffic flow. The microscopic model is based on an optimal velocity function and a reaction time parameter. The corresponding macroscopic formulations in Eulerian or Lagrangian coordinates result in first order convection-diffusion equations. More precisely, the convection is described by the optimal velocity while the diffusion term depends on the reaction time. A linear stability analysis for homogeneous solutions of both continuous and discrete models are provided. The conditions match the ones of the car-following model for specific values of the space discretization. The behavior of the novel model is illustrated thanks to numerical simulations. Transitions to collision-free self-sustained stop-and-go dynamics are obtained if the reaction time is sufficiently large. The results show that the dynamics of the microscopic model can be well captured by the macroscopic equations. For non-zero reaction times we observe a scattered fundamental diagram. The scattering width is compared to real pedestrian and road traffic data.

Introduction

Microscopic and macroscopic approaches for the purpose of vehicular traffic flow modelling have been often developed separately in the engineering community [43,22,[START_REF] Bayen | Lax-Hopf based incorporation of internal boundary conditions into Hamilton-Jacobi equation. Part I: Theory[END_REF][START_REF] Lebacque | Les modèles macroscopiques du trafic[END_REF]. Similar models can also be used in the description of pedestrian dynamics [40, [START_REF] Chraibi | Force-based models of pedestrian dynamics[END_REF][START_REF] Appert-Rolland | Two-way multi-lane traffic model for pedestrians in corridors[END_REF]. Typically, microscopic models are based on the so-called "follow-the-leader" strategy and they are stated as (finite of infinite) systems of Ordinary Differential Equations (ODEs). They are generally based on speed or acceleration functions which depend on distance spacing, speed, predecessor' speed, relative speed and so on. One of the simplest approach is a speed model solely based on the spacing, firstly proposed by Pipes [START_REF] Bürger | On a diffusively corrected kinematic-wave traffic flow model with changing road surface conditions[END_REF] ẋi (t) = W (∆x i (t)), [START_REF] Appert-Rolland | Two-way multi-lane traffic model for pedestrians in corridors[END_REF] where ∆x i (t) = x i+1 (t)x i (t) denotes the spacing between the vehicle (i) to its predecessor (i + 1) and W (•) stands for the equilibrium (or optimal) speed function depending on the spacing. The microscopic models are discrete in the sense that the vehicles or pedestrians i ∈ Z are individually considered. A macroscopic description consider the flow of vehicles or pedestrians (in the following also referred to as agents) as a continuum in Eulerian or Lagrangian coordinates. For instance in the most classical Eulerian time-space framework, the main variables are the density, the flow and the mean speed. The simplest approach is the scalar hyperbolic equation of the celebrated Lighthill-Whitham-Richards (LWR) model [START_REF] Lighthill | On kinematic waves II: a theory of traffic flow on long, crowded roads[END_REF]39]

∂ t ρ + ∂ x (ρV (ρ)) = 0. ( 2
)
Here ρ is the density, V (•) is the equilibrium speed function which is assumed to depend only on the density. The flow f (ρ) = ρV (ρ) is given by the product of the density times the mean speed. The model is derived from the continuity equation for which the flow is supposed in equilibrium. The microscopic and macroscopic models Eq. ( 1) and Eq. ( 2) well reproduce shock-wave phenomena for Riemann problems. Yet such models are not able to describe the observed transition to scattered flow/density relation (the fundamental diagram) with hysteresis and self-sustained stop-and-go phenomena (see [START_REF] Treiterer | The hysteresis phenomenon in traffic flow[END_REF][START_REF] Kerner | Experimental properties of phase transitions in traffic flow[END_REF]9] and Fig. 1). This is due to the fact that spatially homogeneous regime are always in the equilibrium solutions and determined by the functions W (•) and V (•), respectively. Therefore, the microscopic behavior is modified by introducing reaction and relaxation times.

The simplest following model of this type may be the delayed model by Newell [START_REF] Newell | Nonlinear effects in the dynamics of car-following[END_REF] ẋi (t + τ ) = W (∆x i (t)),

with τ the reaction time (if positive). Applying a Taylor expansion in the l.h.s. of the delayed speed model Eq. ( 3), we obtain the second order 'optimal velocity model' (OVM) introduced by Bando et al. in [4]. The OVM has limit-cycles in stationary states, with self-sustained propagation of non-linear stop-and-go waves, and hysteresis curves in the fundamental flow/density diagram (see [34,[START_REF] Orosz | Traffic jams: dynamics and control[END_REF]). Macroscopic second-order models comprised of systems of hyperbolic equations are also able to reproduce non-linear stop-and-go waves and scattering of the fundamental diagram. One of the first approach is the one by Payne and Whitham [37,[START_REF] Whitham | Linear and Nonlinear Waves[END_REF]. The model can be derived from the microscopic Newell model Eq. ( 3). The main drawback of this model is that, as pointed out by Daganzo [START_REF] Daganzo | Requiem for second-order fluid approximations of traffic flow[END_REF], the speed and the density could yield negative values and are not bounded. Note that this drawback is also observed with follow-the-leader models like the OVM and is referred as collision between the vehicles (see for instance [START_REF] Davis | Modifications of the optimal velocity traffic model to include delay due to driver reaction time[END_REF][START_REF] Orosz | Exciting traffic jams: Nonlinear phenomena behind traffic jam formation on highways[END_REF] or [43, Chap. 15]). Aw and Rascle have corrected this issue by replacing the space derivative of the 'pressure' by a convective derivative [START_REF] Aw | Resurrection of "second order" models of traffic flow[END_REF] (AR model). Nowadays extensions of the AR model such as the ARZ, GARZ or generalized models [21, [START_REF] Zhang | A non-equilibrium traffic model devoid of gas-like behavior[END_REF]6,[START_REF] Flynn | Self-sustained nonlinear waves in traffic flow[END_REF]41,16], as well as two phase models coupled with the LWR model [11,[START_REF] Goatin | The Aw-Rascle vehicular traffic flow model with phase transitions[END_REF]12,[START_REF] Blandin | A general phase transition model for vehicular traffic[END_REF], are used to describe transition to congested traffic with scattered fundamental diagrams and self-sustained non-linear shock waves. A general framework is the generic second order model (GSOM) family introduced in [START_REF] Lebacque | Generic second order traffic flow modelling[END_REF][START_REF] Lebacque | Second order traffic flow modeling: supply-demand analysis of the inhomogeneous Riemann problem and of boundary conditions[END_REF]. Most of the approaches are a posteriori based on the continuous description.

In this article, we derive minimalist macroscopic traffic flow models of first order from a microscopic speed model to describe stop-and-go wave phenomena and scattering of the fundamental diagram. The use of first order models allow us to ensure by construction that the speed and the density remain positive and bounded. The starting point is a OV microscopic model of first order including a reaction time parameter. We show in Sec. 2 that the corresponding macroscopic model results in a convection-diffusion equation. The macroscopic model is discretized using distinct Godunov and Euler-based schemes and the linear stability conditions for the homogeneous solutions of these numerical schemes are provided in Sec. 3. The conditions match the ones of the car-following model for specific values of the spatial discretization step. Simulations are carried out in Sec. 4. Systems with different initial conditions are numerically solved. Further, we compare with data of realistic traffic flow as well as pedestrian flow.

2 Microscopic and macroscopic models

The microscopic follow-the-leader model

The microscopic model we use has been introduced in [START_REF] Tordeux | Collision-free nonuniform dynamics within continuous optimal velocity models[END_REF]. It is based on the Newell model Eq. ( 3). In the remaining of the paper, we assume that W : s → W (s) is Lipschitz continuous, non-decreasing and upper bounded in order to get the well-posedness of Eq. ( 3) supplemented with initial conditions x i (t = 0) = x i,0 for any i ∈ Z. We rewrite the equation as

ẋi (t) = W (x i+1 (t -τ ) -x i (t -τ )), (4)
and apply a Taylor expansion in the argument of W . Neglecting higher-order terms in τ we obtain ẋi (t) = W ∆x i (t)τ W (∆x i+1 (t)) -W (∆x i (t)) .

(

) 5
The model is a system of ordinary differential equations of first order with two predecessors in interaction. It is calibrated by the delayed time τ ∈ R, that is a reaction time if positive and an anticipation time if negative, and the optimal speed function W (•). The function W (•) is supposed to be bounded by V 0 > 0, positive and zero if the spacing is smaller than , > 0 being the vehicle's length or size of the pedestrian. Note that the model admits a minimum principle, say ∆x i (t) ≥ ∀i, t. Thus it is by construction collision-free and it has the same stability condition as the initial microscopic Newell model Eq. ( 3) or as the OVM from [4]. This condition is for all

s ∈ R |τ |W (s) < 1/2. ( 6
)
Note that the condition simply reduces to |τ | < T /2 if one considers the linear fundamental diagram W : s → W (s) := max 0, min 1 T (s -), V 0 , with T > 0. When unstable, the model transits to states with collision-free self-sustained stop-and-go dynamics, see [START_REF] Tordeux | Collision-free nonuniform dynamics within continuous optimal velocity models[END_REF].

Derivation of macroscopic models

In the following, we consider i = 1, . . . , N agents with periodic boundary conditions (i.e. the predecessor of the agent N is the agent 1). The derivation of macroscopic models from microscopic models is useful to fully understand the dynamics. In [2], Aw et al. established the connection between a microscopic car-following model and the second-order AR macroscopic traffic flow model. The rigorous proof, based on a scaling limit where the time and space linearly increase while the speed and the density remain constant, assumes homogeneous conditions. We use here the same methodology considering the local density ρ i (t) around the vehicle or pedestrian (i) and at time t > 0, as the inverse of the spacing

ρ i (t) := 1 ∆x i (t) . (7)
The density could also be normalized by multiplication with . Here, we prefer to keep the unit of one over length as density to ease the comparison with the classical models. Then, the microscopic model reads

ẋi (t) = W 1 ρ i (t) -τ W 1 ρ i+1 (t) -W 1 ρ i (t) =: Ṽ (ρ i+1 (t), ρ i (t)), (8)
for a velocity profile Ṽ . Then,

∂ t 1 ρ i (t) = ∂ t ∆x i (t) = Ṽ (ρ i+2 (t), ρ i+1 (t)) -Ṽ (ρ i+1 (t), ρ i (t)). (9)
In [2] it has been observed that Eq. ( 9) is a semi-discretized version of hyperbolic partial differential equation in Lagrangian coordinates. This requires to consider limits of many vehicles or pedestrians N → ∞ and diminishing length → 0. We introduce the continuous variable y ∈ R such that y i = i∆y as counting variable for the number of agents where ∆y is proportional to . By piecewise constant extension of the given spacing, we construct a density ρ(t, y) such that 1 ρ i (t) = 1 ∆y

y i + ∆y 2 y i -∆y 2 1 ρ(t, z)
dz. The quantity 1 ρ i (t) is the volume average over a cell of length ∆y centered at y i . The r.h.s. of (9) describes the flux across the cell boundaries. Introduce

V : k → V (k) = W ( 1 k ) for any k > 0. Then, it follows that Ṽ (k 1 , k 2 ) = V k 2 1-k 2 τ [V (k 1 )-V (k 2 )] for any (k 1 , k 2 ) ∈ (0, +∞) 2 satisfying V (k 1 ) = V (k 2 ) + 1 τ k 2 .
As for OVM, W is non-decreasing, therefore we observe that V is non-increasing on (0, +∞). We obtain from (9)

∂ t 1 ρ i (t) - ∆y ∆y V ρ i+1 1 -τ ρ i+1 Z i+1 -V ρ i 1 -τ ρ i Z i = 0, (10)
where

Z i := V (ρ i+1 )-V (ρ i ).
Provided that -V is increasing and ρ(t, y) is piecewise constant on each cell, the reconstruction at the cell interface y i± 1 2 is given by cell averages, i.e.,

1 ρ t,y i+ 1 2 = 1 ρ i+1 (t)
. Next, we rescale time t → t∆y and also reaction time τ → τ ∆y to obtain

∂ t 1 ρ i (t) - 1 ∆y V ρ i+1 1 -τ ρ i+1 Z i+1 ∆y -V ρ i 1 -τ ρ i Z i ∆y = 0. ( 11
)
Hence, we observe that in the rescaled time and in the limit ∆y → 0 the microscopic model is an upwind discretization of the following macroscopic equation

∂ t 1 ρ -∂ y V ρ 1 -τ ρ∂ y V (ρ) = 0. ( 12
)
The upwind or Godunov scheme is the most mathematically reasonable discretization provided τ is sufficiently small due to the decreasing behavior of V for suitable OVM functions W. Up to second-order in τ we approximate (12) by Taylor expansion and obtain a convection-diffusion model as

∂ t 1 ρ -∂ y V (ρ) = τ ∂ y (ρV (ρ)) 2 ∂ y ρ . (13)
The relation between the density in Lagrangian coordinates and Eulerian coordinates is given by the coordinate transformation (t, y) → (t, x) where y =

x -∞ ρ(t, x)dx. Note that y counts the number of vehicles/pedestrians up to position x in Eulerian coordinates. In the Eulerian coordinates (t, x), the macroscopic model Eq. ( 12) reads

∂ t ρ + ∂ x ρV ρ 1 -τ ∂ x V (ρ) = 0. ( 14
)
The model could be seen as an extension of the LWR model Eq. ( 2) with a modified speed-

density relationship ρ → V (ρ/(1 -τ ∂ x V (ρ))).
For illustrating the behavior of this modified speed-density mapping, we set 14) with constant inhomogeneity I ∈ {-0.3, 0, 0.3} and V : ρ → max{0, min{2, 1/ρ -1}}.

I := τ ∂ x V (ρ))
3 Linear stability analysis

Linear stability analysis of the continuous macroscopic model

A Taylor expansion up to second-order in terms of τ for equation ( 14) yields

∂ t ρ + ∂ x (ρV (ρ)) = -τ ∂ x (ρV (ρ)) 2 ∂ x ρ . ( 15
)
We also consider an initial condition

ρ(t = 0, x) = ρ 0 (x) for any x ∈ R (16)
where

ρ 0 ∈ L 1 (R)∩BV(R). By defining D(ρ) := x -∞
-τ (ρV (ρ)) 2 ∂ x ρ dy, we obtain an equation similar to the one considered in [8], say

∂ t ρ + ∂ x (ρV (ρ)) = ∂ 2 x D(ρ). (17)
One can verify that D(ρ 0 ) is absolutely continuous on R and that ∂ x D(ρ 0 ) ∈ BV(R). In (15), the l.h.s. is the LWR model with additional diffusion proportional to the reaction time parameter τ and that can be either negative or positive. More precisely the diffusion is negative in deceleration phases where the density get higher upstream, and it is positive in the opposite acceleration phases. This type of diffusion seems to induce an instability of the homogeneous (constant) solutions and the formation of oscillations (i.e. jam waves). The diffusion coefficient (ρV (ρ)) 2 depends on the density and the fundamental diagram. In fluid dynamics the coefficient is a characteristic for the flexibility of the random movement responsible for the diffusion. In traffic flows, comparable diffusion-convection forms have been used in [30,8]. We refer the interested reader to [8] (and references therein) for a proof of existence and uniqueness of the solution to Eq. ( 15)-( 16). In the following we analyze the linear stability of homogeneous solutions for the macroscopic model.

Proposition 3.1 The homogeneous configurations for which ρ(x, t) = ρ e for all x and t are linearly stable for the continuous traffic model Eq. ( 15) if and only if

τ < 0. ( 18
)
Note that a negative τ refers to an anticipation time.

Proof If ε(x, t) = ρ(x, t)ρ e is a perturbation to homogeneous solution ρ e , then

ε t = F (ρ e + ε, ε x , ε xx ) = αε + βε x + γε xx + o(LC(ε, ε x , ε xx )), (19)
with

F (ρ, ρ x , ρ xx ) = -∂ x ρV (ρ) -τ (ρV (ρ)) 2 ∂ x ρ , α = ∂F ∂ρ (ρ e , ρ e , ρ e ) = 0, β = ∂F ∂ρx = -V (ρ e ) - ρ e V (ρ e ), γ = ∂F ∂ρxx = -τ (ρ e V (ρ e )) 2
. The solutions of the linear system are the Ansatz ε = ze λt-ixl where λ ∈ C, x, l ∈ R. We get ε t = λε, ε x = -ilε and ε xx = -l 2 ε. Therefore the characteristic equation of the perturbed system Eq. ( 19) is λ l = τ (lρ e V (ρ e )) 2 + il(V (ρ e ) + ρ e V (ρ e )). The homogeneous solution are stable when ε → 0, i.e. (λ l ) < 0 for all l > 0. This holds only if the diffusion is positive. This is τ < 0.

Therefore the macroscopic model is unstable as soon as the reaction time τ is positive which is the physically reasonable case. An explanation is that the Taylor expansion of the original model in terms of τ does lead to a perturbed equation with different properties. However, the discrete model does not have this stability requirement. Therefore, we show below that for suitable discretization of the model we recover stability.

Linear stability analysis for the discrete schemes

Discretizations of the macroscopic models Eqs. ( 12) and ( 13) in 'Lagrangian' coordinates give the initial microscopic model Eq. [START_REF] Bayen | Lax-Hopf based incorporation of internal boundary conditions into Hamilton-Jacobi equation. Part I: Theory[END_REF]. Our purpose in this section is the discretization of the macroscopic models Eqs. ( 14) and (15) in Eulerian coordinates. We denote dt and dx the time and spatial discretization steps and use the Godunov scheme [START_REF] Godunov | A difference scheme for the numerical computation of a discontinuous solution of the hydrodynamic equation[END_REF] for the discretization of the density

ρ i (t + dt) = ρ i (t) + dt dx f i-1 (t) -f i (t) (20)
where f i denotes the flow at cell boundary and has to be determined. For this aim, we introduce the demand and supply functions from the flow-density fundamental diagram f : ρ → f (ρ) := ρV (ρ) as first proposed in [START_REF] Daganzo | Requiem for second-order fluid approximations of traffic flow[END_REF]26] and that read respectively

∆(ρ) := max k≤ρ f (k) and Σ(ρ) := max k≥ρ f (k) ( 21
)
and we define the Godunov flux as G(x, y) := min{∆(x), Σ(y)}. We are now ready to propose three different strategies to compute the boundary flows f i . The first two methods discretize the linearized model Eq. (15) using a splitting scheme which treats separately the convection and the diffusion terms. The last scheme is a simple discretization of the exact macroscopic model Eq. ( 14).

1. The Godunov/Euler scheme: a Godunov scheme for the convection term and an explicit Euler scheme for the diffusive term of the linearised model Eq. (15):

f (1) i = G(ρ i , ρ i+1 ) + τ dx (ρ i V (ρ i )) 2 (ρ i+1 -ρ i ). ( 22
)
Such a scheme is the one used in [2].

2. The Godunov/Godunov scheme: a Godunov scheme for the convection term and a Godunov scheme for the diffusion term of the Taylor-expanded model Eq. (15):

f (2) i = G(ρ i , ρ i+1 ) + τ dx ρ i V (ρ i ) G(ρ i+1 , ρ i+2 ) -G(ρ i , ρ i+1 ) . (23)
3. The Godunov scheme: a Godunov scheme for the modified convection term in the exact macroscopic model Eq. ( 14):

f (3) i = G ρ i 1 -τ dx V (ρ i+1 ) -V (ρ i ) , ρ i+1 1 -τ dx V (ρ i+2 ) -V (ρ i+1 ) . ( 24
)
Note that this scheme is valid if 1 -τ dx V (ρ i+1 ) -V (ρ i ) > 0 for all ρ i and ρ i+1 . By denoting V 0 = sup x V (x), this inequality holds if

τ < dx /V 0 . ( 25
)
Proposition 3.2 The homogeneous configurations for which ρ i (t) = ρ e for all i and t are linearly stable for the discrete traffic model Eq. 20 if and only if

α 2 + β 2 + γ 2 + ξ 2 -2αγ -2βξ + 2f (c l ) < 1, ∀l = 1, . . . , N -1, ( 26
)
with c l = cos(2πl/N ) and f (x) = (αβ + αξ + βξ -3γξ)x + 2(αγ + βξ)x 2 + 4γξx 3 , and α = ∂F ∂ρ i , β = ∂F ∂ρ i+1 , γ = ∂F ∂ρ i+2 and ξ = ∂F ∂ρ i-1 the partial derivatives of the model in equilibrium.

Proof The perturbations to homogeneous solution are the variables ε i (t) = ρ i (t)ρ e . The perturbed system is

ε i (t + dt) = ρ i (t + dt) -ρ e = F (ρ i (t), ρ i+1 (t), ρ i+2 (t), ρ i-1 (t)) -ρ e = α ε i (t) + β ε i+1 (t) + γ ε i+2 (t) + ξ ε i-1 (t) + o(LC(ε i , ε i-1 , ε i+1 , ε i+2 )), (27)
with α = ∂F ∂ρ i , β = ∂F ∂ρ i+1 , γ = ∂F ∂ρ i+2 and ξ = ∂F ∂ρ i-1 at (ρ e , ρ e , ρ e , ρ e ). General conditions for the global stability of the discrete schemes can be obtained for a system of N cells with periodic boundary conditions. The linear perturbed system is ε (t+dt) = M ε (t), with ε = T (ε 1 , . . . , ε N ) and M a sparse matrix with (ξ, α, β, γ) on the diagonal. If M = P DP -1 with D a diagonal matrix, then ε (t) = P D t/ dt P -1 ε (0) → 0 if all the coefficients of D are less than one excepted one equal to 1. M is circulant therefore the eigenvectors of M are z(ι 0 , ι 1 , . . . , ι m-1 ) with ι = exp i 2πl N and z ∈ Z, and the eigenvalues are λ l = α + βι l + γι 2 l + ξι -1 l . The system is linearly stable if |λ l | < 1 for all l = 1, . . . , N -1. This is

λ 2 l = α 2 + β 2 + γ 2 + ξ 2 -2αγ -2βξ + 2f (c l ) < 1, ∀l = 1, . . . , N -1, (28)
with c l = cos(2πl/N ) and f (x) = (αβ + αξ + βξ -3γξ)x + 2(αγ + βξ)x 2 + 4γξx 3 .

These conditions Eq. ( 26) are applied to the different numerical schemes Eq. ( 22), Eq. ( 23) and Eq. ( 24) with optimal speed V (ρ) = 1 T (1/ρ -), T > 0 being the vehicles time gap and ≥ 0 the vehicle' size. For such speed function, the Godunov scheme is simply G(x, y) = 1

T (1y ).

Lemma 3.3 The homogeneous configurations are linearly stable for the Godunov-Euler scheme Eqs. [START_REF] Godunov | A difference scheme for the numerical computation of a discontinuous solution of the hydrodynamic equation[END_REF](21)(22)

if 2τ < T dx ρ 2 e , (29)
and if dt is sufficiently small.

Proof Mixed with the scheme for the density Eq. ( 20), Godunov/Euler scheme Eq. ( 22) is

F 1 (ρ i , ρ i+1 , ρ i+2 , ρ i-1 ) = ρ i + dt T dx (ρ i+1 -ρ i ) + τ T dx ρ i -ρ i-1 ρ 2 i-1 - ρ i+1 -ρ i ρ 2 i , (30)
where

α = 1 -A + 2B, β = A -B, γ = 0 and ξ = -B with A = dt T dx and B = dt τ (T dx ρe) 2 . If τ < 1 2 T dx ρ 2 e , then α > 0 for dt < T dx + 2τ T dx ρ 2 e , (31)
and for all dt ≥ 0 if τ ≥ 1 2 T dx ρ 2 e . Moreover 1α > 0 if τ < 1 2 T dx ρ 2 e while β is positive only if τ < T dx ρ 2 e and the sign of ξ is the one of -τ . The stability conditions are distinguished according to the sign of τ .

• If τ < 0 and Eq. ( 31) holds, then [START_REF] Newell | Nonlinear effects in the dynamics of car-following[END_REF] holds and 1α > 0 on τ < 0. Therefore the system is stable for all τ < 0.

f (x) = α(1 -α)x + 2βξx 2 is strictly convex and is maximal on [-1, 1] for x = -1 or x = 1. Therefore the model is stable if f (-1) < f (1); this is simply -α(1 -α) < α(1 -α) that is always true since α > 0 if
• Several cases have to be distinguished for τ > 0. We assume in the following that Eq. ( 31) holds.

-For 0 < τ < 1 2 T dx ρ 2 e , we have α, 1α, β > 0, ξ < 0 and f (x) = α(1α)x + 2βξx 2 is strictly concave and maximal for

x 0 = -α(1-α) 4βξ > 0. The model is stable if x 0 > 1, this is dt < T dx - 2τ ( ρ e ) 2 . ( 32
)
This condition is more restrictive than Eq. ( 31).

-For 1 2 T dx ρ 2 e < τ < T dx ρ 2 e , we have α, β > 0, 1α, ξ < 0 and f (-1) > f (1) therefore the model is unstable. More precisely, f maximal for x 0 < -1, i.e. the unstable solution have shortest wavelength if

dt < τ - 1 2 T dx ρ 2 e 2τ T dx ρ e -2 . ( 33
)
This condition is also more restrictive than Eq. ( 31). -For τ > T dx ρ 2 e , we have α > 0, 1α, β, ξ < 0 and the system is unstable for all dt with shortest wavelength since f (•) strictly convex and f (-1) > f (1).

The stability conditions for the Godunov/Euler splitting scheme Eq. ( 22) are summarised in Fig. 3. The same conditions as the continuous macroscopic model are obtained for dx → 0 and dt → 0 such that dt / dx → 0.

Stable dt < T dx -2τ T dxρ 2 e Stable dt < T dx -2τ ( ρe) 2 Unstable dt < τ -1 2 T dxρ 2 e 2τ T dxρe 2 Shortest Wavelength

Unstable

For all dt > 0 Shortest Wavelength

and if dt is sufficiently small.

Proof The Godunov numerical schemes Eq. ( 23) and Eq. ( 24) are respectively

F 2 (ρ i , ρ i+1 , ρ i+2 , ρ i-1 ) = ρ i + dt T dx ρ i+1 -ρ i + τ T dx ρ i+1 -ρ i+2 ρ i - ρ i -ρ i+1 ρ i-1 , (35)
and

F 3 (ρ i , ρ i+1 , ρ i+2 , ρ i-1 ) = ρ i + dt T dx ρ i+1 1 -τ T dx 1 ρ i+2 -1 ρ i+2 - ρ i 1 -τ T dx 1 ρ i+1 -1 ρ i . ( 36
)
By construction, both give α = 1-A(1+B), β = A(1+2B), γ = -AB and ξ = 0 with A = dt T dx and B = τ T dx ρe . As expected, the stability conditions of these two schemes are the same. If τ > -T dx ρ e , then α > 0 for

dt < T dx + τ T dx ρe , (37)
and for all dt ≥ 0 if τ ≤ -T dx ρ e . β is positive only if τ > -1 2 T dx ρ e . Moreover 1β > 0 if Eq. (37) holds while the sign of γ is the one of -τ . Here again, the stability conditions are distinguished according to the sign of τ .

• If τ < 0 and Eq. (37) holds,

f (x) = β(1 -β)x + 2αγx 2 is strictly convex is maximal on [-1, 1] for x = -1 or x = 1. Therefore the model is stable if f (-1) < f (1); this is τ > - 1 2 T dx ρ e and dt < T dx + 2 τ T dx ρe . ( 38
)
The condition for dt is weaker than (37) since τ is negative. If τ ≤ -1 2 T dx ρ e then the system is unstable at the shortest wave-length frequency. A sufficiently condition for that the finite system produces the shortest frequency is simply N ≥ 2. Note that no condition holds on dt if τ ≤ -T dx ρ e .

• If τ > 0 and Eq. ( 37) holds then

f (x) = β(1 -β)x + 2αγx 2 is concave and is maximum at arg sup x f (x) = x 0 = -β(1-β) 4αγ > 0. We know that λ 2 0 = α 2 + β 2 + γ 2 -2αβ + f (1) = 1 (case l = 0). Therefore the model is stable if x 0 > 1; this is τ < 1 2 T dx ρ e and dt < T dx - 2τ ρ e . ( 39
)
The condition for dt is stronger than Eq. ( 37). If τ ≥ 1 2 T dx ρ e then the system is unstable at the frequency cos -1 (x 0 ) that is reachable in the finite system if N > 2π/ cos -1 (x 0 ). We have x 0 → 1/2 + T dx ρ e /(4τ ) as dt → 0, going from 1 to 1/2 according to τ (long wave).

The stability conditions for the Godunov/Godunov and Godunov schemes Eq. ( 23) and Eq. ( 24) are summarised in Fig. (4). The same conditions as the microscopic model are obtained at the limit dt → 0 for dx = 1/ρ e , i.e. a space step equal to the mean spacing. 23) and Eq. [START_REF] Kerner | Experimental properties of phase transitions in traffic flow[END_REF]. Note that we have the additional condition τ < dx /V 0 , with V 0 = sup x V (x), for the simple Godunov scheme Eq. ( 24).

on the speed and the density

The microscopic model Eq. ( 5) is collision-free: the spacing remains by construction bigger than the vehicle length > 0, and the speed is positive and bounded. We check whether this property also occurs with the numerical schemes F j , j = 1, 2, 3 of the macroscopic models Eqs. (30), [START_REF] Orosz | Traffic jams: dynamics and control[END_REF] and [START_REF] Orosz | Exciting traffic jams: Nonlinear phenomena behind traffic jam formation on highways[END_REF]. The models are from the first order therefore the speed is necessary positive and bounded if the optimal velocity functions are so defined. Moreover, the density remains bounded in [1, ρ M ], with ρ M = 1/ , if

F j (0, a, b, c) ≥ 0 and F j (ρ M , a, b, c) ≤ ρ M for all (a, b, c). (40)
It is easy to check that such property holds only if τ ≤ 0 for the Godunov/Euler Eq. (30), while it holds for τ ≥dx ρ e /W with the Godunov/Godunov scheme Eq. ( 35), and for τ < dx /V 0 with the simple Godunov scheme Eq. [START_REF] Orosz | Exciting traffic jams: Nonlinear phenomena behind traffic jam formation on highways[END_REF]. As the microscopic model, Eq. ( 35) and Eq. ( 36) are able to describe macroscopically unstable homogeneous solutions with large waves by ensuring that speed and density remain positive and bounded. The relation between instability and selfsustained traffic waves (or jamiton) are notably described in [START_REF] Flynn | Self-sustained nonlinear waves in traffic flow[END_REF]34,[START_REF] Orosz | Traffic jams: dynamics and control[END_REF]41] with microscopic and macroscopic second order models. In the next section, we analyse by simulation the unstable solutions we get with the first order models for different initial conditions.

Simulation results

In this section numerical simulations of the microscopic model Eq. ( 5) and of the simple Godunov scheme Eq. ( 36) macroscopic model are compared. The car-following model Eq. ( 5) is simulated using an explicit Euler scheme. A ring (periodic boundaries) with a length 101 and 50 vehicles is considered. The optimal speed functions are W (∆) = max{0, min{2, ∆ -1}} and V (ρ) = W (1/ρ) corresponding to a triangular fundamental diagram, while the reaction time is τ = 1. The values of the parameters are set to obtain unstable homogeneous solutions. The time step is dt = 0.01. The space step for the Godunov scheme is the mean spacing dx = 101/50 = 2.02 in order to match the stability conditions of both microscopic and macroscopic model (see Eq. ( 6) and Fig. 4) and to hold the CFL conditions (see Eq. ( 25) and Fig. 4). Three experiments are carried out with different initial conditions. In the first one, the initial configuration is a jam. The initial condition is random in the second experiment while it is a perturbed homogeneous configuration in the last one.

Trajectories

In Figs. 5, 6 and 7, the trajectories of the microscopic model and the time series for the density by cell for the discrete macroscopic model (gray levels) are plot for respectively the jam, random and perturbed initial conditions. The jam stationary propagates within the first experiment in Fig. 5. Both microscopic and macroscopic models rigorously describe the same dynamics. The dynamics obtained does not perfectly coincide for the random and perturbed initial conditions (see Figs. 6 and5). Yet most of the dynamics seems to be well recaptured and notably the self-sustained emergence of traffic stop-and-go waves. Note that the waves propagate backward with the speed -/T that is close to the value empirically observed (see [32]).

Fundamental diagram

The fundamental diagram is the plot of the flow or the mean speed as a function of the density. It generally refers to spatial performances [15], that have to be distinguished from temporal ones [START_REF] Wardrop | Some theoretical aspects of road traffic research[END_REF]. Here to deal with spatial performances, we measured the spatial speed and the density and express the flow as the product of the density by the speed. The density for the microscopic model is the inverse of the spacing (see Eq. ( 7)) while the speed in the macroscopic model is (see Eq. ( 9)) The sequences obtained for the perturbed initial conditions (see Fig. 7) are presented in Fig. 8. The performances are instantaneous ones in the sense that they correspond to instantaneous measurements for a vehicle (microscopic model) and a cell (macroscopic) in the system. The variability in such diagram is larger than the one of the aggregated fundamental diagram plotted in Fig. 1 where the performances were averaged over time intervals. Both microscopic and macroscopic systems converge to limit-cycles with self-sustained stopand-go waves resulting in hysteresis curves in the microscopic fundamental diagram. Such phenomenon generate scattering of the fundamental diagram for which some bounds can be calculated [START_REF] Zhang | A non-equilibrium traffic model devoid of gas-like behavior[END_REF]11,[START_REF] Goatin | The Aw-Rascle vehicular traffic flow model with phase transitions[END_REF]12,41,16]. The bounds V + and V -for the fundamental diagrams can here intuitively been determined from the microscopic model. The upper bound V + corresponds to the sequence of a vehicle moving at maximal speed V 0 behind a stopped vehicle:

Ṽ (ρ i , ρ i+1 ) = V ρ 1 -τ ρ(V (ρ i+1 ) -V (ρ i )) . (41)
V + (ρ) = Ṽ (ρ, 1/ ) = V ρ 1 + τ ρV (ρ) . (42)
Due to the reaction time, the distance tends to be smaller and the fundamental diagram is 'over-estimated'. Oppositely, the lower bound V -corresponds to the sequence of a stopped vehicle following a predecessor moving at the maximal speed V 0 :

V -(ρ) = Ṽ (ρ i , 0) = V ρ 1 -τ ρ(V 0 -V (ρ)) . ( 43
)
Here the reaction time induces a delay in the acceleration and an under-estimation of the fundamental diagram.

As in [41,16], the bounds Eqs. ( 42) and (43) obtained with the macroscopic model are compared to real instantaneous pedestrians and road traffic data in Figs. 9 and 10. The pedestrians data comes from a laboratory experiment with participants in a ring geometry [18]. Several experiments have been carried out with different density levels. The road traffic data are real measurement of trajectories on an American highway [START_REF]NGSIM: Next Generation Simulation[END_REF]. The speed, density and the flow are measured as previously (i.e. the density is the inverse of the spacing while the flow is the product of the density by the speed). A triangular fundamental diagram with 3 parameters V (ρ) = min V 0 , 1 T (1/ρ -1) is used again. The parameters are the ones of an estimation by least squares for the pedestrians V 0 = 0.9 m/s, = 0.3 m and T = τ = 1 s, see Fig. 9, while V 0 = 15 m/s, = 5 m and T = τ = 2 s for the vehicles, see Fig. 10. The bounds present a reasonable agreement with the data, even if no clustering of measurements are observed around them.

Conclusion

Starting from a speed following model, we derive a first order convection-diffusion continuum traffic flow model that we discretised using Godunov and Euler schemes. Simulation results shown that discrete macroscopic models can recapture the dynamics of the microscopic model, if specific values for the space discretization are chosen. More precisely, the linear stability conditions of the homogeneous solutions for the macroscopic models match the ones of the microscopic model for specific values of the space discretization and sufficiently small time steps. For unstable conditions, i.e. for large reaction times, the dynamics obtained describe selfsustained stop-and-go waves, with hysteresis cycles and a large scattering of the fundamental flow/density diagram. Such characteristics are observed in real data [START_REF] Treiterer | The hysteresis phenomenon in traffic flow[END_REF][START_REF] Kerner | Experimental properties of phase transitions in traffic flow[END_REF][START_REF] Kerner | Phase transitions in traffic flow[END_REF]9,[START_REF] Zhang | Universal flow-density relation of single-file bicycle, pedestrian and car motion[END_REF] as well as for second order models [START_REF] Zhang | A non-equilibrium traffic model devoid of gas-like behavior[END_REF]11,[START_REF] Goatin | The Aw-Rascle vehicular traffic flow model with phase transitions[END_REF]6,[START_REF] Flynn | Self-sustained nonlinear waves in traffic flow[END_REF]41,12,16]. Here it is achieved with first order models ensuring by construction that the models are physical and 'collision-free' (i.e. bounded and positive speed as well as density). Further investigations are necessary to understand the impact of the shape of the optimal velocity function on the characteristics of the waves. The macroscopic model corresponding to the follow-the-leader model is a first order elliptic convection-diffusion equation, for which the convection part is calibrated by the optimal velocity function (i.e. the fundamental diagram), while the diffusion is proportional to the reaction time parameter. More precisely the diffusion is negative in deceleration phases where the density get higher, and it is positive in acceleration phases where the density decreases. Such mechanism seems to be responsible for the appearance of oscillations and self-sustained non-linear stopand-go waves in the system. This observation remains to be confirmed rigorously, yet it could give us a way to explain the wave formations.

Figure 1 :

 1 Figure 1: Empirical fundamental diagrams. Left, [23, Figure 1] and right, [48, Figure 5].

Figure 2 :

 2 Figure 2: Illustration for the fundamental diagram V : (ρ, I) → V ρ/(1 -I) obtained in the macroscopic model (14) with constant inhomogeneity I ∈ {-0.3, 0, 0.3} and V : ρ → max{0, min{2, 1/ρ -1}}.

Figure 3 :

 3 Figure 3: Summary of the stability conditions for the Godunov/Euler splitting scheme Eq. (22).

Figure 4 :

 4 Figure4: Summary of the stability conditions for the Godunov/Godunov and Godunov schemes Eq. (23) and Eq.[START_REF] Kerner | Experimental properties of phase transitions in traffic flow[END_REF]. Note that we have the additional condition τ < dx /V 0 , with V 0 = sup x V (x), for the simple Godunov scheme Eq. (24).

Figure 5 :Figure 6 :

 56 Figure 5: The trajectories of the microscopic model (cyan curves) and the time series for the density by cell for the discrete macroscopic model (gray levels) for jam initial conditions.

Figure 7 :

 7 Figure 7: The trajectories of the microscopic model (cyan curves) and the time series for the density by cell for the discrete macroscopic model (gray levels) for perturbed initial conditions.

Figure 8 :

 8 Figure 8: Sequence of speed and flow / density relation for the perturbed initial conditions (see Fig. 7). Left, for one vehicle (microscopic model) and right, for one cell (macroscopic model).

Figure 9 :Figure 10 :

 910 Figure9: Instantaneous speed/density and flow/density measurements for real pedestrian flows[18] and the bounds Eqs. (42) and (43) for V 0 = 0.9 m/s, = 0.3 m and T = τ = 1 s.

1 3 3 3 . 1

 1331 The microscopic follow-the-leader model . . . . . . . . . . . . . . . . . . . . . . . 3 2.2 Derivation of macroscopic models . . . . . . . . . . . . . . . . . . . . . . . . . . . Linear stability analysis 5 Linear stability analysis of the continuous macroscopic model . . . . . . . . . . . 5 3.2 Linear stability analysis for the discrete schemes . . . . . . . . . . . . . . . . . . 6 3.3 Bounds on the speed and the density . . . . . . . . . . . . . . . . . . . . . . . . . 10 4 Simulation results 11 4.1 Trajectories . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 4.2 Fundamental diagram . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 5 Conclusion 14


Marc-Paul Errera
email: marc.errera@onera.fr

Roch Roukoz

El Khoury

A numerical efficient approach for coupled aerothermal simulations in fluid-structure systems

  L'archive ouverte pluridisciplinaire

INTRODUCTION

Optimizing heat transfer between fluids and solids is a crucial issue in a wide range of engineering applications, especially in the aerospace industry in which heat transfer enhancement is considered as a key parameter for developing more efficient and effective thermal devices. For instance, the tendency toward higher turbine inlet temperatures is continuing and it continues to raise considerable problems for the design of turbine cooling circuits.

In many situations it is not adequate to consider heat transport in the fluid only, with some idealized boundary conditions at the fluid-solid interface. Instead, the fully coupled problem must be considered, including heat conduction in the solid body. Conjugate heat transfer (CHT) analysis addresses such coupled problems as it is able to study processes which involve variations of temperature within the solids and the fluid simultaneously, due to their mutual thermal interaction.

In partitioned procedures, systems are spatially decomposed and the solution is separately advanced in time over each partition [START_REF] Felippa | Staggered transient analysis procedures for coupled dynamic systems: formulation[END_REF][START_REF] Piperno | Partitioned procedures for the transient solution of coupled aeroelastic problems[END_REF]. This strategy allows the direct use of designed solvers for different fields but may suffer from stability issues and performance limitations.

As has often been stressed, the nature of instabilities derived from a 1D model can give insight into the potential instabilities in 2D/3D flows computations. As a result, the behavior of interface conditions in CHT is often studied using a normal mode analysis, for instance by Roe et al. [START_REF] Roe | Combined interface boundary method for coupled thermal simulations[END_REF], Henshaw & Chand [START_REF] Henshaw | A composite grid solver for conjugate heat transfer in fluid-structure systems[END_REF], Kazemi-Kamyab et al. [START_REF] Kazemi-Kamyab | Accuracy and stability analysis of a secondorder time-accurate loosely coupled partitioned algorithm for transient conjugate heat transfer problems[END_REF], Joshi et al. [START_REF] Joshi | Stability analysis of a partitioned fluid-structure thermal coupling algorithm[END_REF], even if other means of investigation exist [START_REF] Lindström | A stable and high-order conjugate heat transfer problem[END_REF]. This is 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. Recently, it was shown, using a 1D thermal problem, that a numerical transition can be identified [START_REF] Errera | Optimal solutions of numerical interface conditions in fluidstructure thermal analysis[END_REF]. This transition was derived from a normal mode stability analysis, and results in an unconditionally stable coupling procedure that provides high convergence rates. The goal of this paper is to present the implementation of stable and efficient numerical schemes at the fluid-structure interface for aerothermal simulations on the basis of this new concept. This paper is composed of two parts. In the first part, the theoretical study is presented. The numerical interface treatment is briefly described for a Dirichlet-Robin transmission procedure. The pros and cons of this specific scheme are then presented along with the appropriate numerical treatments that optimize it. In order to further enhance the numerical efficiency of the coupled approach, we put forward another complementary numerical alternative based on the physical properties of the local coupled problem. These general numerical concepts are then illustrated in the second part of the paper via a CHT 3D industrial test case with emphasis put on the stability and convergence properties of the coupling schemes.

COUPLING MODEL

The partitioned CHT procedure

The partitioned CHT strategy presented hereafter is motivated by the desire to obtain a fast and stable steady solution of the fluid-structure thermal problem. It takes advantage of the significant discrepancy that exists between the characteristic time scales of the phenomena, namely a fast transient process in the fluid and a very slow one in the structure. Therefore, since the physical time-dependent response is not sought, the steady fluid-structure problem is solved by iteratively coupling a transient fluid solver to a steady state solid solver to achieve temperature and flux convergence at the interface. This implies that only the final fluid-solid steady state result is relevant.

Partitioned techniques are very popular as they allow the direct use of a specialized solver for each subdomain, offering significant benefits in terms of efficiency and code reuse.

The basic conventional serial staggered (CSS) algorithm was used to sequentially execute the CFD and the conduction solvers. The generic cycle is described in Fig. 1 where Ω f is the fluid sub-domain and Ω s is the solid sub-domain. , + n n t t corresponds to the current coupling step at which the two solvers exchange their interface values. This procedure is repeated until a global steady-state solution is obtained.

Interface coupling prototype

The transient energy equation in the fluid domain is discretized using a finite volume method (FVM) while the steady state heat conduction equation within the solid is solved using a finite element method (FEM).

An implicit backward Euler method is employed for the time derivative in the fluid and the heat diffusion terms in both domains are computed using a second order central difference formulation (see details in [START_REF] Errera | Optimal solutions of numerical interface conditions in fluidstructure thermal analysis[END_REF] ).

If the convection term is assumed to be negligible in the first fluid cell, the CHT problem may be modelled by a simple one-dimensional thermally coupled problem [START_REF] Errera | Optimal solutions of numerical interface conditions in fluidstructure thermal analysis[END_REF]. The equations in both domains are then iteratively solved using the CSS algorithm, with interface data being exchanged at each coupling iteration.

A general Robin transmission condition, considered on the solid side, leads to the equation :

f f f s f s T Q T Q α α + = + ˆ ( 1
)
where Q is the heat flux, T the temperature and f α is a coupling coefficient, with (^) denoting the sought values. The same type of interface condition can be defined on the fluid side :

s s s f s f T Q T Q α α + = + ˆ (2)
Applying the right hand side of Eq.( 2) as a boundary condition on the fluid domain, may technically be challenging in most research and commercial CFD packages. To overcome this hurdle, the convective boundary condition used in the CFD solver may be expressed in terms of a reference temperature f T and a convective heat transfer coefficient ext h (W.m -2 .K -1 ). This is obtained easily from Eq. (2) :

) ( ˆf ref ext f T T h Q - = (3) with      + = = s s s ref s ext T Q T h α α (4) n ( 4
Ω ∂ → Ω ∂ (3) Transfer n s n f Ω ∂ → Ω ∂ +1

Stability analysis

The Godunov & Ryabenkii theory [START_REF] Godunov | The theory of difference schemes. An introduction[END_REF][START_REF] Gustafsson | The Godunov-Ryabenkii condition : the beginning of a new stability theory[END_REF] is used to analyze the stability of the coupled fluid-solid procedure because it includes the effects of boundary conditions on the numerical stability of the problem.

The first step consists in introducing the normal mode solution in the discrete model equations, composed of the transient energy equation in the fluid and the steady heat conduction in the solid with Robin conditions on both sides of the interface. The second step consists in considering the eigen-solutions of the discrete problem. After elementary transformations (detailed in [START_REF] Errera | Optimal solutions of numerical interface conditions in fluidstructure thermal analysis[END_REF][START_REF] Errera | Efficient coupling procedures in steady and unsteady thermal analysis[END_REF]), the temporal amplification factor ) (z g of the coupled discrete thermal problem can be expressed as follows

) )( ( ) )( ( ) )( ( ) ( ) ( f s s f s s f f f f s s f f s f K K K K K K K z g α α α α κ α α α α + + - - + + + + = (5)
where f K and s K denote the fluid and solid thermal conductances, respectively defined by size grid

ty conductivi fluid x k K f f f = ∆ = (6)
The same applies to the conductance of the solid domain, but this time it worth noting that as only a steady state is sought at every coupling step, the characteristic size to be considered is the width of the solid domain s Λ size solid

ty conductivi solid k K s s s = Λ = (7)
The temporal amplification factor ) (z g depends upon the complex function ) , (

f f f D z κ κ =
. It can be shown that this function, solution of a quadratic equation obtained from the interior scheme in the fluid domain, introduces the mesh Fourier number f D defined by

2 f f f x t a D ∆ ∆ = (8)
where f a is the fluid thermal diffusivity (m 2 /s), t ∆ is the time step. The mesh Fourier number introduced in Eq. ( 8) characterizes heat conduction in the layer mesh of the transient domain. The stability of the coupled 1D CHT model is thus characterised by the modulus of its temporal amplification factor g(z).

EFFICIENT CHT TREATMENTS

Optimal coefficient

A general Robin-Robin interface condition simultaneously involves two coefficients. It is thus useful to initially limit the interfacial scheme to a Robin condition on the solid side and a Dirichlet condition on the fluid side (prescribed temperature). This is the most commonly used configuration in literature, because it is simple to implement and can deal with most situations likely to arise. But in order to go further, we will not limit ourselves to this scheme, as we shall see later. In other words, against all expectations, the existence of a transition value for f α is highlighted. At this transition value, the shape of the curve changes, resulting in the lowest amplification factor. This value, denoted

) (opt f α , is given by ) 1 ( 2 ) ( f f opt f D K - = α (9)
where f D is a dimensionless Fourier number defined by As shown in Fig. 2, two diametrically opposite behaviors can be identified :

= f D -stable min f f α α ≥ ∀ 1000 = f D -unconditionally stable f f f f D D D D 2 1 1 + + + = (10) f D describes the interval [ [ ∞ + , 0 as f D ranges over [ [
( (

) ) (opt f f α α > 2
: the amplification factor is an increasing function of f α . It is also a branch of hyperbola that tends to one as f α tends to infinity.

At the intersection

) (opt f f α α =
, the two branches of the amplification factor join up. At this remarkable value, the amplification factor attains its absolute minimum.

Numerical Biot number and stability

The expression of the mathematical expression of the optimal coefficient was provided recently at ONERA [START_REF] Errera | Optimal solutions of numerical interface conditions in fluidstructure thermal analysis[END_REF]. But it is now fundamental to identify all the stability zones involved in the Dirichlet-Robin transmission procedure. To the best of our knowledge, this has never been presented before.

The stability condition 1 ) , , ( < ∞ = s f z g α α applied to Eq.( 5) leads to the following lower bound stability

2 ) 1 ( 2 min s f f f K D K - - = α (11)
which can also take the following form

[ ] 1 ) 1 ( 2 ) ( min - - = ∆ f s f D Bi K α (12)
In Eq. ( 12), we have introduced the dimensionless number Bi may be regarded as a non-conventional local Biot number involved in a transient CFD calculation. Let us recall that the conventional Biot number needs a heat transfer coefficient that can be defined only at steady state. The new dimensionless number introduced in Eq. ( 13), is defined at any coupling iteration of a CHT computation and results from a balance between the unsteady fluid and steady solid domain properties.

From Eq. ( 12), it can be seen that two zones can be considered : --------------- -----------------∞ c o n d i t i o

1 ) 1 ( ) ( ≤ - ∆ f D Bi ( 14
- ) (opt f α -
n g U N S T A B L E 1 opt g 1 1 ) 1 ( ) ( ≥ - ∆ f D Bi g 0 g opt g 1 1 ) 1 ( ) ( ≤ - ∆ f D Bi Table 2.

Stability zones

The second zone is Then, the CHT procedure is inherently prone to instability and a relaxation coupling coefficient must be used such that min f f α α > . This lower stability limit may be very large. For example in Fig. 2 The main results of this stability analysis are summarized in Tab. 2.

As we can see, the unconditional stability of the coupling procedure can always be ensured by increasing When this criterion is not satisfied, a certain amount of relaxation is needed and the optimal amount of relaxation is provided by the choice of the optimal coefficient (Eq. ( 9).

Practical consequences

The theoretical results proposed in this paper can have the widest practicable effect in reducing the CPU time and in ensuring unconditional stability. Indeed, substituting Eq. ( 9) back into Eq. ( 5), we obtain the following amplification factor 1 ) , (

) ( ) ( ) ( < + = opt f s opt f opt f K z g α α α (16)
As a result, the optimal coefficient always provides a stable procedure and the lowest amplification factor.

But special attention has to be given to the conditions ( 14) and [START_REF] Heselhaus | A hybrid coupling scheme and stability analysis for coupled solid/fluid turbine blade temperature calculations[END_REF]. If condition [START_REF] El Khoury | Development and validation of coupling methods for the modelling of thermal systems[END_REF] holds, this means that the "dynamic" thermal resistance of the fluid domain at the shared interface is greater than the resistance offered by the whole solid domain. A Dirichlet-Robin condition is therefore appropriate.

If condition [START_REF] Heselhaus | A hybrid coupling scheme and stability analysis for coupled solid/fluid turbine blade temperature calculations[END_REF] holds, and at extreme cases where

1 ) 1 ( ) ( >> - ∆ f D Bi or s opt f K >> ) (
α , imply that the solid thermal gradients are not negligible or that the thermal fluid conductance is large compared to that of the solid. Here, a Dirichlet condition imposed on the fluid does not provide the most efficient solution, even though the optimal coefficient is able to stabilize the procedure. A coupling procedure that is more suited for the physics and the numerical setup at hand would be to impose a Neumann condition on the fluid side.

Conventional and numerical Biot

In this paper a "numerical" local Biot number has been introduced and particular emphasis has been placed on its key role in the stability process. This new local Biot number takes into account the thermal and dynamic response of the boundary layer and directly participates in the stability of the coupled process. It is worth noting that diffusion, on both sides of the interface, dominates and guides the coupling process, as long as transients effects prevail.

On a different level, the conventional Biot Number is a criterion which gives a direct indication of the relative importance of conduction and convection. It measures how is the resistance to heat flow within the solid relative to the resistance presented by the convection processes at the surface. As a result, this number can be a key parameter, at steady state and only at steady state, to determine the stability of the fluid-solid equilibrium. But this parameter, not defined during the transients, cannot be used to set up a numerical CHT procedure, as long as a transient fluid state is involved in this procedure.

Efficient strategy at the interface

Realistic and accurate thermal conjugate heat transfer solutions require high-resolution CFD meshes. Consequently the placement of the first fluid grid point should fall in the near-wall boundary-layer region, in order to ensure y+ values close to unity. This is mandatory when low Reynolds RANS turbulence models are used, due to the sensitivity of their thermal wall fluxes and temperatures to the boundary layer resolution. In such cases, it is safe to neglect the effects of convection on the first fluid mesh element, which allows us to consider it as a small perturbation, that does not affect the nature of the diffusion equation.

Another important point concerns the canonical prototype proposed in this paper. The relevance of this predictive model problem, particularly the so-called "optimal coefficient" has been confirmed in many CHT problems, but it is clear that this model does not contain

1 ) 1 ( ) ( ≥ - ∆ f D Bi (15)
velocity component normal to the wall and thus, this approach is valid if the Péclet number which measures the ratio of convection to diffusion is small. In other words, a good representation of the flow structure is here necessary. This also implies that the grid must be refined enough so that convection may thus be considered as a small perturbation.

The two preceding considerations imply that the ratio

) (∆
Bi is likely to be very high as it directly depends on the first fluid cell size. This is becoming a tendency, with the ever growing computational resources dedicated to CFD computations.

It is therefore essential to develop an efficient strategy for coupled CHT calculations that are shifting towards very fine grids at the interface. In such configurations, the widely used Dirichlet condition on the fluid side, quickly show its limitations. Of course, the dynamic effect of the fluid domain, represented by

) 1 ( f D -
, is capable of stabilizing any CHT calculation, but the numerical performance of such coupling will be far from ideal, as it does not take into account the physics on either side of the interface.

For that matter, a physics-based coupling strategy is proposed below. It relies on the stability analysis of the model problem, while emphasizing on criteria ( 14) and ( 15) that can be locally computed, at every interface mesh cell: ____________________________________________ local computation of

) (opt f α in the fluid if ( 1 ) 1 ( ) ( ≤ - ∆ f D Bi ) then //Dirichlet-Robin procedure ∞ = s α //Dirichlet on the fluid ) (opt f f α α = //optimal Robin on the solid else // 1 ) 1 ( ) ( > - ∆ f D Bi if ( s opt f K >> ) ( α )

TEST CASE

First CHT tests

The formulation of the optimal Robin coefficient was published relatively recently [START_REF] Errera | Optimal solutions of numerical interface conditions in fluidstructure thermal analysis[END_REF], and consequently very few CHT computations that take advantage of it are found in literature. However, some interesting results have been achieved.

For instance, in a steady CHT computation of an effusion cooling system [START_REF] Errera | Temporal multiscale strategies for conjugate heat transfer problems[END_REF], the CPU time necessary to converge was divided by a factor of ten as compared to a conventional method.

In another work devoted to test systematically the values of various coupling coefficients, it has been shown that the optimal coefficient can be applied in an efficient manner as a tool for predicting and obtaining excellent stability properties [START_REF] Errera | Comparative study of coupling coefficients in Dirichlet-Robin procedure for fluid-structure aerothermal simulations[END_REF]. Indeed, the CHT results obtained by employing the optimal coefficient have always proved to be stable and oscillation-free. This would suggest that the one-dimensional normal mode analysis can provide relevant coefficients directly applicable to industrial 3D CHT problems.

We now present results to illustrate the importance of the interface treatment in CHT problems, with emphasis put on the optimal treatment and the strategy proposed in this paper. The details about this computation can be found in [START_REF] El Khoury | Development and validation of coupling methods for the modelling of thermal systems[END_REF].

3D industrial test case

The considered 3D industrial setup represents a cell of a reheat oven and puts into play various interface configurations, ranging from metallic parts to thick insulating components, simultaneously interacting with the fluid flow. This provides a realistic setup with a large number of heteroclite coupled interfaces. The heterogeneous properties of this set-up are ideal to study the nature of the coupling procedures. Radiative heat transfer is not considered in this work, because it is outside the scope of the current work.

The walls of the cell are made out of, thick, high temperature insulating material. Molybdenum electrical heaters are positioned above and below the steel element. Turbulent air, flows from the inlet and interacts with the walls of the chamber and its metallic elements before exiting through the outlet. The upper wall has an imposed external temperature, while the external faces of the lower and side walls are assigned a convection boundary condition. The near-wall fluid mesh density that was used in the CHT calculations is refined enough (y + ~1) to correctly capture the flow boundary layer. Four coupling coefficients have been used :

(0) ) ( 0 opt f f stnrd f K H α α α > = = =
(see [START_REF] Heselhaus | A hybrid coupling scheme and stability analysis for coupled solid/fluid turbine blade temperature calculations[END_REF]) ). This latter coefficient provides an over-stable solution. When the coefficient is equal (

) ( 1 opt f f α α α > = (2) ) ( 2 opt f f α α α = = (1)
) ( 2 opt f f α α α = =
) to the optimal value, the convergence is, not surprisingly, always monotonic (for all the Fourier numbers considered) and provides the best numerical behavior.

These convergence behaviors were predictable from the model problem. As expected, ) (opt f α appears to be the optimal choice, just at the intersection of the Neumann and the Dirichlet branches [START_REF] Errera | Comparative study of coupling coefficients in Dirichlet-Robin procedure for fluid-structure aerothermal simulations[END_REF]. This coefficient provides a monotonic, stable and fast convergence of the coupled simulations.

OPTIMIZED TRANSMISSION SCHEME

The preceding results clearly show that for a given transmission condition scheme (in this case Dirichlet-Robin), it is possible to optimize the numerical procedure in terms of stability and convergence. It is already a significant achievement. However, as stated before and explicitly described in Tab. 2, it is useful to know if the algorithm proposed in Tab. 3 can be used.

In other words, is the Dirichlet-Robin condition the most appropriate scheme ? Of course, we must recognize that the Dirichlet-Robin condition and the resulting optimal coefficient (Eq. ( 9) always gives a monotonic and stable behavior. But can we proceed a step further ?

Indeed, the computational domain in the CHT case contains may disparate zones and it is instructive to consider the values of the coupling coefficients computed by the model problem. In some of these zones, the coupling coefficient needed for stability is clearly very high. This means that stability is achieved at the expense of computational efficiency (CPU time) and obviously, a Dirichlet condition imposed on the fluid side is not the best choice. The problem arises when the numerical Biot number ) (∆ Bi is too high. In this case, a natural physics-based approach is to impose a Neumann condition on the fluid side, as identified in the algorithm of Tab. 3. The first tests with such an approach show that convergence occurs 3.2 times faster than the standard coupling scheme.

In the preceding test case, the ratio between the thermal conductances is approximately It was therefore of great importance to use the strategy defined in Section 3.5. Fig. 5 shows the residual of the Neumann-Dirichlet (N-D) procedure as well the convergence history of the optimal Dirichlet-Robin condition. These CHT computations highlight the importance of a relevant transmission procedure. Let us recall that the N-D procedure works with no relaxation parameters and shows the fastest convergence properties. Clearly, the strategy proposed in this paper (Tab. 3) seems to be the most physics-based approach, while always being immune to oscillations.

CONCLUSION

We have presented the stability analysis of the CHT problem, with emphasis put on the Dirichlet-Robin procedure. It has been shown that the Fourier and the numerical Biot numbers were the key parameters predicting the numerical behavior of the coupling process. They are used to obtain the optimal Dirichlet-Robin procedure that can be regarded as an adaptive procedure to always obtain the fastest rate of convergence and the best stability properties. This is certainly a first essential step and further improvement of the CHT computations were proven to be attainable. The local Biot number introduced in this paper can also be regarded as a good criterion to choose a physics-based approach. With this insight, the most appropriate interface transmission condition can be implemented as opposed to being limited to a single procedure. As a result, a strategy based on two complementary interface schemes has been presented.

This strategy has been applied to a 3D CHT test case exhibiting disparate coupling surfaces making it an excellent candidate to testing the new coupling strategy. This numerical example has clearly illustrated that a single interfacial approach can provide effective and practical solutions directly applicable to aerothermal coupled computations. But a higher numerical efficiency and substantial computing time savings could be attained when using the proposed adaptive algorithm.

Figure 1

 1 Figure 1. CSS Algorithm The interval [ ]

1

 1

  can observe that this function is defined and continuous and that each curve is composed of two half-lines with a singular point similar to a cusp at the intersection. At this point,

Figure 2

 2 Figure 2. ) ( max z g for 2 different Fourier numbers 10 =

:

  the amplification factor is a decreasing function of f α . It is a branch of hyperbola that depends directly upon the Fourier number f D . The greater f D , the smaller the amplification factor.

  Neumann transmission condition can be used, with no relaxation, without affecting the stability of the coupled problem. In practice however, a Dirichlet-Robin transmission is preferred in such cases, in order to avoid an ill posed thermal problem on the solid side.

  ,

  that Condition (14) holds.

Fig. 2

 2 Fig. 2 illustrates, for two different Fourier numbers ( . 10 = f D and

  then //Neumann-Robin procedure

Figure 3

 3 Figure3shows the streamlines and the temperature contours of the flow in a central cross section. The complex geometrical setup of the problem generates

Figure 3 .

 3 Figure 3. Temperature contours & flow streamlines

  As might be anticipated, the convergence history largely depends on the value of the coupling coefficient as shown in Fig.4. For low values of the coefficient 3 α α = f , regardless of the fluid Fourier number employed, the CHT computations are unstable. In other words, oscillations that grow without bound are obtained.

Figure 4 .

 4 Figure 4. Convergence vs coupling iterations (D-R) This tendency can be observed in Fig. 4 at the beginning of the computation before divergence. Larger values

Figure 5 .

 5 Figure 5. Convergence vs coupling iterations (D-R & N-D)

Table 1 .

 1 Normalized Fourier number

	1 0 . It is a rapidly increasing function as can be seen ,
	in Tab. 1.
	D	f	1.	2.	4.2 10. 100. 179. 1000.
	f D 0.27 0.38 0.50 0.64 0.87 0.90 0.95

The proposed strategy, involving two complementary interface treatments, should be extended to cover the most general coupling approach, based on the Robin-Robin transmission procedure, in order to cope with the increasing complexity of CHT simulations.