Nadezhda Kononova

Vyacheslav Shevchenko

Alexander Kokh

Tatiana Nabeeva

David Chapron

Alain Maillard

Asset Bolatov

Bolat Uralbekov

Dy, Ho, Er, Tm Lu Nadezhda Kononova

V S Sobolev

Synthesis of New Isostructural Orthoborates NaBaR(BO 3 ) 2 with R = Tb

Keywords: rare-earth orthoborates, NaBaR(BO 3 ) 2, solid-state synthesis, Raman spectra, luminescence

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The rare-earth borates have potential applications for laser medium, luminescent and nonlinear optical materials [START_REF] Ye | Growth of nonlinear optical crystal Y 0.57 La 0.72 Sc 2.71 (BO 3 ) 4[END_REF][START_REF] Sasaki | Recent development of nonlinear optical borate crystals: key materials for generation of visible and UV light[END_REF][START_REF] Jubera | Luminescent properties of Eu 3+ -activated lithium rare earth borates and oxyborates[END_REF] . Many one-and two cation borates are currently used in photonics [START_REF] Becker | Borate materials in nonlinear optics[END_REF] . A current area of research is for new materials for photonic dilates using more complex borate crystal compounds. Orthoborate crystals containing three cations are among the most attractive to investigate because of their potential application as a light emitting phosphor for phosphor-converted white light emitting diodes [START_REF] Lian | Crystal structure refinement and luminescence properties of Ce 3+ singly doped and Ce 3+ /Mn 2+ co-doped KBaY(BO 3 ) 2 for n-UV pumped white-light-emitting diodes[END_REF] and as a green emitting phosphor for three-dimensional plasma panel [START_REF] Han | Observation of efficient energy transfer from host to rare-earth ions in KBaY (BO 3 ) 2 :Tb 3+ phosphor for plasma display panel[END_REF] .

In our earlier study of solid phase formation in the system M 2 O 3 -BaO-Na 2 O-B 2 O3 (M = Sc, Y) we identified two new borates: NaBaSc(BO 3 ) 2 and NaBaY(BO 3 ) 2 7-8 . Further investigation in this family resulted in the discovery of new ytterbium-bearing borate NaBaYb(BO 3 ) 2

9

. All discovered compounds of NaBaR(BO 3 ) 2 (R = Sc, Y, and Yb) line are isostructural and have eitelite-like, Na 2 Mg(CO 3 ) 2 , structure with space group.

In this study, we report the synthesis of new rare-earth compounds in NaBaR(BO 3 ) 2 orthoborate family (R = Tb, Dy, Ho, Er, Tm, Lu) using high-temperature solid-state reaction, and continue previously described orthoborates NaBaR(BO 3 ) 2 with R = Sc, Y and Yb. The photoluminescence properties of NaBaTb(BO 3 ) 2 are also reported.

Material and Methods

Samples preparation and Synthesis

High purity chemicals of BaCO 3 , Na 2 CO 3 ,

Н 3 BO 3 , La 2 O 3 , Ce 2 O 3 , Nd 2 O 3 , Pr 2 O 3 , Eu 2 O 3 , Gd 2 O 3 , Tb 2 O 3 , Dy 2 O 3 , Ho 2 O 3 , Er 2 O 3 , Tm 2 O 3 , Lu 2 O 3 (
analytical reagent) were purchased and used as starting reactants. Samples of regents were weighed at the stoichiometric ratio, mixed thoroughly and ground into fine powder in an agate mortar. Resulting mixtures with an excess of 3 mass % H 3 BO 3 were heated to 700 0 C in platinum crucible and kept at this temperature for 24 h to decompose the carbonates and boron acid.

Multi-stage high temperature solid-state synthesis was used to prepare rare-earth borates. During the first stage, the starting mixture was annealed at 700 o C for approximately 1 day. The conditions of the second stage were achieved by increasing the temperature in steps of 50°С. After each heating step, intermediate reaction products were sampled and identified by gradually grinding following X-ray powder diffraction analysis to identify crystal forms. Phase equilibrium, indicated by unchanging reflections and intensities in the X-ray diffraction patterns between temperature steps, was served as a criterion for completion of the solid-state reaction. Heating to higher temperatures was observed to result in partial or complete melting.

Samples characterization

X-ray powder diffraction (XRD) data was obtained by a DRON-3 diffractometer with 2θ ranging from 10 0 to 80 0 (CuKα radiation; U = 30-40 kV; I = 25 mA, speed meter -1°/min, belt speed -1 cm/min) using silicon as an external standard. Unit cell parameters for new borates are calculated based on X-ray powder diffraction patterns. In addition, Raman spectra were collected using Lab RAM ARAMIS Horiba Jobin Yvon Raman spectrometer producing light at 633 and 785 nm. The photoluminescence properties of NaBaTb(BO 3 ) 2 were studied with Cary Eclipse Fluorescence spectrophotometer at room temperature.

Results and Discussion

Table 1 summarizes results of phases identified after high temperature solid-state synthesis of rare-earth borates containing sodium and barium. In some cases, we obtained simple and complex rare-earth borates. For example, a Ce 3+ solid state reaction resulted in formation of CeBO 3 simple borate, while Nd 3+ produced Nd 2 Ba 3 (BO 3 ) 4 and NdNa 3 (BO 3 ) 2 complex borates. Using La 3+ , Pr 3+ , Eu 3+ produced La 2 Ba 3 (BO 3 ) 4 , Pr 2 Ba 3 (BO 3 ) 4 , Eu 2 Ba 3 (BO 3 ) 4 , respectively. These latter reactions did not result in a compound isostructural to previously obtained NaBaSc(BO 3 ) 2 , probably because the ionic radii of cerium subgroup rare-earth elements are too large.

For Gd 3+ we observed intermediate products of reaction:

( ) Gd BaO Na O B O NaBaBO GdBO 2 2 2 2 2 1 2 2 3 3 3 + + + = +
Figure 1 shows observed X-ray powder diffraction patterns of the new NaBaR(BO 3 ) 2 (R = Tb, Dy, Ho, Er, Tm and Lu) together with previously prepared orthoborates of NaBaYb(BO 3 ) 2 for comparison. All prepared borates are isostructural as evidenced by the presence of the same diffraction peaks in X-ray patterns. The cell parameters of the new synthesized compounds NaBaR(BO 3 ) 2 were calculated on the basis of the X-ray powder diffraction data using Cellcalc software [START_REF] Miura | CellCalc: A unit cell parameter refinement program on Windows computer[END_REF] with space group . The cell parameters for all new compounds are given in Table 2. Based on these results, it can be seen that the lattice parameters and volume of NaBaR(BO 3 ) 2 decrease with the decreasing of the effective ionic radius of rare earth elements according to Shannon [START_REF] Shannon | Revised effective ionic radii and systematic studies of interatomic distances in halides and chalcogenides[END_REF] .

Using these cell parameters, it can be concluded that synthetic orthoborates are isostructural to NaBaSc(BO 3 ) 2 , NaBaY(BO 3 ) 2 , NaBaYb(BO 3 ) 2, crystallographic data of which have been fully described elsewhere [START_REF] Seryotkin | Synthesis and crystal structure of new layered BaNaSc(BO 3 ) 2 and BaNaY(BO 3 ) 2 orthoborates[END_REF][START_REF] Svetlyakova | Synthesis, growth and some physical properties of new orthoborates ScBaNa(BO 3 ) 2 and YBaNa(BO 3 ) 2[END_REF][START_REF] Svetlyakova | Search for compounds of the NaBaR(BO 3 ) 2 family (R= La, Nd, Gd, and Yb) and the new NaBaYb (BO 3 ) 2 orthoborate[END_REF] . The main features of their crystal structure are anion-cation (Ba,Na)(BO 3 ) layers, which are combined into the base building packages of two types: {R 3+ [Ba 2+ (BO 3 ) 3-] 2 } + and {R 3+ [Na + (BO 3 ) 3-] 2 } -. The fundamental building units of these orthoborates are RO 6 octahedra and BO 3 triangles. Furthermore, how it was determined in ref. [START_REF] Seryotkin | Synthesis and crystal structure of new layered BaNaSc(BO 3 ) 2 and BaNaY(BO 3 ) 2 orthoborates[END_REF] the new compounds are isotypic to mineral eitelite, Na 2 Mg(CO 3 ) 2 [START_REF] Pabst | Synthesis, properties, and structure of K 2 Ca(CO 3 ) 2 , buetschliite[END_REF] with the same space group . XRD patterns depict a marked shift of the X-ray peaks in the transition from Sc-to Tb-compounds due to the differences between the ion radii of studied rare-earth elements, from 0.745 to 0.923, respectively.

Based on XRD results of measured intermediate solidstate reaction products, the reaction of the borates formation can be represented as follows:

( ) R O BaO Na O B O NaBaBO RBO 2 2 2 2 2 2 3 2 2 3 3 3 + + + = + ( ) ( ) NaBaBO RBO NaBaR BO 3 3 3 3 2 + =
Further interaction between the products of this reaction was not observed, since annealing of the sample higher than 950 0 C led to complete melting of solids.

Raman spectroscopy

The Raman spectra of the NaBaR(BO 3 ) 2 crystals were recorded with exciting laser at 633 nm and 785 nm for R = Tb, Dy, Y, Tm, Yb, Lu, Sc and Er, Ho, respectively. At 633 nm the laser was found to excite the energy levels of Er or Ho ions, interfering with the emission spectrum.

Results of Raman spectroscopy presented in Figure 2 and Figure 3 show that synthetic orthoborates are isostructural, evidenced by presence of the same modes in Raman spectra. The Raman spectra show typical vibration modes of the triangular BO 3 group. It is well known that isolated BO 3 ion with D 3h symmetry characterized by four fundamental modes of which three Raman-active v' 1 (A' 1 ) (~950cm -1 ), v' 3 (E') (~1250-1400 cm -1 ), v' 4 (E' 1 ) (~600 cm -1 ) and one infrared active [START_REF] Kasprowicz | Vibrational properties of nonlinear optical Bi 2 ZnOB 2 O 6 single crystals doped with Pr 3+ :μ-Raman investigations[END_REF] .

The peak at ~1200 cm -1 in the Raman spectra of synthesized orthoborates are likely due to the asymmetric stretching mode v' 3 ((E') strong peaks at ~950 cm -1 attributed to the symmetric stretching modes v' 1 (A' 1 ), while relatively high peak at ~600 cm -1 are from bending vibration mode of BO 3 group. Furthermore, the modes in the region of 420-125 cm -1 can be assigned to BO 3 translation and vibration 14 .

Photoluminescence properties

Previous studies showed that rare-earth elements introduced into borate structure may present photoluminescence properties [START_REF] Liu | Spectroscopic properties of rare earths in optical materials[END_REF] . Among synthesized NaBaR(BO) 3 borate family, photoluminescence properties of as-prepared NaBaTb(BO) 3 have been investigated in this study, because Tb 3+ in host materials shows an intense characteristic green emission and relative suitable decay time, which can be used for PDP and LEDs [START_REF] Han | Observation of efficient energy transfer from host to rare-earth ions in KBaY (BO 3 ) 2 :Tb 3+ phosphor for plasma display panel[END_REF][START_REF] Lian | Synthesis, crystal structure, characterization and luminescent properties of KBaTbB 2 O 6[END_REF] .

Figure 4 shows the PLE and PL spectra of NaBaTb(BO 3 ) 2 powder sample. As shown in Figure 4, there is a broad excitation band ranging from 200 to 300 nm with a maximum at about 250 nm and shoulder peaks at ~270 and 292 nm, which can be ascribed to 4f → 5d transitions of Tb 3+ . The PL spectrum excited at 260 nm is composed of several narrow emission peaks at ~490, 550, 580 and 620 nm due to the well-known 5 D 4 → 7 F J (J=6, 5, 4, 3) transitions [START_REF] Liu | Spectroscopic properties of rare earths in optical materials[END_REF] .

The high level emission ( 5 D 3 ) is quenched by cross relaxation in favour of the 5 D 4 emission due to high concentration of Tb 3+ in the compound [START_REF] Blasse | Energy Transfer[END_REF] . The shortest estimated Tb-Tb distance (~5.36 Å) suggests that exchange interaction becomes ineffective and luminescence is predicted to become efficient 17 , making the NaBaTb(BO 3 ) 2 compound a promising green-emitting phosphor material.

Conclusion

Six new rare-earth orthoborates belonging to NaBaR(BO 3 ) 2 crystal family were synthesized with rare earth elements Tb, Dy, Ho, Er, Tm and Lu, while several other trivalent (3+) cations: La, Ce,Pr, Nd, Eu, Gd were not. Based on calculated lattice parameters and volumes these six new borates are isostructural to previously described orthoborates NaBaR(BO 3 ) 2 with R = Sc, Y, Yb and fall within the trigonal system space group , and are isotypic with mineral eitelite, Na 2 Mg(CO 3 ) 2. Moreover, NaBaR(BO 3 ) 2 lattice parameters appear to decrease with the decreasing of the radius of rare earth elements. A strong green emission centered at 540-550 nm under UV excitation was observed for synthesized NaBaTb(BO 3 ) 2 showing promise as a new phosphor material.

Figure 1 :

 1 Figure 1: X-ray powder diffraction patterns of the new NaBaR(BO 3 ) 2 orthoborates (R = Tb, Dy, Ho, Er, Tm and Lu) and NaBaYb(BO 3 ) 2 as a reference.

Figure 2 :

 2 Figure 2: Raman spectra of the new NaBaR(BO 3 ) 2 (R = Tb, Dy, Y, Tm, Yb, Lu and Sc) at 633 nm excitation.

Figure 3 :

 3 Figure 3: Raman spectra of the new NaBaR(BO3)2 (R = Ho and Er) at 785 nm excitation.

Figure 4 :

 4 Figure 4: The excitation and emmission spectra of NaBaTb(BO 3 ) 2

Table 1 :

 1 Results of solid-state synthesis of NaBaR(BO 3 ) 2 orthoborates (R = Sc, Y, cerium (La-Eu) subgroup and yttrium (Gd-Lu) subgroup of REEs.

	№	R	Atomic number	R3+(coordination	Max annealing	Determined or known compounds after
				number -6), Å 10	temperature, о С	X-ray research
	1	Sc	21	0.745	900	NaBaSc(BO 3 ) 2 6,7
	2	Y	39	0.900	900	NaBaY(BO 3 ) 2	6,7
	3	La	57	1.032	1000	La 2 Ba 3 (BO 3 ) 4 , Na 2 B 2 O 4
	4	Ce	58	1.010	900	CeBO 3 , CeO 2 , BaNaBO 3
	5	Pr	59	0.990	900	Pr 2 Ba 3 (BO 3 ) 4 , BaNaBO 3
	6	Nd	60	0.983	900	Nd 2 Ba 3 (BO 3 ) 4 , NdNa 3 (BO 3 ) 2
	7	Eu	63	0.947	900	Eu 2 Ba 3 (BO 3 ) 4 , Na 2 B 6 O 10
	8	Gd	64	0.938	950	GdBO 3 , BaNaBO 3
	9	Tb	65	0.923	800	NaBaTb(BO 3 ) 2 NEW
	10	Dy	66	0.912	800	NaBaDy(BO 3 ) 2 NEW
	11	Ho	67	0.901	850	NaBaHo(BO 3 ) 2 NEW
	12	Er	68	0.890	800	NaBaEr(BO 3 ) 2 NEW
	13	Tm	69	0.880	800	NaBaTm(BO 3 ) 2 NEW
	14	Yb	70	0.868	850	NaBaYb(BO 3 ) 2 8
	15	Lu	71	0.861	800	NaBaLu(BO 3 ) 2 NEW

Table 2 :

 2 Lattice parameters of NaBaR(BO 3 ) 2 orthoborates.

	№	R	a, Å	c, Å	V, Å 3	ρ, g/cm 3
	1	Sc	5.23944(12)	34.5919(11)	822.38(4)	3.912
	2	Y	5.3338(2)	35.8303(19)	882.78(7)	4.140
	3	Tb	5.367(5)	36.30(7)	905(2)	4.806
	4	Dy	5.360(3)	36.16(4)	899(1)	4.874
	5	Ho	5.358(2)	35.87(4)	891(1)	4.946
	6	Er	5.349(2)	35.62(4)	882(1)	5.024
	7	Tm	5.343(2)	35.59(3)	880(1)	5.057
	8	Yb	5.337(3)	35.58(4)	878(1)	5.116
	9	Lu	5.336(3)	35.54(4)	876(1)	5.148

Acknowledgements

The authors would like to thank Dr. Kh. Abdullin for help with the luminescence analysis.


Eric Chénier
email: eric.chenier@univ-mlv.fr

Robert Eymard
email: robert.eymard@univ-mlv.fr

Raphaèle Herbin
email: herbin@cmi.univ-mrs.fr

Ouardia Touazi
email: ouardia.touazi@univ-mlv.fr

Collocated finite volume schemes for the simulation of natural convective flows on unstructured meshes

Keywords: cell centered finite volumes, incompressible Navier-Stokes, Boussinesq hypothesis, stabilization methods, unstructured meshes, natural convection between cylinders

We describe here a collocated finite volume scheme which was recently developed for the numerical simulation of the incompressible Navier-Stokes equations on unstructured meshes, in 2 or 3 space dimensions. We recall its convergence in the case of the linear Stokes equations, and we prove a convergence theorem for the case of the Navier-Stokes equations under the Boussinesq hypothesis. We then present several numerical studies. A comparison between a cluster-type stabilization technique and the more classical Brezzi-Pitkäranta method is performed, the numerical convergence properties are presented on both analytical solutions and benchmark problems and the scheme is finally applied to the study of the natural convection between two eccentric cylinders.

Introduction

Numerical schemes for the incompressible Navier-Stokes equations have been extensively studied: see [START_REF] Girault | Finite element methods for the Navier-Stokes equations: Theory and algorithms[END_REF][START_REF] Gunzburger | Finite element methods for viscous incompressible flows, A guide to thoery, practice, and algorithms[END_REF][START_REF] Glowinski | Numerical Methods for fluids (Part3)[END_REF]26] and references therein. An advantage of the finite volume schemes is that the unknowns are approximated by piecewise constant functions: this makes it easy to take into account additional nonlinear phenomena or the coupling with algebraic or differential equations, for instance in the case of reactive flows; in particular, one can find in [START_REF] Patankar | Numerical Heat Transfer and Fluid Flow[END_REF] the presentation of the classical finite volume scheme on rectangular meshes, which has been the basis of many industrial applications. Proofs of the convergence of the so-called "MAC scheme" [START_REF] Harlow | Numerical calculation of time dependent viscous incompressible flow of fluids with free surface[END_REF] were performed for the Stokes and Navier-Stokes equations, see [START_REF] Ph | Convergence of a finite volume scheme on a MAC mesh for the Stokes problem with right hand side in H -1[END_REF] and references therein. However, the use of rectangular grids makes an important limitation to the type of points satisfying the usual orthogonality condition (see Figure 1). The following notations are

K L x L x K d K|L d K,K|L n KL d K,σ m K|L mσ Figure 1.
Notations for an admissible mesh used. The size of the discretization is defined by: h D = sup{diam(K), K ∈ M}. The set of interior (resp. boundary) edges is denoted by E int (resp. E ext ), that is E int = {σ ∈ E; σ ⊂ ∂Ω} (resp. E ext = {σ ∈ E; σ ⊂ ∂Ω}). For any edge σ ∈ E, we denote by m σ the length or area of σ. For all K ∈ M and σ ∈ E K (the set of edges of K), we denote by n K,σ the unit vector normal to σ outward to K and we denote by d K,σ the Euclidean distance between x K and σ. For all K ∈ M, we denote by m K the surface or volume of K, N K the subset of M of the neighboring control volumes. For all K ∈ M and L ∈ N K , we denote by K|L ∈ E K ∩ E L the common edge of K and L; we then set n KL = n K,K|L , and we denote by d K|L = d K,K|L + d L,K|L the Euclidean distance between x K and x L . We shall measure the regularity of the mesh through the function regul(D) defined by regul(D) = inf dK,σ

diam(K) , K ∈ M, σ ∈ E K ∪ d K,K|L d K|L , K ∈ M, L ∈ N K ∪ 1 card(EK ) , K ∈ M . (1)
We then denote by H D (Ω) ⊂ L 2 (Ω) the space of functions which are piecewise constant on each control volume K ∈ M. For any u ∈ H D (Ω) and for all K ∈ M, we denote by u K the constant value of u in K and we define the discrete finite volume Laplacian ∆ D u ∈ H D (Ω) (corresponding to the homogeneous Dirichlet boundary conditions), such that for all K ∈ M, the constant value ∆ K u of ∆ D u in K is given by:

∆ K u = 1 m K L∈NK m K|L d K|L (u L -u K ) + σ∈EK ∩Eext m σ d K,σ (0 -u K ) . (2)
Note that this approximation is indeed obtained by a finite volume principle; in fact, for a continuous function u ∈ C 2 (Ω), integrating ∆u over a control volume K and using the Gauss divergence theorem and decomposing the boundary ∂K of K into its edges (∂K = ∪ σ∈EK σ) yields:

K ∆u(x)dx = ∂K ∇u(x) • n K (x)dγ(x) = σ∈EK σ ∇u(x) • n K,σ dγ(x),
where n K (resp. n K,σ ) denotes the outward unit normal vector to ∂K (resp. unit normal vector to σ outward to K) and dγ denotes the integration sign on the boundary. Under the orthogonality condition on the mesh, the quotients u(xL)-u(xK ) d K|L and 0-u(xK ) dK,σ are consistent finite difference approximations of ∇u(x) • n K,σ on σ = K|L and σ ∈ E K ∩ E ext . This motivates the choice of the approximate Laplacian [START_REF] Armfield | Ellipticity, accuracy and convergence of the discrete Navier-Stokes equations[END_REF].

Following the continuous framework, we define the following inner product of two elements v, w of the space H D (Ω):

[v, w] D = - Ω w(x)∆ D v(x)dx = - Ω v(x)∆ D w(x)dx = 1 2 K∈M L∈NK m K|L d K|L (v L -v K )(w L -w K ) + K∈M σ∈EK ∩Eext m σ d K,σ v K w K . (3)
The corresponding norm of w ∈ H D (Ω), denoted by

w D = ([w, w] D ) 1/2 ,
satisfies the discrete Poincaré inequality (see e.g. [9]), given by:

w L 2 (Ω) ≤ diam(Ω) w D , ∀w ∈ H D (Ω). ( 4
)
For any continuous function ϕ, we define the interpolation P D ϕ ∈ H D (Ω) by

(P D ϕ) K = ϕ(x K ), for all K ∈ M. (5)
Similarly, for u = (u

(i) ) i=1,...,d ∈ (H D (Ω)) d , v = (v (i) ) i=1,...,d ∈ (H D (Ω)) d and w = (w (i) ) i=1,...,d ∈ (H D (Ω)) d , we define: u D = d i=1 [u (i) , u (i) ] D 1/2 , [v, w] D = d i=1 [v (i) , w (i) ] D ,
and, for all vector field ϕ given by d continuous components, we define

P D ϕ ∈ H D (Ω) d by (P D ϕ) K = ϕ(x K ), for all K ∈ M.
Next, we wish to define a discrete divergence operator which approximates the continuous divergence term ∇ • (z u), where z is a scalar function from Ω to R (say a concentration or a temperature) and u is a vector function, from Ω to R d (typically a velocity). Again using the Gauss divergence theorem, we get that

K div(z u) = ∂K z u • n K dγ(x) = σ∈EK σ z u • n K,σ dγ(x).
Assuming the homogeneous Dirichlet boundary condition for u, a natural (centred) finite volume approximation of the divergence of z u is therefore div 0 D (z, u) ∈ H D (Ω), such that for all K ∈ M, its constant value div 0 K (z, u) on K is given by:

div 0 K (z, u) = 1 m K L∈NK m K|L d K|L n KL • (d L,K|L u K + d K,K|L u L ) z K + z L 2 , (6)
where z ∈ H D (Ω) and u ∈ H D (Ω) d are the sought approximations of z and u. However, since the approximations of the velocity and the pressure fields are collocated, we are led to define a modified discrete divergence operator which involves a family of stabilization parameters, denoted by λ, and the pressure field. Hence, for a velocity field u ∈ H D (Ω) d , for a pressure field p ∈ H D (Ω), for a scalar field z ∈ H D (Ω), we define the discrete divergence operator div D (z, u, λ, p) ∈ H D (Ω), such that for all K ∈ M, its constant value div K (z, u, λ, p) on K is given by:

div K (z, u, λ, p) = 1 m K L∈NK Φ KL (u, λ, p) z K + z L 2 , with Φ KL (u, λ, p) = m K|L d K|L n KL • (d L,K|L u K + d K,K|L u L ) -λ K|L (p L -p K ) , ∀L ∈ N K , (7)
where λ = (λ σ ) σ∈Eint is a given family of nonnegative real numbers. Our first choice of stabilization consists in choosing λ as

λ σ = λh α D , ∀σ ∈ E int , (8)
for a given λ > 0 and α ∈ (0, 2); this yields a stabilization term à la Brezzi-Pitkäranta [START_REF] Brezzi | On the stabilization of finite element approximations of the Stokes equations[END_REF]. Our second choice of stabilization consists in choosing a strictly positive value λ and introducing a partition of M, denoted by C; the elements of C, called clusters in this paper, are subsets of M: examples are shown in Figures 3 and11. We then set

λ K|L = λ if there exists G ∈ C with {K, L} ⊂ G, λ K|L = 0 otherwise. (9)
Several algorithms can be used for partitioning M. One can initialize the clusters by all the control volumes neighboring some initial given control volumes as in the example presented in Figure 3, or by all control volumes having a common vertex as in the example presented in Figure 11. Then the isolated remaining control volumes are included in the cluster which already contains the greatest number of its neighbors. Note that in [START_REF] Chénier | Numerical results using a collocated finite-volume scheme on unstructured grids for incompressible flows[END_REF] as well as in [START_REF] Brezzi | On the stabilization of finite element approximations of the Stokes equations[END_REF], the velocity field between the two control volumes K and L is interpolated in order to ensure an optimal error estimate (see [START_REF] Eymard | On a stabilized collocated Finite Volume scheme for the Stokes problem[END_REF]). For all v ∈ H D (Ω) d , the approximation div D (v, u, λ, p) of div(v ⊗ ū) is obtained by replacing in [START_REF] Chénier | Numerical results using a collocated finite-volume scheme on unstructured grids for incompressible flows[END_REF] the scalars z K and z L by the vectors v K and v L , that is:

div K (v, u, λ, p) = 1 m K L∈NK Φ KL (u, λ, p) v K + v L 2 . (10)
Finally, let us define a discrete gradient operator as the adjoint of the non stabilized divergence operator given in [START_REF] Brezzi | On the stabilization of finite element approximations of the Stokes equations[END_REF]; for any discrete field p ∈ H D (Ω), we define its discrete gradient

∇ D p ∈ H D (Ω) d : Ω ∇ D p(x) • v(x)dx = - Ω p(x)div 0 D (1, v)(x)dx, ∀v ∈ H D (Ω) d .
This implies the following definition for the constant value ∇ K p of ∇ D p in K ∈ M:

∇ K p = 1 m K L∈NK d L,K|L d K|L m K|L n KL (p L -p K ). (11)
Defining

p σ = d K,K|L d K|L p K + d L,K|L
d K|L p L if σ = K|L (note that this interpolation is not natural, contrary to that of (6) or ( 7)), and p σ = p K if σ ∈ E ext ∩ E K , and using the fact that σ∈EK m σ n K,σ = 0, one notices that ∇ K p may also be written

∇ K p = 1 m K σ∈EK m σ p σ n K,σ . (12)
Using the above definitions, we are now able to write the scheme for the linear Stokes equations in the next section. The extension to a nonlinear case is done in section 3.

The finite volume scheme in the linear case

Consider the generalized stationary Stokes equations on Ω ⊂ R d obtained after a time discretization of the unsteady Stokes problem:

η ū -ν∆ ū + ∇p = f , div ū = 0. ( 13
)
where ν is the kinematic viscosity, η is a frequency depending on the time scheme, ū = (ū i ) i=1,...,d denotes the velocity field, d is the space dimension, and p the pressure field, with the following assumptions:

Ω is a polygonal open bounded connected subset of R d , d = 2 or 3 (14) ν ∈ (0, +∞), η ∈ [0, +∞), f ∈ L 2 (Ω) d . (15)
Assuming the homogeneous Dirichlet boundary condition for ū on ∂Ω, the weak formulation of ( 13) is given by:

             ū ∈ H 1 0 (Ω) d , p ∈ L 2 (Ω) with Ω p(x)dx = 0, η Ω ū(x) • v(x)dx + ν Ω ∇ ū(x) : ∇v(x)dx - Ω p(x)div v(x)dx = Ω f (x) • v(x)dx, ∀v ∈ H 1 0 (Ω) d , div ū(x) = 0 for a.e. x ∈ Ω. ( 16
)
Let D be an admissible discretization of Ω. For a given family of nonnegative real numbers λ = (λ σ ) σ∈Eint , we look for (u, p) such that

                   (u, p) ∈ H D (Ω) d × H D (Ω) with Ω p(x)dx = 0, η Ω u(x) • v(x)dx + ν [u, v] D - Ω p(x)div 0 D (1, v)(x)dx = Ω f (x) • v(x)dx, ∀v ∈ H D (Ω) d , div D (1, u, λ, p)(x) = 0, for a.e. x ∈ Ω. ( 17
)
System ( 17) is equivalent to finding the family of vectors (u K ) K∈M ⊂ R d , and scalars (p K ) K∈M ⊂ R solution of the system of equations obtained by writing for each control volume

K of M:      ηu K -ν∆ K u + ∇ K p = 1 m K K f (x)dx div K (1, u, λ, p) = 0, (18)
supplemented by the relation

K∈M m K p K = 0. ( 19
)
Thanks to definition (2) and to the relation [START_REF] Eymard | A staggered finite volume scheme on general meshes for the Navier-Stokes equations in two space dimensions[END_REF], we see that [START_REF] Glowinski | Numerical Methods for fluids (Part3)[END_REF] can be written in a conservative form, that indicates this scheme is indeed a finite volume scheme. We can then prove (see [START_REF] Eymard | Convergence analysis of a collocated finite volume scheme for the incompressible Navier-Stokes equations on general 2 or 3D meshes[END_REF] for details) that if (u, p) ∈ H D (Ω) d × H D (Ω) is a solution to [START_REF] Girault | Finite element methods for the Navier-Stokes equations: Theory and algorithms[END_REF], then the following inequalities hold:

ν u D ≤ diam(Ω) f (L 2 (Ω)) d , (20)
and

1 2 K∈M L∈NK m K|L d K|L λ K|L (p L -p K ) 2 ≤ diam(Ω) ν f 2 (L 2 (Ω)) d . (21)
We now give some consequences of these a priori estimates, restricting our study to the case where there exists some λ > 0 and α ∈ (0, 2) such that ( 8) holds (note that under some geometrical assumptions on the mesh, convergence may also be proved for the stabilization by clusters, see [START_REF] Eymard | Numerical analysis of a clustered collocated finite volume scheme for the incompressible Stokes Navier-Stokes equations on general 2 or 3D meshes[END_REF] and [START_REF] Eymard | On the stability of colocated clustered finite volume simplicial discretizations for the 2D Stokes problem[END_REF]). We then get the uniqueness of a discrete solution to [START_REF] Girault | Finite element methods for the Navier-Stokes equations: Theory and algorithms[END_REF], and, using the Nečas technique, we can prove the following estimate on the pressure:

p L 2 (Ω) ≤ C 1 f (L 2 (Ω)) d , (22)
where C 1 only depends on d, Ω, η, ν, λ, α and on any θ > 0 such that regul(D) ≥ θ, and not on h D . From these estimates, and thanks to some estimates on the translates which may be obtained in a similar way as in [9], we obtain compactness properties of sequences of approximate solutions, as the mesh size tends to 0. Passing to the limit in the scheme, we then get the following convergence result [START_REF] Eymard | Convergence analysis of a collocated finite volume scheme for the incompressible Navier-Stokes equations on general 2 or 3D meshes[END_REF]:

Theorem 2.1 (Convergence) Under hypotheses ( 14)- [START_REF] Eymard | Numerical analysis of a clustered collocated finite volume scheme for the incompressible Stokes Navier-Stokes equations on general 2 or 3D meshes[END_REF], let ( ū, p) be the unique weak solution of the Stokes problem [START_REF] Eymard | Convergence analysis of a collocated finite volume scheme for the incompressible Navier-Stokes equations on general 2 or 3D meshes[END_REF]. Let λ ∈ (0, +∞), α ∈ (0, 2) and θ > 0 be given and let D be an admissible discretization of Ω such that regul(D) ≥ θ (see definitions in section 2.1). Let (u D , p D ) ∈ H D (Ω) d × H D (Ω) be the unique solution to [START_REF] Girault | Finite element methods for the Navier-Stokes equations: Theory and algorithms[END_REF] with λ defined by [START_REF] Chiang | Topological Flow Structures in Backward-Facing Step Channels[END_REF]. Then u D converges to ū in L 2 (Ω) d and p D weakly converges to p in L 2 (Ω) as h D tends to 0.

Note that the proof of the strong convergence of p D to p is also a straightforward consequence of the error estimate which we now state (also given in [START_REF] Eymard | Convergence analysis of a collocated finite volume scheme for the incompressible Navier-Stokes equations on general 2 or 3D meshes[END_REF]), under additional regularity hypotheses.

Theorem 2.2 (Error estimate) Under hypotheses ( 14)-( 15), we assume that the weak solution ( ū, p) of the Stokes problem [START_REF] Eymard | Convergence analysis of a collocated finite volume scheme for the incompressible Navier-Stokes equations on general 2 or 3D meshes[END_REF] i.e. satisfying ( 16) is such that ( ū, p) ∈ H 2 (Ω) d × H 1 (Ω). Let λ ∈ (0, +∞) and α ∈ (0, 2) be given, let D be an admissible discretization of Ω and let θ > 0 such that regul(D (m) ) ≥ θ (see definitions in section 2.1). Let (u D , p D ) ∈ H D (Ω) d × H D (Ω) be the solution to [START_REF] Girault | Finite element methods for the Navier-Stokes equations: Theory and algorithms[END_REF] with λ defined by [START_REF] Chiang | Topological Flow Structures in Backward-Facing Step Channels[END_REF]. Then there exists C 2 , which only depends on d, Ω, ν, η and θ such that

u D -ū 2 L 2 (Ω) ≤ C 2 ǫ(λ, h D , p, ū), (23)
λ h α D |p D | 2 D ≤ C 2 ǫ(λ, h D , p, ū) (24)
p D -p 2 L 2 (Ω) ≤ C 2 ǫ(λ, h D , p, ū). ( 25
)
where

ǫ(λ, h D , p, ū) = max λh α D , h 2-α D λ p 2 H 1 (Ω) + ū 2 H 2 (Ω) . (26)
Hence, for α = 1 we get an order 1/2 for the convergence of the scheme. In fact, this result is not sharp, and the numerical results show a much better order of convergence.

The Navier-Stokes equations under the Boussinesq approximation

Let us consider in this section the case of the steady equations obtained after a time discretization of the Navier-Stokes and energy equations, under the Boussinesq hypotheses

η ū -ν∆ ū + ∇p + div( ū ⊗ ū) = f + T w in Ω, div ū = 0 in Ω, η T -κ∆ T + div( T ū) = g in Ω. (27)
with homogeneous Dirichlet boundary conditions on ū and T , assumptions ( 14)-( 15) on Ω, ν, η and f and with κ the thermal diffusivity, w and g such that:

κ ∈ (0, +∞), w ∈ R d , g ∈ L 2 (Ω). ( 28
)
Then we say that ( ū, p, T ) is a weak solution of [START_REF] Nicolaïdes | Analysis and Convergence of the MAC Scheme I. The Linear Problem[END_REF] if

                                 ū ∈ H 1 0 (Ω) d , p ∈ L 2 (Ω) with Ω p(x)dx = 0, T ∈ H 1 0 (Ω), Ω (η ū(x) + div( ū ⊗ ū)(x)) • v(x)dx + ν Ω ∇ ū(x) : ∇v(x)dx - Ω p(x)divv(x)dx = Ω (f (x) + T (x)w) • v(x)dx, ∀v ∈ H 1 0 (Ω) d , Ω (η T (x) + div( T ū)(x)) S(x)dx + κ Ω ∇ T (x) • ∇ S(x)dx = Ω g(x) S(x)dx, ∀ S ∈ H 1 0 (Ω), div ū(x) = 0 for a.e. x ∈ Ω. ( 29
)
We now give the finite volume scheme for this problem. Under hypotheses ( 14)-( 15) and ( 28), let D be an admissible discretization of Ω. For a given family of nonnegative real numbers λ = (λ σ ) σ∈Eint , we look for (u, p, T ) such that

                               (u, p, T ) ∈ H D (Ω) d × H D (Ω) × H D (Ω) with Ω p(x)dx = 0, ν [u, v] D + Ω (ηu(x) + div D (u, u, λ, p)(x)) • v(x)dx - Ω p(x)div 0 D (1, v)(x)dx = Ω (f (x) + T (x)w) • v(x)dx, ∀v ∈ H D (Ω) d div D (1, u, λ, p) = 0, κ[T, S] D + Ω (ηT (x) + div D (T, u, λ, p)(x))S(x)dx = Ω g(x)S(x)dx, ∀S ∈ H D (Ω). ( 30
)
Note that the relation div

D (1, u, λ, p) = 0 implies Ω v(x) • div D (u, u, λ, p)(x)dx = - Ω u(x) • div D (v, u, λ, p)(x)dx, Ω S(x)div D (T, u, λ, p)(x)dx = - Ω T (x)div D (S, u, λ, p)(x)dx. ( 31
)
which shows that the discrete versions of the operators behave as the continuous ones. Again, system ( 30) is equivalent to finding the family of vectors (u K ) K∈M ⊂ R d , and scalars (p K , T K ) K∈M ⊂ R such that, for each control volume K of M:

             ηu K -ν∆ K u + div K (u, u, λ, p) + ∇ K p = 1 m K K f (x)dx + T K w, div K (1, u, λ, p) = 0, ηT K -κ∆ K T + div K (T, u, λ, p) = 1 m K K g(x)dx, (32)
supplemented by the relation:

K∈M m K p K = 0.
One should notice that, thanks to [START_REF] Tai | Parallel computation of unsteady three-dimensional incompressible viscous flow using an unstructured multigrid method[END_REF], one has two ways of implementing the scheme [START_REF] Touazi | Etude numérique d'un schéma de volumes finis sur des maillages non structurés pour des écoulements de fluide visqueux à densité constante ou faiblement variable[END_REF]. These two ways are not equivalent from the point of view of the convergence properties of the numerical algorithms described in Section 4.1.

Theorem 3.1 (Convergence of the scheme) Under hypotheses ( 14)-( 15) and [START_REF] Patankar | Numerical Heat Transfer and Fluid Flow[END_REF], let (D (m) ) m∈N be a sequence of admissible discretizations of Ω such that h D m tends to 0 as m → ∞ and such that there exists θ > 0 with regul(D (m) ) ≥ θ, for all m ∈ N (see definitions in section 2.1). Let λ ∈ (0, +∞) and α ∈ (0, 2) be given. Then, for all m ∈ N, there exists at least one solution to [START_REF] Peyret | Computational methods for fluid flow[END_REF] with D = D (m) and λ = λ (m) defined by ( 8), that we denote

(u (m) , p (m) , T (m) ) ∈ (H D (m) (Ω)) d × H D (m) (Ω) × H D (m) (Ω).
Moreover, there exists ( ū, p, T ) satisfying ( 29) and a subsequence of (D (m) ) m∈N , again denoted (D (m) ) m∈N , such that the corresponding subsequence of solutions (u (m) , T (m) ) m∈N converges to ( ū, T ) in L 2 (Ω) and (p (m) ) m∈N weakly converges to p in L 2 (Ω).

Proof Several keypoints of the proof of this convergence result may be found in [START_REF] Eymard | Convergence analysis of a collocated finite volume scheme for the incompressible Navier-Stokes equations on general 2 or 3D meshes[END_REF], where we showed the convergence of the scheme for the isothermal Navier-Stokes equations. The main differences with [START_REF] Eymard | Convergence analysis of a collocated finite volume scheme for the incompressible Navier-Stokes equations on general 2 or 3D meshes[END_REF] are the fact that we have to handle the additional convection-diffusion equation on the temperature, and that we do not use the Bernoulli pressure thanks to the definition (7) of the discrete divergence, which includes the stabilization pressure term. Let us then notice that, thanks to (31), we get

Ω u(x) • div D (u, u, λ, p)(x)dx = 0 and Ω T (x)div D (T, u, λ, p)(x)dx = 0. ( 33
)
Hence, taking u = v in the first equation of ( 30) and S = T in the last one, we get the estimate

κ T D ≤ diam(Ω) g L 2 (Ω) , (34)
and therefore

ν u D ≤ diam(Ω) f L 2 (Ω) d + |w| diam(Ω) κ g L 2 (Ω) , (35)
and, thanks to [START_REF] Chiang | Topological Flow Structures in Backward-Facing Step Channels[END_REF],

λ h α D 1 2 K∈M L∈NK m K|L d K|L (p L -p K ) 2 ≤ diam(Ω) ν f (L 2 (Ω)) d + |w| diam(Ω) κ g L 2 (Ω) 2 . ( 36
)
We then deduce the existence of C 1 , only depending on d, Ω, η, ν, λ, κ, w and θ, and not on h D , such that

p L 2 (Ω) ≤ C 1 f L 2 (Ω) d + g L 2 (Ω) + f L 2 (Ω) d + g L 2 (Ω) 2 (37)
in a similar way as in [START_REF] Eymard | Convergence analysis of a collocated finite volume scheme for the incompressible Navier-Stokes equations on general 2 or 3D meshes[END_REF]. These estimates are sufficient to allow the application of the topological degree method (see theorem 4.3 of [START_REF] Eymard | Convergence analysis of a collocated finite volume scheme for the incompressible Navier-Stokes equations on general 2 or 3D meshes[END_REF]). Hence, for all m ∈ N, we prove the existence of at least one solution (u (m) , p (m) , T [START_REF] Peyret | Computational methods for fluid flow[END_REF] with D = D (m) and λ = λ (m) defined by [START_REF] Chiang | Topological Flow Structures in Backward-Facing Step Channels[END_REF]. These estimates also prove the existence of ū ∈ H 1 0 (Ω) d , p ∈ L 2 (Ω) with Ω p(x)dx = 0, of T ∈ H 1 0 (Ω), and of a subsequence of (D (m) ) m∈N , again denoted (D (m) ) m∈N , such that (u (m) ) m∈N converges in L 2 (Ω) to ū, (p (m) ) m∈N weakly converges in L 2 (Ω) to p and (T (m) ) m∈N converges in L 2 (Ω) to T . It remains now to prove that ( ū, p, T ) satisfies [START_REF] Petrone | Stability of free convection in air-filled horizontal annuli: Influence of the radius ratio[END_REF]. Note that div ū(x) = 0 for a.e. x ∈ Ω is obtained in the same way as in [START_REF] Eymard | Convergence analysis of a collocated finite volume scheme for the incompressible Navier-Stokes equations on general 2 or 3D meshes[END_REF]. The only new points to check concern the right hand side f + T w instead of f , the additional equation and the nonlinear terms. Since (T (m) ) m∈N converges in L 2 (Ω) to T , the convergence to f + T w of its discrete counterpart is clear. Let ϕ ∈ C ∞ c (Ω) d be given. Let us prove that

(m) ) ∈ H D (m) (Ω) d × H D (m) (Ω) × H D (m) (Ω) to
Ω P D (m) ϕ(x) • div D (u (m) , u (m) , λ (m) , p (m) )(x)dx converges to Ω ϕ(x) • div( ū ⊗ ū)(x)dx as m → ∞.
For a given discretisation D and for (u, p) where we use the notation Φ KL defined in [START_REF] Chénier | Numerical results using a collocated finite-volume scheme on unstructured grids for incompressible flows[END_REF], and where the function G D (p, ϕ) : Ω → R d is defined, for all K ∈ M and σ ∈ E K , by the constant value G K,σ (p, ϕ) in the cone D K,σ with basis σ and vertex x K , such that

∈ H D (Ω) d × H D (Ω), defining λ by (8), let us write Ω P D ϕ(x) • div D (u, u, λ, p)(x)dx = T 1 (D, u, p) + T 2 (D, u, p) with T 1 (D, u, p) = 1 2 K∈M L∈NK Φ KL (u, 0, 0) u L + u K 2 • (ϕ(x K ) -ϕ(x L )) , T 2 (D, u, p) = Ω u(x) • G D (p, ϕ)(x)dx,
G K,K|L (p, ϕ) = 1 2 λh α D m(D K,K|L ) m K|L d K|L (p L -p K )(ϕ(x L ) -ϕ(x K )), ∀L ∈ N K , and G K,σ (p, ϕ) = 0, ∀σ ∈ E K ∩ E ext .
We remark that, thanks to (36), G D (m) (p (m) , ϕ) converges to 0 in L 2 (Ω) d as as m → +∞. Hence we get that lim

m→∞ T 2 (D (m) , u (m) , p (m) ) = 0.
In the same way as in [START_REF] Eymard | Convergence analysis of a collocated finite volume scheme for the incompressible Navier-Stokes equations on general 2 or 3D meshes[END_REF], we have

lim m→∞ T 1 (D (m) , u (m) , p (m) ) = Ω ϕ(x) • div( ū ⊗ ū)(x)dx.
Hence, the expression [START_REF] Eymard | Convergence analysis of a collocated finite volume scheme for the incompressible Navier-Stokes equations on general 2 or 3D meshes[END_REF], and the convergence study of the corresponding right hand side is straightforward. This completes the proof of the theorem.

Ω P D (m) ϕ(x)•div D (u (m) , u (m) , λ (m) , p (m) )(x)dx converges to Ω ϕ(x)• div( ū ⊗ ū)(x))dx as m → ∞. For a given function ϕ ∈ C ∞ c (Ω), the convergence proof of Ω P D (m) ϕ(x)div D (T (m) , u (m) , λ (m) , p (m) )(x)dx to Ω ϕ(x)div( T ū)(x))dx as m → ∞ follows similar arguments. The convergence of [T (m) , P D (m) ϕ] D (m) to Ω ∇ T • ∇ϕ is proven in Proposition 2.4 of
Remark 1. Note that it is possible to prove a convergence theorem, similar to Theorem 3.1, in the case of nonhomogeneous Dirichlet boundary and g = 0. We then consider the upstream approximation of div( T ū), given by

div ups K (T, u, λ, p) = 1 m K L∈NK Φ KL (u, λ, p) + T K -Φ KL (u, λ, p) -T L ,
where Φ KL (u, λ, p) is defined in [START_REF] Chénier | Numerical results using a collocated finite-volume scheme on unstructured grids for incompressible flows[END_REF], and where, for all s ∈ R, we set s + = max(s, 0) and s -= max(-s, 0). In this case, it is possible to derive an L ∞ estimate on T , which allows for similar inequalities to (35) and (36), and then it is easy to conclude the convergence proof, following the methods used in [9] in the case of the approximation by the finite volume method of a nonhomogeneous Dirichlet elliptic problem. Note that the convergence of the scheme, keeping the centered approximation (7) of div K (T, u, λ, p) instead using div ups K (T, u, λ, p), remains an open problem, because the lack of a bound on the velocities prevents from getting the stability of the resulting centred convection diffusion scheme. Nevertheless, the centred approximation remains stable in many cases, including the numerical tests given in section 4.5.

Numerical results

Practical implementation

Let us first comment the nature of the different problems to be solved. The scheme ( 18) is a coupled linear problem with a symmetric matrix, but the scheme (32) is a coupled nonlinear problem, the Jacobian matrix of which is nonsymmetric (taking the opposite of the fluid mass conservation equation). Our numerical choices are as follows:

1. The nonlinear problem [START_REF] Touazi | Etude numérique d'un schéma de volumes finis sur des maillages non structurés pour des écoulements de fluide visqueux à densité constante ou faiblement variable[END_REF] is solved by a coupled under-relaxed Newton iteration method (in order to avoid large variations of the unknown in one iteration which could lead to the divergence of the method). This implies that at each iteration, a linear system involving the Jacobian matrix of the full nonlinear system, must be solved, with 4 unknowns per control volume in 2D and 5 in 3D. Several advantages are gained from this procedure: there is no approximate Poisson problem to solve in order to compute the pressure (it is well known that such a procedure induces the appearance of boundary layers), the method is completely systematic and the programming easy; furthermore steady problems may be solved directly, without using artificial transient problems.

2. The coupled linear systems are solved in our prototype, using either a direct solver (Gaussian elimination without pivoting, for small tests) or a GMRES or BiCGSTAB method for large ones with a incomplete LU preconditioning. These methods have been shown to be efficient in all the test cases which are presented in this paper.

Another advantage of the numerical method that we have used is that the accuracy which is requested in the linear system resolution can be linked with the residual at each Newton iteration, which avoids wasting time in precise linear resolution when the Newton iteration is still far for convergence.

Study of the stabilization methods

We implemented the two stabilization techniques ( 8) and ( 9) for the steady and transient Navier-Stokes equations with or without the energy equation [START_REF] Chénier | Numerical results using a collocated finite-volume scheme on unstructured grids for incompressible flows[END_REF]. Let us first show on an example the importance of introducing a stabilization term in the numerical scheme. We consider the case Ω = (0, 1) × (0, 1), η = 0, ν = 1 and g = 0 (it corresponds to an isothermal case), and f chosen such that (27) holds with T (x, y) = 0, ū(1) (x, y) = -∂ y Φ(x, y) and ū(2) (x, y) = ∂ x Φ(x, y), computed from the potential Φ(x, y) = 600[x(1x)y(1y)] 2 , the pressure being given by p(x, y) = 100(x 2 + y 2 -2 3 ), for all (x, y) ∈ Ω. We then obtain that the numerical solution obtained for the velocity is accurate whatever the method and the magnitude of stabilization. On the contrary, it is clear from Figure 2 that the pressure field suffers from severe oscillations when λ = 0 or is close to 0; these results are obtained using the cluster stabilization (9) (the clusters are shown on Figure 3), but a similar behaviour on the pressure is obtained with the Brezzi Pitkäranta stabilization.

Let us compare the two stabilization terms respectively given by ( 8) and ( 9). An advantage of ( 8) is that it leads to an easier mathematical analysis of the convergence of the scheme, but a drawback is that it yields some redistribution of the fluid mass over the whole domain. Moreover, in order to obtain convergence, one needs to let the stabilization parameter tend to 0 with the size of the mesh h D (hence the expression of the stabilization term in ( 8)). On the contrary, in the stabilization term given by ( 9), there is no need to let λ tend to 0 with the size of the mesh, which means that the presence of a finite stabilization does not decrease the quality of the approximation. Let us show this property on the backward facing step example at Reynolds number 800 (the complete data are presented, for example, in [START_REF] Chiang | Topological Flow Structures in Backward-Facing Step Channels[END_REF]). We show in Figure 4 the computed isovalue lines of the horizontal component of the velocity u (1) , obtained for a mesh of 6300 cells. Roughly speaking, the isovalue lines u (1) = 0 (dashed in Figure 4) separates the two recirculation vortices; the determination of this isoline therefore allows the evaluation of the reattachment lengths of both vortices. We then study the effect of the stabilization term on these reattachment lengths. In order to do so, we consider a given mesh and decrease the stabilization term. This effect is shown in Figure 5, where we draw the values of both reattachment lengths computed with the same mesh, as a function of the value of the "weight" of the stabilization, that is the global value of λh α D in [START_REF] Chiang | Topological Flow Structures in Backward-Facing Step Channels[END_REF], and λ in (9). We see that, as the stabilization value increases, the reattachment lengths obtained with the Brezzi-Pitkäranta stabilization decrease much more than those given by the cluster stabilization. This shows that, in the case of problems requesting a high stabilization level, the numerical results expected using the cluster stabilization remains more accurate. The adaptation of the convergence proofs to the cluster stabilization is the object of recent works [START_REF] Eymard | Numerical analysis of a clustered collocated finite volume scheme for the incompressible Stokes Navier-Stokes equations on general 2 or 3D meshes[END_REF], [START_REF] Eymard | On the stability of colocated clustered finite volume simplicial discretizations for the 2D Stokes problem[END_REF].

Numerical confirmation of the convergence study

Let us first emphasize that the order of convergence of a numerical scheme must be distinguished from the consistency order, obtained by introducing in the numerical scheme an analytical solution taken in some points: we noticed above that for example, we do not require that the approximate gradient of pressure be consistent. Indeed, the definition ( 11) is not consistent in the general case and we only require that it is the transpose of the discrete divergence approximate. Moreover, as mentioned above, the orders of convergence that were mathematically proven (see Theorem 2.2) are hopefully much lower that the observed numerical orders of convergence. Indeed, let us first consider the well-known transient Green-Taylor analytical solution (isothermal case), in the square domain Ω = (-0.5, 0.5) × (-0.5, 0.5), given by ū(1) (x, y, t) =cos(πx) sin(πy) exp(-2π 2 t/Re) ū(2) (x, y, t) = sin(πx) cos(πy) exp(-2π 2 t/Re)

p(x, y, t) = - cos(2πx) + cos(2πy) 4 exp(-4π 2 t/Re),
where Re is the Reynolds number. The initial and boundary conditions are deduced from the exact solution. We consider the case Re = 10, and we compare the approximate solutions given by the collocated scheme using the cluster stabilization (9), discretized by the Crank-Nicholson method, with the analytical ones at time t = 0.3, for rectangular meshes with respectively 400, 1600, 6400 and 25 600 control volumes, keeping the time step equal to 0.001. We then plot in Figure 6 the L 2 norms of the errors e(u) and e(p) between the computational and the analytical solution, for these different meshes. We observe that the order of convergence, given by the slope of fitted lines on these curves, is close to 2 for the velocity, and 1 for the pressure.

We study in a similar way the error as a function of the time steps, on a 160 × 160 rectangular mesh. The L 2 norms of the errors e(u) and e(p) are plotted in Figure 7, showing similar order of convergence as in space. -

-2 -1.9 -1.8 -1.7 -1.6 -1.5 -1.4 -1.3 -1.2 -1.1 -1 -0.9 log Since there is no straightforward analytical solution with null right-hand-side in the case of non isothermal problem, we considered a steady problem (η = 0) where a divergence free analytical solution is arbitrarily chosen, and the right-hand-sides f and g are the residual of the continuous Boussinesq and energy equations without source terms. Hence we define the scalar potential ψ as in [START_REF] Touazi | Etude numérique d'un schéma de volumes finis sur des maillages non structurés pour des écoulements de fluide visqueux à densité constante ou faiblement variable[END_REF], by ψ(x, y) = sin 2 (πx) sin 2 (πy), the temperature and the pressure fields are respectively given by T (x, y) = sin 2 (πx) sin 2 (πy) p(x, y) = sin 2 (πx) sin 2 (πy) and the velocity field is given by ū(1) (x, y) = ∂Ψ ∂y = 2π cos(πy) sin(πy) sin 2 (πx)

-2 -1.9 -1.8 -1.7 -1.6 -1.5 -1.4 -1.3 -1.2 -1.1 -1 -0.9 log 10 (||e(u)||) L2 log 10 (∆t) 4.2 -4 -3.8 -3.6 -3.4 -3.2 -3 -2.8 -2.6 -2.4 -2.2
-2 -1.9 -1.8 -1.7 -1.6 -1.5 -1.4 -1.3 -1.2 -1.1 -1 -0. -1.2 -1.1 -1 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3
ū(2) (x, y) = - ∂Ψ ∂x = -2π cos(πx) sin(πx) sin 2 (πy).
Choosing the cluster stabilization (9), we set λ = 0.001, ν = κ = 1, w is the unit vertical vertor pointing upward. We proceed in the same way as we did for the Green-Taylor analytical solution, computing the order of convergence for variable space steps. We found an order 2 convergence for the velocity, the temperature, and an order 1 for the pressure in the case of rectangular meshes or Delaunay meshes (see Figure 8 for example of such a mesh).

Comparison with other methods

Let us first again consider the same backward facing step problem as the one used in Section 4.2 for the study of the effects of the stabilization method. We apply the collocated scheme, using a mesh with 109 050 triangles. In table I, we provide the first and second reattachment lengths in the case Re= 800 given in the literature, together with ours. We observe that our results show an excellent agreement with that obtained by other methods.

Let us now turn to the 3D lid driven cavity test case. The considered domain is the unit cube, and we consider two values for the Reynolds number (Re= 400 and Re= 3200). We compare the profiles of some components of the velocity along median lines with that given in [START_REF] Tai | Parallel computation of unsteady three-dimensional incompressible viscous flow using an unstructured multigrid method[END_REF] with 97 336 control volumes (finite volume method on unstructured meshes). We use 1st vortex 2nd vortex Armaly et al. (experimental results) [START_REF] Armaly | Experimental and theoretical investigation of backward-facing step flow[END_REF] 14.2 20.0 Chiang et al. [START_REF] Chiang | Topological Flow Structures in Backward-Facing Step Channels[END_REF] 12.3 20.8 Kim & Moin [START_REF] Kim | Application of a fractionnal-step method to incompressible Navier-Stokes equations[END_REF] 12.0 Kaiktsis et al. [START_REF] Kaiktsis | Onset of three-dimensionality, equilibria, and early transition in flow over a backward-facing step[END_REF] 11.9 Centred staggered scheme [START_REF] Eymard | A staggered finite volume scheme on general meshes for the Navier-Stokes equations in two space dimensions[END_REF] 13.3 20.6 Collocated scheme 12.0 19.9

Table I. First and second reattachment lengths for Re= 800 a simple 36 × 36 × 36 parallelepipedic mesh with refined sides along the boundaries of the domain. The stabilization parameter is taken equal to 10 -5 . We observe in Figure 9 that the results obtained by the collocated scheme and those which are given in the reference [START_REF] Tai | Parallel computation of unsteady three-dimensional incompressible viscous flow using an unstructured multigrid method[END_REF] are very close.

Natural convection between two cylinders

We consider the natural convective flows between two horizontal cylinders where the surface of the inner cylinder is heated. Such a configuration is, for instance, a model for thermal storage systems. An important industrial problem is therefore to control the heat transfer, either to enhance or to reduce it. Major published works have considered fluid flows and heat transfer between concentric cylinders [START_REF] Petrone | Stability of free convection in air-filled horizontal annuli: Influence of the radius ratio[END_REF]. A simple way to modify the thermal performances is to change the geometry and for instance by introducing an eccentricity in the annular space.

The geometrical characteristics of the cavity are defined by the inner C i and outer C o cylinders of radii r * i and r * o , by the eccentricity ǫ * and by the eccentricity angle φ ǫ (Fig. 10). In our study, uniform temperatures T * (C i ) = T * i and T * (C o ) = T * o are imposed on the no-slip cylinder walls (which means homogeneous Dirichlet boundary conditions for the velocity). The flow is assumed two-dimensional, steady, laminar and incompressible with constant physical properties except for the density when formulating the buoyancy effect. We then scale the lengths by the annular gap d * = r * or * i , the velocity components by the ratio between the thermal diffusivity and d * , we introduce the dimensionless difference temperature T = (T * -T

* r )/(T * i -T * o ) φǫ r * o r * i δ * ǫ * δ * Co Ci ξ ξ Figure 10. Geometry.
where T * r = (T * i + T * o )/2 is the reference temperature. The dimensionless Navier-Stokes and energy equations are then given by [START_REF] Nicolaïdes | Analysis and Convergence of the MAC Scheme I. The Linear Problem[END_REF] in which f = 0, η = 0, g = 0, ν = Pr, κ = 1 and w = Ra Pr e (2) , where Ra and Pr are the Rayleigh and Prandtl numbers and e (2) is the upward vertical unit vector. The boundary conditions for the temperature are then nonhomogeneous Dirichlet boundary conditions ( T = 1/2 on C i and T = -1/2 on C o ). We denote the radius ratio by R = r * o /r * i , the dimensionless eccentricity by ǫ = ǫ * /d * and the eccentricity angle by φ ǫ .

Voronoï tessellations are used to mesh the annular space. This choice meets the definition of an admissible mesh, since the straight line joining the centers of two adjacent cells is orthogonal to their common edge (see section 2.1). Since the cell centers are mainly located at the vertices of equilateral triangles, most of the resulting Voronoï cells have hexagonal shapes. To adjust the mesh to the natural convective flows, cells with a smaller hexagonal size were used along the cylinder walls, in a layer of thickness δ = δ * /d * , and in an angular sector 2ξ around the inner cylinder where a thermal plume and higher velocity occur (Fig. 10). The cells in the boundary layer and in the angular sector are twice smaller than the cells which define the basic hexagonal cell size. An enlargement of such a mesh is presented in figure 11a.

We use the cluster stabilization method (9), initializing the clusters by control volumes having one common vertex. Figure 11b illustrates the partition into clusters of the mesh shown in Figure 11a.

The heat transfer is measured by the mean Nusselt number at the inner cylinder

N u i = 1 2πr i Ci ∇ T (x) • n ∂Ω (x) ds(x),
where r i = r * i /d * and n ∂Ω (x) is the normal outward unit vector at point x ∈ C i . For coaxial cylinders and moderate Rayleigh number values, the flow field consists of two-dimensional crescent-shape cells, symmetrical with respect to the vertical straight line containing the centers of the cylinders. When the eccentricity is non zero and φ ǫ = 0 o , φ ǫ = 180 o , the vertical symmetry of annular space is broken and a fluid flow occurs around the inner cylinder. This flow rate is computed by ψ Ci with ψ the stream function defined by u The effects of the eccentricity and eccentricity angle are studied for a radius ratio R = 2, a Rayleigh number value Ra = 10 4 and a Prandtl number value Pr = 0.7. The angular sector and the thickness of the boundary layer, where the cell size is half the basic hexagonal cell size, are equal to 2ξ = 270 o and δ = 0.2. The basic hexagonal cell size 0.03 is classically adopted for the simulations.

In figures 12-14 are drawn the isotherms as a function of φ ǫ , for ǫ = 0, 0.5 and 0.95. Continuous and dashed lines correspond to positive and negative temperature values respectively. As foreseen, the isotherms are no longer symmetric when φ ǫ = 0 o or = 180 o and ǫ > 0. For φ ǫ ≤ 90 o , a thermal plume develops above the hot cylinder. On the other hand, the temperature field is rather convected from the lateral region of the inner cylinder when this latter one is located in the upper region of the outer cylinder, namely for φ ǫ > 90 o . When the cylinders are close, the Rayleigh number based on the smallest gap between the two walls is small, and a conduction thermal regime is therefore locally expected. Figure 15 illustrates the mean Nusselt number at the inner cylinder, N u i , as a function of φ ǫ and ǫ. The heat transfer is weakly affected by the eccentricity angle and the eccentricity if ǫ < 0.5. Conversely, the mean Nusselt number is doubled between ǫ = 0.5 and ǫ = 0.95. The symmetry breaking of the geometry with respect to the vertical line containing the center of the outer cylinder generates a fluid flow circulation around the inner cylinder for which the dimensionless flow rate is measured by ψ Ci (Fig. 16). The circulation is maximal for large eccentricities and when the inner cylinder is in the upper region of the outer cylinder. It is however clear that if both cylinders are too close, the fluid flows with difficulty around the hot cylinder and even cannot pass when they are side by side. In this limit case, the regularity of the boundary of the domain decreases, the boundary condition is no longer given by the trace of a regular function and the convergence of the scheme is far from being proven (the distance between the two cylinders should be discretized by a sufficient number layers of control volumes, which is maybe not the case in our computations near the limit case). As it can be observed in figures 15 and 16, some values seem to exhibit spurious variations as a function of φ ǫ and ǫ. For example ψ Ci decreases sharply in the vicinity of φ ǫ ≈ 140 o and ǫ ≈ 0.5. One could also wonder whether this behavior might be explained by the occurrence of bifurcations, and a stability analysis of the steady flow should be conducted to confirm the appearance of transitions in the fluid flow.

Conclusions and perspectives

Collocated schemes for the incompressible Navier-Stokes equations were analyzed and implemented on unstructured meshes. The collocated scheme, for which the proof of the convergence of the scheme has been completed in some cases, can handle some cases with a complex geometry, as the problem of natural convection between two cylinders with eccentricity, using unstructured meshes such as Voronoï tesselations or triangular meshes. However, the Brezzi-Pitkäranta stabilization method was found to be more sensitive to high Reynolds numbers on tests such as the backward step. The cluster stabilization was then implemented and improved dramatically the stability of the collocated scheme (in fact, clustering leads to large scale grids for the pressure). Work is in progress concerning the convergence analysis of the scheme with clustering, and of the viscous tensor term in the compressible case. The case of non constant viscosity also needs to be addressed.

Copyright c 2006

 2006 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Fluids 2006; 00:0-0 Prepared using fldauth.cls

Figure 2 .Figure 3 .

 23 Figure 2. Pressure field for different values of λ (from the left to the right: 0, 10 -6 , 10 -4 , 10 -2 )

Figure 4 .

 4 Figure 4. Isovalue lines of u(1) for the backward facing step (the step is visible at the lower left part of the figure) with Re= 800. The dashed lines correspond to the isovalue u (1) = 0 line.

Figure 5 .

 5 Figure 5. Reattachment lengths for different stabilization parameters. The dashed lines correspond to the Brezzi-Pitkäranta stabilization, the continuous ones to the cluster stabilization.

Figure 6 .

 6 Figure 6. L 2 norm of error (left: velocity, right: pressure) for different space steps

Figure 7 .

 7 Figure 7. L 2 norm of error (left: velocity, right: pressure) for different time steps

Figure 8 .

 8 Figure 8. Example of a Delaunay mesh

Figure 9 .

 9 Figure 9. 3D Lid driven cavity

  (1) = ∂ψ ∂y , u(2) = -∂ψ ∂x and ψ Co = 0. Copyright c 2006 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Fluids 2006; 00:0-0 Prepared using fldauth.cls (a) Mesh (b) clusters

Figure 11 .

 11 Figure 11. Voronoï tessellations based on hexagonal cells and partition of the computational domain into clusters.

Figure 12 .

 12 Figure 12. Isotherms of the steady solutions, ǫ = 0.

Figure 13 .

 13 Figure 13. Isotherms of the steady solutions, ǫ = 0.5 and φǫ = 0 o , 45 o , 90 o , 135 o , 180 o .

Figure 16 .

 16 Figure 16. Dimensionless flow rate around the inner cylinder, ψC i as a function of ǫ and φǫ.

Copyright c 2006 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Fluids 2006; 00:0-0 Prepared using fldauth.cls


Numerical results using a colocated finite-volume scheme on unstructured grids for incompressible fluid flows

Eric Chénier, R. Eymard, O.