D Ameur
email: ameur@l3m.univ-mrs.fr

I Raspo

U M R M2p2

NUMERICAL SIMULATION OF THE POISEUILLE-RAYLEIGH-BÉNARD INSTABILITY FOR A SUPERCRITICAL FLUID IN A MINI-CHANNEL

The Poiseuille-Rayleigh-Bénard problem, involving the onset of thermoconvective structures in channels heated from below, was the subject of many theoretical, numerical and knowledge, this problem was never studied for supercritical fluids. The objective of this paper is to study the influence of the specific properties of such fluids on thermoconvective instability phenomena compared with those observed in the perfect gas case. The effect of the distance to the critical point is also investigated. The numerical approach used is based on the Navier-Stokes equations in the framework of the low Mach number approximation.

KEY WORDS: thermoconvective instability, convective instability, absolute instability, mixed convection flows, spectral method, low Mach number approximation Moreover, they are very sensitive to temperature and pressure variations. These tunable properties motivated the use of SCF in many industrial applications, such as supercritical fluid extraction [START_REF] Reverchon | Supercritical fluid extraction and fractionation of natural matter[END_REF], particle generation [START_REF] Cocero | Encapsulation and coprecipitation processes with supercritical fluids: fundamentals and applications[END_REF], in particular in pharmaceutical industry, and preparation of metallic films using supercritical fluid deposition [START_REF] Erkey | Preparation of metallic supported nanoparticles and films using supercritical fluid deposition[END_REF].

The instability onset and the development of the associated patterns in PRB flows were extensively studied for incompressible fluids and perfect gas from many years because the PRB configuration is relevant for several technological processes such as the cooling of microelectronic equipments [START_REF] Incropera | Convective heat transfer in electronic equipment cooling[END_REF] or the growth of thin crystal films from chemical vapor deposition (CVD) [START_REF] Evans | A study of traveling wave instabilities in a horizontal channel flow with applications to chemical vapor deposition[END_REF][START_REF] Evans | Unsteady three-dimensional mixed convection in a heated horizontal channel with applications to chemical vapor deposition[END_REF][START_REF] Evans | Thermally unstable convection with applications to chemical vapor deposition channel reactors[END_REF][START_REF] Spall | Unsteady mixed convection in horizontal ducts with applications to chemical vapor deposition processes[END_REF]. A detailed review of various studies reported in literature was performed by [START_REF] Nicolas | Bibliographical review on the Poiseuille Rayleigh Bénard flows: the mixed convection flows in horizontal rectangular ducts heated from below[END_REF]. These studies showed that when the base flow becomes unstable, two kinds of thermoconvective structures may appear: transversal rolls at low

Reynolds number (about Re<10) and longitudinal rolls at higher Reynolds number. The transversal rolls axes are perpendicular to the mean flow direction and these patterns can be considered as a quasi two-dimensional structure, whereas the longitudinal rolls axes are parallel to the mean flow and the three velocity components are excited. Results of linear stability analysis for incompressible flows showed that the transversal rolls are due to a convective or an absolute instability depending on the value of the Rayleigh number [START_REF] Müller | Transversal convection patterns in horizontal shear flow[END_REF][START_REF] Carrière | Convective versus absolute instability in mixed Rayleigh-Bénard-Poiseuille convection[END_REF]. When the flow is convectively unstable, the initial perturbation locally grows, it moves downstream and it finally leaves the system.

Hence, the rolls cross the channel during time and the flow becomes stable again once they moved away. On the other hand, in the case of absolute instability, the initial perturbation locally grows but it also expands in the whole system so that the rolls appear throughout the channel all the time although they move downstream: when a roll moves away, another one is locally generated in its place by the perturbation. The convective or absolute nature of the instability depends on the value of the Rayleigh number Ra, as is shown by the stability diagram in the plane Ra-Re (Fig. 1). The critical Rayleigh numbers for the onset of convective instability, Rac1, and for the onset of absolute instability, Rac2, are computed from the formulae of [START_REF] Müller | Transversal convection patterns in horizontal shear flow[END_REF] for a Prandtl number Pr=0.6567 (corresponding to the value for CO2 as a perfect gas). It can be seen that they both increase with the Reynolds number. On the other hand, [START_REF] Carrière | Convective versus absolute instability in mixed Rayleigh-Bénard-Poiseuille convection[END_REF] showed that the longitudinal rolls can only be convectively unstable and that the critical Rayleigh number is independent of Re and Ra. These theoretical predictions were confirmed by experiments.

m was never studied for fluids near their liquid-gas critical point while it is relevant for supercritical fluid deposition processes and it can be expected that the fluid response to the bottom boundary heating may be quite different in this special zone of the phase diagram. Indeed, as was already said, near the critical point, thermodynamic properties, such as density, and transport coefficients exhibit large variations with temperature and pressure. More precisely, on the critical isochore, several properties, such as isothermal compressibility, thermal expansion coefficient, specific heats or thermal conductivity, diverge as the critical temperature is approached, while other properties, such as thermal diffusivity, tend to zero. These critical behaviors cause some peculiar phenomena in supercritical fluids, such as a fast heat transport by thermoacoustic effects (also called Piston effect) [START_REF] Onuki | Fast adiabatic equilibration in a singlecomponent fluid near the liquid-vapor critical point[END_REF][START_REF] Boukari | Critical speeding up in pure fluids[END_REF][START_REF] Zappoli | Anomalous heat transport by the piston effect in supercritical fluids under zero gravity[END_REF]. Natural convection in SCF was extensively studied in the past decade [START_REF] Zappoli | Thermoacoustic and buoyancy-driven transport in a square side-heated cavity filled with a near-critical fluid[END_REF][START_REF] Chiwata | Thermal plumes and convection in highly compressible fluids[END_REF][START_REF] Amiroudine | Direct numerical simulation of instabilities in a two-dimensional near-critical fluid layer heated from below[END_REF][START_REF] Furukawa | Convective heat transport in compressible fluids[END_REF][START_REF] Furukawa | Convection in a very compressible fluid: comparison of simulations with experiments[END_REF][START_REF] Amiroudine | Piston-effect-induced thermal oscillations at the Rayleigh-Bénard threshold in supercritical He[END_REF][START_REF] Raspo | Unsteady two-dimensional convection in a bottom heated supercritical fluid[END_REF][START_REF] Accary | Three-dimensional Rayleigh-Bénard instability in a supercritical fluid[END_REF]Accary et al., 2005a;Accary et al., 2005b;[START_REF] Gorbunov | Conditions of Rayleigh-Bénard convection onset and heat transfer in a near-critical medium[END_REF], especially in the Rayleigh-Bénard configuration [START_REF] Amiroudine | Direct numerical simulation of instabilities in a two-dimensional near-critical fluid layer heated from below[END_REF][START_REF] Furukawa | Convective heat transport in compressible fluids[END_REF][START_REF] Furukawa | Convection in a very compressible fluid: comparison of simulations with experiments[END_REF][START_REF] Amiroudine | Piston-effect-induced thermal oscillations at the Rayleigh-Bénard threshold in supercritical He[END_REF][START_REF] Raspo | Unsteady two-dimensional convection in a bottom heated supercritical fluid[END_REF][START_REF] Accary | Three-dimensional Rayleigh-Bénard instability in a supercritical fluid[END_REF]Accary et al., 2005a;Accary et al., 2005b;[START_REF] Gorbunov | Conditions of Rayleigh-Bénard convection onset and heat transfer in a near-critical medium[END_REF]. These works showed, in particular, that the adiabatic temperature gradient, which is generally relevant only at large length scale, like for the study of atmospheric flows, must be taken into account for the convection onset in supercritical fluids in cavities of a few millimeters in height. We can therefore suppose that the specific properties of SCF may also modify the instability onset in the PRB configuration.

In the present paper, the stability of a supercritical fluid, modeled by the Peng-Robinson equation of state, in the PRB configuration is investigated using 2D direct numerical simulations. The first section is devoted to the mathematical formulation of the treated problem, namely the Navier-Stokes equations in the framework of the low Mach number approximation and the boundary conditions. Then, the numerical method, based on a Chebyshev collocation approximation, is described in details. Results obtained for a

Reynolds number equal to 50.57 in a channel of height H presented in the third section. First, we studied the influence of the temperature increase on the onset and the development of convection for a fluid set at 1K above its critical temperature. The relation between the Rayleigh number and the instability characteristics is thus investigated. In a second part, we are focusing on the influence of the proximity to the critical point on the instability onset.

PHYSICAL PROBLEM AND MODELLING

The physical model consists of a 2D channel of height =1mm, with an aspect ratio =10 or 15. The carried fluid is supercritical CO2 subject to gravity. Initially, the fluid is at a uniform temperature slightly above the critical temperature, Ti (1+ )Tc (with <<1), and at a mean density equal to the critical density c , and it flows through the channel according to a Poiseuille profile. Then, the temperature of the bottom wall is gradually increased up to Tch =Ti + (with ranging from about ten to some hundreds milliKelvin) from a distance from inlet (see Fig. 2).

The supercritical fluid is modeled by the Peng-Robinson equation of state. This equation implicitly accounts for the divergence of the thermal expansion coefficient , of the isothermal compressibility and of the specific heat at constant pressure CP near the liquidgas critical point. The divergence of the thermal conductivity is modelled by the formula

0.5 [1 1 ] b c T T
. The physical parameters of CO2 were used:

Tc =304.13K, c 467.8kg.m -3 , b 0.04412W.m -1 .K -1 , =0.75.
The evolution of the flow is governed by the time-dependent 2D Navier-Stokes equations coupled with the energy and the Peng-Robinson equations. These equations are solved in the framework of the low Mach number approximation [START_REF] Paolucci | On the filtering of sound from the Navier-Stokes equations[END_REF]: the pressure is thus split into a thermodynamic part, Pth , which is constant in space and appears in the energy equation and in the equation of state, and a dynamic part, Pdyn , involved in the momentum equation. However, the basic approximation of [START_REF] Paolucci | On the filtering of sound from the Navier-Stokes equations[END_REF] is modified as proposed in Accary et al. (2005c) to account for the stratification of the fluid near the critical point since for time and b for thermal conductivity.

The specific heat at constant volume CV and the dynamic viscosity were fixed to their background values, CVb 472.313J.kg -1 .K -1 and b =3.2702×10 -5 Pa.s. The dimensionless governing equations are therefore:

. 0 t V (1) 1 1 . . 3 dyn i y Pr P t Ra Fr V V V V V e
(2)

0 0 0 . 1 . . T P Pr T T T t T Pr Ra V V (3) 2 2 2 1 1 2 th hyd a T T P P b b b (4)
with a and b the dimensionless energy parameter and covolume respectively and the Soave function defined by:

2 1.487422 , 0.253076 , 1 1 a b T m T
where m is computed from the acentric factor ( =0.225 for CO2) by the formula 2 0.37464 1.54226 0.26992 m .

In Eqs. ( 1)-( 4), Phyd is the hydrostatic pressure introduced by the modification of the Low Mach number approximation (Accary et al., 2005c), ey is the unit vector in the y-direction and 0 and Pr0 are respectively the specific heats ratio and the Prandtl number for the perfect gas ( 0=1.4, Pr0=0.6567). The other dimensionless numbers introduced are the Prandtl number

Pr, the Rayleigh number Ra, the Mach number Ma and the Froud number Fr which are defined by:

2 2 3 0 , = , , ref ref P b c P b b b c U U C g C T H Pr Ra Ma Fr g H R T .
In the above formula, the physical parameters and CP are calculated for the initial state (Ti i from the equation of state. The initial condition for the dimensionless variables in

=[0, ]×[0,1] is: 2 , 1.5 1 2 1 , , 0 i i Pr u x y Re y v x y Ra ( 5
)
where Re is the Reynolds number defined by , 1

K y i K e x y K e (7) , i thi hyd P x y P P y (8) with 2 2 0 1 2 2 2 1 2 1 1 , 1 1 2 i a b T Ma K K K Fr b b b and: 2 2 1 2 2 1 , 1 1 1 2 1 K y i thi hyd K a T e P P y K K b b b e
Equations ( 1)-( 4) are solved with the following boundary conditions:

-On the channel walls, the no-slip condition is prescribed for the velocity. The top wall ( y=1) is kept at the initial temperature Ti. On the bottom heated wall (y=0), in order to avoid a discontinuity of the temperature profile, the following boundary condition is imposed for

0 x : , 0, 2 2 2 18 2 i T x t T T th x th th th (9)
This boundary condition allows a continuous transition between the cold entry zone for 0

x 1 and the hot zone for 1 x .

-At the channel inlet (x=0), temperature is kept at its initial value Ti and a parabolic profile is imposed for velocity:

2 0, , 1 0,

NUMERICAL METHOD

Time scheme and space approximation

Equations ( 1)-( 4) are discretized in time with a second order semi-implicit scheme: the convective terms are evaluated by an Adams-Bashforth scheme, a second order backward Euler scheme is used for the discretization of time derivatives and the diffusive terms as the energy source term . T P T

V are implicitly treated. The outlet boundary condition ( 11) is also discretized using the same second order scheme, leading to a Dirichlet boundary condition for the variables at the current time step n+1.

The space approximation is performed using a Chebyshev-collocation method with Gauss-Lobatto points. For the computation of the convective terms, the derivatives are calculated in the spectral space and the products are performed in the physical one; the connection between the spectral and the physical spaces is realized through a FFT algorithm. On the other hand, the spectral differentiation matrices are used for the derivatives in the diffusive terms.

Numerical algorithm

At a given time step, the discretized equations are coupled because of the implicit treatment of the energy source term involving .V. This term must be implicitly evaluated for supercritical fluids. However, it is possible to completely uncouple the solution of the energy equation and the computation of the dynamic field by calculating the velocity divergence from the sole knowledge of the thermodynamic variables [START_REF] Ouazzani | A new numerical algorithm for low Mach number supercritical fluids[END_REF][START_REF] Raspo | Un algorithme faible nombre de Mach pour la simulation des écoulements de fluides supercritiques par des méthodes spectrales[END_REF]. The procedure consists in taking the total derivative of the equation of state (4) written in the form F(Pth )=0. Then, using the energy equation ( 3) and the continuity equation ( 1), the following formula is obtained for the velocity divergence: In the case of an open system, as the one considered here, the thermodynamic pressure Pth is also constant in time. Therefore, the time derivative dPth/dt disappears in Eq. ( 12). Otherwise, this derivative can be evaluated using a second order backward Euler scheme as in [START_REF] Raspo | Un algorithme faible nombre de Mach pour la simulation des écoulements de fluides supercritiques par des méthodes spectrales[END_REF].

Thanks to Eq. ( 12), the discretized equations are solved in two successive steps: first, the thermodynamic variables are computed through an iterative algorithm [START_REF] Ouazzani | A new numerical algorithm for low Mach number supercritical fluids[END_REF][START_REF] Raspo | Un algorithme faible nombre de Mach pour la simulation des écoulements de fluides supercritiques par des méthodes spectrales[END_REF] and then the Navier-Stokes equations are solved using the modified projection method developed in [START_REF] Hugues | An improved projection scheme applied to pseudospectral methods for the incompressible Navier-Stokes equations[END_REF] and extended to variable density flows. The two steps are described in detail below.

-1 st step: Computation of ( .V) At time step n+1, the discretized energy equation can be written as a Helmholtz equation with coefficients involving n+1 and n+1 . In order to take advantage of the efficiency of the full diagonalization technique developed for the solution of Helmholtz problems with time-independent coefficients [START_REF] Haldenwang | Chebyshev 3D spectral and 2D pseudospectral solvers for the Helmholtz equation[END_REF], the density and the thermal conductivity are split into a constant part, namely their value at t=0, and a time-dependent part. So, the terms involving the time-dependent parts are treated as source terms and the discretized energy equation can be written as a Helmholtz equation with constant coefficients. The complete diagonalization of the Helmholtz operator with constant coefficients is performed once for all in a preprocessing stage. Then, at each time step, the solution of the Helmholtz problem is reduced to matrix products, leading to a very efficient solution technique on supercomputers.

More precisely, the computation of T, and .V is performed through the following iterative algorithm:

1. The variables T k-1 , k-1 and ( .V) k-1 are initialized at their values at the previous time step n;

2. The temperature T k is obtained by the solution of the Helmholtz equation:

1 1 1 1 0 0 1 1 1 1 1 0 , 1 3 3 . 2 2 4 1 . 2 . k k k k k k i i i i 0 0 k n n k k k th n n Pr T Pr T T T Pr Ra t t Pr Ra P T T T V T t AB T V
with the boundary conditions of the problem. In the above equation, t is the time-step and the notation AB(.) means an Adams-Bashforth evaluation of the quantity: 4. The density k is computed from the equation of state (4); 5. The velocity divergence ( .V) k is computed by Eq. ( 12).

The steps 2 to 5 are repeated until convergence is achieved on temperature and density. The convergence criterion used is Max <10 -11 , with Res =Max(( kk-1 )/ k-1 ) for and the maximum number of iterations is fixed to 150.

-2 nd step: Computation of (V, Pdyn)

When starting the solution of the Navier-Stokes equations, temperature, density and velocity divergence at the current time step n+1 are known. It is therefore possible to use a projection-type algorithm such as those developed for the solution of incompressible Navier-Stokes equations, with some modifications. In this work, we have used the modified projection method proposed by [START_REF] Hugues | An improved projection scheme applied to pseudospectral methods for the incompressible Navier-Stokes equations[END_REF] and analysed in details in [START_REF] Raspo | A spectral projection method for the simulation of complex three -dimensional rotating flows[END_REF]. It was shown in particular that the modification introduced by

Hugues and Randriamampianina, namely the computation of a preliminary pressure, improves the accuracy on pressure and allows to reduce the slip velocity compared to other projection methods [START_REF] Raspo | A spectral projection method for the simulation of complex three -dimensional rotating flows[END_REF]. The modified projection method consists in solving the Navier-Stokes equations by three successive steps as follows:

Computation of a preliminary pressure First, a preliminary pressure

W V V V V V e V
where W n+1 represents the velocity boundary condition at time step n+1. At the channel outlet, the discretization of Eq. ( 11) leads to:

1 1 2 4 3 2 n n n n t Pr Re t Ra x V V V W

Computation of a predicted velocity V *

The predicted velocity field V * is computed implicitly from the momentum equation with the gradient of the preliminary pressure instead of that of the actual pressure

V V V V V V V e V W
Here again, we have to solve Helmholtz equations with coefficients involving n+1 for each velocity component. As for the energy equation, the density n+1 is split into a constant part and a time-dependent part and Helmholtz equations with constant coefficients are solved iteratively for velocity components. The convergence is achieved when Max(Resu, Resv) <10 -13 . This criterion is fulfilled after 3 or 4 iterations.

Correction step

The converged velocity field V * is then corrected by taking into account the pressure gradient at the current time step n+1 so that the final velocity field satisfies the continuity equation ( 1). This correction is performed through the solution of the following Poisson problem for the intermediate variable [START_REF] Le | Direct numerical simulation of turbulent flow over a backwardfacing step[END_REF].

The actual velocity field and pressure at the current time step n+1 are finally calculated in by the formulae:

1 * 1 1 n n V V , 1 1 3 2 n n dyn dyn P P t .

RESULTS AND DISCUSSION

Simulations were performed for a Reynolds number Re=50.57. The passage from the cold entry zone to the heated zone of the channel induces a perturbation of the flow. If the flow is stable, this perturbation decreases as it moves downstream and it leaves the system.

On the other hand, when the flow is unstable, this perturbation grows and it gives rise to thermoconvective patterns. Instantaneous temperature fields at several times are presented in order to show the spatio-temporal evolution of these structures.

First, the fluid was set at 1K above its critical temperature, which corresponds to a dimensionless distance to the critical point 3.288×10 -3 , and several values of the bottom wall heating were considered. These results were obtained in a channel with an aspect ratio =10. Then, the influence of the proximity to the critical point was investigated in order to study the effect of the Prandtl number variations on the thermoconvective instability. In this case, in order to track further the different structures in the channel, we increased the aspect ratio to =15.

Thermoconvective instabilities at 1K from the critical point

Computations were carried out for ranging from 0.02K to 0.1K. Because of the divergence of the thermal expansion coefficient and of the isobaric specific heat Cp near the critical point, these small temperature increases induce large Rayleigh numbers ranging from 1.05×10 6 to 5.23×10 6 . Figure 3 shows the temperature fields at different calculation times for several values of . In all the cases, instability appears in the form of thermal plumes which develop on the hot boundary layer. These structures are similar to those previously obtained for a supercritical fluid in the Rayleigh-Bénard configuration [START_REF] Amiroudine | Direct numerical simulation of instabilities in a two-dimensional near-critical fluid layer heated from below[END_REF][START_REF] Chiwata | Thermal plumes and convection in highly compressible fluids[END_REF][START_REF] Furukawa | Convective heat transport in compressible fluids[END_REF][START_REF] Furukawa | Convection in a very compressible fluid: comparison of simulations with experiments[END_REF][START_REF] Amiroudine | Piston-effect-induced thermal oscillations at the Rayleigh-Bénard threshold in supercritical He[END_REF][START_REF] Raspo | Unsteady two-dimensional convection in a bottom heated supercritical fluid[END_REF][START_REF] Accary | Three-dimensional Rayleigh-Bénard instability in a supercritical fluid[END_REF]Accary et al., 2005a;Accary et al., 2005b). However, in a closed cavity, the bottom wall heating induces a piston effect that homogeneously increases the bulk temperature and that gives rise to a second unstable thermal boundary layer on the top cold wall. This is not the case in an open channel as it can be seen in Fig. 3. The couple ,

ins ins t x
(where ins t is the time corresponding to the beginning of the isotherms deformation and ins x the axial position of the first plume) which characterizes the instability onset is respectively (4 1.22s), (3 , 0.71s), (1 , 0.33s) for 0.02K, 0.04K and 0.1K. So the thermoconvective instability appears earlier and closer to the beginning of the heated zone when the Rayleigh number is increased. Furthermore, as the heating increases, the disturbances grow faster and the thermal plumes are larger.

Figure 3 also shows that the number and the development of the thermoconvective structures are different according to the value of . While instability appears in the form of one or some plumes which move downstream during time for about 10 -2 K (cases A and B), many structures develop in the whole channel for the largest heating (case C). In addition, in this last case, as the structures move downstream, new thermal plumes continuously appear near the beginning of the heated zone. This time evolution of the patterns may indicate that the structures are due to a convective instability in cases A and B and to an absolute instability in case C. Indeed, according to the previous studies for incompressible fluids, when the flow is convectively unstable, the initial perturbation locally grows, it moves downstream and it finally leaves the system. Hence, the thermoconvective structures cross the channel during time and the flow becomes stable again once they moved away. On the other hand, in the case of absolute instability, the initial perturbation locally grows but it also expands in the whole system so that the thermoconvective structures appear throughout the channel all the time although they move downstream: when a structure moves away, another one is locally generated in its place by the perturbation. However, results reported in literature for incompressible fluids showed that longitudinal rolls can never be absolutely unstable [START_REF] Carrière | Convective versus absolute instability in mixed Rayleigh-Bénard-Poiseuille convection[END_REF]. Therefore, if the patterns in Nevertheless, the above results of linear stability analysis cannot be directly applied to our case since they were obtained with the assumption of incompressibility and constant physical properties, which is totally false for supercritical fluids.

Simulations were performed for a temperature increase 0.27K, corresponding to a Rayleigh number Ra=5.23×10 5 . Because of the strong decrease of the thermal expansion coefficient and of the isobaric specific heat as the distance to the critical point increases (see table 1), the Rayleigh number is much smaller than the smallest value that we considered at 1K from Tc (case A). Figure 5 shows the temperature fields obtained at two calculation times.

It reveals that the thermoconvective structures develop much faster compared to case A of Fig. 3 and that they are also much larger. Results obtained in section 4.1 for Ti =Tc +1K

showed that, when the flow is unstable, the disturbances grow faster and the thermal plumes are larger as Ra increases. Therefore, the fast and strong growth of the perturbations in Fig. 5 means that the value Ra=5.23×10 5 is much larger than the critical value Rac1 of the instability threshold for Pr=6.28. Two conclusions can be drawn from these results:

-first, contrary to what can be expected, the flow is convectively unstable for smaller values of Ra when moving away from the critical point. This is due to the decrease of the Prandtl number;

-secondly, the threshold of convective instability depends on the Prandtl number, since it depends on the distance to the critical point. Therefore, according to the previous studies [START_REF] Müller | Transversal convection patterns in horizontal shear flow[END_REF][START_REF] Carrière | Convective versus absolute instability in mixed Rayleigh-Bénard-Poiseuille convection[END_REF], this result confirms that the thermoconvective structures observed are likely transversal rolls, since the critical Rayleigh number for longitudinal rolls is independent of Pr. However, experiments on incompressible fluids and perfect gas showed that transversal rolls only exist for small values of Re (Re<10).

Therefore, the present results would indicate that, in a supercritical fluid, the transition between transversal rolls and longitudinal rolls occurs for larger Reynolds numbers than for perfect gas or incompressible fluids.

CONCLUSION

In this paper, the Poiseuille-Rayleigh-Bénard problem for a supercritical fluid, modeled by the Peng-Robinson equation of state, was studied for a Reynolds number Re=50.57. The temperature fields obtained in a channel of height =1mm revealed the existence of thermoconvective instabilities appearing in the form of thermal plumes, similar to those previously observed in supercritical fluids in the Rayleigh-Bénard configuration [START_REF] Chiwata | Thermal plumes and convection in highly compressible fluids[END_REF][START_REF] Amiroudine | Direct numerical simulation of instabilities in a two-dimensional near-critical fluid layer heated from below[END_REF][START_REF] Furukawa | Convective heat transport in compressible fluids[END_REF][START_REF] Furukawa | Convection in a very compressible fluid: comparison of simulations with experiments[END_REF][START_REF] Amiroudine | Piston-effect-induced thermal oscillations at the Rayleigh-Bénard threshold in supercritical He[END_REF][START_REF] Raspo | Unsteady two-dimensional convection in a bottom heated supercritical fluid[END_REF][START_REF] Accary | Three-dimensional Rayleigh-Bénard instability in a supercritical fluid[END_REF]Accary et al., 2005a;Accary et al., 2005b). These thermal plumes develop on the hot boundary layer and then move downstream during time. Thus, contrary to the perfect gas case, for which the convective structures occupy the whole height of the channel, instability develops here on the fine thermal boundary layer along the heated zone. The temporal evolution of the patterns indicates that the instability is convective or absolute depending on the value of the Rayleigh number. As a consequence, since longitudinal rolls can never be absolutely unstable [START_REF] Carrière | Convective versus absolute instability in mixed Rayleigh-Bénard-Poiseuille convection[END_REF], the thermoconvective structures obtained are likely transversal rolls.

Moreover, the results obtained at different distances from the critical point showed that the critical Rayleigh number for the convective instability threshold depends on the Prandtl number. This variation of the critical Rayleigh number with the Prandtl number is also characteristic of transversal rolls according to previous studies [START_REF] Müller | Transversal convection patterns in horizontal shear flow[END_REF][START_REF] Carrière | Convective versus absolute instability in mixed Rayleigh-Bénard-Poiseuille convection[END_REF].

Obviously, the results presented here are preliminary results and further simulations are necessary to confirm the type of structures that we obtained. In particular, simulations will be performed to determine the convective instability threshold for various values of the Reynolds number as well as the Reynolds number for the transition between transversal rolls and longitudinal rolls. Three-dimensional computations will also be carried out.

Figure captions

  -Bénard (PRB) flows are mixed convection flows in horizontal channels heated from below. When the Rayleigh number exceeds a critical value, thermoconvective structures develop in the channel. This paper focuses on these thermal instabilities in the case of fluids near their gas-liquid critical point. More precisely, a supercritical fluid (SCF) is a fluid for which temperature and pressure are larger than those of the critical point. In this zone of the phase diagram, the physical properties of the fluid (density, viscosity, diffusivity) are intermediate between those of liquids and gases.

  it was proposed inAccary et al. (2005c), the stratification of the fluid is taken into account leading to the initial condition for density and pressure:

  thermal conductivity k is updated;

V

  It must noted that the above Poisson-Neumann problem must fulfill a compatibility condition to be solvable. Because of the inhomogeneous Dirichlet condition for the predicted velocity V * at the channel outlet, this compatibility condition is not fulfilled. A renormalization of V * at the outlet boundary was therefore performed before the solution of the Poisson-Neumann problem as proposed in

  Fig. 4. In this figure, the critical Rayleigh number Rac1 for convective instability onset was

Figure 1 :

 1 Figure 1: Linear stability diagram of the PRB flow for incompressible fluids showing the

Figure 2 :

 2 Figure 2: Channel geometry with top and bottom thermal boundary conditions.

Figure 3 :

 3 Figure 3: Temperature fields for Ti =Tc +1K (Pr=31.72) and several values of Ra.

Figure 4 :

 4 Figure 4: Linear stability diagram of the PRB flow for incompressible fluids showing the

Figure 5 :

 5 Figure 5: Temperature field for Ti =Tc +5K (Pr=6.28) and Ra=5.23×10 5 .

  Figure 1

  Ma 2 /Fr (where Ma and Fr are respectively the Mach and the Froud numbers) is not in O(Ma 2 ). We chose as characteristic variables Tc for temperature, c for density, c

				c for
	pressure (with	188.92J.kg -1 .K -1 the perfect gas constant),	for length, ref U	g T H
	for velocity (with	the Earth gravity),

ref H U

Table titles Table 1 :

 titles1 Thermodynamic properties and Prandtl number at 1K and 5K above the critical point.

	Table 1
	i T T c	(K)	(K -1 )	Cp (J/kg/K)	Pr
	1		0.82199	42803.48	31.72
	5		0.15334	8478.77	6.28

Acknowledgements.

for its financial support (White Program ANR-09-BLAN-0105-01). This work was performed using HPC resources from GENCI-IDRIS (grants 2011-20321).


S Henry

J M Faure

Elaboration of invariant safety properties from fault-tree analysis

Formal verification of PLC programs using model-checking requires to elaborate previously temporal logic formulae that state in a formal way the properties that must be checked. Unfortunately temporal logic is a formalism totally unknown by automation engineers. This explains why PLC programs developers willing to verify the behavior of their programs are unable to use the existing modelchecking tools. Temporal logic formulae elaboration is a too difficult task. In order to overcome this problem and to bridge the gap between PLC programs development and model-checking, this paper proposes a methodology enabling to obtain invariant safety properties from fault-tree analysis. Fault-tree analysis is a quite popular analysis method often used in industry when designing critical systems. Hence using the results of this analysis to elaborate formal properties may contribute to increase the use of formal verification techniques.

I. INTRODUCTION

Programmable Logic Controllers (PLC) are widely employed nowadays to implement industrial control systems. The first PLC systems have been used in the 1980s to control simple processes, mainly in production systems where they have replaced hard-wired switching networks based on relays. The new capabilities of PLCs, the comfort of PLC programming languages and of programming environments as well as economic constraints have led to an increasing use of PLC-based systems in all the industrial fields and even for safetyrelated applications. Examples can be found in transport systems, in critical industries like chemical processes, oil industry, power production and distribution. The growing demand of the society for risk prevention combined with this expanding use of PLC systems explain why PLC programs verification [FAU,01], [FRE,00] is to-day a major industrial concern. PLC program verification does not aim merely at checking the intrinsic properties of the program, e.g. no infinite loop, no locking point, …, regardless the application, but mainly at checking that the program behaves as required. This paper focuses only on this last kind of property: the compliance of a given PLC program with the properties required for the application. Moreover these properties are often ranked in two categories: safety properties (what shall not be done) and liveliness properties (what must be done). In this paper only safety properties will be dealt with.

The LURPA has achieved since ten years several works on PLC programs verification using model-checking. A significant result of these works is a formal semantics of the SFC and ladder diagram languages [LAM, 00], [DES,02]. This formal definition enables to represent any program written in these languages as a state automaton and consequently to translate such programs into modelcheckers compliant languages. The underlying theory of model checking [MAC,93] [SCH,99] is indeed state automata theory. Hence, using a model-checker for PLC program verification implies to have a state model of the program as well as formal properties, in some temporal logic, that expresses in a formal way the application requirements and that must be satisfied by the program (Fig. 1).

Obtaining these formal properties is unfortunately a major problem. The application requirements are indeed expressed in industry in a quite informal way, i.e. some sentences in natural language or drawings., but never with sound mathematical formulae. Moreover, temporal logic [PNU,77] [EME,90] is not a classical logic and is difficult to grasp for automation engineers. It really matters to highlight that the use of temporal logic for expressing properties acts as a brake for the industrial use of model-checking of PLC programs.

Modelling

State automaton

The works presented in this paper are aimed at easing the elaboration of formal properties of PLC programs. Other methods, developed with the same objective, are presented in [FIL,99], [JUL,99], [KLE,01]. Though they may help usefully designers looking for formal properties, they require some experience and are not directly connected to any existing industrial method for dependability or safety analysis. On the opposite the method presented hereafter has been developed in order to link formal properties elaboration to a method commonly employed in industry: fault-tree analysis. This constraint has been put to facilitate the acceptance of industrial users.

In the scope of this paper, only invariant safety properties will be considered. This paper is structured as follows. The section 2 gives an overview on the developed method. A small example used in the rest of the paper to illustrate this method is presented in the section 3. The section 4 shows the design of a fault-tree taking into account PLC program faults. Formal statements are then derived from this fault-tree in section 5. The scope of these statements as well as program invariants elaboration are discussed in the last section.

II. METHOD OVERVIEW

Lots of techniques have been developed in order to assess or to enhance systems dependability when specifying and designing [VIL,88]. Fault-tree analysis [IEC, 90] has been chosen for the following reasons:

• This analysis technique is widely used for critical systems design. • It aims to find the different causes of a given fault.

• It enables to determine how failures combinations may lead to faults.

• Several research works have been developed to facilitate fault-trees design from functional analysis [DUT,95] or to generate automatically fault-trees [GAL,99] [GAU,01].

The objective of the work presented in this paper is to take benefit of the results of a fault-tree analysis to elaborate invariant safety properties. Figure 2 depicts the proposed method. The general idea is to perform from the requirements a fault-tree analysis taking into account PLC programs faults and then to deduce temporal logic formulae from the obtained fault-trees.

III. ILLUSTRATIVE EXAMPLE

This method for properties elaboration will be presented thanks to a small example used in the framework of a benchmarking project by the French working group COSED [COS]. This example, sketched in figure 3, is a pick-and-place manipulator, part of an assembly line located at the Mechanical Engineering Department of the ENS Cachan. The aim of this manipulator is to pick up gearwheels with suction cups and to transfer the gearwheels to gear housings using two pneumatic cylinders. Only automatic operations will be considered.

The controlled system may be decomposed into three parts:

-Horizontal movement:

1 double-acting cylinder 1 dual-solenoid valve 2 magnetic sensors (rightmost and leftmost positions)

-Vertical movement: 1 double-acting cylinder 1 dual-solenoid valve 2 magnetic sensors (upper and lower positions)

-Drawing up: suction cups with a venturi system 1 single-solenoid valve no sensor

The control system is implemented by a control program running on a PLC. The inputs and outputs of this program are given in figure 3.

In order to avoid collisions with mechanical systems located between the two stations, all the horizontal movements must be performed with the drawing up system in the upper position. Hence any downward movement between the two stations shall be forbidden and considered as a fault during safety analysis. The aim of fault-tree analysis [IEC,90] is to determine all the possible causes of a considered fault, called "undesirable event". To reach this objective, fault-tree analysis looks for how combinations of elementary faults (faults that do not come from other faults) can lead to this fault. The result of a fault-tree analysis can be represented in the form of a tree including faults and logical gates. The root of this tree is the "undesirable event", its leaves are the elementary faults. The figure 4 presents an extract of a fault-tree designed for the "Collision with the environment at the picking station" fault. Similar analysis can be carried out for other faults, like: "Collision with the environment at the placing station" or "Drawing up system fault".

In a classical fault tree [NIE, 02], any fault is related to a physical component failure or to an operator fault or error.

When dealing with automatic systems, faults can come from:

• The controlled system. In this case only physical components failures are considered. • The control system. Faults can be then either physical failures of the components of this system (sensors, input/output cards, processor, …), or faults or errors of the control software.

It is the reason why the fault-tree depicted in the figure 4 includes physical components failures and software faults. The fault "Move down with Leftmost Position false" is related to the control system and is a fault of the control software, whilst the fault "Move down with Leftmost Position information true", is also related to the control system but is a failure of an hardware component.

Hence, the obtained fault-tree includes three kinds of faults: control software faults, physical components failures and faults that are combination of these two kinds of faults (Move down without being at the leftmost position for instance). In order to elaborate properties dealing only with control software faults, control program faults shall be decoupled from other faults (coming from physical components failures). This can be achieved by introducing the following constraint, that the fault-tree shall verify: Any control program fault must be combined with physical components faults only using OR logical gates.

Fig. 4. Fault-tree for the "Collision with the environment at the picking station" fault

V. FROM FAULT-TREE TO FORMAL STATEMENTS

A System properties

Assuming that the global system, that includes the controlled system and the control system, can be described in the form of a state automaton, the following property in temporal logic can be derived from the previous fault-tree:

AG ¬ (Collision with the environment at the picking station) This property may be interpreted as follows: "For all paths, in each state, it is impossible to move down without being at the leftmost position and to move to the left without being at the upper position", and be rewritten as a system of two properties: "For all paths, in each state, it is impossible to move down without being at the leftmost position" and "For all paths, in each state, it is impossible to move to the left without being at the upper position": AG (¬Move down without being at the leftmost position) AG (¬Move to the left without being at the upper position)

INPUTS

When a fault is indeed the output of an OR gate, this fault will never occur if none of the faults inputs of this gate occurs.

B Control software properties

In order to obtain program properties, only software faults will be kept in the following. In the case of the studied fault-tree, the two hereafter mentioned faults will be retained:

• Move down with Leftmost Position false

• Move left with Upper Position false.

The root of this fault-tree (the "undesirable event") is indeed the output of a OR logical gate. Each input of this gate is itself output of a OR gate. This structure enables to state that the undesirable event will never occur if none of the control program faults occurs, that leads to the following properties system:

AG ¬( Move down with Leftmost Position false)

AG ¬( Move left with Upper Position false)

This system means by reasoning in the state automaton describing the program, that, for all paths, in each state, the two properties must hold to avoid the fault for which the fault-tree analysis has been made.

C Taking into account PLC inputs-outputs

It is now possible to replace textual expressions by the input and output variables of the control program, that leads to:

AG ¬( S3 ∧ ¬e1 ) (1) AG ¬( S2 ∧ ¬e3 ) (2)
VI. PROGRAM INVARIANTS ELABORATION

The properties system obtained in the previous section has been deduced from a fault-tree related to collision at the picking station. Nevertheless checking invariant properties of a program implies that the properties to check shall hold for all the states of the automaton representing the program, i.e. in our case when the drawing up system is at the placing station or moves from a station to the other.

It matters now to ensure that the above obtained properties are really invariants and, if not, to modify them so as to obtain real invariants.

A. Model-based verification: introducing conditions on inputs

[FRE, 00] describes two approaches for PLC programs verification:

• The first one does not use any plant model. All inputs combinations are possible. • The other one uses a plant model that limits inputs combinations to physically plausible ones. Hence constraints on inputs are introduced, stating for instance that two sensors used to indicate the rightmost and leftmost positions of a cylinder must never give simultaneously a true signal.

The main advantage of this second approach is to reduce the number of reachable states and therefore the checking time. It is the reason why this approach is often used and will be taken in what follows. In the case of the studied example, it is possible to introduce the two following conditions that constraint PLC input variables:

e1 ∧ e2=0 (c1)
the horizontal cylinder can not be at the same time at the rightmost position and at the leftmost position e3 ∧ e4=0 (c2) the vertical cylinder can not be at the same time at the upper position and at the lower position.

B. Scope of the previous properties

The properties system (1), (2) has been deduced considering the "Collision with the environment at the picking station" fault. But a collision with the environment can also occur at the placing station. From the fault-tree of the "Collision with the environment at the placing station" a system of two invariant properties can be derived in a similar way:

AG ¬( S3 ∧ ¬e2 ) (3) AG ¬( S1 ∧ ¬e3 ) (4)
The first attempt to obtain a program invariant is to consider a generic fault "Collision with the environment" and to merge the four previously obtained properties (1), ( 2), ( 3) an (4):

AG ¬( Collision with the environment)

AG ¬[(S3∧¬e1)∨(S3∧¬e2)∨(S1∧¬e3)∨(S2∧¬e3)]
AG ¬(S3∧¬e1)∧¬(S3∧¬e2)∧¬( S1∧¬e3)∧¬( S2∧¬e3)

AG (¬S3∨e1)∧(¬S3∨e2)∧¬(S1∧¬e3)∧¬(S2∧¬e3) AG [¬S3∨(e1∧e2)]∧¬(S1∧¬e3)∧¬(S2∧¬e3)
It is easy to see that this attempt was wrong. The condition c1 [e1 ∧ e2] is always false. The property is true only if ¬S3, i.e. if the drawing up system never goes down.

The properties systems previously obtained are true only in the conditions for which the fault-tree analysis has been performed. The system (1), ( 2) is true only for a collision at the picking station and the system (3), ( 4) is true only for a collision at the placing station. It is not allowed to merge them so as to obtain a program invariant, true whatever the state of the controlled system would be.

C. Program invariants

In order to obtain program invariants (properties true for all states of the program), previous statements are to be rewritten by introducing a condition indicating from which fault-tree analysis they have been deduced, i.e. for which situation of the plant (controlled system) there are true. This implies to add to each present formula coming from a given fault-tree an other present formula stating the situation of the plant considered when designing the faulttree. Hence program invariants will be written in the following way:

AG ¬(present formula #1∧ present formula #2)

The first formula is merely that deduced from the faulttree. The second one indicates the plant situation.

In the case of the studied example, a collision at the picking station does not imply that the drawing up system is strictly at this station, but that it arrives, leaves or is waiting for at this station. Hence this situation must be translated by "Not at the placing station" with a control point of view and will be characterized with PLC inputs by using the ¬e2 condition. The properties system (1) and

(2) becomes consequently:

AG ¬( S3 ∧ ¬e1 ∧ ¬e2) (5) AG ¬( S2 ∧ ¬e3 ∧ ¬e2) (6)
In a similar way, "Collision at the placing station" will be translated by ¬e1 (not to be at the picking station). The properties system (3) and ( 4) becomes consequently:

AG ¬( S3 ∧ ¬e2 ∧ ¬e1) (7) AG ¬( S1 ∧ ¬e3 ∧ ¬e1) (8)
The property ( 5) is identical to the property ( 7) and means that the control program shall never set the "Move down" output if neither the input variable e1 (drawing up system at the leftmost position), nor the variable e2 (drawing up system at the rightmost position) are true. Properties ( 6) and ( 8) mean respectively that no movement to the left (to the right) can occur if he drawing up system is not at the upper position.

These four properties are the searched program invariants.

VII. CONCLUSION

This paper has presented the first results of a work whose objective is to facilitate the elaboration of invariant safety properties using a preliminary fault-tree analysis. The presented method has been tested on other examples, e.g. the control of a water distribution system including several pumps. The main difficulty of this method lies in the scope of the statements directly obtained from the fault-tree. Nevertheless it has been shown in section VI that it is possible to obtain from these statements program invariants by introducing a formula clearly stating the scope of the fault-tree analysis from which the statement has been obtained.

The perspectives of this work are numerous. First of all the robustness of this method has to be evaluated thanks to several examples. The next step will aim at easing dynamic properties elaboration. This implies to design fault-trees that do not include only logical gates, but also gates dealing with time and events sequences. Some works taking into account that kind of specific elements are presented in [HEN,02]. The last step, that must be carried out in strong cooperation with safety engineers and PLC developers, may be the development of a software tool supporting the method and linked to existing fault-tree design commercial software.

Fig. 1 .

 1 Fig. 1. Formal verification using model-checking

Fig. 2 .

 2 Fig. 2. Using fault-tree analysis to elaborate formal properties

Fig. 3 .

 3 Fig. 3. Studied example

¬

  that the considered fault must never happen or, with other words, that for all paths, for each state, this property shall always hold. Using the decomposition provided by the fault-tree it is then possible to obtain another equivalent expression: AG ¬ (Move down without being at the leftmost position ∨ Move to the left without being at the upper position) That becomes when applying the De Morgan theorem: AG (¬Move down without being at the leftmost position ∧ Move to the left without being at the upper position)