Aymeric Sonolet
email: aymeric.sonolet@cea.fr

Guillaume Bois
email: guillaume.bois@cea.fr

Mathis Grosso

Adrien Toutant
email: adrien.toutant@univ-perp.fr

Conservative formulations

Keywords: Two-phase flow, Front-tracking, Conservative formulation, Heat transfer, Sharp discontinuity, Finite volume, One-fluid formulation, Interfacial flux reconstruction

and research institutions in France or abroad, or from public or private research centers.

Introduction

Heat transfer simulations in turbulent bubbly flows take an increasing place in numerous domains, such as naval propulsion, vapour generation, simulation of accidents in nuclear power plant. This article concerns the numerical formulation of heat transfer calculations in incompressible bubbly flows without phase-change nor wall for large scale simulation, and especially the treatment of discontinuous properties. This kind of simulation consists in approximating the local variables of the flow at different times on a grid of points. In the case of a large scale simulation of a two-phase flow, the interface is also simulated. The interface calculation, in order to treat discontinuities is of primary importance. Several methods are actively used and none of them seem to take precedence in the literature. Few methods are purely Eulerian, and does not rely on the interface position [START_REF] Ding | Diffuse interface model for incompressible two-phase flows with large density ratios[END_REF], but they present interface diffusion issues.

The most widespread sharp method is the Volume of Fluids (VOF), introduced by Hirt and Nichols [START_REF] Hirt | Volume of fluid (VOF) method for the dynamics of free boundaries[END_REF]. In this method, a colour function (with a value of one in the liquid and zero in the gas) is convected. It holds great simplicity, and it regularizes the discontinuous properties, allowing a stable resolution. But this regularization is artificial, so the interface needs a reconstruction method.

Various improvements of this method have been proposed [START_REF] Massoni | Proposition de méthodes et modèles eulériens pour les problèmes à interfaces entre fluides compressibles en présence de transfert de chaleur: Some models and Eulerian methods for interface problems between compressible fluids with heat transfer[END_REF], [START_REF] Agarwal | Planar Simulation of Bubble Growth in Film Boiling in Near-Critical Water Using a Variant of the VOF Method[END_REF], [START_REF] Palmore | A volume of fluid framework for interfaceresolved simulations of vaporizing liquid-gas flows[END_REF], [START_REF] Malan | A geometric VOF method for interface resolved phase change and conservative thermal energy advection[END_REF]. To improve the interface tracking, Tanguy et al. [START_REF] Tanguy | Développement d'une méthode level set pour le suivi d'interfaces et applications[END_REF] and Tryggvason et al. [START_REF] Tryggvason | Direct numerical simulations of gas-liquid multiphase flows[END_REF] refer to the level set (LS) method. In this paradigm, the interface is represented by a signed 3D distance function. This function is also convected. Different variations of this method exist, such as the coupling of VOF and LS by Sussman et al. [START_REF] Sussman | A sharp interface method for incompressible two-phase flows[END_REF], or the CLS method by Olsson and Kreiss [START_REF] Olsson | A conservative level set method for two phase flow[END_REF], Olsson et al. [START_REF] Olsson | A conservative level set method for two phase flow II[END_REF], or ACLS by Desjardins et al. [START_REF] Desjardins | An accurate conservative level set/ghost fluid method for simulating turbulent atomization[END_REF], Desjardins and Moureau [START_REF] Desjardins | Methods for multiphase flows with high density ratio[END_REF]. Some efforts have also been made to develop this method on unstructured grids by Balcázar et al. [START_REF] Balcázar | A finitevolume/level-set method for simulating two-phase flows on unstructured grids[END_REF].

Benchmarks of the different methods was done by Boniou et al. [START_REF] Boniou | Comparison of interface capturing methods for the simulation of two-phase flow in a unified low-Mach framework[END_REF] and Tanguy et al. [START_REF] Tanguy | Benchmarks and numerical methods for the simulation of boiling flows[END_REF]. Both VOF and LS are used with various interfacial reconstruction possibilities. The most common reconstruction is the PLIC (Piecewise Linear

Interface Construction) method presented by Youngs [START_REF] Youngs | Time-dependent multi-material flow with large fluid distortion[END_REF]. Higher order reconstructions were proposed by Scardovelli and Zaleski [START_REF] Scardovelli | Interface reconstruction with least-square fit and split Eulerian-Lagrangian advection[END_REF]. A benchmark of different interface reconstruction methods was done by Aulisa et al. [START_REF] Aulisa | Interface reconstruction with least-squares fit and split advection in three-dimensional Cartesian geometry[END_REF]. Those methods are not used in this article, but some numerical schemes are inspired from them. Another popular interface calculation method is the front-tracking (FT) method developed by Glimm et al. [START_REF] Glimm | Three-Dimensional Front Tracking[END_REF][START_REF] Glimm | Simple front tracking[END_REF], and further improved by Unverdi and Tryggvason [START_REF] Unverdi | A front-tracking method for viscous, incompressible, multi-fluid flows[END_REF], Tryggvason et al. [START_REF] Tryggvason | A Front-Tracking Method for the Computations of Multiphase Flow[END_REF][START_REF] Tryggvason | Direct numerical simulations of gas-liquid multiphase flows[END_REF]. Its advantages are the precise description of the front for the calculation of surface tension forces and the control over topological changes; it does not need an interface reconstruction method afterwards. Thus, it holds greater precision. Besides, the phase indicator function and physical properties can be considered sharply, without artificial thickening further than the cell size. The drawbacks mainly consist in the difficult implementation of the algorithm and its computational cost. In this method, the interface localization is known with high precision, but the evolution of discontinuous properties is not synchronized with the evolution of volumetric fields (such as velocity or temperature). This leads to a conservation issue or even instabilities in some schemes, as we shall see respectively in Sec. 2.1 and 3.1.2.

In addition to the choice of the numerical method to deal with the interface a discretization has to be chosen. Some use a modified Finite Element Method to tackle the discontinuity and conservation issues such as Garajeu and Medale [START_REF] Garajeu | An XFEM model for incompressible two-fluid flows with arbitrary high contrasts in material properties[END_REF].

This method has a good geometrical flexibility in three-dimensional cases with a good accuracy. In this article the finite volume (FV) method is considered, it is one of the most common methods used. It has the property to conserve the transported value. The equations are integrated over small volumes (in our case, a regular mesh of hexahedral volumes in 3D), divergence terms are changed into the integral of fluxes over the faces of the cell (thanks to Gauss's theorem).

This is fully illustrated in Sec. 1. The other classical method used is the finite difference method [25, 15, . . . ].

The discontinuity of the fluid properties across the interface can be handled with different methods. Usually, the grid discretization does not correspond to the interface localization. When an interfacial condition has to be satisfied across the interface, the ghost-fluid method is commonly used to tackle the discontinuity of velocity for compressible or incompressible flows with phase change. The ghost-fluid method was introduced by Fedkiw et al. [START_REF] Fedkiw | A Non-oscillatory Eulerian Approach to Interfaces in Multimaterial Flows (the Ghost Fluid Method)[END_REF] and applied to the VOF paradigm. It was later improved in [START_REF] Aslam | A partial differential equation approach to multidimensional extrapolation[END_REF], [START_REF] Desjardins | An accurate conservative level set/ghost fluid method for simulating turbulent atomization[END_REF] and [START_REF] Li | Alternating direction ghost-fluid methods for solving the heat equation with interfaces[END_REF]. It was also successfuly adapted to the LS framework by Tanguy et al. [START_REF] Tanguy | A Level Set Method for vaporizing two-phase flows[END_REF] and Gibou et al. [START_REF] Gibou | A level set based sharp interface method for the multiphase incompressible Navier-Stokes equations with phase change[END_REF], and enhanced by Li et al. [START_REF] Li | Alternating direction ghost-fluid methods for solving the heat equation with interfaces[END_REF]; Terashima and Tryggvason [START_REF] Terashima | A front-tracking/ghost-fluid method for fluid interfaces in compressible flows[END_REF] and Donghong et al. [START_REF] Donghong | A ghost fluid based front tracking method for multimedium compressible flows[END_REF] introduced it in the FT framework. Another method to treat interface conditions is the one-fluid integrated formulation described in Prosperetti and Tryggvason [33], Toutant et al. [START_REF] Toutant | Jump conditions for filtered quantities at an under-resolved discontinuous interface. Part 1: Theoretical development[END_REF] and [START_REF] Morel | Mathematical Modeling of Disperse Two-Phase Flows[END_REF] and used in [START_REF] Bois | DNS of Turbulent Bubbly Flows in Plane Channels Using the Front-Tracking Algorithm of TrioCFD[END_REF].

In this method, fluid properties are expressed as a continuous function of the volume fraction of the fluid in the cell. This method is at the foundation of some formulations presented in this article.

All those methods have to tackle the issue of thermal energy conservation, which is particularly difficult to attain in heat transfer because of the jump of properties at the interface. Several studies and proposition were made, but mostly with a VOF method and with phase change [START_REF] Malan | A geometric VOF method for interface resolved phase change and conservative thermal energy advection[END_REF][START_REF] Zhao | Boiling and evaporation model for liquid-gas flows: A sharp and conservative method based on the geometrical VOF approach[END_REF]. The aim of the article is to present the limitations of various formulations in the finite volume front-tracking framework, and to propose improvements. By formulation, we do not mean the specific discretization used to interpolate the values where needed, but rather the way the local thermal equation is written just before its integration on a cell volume (finite volume method) and in time. As represented on Fig. 1, a complete formulation involves multiple coherent choices and approximations. At each step, several choices are possible. This complexity is mainly due to the discontinuity of the fluid properties. The issue of thermal energy conservation can be compared to the conservation of momentum in simulations where there is a discontinuity of density and viscosity which has been a topic of more publications [START_REF] Malan | A geometric VOF method for interface resolved phase change and conservative thermal energy advection[END_REF][START_REF] Desjardins | An accurate conservative level set/ghost fluid method for simulating turbulent atomization[END_REF][START_REF] Vaudor | A Consistent Mass and Momentum Flux Computation Method Using Rudman-Type Technique With a CLSVOF Solver[END_REF] in particular for studies on droplets or liquid jet fragmentation. Some simulations were made with a jump of thermal properties (O(10 -50) in Dabiri and Tryggvason [START_REF] Dabiri | Heat transfer in turbulent bubbly flow in vertical channels[END_REF]), especially in the field of thermocapillary simulations (Nas and Tryggvason [40], Nas et al. [START_REF] Nas | Pattern formation of drops in thermocapillary migration[END_REF]). But these simulations do not specifically concern heat transfer at the interfaces, nor directly observe the conservation of thermal energy. Ma and Bothe [START_REF] Ma | Direct numerical simulation of thermocapillary flow based on the volume of fluid method[END_REF] propose a solution in the VOF scheme to suppress the energy conservation issue, and analyse a simulation using this scheme for a small discontinuity ratio of heat mass capacity and conductivity (1/2). The paper is organized as follows. In the first part, the classical one-fluid formulation is presented for heat transfer resolution (Sec. 1). Then, in Sec. 2 several 3D simulations are done in this paradigm using TrioCFD-IJK [43] to demonstrate first the improvement of the one-fluid formulation in a purely diffusive case, and second the drawbacks of this formulation on the thermal energy conservation when convection is also involved. In Sec. 2.2 this observation is reproduced with a 1D code FluidDyn [START_REF] Sonolet | Flu1dDyn[END_REF], which aims to test and evaluate easily several formulations. Different tests are made in order to characterize the origin of the energy loss: one on the effect of different space interpolation schemes and one on the evolution of the energy loss depending on the space and time steps. Finally, in Sec. 3, several new formulations are presented in the fronttracking, finite volume paradigm; they are tested on the same case as in Sec. 2 and compared with each other.

1. Mathematical formulation of the Temperature One-Fluid method (TOF)

In this section, a complete formulation and the process to obtain it is presented. This formulation is a way to compute the evolution of the temperature field in a one-fluid paradigm with a staggered grid representation. Few improvements are made from the classical one-fluid formulation presented by Tryggvason et al. [START_REF] Tryggvason | Direct numerical simulations of gas-liquid multiphase flows[END_REF], see Eq. [START_REF] Aulisa | Interface reconstruction with least-squares fit and split advection in three-dimensional Cartesian geometry[END_REF]. Their benefits are illustrated in Sec. 1.8.

Front-tracking method (FT)

The FT method is used in direct numerical simulations of bubbly flows with a growing success. It brings more complexity than LS methods, which explains why it is less widespread, but it can reach more accuracy. FT does not lead to numerical diffusion of the interface position and discontinuities can thus be considered sharply, without smearing. This method consists in using a Lagrangian mesh, made of triangular elements, to represent the interface between liquid and gas phases (see Fig. 2). The markers of this mesh are moved by interpolating the velocity from the Eulerian mesh to the Lagrangian mesh. Then, the phase indicator function, computed on the Eulerian mesh, is calculated according to the volume rate of each phase in each cell. To know where the liquid and gas are, an algorithm of connex components is used. Thus, the indicator is 1 in the bubble, 0 outside, and has an intermediate value in cells crossed by the interface. The treatment of those so-called two-phase cells is of great importance as shall be seen in the rest of the paper. Indeed, it is necessary to prevent the diffusion of the physical characteristics between the two phases. The phase indicator function χ k , with the indices k corresponding to the phase in which the indicator equals one, is governed by the following transport equation:

∂χ k ∂t + v k • ∇χ k = 0 (1)

Hypotheses on the flow dynamics

In each phase. For the flow dynamics, the following phenomena are considered in each phase:

• The fluid is incompressible, the divergence of the velocity field is therefore zero.

• The viscous force follows Newton's law of proportionality to the stress tensor.

• Gravity acts in a constant and uniform manner.

• There is no phase-change and the temperature has no effect on the flow dynamics, so the density and the viscosity are constant and uniform per phase. Such flows have the advantage of making the velocity and interfaces independent of the thermal fields. Thus for the same flow dynamics, depending on the material properties of the bubble phase and the liquid phase, different thermal solutions are obtained.

Naviez-Stokes equations in each phase under these assumptions are written:

∇ • v k = 0, (2)
∂ρ k v k ∂t + ∇ • (ρ k v k v k ) = -∇p + ρ k g + ∇ • (µ k ∇v k ) . (3)
At interfaces. The surface tension force is characterized by a constant and uniform surface tension coefficient σ on the surface of the bubbles. This coefficient is supposed independent from the temperature, because of the nature of the flow we aim to model in future works: a boiling flow. Consequently, Marangoni's effects are not considered.

Consequences. Under those hypotheses, the conditions at the interfaces on the velocity field and the mass field are the following.

• Because of the viscosity, there is continuity of the tangential speed.

• As there is no phase change, the normal speed at the interface is continuous.

• As the flow is also incompressible, there is no volume change of the simulated liquid and gas in the domain.

These relations are written as follows,

v l I = v v I (4) k T k • n k I = σκn v I , (5)
with

T k = -p k I + µ k ∇v k + ∇ v k the stress tensor, κ = ∇ s • n v I the local
interface curvature and n v I the normal vector to the interface directed towards the liquid phase. In this section, we focus only on the description of the evolution of the temperature field. The complete set of equations fully describing the evolution of the two-phase mixture is given in Sec. 2.1 by the equation system [START_REF] Fedkiw | A Non-oscillatory Eulerian Approach to Interfaces in Multimaterial Flows (the Ghost Fluid Method)[END_REF]. Under the hypothesis that

• the thermal energy production due to viscous dissipation is neglected,

235

• the conductive heat flux can be written using the Fourier law q k = -λ k ∇T k , one can write in phase k:

∂T k ∂t + ∇ • (T k v k ) = 1 ρ k c p k ∇ • (λ k ∇T k ) (6)
The corresponding interfacial jump conditions associated to this temperature transport equation are given by the temperature continuity

T k | I = T I (7)
and the equality of interfacial fluxes:

k λ k ∇T k | I • n kI = 0 (8)
The interfacial jump of normal velocity is null in this case because phase change is not considered:

k v k • n kI = 0 (9)
One-fluid variables are defined over the whole domain (over both phases) by

T = k χ k T k λ = k χ k λ k ρc p = k χ k ρ k c p k v = k χ k v k (10)
in order to obtain the following equation using Eq. ( 1) and after applying the interfacial jump conditions (7), ( 8) and ( 9):

∂T ∂t + ∇ • (T v) = 1 ρc p ∇ • (λ∇T ) (11)
This method was introduced by Kataoka [START_REF] Kataoka | Local instant formulation of two-phase flow[END_REF] for the equation of movement and the same method was used on the temperature field by Morel [START_REF] Morel | Mathematical Modeling of Disperse Two-Phase Flows[END_REF].

Numerical resolution of the flow dynamics

The one-fluid Navier-Stokes equations are written as follows:

∇ • v = 0 (12) ∂v ∂t + ∇ • (vv) = - 1 ρ ∇p + 1 ρ ∇ • µ ∇v + ∇ v + g - 1 ρ σκn v I δ I (13)
These equations are discretized and solved in finite volumes. The finite volume form of the one-fluid momentum equation Eq. ( 13) is written:

1 V dV ∂v ∂t dV = 1 V - S v v • ndS + 1 ρ a S µ ∇v + ∇ v -pI • ndS + V gdV - 1 ρ a S I σκn I dS , (14)
with V the volume of a dual cell, S its closed surface, S I the interfacial area contained in the dual cell and ρ a the arithmetic averaged density in a dual cell. As explained in du Cluzeau et al. [START_REF] Du Cluzeau | Analysis and modelling of reynolds stresses in turbulent bubbly up-flows from direct numerical simulations[END_REF], the Navier-Stokes equations are then solved by a projection method (Puckett et al. [START_REF] Puckett | A high-order projection method for tracking fluid interfaces in variable density incompressible flows[END_REF]) using respectively fourth-order and second-order central differentiations for evaluation of the convective and diffusive terms on a fixed, staggered Cartesian grid. Fractional time stepping leads to a third-order Runge-Kutta scheme (Williamson [48]). In the two-step prediction-correction algorithm, a surface tension source is added to the main flow source term and to the evaluation of the convection and diffusion operators in order to obtain the predicted velocity (see Mathieu [START_REF] Mathieu | Etudes physique, expérimentale et numérique des mécanismes de base intervenant dans les écoulements diphasiques en micro-fluidique[END_REF] for further information).

The surface tension is computed geometrically from the mesh curvature based on the surface divergence of the interface normal vector. Then, an elliptic pressure equation is solved by an algebraic multigrid method to impose a divergence-free velocity field. TrioCFD has already been used for two-phase (Toutant et al.

[50, 51, 52, 53], Bois et al. [START_REF] Bois | Dns of a turbulent steam/water bubbly flow in a vertical channel[END_REF], Bois and du Cluzeau [START_REF] Bois | DNS of turbulent bubbly flows in plane channels using the front-tracking algorithm of triocfd[END_REF], Bois [START_REF] Bois | Direct numerical simulation of a turbulent bubbly flow in a vertical channel: Towards an improved second-order reynolds stress model[END_REF], du Cluzeau et al. [START_REF] Du Cluzeau | Analysis and modelling of reynolds stresses in turbulent bubbly up-flows from direct numerical simulations[END_REF]) and single-phase (Chandesris and Jamet [START_REF] Chandesris | Boundary conditions at a planar fluid-porous interface for a poiseuille flow[END_REF], Dupuy et al. [START_REF] Dupuy | Turbulence kinetic energy exchanges in flows with highly variable fluid properties[END_REF]) flow studies.

Finite volume, formulation choices

In this part, details of the numerical schemes are discussed, with a particular emphasis on the interpolation of the discontinuous variables such as physical properties. The finite volume methods consist in integrating the diffusion convection equation [START_REF] Olsson | A conservative level set method for two phase flow II[END_REF] over each control volume. The integrated values are later assimilated to the values at the centre of the cell to calculate the values at the faces (convective and diffusive fluxes). The integration of Eq. ( 11) over a control volume V using Gauss's theorem gives

∂T m ∂t = - 1 V S T v • ndS + 1 V V 1 ρc p ∇ • (λ∇T ) dV (15)
with T m = 1/V V T dV the mean temperature on a control volume. But in this expression, the diffusive term is not expressed as the integral of a flux over a surface. In order to do so, the following approximation is used

V k ∇ • (λ∇T ) dV ≈ I k V ∇ • (λ∇T ) dV (16)
which leads to

V 1 ρc p ∇ • (λ∇T ) dV = 1 ρc p h V ∇ • (λ∇T ) (17)
where I k = 1/V V χ k dV is the volume rate of phase k in the cell, and ρc p h is the harmonic mean of the mass heat capacity defined by:

1 ρc p h = I l ρc p l + I v ρc p v (18)
To the best of our knowledge, this expression is original, the usual mean used being a classical arithmetic mean, ρc p a = I l ρc p l + I v ρc p v . This new expression is used with success in simulations without convection (on a comparison with Panda et al. [START_REF] Panda | Influence of Gas Fraction on Wall-to-Liquid Heat Transfer in Dense Bubbly Flows[END_REF], see Sec. 1.8). Finally, the diffusive term is expressed as

1 ρc p h S λ∇T • ndS (19)
which leads to the following finite-volume integration of the temperature equation:

∂T m ∂t = - 1 V S T v • ndS + 1 V 1 ρc p h S λ∇T • ndS (20)

Discretization, finite difference interpolation

Once the simplification of the diffusion is made, one has to express the convective and diffusive fluxes as functions of the mean temperature over the control volume of several cells. In the framework of the staggered grid finite difference method, the scalar values (pressure, temperature, physical properties, phase indicator, etc.) are represented at the cell centres, and their space derivatives are located at the faces of the elementary volumes. Thus, the velocity is at the faces of those elements whereas the temperature is located at their centre.

Hence, in Eq. 21, because of the staggered representation, the interpolation of the different fluxes has to be chosen

∂T m ∂t = - 1 V f T f v f S f + 1 V ρc p h f q f S f (21)
with the indices f indexing the faces of the cell,

v f = v f •n f and q f = (λ∇T ) f •n f . Convection. To calculate the convective term V ∇ • (T v) dV = S T v • ndS at
the faces (S) it is necessary to interpolate the temperature from the centres of the volumes to the faces. The integral is calculated as:

S T v • ndS ≈ f T f v f S f with v f = v f • n f .
The following numerical schemes are used in the rest of the article to calculate T f :

• The upwind scheme is the lowest order stable scheme for convection. It has the advantage of being non-oscillatory and stable, but it leads to a very large numerical diffusion compared to higher order schemes. It is also used by Boniou et al. [START_REF] Boniou | Comparison of interface capturing methods for the simulation of two-phase flow in a unified low-Mach framework[END_REF] near the interfaces.

• The QUICK scheme uses a stencil of size 3. It takes into account 2 upstream values and one downstream. Thus, to interpolate the value at face i + 1/2 if the velocity component is positive in this direction, the values of cells {i -1; i; i + 1} are used. The interpolation is non-linear as it is corrected by slope limiters. It was introduced by Leonard [START_REF] Leonard | A stable and accurate convective modelling procedure based on quadratic upstream interpolation[END_REF].

• The WENO scheme uses a stencil of size 5 and is non-oscillatory, which

gives it an advantage over the QUICK type scheme. In the case where the velocity is positive at face i + 1/2, the temperature is interpolated at this location using the values of the cells {i -2; i -1; i; i + 1; i + 2}.

Three linear interpolations are calculated with the stencils {i -2; i -1; i}, {i -1; i; i + 1} and {i; i + 1; i + 2}. Then, non-linear 280 coefficients are calculated to weight the results of those 3 interpolations and form the final value. It was introduced by Liu et al. [START_REF] Liu | Weighted Essentially Non-oscillatory Schemes[END_REF], [START_REF] Zhang | ENO and WENO Schemes[END_REF].

Diffusion. The diffusive term is written as the product of the facial harmonic mean of the thermal conductivity λ h f and the temperature gradient ∇ f T . The values of λ h and T m are known at cell centres. Face values λ h f and ∇ f T must be deduced from them. The treatment of λ f is handled in the part on discontinuous properties below. For the assessment of the discontinuous temperature gradient at the face ∇ f T , a centred linear discretization of order 2 is used as described in Eq. ( 22). This scheme is classical for a discretization of a gradient to faces.

It has a stencil of size two.

∇ f T = T mi+1 -T mi ∆x (22)
Discontinuous properties. The interpolation of the thermal conductivity λ to faces is done with the following method. First, the value of the conductivity at the cell centre is chosen such that 1/λ h = I l /λ l + I v /λ v where the subscript h denotes the harmonic mean as stated by Patankar [START_REF] Patankar | Numerical Heat Transfer and Fluid Flow[END_REF]. Then, this variable is evaluated at multiple cell centres, which are marked with the subscript i.

Finally, the value of the thermal conductivity λ can be interpolated from the cell centres with a second order harmonic centred interpolation as in Eq. ( 23).

This scheme is particular in the case of interpolation of the conductivity which is piecewise constant and strongly discontinuous. Actually it corresponds to the harmonic mean of λ using a linear interpolation of the void fraction at the face.

λ h f = 2 λ hi λ hi+1 λ hi + λ hi+1 (23)
In this formulation, 1/ρc p h does not require interpolation because it is needed at cell centres.

Temporal integration

In direct simulations of interfacial flows, time integration is usually explicit due to the natural limitations of the sharp properties' variations that strongly constrains the time step. For more accuracy a third-order Runge-Kutta method with low memory borrowing from Williamson [START_REF] Williamson | Low-storage Runge-Kutta schemes[END_REF] is used, but in this article for the simplicity of the notations all temporal discretizations are written using the Euler explicit scheme. Moreover, all results presented use the Euler scheme when it is not specifically expressed otherwise. The time step is limited by the CFL and Fourier numbers.

Integrating Eq. ( 21) over a time step ∆t, one gets:

T (n+1) m -T (n) m = - 1 V ∆t f T f v f S f dt + 1 V ∆t 1 ρc p h f λ h f ∇ f T S f dt (24)
After discretization with an Euler time scheme, Eq. ( 24) becomes:

T (n+1) m -T (n) m = - 1 V f T (n) f v (n) f S f ∆t + 1 V 1 ρc p (n) h f λ (n) h f ∇ f T (n) S f ∆t (25)

Heat diffusion around a single bubble

In this section a 3D simulation of a single bubble was made. The set-up is the same as that presented by Panda et al. [START_REF] Panda | Influence of Gas Fraction on Wall-to-Liquid Heat Transfer in Dense Bubbly Flows[END_REF]. The whole simulation is without velocity as the bubble is static. Initially the temperature of the bubble is different from the temperature of the surrounding fluid. A thermal property jump is set between the bubble and the fluid. As the simulation is symmetric, the results are shown along the radius of spherical coordinates. The simulation was made in order to test the different possibilities of means (arithmetic and harmonic) in the TOF formulation and their effects. A temperature reference solution was computed in 1D with the semi-analytical method described in [START_REF] Panda | Influence of Gas Fraction on Wall-to-Liquid Heat Transfer in Dense Bubbly Flows[END_REF].

The temperature is initialized at T 0 = 373K and T 1 = 273K. The time step for the 1D case (in spherical coordinates) is chosen to be ∆t 1D = 2.00 × 10 -4 s and is coupled with a grid resolution D b /∆r = 80 which has led to convergence according to Panda et al. [START_REF] Panda | Influence of Gas Fraction on Wall-to-Liquid Heat Transfer in Dense Bubbly Flows[END_REF].

In these simulations the heat transfer properties, i.e. the thermal conductivity λ and the mass heat capacity ρc p are constant and uniform in each phase.

The profiles of temperature at different times are reported on Fig. 3. The grey part corresponds to bubble's inside while the white part is for the surrounding fluid. In Figs. 3a and3c, the ratio of thermal conductivity is λ in /λ out = 200

("in" means inside the spherical grey part, i.e. r < D b or r < 1) and the heat capacity is constant. In Figs. 3b and3d, the ratio of heat capacity is ρc p out /(ρc p in ) = 200 and the thermal conductivity is constant. Those ratios are choosen so that the temperature in the phase r > D b stays almost constant and that the interfacial temperature is close to the bubble temperature.

Figure 3 clearly shows that the harmonic mean of 1/ρc p (see Eq. ( 20)) gives good results, contrary to the arithmetic one.

TOF formulation, a conservation issue

At first, the lack of conservation of the integrated TOF formulation is shown on a 3D homogeneous bubbly flow in Sec. 2.1. Then, in Sec. 2.2, a similar one dimensional simulation with temperature convection and diffusion is made in order to demonstrate that this case induces an equivalent energy loss. It justifies the use of one-dimensional simulations in this paper, in particular to study the effect of the convection interpolation scheme (Sec. 2.2.3) and to characterize mesh convergence (Sec. 2.2.4). Those analyses lead to the conclusion that the use of a better scheme with a finer mesh is not sufficient to achieve a reasonable energy conservation in cases with large property discontinuities.

A homogeneous bubbly flow

This simulation was made in the objective to quantify the conservativeness of the TOF method in the bulk of bubbly flows.

The hypotheses stated in Sec. system

∂χ k ∂t + v • ∇χ k = 0 ∇ • (v) = 0 ∂v ∂t + ∇ • (v v) = - 1 ρ ∇p + g + 1 ρ ∇ • µ(∇v + ∇v T ) + 1 ρ σκn δ I ∂T ∂t + ∇ • (T v) = 1 ρc p ∇ • (λ∇T )                        (26)
with p the pressure, g the gravity acceleration, µ the dynamic viscosity, σ the surface tension coefficient and κ the bubble local curvature. It is important to note that the only choice in the equation system ( 26) that has a consequence over the study of the formulation for the temperature equation is that of the incompressibility.

Set up

Geometry. The simulation is periodic in any direction and the computational domain is described in Tab. 3.

Initial conditions. Initially the velocity field was established without any temperature field calculation. Bubbles were uniformly distributed to create a void fraction of α = 6% (see Tab. 2). The temperature is initially set unity in the bubbles and zero in the liquid.

Parameter values. The set-up was chosen so that the temperature in the bubble remains more or less constant, i. e., with a higher diffusivity and heat capacity in the bubbles than in the liquid. Parameters are constant by phase and reported in Tab. 1. The subscript b stands for "bubble". For the flow dynamics, parameters are identical to the case at Re b = 400 in du Cluzeau et al. [START_REF] Du Cluzeau | Analysis and modelling of bubble-induced turbulence from DNS of swarms[END_REF], with Re b = ρ l u r D b /µ l and u r the mean relative velocity between the liquid and the gas ( l denotes a liquid value). In Tab. 2, the Eötvös number (Eo = ∆ρgD 2 b /σ) refers to the deformation capacity of the bubbles, the Prandtl number (Pr l = c p l µ l /λ l ) characterizes the thermal diffusivity with respect to the dynamic diffusivity and α v is the void fraction of the gas ( v denotes a gas value). Finally, numerical parameters are reported in Tab. 3, the timestep is chosen such that the Fourier number F o = λ max ∆t/(ρc p min ∆x 2 ) and the Courant number CF L = u∆t/∆x are lower than 1.

Liquid

Gas

ρ (kg m -3 ) 1.17 × 10 3 8.75 × 10 1 µ (Pa s) 3.53 × 10 -4 1.35 × 10 -5 λ (W m -1 K -1 ) 5.54 × 10 -2 1.54 × 10 1 c p (J kg -1 K -1 ) 6.03 × 10 1 8.07 × 10 3 σ (N m -1 ) 1.81 × 10 -2 g (m s -2 ) 9.81 D b (m) 1.00 × 10 -3

Results and analysis

To exploit this DNS, statistics are taken in space (average according to the 3 spatial directions) for each phase. The time evolution of the temperature of each phase and of the Nusselt number are illustrated in Fig. 4. This simulation shows that the temperatures converge towards each other; this leads to deacreasing heat flux and Nusselt number with time.

N u(t) = A b q • ndSD b λ l ∆T 0 A b (

A 1D periodic bubbly flow with heat transfer

In order to propose and test various formulations in this one-dimensional configuration, a light research code was written in the form of a python library, called FluidDyn, destined to solve one-dimensional flows using the front-tracking and finite volume frameworks. In the following, simulations are done with Flu-idDyn.

The velocity v is considered constant, the bubbles are marked by 1D discontinuous functions computed as two markers for each Heaviside associated to a bubble. The domain is periodic. The only equation that is no longer trivial in this case is the temperature equation [START_REF] Olsson | A conservative level set method for two phase flow[END_REF]. Its convection and diffusion terms keep the same form as in the 3D case without becoming trivial. The major difference with the 3D case is that the velocity is now constant.

Simulation set-up

Physical parameters. To obtain similar conditions as in the previous 3D analysis, the same physical parameters (λ, c p and ρ) as presented in Tab. 1 are used.

The simulation is carried out with the same length in the velocity direction, L tot = 0.02m and the same void fraction α = 6%; to satisfy those conditions in one dimension, we have:

L tot α = N b D b .
The number of bubbles and the bubbles' diameter still needs to be adjusted. The choice was to preserve the interfacial area, and so the number of bubbles in the 1D simulation is set to

N b = 4,

Analysis of temperature profiles

As the simulation is one dimensional, the velocity is constant. It is therefore possible to obtain a solution for the diffusion-convection problem simply by translation of the solution of the diffusion problem. This translation is x = y-vt. very fine grid (twice finer, with ∆x = 2.00 × 10 -5 m). Without convection, the simulation using harmonic means of properties is almost conservative in energy, as observed in Sec. 1.8. We use this methodology to provide reference solutions for convection and convection/diffusion problems. In practice, we translate the solution of the pure diffusion problem by a distance -v∆t for the graphical comparison of the curves. As one can see on Fig. 7, the TOF method leads to an under-evaluation of the bubbles' temperature, which leads to a large energy loss (see Fig. 6). A loss of energy in this simulation is observed as in the three-dimensional case presented in Sec. 2. It encourages us in considering the test case of Sec. 2.2 as a simple but meaningful test case for three-dimensional simulations. Thus, this test case is used to evaluate new formulations using FluidDyn solver.

In all the following, 1D simulations are performed with the same parameters as presented in Sec. 2.2, except for the number of bubbles which is chosen to be one in order to increase the resolution without increasing the computational cost. Thus, α, L tot and n tot are kept the same, but a I and D b change, leading to a decreasing energy loss over the same time (the exchange surface diminishes)

and an increased resolution (the number of cells per bubble diameter increases from around 8 to 30).

Convective schemes

This study is carried out to validate the impact of using a QUICK scheme rather than a WENO scheme. The set-up is identical to Sec. The superposition of the results on the finer meshes for the three interpolation methods can be seen on Fig. 9. As expected, the WENO scheme gives the best results, the upwind scheme is far too diffusive compared to the other two, and the QUICK is not too diffusive. Its potential oscillatory behaviour is compensated by the diffusion. As the energy loss is almost the same using WENO or QUICK (see Fig. 8), it shows that the energy is not lost only due to the lesser capability of the QUICK scheme. Moreover the temperature profileis are almost identical with both schemes (see Fig. 9); hence the QUICK scheme is satisfactory in those conditions and with the mesh refinement considered. This leads to the conclusion that the enhancements that have to be made to attain greater energy conservation are mainly not a matter of interpolation.

Mesh convergence for a convection-diffusion simulation

The mesh convergence observed in Fig. 10 has a coefficient of approximately -0.38. This gives us a rough idea of the number of cells per diameter the mesh will need in order for the energy loss to be lower than 5% using a logarithmic regression. As an example, on this short simulation time, the mesh requirement would be D b /∆x > 700, which is approximately 30 times the number of cells per diameter originally used in the 3D simulation; in three dimensions, the number of cells should then be multiplied by 27000. Moreover, considering that in this simulation the interfacial area is four times less than in the 3D simulation presented, reaching 5% of energy loss would be even harder in the 3D homogeneous bubbly flow. Finally, it should be noted that the mesh chosen for the 3D simulation gives a result almost converged for the corresponding dynamic simulation (see [START_REF] Du Cluzeau | Analysis and modelling of bubble-induced turbulence from DNS of swarms[END_REF]). Consequently, a mesh convergence on the full three-dimensional simulation is not a satisfactory solution to tackle the energy conservation issue.

In summary, the new formulation based on the one-fluid temperature proposed in Sec. 1.7 relies on harmonic means of heat capacity and conductivity.

However, even though this change gives good results in a pure diffusion case, it fails to eliminate the large energy loss observed in a convection-diffusion case.

Moreover, only few improvements can be expected from the change of the interpolation method as observed in Sec. 2.2.3, and mesh convergence of the energy loss is really slow compared to an acceptable mesh refinement. The conservation issue is much more complex and has to be treated using a different formulation.

In Sec. 3, new temperature and energy formulations are then introduced, different from that presented in Eq. ( 25) and used until now.

New formulations: quest for conservative and accurate formulations

In this section, new and original methods are presented in order to hold greater conservativeness. First, come two formulations starting with the temper-

T m n, v f V l , ρC p l , λ l , I l S I V v , ρC p v , λ v I v , S f v S f l

Improvement of the temperature equation integration in the temperature formulation

The local one-fluid temperature equation used in this part is Eq. ( 11). In Sec. 1, the mean values over the cell of the discontinuous properties λ and ρc p are used. Sec. 2.1 showed that the main drawback of this one-fluid approach (TOF) is that the evolution of the fluxes is inaccurate and the total energy is not naturally conserved in the presence of convection. This section presents alternative methods in order to improve the integration of the temperature Interp.:

T f = f (T m ), q f = λ f ∇ f T , with λ f = f (λ h ) and ∇ f T = f (T m ).

Hyp. (17): diffusion divided by

ρc p h .
EOFm, Eq. ( 43).

Interp.:

h f = f (h m ), q f = λ f ∇ f T , with λ f = f (λ h ) and ∇ f T = f (T m ).
Hyp. [START_REF] Vaudor | A Consistent Mass and Momentum Flux Computation Method Using Rudman-Type Technique With a CLSVOF Solver[END_REF]:

h m = ρc p a × T m .
Single-phase TSP, Eq. ( 31).

Interp.:

T f = f (T k , T I ), q f = f (λ k , T k , q I )
with

T I = f (λ k , T k ) and q I = f (λ k , T k ).

Hyp. (17): diffusion divided by

ρc p h .
ESP, Eq. ( 49).

Interp.:

T f = f (T k , T I ), q f = f (λ f , T k , q I )
with

T I = f (λ k , T k ) and q I = f (λ k , T k ).
Hyp. ( 45 Hyp. [START_REF] Toutant | Jump conditions for filtered quantities at an under-resolved discontinuous interface. Part 1: Theoretical development[END_REF], q constant over a whole face/interface during a time step.

equation [START_REF] Olsson | A conservative level set method for two phase flow II[END_REF]. There are two possibilities to improve the conservativeness of the temperature formulation:

• either keeping the formulation presented in Sec. 1, but calculating the fluxes with an interpolation of single-phase and interfacial values at the faces as shown in Fig. 12,

490

• or by separating the volume integral over both phases and using the discontinuous properties of each phase and interfacial fluxes.

Both methods need the calculation of interfacial values to make single-phase interpolations (interpolation without using values at mixed cells), which is presented immediately in Sec. 3.1.1. The temperature reconstruction is done using the continuity of the temperature (Eq [START_REF] Tanguy | Développement d'une méthode level set pour le suivi d'interfaces et applications[END_REF]) and of the normal heat flux (Eq (8)) at the interface along with the Taylor series of the temperature in each phase k:

T k (x) = T I + n i=1 ∂ i T ∂x i I k (x -x I ) i i! + O((x -x I ) n ) (28)
∂ i T /∂x i | I k is the i th derivative of the temperature evaluated at the interface in the phase k (it is discontinuous across the interface for i = 1). By evaluating the expression in adjacent single-phase cells, we obtain the following first-order

reconstruction            T i-1 = T I -∇T I -d i-1 T i+1 = T I + ∇T I + d i+1 λ I -∇T I -= λ I + ∇T I + (29)
with d i the distance between the centre of the cell i and the interface, and ∇T I ± the space derivative of the temperature immediately before of after the interface.

The equation system [START_REF] Tanguy | A Level Set Method for vaporizing two-phase flows[END_REF] can never be stiff because the distances d i-1 and d i+1 are greater than ∆x/2. One can note that the temperature T i in the mixed element is not used at all in this reconstruction. The resolution of this system is used to deduce T I , ∇T I -and ∇T I + .

Temperature with Single-Phase fluxes (TSP)

The main idea in this section is to compute single-phase fluxes i. e., fluxes interpolated from phase or interfacial values only, but not from values of the other phase nor mean values on mixed cells. Then, the formulation used is still that presented in Eq. ( 20), but the fluxes interpolations are changed in the following ways:

Temperature interpolation at the faces T f . In 1D, a face is always pure thus it can be associated to a given phase. In this method, T f is an interpolation of T m of single-phase cells (with the same phase as the face) and of the interfacial temperature T I . The interpolation used for the temperature is a third order upwind interpolation, using q I = λ∇T | I • n I if needed (on the face i + 1/2). The values used are represented on Fig. 13. Here are the interpolations used near the interface assuming the velocity is positive from left to right:

• for i-3/2 and i-1/2, a linear third order upwind scheme is used, with a 3 points stencil, using the temperature and position of the interface instead of the temperature of the two-phase cell i if needed.

• on i+1/2, a linear third order upwind scheme is used, using the interfacial temperature instead of the temperature of the mixed cell i, the mean temperature of the cell i + 1, and the interface flux q I .

• for i+3/2 and i+5/2, a linear third order upwind scheme is used, with a 3 points stencil similar to the first case, using the temperature and position of the interface instead of the temperature of the two-phase cell i if needed.

525

In those interpolations, the spacing between points is non-uniform due to the moving position of the interface. Besides, at face i + 1/2, the interpolation is based on two temperatures and a gradient instead of 3 temperatures. Heat flux q f discretization at the faces. The discretization used is a second order centred linear interpolation when the cells are single-phase. On the particular 530 faces i -1/2 and i + 1/2, T I and q I are used with the pure values of the temperatures at i ± 1 and i ± 2 in a third order "centred" interpolation (the same number of values is used from each part of the face, but they are not necessarily evenly spaced).

T i-3 T i-2 T i-3 2 T i-1 T i-1 2 T i T i+ 1 2 T i+1 T i+ 3 2 T i+2 T i+3 T I T i-3 T i-2 T i-3 2 T i-1 T i-1 2 T i T i+ 1 2 T i+1 T i+ 3 2 T i+2 T i+3 T I T i-3 T i-2 T i-3 2 T i-1 T i-1 2 T i T i+ 1 2 T i+1 T i+ 3 2 T i+2 T i+3 T I q I T i-3 T i-2 T i-3 2 T i-1 T i-1 2 T i T i+ 1 2 T i+1 T i+ 3 2 T i+2 T i+3 T I
Heat mass capacity ρc p h time integration. The integral over time of 1/ρc p h is calculated using a Cranck-Nicolson second order centred time scheme:

∆t 1 ρc p h dt ≈ ∆t 2 1 ρc p (n+1) h + 1 ρc p (n) h (30)
This change brings more precision when the jump 1/ρc p l -1/ρc p v is large, but it is not mandatory. Simulations were made without this improved time scheme in FluidDyn and the results were not very different for the chosen parameters.

Using all those discretizations, the following equation is obtained: This section presents a method based on the decomposition of the volume integral over each phase separately. A similar method was also introduced very recently by Zhao et al. [START_REF] Zhao | Boiling and evaporation model for liquid-gas flows: A sharp and conservative method based on the geometrical VOF approach[END_REF] in the context of simulations of boiling flows. The concept is to write conservative equations over both sub-volumes of each phase.

T (n+1) m -T (n) m = ∆t V   - f T (n) f v (n) f S f + 1 2 1 ρc p (n+1) h + 1 ρc p (n) h f q (n) f S f 
Splitting the volume integral, Eq. ( 20) becomes

∂T m ∂t = 1 V - S T f v f • ndS + 1 ρc p l S l λ l ∇T l • ndS + 1 ρc p v Sv λ v ∇T v • ndS (32) = 1 V - S T f v f • ndS + S λ ρc p ∇T • ndS + 1 ρc p l - 1 ρc p v S I λ∇T I • n I l dS (33)
using the interfacial flux equality Eq. ( 8). As 1/ρc p is not continuous over a face during a time step when the interface pass through the face, but λ∇T is, the following approximations are made:

v f ∆t n f , v f , q f T f S I n I l S f v ρc p v S f l ρc p l T m
∆t S f λ ρc p ∇T • ndS dt ≈ q (n) f k 1 ρc p k ∆t S f k (t)dt ∆t S I λ∇T • n I l dS dt ≈ q (n) I l ∆t S I dt (34)
The different terms, and especially the part of the surface S f k wetted by phase k, are defined on Fig. 14. Here, the face values of temperature T f and heat transfer q f are computed in the same way as in Sec. 3.1.2. q I l is the normal interfacial heat flux going from the liquid to the gas. After time integration, discretization and using Eq. ( 34) it gives:

T (n+1) m -T (n) m = - ∆t V f T (n) f v (n) f S f + 1 V f k 1 ρc p k ∆t S f k (t)dt q (n) f + 1 V 1 ρc p l - 1 ρc p v q (n) I l ∆t S I dt (35)
3.1.4. Comparison of the temperature formulations (TOF, TSP, TSV)

As one can see on Fig. 15, none of those new formulations is purely conservative, but both new formulations hold better conservation properties. The default of the TSP method is that it gains energy and the temperature is no longer symmetric as it should be (see Fig. 16). TSP only gains half as much energy as the TOF method loses; moreover, the TSV method gains energy too, but less than the TSP formulation. The advantages of those methods lay in the fact that they still use a temperature point of view. But obviously, their conservation is not sufficient given the highly difficult conditions (large ρc p and λ jumps) used in this study. In both cases as one can see in Fig. 16, the diffusion profile is not symmetrical with respect to the bubble centre at x = 0.5D b as it should be. The temperature gradient inside the bubble is almost constant. The peaks observed at the interface of the bubble comes from a flux jump across the interface that should not be observed, thus showing that the diffusion flux is miscalculated at the interface. Hence, the new methods are more conservative but they produce temperature profiles less satisfactory than the TOF method.

Energy formulation

Even if the improvements proposed in Sec. local equations, the jump conditions (Eqs. ( 8), ( 7) and ( 9)) and the one-fluid conventions (Eq. ( 10)), is used:

∂ρc p T ∂t = -∇ • (ρc p T v) + ∇ • (λ∇T ) (36)
3.2.1. Energy One-Fluid formulation (EOF and EOFm)

After space integration, equation ( 36) becomes:

∂ V ρc p T dV ∂t = - S ρc p T v • ndS + S λ∇T • ndS (37)
To discretize it, two new problems arise: one for the computation of ∇T and the other for the computation of ρc p T at the faces. On the one hand, assessing the gradient ∇T at the faces requires an expression to build the temperature at the cell centres; this issue is addressed using Eq. [START_REF] Vaudor | A Consistent Mass and Momentum Flux Computation Method Using Rudman-Type Technique With a CLSVOF Solver[END_REF]. On the other hand, the interpolation of ρc p T to the face is not straightforward because of the discontinuity of ρc p which in practice makes a linear interpolation very inaccurate and impractical. For these two reasons (flux assessment and non-linear interpolation of the discontinuous ρc p T ), we want to keep the cell temperature T m and the cell mass heat capacity ρc p a as separate variables instead of using a cell enthalpy defined as h m = V ρc p T dV . Then, the following approximation is used and it is exact when the cell is single-phase:

h m ≈ 1 V V T dV × 1 V V ρc p dV = ρc p a × T m (38)
With this hypothesis, a straightforward method can be written to discretize the energy equation [START_REF] Zhao | Boiling and evaporation model for liquid-gas flows: A sharp and conservative method based on the geometrical VOF approach[END_REF]:

ρc p (n+1) T (n+1) m -ρc p (n) T (n) m = - ∆t V f (ρc p (n) a T (n) m ) f v (n) f S f + ∆t V f λ (n) h f ∇ f T (n) S f (39)
However, as Eq. ( 38) is no longer valid for cells crossed by the interface, this formulation leads to numerical instabilities. Precautions have to be taken to ensure the convergence of the scheme using the cell temperature T m . A way to take care of the instabilities that occur in the straightforward EOF method is to use a modified EOF method (EOFm) inspired by a method presented by Boniou et al. [START_REF] Boniou | Comparison of interface capturing methods for the simulation of two-phase flow in a unified low-Mach framework[END_REF] for the density and momentum. A convection equation of ρc p is written as follows

ρc p * a -ρc p (n) a = - 1 V ∆t S ρc p v • ndSdt (40)
where ρc p * a is a temporary variable for the time advancement of the heat mass capacity. For simplicity ρc p * a is written as an explicit Euler evolution of ρc p , but it is possible to do this evolution using a Runge-Kutta scheme if the rest of the numerical scheme is also Runge-Kutta. This ρc p * a is obtained consistently with T m during a convection step:

ρc p * a T (n+1) m -ρc p (n) a T (n) m = 1 V ∆t - S ρc p T v • ndS + S λ∇T • ndS dt (41)
Manipulating equation [START_REF] Nas | Pattern formation of drops in thermocapillary migration[END_REF] to explicitly introduce the time increment of temperature and using equation ( 40) leads to:

ρc p * a T (n+1) m -T (n) m = 1 V ∆t -T (n) m S ρc p v • ndS - S ρc p T v • ndS + S λ∇T • ndS dt (42)
The convective flux of ρc p T is interpolated at the faces from the value at the cell centre of ρc p a T m . The convective flux of ρc p at the faces in Eq. ( 40) is calculated using the same interpolation scheme to ensure accuracy. It is an essential precaution to use the same scheme for both interpolations to guarantee the coherency of both fluxes. Besides, the diffusive flux is defined as the product of λ h f and ∇ f T . It is easy now to calculate ∇ f T as a function of T m as the cell temperature is chosen as a main variable of the problem even in this energy formulation. Finally, after time integration and an Euler explicit discretization, Eq. ( 43) is obtained:

ρc p * a T (n+1) m -T (n) m ∆t = - 1 V T (n) m f ρc p (n) a f v (n) f S f - 1 V f (ρc p (n) a T (n) m ) f v (n) f S f + 1 V f λ (n) h f ∇ f T (n) S f (43)
Face interpolations are made using the following schemes:

555

• the face energy (ρc p a T m ) f and mass capacity ρc p a f are interpolated with a QUICK scheme from ρc p a T m and ρc p a , and as stated in [START_REF] Boniou | Comparison of interface capturing methods for the simulation of two-phase flow in a unified low-Mach framework[END_REF], an upwind interpolation is used when the interface crosses the QUICK stencil;

• the thermal conductivity λ h f is interpolated with a harmonic second-order scheme from λ h as in Eq. ( 23);

• the temperature gradient ∇ f T is computed with a second-order central difference scheme from T m as in Eq. ( 22).

It is responsible for the loss of the energy conservation property of the numerical scheme, which was ensured by the formulation [START_REF] Zhao | Boiling and evaporation model for liquid-gas flows: A sharp and conservative method based on the geometrical VOF approach[END_REF]. But in this formulation, it is inevitable due to the fact that keeping ρc p * diffuses this property, and not going through this intermediate leads to a numerical instability. The hypothesis which seems the most flawed is that of the continuity of ρc p T v. Therefore, it is necessary to compute the fluxes (especially the convective one) to the faces more accurately.

Energy with Single-Phase fluxes (ESP)

The first and highest order error lies in the values of the physical characteristics at the interfaces. In 1D, it seems obvious that the physical properties at a face are single-phase and not weighted by a phase indicator function. Actually, a problem arises when the interface passes through a face during the time step.

Indeed, this problem appears clearly if one considers only convection. If we use a pure ρc p while it changes during the time step, an important error on the energy transfer through this face occurs. This error comes from the inequality (for an explicit Euler integration):

∆t S ρc p T v • ndSdt = ∆t f ρc p (n) f T (n) f v (n) f S f (44)
As the temperature is continuous, we can imagine a case where the temperature is constant everywhere and the interface was immediately upstream of the face at step n. In this case, we have a maximum error on this time step which is

(ρ l c p l -ρ v c p v )T vS∆t.
When integrated in time and space, this error is large even though it is located only at the interface. It is thus necessary to correct this error by better estimating the integral ∆t S ρc p T v • ndSdt. To improve this estimation, the effect of the temperature's variation on a face and on a single time step is neglected compared to the variation of the thermal mass capacity.

Then, the integral on a face S f is calculated in the following way

∆t S f ρc p T v • ndSdt ≈ ∆t S f ρc p dSdt × T (n) f v (n) f (45) = ∆t k ρc p k S f k (t) dt × T (n) f v (n) f (46)
with v f = v • n. In the particular case where the simulation is one-dimensional, this integration can be easily calculated as:

∆t ( k ρc p k S f k (t))dt = S f k ρc p k ∆t χ k (t)dt (47)
because χ k is either zero or one on the whole face, and its time evolution is known. Thus, the time integration of the wetted surface S f k (t) is straightforward 575 in this case.

The same time/phase integration is not made on λ because the product λ∇T is continuous. Then, λ∇T is directly estimated from λ∇T I and the neighbouring single-phase values. Interpolations of interfacial values are described in Sec. 3.1.1. Now that the discontinuous terms have been integrated by phase, the use of ρc p * is no longer needed to maintain stability. And so this time the formulation is genuinely conservative. It can directly be written as

ρc p (n+1) a T (n+1) m -ρc p (n) a T (n) m = 1 V ∆t - S ρc p T v • ndS + S λ∇T • ndS dt (48)
and after discretization it gives:

ρc p (n+1) a T (n+1) m -ρc p (n) a T (n) m = ∆t V   - 1 ∆t f k ρc p k ∆t S f k (t)dt T (n) f v (n) f + f q (n) f S f   (49)
The same reasoning as in the temperature formulation can be done in the separation of volumes, but it actually leads to the same equation (Eq. ( 49)).

That is because no hypothesis on mean properties is made (such as in Eq. ( 17)), and as the energy is a conservative variable, its energy balance over a volume equals the sum of both balance over two sub-volumes. So the separation on each phase does not lead to a different formulation.

Comparison of the energy formulations (EOFm and ESP) to the standard temperature formulation (TOF)

In the energy paradigm, the main advantage lays in the inherent conservation properties of the formulation. The problem comes from the difficulty to recover the temperature gradient ∇ f T necessary to the diffusive term. Because of this difficulty, the formulation adapted for stability reasons (EOFm) uses a convected volume heat capacity, which thus is numerically diffused, and which is forgotten at the end of the time step. As a consequence, this methods is no longer conservative. As one can see in Fig. 17, it leads to an augmentation of the global thermal energy in the simulation. This comes from the fact that the gas temperature is over-estimated and that this phase contains most of the energy (see Fig. 18), as the higher value of heat capacity is in the bubbles in our conditions. Moreover, the temperature gradient inside the bubble is clearly different from the reference one as it is almost always positive, which is non physical (it should be negative in the part between the centre of the bubble and

x/D b = 1, so that it is symmetrical). The change from ρc p * to ρc p (n+1) creates an energy increase in the mixed cell. On the contrary, the ESP method is purely conservative. This perfect conservation property comes from the finite volume method used with fluxes taken identicals on a face shared by two cells. This does not mean that there are no errors in this scheme, the errors come from the fluxes estimations as classically encountered with any discretization method.

However, one can see on Fig. 18 that the reconstructed temperature profile is very satisfactory and almost matches the reference temperature. The diffusive flux at the interfaces seems to be very slightly under-estimated compared to the reference one because the temperature is slightly greater than the reference in the bubble and slightly lower than the reference in the liquid far away from the bubble.

Convergence order of the ESP mehtod

The order of convergence of the newly developed method (the ESP method) is compared to the convergence of the temperature one-fluid (TOF) method.

The simulations use the same parameters than in Sec. 2.2.1, reported Tab. 1, but with only one bubble. The reference used to compute the L 2 error is the result of a computation with the TOF method, without convection. This reference is computed with D b /∆x = 170. The error is evaluated after a simulation time of 0.10 s. This corresponds to the travel of a domain length by the bubble. As the meshes used might not coincide with the reference mesh, a linear interpolation of the reference temperature at the tested mesh is done. For the meshes assessed, that the proposed numerical formulation possesses both a higher order and a smaller initial error for the range of resolution considered. The poor order of the TOF method comes form the lack of conservation. The second order observed for the ESP method corresponds to the diffusion method order.

Conclusion

In this article, the energy conservation property of different numerical formulations was evaluated for the temperature equation without phase change in the finite volume/front-tracking framework. At first, the natural method considering the temperature as a one-fluid variable was improved (see Sec. On the contrary, an energy formulation is naturally conservative, provided that the temperature is calculated with a consistent ρc p over time, but it is not the case for the stability. The temperature has to be reconstructed with particular care so as to have a stable numerical scheme and a correct assessment of

6 FV Finite Volume, page 6 LS Level Set, page 5 PLIC 20 χ 2 V 3 Pr 21 n 12 q 2 v 21 D 21 g- 3 I 21 T 17 amean, page 14 bthe bubbles, page 21 f 15 hmean, page 14 Ithe interface, page 12 i 15 k 10 lthe liquid, page 14 m 14 v

 665202321122212132117142115141215101414 Level Set, page 5 CFL Courant-Friedrichs-Lewy, page 17 CLS Conservative Level Set, page 5 DNS Direct Numerical Simulation, page 22 EOFm Energy One-Fluid modified, page 31 ESP Energy Single-Phase, page 31 FT Front-Tracking, page Piecewise Linear Interface Construction, page 6 QUICK Quadratic Upstream Interpolation for Convective Kinematics, page 15 TOF Temperature One-Fluid, page 18 TSP Temperature Single-Phase, page 31 TSV Temperature Separated Volumes, page 31 VOF Volume Of Fluids, page 5 WENO Weighted Essentially Non-Oscillatory, page 15 Greek symbols α Void fraction, or gas volume fraction over the whole domain, page Indicator function, page 10 ∆r Mesh size in Sec. 1.8Cell volume, page 14 m Prandtl number, Pr = c p µ/λ, page Vector normal to the surface of the corresponding volume, directed towards the exterior of that volume, page Heat flux, q = λ∇T , page 12 W.m -Fluid velocity, page 11 m.s -1 A b Area of all the bubbles in the simulation volume, page 22 m 2 a I Interfacial area, page 24 m -1 c p Heat capacity, page 12 J.kg -1 .K -1 CF L Courant number, CF L = v∆t/∆x, page Thermal diffusivity D = λ/(ρc p ), page 18 m 2 .s -1 D b Diameter of the bubbles, page 21 m Eo Eötvös number, Eo = ∆ρgL 2 /σ, page 21 F o Fourier number, F o = D∆t/∆x 2 , page GravityMean indicator function over a cell, I k = 1/V V χ k dV , page 14 L j Length of the domain in the j th direction, page 18 m N b Number of bubbles, page 24 N j Number of cells in the j th direction, page 18 N u Nusselt number, page 22 p Pressure, page 20 Pa Re Reynolds number, Re = ρvL/µ, page TemperatureThe mean relative velocity between both phases, u r = u l -u v , page 21m.s -1 Indices and exponants * Next time value, determined with an integrated convection equation, page 40 (n) Timestep number, page Arithmetic Of At the face f of the cell, page Harmonic At Number of the cell, serve as a 1D index, used as coordinates, page In/of the phase k, ex: χ k is the indicator function of the phase k, page In Mean over a control volume, ex: T m = 1/V V T dV , page In the gas, page 14 tot Over the whole domain, page 24 Dirac function at the interface, page 20 ∇ • () Divergence, page 12 m -1

Figure 1 :

 1 Figure 1: Steps necessary to get a coherent formulation.

Figure 2 :

 2 Figure 2: Lagrangian mesh over the bubble and slices of the Eulerian mesh, coloured respectively by the curvature and integrated phase indicator function.

1. 3 .

 3 Temperature equation in the local one-fluid paradigm 230

3 (d) 1 /

 31 1.3 lead to the following one-fluid equation ρcp h , no λ jump

Figure 3 :

 3 Figure 3: Study of the impact of the use of arithmetic or harmonic means for 1/ρcp and λ. r = r/D b is the non-dimensional radius and Θ is the non-dimensional temperature. The temperature profiles are represented at different simulation times (s) both using TrioCFD-IJK and the semi-analytical solution (1D-FD presented by Panda et al. [59]).

Figure 4 :Figure 5 :

 45 Figure 4: Phase temperature and Nusselt number evolution in a homogeneous bubbly flow as a function of time.

1 512

 1 then leading to D b = 3 × 10 -4 m. The velocity is chosen to be identical to the relative velocity of the three-dimensional homogeneous bubbly flow, v = 0.2m/s. Numerical parameters. The number of cells in the length of the simulation is chosen to be as in 3D i. e., n tot = 512 and ∆x = 3.91 × 10 -5 m and the nondimensional numbers for numerical convergence (the Fourier number F o and the Courant number CF L) are imposed by the same constraints as in 3D simulations (see Tab. 4). n tot ∆x (m) L tot (m) α a I m -1 N b D b (m) v m s -

Figure 6 :

 6 Figure 6: Energy loss as a function of time in the equivalent 1D configuration. The simulation uses a RK3 time scheme and QUICK interpolation, as the 3D simulation.

Figure 7 :

 7 Figure 7: Temperature profile at t = 0.4s for the 1D simulation compared to the semianalytical reference. The simulation uses a RK3 time scheme and QUICK interpolation, as the 3D simulation.

2 . 2 with D b = 1 . 20 ×

 22120 10 -3 m.

Figure 8 :

 8 Figure 8: Comparison of the different convection schemes: effects on the energy.

Figure 9 :

 9 Figure 9: Comparison of the different convection schemes after a simulation time of 0.4s using an Euler time scheme.

T

  Ref, ∆x = 1.50 × 10 -5 TOF, ∆x = 1.50 × -5 TOF, ∆x = 2.00 × -5 TOF, ∆x = 3.00 × -5 TOF, ∆x = 3.91 × -5 TOF, ∆x = 5.00 × -5 TOF, ∆x = 7.02 × -5 TOF, ∆x = 9.01 × -5 (a) Temperature profile for different meshes.

3 × 10 1 4 × 10 1 6 × 10 1 D b /∆x 10 4 6 × 10 3 7 × 10 3 8 × 10 3 9 × 10 3 ∆E(J/m 3 )

 311433333 Regression slope: -0.38 (b) Energy loss as a function of the mesh size ∆x.

Figure 10 :

 10 Figure 10: Convergence of the simulation with convection and diffusion, with decreasing ∆x and ttot = 0.1s and ∆t = 1.02 × 10 -6 s, the Euler time scheme and a QUICK interpolation are used.

Figure 11 :

 11 Figure 11: Scheme of a cell illustrating the known values from which other values are interpolated.

  , and thus f (Y ) is an interpolation using Y values taken from several cells. TOF: temperature one-fluid, EOFm: energy one-fluid modified, TSP: temperature with single-phase fluxes, ESP: energy with single-phase fluxes, TSV: temperature separated volume integral.

495 3 . 1 . 1 .Figure 12 :

 31112 Figure 12: Interpolation scheme for the temperature and the heat flux at the interface.

Figure 13 :

 13 Figure 13: Values used for the interpolation of the temperature near the interface in one dimension with the fluid velocity positive from left to right. The green values are used to interpolate the red values.

 ( 31 ) 3 . 1 . 3 .

 31313 Temperature with Separation of the Volume integral over both phases 535 (TSV)

Figure 14 :

 14 Figure 14: Scheme of a cell illustrating the volume tracking of each phase in a conservative approach.

Figure 15 :

 15 Figure 15: Energy evolution with different temperature formulations.

3 . 1 Figure 16 :

 3116 Figure 16: Final temperature and diffusion flux (at t = 2s) with different temperature formulations. The fluid goes from left to right and the domain is actually much larger. An Euler time scheme and a QUICK interpolation (except at the interface vicinity) are used.

Figure 17 :

 17 Figure 17: Energy evolution with different energy formulations.

Figure 18 :× 10 1 D b /∆x 10 - 3 10 - 2 Error L 2 ESP,

 18110322 Figure 18: Final temperature and diffusion flux (at t = 2s) with different energy formulations.

Figure 19 :

 19 Figure 19: L 2 error compared to a non-convected reference solution.

  1.8) to tackle the interface condition issue that appears without this modification (see Figs.3a and 3b) in a pure diffusion case. But this improved formulation still does not conserve the thermal energy in the fluid. An error of 50% is obtained. Several new and original formulations are proposed and tested to improve the energy conservation and the scheme accuracy. We show that as long as the equation is written as the time derivative of the temperature (and not the energy), it does not hold energy conservation in a full convection-diffusion simulation (see Sec. 2.1). With a temperature evolution equation, conservation is not attained because of the heat capacity discontinuity between the phases which prevents the possibility to write the equation in a conservative way (only with fluxes).

  diffusive fluxes. Finally, a conservative formulation (Energy with Single-Phase fluxes, ESP) was proposed, with good interfacial temperature prediction (see Sec. 3.2.2) by discarding some approximations. Discontinuous properties are integrated over space and time on the wetted part of the face. Besides, the temperature used in the convective flux and in the heat diffusion flux are not calculated using the temperature of the mixed cell. They are reconstructed from the interface values and values located in the same phase. This new numerical method will be implemented and tested in three-dimensional cases in future works, opening a way to conservative simulations of bubbly flows with highly discontinuous thermal properties in the front-tracking framework.

Table 1 :

 1 Physical parameters

	Re b	Eo	ρ l ρv	µ l µv	Pr l	Pr v	ρ l cp l ρvcp v	α v
	400 0.636 13.4 26.2 0.384 0.00705	0.1	6%

Table 2 :

 2 Dimensionless quantities

	512 128 128 57	0.02	0.005	0.005	25.6

n x n y n z N b L x (m) L y (m) L z (m) D b /∆x

Table 3 :

 3 Numerical values

Table 4 :

 4

Numerical values used in the 1D simulation aimed to reproduce the energy loss observed in a 3D homogeneous bubbly flow.

Table 5 :

 5 Panel of methods used. In this table f refers to an arbitrary interpolation function

  should be noted that ρc p * is an intermediate variable. It is forgotten at the end of the time step as ρc p n+1 is recomputed from the front position and from the colour function I n+1 , in the classical way used for any other physical property. This replacement of ρc p *

	a by ρc p	n+1 a

Acknowledgement

The authors also would like to convey their sincere thanks to GENCI and the TGCC for providing the necessary computational resources for performing the DNS calculations. This work was granted access to the HPC resources of CINES under the allocations A0102b7712 and A0122b7712 made by GENCI. Also, we would like to acknowledge the CEA/DEN for the HPC resources allocated on IRENE super-computers at TGCC.


Keywords: . Data Clustering: Algorithms and Applications, 29 12/91 Data Classifications(s) Clustering: motivation Empirical procedures Automatic procedures Experiments Data Clustering: Algorithms and Applications, 29 13/91 Data Classifications(s) Clustering: motivation Empirical procedures Automatic procedures Experiments

Classifications(s) Clustering: motivation Empirical procedures Automatic procedures Experiments To go further Today's data (1/2) Today, it is easy to collect many features, so it favors data variety and/or mixed data missing data uncertainty (or interval data) Classifications(s) Clustering: motivation Empirical procedures Automatic procedures Experiments To go further Large data sets (n) 2 2

Mixed, missing, uncertain

Coding for data x

A set of n individuals x = {x 1 , . . . , xn} with x i a set of (possibly non-scalar) d variables

x i = {x i1 , . . . , x id } where x ij 2 X j
A n-uplet of individuals x = (x 1 , . . . , xn)

with x i a d-uplet of (possibly non-scalar) variables

x i = (x i1 , . . . , x id ) 2 X where X = X 1 ⇥ . . . . . . X d

We will pass from a coding to another, depending on the practical utility (useful for some calculus to have matrices or vectors for instance)

11/91

Genesis of "Big Data"

The Big Data phenomenon mainly originates in the increase of computer and digital resources at an ever lower cost

More data for what?

Opportunity to improve accuracy of traditional questionings

Here is just illustrated the e↵ect of n

In further lessons will be illustrated the e↵ect of d (be patient) Supervised classification (1/3)

Data: learning dataset D = {x O , z} n individuals: x = {x1, . . . , xn} = {x O , x M }, each x i 2 X Observed individuals x O Missing individuals x M Partition in K groups G1, . . . , G K : z = (z1, . . . , zn) 6 , z i = (z i1 , . . . , z iK ) 0 x i 2 G k , z ih = I {h=k} Aim: estimation of an allocation rule r from D r : X ! {1, . . . , K } x O n+1 7 ! r (x O n+1 ).

Other possible coding for z

Sometimes it can be more convenient in formula to code z = (z 1 , . . . , zn) a where

x i 2 G k , z i = k
a Sometimes it would be more convenient to use a set also: z = {z 1 , . . . , zn}

6 Sometimes it would be more convenient to use a set also: z = {z 1 , . Semi-supervised classification (1/3) We will gradually answer along the lessons: Partition: retrieve less that 19 999$ (group of "low income"), between 20 000$ and 39 999$ group of "average income"), more than 40 000$ (group of "high income")

Data: learning dataset D = {x O , z O } n individuals: x = {x1, . . . , xn} = {x O , x M } belonging to a space X Observed individuals x O Missing individuals x M Partition: z = {z1, . . . , zn} = {z O , z M } Observed partition z O Missing partition z M Aim: estimation of an allocation rule r from D r : X ! {1, . . . , K } x O n+1 7 ! r (x O n+1 ). Idea: x is cheaper than z so #z M #z O 22/91
-1 -0.5 0 0.5 1 1.5 2 2.5 -1.5 -1 -0.5 0 0.5 1 1.5 1st MCA axis 2nd MCA axis -1 -0.5 0 0.5 1 1.5 2 2.5 -1.5 -1 -0.5 0 0.

Four main clustering structures: partition

Each object belongs to exactly one cluster Definition: a set of K non-empty parts of x: P = (G 1 , . . . , G K ):

For all k 6 = k 0 , G k \ G k 0 = ; G1 [ . . . [ G K = x Example: x = {a, b, c, d, e} P = {{a, b} | {z } G 1 , {c, d, e} | {z } G 2 } Notation: z = (z 1 , . . . ,

Four main clustering structures: partition with outliers

Each object belongs to one cluster or less Definition: a set of K non-empty parts of x: P = (G 1 , . . . , G K ):

For all k 6 = k 0 , G k \ G k 0 = ; G1 [ . . . [ G K ⇢ x Example: x = {a, b, c, d, e} P = {{a, b} | {z } G 1 , {d, e} | {z } G 2 },

Four main clustering structures: multiclass partition

Each object belongs to one cluster or more Definition: a set of K non-empty parts of x: P = (G 1 , . . . , G K ):

G1 [ . . . [ G K = x Example: x = {a, b, c, d, e} P = {{a, b, c} | {z } G 1 , {c, d, e} | {z } G 2 } Not considered
◆ : H ! R + is an index of a hierarchy H if ◆ is a non-decreasing mapping: (G1, G2) 2 H 2 , G1 ⇢ G2, G1 6 = G2 ) ◆(G2) < ◆(G2) For all x i 2 x, ◆({x i }) = 0
The couple (H, ◆) is an indexed hierarchy and can be visualized by a dendrogram Example (continuated):

H = {{a} |{z} 0 , {b} |{z} 0 , {c} |{z} 0 , {d} |{z} 0 , {e} |{z} 0 , {d, e} | {z } 1 , {a, b} | {z } 2 , {c, d, e} | {z } 2.5 , {a, b, c, d, e} | {z } 5 , } 44

Clustering is the cluster building process

According to JSTOR, data clustering first appeared in the title of a 1954 article dealing with anthropological data Need to be automatic (algorithms) for complex data: mixed features, large data sets, high-dimensional data. . .

X = R d
Euclidean distance (or L 1 norm) with metric M:

M (x 1 , x 2 ) = kx 1 x 2 k M = p (x 1 x 2 ) 0 M(x 1 x 2 )
with a metric M: symmetric and positive definite d ⇥ d matrix Manhattan distance or city-block distance or L 1 norm:

(x 1 , x 2 ) = d X j=1 kx 1j x 2j k . . . 53/
X = {0, 1} d thus x ij 2 {0, 1}
The Hamming distance counts the dimensions on which both vectors are di↵erent:

(x 1 , x 2 ) = d X j=1 I {x 1j 6 =x 2j }
The matching coe cient counts dimensions on which both vectors are non-zero:

(x 1 , x 2 ) = d X j=1 I {x 1j =x 2j =1}
. . .

Similarities, dissimilarities, distances: similarities proposals

It is easy to construct similarities from dissimilarities or distances A typical example:

s(x 1 , x 2 ) = exp( (x 1 , x 2 )) Thus smax = 1
The Gaussian similarity is an important case ( > 0):

s(x 1 , x 2 ) = exp ✓ kx 1 x 2 k 2 2 2 ◆ 56/

Within-cluster inertia criterion

Select the partition z minimizing the criterion

W M (z) = n X i=1 K X k=1 z ik kx i xk k 2 M Look for

Between-cluster inertia criterion

Select the partition z maximizing the criterion

B M (z) = K X k=1 n k kx xk k 2 M
Look for clusters far from each other

x the mean (or center) of the whole data set x:

x = 1 n n X i=1 x i It is in fact equivalent to minimize W M (z) since B M (z) + W M (z) =

Fuzzy within-cluster inertia criterion

Select the fuzzy partition c minimizing the criterion

W fuzzy M,m (c) = n X i=1 K X k=1 c m ik kx i xk k 2 M Fuzzy partition c = {c 1 , . . . , c K }, c k = (c 1k , . . . , c nk ), c ik 2 [0, 1] defined by c ik = 1 P K k 0 =1 ⇣ kx i xk k M kx i xk 0 k M ⌘ 2 m 1
m 1 determines the level of cluster fuzziness: m large produces fuzzier clusters; m = 1 (in the limit) produces hard clustering (m converges to 0 or 1) xk is the fuzzy mean (or center) of the kth cluster Feature space: a non-linear mapping to a higher-dimensional space

: X 7 ! Y x i ! (x i )
Expected property: linearly separable clusters in the feature space Di culty: define

The "kernel trick": not necessary to define , define just an inner product

(x i , x i 0 ) = (x i ) 0 (x i 0 )
 is so-called a kernel Kernel clustering: fuzzy within-cluster rewriting W fuzzy with the data set (x) = { (x 1 ), . . . , (xn)} (we drop M and m)

W fuzzy (c) = n X i=1 K X k=1 c m ik k (x i ) ¯ k (x)k 2 with ¯ k (x) = 1 n k P n i=1 c m ik (x i ) P n i=1 c m ik and c ik = 1 P K k 0 =1 ⇣ k (x i ) ¯ k (x)k k (x i ) ¯ k 0 (x)k ⌘ 2 m 1
After some algebra, only appears inside an inner product since:

k (x i ) ¯ k (x)k 2 = (x i ) 0 (x i ) 2 P n i=1 c m ik (x i ) 0 (x i 0 ) P n i 0 =1 c m i 0 k + P n i 0 =1 P n i 00 =1 c m i 0 k c m i 00 k (x i 0 ) 0 (x i 00 ) ( P n i 0 =1 c m i 0 k ) 2 = (x i , x i ) 2 P n i=1 c m ik (x i , x i 0 ) P n i 0 =1 c m i 0 k + P n i 0 =1 P n i 00 =1 c m i 0 k c m i 00 k (x i 0 , x i 00 ) ( P n i 0 =1 c m i 0 k ) 2
Consequently, need only to define  for performing the (fuzzy) K -means. . .

(x i , x i 0 ) = (x 0 i x i 0 + a) b Gaussian (x i , x i 0 ) = exp( kx i x i 0 k/(2 2 )) Sigmoid (x i , x i 0 ) = tanh(ax 0 i x i 0 + b)
S ii 0 = s ii 0 if i 6 = i 0 , S ii = 0
Similarities can be viewed as proximities between individuals (nodes) in a graph

If some disconnected sub-graphs, S can be reordered as block diagonal

We "can see" the clusters thus clustering is expected to be easy

Properties of L

1 For every vector u = (u 1 , . . . , u 0 n ) 2 R n , we have

u 0 Lu = 1 2 n X i=1 n X i 0 =1 s ii 0 (u i u i 0 ) 2 2 L is
y i = 1 G 1 (i), . . . , 1 G K (i)
From the previous property, y i s in the same clusters are very similar. . . Spectral clustering: some last elements

Other Laplacian choices are possible, as the normalized one Single linkage criterion

L norm = D 1/2 LD 1/2 = I D 1/2 SD 1/2
(G 1 , G 2 ) = min{ (x i , x i 0 ) : x i 2 G 1 , x i 0 2 G 2 } 74/
Complete linkage criterion

(G 1 , G 2 ) = max{ (x i , x i 0 ) : x i 2 G 1 , x i 0 2 G 2 } 75/
Average linkage criterion

(G 1 , G 2 ) = 1 n 1 n 2 X x i 2G 1 X x i 0 2G 2 (x i , x i 0 ) 76/91 Data Classifications(s)
Clustering:

Ward linkage criterion

(G 1 , G 2 ) = n 1 n 2 n 1 + n 2 2 (µ 1 , µ 2 )
P = {G 1 , . . . , G K , G K +1 , G K +2 } W (P) = K X k=1 X x i 2G k (x i , µ k ) + X x i 2G K +1 (x i , µ K +1 ) + X x i 2G K +2 (x i , µ K +2 )
Partition at a the next iteration:

P + = {G 1 , . . . , G K , G K +1 [ G K +2 | {z } G + K +1 } W (P + ) = K X k=1 X x i 2G k (x i , µ k ) + X x i 2G K +1 [G K +2 (x i , µ + K +1 )
After some algebra and using the Koenig-Huygens theorem, we obtain Semi-supervised clustering 20

W (P + ) W (P) = n K +1 n K +2 n K +1 + n K +2 (µ K +1 , µ K +2 )
The user has to provide any external information he has on the partition Pair-wise constraints:

  If we have a look at the true partition, indeed visualizing was not enough. . . Need to define more "automatic" methods acting directly in the space X

  Number of possible partitions z: #{z} ⇡ K n K ! For instance : n = 40, K = 3, #{z} ⇡ 2.10 18 (two billion of billion), thus about 64 000 years for a computer calculating one million of partitions per second Thus impossible to minimize directly W M (z) But equivalent to minimize on z and µ = {µ 1 , . . . , µ K } with µ k 2 R d Alternating optimization between the partition z and the center µ of clusters Start from a random µ among x (avoid to start from a random z) Each iteration decreases the criterion WM (z, µ) (so also W M (z)) Converge into a finite (and often low) number of iterations Memory complexity O((n + K )d), time complexity O(nb.iter.Kdn)

  = P n k=1 c ik indicates the fuzzy number of individuals in cluster k Optimize W fuzzy M,m (c) (through W fuzzy M,m (c, µ)) by the fuzzy K -) space: non-linearly separable clusters

  symmetric and positive semi-definite 3 The smallest eigenvalue of L is 0, the corresponding eigenvector is the constant one vector 1 = (1, . . . , 1) 0 2 R n 4 L has n non-negative, real-valued eigenvalues 0 = 1  2  . . .  n Spectral clustering: specific spectral structure of L Linking clusters and L 1 The multiplicity K of the eigenvalue 0 of L equals the number of clusters G 1 , . . . , G K (or the connected components in the graph) 2 The eigenspace of eigenvalue 0 is spanned by the indicator vectors 1 G 1 , . . . , 1 G K of those components, where 1 G k (i) = 1 if x i 2 G k , zero otherwise [A. Singh (2010). Spectral Clustering. Machine Learning.] Spectral clustering: eigenvectors as new data to cluster Define a new data set y = {y 1 , . . . , yn} 2 R K where

  . . . whereas y i s in di↵erent clusters are very di↵erent Thus just run a K -means (for instance) to finish the clustering job! [A.K. Jain (2008). Data Clustering: 50 Years Beyond K-Means]

  obtained by cuting the dendrogram A dissimilarity matrix between pairs of individuals is enough Memory complexity O(n 2 ln n), space complexity O(n 2 ) (do not scale well. . . ) Suboptimal optimisation of W M (•) (see next slide) and K -means Two successive Ward hierarchical algorithm iterations minimize W Partition at a given iteration:

A

  must-link constraint specifies that the point pair connected by the constraint belong to the same cluster A cannot-link constraint specifies that the point pair connected by the constraint do not belong to the same cluster Attempts to derive constraints from domain ontology and other external sources into clustering algorithms include the usage of WordNet ontology, gene ontology, Wikipedia, etc. to guide clustering solutions 20 O. Chapelle et al. (2006), A.K. Jain (2008). Data Clustering: 50 Years Beyond K-Means. quite recent: blogs, Web pages, retail chain, credit card transaction streams, network packets received by a router and stock market, etc. As the data gets modified, clustering must be updated accordingly: ability to detect emerging clusters, etc. Often all data cannot be stored on a disk This imposes additional requirements to traditional clustering algorithms to rapidly process and summarize the massive amount of continuously arriving data Data stream clustering a significant challenge since they are expected to involved singlevision unifiée pour une utilisation éclairée Evaluation d'une méthode de clustering 91/91

  But more complex with mixed, missing or uncertain data. . .

	Data Data Data Data Data Data Data Data		Classifications(s) Classifications(s) Classifications(s) Classifications(s) Classifications(s) Classifications(s) Classifications(s) Classifications(s)			Clustering: motivation Clustering: motivation Clustering: motivation Clustering: motivation Clustering: motivation Clustering: motivation Clustering: motivation Clustering: motivation	Empirical procedures Empirical procedures Empirical procedures Empirical procedures Empirical procedures Empirical procedures Empirical procedures Empirical procedures			Automatic procedures Automatic procedures Automatic procedures Automatic procedures Automatic procedures Automatic procedures Automatic procedures Automatic procedures			Experiments Experiments Experiments Experiments Experiments Experiments Experiments Experiments	To go further To go further To go further To go further To go further To go further To go further To go further
									Semi-supervised classification (2/3) Semi-supervised classification (3/3) Unsupervised classification (1/3) Unsupervised classification (2/3) Unsupervised classification (3/3) Traditional solutions (1/3) Traditional solutions (2/3) Traditional solutions (3/3)
	6															6
	5	No mixed, missing or uncertain data:			5
	-4 -2 0 2 4 -3 -2 -1 0 1 2 3 4 -4 -5	Mixed, missing, uncertain Individuals x O ? 0.5 0.3 0.1 green red 0.3 0.6 {red,green} 0.9 [0.25 0.45] red # # # continuous continuous categorical •? Aim: estimation of the partition z 7 and the number of groups K Partition z O , 5 0 ? ? , 3 1 0 0 , 3 ? ? ? , ? 0 0 1 , # integer -5 -4 -3 -2 -1 0 1 2 3 4 ! Also known as: clustering Mixed, missing, uncertain Individuals x O Partition z O ? 0.5 red 5 ? ? ? 0.3 0.1 green 3 ? ? ? 0.3 0.6 {red,green} 3 ? ? ? 0.9 [0.25 0.45] red ? ? ? ? # # # -2 0 2 4 Supervised classification 9 Two main model-based frameworks 8 Generative models: linear/quadratic discriminant analysis Generative models Predictive models: logistic regression, support vector machines (SVM), k nearest Group G2 or G3 G1 ??? G3 , Group , ??? , ??? , ??? , ??? Model p(x, z) Thus direct model for p(x) = P neighbourhood, classification trees. . . z p(x, z) Easy to take into account some missing z and x Semi-supervised classification 10 ! Generative models: mixture models Predictive models Predictive models: low density separation (transductive SVM), graph-based methods. . . Model p(z|x) or sometimes 1 {p(z|x)>1/2} or also ranking on p(z|x) Avoid asumptions on p(x), thus avoids associated error model Unsupervised classification 11 # continuous continuous categorical -4 di cult to take into account some missing z and x Generative models: k-means like criteria, hierarchical clustering, mixture models integer Predictive models: -
	-6 -6	-6	-4	-4	-2	-2	0	0	2	2	4	4	6	6	8	-6 -6	-4	-4	-2	-2	0	0	2	2	4	4	6	6	8
					{x O , z O } and x O n+1 x O									r and ẑn+1 {x O , ẑ}
			23/91 7 We will see other clustering structures further, here to fix ideas with the main one 24/91 25/91 26/91 27/91 8 Sometimes presented as distance-based instead of model-based: see some links in my 3th lesson 28/91 29/91

Data: learning dataset D = x O , so z O = ; 9 Govaert et al., Data Analysis, Chap.6, 2009 10 Chapelle et al., Semi-supervised learning, 2006 11 Govaert et al., Data Analysis, Chap.7-9, 2009

Missing/uncertain data: multiple imputation is possible but it should ideally take into account the classification purpose at hand Mixed data: some heuristic methods with recoding

  QuestionsDoctors say two stages of the disease are well separated with the variables We forget the classes (Stages of the disease) to retrieve clustering by visualizing

	Data Data Data	Classifications(s) Classifications(s) Classifications(s)	Clustering: motivation Clustering: motivation Clustering: motivation	Empirical procedures Empirical procedures Empirical procedures	Automatic procedures Automatic procedures Automatic procedures	Experiments Experiments Experiments	To go further To go further To go further
			Interdisciplinary endeavour Limit of visualization (1/3) Limit of visualization (2/3)
		Perform PCA for continuous variables, MCA for categorical ones
		Prostate cancer data 14 Discard missing values since not convenient for traditional PCA and MCA
		Individuals: 506 patients with prostatic cancer grouped on clinical criteria into
		two Stages 3 and 4 of the disease
				46/91	45/91		/91

Variables: d = 12 pre-trial variates were measured on each patient, composed by eight continuous variables (age, weight, systolic blood pressure, diastolic blood pressure, serum haemoglobin, size of primary tumour, index of tumour stage and histolic grade, serum prostatic acid phosphatase) and four categorical variables with various numbers of levels (performance rating, cardiovascular disease history, electrocardiogram code, bone metastases) Some missing data: 62 missing values (⇡ 1%)

14

Byar DP, Green SB (1980): Bulletin Cancer, Paris 67:477-488

  Clustering of clustering algorithms 15 Jain et al. (2004) hierarchical clustered 35 di↵erent clustering algorithms into 5 groups based on their partitions on 12 di↵erent datasets. It is not surprising to see that the related algorithms are clustered together. For a visualization of the similarity between the algorithms, the 35 algorithms are also embedded in a two-dimensional space obtained from the 35x35 similarity matrix. Put in the same cluster individuals x 1 and x 2 which are similar Put in di↵erent clusters individuals x 1 and x 2 which are dissimilar

	Data Data Data Data	Classifications(s) Classifications(s) Classifications(s) Classifications(s)	Clustering: motivation Clustering: motivation Clustering: motivation Clustering: motivation	Empirical procedures Empirical procedures Empirical procedures Empirical procedures	Automatic procedures Automatic procedures Automatic procedures Automatic procedures	Experiments Experiments Experiments Experiments	To go further To go further To go further To go further
		Outline Dissimilarities, distances, similarities: definitions Similarities, dissimilarities, distances: examples for quantitative data
		Data General idea
		Classifications(s)
		Clustering: motivation
		Empirical procedures
		Automatic procedures
		Experiments For any x1 2 X , (x1, x1) = smax with smax Symmetry: for any (x1, x2) 2 X 2 , s(x1, x2) = s(x2, x1) s(x1, x2) for any (x1, x2) 2 X 2 , x1 6 = x2
		To go further Often dissimilarities (or similarities) are enough for clustering
		50/91	49/91	51/91	52/91

15 A.K. Jain (2008). Data Clustering: 50 Years Beyond K-Means. Definition of a dissimilarity : X ⇥ X ! R + Separation: for any (x1, x2) 2 X 2 , (x1, x2) = 0 , x1 = x2 Symmetry: for any (x1, x2) 2 X 2 , (x1, x2) = (x2, x1) Definition of a distance : X ⇥ X ! R + : it is a dissimilarity with a new property Triangular inegality: for any (x1, x2, x3) 2 X 3 , (x1, x3)  (x1, x2) + (x2, x3)

Definition of a similarity (or a nity) s : X ⇥ X ! R +

  Spectral clustering: the idea 17Build a similarity matrix S n ⇥ n from x:

	Data	Classifications(s)	Clustering: motivation	Empirical procedures	Automatic procedures	Experiments	To go further
					66/91

  , . . . , dn) is the diagonal degree matrix with d i = P n i 0 =1 S ii 0 Have in mind S if easier ; the structure of L is similar

	Data	Classifications(s)	Clustering: motivation	Empirical procedures	Automatic procedures	Experiments	To go further
			Spectral clustering: Laplacian matrix
		For technical reasons, we study (unnormalized graph) Laplacian L instead of S
				L = D S
		where D = diag(d 1
				67/91

17 Figures from [A. Jain at al.. Data Clustering: A Review.]

  Popularity of K -means and hierarchical clusteringEven K -means was first proposed over 50 years ago, it is still one of the most widely used algorithms for clustering for several reasons: ease of implementation, simplicity, e ciency, empirical success. . . and model-based interpretation (see later)

	Data Data Data Data Data	Classifications(s) Classifications(s) Classifications(s) Classifications(s) Classifications(s)	Clustering: motivation Clustering: motivation Clustering: motivation Clustering: motivation Clustering: motivation	Empirical procedures Empirical procedures Empirical procedures Empirical procedures Empirical procedures	Automatic procedures Automatic procedures Automatic procedures Automatic procedures Automatic procedures	Experiments Experiments Experiments Experiments Experiments	To go further To go further To go further To go further To go further
				Outline Outline Packages
		Data Data
		Classifications(s) Classifications(s)
		Clustering: motivation Clustering: motivation R (statistical community): have a look at the CRAN library
		Python (machine learning community): have a look at the scikit-learn library
		Empirical procedures Empirical procedures
		Automatic procedures Automatic procedures
		Experiments Experiments
		To go further To go further
		81/91 87/91			79/91 80/91	88/91

Thus, it is like performing K -means under (partial) partition P constraints

[J. Suykens and C. Alzate (2011). Kernel spectral clustering: model representations, sparsity and out-of-sample extensions.]

Data set used for illustration 82/91