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

Robert Eymard

Raphaèle Herbin

A collocated finite volume scheme to solve free convection for general non-conforming grids

Keywords: Collocated finite volume schemes, general non-conforming grids, Navier-Stokes equations, Boussinesq assumption PACS: 65C20, 76R10

We present a new collocated numerical scheme for the approximation of the Navier-Stokes and energy equations under the Boussinesq assumption for general grids, using the velocity-pressure unknowns. This scheme is based on a recent scheme for the diffusion terms. Stability properties are drawn from particular choices for the pressure gradient and the non-linear terms. Convergence of the approximate solutions may be proven mathematically. Numerical results show the accuracy of the scheme on irregular grids.

Introduction

Finite volume methods have been widely used in computational fluid dynamics for a long time; they are well adapted to the discretization of partial differential equations under conservative form; one of their attractive features is that the resulting discretized equation has a clear physical interpretation [START_REF] Patankar | Numerical Heat Transfer and Fluid Flow[END_REF]. In the framework of incompressible fluid flows, two strategies are often opposed, namely staggered and collocated schemes. The staggered strategy, which has become very popular since Patankar's book [START_REF] Patankar | Numerical Heat Transfer and Fluid Flow[END_REF], remains mainly restricted to geometrical domains with parallel and orthogonal boundary faces. Therefore, for computations on complex domains with general meshes, the collocated strategy which consists in approximating all unknowns on the same set of points (called collocation points but also cell centers or simply centers), is often preferred, even though the pressure-velocity coupling demands some cure for the stabilization of the well-known checkerboard pressure modes; to this purpose, various pressure stabilization procedures, based on improvements of the Momentum Interpolation Method proposed by Rhie and Chow [START_REF] Rhie | Numerical study of the turbulent flow past an airfoil with trailing edge separation[END_REF], are frequently used [START_REF] Nägele | On the influence of different stabilisation methods for the incompressible Navier-Stokes equations[END_REF].

In [START_REF] Chénier | Numerical results using a colocated finitevolume scheme on unstructured grids for incompressible flows[END_REF][START_REF] Touazi | Simulation of natural convection with the collocated clustered finite volume scheme[END_REF], a collocated finite volume scheme for incompressible flows is developed on so called "admissible" unstructured meshes, that are meshes satisfying the two following conditions: the straight line joining the centers of two adjacent control volumes is perpendicular to the common edge, and the neighboring control volumes and the associated centers are arranged in the same order, with respect to the common edge. Rectangular or orthogonal paral-Robert.Eymard@univ-paris-est.fr (Robert Eymard), Raphaele.Herbin@latp.univ-mrs.fr (Raphaèle Herbin). lelepipedic meshes, triangular (2D) or tetrahedral (3D) Delaunay meshes, and Voronoi meshes fulfill these requirements. Under this assumption, the isotropic diffusion fluxes can be consistently approximated by a two-point finite difference scheme. Using this approximation for a pure diffusion problem yields a symmetric "M-matrix" (which ensures monotony); the stencil is limited to the control volume itself and its natural neighbors and it leads to the classical 5-and 7-point schemes on rectangles and orthogonal parallelepipeds. Unfortunately, although the use of such grids considerably widens the variety of geometric shapes which can be gridded, it is far from solving all the critical needs resulting from actual problems:

• For complex 3D domains, it is well known that the use of a large number of flat tetrahedra produces high discretization errors; generalized hexahedric meshes are often preferred: these are made of 3D elevations of quadrangular meshes, for which the faces of the control volumes may no longer be planar.

• To our knowledge, there is yet no mesh software able to grid any geometrical shape in 3D using Voronoi or Delaunay tessellations while respecting the boundaries and the local refinement requirements.

• In compressible flows, the approximation of the full tensor by the usual two point scheme is no longer consistent even on admissible meshes; multi-point approximations are therefore required.

• Boundary layers are classically meshed with refined grids, so that the discretization scheme should be able to deal with non-conforming meshes.

Whereas there is no real difficulty to discretize the convective terms for general non-conforming grids, writing accurate diffusion approximations, particularly relevant for low Reynolds (Péclet) flows, is still a challenge on such meshes.

In the early 80's, Kershaw [START_REF] Kershaw | Differencing of the diffusion equation in lagrangian hydrodynamic codes[END_REF] first proposed a nine-point scheme on structured quadrilateral grids by using the restrictive assumption of a smooth mapping between the logical mesh and the spatial coordinates. Since then, numerous works have been published to efficiently solve the diffusion equations in general geometry (see [START_REF] Breil | A cell-centered diffusion scheme on two-dimensional unstructured meshes[END_REF] for a review of recent papers). The drawbacks of the actual schemes for diffusion are often linked with one or several of these key points:

• a non-local stencil (quite dense matrices);

• cell-centered but also face-centered unknowns (large matrices);

• non-symmetric definite positive matrices (loss of the energy balance);

• loss of the convergence or of the accuracy on some particular grids;

• loss of monotony for solutions in purely diffusive problems (the resulting matrix is not an "M-matrix").

We focus in this paper on the approximation of the Navier-Stokes and energy equations under the Boussinesq assumption, using a new scheme for diffusion terms. This scheme is shown to provide a cell-centered approximation with a quite reduced stencil, leading to symmetric definite positive matrices and allowing a mathematical proof of convergence. Although the diffusion matrix may not be shown to be an M-matrix in the general case, the maximum principle is nevertheless preserved in our numerical three-dimensional simulations.

In this scheme, the discrete pressure gradient and the non-linear contributions are approximated so that the discrete kinetic and energy balances mimic their continuous counterparts. Indeed, the pressure gradient is chosen as the dual operator of the discrete divergence, and the discretization is such that there is no contribution of the non-linear velocity transport in the increase of kinetic energy. In order to suppress the pressure checkerboard modes, the mass balance is stabilized by a pressure term which only redistributes the fluid mass within subsets of control volumes, the characteristic size of which is two or three times the local mesh size.

The remainder of this paper is divided into four sections. In Section 2, the continuous formulation is presented in the framework of free convection. Section 3 presents the discrete scheme for general non-conforming meshes, with illustrations in the simplified case of uniform rectangular grids, and some mathematical properties. The fourth section is devoted to the numerical validation, first with analytical solutions and then with a classical natural convection problem.

Continuous formulation

Let d be the dimension of the space (d = 2 or 3) and let Ω ⊂ R d be an open polygonal connected domain. For x ∈ Ω, our aim is to compute an approximation of the velocity u(x) = d i=1 u (i) (x)e i , the pressure p(x) and the temperature T (x), solution of the steady and dimensionless Navier-Stokes and energy equations under the Boussinesq approximation:

-Pr∆u + ∇p + (u • ∇)u -Ra Pr T e 3 = f (x) in Ω, (1a)
-∆T + (u • ∇)T = g(x) in Ω, ( 1b
) divu = 0 in Ω, ( 1c
)
where e 3 indicates the vertical upward direction, f (x) = d i=1 f (i) (x)e i and g(x) are dimensionless regular functions modeling source or sink in the momentum or heat balances; Pr and Ra denote the Prandtl and Rayleigh numbers respectively. We consider the case of the homogeneous Dirichlet boundary conditions for the velocity and of the mixed Dirichlet-Neumann boundary con-ditions for the temperature. These boundary conditions read as follows:

                           u(x) = 0 x ∈ Γ, T (x) = T b (x) x ∈ Γ 1 , -∇T (x) • n(x) = q b (x) x ∈ Γ 2 , (2)
where Γ 1 , Γ 2 are subsets of the boundary Γ of the domain Ω such that Γ 1 ∩ Γ 2 = ∅ and Γ 1 ∪ Γ 2 = Γ, and n(x) is the outward unit normal vector to the boundary. We assume that T b is the trace on Γ 1 of a function again denoted

T b such that T b ∈ H 1 (Ω)
, and we define the functional space

H 1 Γ 1 ,0 (Ω) = {T ∈ H 1 (Ω); T (x) = 0 on Γ 1 }.
Then a weak formulation of equations (1a-1c) with boundary conditions (2) reads:

Find u ∈ H 1 0 (Ω) d , p ∈ L 2 (Ω) with Ω p(x)dx = 0, and T with T -T b ∈ H 1 Γ 1 ,0 (Ω) such that Pr Ω ∇u : ∇vdx - Ω p divvdx + Ω div(u ⊗ u) • vdx -Ra Pr Ω T e 3 • vdx = Ω f (x) • vdx, ∀v ∈ H 1 0 (Ω) d , (3a)
Ω ∇T • ∇θdx + Ω div(uT )θdx = Ω g(x)θdx - Γ 2 q b (x)θ(x)dx, ∀θ ∈ H 1 Γ 1 ,0 (Ω). (3b) divu(x) = 0 for a.e. x ∈ Ω, (3c)
Although our discretization scheme belongs to the finite volume family, we shall also be using the weak form (3a-3c) in our discretization. Indeed, the discretization of the diffusive terms -Pr∆u in (1a) and -∆T in (1b) is obtained by the construction of a discrete gradient which is then used to approximate the term Pr Ω ∇u : ∇vdx in (3a) and Ω ∇T • ∇θdx in (3b).

Numerical scheme

In this section we present the discretization scheme for Problem (1)-( 2) under its weak form [START_REF] Chénier | Numerical results using a colocated finitevolume scheme on unstructured grids for incompressible flows[END_REF]. The next paragraph is devoted to the notations for general discretization meshes and to the description of the discrete degrees of freedom. We then describe the approximation of the diffusive terms (Sec. 3.2).

Because of the collocated choice of the unknowns, a stabilization is needed.

The stabilization we choose is imposed on the mass flux (rather than the overall balance) and also appears in the momentum and energy equations through the convective contributions: this is described in Section 3.3. It also involves the choice of some coefficients which are defined in Section 3.4. The complete discrete problem is finally given in Section 3.5, and some of its mathematical properties sketched in Section 3.6.

Mesh and discrete spaces

We denote by D = (M, E, P) a space discretization, where (see Fig. 1):

• M is a finite family of "control volumes", i.e. non empty connected open disjoint subsets of Ω such that Ω = ∪ K∈M K. For any control volume K ∈ M, we denote by ∂K = K \ K its boundary, m K > 0 its measure (area if

d = 2, volume if d = 3
) and h K its diameter (that is the largest distance between any two points of K).

• E is a finite family "edges" (d = 2) or "faces" (d = 3) of the mesh; these are assumed to be non empty open disjoint subsets of Ω, which are included in a straight line if d = 2 or in a plane if d = 3, and with non zero measure.

We assume that, for all K ∈ M, there exists a subset E K of E such that Note that there are cases in which K/L includes two or more edges or faces, see for instance the third mesh for the unit cube, section 4.1, and Figure 5b. We also assume that, if σ ∈ E ext , then either σ ⊂ Γ 1 or σ ⊂ Γ 2 . For all σ ∈ E, we denote by x σ and m σ the barycenter and the measure of σ. For all K ∈ M and σ ∈ E K , we denote by n K,σ the unit vector normal to σ outward to K.

∂K = ∪ σ∈E K σ.
• P is a family of collocation points P = (x K ) K∈M of Ω which is chosen such that for all K ∈ M and for all x ∈ K, the property [x K , x] ⊂ K holds.

Note that this choice is possible for quite general polygons, including those with re-entrant corners, see Fig. 1. The Euclidean distance d K,σ between x K and the hyperplane including σ is thus positive. We also denote by C K,σ the cone with vertex x K and basis σ.

Next, for any σ ∈ E int , we choose some real coefficients (β L σ ) L∈M such that the barycenter x σ of σ is expressed by

x σ = L∈M β L σ x L , L∈M β L σ = 1. (4)
In three space dimensions, it is always possible to restrict the number of nonzero coefficients β L σ to four (in practice, the scheme has been shown to be robust with respect to the choice of these four control volumes, taken close enough to the considered edge).

Note that in the case of an uniform rectangular grid which is depicted in Figure 2, an obvious choice for the coefficients β L σ is obtained by noticing that x i+1/2,j = (x i,j + x i+1,j )/2 and x i,j+1/2 = (x i,j + x i,j+1 )/2. Thus for any edge σ, only two coefficients β L σ need to be nonzero.

We now define the finite dimensional space

R M × R E (where an element v ∈ R M × R E is defined by the family of real values ((v K ) K∈M , (v σ ) σ∈E )
) and the following subspaces:

• X D = u ∈ R M × R E , ∀σ ∈ E int , u σ = L∈M β L σ u L (the dimension of X D
is the number of control volumes plus that of boundary edges),

• X D 0 = u ∈ X D , ∀σ ∈ E ext , u σ = 0 (the dimension of X D 0 is the number of control volumes), • X D Γ 1 ,0 = θ ∈ X D , ∀σ ∈ E ext ∩ Γ 1 , θ σ = 0 (the dimension of X D Γ 1 ,0
is the number of control volumes plus that of boundary edges on Γ 2 ).

Discretization of diffusive terms

Let us first define a discrete gradient for the elements of X D on cell K ∈ M.

We set, for any u ∈ X D and K ∈ M:

∇ K u = 1 m K σ∈E K m σ (u σ -u K )n K,σ . (5)
Note that this is a centered gradient.

As an illustration, consider the case of the two dimensional uniform rectangular grid depicted in Fig. 2, let u ∈ X D , and choose the natural choice

β L σ = 1/2 if
σ is a side of L and 0 otherwise. Then with the (natural) notations of Fig. 2, one has:

∇ K i,j u =         1 2h x (u i+1,j -u i-1,j ) 1 2h y (u i,j+1 -u i,j-1 )         ;
if we apply this formula to the element 1 K i,j ∈ X D with all of components equal to zero except for the one associated with K i,j which is equal to 1, we get that the vector ∇ L 1 K i,j is zero for all control volumes L except those neighboring K i,j , as shown on Fig. 3a. Considering the checkerboard solutions on uniform rectangular grids, the above expression shows that this discrete gradient may vanish for some non-constant functions.

An approximation of the diffusive terms Pr Ω ∇u : ∇vdx in (3a) and Ω ∇T • ∇θdx in (3b) using this discrete gradient (5) would yield a non-coercive form.

Thus, we shall work with a modified gradient, defined in ( 6)-( 7) below.

To this end, for all σ ∈ E K we first define R K,σ u ∈ R which may be seen as a consistency error on the normal flux, by:

R K,σ u = √ d d K,σ (u σ -u K -∇ K u • (x σ -x K )) .
One may note that R K,σ u = 0 if u K and u σ are the exact values of a linear function at points x K and x σ , for all K and σ ∈ E K . We then give the following expression for a stabilized discrete gradient of u ∈ X D in each cone C K,σ :

∇ K,σ u = ∇ K u + R K,σ u n K,σ . (6)
In the case of the uniform rectangular grid given in Fig. 2 and for K = K i,j and σ = σ i+1/2,j , we find:

R K i,j ,σ i+1/2,j = √ 2 2h x (u i+1,j + u i-1,j -2u i,j ).
The stabilized discrete gradient (6) applied to the above defined element

1 K i,j
of X D provides nonzero contributions on the triangular subcells C K,σ of the cell K i,j and of its neighbors (Fig. 3b).

The global discrete gradient is chosen as the function ∇ D u:

∇ D u(x) = ∇ K,σ u, for a.e. x ∈ C K,σ , ∀K ∈ M, ∀σ ∈ E K . (7)
We then plan to approximate the term Ω ∇u(x) • ∇v(x)dx by:

Ω ∇ D u(x) • ∇ D v(x)dx = K∈M σ∈E K m σ d K,σ d ∇ K,σ u • ∇ K,σ v, ∀u, v ∈ X D . (8)
In fact, it is shown in [START_REF] Eymard | A new colocated finite volume scheme for the incompressible Navier-Stokes equations on general non matching grids[END_REF][START_REF] Eymard | Discretization schemes for anisotropic 29 diffusion problems on general nonconforming meshes[END_REF] that this expression defines a symmetric inner product on X D , and provides a good approximation for Ω ∇u(x) • ∇v(x)dx; this approximation may be seen as a low degree discontinuous Galerkin method.

If one seeks a finite volume interpretation of this scheme, it is possible, expressing u σ and v σ for all σ ∈ E int thanks to the relations (4), to show that

Ω ∇ D u(x) • ∇ D v(x)dx = K∈M   L∈N K F K,L (u)v K + σ∈E K ∩Eext F K,σ (u)(v K -v σ )   , (9)
where for any K ∈ M, N K is the subset of cells playing a part in the barycenter expression of x σ , for all edges of the cell K and of the neighbors of K, i.e.

N K = {M ∈ M; β M σ ′ = 0, ∀σ ′ ∈ E L , ∀L ∈ M σ , ∀σ ∈ E K }; F K,L (u) is a linear function of the unknowns (u L ) L∈M which is such that F K,L (u) = -F L,K (u).
In the general case, the expression of F K,L (u) is rather complicated (see [START_REF] Eymard | Discretization schemes for anisotropic 29 diffusion problems on general nonconforming meshes[END_REF]).

In the case of the uniform rectangular grid of Fig. 3, this expression simplifies into the usual two point flux; for instance, the flux from K i,j to K i+1,j reads:

F K i,j ,K i+1,j (u) = h y u i,j -u i+1,j h x .
More generally, a two point flux is also obtained in the two or three dimensional non uniform rectangular cases. Indeed, locating x K at the center of gravity of the cell K, the relation

(x σ -x K )/d K,σ = n K,σ holds.
It is then possible to write

x σ = (d L,σ x K + d K,σ x L )/(d L,σ + d K,σ ) and u σ (v) = (d L,σ u K + d K,σ u L )/(d L,σ + d K,σ ) for all σ such that σ ⊂ K/L,
and for all u ∈ X D . Using the identity

σ∈E K m σ (x σ -x K )n t K,σ = m K I
where t designates the transposition and I the identity matrix, we obtain [5, Lemma 2.1]:

Ω ∇ D u(x) • ∇ D v(x)dx = σ∈E int ,σ⊂K/L m σ d K,σ + d L,σ (u L -u K )(v L -v K ) + σ∈Eext∩E K m σ d K,σ (u σ -u K )(v σ -v K ).
Then the previous relation leads to define N K as simply the set of the natural neighbors of K, and to define the fluxes by the natural two-point difference scheme, in the same manner as in [START_REF] Chénier | Numerical results using a colocated finitevolume scheme on unstructured grids for incompressible flows[END_REF][START_REF] Touazi | Simulation of natural convection with the collocated clustered finite volume scheme[END_REF]:

F K,L (u) = m σ d K,σ + d L,σ (u K -u L ) for σ ∈ E int , σ ⊂ K/L F K,σ (u) = m σ d K,σ (u K -u σ ) for σ ∈ E ext ∩ E K . ( 10
)
The classical and cheap 5-and 7-point schemes on rectangular or orthogonal parallelepipedic meshes is then recovered. An advantage can then be taken from this property, by using meshes which consist in orthogonal parallelepipedic control volumes in the main part of the interior of the domain, as illustrated by the cone-shaped cavity (Fig. 5c) in Section 4.1.

Note that the approximation of -K ∆udx is obtained by letting v K = 1,

v L = 0 for L = K and v σ = 0 for σ ∈ E ext in (9): - K ∆udx ≃ L∈N K F K,L (u) + σ∈E K ∩Eext F K,σ (u),
so that we may define an approximate Laplace operator ∆ D by the constant values ∆ K u on the cells K:

-∆ K u = 1 m K   L∈N K F K,L (u) + σ∈E K ∩Eext F K,σ (u)   . (11)
In the case of the uniform rectangular grid of Fig. 3, this discrete Laplacian leads to the usual five point formula:

-∆ K i,j u = 1 h x 2 (2u i,j -u i+1,j -u i-1,j ) + 1 h y 2 (2u i,j -u i,j+1 -u i,j-1 ) .
In the general case, the stencil of the discrete operator on cell K is defined by N K (see relation [START_REF] Touazi | Simulation of natural convection with the collocated clustered finite volume scheme[END_REF]) and therefore depends on the way the barycenters x σ are computed. For general grids, the equation for a given cell usually concerns the unknowns associated to itself, its neighbors, the neighbors of its neighbors and possibly some additional adjacent cells. The resulting matrix is usually not an "M-matrix", except on particular meshes such as conforming orthogonal parallelepipeds, in which case we obtain the usual two point flux scheme, as previously pointed out.

Pressure-velocity coupling, mass balance and convective contributions

For all u ∈ (X D 0 ) d , we define a discrete divergence operator by:

div K u = 1 m K σ∈E K m σ u σ • n K,σ , ∀K ∈ M. ( 12
)
where

u σ = d i=1 u (i) σ e i . Notice that div K u = d i=1 (∇ K u (i) ) (i) ,
with ∇ K u (i) defined by ( 5).

In the case of the uniform rectangular mesh of Fig. 2, this operator reads:

div K i,j u = 1 2h x u (1) i+1,j -u (1)
i-1,j + 1 2h y u

i,j+1u

i,j-1 .

We then define the function div D u by the relation div D u(x) = div K u, for a.e. x ∈ K, ∀K ∈ M.

The discrete gradient operator used for the pressure gradient is defined as the dual of this divergence operator. More precisely, we mimic at the discrete level the (formal) equality

Ω p divv dx = -Ω ∇p • v dx. The discrete equivalent of Ω p divv dx reads L∈M m L p L div L v with v ∈ (X D ) d
and p ∈ X D ; we then define the discrete pressure gradient which we denote

∇ K p (on cell K) , such that L∈M m L ∇ L p • v L = - L∈M m L p L div L v, ∀ v ∈ (X D ) d .
¿From the definition of the divergence (12) and of X D , we thus seek ( ∇ L p) L∈M such that:

L∈M m L ∇ L p • v L = - L∈M p L σ∈E L ∩E int m σ M ∈M β M σ v M • n L,σ . (13)
Taking for v the element of (X D ) d with components v (j) L = 1 if j = i and L = K, and 0 otherwise , we thus get:

m K ∇ K p = - L∈M p L σ∈E L ∩E int m σ β K σ n L,σ = σ∈E int ,σ⊂L/M m σ β K σ (p M -p L )n L,σ . (14)
Remark that ∇ K p is neither constructed with the discrete gradient [START_REF] Eymard | Discretization schemes for anisotropic 29 diffusion problems on general nonconforming meshes[END_REF] nor with the stabilized one (6); its expression is the dual form of the divergence [START_REF] Tric | A first incursion into the 3D structure of natural convection of air in a differentially heated cubic cavity, from accurate numerical solutions[END_REF].

In the case of non uniform rectangles (d = 2) or parallelepipeds (d = 3)

with collocated points at the gravity center of the cells, and for all σ ∈ E int with σ ⊂ K/L, we only need two non-zero coefficients β M σ :

β K σ = d L,σ d K,σ +d L,σ and β L σ = d K,σ d K,σ +d L,σ . Therefore, relation (14) reduces to m K ∇ K p = σ∈E K ,σ⊂K/L m σ d L,σ d K,σ + d L,σ (p L -p K )n K,σ . Using p K σ∈E K m σ n K,σ = 0, m K ∇ K p = σ∈E K ,σ⊂K/L m σ d L,σ p L + d K,σ p K d K,σ + d L,σ n K,σ + σ∈E K ∩Eext m σ p K n K,σ
For a uniform grid and a control volume without boundary faces, the above expression resumes to

m K ∇ K p = σ∈E K ,σ⊂K/L m σ p L + p K 2 n K,σ ,
which provides, in the particular case of Figure 2, the usual formulation

h x h y ∇ K p = h y p i+1,j -p i-1,j 2 + h x p i,j+1 -p i,j-1 2 .
As recalled in the introduction of this paper, a pressure stabilization method is implemented in the mass conservation equation in order to prevent from oscillations of the pressure, as for instance in [START_REF] Brezzi | On the stabilization of finite element approximations of the Stokes equations, Efficient solutions of elliptic systems[END_REF] in the finite element setting, [START_REF] Nägele | On the influence of different stabilisation methods for the incompressible Navier-Stokes equations[END_REF][START_REF] Rhie | Numerical study of the turbulent flow past an airfoil with trailing edge separation[END_REF] in the finite volume setting. The originality of our approach is that we directly include the stabilizing diffusive pressure flux in the approximated mass flux, so that it will appear not only (as usual) in the mass equation, but also in the momentum equation through the non-linear convective term. From the mathematical point of view, this helps in obtaining simple estimates on the velocity and pressure, but more importantly, it ensures that the contribution of the discrete non-linear convective term to the kinetic (and thermal) energy balance is zero, just as in the continuous case. Let us define the stabilized mass flux across σ ⊂ K/L by

Φ λ K,σ (u, p) = m σ (u σ • n K,σ + λ σ (p K -p L )) , (15)
where (λ σ ) σ∈E int is a given family of positive real numbers, the choice of which is discussed below. Note that the quantity λ σ (p Kp L ) may be seen as a numerical pressure diffusion flux, and that the overall numerical flux remains

conservative, that is, if σ ⊂ K/L then Φ λ K,σ (u, p) + Φ λ L,σ (u, p) = 0.
In the case of the uniform rectangular mesh of Fig. 2, the expression of this flux through a vertical edge σ i+1/2,j reads:

Φ λ K i,j ,σ i+1/2,j u = h y 2 u (1) i+1,j + u (1)
i,j + h y λ i+1/2,j (p i,jp i+1,j ).

We then use the modified flux, in order to define a stabilized centered transport operator which is defined, for all u ∈ (X D 0 ) d , w ∈ X D and K ∈ M, by

div λ K (w, u, p) = 1 m K σ∈E K ,σ⊂K/L Φ λ K,σ (u, p) w K + w L 2 .
An interesting remark is that, in the case where the mass balance equation in the control volume K is satisfied, that is:

div λ K (1, u, p) = 1 m K σ∈E K ∩E int Φ λ K,σ (u, p) = 0, then σ∈E K ∩E int Φ λ K,σ (u, p)w K = 0,
so that the following relation also holds:

div λ K (w, u, p) = 1 m K σ∈E K ,σ⊂K/L Φ λ K,σ (u, p) w L -w K 2 .
We shall use this latter form in the practical implementation, in particular in the discretization of the non-linear convection term. Indeed, it is more efficient when computing the Jacobian matrix of the momentum equation, since it avoids summing up values of the same amplitude.

In the particular case of the uniform rectangular grid of Fig. 2, the summation σ∈E K ,σ⊂K/L involves the four edges between the control volume

K = K i,j and its neighbors L = K i+1,j , K i-1,j , K i,j+1 , K i,j-1 . If L = K i-1,j
the expression which appears in the summation reads:

h y   - u (1) i,j + u (1) i-1,j 2 + λ i-1/2,j (p i,j -p i-1,j )   w i,j + w i-1,j 2 .
When the local grid Reynolds (or Péclet) number is much larger than 1, an upwind scheme must be applied that consists in substituting div λ K (w, u, p) by

div λ,up K (w, u, p) = 1 m K σ∈E K ,σ⊂K/L max(Φ λ K,σ (u, p), 0)w K + min(Φ λ K,σ (u, p), 0)w L .
In both the centered or upwind cases, the functions div λ D (w, u, p) and div λ,up D (w, u, p) are defined by their constant values in each control volume. For u, w ∈ (X D 0 ) d , we also define the centered vector divergence operator div λ D (w, u, p) such that the i-th component of div λ D (w, u, p) is equal to div λ D (w (i) , u, p), for i = 1, . . . , d; a similar expression applies for the upwind vector divergence operator div λ,up D (w, u, p).

Choice for the parameters (λ σ ) σ∈E int

Different strategies can be applied to define the parameters (λ σ ) σ∈E int . Amongst all of them we applied the "cluster stabilization method" [START_REF] Chénier | Numerical results using a colocated finitevolume scheme on unstructured grids for incompressible flows[END_REF][START_REF] Touazi | Simulation of natural convection with the collocated clustered finite volume scheme[END_REF]; it consists in constructing a partition of M, denoted G, and setting λ σ = λ > 0 if there exists G ∈ G (such G ⊂ M is called a cluster) with σ ⊂ K/L, K and L belonging to G, and λ σ = 0 otherwise. Here is an example of an algorithm creating a cluster partition:

(1) for all cells K ∈ M, initialize a new cluster if K and its neighboring cells do not already belong to a cluster;

(2) for any remaining isolated cell L, connect it to the closest cluster having the largest number of common edges with L.

This algorithm is now applied to the mesh of figure (4a), the cells being described from left to right, from the lower to the upper row. Figures (4b) and (4c) illustrate the first cluster and the set of the clusters at the end of the first step of the algorithm. Figure 4d shows the clusters after the isolated (unnumbered) cells of figure 4c have been connected. The choice of the value λ depends on the physical problem and it must be chosen both large to avoid the appearance of spurious modes and small to preserve an accurate approximation of the mass equation. In [START_REF] Touazi | Simulation of natural convection with the collocated clustered finite volume scheme[END_REF], comparisons with other stabilization methods were performed and the results showed that solutions are little sensitive to the value of λ. For the natural convection example presented in section 4.2, λ = 10 -8 .

Resulting discrete equations

We denote by T b,D the element T ∈ X D such that T K = 0 for all K ∈ M, T σ = 0 for all σ ∈ E int and all σ ∈ E ext with σ ⊂ Γ 2 , and, for all σ ∈ E ext with

σ ⊂ Γ 1 , T σ = 1 m σ σ T b (x)ds(x). (16)
Let H M (Ω) ⊂ L 2 (Ω) denote the set of functions which are constant in each K ∈ M; for any function q ∈ H M (Ω), we shall denote by q K its constant value on K ∈ M. We then define the mapping P M :

X D → H M (Ω) by v ∈ X D → P M v with P M v(x) = v K for a.e.
x ∈ K and all K ∈ M. We also define the mapping P E :

X D → L 2 (Γ) by v ∈ X D → P E v with P E v(x) = v σ
for a.e. x ∈ σ and all σ ∈ E ext .

Let us then use the previously defined discrete operators to formulate a discrete approximation to problem (3):

Find u = (u (i) ) i=1,d ∈ (X D 0 ) d , p ∈ H M (Ω) with Ω p(x)dx = K∈M m K p K = 0 and T -T b,D ∈ X D Γ 1 ,0 such that: Pr Ω ∇ D u : ∇ D v dx - Ω p div D v dx + Ω div λ D (u, u, p) • P M v dx -Ra Pr Ω P M T e 3 • P M v dx = Ω f • P M v dx, ∀v ∈ (X D 0 ) d , ( 17
) Ω ∇ D T • ∇ D θ dx + Ω div λ D (T, u, p)P M θ dx = Ω g P M θ dx - Γ 2 q b P E θ ds, ∀θ ∈ X D Γ 1 ,0 , (18)
div λ D (1, u, p) = 0 a.e. in Ω. ( 19
)
We then deduce from (17) the d discrete momentum balances over the control volume K, letting v (i) = 1 in K, and 0 otherwise; these equations read, in vector form:

-Pr m K ∆ K u + σ∈E int ,σ⊂M/L m σ β K σ (p M -p L ) n L,σ + σ∈E K ,σ⊂K/L Φ λ K,σ (u, p) u K + u L 2 -Ra Pr m K T K e 3 = K f dx (20)
(where -∆ K u is the vector valued discrete Laplace operator defined by [START_REF] Touazi | Simulation of natural convection with the collocated clustered finite volume scheme[END_REF] for each of its components). Similarly, we deduce from (18) the discrete energy balance over the control volume K, letting θ = 1 in K, and 0 otherwise; this equation reads:

-m K ∆ K T + σ∈E K ,σ⊂K/L Φ λ K,σ (u, p) T K + T L 2 = K gdx. ( 21
)
Recall that, for all K ∈ M, and all σ ∈ E K such that σ ⊂ Γ 1 , the Dirichlet boundary condition ( 16) is given. We deduce from (18) the relation imposed by the Neumann boundary condition for the thermal flux, letting θ σ = 1 and 0 otherwise, for some σ ∈ E K with σ ⊂ Γ 2 :

F K,σ (T ) = σ q b (x)ds(x). ( 22
)
Note that the above relation is natural, accounting for the fact that F K,σ (T )

approximates the heat flux at the edge σ. Finally, we write (19) in a given control volume K:

σ∈E K ∩E int Φ λ K,σ (u, p) = 0. ( 23
)
The stencil of the scheme (20-23) is always determined by that of the diffusion operator.

As previously mentioned, in the case of orthogonal quadrilateral or parallelepiped grids, the diffusion flux through an edge σ gives the classical two-point difference scheme [START_REF] Rhie | Numerical study of the turbulent flow past an airfoil with trailing edge separation[END_REF] and the barycenter coordinates are simply calculated by an arithmetic average,

x σ = (d L,σ x K + d K,σ x L )/(d L,σ + d K,σ ).
Hence, in this case, equations (20)-( 23) then read, for a given K ∈ M:

Q (K) ∆u + Q (K) ∇p + Q (K) u•∇u + Q (K) T = K f dx, E (K) ∆T + E (K) u•∇T = K gdx, M (K) divu = 0, with: Q (K) ∆u = Pr   σ∈E K ,σ⊂K/L m σ u L -u K d K,σ + d L,σ + σ∈E K ∩Eext m σ 0 -u K d K,σ   Q (K) ∇p = σ∈E K ,σ⊂K/L m σ d L,σ d K,σ + d L,σ (p L -p K )n K,σ Q (K) u•∇u = σ∈E K ,σ⊂K/L m σ d L,σ u K + d K,σ u L d K,σ + d L,σ • n K,σ + λ σ (p K -p L ) u K + u L 2 Q (K) T = -Ra Pr m K T K e 3 , E (K) ∆T = σ∈E K ,σ⊂K/L m σ T L -T K d K,σ + d L,σ + σ∈E K ∩Eext m σ T σ -T K d K,σ E (K) u•∇T = σ∈E K ,σ⊂K/L m σ d L,σ u K + d K,σ u L d K,σ + d L,σ • n K,σ + λ σ (p K -p L ) T K + T L 2 , M (K) divu = σ∈E K ,σ⊂K/L m σ n K,σ d L,σ u K + d K,σ u L d K,σ + d L,σ + λ σ (p K -p L ) ,
and with the Dirichlet boundary (16) applied to σ ⊂ Γ 1 and the Neumann boundary condition on σ ⊂ Γ 2 being reduced to m σ (T σ -T K )/d K,σ = σ q b (x)ds(x).

Some mathematical properties

The system of discrete equations (20-23) is a system of non-linear equations.

The mathematical proof of the existence of at least one solution can be shown in the particular case T b = 0 and q b = 0, which we consider in this section.

Indeed, in this case, we can show some a priori bounds on T and u. We first let θ = T in (18). Using the relation

Ω div λ D (T, u, p) P M T dx = 0,
which results from (19), we get

∇ D T 2 L 2 (Ω) d = Ω g P M T dx.
Thanks to a discrete Poincaré inequality which follows from [5, Lemma 5.3],

we get that there exists C T , only depending on the regularity of the mesh and on g, but not on the size of the mesh, such that

∇ D T L 2 (Ω) d ≤ C T .
We then let v = u in (17). We get, thanks to (15) and ( 19),

Pr ∇ D u 2 (L 2 (Ω) d ) d + σ∈E int ,σ⊂K/L m σ λ σ (p L -p K ) 2 = Ω (f + Ra Pr P M T e 3 ) • P M u dx.
Again using the Poincaré inequality, we conclude that there exists C u , only depending on the regularity of the mesh, on Ra, Pr, f and g, but not on the size of the mesh, such that

∇ D u (L 2 (Ω) d ) d ≤ C u .
Hence, using the topological degree method, we can prove the existence of at least one solution. Moreover, these inequalities are then sufficient to get compactness properties, which show that, from a sequence of discrete solutions with the space step tending to zero, we can extract a converging subsequence, for suitable norms. Then we can prove that the limit of this subsequence has a sufficient regularity, in relation with the weak sense provided by [START_REF] Chénier | Numerical results using a colocated finitevolume scheme on unstructured grids for incompressible flows[END_REF]. It is then possible to pass to the limit on (18), ( 17) and (19), using test functions which are interpolation of regular ones. We then get that the limit of the converging subsequence satisfies (3).

Numerical validation

Numerical implementation is performed for three dimensional domains. We first validate our results on known analytical solutions which allow us to compute the scheme's order of convergence. We then turn to a natural convection case which is referenced in the literature.

The domains considered are either the unit cube or a circular centered cone.

In both cases, we use structured (rectangular and non rectangular) meshes, and we denote by N the number of cells in each of the three space directions.

The set of non-linear equations (20-23) is solved by an under-relaxed Newton method where the unknowns are the velocity u K , the pressure p K and the temperature T K and T σ for all K ∈ M and σ ∈ E ext ∩ Γ 2 . The solutions of the linear systems are computed with a parallel Generalized Minimal RESidual method provided by the scalable linear solvers package HYPRE with a preconditioning based on the block Jacobi iLU factorization carried out by the Euclid library [START_REF]The Regents of the University of California[END_REF].

Analytical solutions

We consider two closed cavities, cubic or cone-shaped, in which the fluid flow For any regular function ψ (ψ = (u (i) , i = 1, • • • , 3, p or T ), the scheme's relative accuracy for the usual L ∞ , L 2 and H 1 norms is measured by

ǫ ∞ (ψ) = max K∈M |ψ K -ψ ref (x K )| max K∈M |ψ ref (x K )| , ǫ 2 (ψ) = K∈M m K (ψ K -ψ ref (x K )) 2 K∈M m K (ψ ref (x K )) 2 1/2 , ǫ H 1 (ψ) = K∈M m K |∇ K ψ -∇ψ ref (x K )| 2 K∈M m K |∇ψ ref (x K )| 2 1/2 , (24)
where | • | denotes the usual Euclidean inner product in R 3 . For each of the above relative error, the scheme's order of convergence is defined by the mean slope of the logarithm of the relative error as a function of the logarithm of the largest cell diameter max K∈M h K , the slope being calculated by a least square method.

Three different meshes are studied for the unit cubic enclosure. Except for the last mesh, the points x K are located at the gravity center of the cells.

(1) The first mesh is an uniform mesh consisting of orthogonal parallelepipeds of size 1/N 3 .

(2) The second one (Fig. 5a) is constructed by a smooth mapping between the uniform mesh and the spatial coordinates [START_REF] Breil | A cell-centered diffusion scheme on two-dimensional unstructured meshes[END_REF]. The vertices x s (i, j, k) =

x (l) s (i, j, k) l=1,••• ,3
of the elementary distorted cubes are defined by:

∀(i, j, k) ∈ N([1, n 1 + 1]) × N([1, n 2 + 1]) × N([1, n 3 + 1]),
x (1) s (i, j, k) = 1cos π(i -1) 2N

x (2) s (i, j, k) = (j -1)

N + 0.1 sin 2π(j -1) N sin 2π(k -1) N x (3) s (i, j, k) = (k -1) N + 0.1 sin 2π(j -1) N sin 2π(k -1) N .
(3) The third mesh (Fig. 5b) is constructed from the first one in the follow-ing way: the points x K remains at the gravity centers of the basic mesh whereas the vertices x s of the cells are randomly displaced in each space direction at most of 0.45/N. Unlike the previous mesh, which consists in hexahedra with plane faces, the four edges of a face are now not included in a same plane, with the exception of edges which belong to the boundaries of the cubic domain. Since the consistency of the discrete gradient defined in (5) (and therefore of that defined in ( 7)) holds if the faces σ are plane, we replace each of the non planar faces by two triangular faces.

The second enclosure, the cone-shaped cavity, is bounded by the lateral surface ((x (1) -0.5) 2 + (x (2) -0.5) 2 = ((6 -5x (3) )/12) 2 for x (3) ∈ [0, 1] and by two plane discs (x (1) -0.5) 2 + (x (2) -0.5) 2 ≤ 1/4 for x (3) = 0 and (x (1) -0.5) 2 + (x (2) -0.5) 2 ≤ 1/12 2 for x (3) = 1. To make the most of the scheme which locally reduces into a 7-point difference scheme on parallelepiped cells, the mesh consists of cubes which were cut to match the lateral curved boundary.

Thus, the mesh error tends quadratically to zero with respect to the mesh size. To avoid too large differences of volume sizes between adjacent cells that may deteriorate the numerical accuracy, the boundary cells having a volume less than 0.1/(n i ) d are merged into adjacent cells.

We are first interested in the Poisson problem for a scalar variable (Eq. ( 21)

with u = 0) where Γ 1 = Γ and Γ 2 = ∅. It was first checked that the errors obtained with a linear analytical solution on the different meshes and cavities are of the order of the computer accuracy, even for the coarsest grids. The next analytical test consists in choosing the reference solution 1) , x (2) , x (3) ) = sin(πx (1) ) cos(πx (2) ) cos(πx (3) ) with appropriate Dirichlet boundary conditions (Fig. 6a-c). The orders of convergence are reported in table [START_REF] Breil | A cell-centered diffusion scheme on two-dimensional unstructured meshes[END_REF]. The accuracy of the scheme is close to 2 when considering the L 2 -norm and it slightly decreases with the L ∞ -norm but always remains larger than 1.50. The order of convergence for the gradients (H 1 -norm) is larger than 1.

T ref (x (
Next, we examine the convergence behavior of the isothermal Navier-Stokes equations by setting u ref (x) = ∇ ∧ d i=1 (4x (1) (x (1) -1)) 3 (4x (2) (x (2) -1)) 4 (4x (3) (x (3) -1)) 5 e i and p ref (x) = cos(πx (1) ) cos(πx (2) ) cos(πx (3) ) in (Eq. 1a)

with Pr = 1 and Ra = 0 (note that the dimensionless writing of the equations is meaningless because the current reference velocity is related to the thermal diffusivity which has no reason to appear for isothermal problems.

Another velocity reference should be used, based on the viscous diffusivity).

Table [START_REF] Brezzi | On the stabilization of finite element approximations of the Stokes equations, Efficient solutions of elliptic systems[END_REF] indicates that the convergence rates of the velocity components are larger than 1.90 on the three finer meshes when the relative error is based on the L 2 -norm and first order accurate for the pressure for distorted meshes. In accordance with the diffusion problem when the L ∞ -norm is used, the orders of convergence slightly decrease for the velocity but a convergence rate larger than 1.60 is still measured. The convergence rates of the gradients are better than the expected first order. Unsurprisingly, the L ∞ and H 1 -norms of the pressure do not tend to zero with the mesh size because it simply appears in the momentum equation as Lagrangian multiplier of the mass equation. Thus the only guaranteed convergence for the pressure is based on the L 2 -norm.

Natural convection problem

We consider an air filled unit-cubic enclosure with isolated walls except the two face to face vertical isothermal surfaces at x (1) = 0 and 1. The governing fluid flow equations are solution of system (1) and (2) with f (x) = 0, g(x) = 0, q b (x) = 0, T (0, x (2) , x (3) ) = -0.5 and T (1, x (2) , x (3) ) = 0.5. The Prandtl and Rayleigh numbers are fixed to Pr = 0.71 and Ra = 10 7 and the stabilization parameter is chosen equal to λ = 10 -8 in the mass equation. Because very small boundary layers take place along the walls, the vertices are located at the Gauss-Lobatto points and the collocation points x K at the gravity centers of the cells. To also study the effect of non-cubic meshes, the coordinates of the previous defined vertices are randomly moved, for each direction, of a magnitude at most equal to 0.45 the size of the cell in this direction.

Table [START_REF] Eymard | Discretization schemes for anisotropic 29 diffusion problems on general nonconforming meshes[END_REF] presents the maxima of the velocity components, the average Nusselt number on the isothermal walls and their relative differences with respect to reference data [START_REF] Tric | A first incursion into the 3D structure of natural convection of air in a differentially heated cubic cavity, from accurate numerical solutions[END_REF]. For the cubic meshes, the relative differences with respect to reference data seem to tend to zero: for the finer mesh, our results depart from less than 1% with respect to the reference values. Although the solutions are less accurate for the shaken meshes, a convergent behavior is also observed.

The relatively large differences measured on u (2) ∞ , in comparison with the other velocity components, are probably the result of the flow shape in which the main motion occurs in the (e 1 , e 3 )-planes with a small secondary flow in the transverse planes. It is also interesting to note the good accuracy of the average Nusselt numbers. This accuracy is essentially obtained thanks to the use of the stabilized mass flux Φ λ K,σ (u, p) in the heat transport expression which ensures the conservation of the average heat flux balance on the boundaries of the cavity.

Conclusion

In this paper we presented a new scheme which is well suited for the simulation of incompressible viscous flows on irregular and non-conforming grids.

This possibility seems to open a large field of new applications (grid refinement as a function of an a posteriori error computation, free boundaries, . . . ). We emphasize that the convergence of the scheme may be proven mathematically, and that the obtained numerical results are accurate. Although we presented this scheme in the steady case, its extension to transient regimes is straightforward. In this latter case, one should consider optimizing the linear solving step by using suitable projection algorithms.

List of Figures

1 Cell K and neighbors for a mesh of a two-dimensional (d = 2)

domain Ω: σ an interior edge (σ ⊂ K/L), σ ′ a boundary edge (σ ′ ∈ E K ), m K the measure of the cell K, x K the collocation point, x σ the barycenter of σ of measure m σ , d K,σ the Euclidean distance between x K and σ, n K,σ the unit outward vector normal to σ and C K,σ the cone with vertex x K and basis σ (similar notations apply for the edge σ ′ ). 37 2 Cell K and neighbors for a mesh of a two-dimensional (d = 2)

domain Ω meshed with an uniform rectangular grid with step size h x in the x direction and h y in the y direction: σ = σ i,j-1/2

and σ ′ = σ i+1/2,j two interior edges with σ i,j-1/2 ⊂ K i,j /K i,j-1 x σ

and σ i+1/2,j ⊂ K i,j /K i+1,j , m K = h x h
d K , σ hK σ m σ ′ σ ′ x σ ′ x K m σ L N d K , σ ′ n K,σ n K,σ ′ Γ M Ω K C K,σ
C K,σ ′ Fig. 1.

n K,σ ′ d K,σ ′ = hx 2 n K,σ K = K i,j m σ = h x σ = σ i,j-1/2 M = K i-1,j h y L = K i,j-1 m σ ′ = h y
x σ = x i,j-1/2 σ ′ = σ i+1/2,j d K,σ = hy 2

x K = x i,j

x σ ′ = x i+1/2,j i i -1 j -1 j C K,σ C K,σ ′ h K = h 2 x + h 2 y

  The set E is assumed to be partitioned into external and interior edges (d = 2) or faces (d = 3): E = E int ∪ E ext , with σ ⊂ ∂Ω for any σ ∈ E ext and σ ⊂ Ω \ ∂Ω for any σ ∈ E int . Any boundary edge σ is assumed to belong to a set E K for one and only one K ∈ M; any interior edge σ is assumed to belong to exactly two sets E K and E L with K = L, and in this case σ is included in the common boundary of K and L, denoted K/L.

  and the heat transfer are known a priori. Let p ref , u ref and T ref be some known pressure, divergence free velocity and temperature fields; we then compute f and g by Eqs. (1a,1b) where we have set u ≡ u ref , p ≡ p ref and T ≡ T ref .

  y the measure of any cell K, x K the collocation point (center of gravity), x σ the barycenter of σ of measure m σ , d K,σ = h x /2 the Euclidean distance between x K and σ, n K,σ the unit outward vector normal to σ and C K,σ the cone with vertex x K and basis σ. 38 3 Gradient of the element 1 K i,j of X D for an uniform rectangular grid with step size h x in the x direction and h y in the y direction; (a) centered gradient, (b) stabilized gradient. 39 4 Example of cluster construction; (a) Mesh, (b) First cluster, (c) All the clusters are constructed but isolated cells (unnumbered cells) still remain, (d) The isolated cells are connected to neighboring clusters. 40 5 Example of meshes applied for convergence analysis, N = 20 and x ≡ x (1) , y ≡ x (2) , z ≡ x (3) ; (a) Smooth mapping, (b) Random meshes, (c) Truncated cone mesh. 41 6 Solution accuracy of the Poisson problem as a function of log(max K∈M h K ), for 10 ≤ N ≤ 100 and different geometries and meshes; (a) log(ǫ 2 (T )), (b) log(ǫ ∞ (T )), (c) log(ǫ H 1 (T )

Fig. 2 .

 2 Fig. 2.

Fig. 3 .

 3 Fig. 3.

  Fig. 4.

  Fig. 6.

Table 2

 2 ). 42

			Cubic enclosure		Truncated conic
		cubic meshes smooth meshes random meshes	enclosure
	L 2 -norm	2 (2)	1.96 (1.97)	1.87 (1.93)	2.09 (2.03)
	L ∞ -norm	1.99 (2)	1.81 (1.68)	1.74 (1.88)	1.92 (1.58)
	H 1 -norm	2 (2)	1.50 (1.28)	1.16 (1.07)	1.55 (1.46)
	Table 1

dx 2 dx 3 and their relative differences to reference values [START_REF] Tric | A first incursion into the 3D structure of natural convection of air in a differentially heated cubic cavity, from accurate numerical solutions[END_REF],

∞ , Nu} for the natural convection problem with Ra = 10 7 and Pr = 0.71. 36


Klaus Werner

Iurii Karpenko

Tanguy Pierog

Marcus Bleicher

Kirill Mikhailov

Femtoscopy within a hydrodynamic approach based on flux tube initial conditions

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.

EPOS is a multiple scattering model in the spirit of the Gribov-Regge approach [1]. Here, one does not refer to simply multiple hard scatterings, the elementary processes correspond to complete parton ladders, which means hard scatterings plus initial state radiation. In this case, this elementary process carries an important fraction of the available energy, and therefore we treat very carefully the question of energy sharing in the multiple scattering process. Open and closed ladders have to be considered, in order to have a consistent quantum mechanical treatment. The corresponding graphs are squared, and we employ cutting rule techniques and Markov chains to obtain finally partial cross sections. The cut parton ladders are identified with longitudinal color fields or flux tubes, treated via relativistic string theory.

In case of very high energy pp collisions (at the LHC) or heavy ion scatterings already at RHIC, many flux tubes overlap and produce high energy densities. Let us consider the energy density at an early time in an AuAu scattering at RHIC, as obtained from an EPOS simulation [1]. The energy density shows a very bumpy structure concerning the x-y dependence (transverse coordinates), whereas the variation with the space-time rapidity s is small. There are in particular peaks in the x-y plane, which show up at the same position at different values of η s . So we have sub-flux- 1 The article is published in the original.

tubes which exhibit a long range structure in the longitudinal variable η s , with a typical transverse width of the order of a fermi. This is exactly the width we obtain if we compute the initial energy density in protonproton scattering at the LHC. This means, if a hydrodynamic treatment is justified for Au-Au collisions at RHIC, it is equally justified for pp scattering at the LHC, provided the energy densities are high enough, which seems to be the case.

We are therefore going to employ a new tool for treating very high energy hadronic interactions including a hydrodynamic evolution (even in pp), as explained in [1]. In [1], we test the approach by investigating all soft observables of heavy ion physics, in case of AuAu scattering at 200 GeV. Let us consider proton-proton scattering now. In [2], we compare three different scenarios: the full calculations, including hydro evolution and hadronic cascade (full), the calculation without hadronic cascade (no casc), and the calculation without hydro and without cascade (base). Concerning particle spectra, we refer to [2], where we compare the corresponding calculations with experimental data, for pp scattering at 900 GeV.

The space-time evolution of the "full" hydrodynamic approach will be completely different compared to the "base" approach, where particles are directly produced from breaking strings, as can be seen from Fig. 1, where we plot the distribution of formation points of π + as a function of the radial distance r = Bogolyubov Institute for Theoretical Physics, Kiev c KIT, Institut für Kernphysik, Karlsruhe, Germany d FIAS, JW Goethe Universität, Frankfurt am Main, Germany e Institute for Theoretical and Experimental Physics, Moscow e-mail: werner@subatech.in2p3.fr, karpenko@bitp.kiev.ua, pierog@ik.fzk.de, bleicher@th.physik.uni-frankfurt.de, kmikhail@itep.ru Abstract-In pp scattering at LHC energies, large numbers of elementary scatterings will contribute significantly, and the corresponding high multiplicity events will be of particular interest. Elementary scatterings are parton ladders, identified with color flux tubes. In high multiplicity events, many of these flux tubes are produced in the same space region, creating high energy densities. We argue that there is good reason to employ the successful procedure used for heavy ion collisions: matter is assumed to thermalize quickly, so that the energy from the flux tubes can be taken as initial condition for a hydro-dynamic expansion. This scenario gets spectacular support from very recent results on Bose-Einstein correlations in pp scattering at 900 GeV at the LHC.

(in the pp center-of-mass system (cms)) for the three scenarios. Only particles with space-time rapidities around zero are considered. All calculations in this paper refer to high-multiplicity events in pp scattering at 900 GeV, with a mean dn/dη(0) equal to 12.9. The "base calculation" (dotted line) gives as expected a steeply falling distribution as a function of r. In the two cases involving a hydrodynamic evolution, particle production is significantly delayed, even more in the case of the full calculation, with hadronic cascade. The bump in the two latter scenarios is due to particles being produced from the fluid, the small p t contribution is due to corona particles. This particular space-time behavior of the hydrodynamic expansions should clearly affect Bose-Einstein correlations-what we are going to investigate in the

x 2 y 2 +
following. There is a long history of so-called femtoscopic methods [4], where the study of two-particle correlations provides information about the source function S(P, r'), being the probability of emitting a pair with total momentum P and relative distance r'. Under certain assumptions, the source function is related to the measurable two-particle correlation function CF(P, q) as CF(P, q) = ∫d 3 r'S(P, r')|Ψ(q', r')| 2 , with q being the relative momentum, and where Ψ is the outgoing two-particle wave function, with q' and r' being relative momentum and distance in the pair center-ofmass system. The source function S can be obtained from our simulations, concerning the pair wave function, we follow [5], some details are given in [1].

Here, we investigate π + π -correlations. We calculate the correlation function with Bose-Einstein (BE) quantum statistics included, but no Coulomb corrections. Weak decays are not carried out. We compute correlation functions for different k T interval defined as (in MeV): KT1 = [100, 250], KT3 = [400, 550], KT5 = [700, 1000], where k T of the pair is defined as k T = (|p t (pion 1) + p t (pion 2)|)/2. In Fig. 2a, we show the results for KT1. We compare the three different scenarios: "full calculation" (solid line), "calculation without hadronic cascade" (dashed), and "calculation without hydro and without cascade" (dotted), and data from ALICE [6]. The data are actually not Coulomb corrected, because the effect is estimated to be small compared to the statistical errors. We consider here the high-multiplicity class, with dn/dη(0) = 11.2, close to the value of 12.9 from our simulated highmultiplicity events. We compare with the real data (not polluted with simulations), normalized via mixed events, and we do the same with our simulations. Despite the limited statistics, in particular at large k T , we see very clearly that the "full" scenario, including hydro evolution and hadronic cascade, seems to fit the data much better than the two other ones. Usually people like to extract radii from these distributions, so when we make a fit of the form CF -1 = λexp(-R|q|), in the |q| range from 0.05 to 0.70. We obtain the radii given in the figure. So the radii are very different, varying from 0.69 fm (base approach) to 1.80 fm (full model), which is understandable from Fig. 1. We prefer an exponential fit rather than a Gaussian, simply because the former one works, the latter one does not. We do not want to give a precise meaning to R, it simply characterizes the distribution.

Normalizing by mixed events is something one can easily do experimentally (this is why we compare with these data), but it is clear that one has still unwanted correlations, like those due to energy-momentum conservation, which is not an issue in mixed events. Doing simulations, life is easier. We can take simulations without Bose-Einstein correlations as base line, rather than mixed events. This is referred to as "real/bare" normalization (to be distinguished from the "real/mixed" case discussed earlier). The corresponding results are shown in Fig. 2b, the solid line (full calculation) is now completely horizontal away from the peak region, the radius from the exponential fit is 2.11 fm instead of 1.83 fm for the "mixed" normalization. For the other k T regions, the situation is similar, the final results for all the three k T regions for the full calculation are shown in Fig. 3a, together with the radii from the exponential fit: they are almost identical, around 2 fm.

We get to the same conclusion as outlined in [6]: the radii are k T -independent, contrary to what has been observed in AuAu scattering. How can it be that our hydrodynamic scenario gives a strong k T dependence in AuAu, but not in pp? To answer this question, we compute the "true" correlation function (real/bare normalization) for the calculation without hydro and without cascade (just string decay). The results are shown in Fig. 3b. Surprisingly, here we get a strong k T dependence of the radii, but the "wrong" way: we have 0.64 for KT1 and 1.63 fm for KT5! Actually such behavior is quite normal: the distribution is broader for high-p t particles, because high-p t resonances live longer and can move further out before decaying. This effect is in principle also present in AuAu scattering, but it is much more visible for the small pp system. So in pp we have two competing effects:

• radii increase with k T , due to the bigger size of the source of the high-p t particles compared to the low-p t ones,

• radii decrease with k T , as in AuAu (see [1]), in case of collective flow, due to the p-x correlation.

This p-x correlation exists indeed for the case of hydrodynamic evolutions, and is much smaller in the basic scenario. So in the hydro scenarios, the two competing effects roughly cancel, the radii are k T -independent.

To summarize: we employ a hydrodynamic approach to pp scattering at 900 GeV. A very interesting application is Bose-Einstein correlations. We have shown that, as in heavy-ion scattering, the hydrodynamic expansion leads to momentum-space correlations, which clearly affect the correlation functions. To see the signal is nontrivial due to the fact that in addition to the x-p correlations (which leads to decreasing radii with k T ), there is a second effect which works the other way round: the single particle source size is p t -dependent, which is an important effect in pp, not so in heavy-ion scattering. In this sense we can interpret the k T independence of the radii as a real flow effect. Our simulation does not only reproduce the k T independence, but also the whole correlation functions, which is not at all reproduced from the "base scenario" without hydro and without cascade. So the correlation data provide a very strong evidence for a collective hydrodynamic expansion in pp scattering at the LHC.

The research was carried out within the scope of the EUREA: European Ultra Relativistic Energies Agreement (European Research Group: "Heavy Ions at Ultrarelativistic Energies"). Iu. K. acknowledges partial support by the National Academy of Sciences of Ukraine (Agreement-2012) and by the State Fund for Fundamental Researches of Ukraine (Agreement-2012). T. P. and K. W. acknowledge partial support by a PICS (CNRS) with KIT (Karlsruhe). K. M. acknowledges partial support by the RFBR-CNRS grants no. 08-02-92496-NTsNIL_a and no. 10-02-9311 l-NTsNIL_a.

Femtoscopy within a Hydrodynamic Approach Based on Flux Tube Initial Conditions 1 K.

 1 Werner a , Iu. Karpenko a, b , T. Pierog c , M. Bleicher d , and K. Mikhailov e a SUBATECH, University of Nantes-IN2P3/CNRS-EMN, Nantes, France b

Fig. 1 .

 1 Fig. 1. (Color online) The distribution of formation points of π+ as a function of the radial distance in a high-multiplicity event from pp scattering at 900 GeV, for the following scenarios: the full calculation (solid line), a calculation without hadronic cascade (dashed), and a calculation without hydro and without cascade (dotted).

Fig. 2 .

 2 Fig. 2. (Color online) The correlation functions CF for π + π -pairs as obtained from our simulations, for the three different scenarios, for k T bin KT1, compared to data (points). (a) Normalization via mixed events; (b) normalization via a simulation without BE correlation ("bare").

  π + π + BE correlations real/mixed Base: R = 0.69 No casc: R = 1.39 Full model: R = 1π + π + BE correlations real/mixed Base: R = 0.64 No casc: R = 1.47 Full model: R = 2.11 KT1 EPOS 2.05

Fig. 3 .

 3 Fig. 3. (Color online) Correlation function, normalized by using a simulation without BE correlation ("bare"), for three k T intervals. (a) Full model; (b) calculation without hydro and without cascade.