Pierre-Henri Maire
email: maire@celia.u-bordeaux1.fr

Jérôme Breil
email: breil@celia.u-bordeaux1.fr

A high-order finite volume cell-centered scheme for anisotropic diffusion on two-dimensional unstructured grids

Keywords: Anisotropic diffusion, isotropic diffusion, cell-centered scheme, high-order finite volume method, two-dimensional unstructured grid, cylindrical geometry PACS: 47.11.Df, 47.10.ab, 47.40.Nm 2000 MSC: 76N15, 65M06

In this paper, we describe a high-order cell-centered finite volume method for solving anisotropic diffusion on two-dimensional unstructured grids. The resulting numerical scheme, named CCLAD (Cell-Centered LAgrangian Diffusion), is characterized by a local stencil and cell-centered unknowns. It is devoted to the resolution of diffusion equation on distorted grids in the context of Lagrangian hydrodynamics wherein a strong coupling occurs between gas dynamics and diffusion. The space discretization relies on the introduction of two half-edge normal fluxes and two half-edge temperatures per cell interface using the partition of each cell into sub-cells. For each cell, the two half-edge normal fluxes attached to a node are expressed in terms of the half-edge temperatures impinging at this node and the cell-centered temperature. This local flux approximation can be derived through the use of either a sub-cell variational formulation or a finite difference approximation, leading to the two variants CCLADS and CCLADNS. The elimination of the half-edge temperatures is performed locally at each node by solving a small linear system which is obtained by enforcing the continuity condition of the normal heat flux across sub-cell interface impinging at the node. The accuracy and the robustness of the present scheme is assessed by means of various numerical test cases.

Introduction

We present a finite volume scheme to solve diffusion equations on two-dimensional unstructured grids. This scheme is the extension to anisotropic diffusion of the work initially described in [START_REF] Breil | A cell-centered diffusion scheme on two-dimensional unstructured meshes[END_REF]. In deriving this numerical method, we aim at developing numerical modeling of physical phenomena encountered in plasma physics. More precisely, we are concerned by heat transfer within laserheated plasma flows such as those obtained in the domain of direct drive Inertial Confinement Fusion [START_REF] Atzeni | The physics of inertial fusion[END_REF]. Let us emphasize that for such flows, the heat conduction equation is strongly coupled to the gas dynamics equations describing the plasma motion. These latter equations, otherwise called Lagrangian hydrodynamics equations, are solved employing a Lagrangian numerical method wherein the computational grid is moving with the fluid. In addition, the thermodynamic variables, that is, the density, the pressure, the specific internal energy and the temperature are located at the cell center, refer to [START_REF] Maire | A high-order one-step sub-cell force-based discretization for cell-centered lagrangian hydrodynamics on polygonal grids[END_REF][START_REF] Maire | A unified sub-cell force-based discretization for cell-centered lagrangian hydrodynamics on polygonal grids[END_REF][START_REF] Maire | A high-order cell-centered Lagrangian scheme for two-dimensional compressible fluid flows on unstructured meshes[END_REF]. Bearing in mind this coupling between Lagrangian hydrodynamics and diffusion, leads to the following requirements concerning the diffusion scheme under consideration:

• It should be a finite volume scheme wherein the primary unknown, i.e., the temperature, is located at the cell center.

• It should be sufficiently accurate and robust scheme to handle highly distorted grids which result from the fluid motion.

The present work aims at describing a finite volume scheme which fulfils the previous requirements.

On that account, we denominate it using the acronym CCLAD which stands for Cell-Centered LAgrangian Diffusion. Before describing the main features of CCLAD scheme, let us briefly give an overview of the existing cell-centered diffusion schemes. It is well known, see [START_REF] Eymard | Finite Volume methods[END_REF], that the standard finite volume algorithms, such as the five-point scheme, behave poorly on highly skewed quadrilateral grids. In this situation, the five-point scheme produces a numerical solution wherein the diffusion front is aligned with the grid distortions. This undesirable behavior is due to the crude finite difference approximation used for discretization of the face fluxes.

Kershaw, in his pioneering work [START_REF] Kershaw | Differencing of the Diffusion Equation in Lagrangian Hydrodynamic Codes[END_REF] has proposed a nine-point scheme on structured quadrilateral grids, which partially resolves the above mentioned difficulties. His scheme consists of a cell-centered variational method based on a smooth mapping between the logical mesh coordinates and the spatial coordinates. This algorithm reduces to the classical five-point scheme on an orthogonal grid. In addition, it leads to a diffusion matrix, which is symmetric positive definite. Although this method is restricted to structured quadrilateral grids, it has been successfully used in many Lagrangian codes devoted to the numerical simulation of Inertial Confinement Fusion, see for instance [START_REF] Moses | Radiation diffusion in DRACO using Kershaw difference scheme[END_REF]. Note that Kershaw's scheme has been recently extended to unstructured grids [START_REF] Rebourcet | Some remarks on Kershaw's legacy diffusion scheme[END_REF] and to anisotropic diffusion [START_REF] Rebourcet | An old scheme adapted to new problems: Kershaw's Finite Differences diffusion scheme adapted to new problems[END_REF]. However, the underlying assumption of a smooth mapping used by Kershaw is too restrictive. As it has been shown in [START_REF] Morel | A Cell-Centered Lagrangian-Mesh Diffusion Differencing Scheme[END_REF], a mesh refinement with Kershaw's scheme does not give a convergent solution unless the mesh becomes smooth as it is refined. Moreover, it appears that the normal flux continuity across cell interfaces is not ensured.

These drawbacks have motivated the work of Morel and his co-authors. In [START_REF] Morel | A Cell-Centered Lagrangian-Mesh Diffusion Differencing Scheme[END_REF] they developed a cell-centered diffusion scheme, which treats rigorously material discontinuities and gives a second order accuracy regardless of the smoothness of the mesh. However, this scheme has two disadvantages: there are cell-edge unknowns in addition to the cell-centered unknowns and the diffusion matrix is asymmetric.

A significant improvement was provided by Shashkov and Steinberg. In [START_REF] Shashkov | Solving Diffusion Equations with Rough Coefficients in Rough Grids[END_REF], [START_REF] Shashkov | Support-Operator Finite-Difference Algorithms for General Elliptic Problems[END_REF] they derived an algorithm using the Support Operators Method (SOM), also named mimetic finite difference method. This method, see [START_REF] Shashkov | Conservative Finite-Difference Methods on General Grids[END_REF], constructs discrete analogs of the divergence and flux operators that satisfy discrete analogs of important integral identities relating the continuum operators. By this way, the discrete flux operator is the negative adjoint of the discrete divergence in an inner scalar product weighted by the inverse conductivity. This SOM diffusion scheme gives the second order accuracy on both smooth and non smooth meshes either with or without material discontinuities. It has a non local stencil and a dense symmetric positive definite matrix representation for the diffusion operators. The introduction of both cell-centered and face-centered unknowns in [START_REF] Morel | A local support-operators diffusion discretization scheme for quadrilateral r-z meshes[END_REF] leads to a variant of this scheme, which has a local stencil. Many extensions of this algorithm have been recently developed. One can find in [START_REF] Hyman | Mimetic finite difference methods for diffusion equations[END_REF][START_REF] Kuznetsov | The mimetic finite difference method on polygonal meshes for diffusion-type problems[END_REF][START_REF] Lipnikov | Mimetic finite difference methods for diffusion equations on non-othogonal non-conformal meshes[END_REF][START_REF] Gyrya | High-order mimetic finite difference method for diffusion problems on polygonal meshes[END_REF][START_REF] Lipnikov | The mimetic finite difference discretization of diffusion problem on unstructured polyhedral meshes[END_REF], developments that take into account nonisotropic materials, polygonal and non-conformal meshes, and also polyhedral meshes. This method has been also recently applied to solve the three-dimensional diffusion equation in multi-material domains containing mixed cells [START_REF] Garimella | Solution of the diffusion equation in multi-material domains by sub-division of elements along reconstructed interfaces[END_REF].

In [START_REF] Moulton | Workshop: Compatible Spatial Discretizations for Partial Differential Equations[END_REF], the authors present mimetic preconditioners for mixed discretizations of the diffusion equation. In this paper, SOM is used with two fluxes per edge in order to construct the local flux discretization. Recently, Lipnikov, Shashkov and Yotov developed a local flux mimetic finite difference method in [START_REF] Lipnikov | Local flux mimetic finite difference methods[END_REF][START_REF] Lipnikov | Local flux mimetic finite difference methods[END_REF], which is very similar to our derivation. They also use two degrees of freedom per edge to approximate the flux. They obtained a symmetric, cell-centered finite difference scheme. Moreover, they demonstrated theoretically the second-order convergence for the temperature in the case of simplicial meshes.

In [START_REF] Aavatsmark | Discretization on unstructured grids for inhomogeneous, anistropic media. Part I: derivation of the methods[END_REF], [START_REF] Aavatsmark | Discretization on unstructured grids for inhomogeneous, anistropic media. Part II: discussion and numerical results[END_REF] and [START_REF] Aavatsmark | An introduction to multipoint flux approximations for quadrilateral grids[END_REF] Aavatsmark and co-authors have proposed an alternative approach named Multi-Point Flux Approximation (MPFA). There, the flux is approximated by a multi-point flux expression based on transmissibility coefficients. These coefficients are computed using continuity of the flux and the temperature across the cell interfaces. This method has only cell-centered unknowns and a local stencil. In [START_REF] Klausen | Relationships among some locally conservative discretization methods wich handle discontinuous coefficients[END_REF], Klausen and Russel present the relationships between the Mixed Finite Element Method (MFEM), the Control Volume Mixed Finite Element Method (CVMFEM), the SOM and the MPFA. The latter can be applied in the physical space to quadrilateral and to unstructured grids. For quadrilaterals, which are not parallelograms the MPFA provides a second order scheme [START_REF] Aavatsmark | Numerical convergence of the MPFA O-method for general quadrilateral grids in two and three dimensions[END_REF] but the diffusion matrix is non-symmetric. In [START_REF] Aavvatsmark | Convergence of a symmetric MPFA method on quadrilateral grids[END_REF], the authors develop a MPFA method for quadrilateral grids in the reference space and its relationship to the MFEM. This approach yields a system of equations with a symmetric matrix. It shows a second-order convergence on smooth distorted grids. However for rough grids the reference space method suffers from a reduction or loss of convergence.

The relation between the finite volume and the MFEM is also studied by Thomas and Trujillo in [START_REF] Thomas | Mixed finite volume methods[END_REF]. These authors use a sub-triangulation, identical to the one used in the present paper. They are also able to eliminate auxiliary unknowns. However, the degrees of freedom for the scalar unknown are located on the vertices of the mesh.

Another class of finite volumes schemes for solving anisotropic diffusion equations on two and three-dimensional unstructured grids has been developed by Hermeline [START_REF] Hermeline | A finite volume method for the approximation of diffusion operators on distorded meshes[END_REF][START_REF] Hermeline | A finite volume method for approximating 3d diffusion operators on general meshes[END_REF]. This type of scheme, termed as Discrete Duality Finite Volume (DDFV) scheme, also arises from the construction of discrete analogs of the divergence and flux operators which fulfil the discrete counterpart of vector calculus identities. However, this method requires to solve the diffusion equation not only over the primal grid but also over a dual grid. Namely, there are both cell-centered and vertex centered unknowns. Compared to a classical cell-centered finite volume scheme, the DDFV method necessitates twice as much degrees of freedom over quadrangular grids. Let us point out that this method might be difficult to use in the context of the coupling between Lagrangian hydrodynamics and diffusion.

We also mention the papers [START_REF] Crumpton | Discretisation and multigrid solution of elliptic equations with mixed derivatives terms and strongly discontinuous coefficients[END_REF] and [START_REF] Edwards | Finite volume discretization with imposed flux continuity for the general tensor pressure equation[END_REF] where local expressions for the diffusion flux has been derived in the context of finite volume scheme for the diffusion equation. In the same framework, Le Potier has derived a local flux approximation in [START_REF] Potier | Schéma volumes finis pour des opérateurs de diffusion fortement anisotropes sur des maillages non structurés[END_REF] and [START_REF] Potier | A finite volume method for the approximation of highly anisotropic diffusion operators on unstructured meshes[END_REF] which is very similar to the MPFA symmetric method and to our method.

We conclude this non-exhaustive review by quoting recent works concerning the development of monotone finite volume methods for diffusion equation which preserve the positivity of the solution. It is well known that high-order linear methods, such as the multi-point flux approximation, mixed finite element and mimetic finite difference methods, are not monotone on strongly anisotropic meshes or for diffusion problems with strongly anisotropic coefficients. On the other hand, the finite volume method with linear two-point flux approximation is monotone but not even firstorder accurate in these cases. This flaw has been corrected by constructing monotone finite volume schemes wherein the discretization is based on a non-linear two-point flux approximation [START_REF] Potier | Schéma volumes finis monotone pour des opérateurs de diffusion fortement anisotropes sur des maillages de triangles non structurés[END_REF][START_REF] Lipnikov | Monotone finite volume schemes for diffusion equations on unstructured triangular and shape-regular polygonal meshes[END_REF][START_REF] Lipnikov | Interpolation-free monotone finite volume method for diffusion equations on polygonal meshes[END_REF]. The drawback of these methods lies in the fact that they require the solution of a global non-linear problem by means of an iterative procedure such as a fixed point algorithm, even in the case of a linear diffusion equation.

Finally, it seems that the diffusion scheme derived from the SOM has the best combination of ideal properties of any previous finite-difference scheme. The only drawback lies in the fact that there are both cell-centered and face-centered unknowns. For instance, using a bi-dimensional computational domain D paved with C D triangular (quadrangular) cells one has asymptotically 2.5 C D (3 C D ) unknowns. In addition, the treatment of the supplementary face-centered unknowns leads to a more complicated algorithm than usual when coupling the diffusion scheme to hydrodynamics.

This disadvantage has motivated us to propose the CCLAD scheme, which retains as well as possible the good properties of the SOM diffusion scheme.

The main feature of CCLAD scheme lies in the introduction of two half-edge normal fluxes and two half-edge temperatures per edge. For each cell, the two half-edge normal fluxes impinging on a node are expressed in terms of the two-half edge temperatures and the cell-centered temperature. This discretization of the half-edge fluxes is derived using either a local variational formulation written for each cell corner or a finite difference approximation based on a Taylor expansion. The former numerical expression leads to the version of CCLAD scheme denoted CCLADS whereas the latter numerical expression yields the version called CCLADNS. For both versions, the halfedge temperatures are eliminated locally by invoking the continuity of the temperature and the normal flux across each edge. This elimination procedure of the half-edge temperatures in terms of the cell-centered temperature is achieved by solving a cyclic tridiagonal linear system at each node. Collecting the contribution of each node allows to construct easily the global sparse diffusion matrix. Let us emphasize that the node-based construction of CCLAD scheme is particularly well adapted to cope with unstructured grids and also provide a natural treatment of the boundary conditions. In summary, CLADD scheme, in its two versions CCLADS and CCLADNS, has the following properties.

• It is characterized by cell-centered unknowns and a local stencil on two-dimensional unstructured grids; it reduces to a nine-point scheme on quadrangular distorted grids.

• For rectangular grids, in the case of isotropic diffusion, it reduces to the standard five-point scheme and the treatment of discontinuous conductivity coefficients is equivalent to the well known harmonic averaging procedure.

In addition, since CCLADS construction is based on a local variational formulation, it provides a positive definite (resp. symmetric) representation of the discrete diffusion operator if the conductivity tensor is positive definite (resp. symmetric). In its semi-discrete version, this scheme satisfies a L 2 stability property. For triangular grids, it preserves linear solutions and it is characterized by a second-order accuracy. For smooth distorted quadrangular grids it exhibits an accuracy which is almost second-order. On the other hand, CCLADNS does not lead in general to a symmetric positive definite discrete diffusion operator even if the conductivity tensor is symmetric positive definite. Its semi-discrete version does not satisfy a L 2 stability property. However, since the numerical flux approximation is based on a Taylor expansion, it is characterized by a second-order accuracy on general unstructured grids. It is interesting to note that CCLADS scheme, like the MPFA reference space method [START_REF] Aavvatsmark | Convergence of a symmetric MPFA method on quadrilateral grids[END_REF] suffers from a reduction or the loss of convergence on quadrangular random grids i.e., the grids with perturbations of order h, where h is the mesh size parameter. However, such grids are seldom encountered in real life simulations wherein the numerical method solving the diffusion equation is coupled to Lagrangian hydrodynamics. Moreover, let us point out that in the context of Arbitrary Lagrangian Eulerian (ALE) computations, the rezoning procedure inherent to the ALE algorithm produces smooth grid for which our scheme exhibits an almost second-order convergence.

The remainder of this paper is organized as follows. In Section 2 we first give the problem statement introducing the governing equations, the notation and assumptions and our motivation regarding the underlying physical model, which is strongly linked to plasma physics. This is followed by Section 3 which is devoted to the space discretization of CCLADS scheme. In this section, we derive the half-edge fluxes approximation by means of a sub-cell-based variational formulation. We also describe the elimination of the half-edge temperatures in terms of the cell-centered unknowns to achieve the construction of the global discrete diffusion operator. After stating the main properties of the semi-discrete scheme, we conclude this section by giving indications related to boundary conditions implementation. In Section 4 we present the space discretization of CCLADNS scheme wherein the half-edge fluxes approximation results from a finite difference approximation. We pursue in Section 5 by describing the extension of CCLAD scheme to cylindrical geometry. The time discretization is developed in Section 6. Finally, the robustness and the accuracy of CCLADS and CCLADNS schemes are assessed using various representative test cases in Section 7.

Problem statement

Governing equations

Let D be an open set of the d-dimensional space IR d , occupied by a thermally and possibly electrically conductive material. Let x denotes the vector position of an arbitrary point inside the domain D and t > 0 the time. The unsteady thermal state of the domain under consideration is described by means of the specific internal energy field ε = ε(x, t). In the absence of fluid motion, heat propagation throughout the domain is governed by a partial differential equation, which is known as the heat conduction equation and writes as

ρ ∂ε ∂t + ∇ • q = ρr. (1)
Here, ρ is a positive real valued function, which stands for the mass density of the material. The source term, r, corresponds to the specific heat supplied to the material. Specific internal energy is expressed in terms of the mass density and the temperature, T , through the use of an equation of state written under the form ε = ε(ρ, T ). Introducing the specific heat capacity at constant volume C v = ( ∂ε ∂T ) ρ leads to rewrite the heat conduction equation as

ρC v ∂T ∂t + ∇ • q = ρr. (2)
In writing this equation, we have assumed that either ε does not depend on density or ρ does not depend on time. Let us note that C v is a positive real valued function, it is also an intrinsic material properties which may depend on temperature. In the sequel, we assume that ρ, C v , and r are known functions. The vector valued function, q, stands for the heat flux. It is expressed in terms of the temperature gradient by means of the generalized Fourier law q = -K∇T.

The second-order tensor K is the conductivity tensor which is also an intrinsic property of the material under consideration. Let us point out that this representation for the flux based on the linear relation between the flux and the temperature gradient fails if the particle mean free path becomes comparable to the temperature scale length. This occurs very often in laser driven plasma flows. There more complicated models are employed related either to local flux limitation (flux limiters) or to the non-local relation between the flux and the temperature gradient, refer to [START_REF] Atzeni | The physics of inertial fusion[END_REF].

In the sequel, we will study uniquely classical heat transfer phenomena wherein the heat flux is described using the above Fourier law.

According to the Second Law of thermodynamics, Fourier law has to obey the constraint

q • ∇T ≤ 0. ( 3
)
This requirement on the constitutive law of the heat flux simply states that heat flux direction is opposite to temperature gradient, that is, heat flows from hot region to cold region. Mathematically speaking, this thermodynamic constraint amounts to state that the conductivity tensor is a positive definite tensor

Kφ • φ > 0, ∀φ ∈ IR d . (4)
Depending on the material under consideration, the conductivity tensor may be either symmetric, i.e., K t = K or may satisfy K(-B) = K t (B) if the material is in the magnetic field B, [START_REF] De Groot | Non-equilibrium thermodynamics[END_REF]. For an isotropic material, heat flux does not depend on the space directions and thus conductivity tensor reduces to

K = κ I d , (5)
where κ denotes the heat conductivity, which is a positive scalar valued function and I d denotes the unit tensor.

Being given the material properties, ρ, C, K and the heat supply r, we want to solve the heat conduction equation (2) on the computational domain D. To do so, we need to define initial and boundary conditions. The initial condition is prescribed as

T (x, 0) = T 0 (x), ∀x ∈ IR d , (6)
where T 0 denotes the initial temperature field. Regarding the boundary conditions, three types can be imposed on the boundary of the domain, ∂D. These are: Dirichlet, Neumann and Robin boundary conditions, they consists in specifying respectively the temperature, the flux and a weighted combination of the temperature and the flux. Introduction of the partition ∂D = ∂D D ∪∂D N ∪∂D R of the boundary domain, allows to write the boundary conditions as

T (x, t) = T ⋆ (x, t), ∀x ∈ ∂D D , (7a)
q(x, t) • n = q ⋆ N (x, t), ∀x ∈ ∂D N , (7b) αT (x, t) + βq(x, t) • n = q ⋆ R (x, t), ∀x ∈ ∂D R . (7c)
Here, T ⋆ , q ⋆ N and q ⋆ R denote respectively the prescribed temperature and fluxes. In addition, α = α(x, t) and β = β(x, t) are real valued specified functions.

Comment 1. Let us remark that the normal component of the heat flux at the interface between two media 1 and 2 with different properties is continuous,

(K∇T ) 1 • n 12 = (K∇T ) 2 • n 12 ,
where n 12 is the unit normal to the interface. The temperature itself is also continuous.

Underlying physical models

In deriving a numerical method to solve the heat conduction equation ( 2), we aim at developing numerical modeling of physical phenomena encountered in plasma physics. More precisely, we are concerned by heat transfer within laser-heated plasma flows such as those obtained in the domain of direct drive Inertial Confinement Fusion, refer to [START_REF] Atzeni | The physics of inertial fusion[END_REF]. In this context, the energy released by the laser is transferred throughout the plasma flows by means of electron heat conduction. Omitting the pressure work term, which results from coupling to hydrodynamics, the electron temperature, T e , is governed by a heat conduction equation similar to [START_REF] Aavatsmark | An introduction to multipoint flux approximations for quadrilateral grids[END_REF]. In the classical regime, the electron heat flux, q e is given by the Spitzer-Härm law: q e = -κ e ∇T e , where the electron thermal conductivity, κ e , depends on the electron temperature as a power law, i.e., κ e (T e ) ∼ T e , refer to [START_REF] Dautray | La fusion thermonucléaire par laser[END_REF][START_REF] Ya | Physics of Shock Waves and High-Temperature Hydrodynamic Phenomena[END_REF]. This corresponds to an isotropic nonlinear heat conduction equation. However, in presence of magnetic fields, this isotropic model for heat conduction is not valid anymore. Such a situation occurs frequently for laser driven plasma wherein the density and pressure gradients are not colinear. In this particular case, a self-generated magnetic field, B, is created by the rotational component of the ambipolar electric field, E = -1

eNe ∇P e , where e is the electron charge, N e is the electron density per unit volume and P e is the electron pressure, refer to [START_REF] Dautray | La fusion thermonucléaire par laser[END_REF]. Knowing that the magnetic field is governed by the Faraday law, i.e., ∂B ∂t + ∇ × E = 0, we deduce that the time evolution of the magnetic field is governed by the following equation

∂B ∂t = 1 e ∇( 1 N e ) × ∇P e . (8)
The magnetic field dramatically modifies electron heat transport leading to a anisotropic electron conductivity. It implies not only a reduction of the magnitude of the heat flux but also its rotation.

Using plasma kinetic theory, Braginskii [START_REF] Braginskii | Transport Process in a Plasma[END_REF] has obtained the following expression of the electron heat flux with magnetic field

q e = -κ (∇T e • b)b -κ ⊥ [∇T e -(∇T e • b)b] -κ ∧ b × ∇T e , (9)
where b = B |B| is unit vector corresponding to the direction of the magnetic field and κ , κ ⊥ and κ ∧ are scalar conductivities given in [START_REF] Braginskii | Transport Process in a Plasma[END_REF]. Let us consider a two-dimensional plasma flow in planar geometry. Let (e x , e y , e z ) be the orthonormal basis of IR 3 and suppose that the two-dimensional flow is contained in the frame (x, y) equipped with the orthonormal basis (e x , e y ). By virtue of (8), it is obvious that the self-generated magnetic field is transverse to the two-dimensional flow, that is, B = Be z . Setting b = B |B| , we have b = be z , where b 2 = 1. Bearing this in mind, the electron heat flux expression (9) collapses to Developing the above equation over the Cartesian frame (x, y) allows to write the electronic heat flux

q e = -κ ⊥ ∇T e -κ ∧ b × ∇T e . ( 10
)
q e = -K e ∇T e , (11)
where the second-order tensor K e corresponds to the electronic conductivity defined by

K e = κ ⊥ -bκ ∧ bκ ∧ κ ⊥ . ( 12
)
The Braginskii transport coefficients κ , κ ⊥ and κ ∧ can be expressed in terms of the Spitzer-Härm conductivity, κ e , as κ = κ e , κ ⊥ = κ e f ⊥ (Ω e τ ei ), κ ∧ = κ e f ∧ (Ω e τ ei ), where f ⊥ , f ∧ are the functions describing the magnetization of the heat flux. In addition, Ω e ∼| B |, is the electron cyclotron frequency and τ ei the electron-ion collision frequency. Note that Ω e has the dimension of the reciprocal of time, thus parameter Ω e τ ei is dimensionless; it describes the effect of the magnetic field on the electron heat conductivity as a ratio between the electron gyration time in the magnetic field and the electron collision time. Bearing this in mind we have displayed in Figure 1 the normalized Braginskii transport coefficients with respect to the normalized parameter Ω e τ ei knowing that

f ⊥ (x) = 1 1+x 2 , f ∧ (x) = x 1+x 2 .
For a weak magnetic field, that is, Ω e τ ei ∈ [0, 1], we have κ ⊥ > κ ∧ , whereas for a strong magnetic filed κ ⊥ < κ ∧ . In the limit | B |→ 0, the normalized parameter Ω e τ ei also tends to zero and the Braginskii coefficients behave as follows: κ ⊥ → κ e and κ ∧ → 0. In this regime, the anisotropic conductivity tensor K e recovers the isotropic Spitzer-Härm conductivity, i.e., K e → κ e I 2 .

We conclude this section by remarking that K e is not symmetric and transforms as K e (-b) = K t e (b). This non-symmetry of the conductivity tensor is a consequence of the presence of the magnetic field, this behavior is known as the Righi-Leduc effect, refer to [START_REF] De Groot | Non-equilibrium thermodynamics[END_REF] chapter XI. It rotates the heat flux vector without changing its absolute value. This property corresponds to the following result: for all arbitrary vector φ

K e φ • φ = κ ⊥ | φ | 2 .
Since κ ⊥ is positive, the above result shows that K e is a positive definite tensor which satisfies the thermodynamic requirement [START_REF] Aavatsmark | Discretization on unstructured grids for inhomogeneous, anistropic media. Part II: discussion and numerical results[END_REF]. It is interesting to note that for an arbitrary vector φ, K e φ can be decomposed as follows

K e φ = κ ⊥ φ + bκ ∧ R π 2 φ, (13)
where R π 2 denotes the counterclockwise rotation through the angle π 2 . The above equation states that the anisotropic conductivity tensor acts as an isotropic conductivity tensor supplemented by a rotation tensor which follows directly from the magnetic field. Finally, computing the divergence of the anisotropic heat flux [START_REF] Breil | Modelling of the magnetic field effects in hydrodynamic codes using a second order tensorial diffusion scheme[END_REF] yields

∇ • q e = -∇ • (κ ⊥ ∇T e ) + A • ∇T e , (14)
where

A = [-∂ ∂y (bκ ∧ ), ∂ ∂x (bκ ∧ )] t .
Under this form the anisotropic diffusion operator appears as the sum of an isotropic diffusion operator plus an advection operator characterized by the velocitylike vector A. This decomposition suggests to solve the anisotropic heat conduction equation discretizing separately the isotropic diffusion operator and the advection operator. However, such a splitting strategy may suffer from a lack of robustness in case of strong magnetic fields, refer to [48]. That is why, we prefer to develop a computational method devoted to the discretization of the whole anisotropic diffusion operator. This computational method, which will be presented in the sequel, is the natural extension to anisotropic heat conduction of the finite volume scheme that has been initially derived in [START_REF] Breil | A cell-centered diffusion scheme on two-dimensional unstructured meshes[END_REF].

Notation and assumptions

Our motivation is to describe a finite volume scheme that solves the anisotropic heat conduction equation on two-dimensional unstructured grids. Before we proceed any further, let us introduce the notation and the assumptions required for the present work. Let D be an open set of the two-dimensional space IR 2 . We aim at constructing a numerical scheme to solve the following initial-boundary-value problem

ρC v ∂T ∂t + ∇ • q = ρr, (x, t) ∈ D × [0, T], (15a)
T (x, t) = T 0 (x), x ∈ D, (15b)
T (x, t) = T ⋆ (x, t), x ∈ ∂D D , (15c)
q(x, t)

• n = q ⋆ N (x, t), x ∈ ∂D N , (15d)
where T > 0 denotes the final time. Equation (15a) is a partial differential parabolic equation of second order for the temperature T , wherein the conductive flux, q, is defined according to

q = -K∇T. ( 16
)
We suppose that the second-order tensor K is positive definite to ensure thermodynamic consistency.

In addition, we make the assumption that there exists h > 0 such that This condition states that the diffusion operator in (15a) is strongly elliptic [START_REF] Pironneau | Optimal Shape Design for Elliptic System[END_REF], this ensures the well-posedness of problem [START_REF] Edwards | Finite volume discretization with imposed flux continuity for the general tensor pressure equation[END_REF] by means of the Lax-Milgram theorem [START_REF] Raviart | Introduction à l'analyse numérique des équations aux dérivées partielles[END_REF].

Kφ • φ ≥ h | φ | 2 , ∀φ ∈ D. ( 17
) p p + p - n + pc n - pc ω pc ω c
Having defined the problem we want to solve, let us introduce some notation necessary to develop the discretization scheme. Let ∪ c ω c denotes a partition of the computational domain D into polygonal cells ω c . The counterclockwise ordered list of vertices (points) of cell c is denoted by P(c). In addition, p being a generic point, we define its position vector denoted as x p and the set C(p) which contains all the cells surrounding point p. Being given p ∈ P(c), p -and p + are the previous and next points with respect to p in the ordered list of vertices of cell c. Let ω c be a generic polygonal cell, for each vertex p ∈ P(c), we define the sub-cell ω pc by connecting the centroid of ω c to the midpoints of edges [p -, p] and [p, p + ] impinging at node p, refer to Figure 2. In two dimensions the sub-cell, as just defined, is always a quadrilateral regardless of the type of cells that compose the underlying grid. The boundaries of the cell ω c and the sub-cell ω pc are denoted respectively ∂ω c and ∂ω pc . Finally, considering the intersection between the cell and subcell boundaries, we introduce half-edge geometric data. As the name implies, a half-edge is a half of an edge and is constructed by splitting an edge down its length. More precisely, we define the two half-edges related to point p and cell c as ∂ω - pc = ∂ω pc ∩ [p -, p] and ∂ω + pc = ∂ω pc ∩ [p, p + ]. The unit outward normal and the length related to half-edge ∂ω ± pc are denoted respectively n ± pc and l ± pc . To proceed with the construction of numerical scheme, let us integrate (15a) over ω c and make use of the divergence formula. This leads to the weak form of the heat conduction equation

d dt ωc ρC v T (x, t) dv + ∂ωc q • n ds = ωc ρr(x, t)dv, (18)
where n denotes the unit outward normal to ∂ω c . We shall first discretize this equation in the spatial variable x. The physical data, ρ, C v and r are supposed to be known functions with respect to space and time variables. We represent them using a piecewise constant approximation over each cell ω c . The piecewise constant approximation of any variable will be denoted using subscript c. The tensor conductivity K space approximation is also constructed using a piecewise constant representation over each cell, which is denoted by K c . Concerning the unknown temperature field, the discretization method we are going to use is the finite volume method for which the finite dimensional space to which the approximate solution belongs is also the space of piecewise constant functions. Bearing this in mind, [START_REF] Gyrya | High-order mimetic finite difference method for diffusion problems on polygonal meshes[END_REF] rewrites

m c C vc d dt T c + ∂ωc q • n ds = m c r c , (19)
Here, m c denotes the mass of the cell, that is, m c = ρ c | ω c | where | ω c | stands for the volume of the cell. Let us point out that T c = T c (t) is nothing but the mean value of the temperature over ω c

T c (t) = 1 | ω c | ωc T (x, t) dv.
To define completely the space discretization it remains to discretize the surface integral in the above equation. To do so, let us introduce the following piecewise constant approximation of the normal heat flux over each half-edge

q ± pc = 1 l ± pc ∂ω ± pc q • n ds. ( 20
)
The scalar q ± pc stands for the half-edge normal flux related to the half-edge ∂ω ± pc . Knowing that ∂ω c = ∪ p∈P(c) ∂ω ± pc , the discretized heat conduction equation writes as

m c C vc d dt T c + p∈P(c) l - pc q - pc + l + pc q + pc = m c r c . (21)
We conclude this paragraph by introducing as auxiliary unknowns the half-edge temperatures T ± pc defined by

T ± pc = 1 l ± pc ∂ω ± pc T (x, t) ds. (22)
In writing this equation, we have also assumed a piecewise constant approximation of the temperature field over each half-edge. By virtue of Comment 1, the piecewise constant approximations of the normal heat flux and temperature along each edge are defined such that these half-edgebased quantities are continuous across each edge. To exhibit these continuity conditions, let us consider two neighboring cells, denoted by subscripts c and d, which share a given edge, refer to Figure 3. This edge corresponds to the segment [p, p + ], where p and p + are two consecutive points in the counterclockwise numbering attached to cell c. It also corresponds to the segment [r -, r], where r -and r are two consecutive points in the counterclockwise numbering attached to cell d. Obviously, these four labels define the same edge and thus their corresponding points coincide, i.e., p ≡ r, p + ≡ r -. The sub-cell of cell c attached to point p ≡ r is denoted ω pc , whereas the sub-cell of cell d attached to point r ≡ p is denoted ω rd . This double notation, in spite of its heaviness, allows to define precisely the half-edge fluxes and temperatures at the half-edge corresponding to the intersection of the two previous sub-cells. Namely, viewed from sub-cell ω pc (resp. ω rd ), the half-edge flux and temperature are denoted q + pc and T + pc (resp. q - rd and T - rd ). Bearing this notation The half-edge fluxes, q ± pc , q ± rd and temperatures, T ± cp , T ± rd are displayed using blue color.

T + pc = T - rd ω pc p - n - pc r + p + = r - q + pc + q - rd = 0 ω rd p = r n + rd n + pc n - rd T + rd q - pc T - pc q + rd T c ω c T d ω d
in mind, continuity conditions at the half-edge (ω pc ∪ ω rd ) for the half-edge fluxes and temperatures write explicitly as

q + pc + q - rd = 0, (23a)
T + pc = T - rd . (23b)
The continuity condition for the heat flux follows from the definition of the unit outward normals related to (ω pc ∪ ω rd ), i.e., n + pc = -n - rd . To achieve the space discretization of [START_REF] Hyman | Mimetic finite difference methods for diffusion equations[END_REF], it remains to construct a consistent approximation of the half-edge normal flux, that is, to define a numeric half-edge-based flux function h ± pc such that

q - pc = h - pc (T - pc -T c , T + pc -T c ), q + pc = h + pc (T - pc -T c , T + pc -T c ). (24)
Here, h ± pc denotes a real valued function which is continuous with respect to its arguments. Let us note that we have expressed this function in terms of the temperature difference T c -T ± pc since the heat flux is proportional to the temperature gradient. The next steps in the design of our finite volume scheme will be the following:

• Definition of an approximation of the half-edge numerical fluxes in terms of the half-edge temperatures and the cell-centered temperature.

• Elimination of the half-edge temperatures through the use of the continuity condition (23) across sub-cell interface.

These tasks will be the main topics of the next sections.

Space discretization of the CCLADS scheme

In this section, we present the space discretization of the CCLADS scheme, wherein the half-edge fluxes approximation results from a local variational formulation. Before proceeding any further, we start by giving a useful and classical result concerning the representation of a vector in terms of its normal components. This result leads to the expression of the standart inner product of two vectors, which will be one the tools utilized to derive the sub-cell variational formulation.

Expression of a vector in terms of its normal components

Here, we recall briefly the methodology which has been thoroughly exposed by Shashkov in [START_REF] Shashkov | Conservative Finite-Difference Methods on General Grids[END_REF][START_REF] Morel | A local support-operators diffusion discretization scheme for quadrilateral r-z meshes[END_REF]. Let φ be an arbitrary vector of the two-dimensional space IR 2 and φ pc its piecewise constant approximation over the sub-cell ω pc . Let φ ± pc be the half-edge normal components of φ pc , that is,

φ pc • n - pc = φ - pc , φ pc • n + pc = φ + pc .
Introducing the corner matrix

J pc = [n - pc , n + pc ]
, the above 2 × 2 linear system rewrites

J t pc φ = φ - pc φ + pc ,
where superscript t denotes the transpose matrix. Provided that n - pc and n + pc are not colinear, the above system has always a unique solution written under the form

φ pc = J -t pc φ - pc φ + pc . (25)
This equation allows to express any vector in terms of its normal components on two non-colinear unit vectors. This representation allows to compute the inner product of two vectors φ pc and ψ pc as follows

φ pc • ψ pc = (J t pc J pc ) -1 ψ - pc ψ + pc • φ - pc φ + pc . ( 26
)
The 2 × 2 matrix H pc = J t pc J pc is defined by

H pc = n - pc • n - pc n - pc • n + pc n + pc • n - pc n + pc • n + pc = 1 -cos θ pc -cos θ pc 1 , (27)
where θ pc denotes the measure of the angle between the two half-edges of sub-cell ω pc impinging at point p, refer to Figure 4. This matrix admits an inverse provided that θ pc = kπ, where k is an integer. Under this condition, H -1 pc is readily obtained

H -1 pc = 1 sin 2 θ pc 1 cos θ pc cos θ pc 1 .
This matrix, which is symmetric definite positive, represents the local metric tensor associated to the sub-cell ω pc . Let us remark that we have recovered exactly the expressions initially derived in [START_REF] Morel | A local support-operators diffusion discretization scheme for quadrilateral r-z meshes[END_REF].

Half-edge fluxes approximation based on a local variational formulation 3.2.1. Sub-cell-based variational formulation

We construct an approximation of the half-edge fluxes by means of a local variational formulation written over the sub-cell ω pc . Contrary to the classical cell-based variational formulation used in the context of Mimetic Finite Difference Method [START_REF] Hyman | Mimetic finite difference methods for diffusion equations[END_REF], the present sub-cell-based variational formulation leads to a local explicit expression of the half-edges fluxes in terms of the half-edges temperatures and the mean cell temperature. The local and explicit feature of the half-edge fluxes expression is of great importance, since it allows to construct a numerical scheme with only one unknown per cell.

Our starting point to derive the sub-cell-based variational formulation consists in writing the partial differential equation satisfied by the flux. From the heat flux definition [START_REF] Eymard | Finite Volume methods[END_REF], it follows that q satisfies K -1 q + ∇T = 0.

Let us point out that the present approach is stronlgy linked to the mixed formulation utilized in the context of mixed finite element discretization [START_REF] Thomas | Mixed finite volume methods[END_REF][START_REF] Aavvatsmark | Convergence of a symmetric MPFA method on quadrilateral grids[END_REF][START_REF] Lipnikov | Local flux mimetic finite difference methods[END_REF]. Dot-multiplying this equation by an arbitrary vector φ ∈ D and integrating over the cell ω pc yields

ωpc φ • K -1 q dv = - ωpc φ • ∇T dv, ∀φ ∈ D. ( 29
)
Integrating by part the right-hand side and applying the divergence formula leads to the following variational formulation

ωpc φ • K -1 q dv = ωpc T ∇ • φ dv - ∂ωpc T φ • n ds, ∀φ ∈ D. (30)
This sub-cell-based variational formulation is the base to construct a local and explicit numerical approximation of the half-edge fluxes. Replacing T by its piecewise constant approximation T c in the first integral of the right-hand side and applying the divergence formula to the remaining volume integrals leads to

ωpc φ • K -1 q dv = T c ∂ωpc φ • n ds - ∂ωpc T φ • n ds.
Partitioning the sub-cell boundary as ∂ω pc = (∂ω pc ∩ ∂ω c ) ∪ (∂ω pc ∩ ω c ) in the latter equation yields

ωpc φ•K -1 q dv = T c ∂ωpc∩∂ωc φ•n ds+T c ∂ωpc∩ωc φ•n ds- ∂ωpc∩∂ωc T φ•n ds- ∂ωpc∩ωc T φ•n ds.
Replacing T by T c in the fourth surface integral of the right-hand side and noticing that the second integral is equal to the last one allows to write the sub-cell-based variational formulation under the form

ωpc φ • K -1 q dv = T c ∂ωpc∩∂ωc φ • n ds - ∂ωpc∩∂ωc T φ • n ds. (31)
At this point it is interesting to remark that this sub-cell-based formulation is a sufficient condition to recover the classical cell-based variational formulation. Since the set of sub-cells of ω c is a partition of this cell, we have The half-edge fluxes, q ± pc , and temperatures, T ± pc are displayed using blue color.

ω c = p∈P(c) ω pc , ∂ω c = p∈P(c) (∂ω pc ∩ ∂ω c ). q + pc T + pc p T c ω pc n - pc n + pc p + T - pc q - pc ω c θ pc p -
Thus, by summing ( 31) over all the sub-cells of ω c , we obtain

ωc φ • K -1 q dv = T c ∂ωc φ • n ds - ∂ωc T φ • n ds. (32)
This last equation corresponds to cell-based variational formulation of the partial differential equation [START_REF] Lipnikov | Monotone finite volume schemes for diffusion equations on unstructured triangular and shape-regular polygonal meshes[END_REF]. This form is used in the context of Mimetic Finite Difference Method [START_REF] Hyman | Mimetic finite difference methods for diffusion equations[END_REF] Returning to the sub-cell based variational formulation, we discretize the right-hand side of (31) by introducing the half-edge normal components of φ and the piecewise constant approximation of the half-edge temperatures as follows

ωpc φ • K -1 q dv = -[l - pc (T - pc -T c )φ - pc + l + pc (T + pc -T c )φ + pc ]. (33)
Assuming a piecewise constant representation of the test function allows to compute the volume integral in the left-hand side thanks to the quadrature rule

ωpc φ • K -1 q dv = w pc φ pc • K -1 c q pc , (34)
Here, K c denotes the piecewise constant approximation of the conductivity tensor and φ pc , q pc are the piecewise constant approximations of vectors φ and q, refer to Figure 4. In addition, w pc denotes some positive corner volume related to sub-cell ω pc , which will be determined later. Note that the corner volumes associated to the same cell ω c must satisfy the consistency condition p∈P(c)

w pc =| ω c | . (35)
Namely, the corner volumes of a cell sums to the volume of the cell. This is the minimal requirement to ensure that constant functions are exactly integrated using the above quadrature rule. Now, combining ( 34) and ( 33) and using the expression of the vectors q and φ in terms of their half-edge normal components leads to the following variational formulation

w pc (J t pc K c J pc ) -1 q - pc q + pc • φ - pc φ + pc = - l - pc (T - pc -T c ) l + pc (T + pc -T c ) • φ - pc φ + pc . ( 36
)
Knowing that this variational formulation must hold for any vector φ pc , this implies

q - pc q + pc = - 1 w pc (J t pc K c J pc ) l - pc (T - pc -T c ) l + pc (T + pc -T c ) . ( 37
)
This equation constitutes the approximation of the half-edge normal fluxes over a sub-cell. This local approximation is coherent with expression of the constitutive law ( 16) in the sense that the numerical approximation of the heat flux is equal to a tensor times a numerical approximation of the temperature gradient. This tensor can be viewed as an effective conductivity tensor associated to the sub-cell ω pc . Thus, it is natural to set

K pc = J t pc K c J pc . (38)
Let us emphasize that this corner tensor inherits all the properties of the conductivity tensor K c . Namely, K c being positive definite, K pc is also positive definite. This comes from the fact that

K pc φ • φ = K c (J pc φ) • (J pc φ), ∀φ ∈ IR 2 .
Using a similar argument, note that if K c is symmetric, K pc is also symmetric. Recalling that

J pc = [n - pc , n + pc ],
we readily obtain the expression of the corner tensor K pc in terms of the unit normal n ± pc

K pc = K c n - pc • n - pc K c n + pc • n - pc K c n - pc • n + pc K c n + pc • n + pc . (39)
Let us remark that in the isotropic case, i.e, K c = κ c I d , the corner tensor collapses to

K pc = κ c H pc , (40)
where κ c denotes the piecewise constant scalar conductivity over cell ω c and H pc is the second-order tensor defined by [START_REF] Lipnikov | The mimetic finite difference discretization of diffusion problem on unstructured polyhedral meshes[END_REF]. We conclude by claiming that a sub-cell-based variational formulation has allowed to construct the following numerical approximation of the half-edge normal fluxes

q - pc q + pc = - 1 w pc K pc l - pc (T - pc -T c ) l + pc (T + pc -T c ) . (41)
Here, w pc is a positive volume weight, which will be determined later, and the corner conductivity tensor, K pc is expressed by [START_REF] Pironneau | Optimal Shape Design for Elliptic System[END_REF].

Comment 2. It is interesting to remark that the corner tensor K pc is a linear function with respect to the piecewise constant approximation of the conductivity tensor K c . This follows directly from [START_REF] Moulton | Workshop: Compatible Spatial Discretizations for Partial Differential Equations[END_REF]. In addition, the corner tensor corresponding to the transpose of K c is the transpose of K pc , i.e., K pc (K t c ) = K t pc (K c ).

Fundamental inequality satisfied by the discrete approximation of the half-edge fluxes

The goal of this paragraph is to show that the discrete approximation of the half-edges normal fluxes [START_REF] Potier | Schéma volumes finis monotone pour des opérateurs de diffusion fortement anisotropes sur des maillages de triangles non structurés[END_REF] derived from the sub-cell-based variational formulation satisfies a discrete version of the fundamental inequality (3), which follows from the Second Law of thermodynamics. This discrete analogous of the fundamental inequality states that for half-edge fluxes defined according to [START_REF] Potier | Schéma volumes finis monotone pour des opérateurs de diffusion fortement anisotropes sur des maillages de triangles non structurés[END_REF] the following inequality holds

c∈C(p) (l - pc q - pc + l + pc q + pc )T c ≥ 0, ( 42
)
where C(p) denotes the set of cells surrounding point p. To prove this inequality, let us introduce I p as being the nodal quantity defined by

I p = c∈C(p) (l - pc q - pc + l + pc q + pc )T c . (43)
We shall prove that I p is always non-negative using the sub-cell variational formulation derived in Section 3.2. From (36) it follows that for all vector φ the following identity holds

w pc K -1 pc q - pc q + pc • φ - pc φ + pc = - l - pc (T - pc -T c ) l + pc (T + pc -T c ) • φ - pc φ + pc ,
where w pc is the positive volume weight and K pc the definite positive corner conductivity tensor. Applying this identity for φ = q and rearranging the right hand-side yields

w pc K -1 pc q - pc q + pc • q - pc q + pc = (l - pc q - pc + l + pc q + pc )T c -(l - pc q - pc T - pc + l + pc q + pc T + pc ).
We notice that the left-hand side of the above equation is always non-negative since K pc is positive definite. Summing the above equation over all the cells surrounding point p leads to c∈P(c)

w pc K -1 pc q - pc q + pc • q - pc q + pc = c∈P(c) (l - pc q - pc + l + pc q + pc )T c - c∈P(c) (l - pc q - pc T - pc + l + pc q + pc T + pc ). ( 44
)
It is interesting to mention that the above equation is the discrete analogous of the following integral identity

- ωp ∇T • q dv = ωp T ∇ • q dv - ∂ωp T q • n ds, (45)
where ω p denotes the dual cell which results from the union of the sub-cells surrounding point p, i.e., ω p = c∈C(p) ω pc . Returning to [START_REF] Raviart | Introduction à l'analyse numérique des équations aux dérivées partielles[END_REF], we observe that the second term in the right-hand side vanishes due to the continuity condition of the fluxes at the half-edges impinging at point p, refer to (23a). Finally, (44) turns to Let us emphasize that inequality [START_REF] Potier | Schéma volumes finis pour des opérateurs de diffusion fortement anisotropes sur des maillages non structurés[END_REF] follows directly from the fact that we have used a variational formulation to derive the numerical approximation of the flux. This inequality will be of great importance to prove several crucial results regarding the properties of the finite volume scheme CCLADS.

I p = c∈C(p) w pc K -1 pc q - pc q + pc • q - pc q + pc ≥ 0. ( 46
) p + p θ pc T c T + pc n + pc t + pc p - n - pc t - pc T - pc

Volume weight computation

In this paragraph, we aim at deriving practical formulas to compute the volume weight, w pc , present in the flux approximation [START_REF] Potier | Schéma volumes finis monotone pour des opérateurs de diffusion fortement anisotropes sur des maillages de triangles non structurés[END_REF]. To begin with, let us consider a triangular cell, ω c , characterized by its counterclockwise ordered vertices p -, p and p + , refer to Figure 5. We state that the flux approximation (41) preserves linear fields over triangular cells provided that the volume weight is such that

w tri pc = 1 3 | ω c | . ( 47
)
To prove this result, let us consider T h = T h (x) a piecewise linear approximation of the temperature field, i.e.,

T h (x) = T c + (∇T ) c • (x -x c ), ∀x ∈ ω c . (48)
Here, x c = 1 3 (x p -+ x p + x p + ) is the centroid of ω c and T c = T h (x c ) denotes the mean temperature of the cell. In addition, (∇T ) c corresponds to the uniform temperature gradient of the cell. Using the piecewise constant approximation of the conductivity tensor, K c , this gradient is rewritten (∇T ) c = -K -1 c q c , where q c is the piecewise constant approximation of the flux. With this notation, (48) transforms into

T h (x) = T c -K -1 c q c • (x -x c ), ∀x ∈ ω c . (49)
Expressing the two vectors q c and (x-x c ) in terms of their half-edge normal components by means of (25) yields

T h (x) = T c -K -1 pc q - pc q + pc • (x -x c ) - pc (x -x c ) + pc , ∀x ∈ ω c , (50)
where K pc = J t pc K c J pc . Since this equation holds for all points in ω c , we apply it to x - pc and x + pc given by

x - pc = 2x p + x p - 3 , x + pc = 2x p + x p + 3 . (51)
This results in

T h (x ± pc ) -T c = -K -1 pc q - pc q + pc • (x ± pc -x c ) - pc (x ± pc -x c ) + pc .
Knowing that

x - pc -x c = 1 3 (x p -x p + ) = - 2 3 l + pc t + pc , x + pc -x c = 1 3 (x p -x p -) = 2 3 l - pc t - pc ,
where t ± pc are the half-edge unit tangent vectors such that n ± pc × t ± pc = e z , refer to Figure 5, using

t - pc • n + pc = sin θ pc and t + pc • n - pc = -sin θ pc leads to T h (x - pc ) -T c = - 2 3 l + pc sin θ pc K -1 pc q - pc q + pc • 1 0 , T h (x + pc ) -T c = - 2 3 l - pc sin θ pc K -1 pc q - pc q + pc • 0 1 .
Rearranging the above equations allows to express the half-edge normal components of the flux as

q - pc q + pc = - 3 2l - pc l + pc sin θ pc K pc l - pc (T h (x - pc ) -T c ) l + pc (T h (x + pc ) -T c ) . (52)
In writing this equation we have obtained an expression of the half-edge fluxes which is exact for a linear approximation of the temperature field over a triangular cell. The comparison between this formula and the general formula obtained previously shows that the volume weight is given by w pc = 2 3 l - pc l + pc sin θ pc , which is nothing but one third of the cell volume. In addition, this comparison reveals that the piecewise constant half-edge approximations of the temperature have a clear geometrical interpretation since T ± pc = T h (x ± pc ), refer to Figure 5. Having defined the volume weight for triangular cells, we conclude this paragraph by giving some indications about the volume weight definition for other types of cells. For quadrangular cells, according to [START_REF] Hyman | Mimetic finite difference methods for diffusion equations[END_REF], a reasonable choice is to set

w quad pc = l - pc l + pc sin θ pc . (53)
This results in a corner volume equal to the half of the area of the triangle formed by points p -, p and p + , refer to Fig 5 . Unfortunately this choice does not allow to preserve linear solution on quadrangular grids, except on grids made of parallelograms, refer to Comment 7. However, the numerical results obtained on quadrangular grids with this choice appeared to be quite satisfactory as we shall see in the section devoted to the numerical results. For general polygonal cells, two possible choices are obtained setting

w poly1 pc = 1 | P(c) | | ω c |, (54a)
w poly2 pc =| ω pc |, (54b)
where | P(c) | is the total number of sub-cells in cell c. Comment 3. Let us point out that all these volume weights must be positive and consequently they have been defined assuming that we are dealing with valid cells, that it, is convex cells. The occurence of non-convex cells can be treated using the rough remedy which consists in replacing the original volume weight by its absolute value.

c + 1 c c -1 T c ω p,c-1 ω c-1 T c-1 ω c n c c n c-1 c-1 ω pc n c-1 c T c T c+1 T c-1 n c c+1 p

Elimination of the half-edge temperatures

From ( 41), it appears that the numerical approximation of the half-edge fluxes at a corner depends on the difference between the mean cell temperature and the half-edges temperatures. The mean cell temperature is the primary unknown whereas the half-edge temperatures are auxiliary unknowns, which can be eliminated by means of continuity argument (23a). Namely, we use the fact that the half-edge normal fluxes are continuous across each half-edges impinging at a given point. This local elimination procedure, which will be describe below, yields a linear system satisfied by the half-edge temperatures. We will show that this system admits always a unique solution which allows to express the half-edge temperatures in terms of the mean temperatures of the cells surrounding the point under consideration. Therefore, this local elimination procedure results in a finite volume discrete scheme with one unknown per cell.

Local notation around a point

To derive the local elimination procedure, we shall introduce some convenient notation. Let p denotes a generic point which is not located on the boundary ∂D. The treatment of boundary points is postponed to Section 3.6, which is devoted to boundary conditions implementation. Let C(p) be the set of cells that surround point p. The edges impinging at point p are labelled using the subscript c ranging from 1 to C p , where C p denotes the total number of cells surrounding point p.

The cell (sub-cell) numbering follows the edge numbering, that is, cell ω c (sub-cell ω pc ) is located between edges c and c+1, refer to Figure 6. The unit outward normal to cell ω c at edge c is denoted by n c c whereas the unit outward normal to cell ω c at edge c + 1 is denoted by n c+1 c . Assuming the continuity of the half-edge temperatures leads to denote by T c the unique half-edge temperature of the half-edge c impinging at point p. Note that we have omitted the dependency on point p in the indexing each time this is possible to avoid too heavy notion. With this notation, the expression of the half-edge fluxes (41) turns into

q c c q c c+1 = - 1 w pc K pc l c (T c -T c ) l c+1 (T c+1 -T c ) , ∀c ∈ C(p). (55)
Here, q c c (resp. q c c+1 ) denotes the half-edge normal flux at edge c (resp. c + 1) viewed from cell c. In addition l c denotes the half of the length of edge c. In writing these equations, we assume a periodic numbering around the point p. According to [START_REF] Pironneau | Optimal Shape Design for Elliptic System[END_REF], the sub-cell conductivity tensor is defined as

K pc = K c n c c • n c c K c n c c+1 • n c c K c n c c • n c c+1 K c n c c+1 • n c c+1 , ∀c ∈ C(p), (56)
where K c is the piecewise constant approximation of the conductivity tensor in cell c. Combining ( 55) and ( 56) yields the explicit expressions

q c c = -α c [l c (K c n c c • n c c )(T c -T c ) + l c+1 (K c n c c+1 • n c c )(T c+1 -T c )], (57a)
q c c+1 = -α c [l c (K c n c c • n c c+1 )(T c -T c ) + l c+1 (K c n c c+1 • n c c+1 )(T c+1 -T c )], (57b)
where we have introduced the inverse of the volume weight setting α c = 1 wpc . Shifting index c, i.e., c → c -1, in (57b) leads to the following expression for the half-edge normal flux at edge c viewed from cell c -1

q c-1 c = -α c-1 [l c-1 (K c-1 n c-1 c-1 • n c-1 c )(T c-1 -T c-1 ) + l c (K c-1 n c-1 c • n c-1 c )(T c -T c-1 )]. ( 58
)
3.3.2. Linear system satisfied by the half-edge temperatures Bearing this in mind, we are now in position to proceed with the elimination of the half-edge temperatures by writing the continuity of the half-edge normal fluxes at each edge c. This continuity condition at edge c reads as

l c q c-1 c + l c q c c = 0, ∀c ∈ C(p). (59)
Let us remark that this continuity condition provides C p equations for the C p auxiliary unknowns T c . Substituting (58) and (57a) into the continuity condition yields

α c-1 l c-1 l c (K c-1 n c-1 c-1 • n c-1 c )T c-1 + [α c-1 l 2 c (K c-1 n c-1 c • n c-1 c ) + α c l 2 c (K c n c c • n c c )]T c + α c l c l c+1 (K c n c c+1 • n c c )T c+1 = α c-1 l c [l c-1 (K c-1 n c-1 c-1 • n c-1 c ) + l c (K c-1 n c-1 c • n c-1 c )]T c-1 + α c l c [l c (K c n c c • n c c ) + l c+1 (K c n c c+1 • n c c )]T c .
To write this equation under a more concise form, let us introduce T = (T 1 , . . . , T Cp ) t as the vector of the cell-centered temperatures around point p and T = (T 1 , . . . , T Cp ) t as the vector of the halfedge temperatures around point p. The continuity condition (59) amounts to write that T satisfies the following C p × C p linear system

MT = ST . (60)
Let us remark that M is a tridiagonal cyclic matrix. This cyclic form is natural consequence of the periodic numbering we have used in solving continuity equations (59). The non-zero terms corresponding to the cth row of this matrix write as

     M c,c-1 = α c-1 l c-1 l c (K c-1 n c-1 c-1 • n c-1 c ), M c,c = α c-1 l 2 c (K c-1 n c-1 c • n c-1 c ) + α c l 2 c (K c n c c • n c c ), M c,c+1 = α c l c l c+1 (K c n c c+1 • n c c ). ( 61
)
From the first equation it follows that

M c+1,c = α c l c l c+1 (K c n c c • n c c+1
). The comparison of this term with M c,c+1 shows that M is symmetric if and only if the conductivity tensor, K c is also symmetric. Regarding S, it is a bidiagonal cyclic matrix, the non-zero terms corresponding to the cth row are:

S c,c-1 = α c-1 l c [l c-1 (K c-1 n c-1 c-1 • n c-1 c ) + l c (K c-1 n c-1 c • n c-1 c )], S c,c = α c l c [l c (K c n c c • n c c ) + l c+1 (K c n c c+1 • n c c )]. (62)
We remark that M can be decomposed as

M = LNL. (63)
N c,c-1 = α c-1 (K c-1 n c-1 c-1 • n c-1 c ), N c,c = α c-1 (K c-1 n c-1 c • n c-1 c ) + α c (K c n c c • n c c ), N c,c+1 = α c (K c n c c+1 • n c c
). Let us write this matrix explicitly in the particular case C p = 4, which corresponds to a point surrounding by 4 cells

N =   α 4 (K 4 n 4 1 •n 4 1 )+α 1 (K 1 n 1 1 •n 1 1 ) α 1 (K 1 n 1 2 •n 1 1 ) 0 α 4 (K 4 n 4 4 •n 4 1 ) α 1 (K 1 n 1 1 •n 1 2 ) α 1 (K 1 n 1 2 •n 1 2 )+α 2 (K 2 n 2 2 •n 2 2 ) α 2 (K 2 n 2 3 •n 2 2 ) 0 0 α 2 (K 2 n 2 2 •n 2 3 ) α 2 (K 2 n 2 3 •n 2 3 )+α 3 (K 3 n 3 3 •n 3 3 ) α 3 (K 3 n 3 4 •n 3 3 ) α 4 (K 4 n 4 1 •n 4 4 ) 0 α 3 (K 3 n 3 3 •n 3 4 ) α 3 (K 3 n 3 4 •n 3 4 )+α 4 (K 4 n 4 4 •n 4 4 )   .
A closer inspection of the above matrix reveals an interesting block-structure. Namely, N can be decomposed as

N = 4 c=1 α c N c ,
where N c are 4 × 4 sparse matrices which read

N 1 =     K 1 n 1 1 • n 1 1 K 1 n 1 2 • n 1 1 0 0 K 1 n 1 1 • n 1 2 K 1 n 1 2 • n 1 2 0 0 0 0 0 0 0 0 0 0     , N 2 =     0 0 0 0 0 K 2 n 2 2 • n 2 2 K 2 n 2 3 • n 2 2 0 0 K 2 n 2 2 • n 2 3 K 2 n 2 3 • n 2 3 0 0 0 0 0     , N 3 =     0 0 0 0 0 0 0 0 0 0 K 3 n 3 3 • n 3 3 K 3 n 3 4 • n 3 3 0 0 K 3 n 3 3 • n 3 4 K 3 n 3 4 • n 3 4     , N 4 =     K 4 n 4 1 • n 4 1 0 0 K 4 n 4 4 • n 4 1 0 0 0 0 0 0 0 0 K 4 n 4 1 • n 4 4 0 0 K 4 n 4 4 • n 4 4     .
This decomposition extends readily to the general case as

N = Cp c=1 α c N c , (64)
where N c is the C p × C p defined by

N c =           0 • • • 0 0 • • • 0 . . . . . . . . . . . . . . . 0 • • • K c n c c • n c c K c n c c+1 • n c c . . . 0 0 • • • K c n c c • n c c+1 K c n c c+1 • n c c+1 . . . 0 . . . . . . . . . . . . . . . 0 • • • 0 0 • • • 0           . ( 65
)
The non-zero terms of this matrix consist of 2 × 2 block which is located at the intersection of the cth and c + 1th columns and rows. In addition, this 2 × 2 block coincides with the sub-cell conductivity matrix K pc defined by [START_REF] Ya | Physics of Shock Waves and High-Temperature Hydrodynamic Phenomena[END_REF].

Bearing this decomposition in mind, we claim that M is a positive definite matrix, that is,

MT • T > 0, ∀T ∈ IR Cp . (66)
To prove this fundamental result we proceed in two steps. First, we note that M is positive definite if and only if N is positive definite. This first result follows from

MT • T =(LNL)T • T , thanks to (63) =L[(NL)T ] • T =N(LT ) • (LT ),
since L is symmetric, i.e., L t = L. Second, we prove that N is positive definite by computing NT • T using the decomposition (64)

NT • T = Cp c=1 α c N c T • T = Cp c=1 α c K pc T c T c+1 • T c T c+1 .
Since α c is positive and K pc is positive definite, the right-hand side of the last equation is always positive, which ends the proof. This shows that matrix M inherits the properties of the corner conductivity matrix K pc which results from the discretization of the half-edge normal fluxes. Particularly, M is symmetric provided that K pc is also symmetric. This properties transfer from K pc to M also emphasizes that the first step in the numerical scheme derivation, i.e., the half-edge normal fluxes discretization, must be performed with great care. We state that the linear system (60) has always a unique solution which provides the half-edge temperatures in terms of the mean cell temperatures surrounding a given point p

T = (M -1 S)T , (67)
where the matrices M and S are defined respectively by ( 61) and (62). The solution of the corresponding tridiagonal cyclic system is easily obtained by using the numerical algorithm proposed in [START_REF] Press | Numerical Recipes in Fortran 77[END_REF].

Comment 4. It is important to note that the solution of the linear system (67) has to preserve uniform temperatures field. Thus, matrix M -1 S satisfies the following property

(M -1 S)1 Cp = 1 Cp , (68)
where 1 Cp ∈ IR Cp is vector whose all entries are equal to 1. This means that 1 Cp is the eigenvector of M -1 S associated to the eigenvalue 1.

It is interesting to mention that matrices M et S exhibit a linear dependency on the conductivity tensor K. This result is a direct consequence of the definitions (61) and (62). In addition, the matrix M associated to the transpose of K corresponds to the transpose of the matrix M associated to K, i.e.,

M(K t 1 , . . . , λK t Cp ) = M t (K 1 , . . . , K Cp ). (69)
By way of contrast, the matrix S associated to the transpose of K does not correspond the transpose of the matrix S associated to K since S is a bidiagonal matrix. In this case, we need to introduce the new matrix S defined as

S(K 1 , . . . , K Cp ) = S(K t 1 , . . . , λK t Cp ). ( 70
)
From this definition and (61), the non-zero entries of S write explicitly

Sc,c-1 = α c-1 l c [l c-1 (K c-1 n c-1 c • n c-1 c-1 ) + l c (K c-1 n c-1 c • n c-1 c )], Sc,c = α c l c [l c (K c n c c • n c c ) + l c+1 (K c n c c • n c c+1 )]. (71)
Here, we have used the fact that

K t c n • m = K c m • n for all vectors (n, m) ∈ IR 2 × IR 2 .
Note that S = S if and only if the conductivity tensor is symmetric.

Maximum principle for the half-edge temperatures in the case of a symmetric positive definite conductivity tensor

In this paragraph, we aim at deriving sufficient conditions such that the half-edges temperatures satisfy a maximum principle, in the case of a symmetric positive definite conductivity tensor. More precisely, suppose that the mean cell temperatures are such that 0 < T c ≤ Θ for all c ∈ C(p), we want to exhibit conditions related to matrices M and S so that the half-edge temperatures satisfy also 0 < T c ≤ Θ for all c ∈ C(p), where Θ > 0 is a given temperature. To derive these sufficient conditions, we need to introduce the notion of M-matrix. Following the definition of [START_REF] Varga | Matrix Iterative Analysis[END_REF], a real matrix A with non-positive off-diagonal entries is an M-matrix if A is non-singular and the entries of A -1 are non-negative. A useful characterization of M-matrix is given by the following proposition: if A is a symmetric positive definite matrix with non-positive off-diagonal entries, then it is an Mmatrix. This latter kind of matrix is sometimes called Stieltjes matrix [START_REF] Varga | Matrix Iterative Analysis[END_REF]. Bearing this in mind, we are going to find a sufficient condition to ensure that M is an M-matrix. First, let us point out that M is symmetric definite positive since by assumption the conductivity tensor, K c , is symmetric definite positive. Thus, observing the off-diagonal entries of M given by (61), a sufficient condition to ensure that M is an M-matrix consists in prescribing

(K c n c c+1 • n c c ) ≤ 0, ∀c ∈ C(p). (72)
Since K c is symmetric positive definite, Schwartz inequality leads to

| (K c n c c+1 • n c c ) | (K c n c c+1 • n c c+1 ) 1 2 (K c n c c • n c c ) 1 2 ≤ 1.
Hence, there exists a unique

ν c ∈ [0, π] such that cos ν c = - (K c n c c+1 • n c c ) (K c n c c+1 • n c c+1 ) 1 2 (K c n c c • n c c ) 1 2
, ∀c ∈ C(p).

Defining the symmetric positive definite tensor U c as being the square root of K c , i.e., K c = U 2 c , the angle ν c rewrites as

cos ν c = - (U c n c c+1 • U c n c c ) (U c n c c+1 • U c n c c+1 ) 1 2 (U c n c c • U c n c c ) 1 2
, ∀c ∈ C(p).

This shows that πν c is the measure of the angle between the vectors U c n c c+1 and U c n c c , which are the images of the unit normal vectors n c c+1 and n c c in the transformation associated to the local metric defined by tensor K c . With this notation, sufficient condition (72), turns into the more explicit form

ν c ∈ [- π 2 , π 2 ], ∀c ∈ C(p). (74)
Assuming this condition, we have M -1 c,d ≥ 0 for all (c, d) ∈ C 2 (p). Recalling that T = (M -1 S)T , to ensure the maximum principle for T it remains to exhibit a sufficient condition such that the entries of S are non-negative. Introducing the angle ν c in the expressions of S entries (62), this condition readily writes as

cos ν c ≤ min( l c l c+1 (K c n c c • n c c ) 1 2 (K c n c c+1 • n c c+1 ) 1 2 , l c+1 l c (K c n c c+1 • n c c+1 ) 1 2 (K c n c c • n c c ) 1 2
), ∀c ∈ C(p).

This condition results in a limitation of the permitted values of ν c , this limitation expressing in terms of the cell aspect ratio with respect to the local metric defined by K c . Finally, we claim that the maximum principle for the half-edge temperatures holds provided the sufficient conditions (74) and (75) are satisfied. Indeed, assuming (74) and (75) implies that the entries of M -1 S are non-negative. In addition, by virtue of Comment 68 we have

d∈C(p) (M -1 S) c,d = 1, ∀c ∈ C(p).
Thus, the entries of M -1 S satisfy the following inequality

0 ≤ (M -1 S) c,d ≤ 1, ∀(c, d) ∈ C(p) 2 .
Therefore, each component of T = M -1 ST is a convex combination of the components of T since

T c = c∈C(p) (M -1 S) c,d T d .
Using the above arguments it is clear that if T d ∈ [0, Θ] then T c ∈ [0, Θ], which ends the proof.

We conclude by stating the following maximum principle for the half-edge temperatures: being given a symmetric positive definite conductivity tensor and assuming that the following geometric conditions hold for all c ∈ C(p)

ν c ∈ [- π 2 , π 2 ], cos ν c ≤ min( l c l c+1 (K c n c c • n c c ) 1 2 (K c n c c+1 • n c c+1 ) 1 2 , l c+1 l c (K c n c c+1 • n c c+1 ) 1 2 (K c n c c • n c c ) 1 2
),

where ν c is the angle defined by (73). If the mean cells temperature are such that T c ∈ [0, Θ] for all c ∈ C(p) then the half-edge temperatures satisfy T c ∈ [0, Θ] for all c ∈ C(p).

Let us notice that these geometric conditions are quite difficult to use since they also depend on the local value of the conductivity tensor. We will see in the next paragraph that they are easier to interpretate in the case of an isotropic conductivity.

The case of isotropic conductivity

It is interesting to investigate the particular case of an isotropic conductivity, i.e., K c = κ c I d , where κ c > 0 denotes the piecewise constant approximation of the scalar conductivity κ in cell c. If θ c denotes the measure of the angle between edges c and c + 1 then n c c • n c c+1 = -cos θ c . It follows from (61) that the cth row of matrix M reduces to

     M c,c-1 = -α c-1 κ c-1 l c-1 l c cos θ c-1 , M c,c = α c-1 κ c-1 l 2 c + α c κ c l 2 c , M c,c+1 = -α c κ c l c l c+1 cos θ c . (76)
It turns out that M is a symmetric matrix. Using the result of the previous paragraphs, we know that M is definite positive. However, it is quite instructive to demonstrate this directly. Recalling that θ c = kπ, where k is an integer, implies that M c,c >| M c,c-1 | + | M c,c+1 | for all c ∈ C(p), thus M is strictly diagonally dominant. Noticing that all the diagonal entries of M are positive real numbers allows to claim that this matrix is nonsingular and all its eigenvalues are positive [START_REF] Varga | Matrix Iterative Analysis[END_REF]. In addition, the off-diagonal entries are negative provided that θ c ∈ [-π 2 , π 2 ] for all c ∈ C(p). Namely, M is an M-matrix provided that the previous angular condition holds. Let us note that we have recovered the angular condition (74). Indeed by setting K c = κ c I d into the definition of ν c (73), we find that ν c coincides with θ c . Bearing this in mind, we deduce from the previous paragraph the isotropic version of the maximum principle for the half-edge temperatures: under the geometric conditions if the mean cell temperatures are such that 0 < T c ≤ Θ for all c ∈ C(p), then the half-edge temperatures satisfy the inequality 0 < T c ≤ Θ, for all c ∈ C(p), where Θ > 0 is a given temperature. Note the second condition ensure that the entries of S are non-negative. This condition could have been obtained proceeding directly with S entries knowing that in the isotropic case they simply write

∀c ∈ C(p),      θ c ∈ [- π 2 , π 2 ], cos θ c ≤ min( l c+1 l c , l c l c+1 ), (77)
S c,c-1 = α c-1 κ c-1 l c (l c -l c-1 cos θ c-1 ), S c,c = α c κ c l c (l c -l c+1 cos θ c ). ( 78
)
To investigate further the impact of the geometric conditions, being given the half-edge lengths l c and l c+1 , let us introduce the limit angle θ lim c = arccos[min( l c+1 lc , lc l c+1 )]. Condition (77) is equivalent to the requirement that

θ c ∈ [-π 2 , -θ lim c ] ∪ [θ lim c , π 2
]. In Figure 7, we have plotted the limit angle variation in terms of the sub-cell aspect ratio ζ = lc l c+1 . Note that we have only displayed the positive value knowing that the negative value is obtained by symmetry about zero. The permitted values of the sub-cell angle θ c to ensure that the geometric conditions hold is the domain delimited at the bottom by the curve θ c = θ lim c (ζ) and at the top by the straight line θ c = π 2 . This graph shows that for moderate aspect ratios the geometric condition (77) is not too much restrictive. However, for high aspect ratios, this conditions becomes very restrictive in the sense that the permitted angles are closed to right angles.

Comment 5. We note that in the case of a rectangular grid, the present scheme recovers the well known five-point scheme. In this case, each vertex of the grid is surrounded by four rectangular grids, thus θ c = π 2 for all c = 1 . . . 4. Defining the volume weights w pc according to (53) yields α c = 1 lcl c+1 . We obtain that M reduces to the diagonal matrix defined as

M c,d = l c ( κ c-1 l c-1 + κ c l c+1 )δ c,d .
Regarding S, its cth row writes as

S c,c-1 = l c κ c-1 l c-1 , S c,c = l c κ c l c+1
.

The half-edge temperatures are readily deduced from the solution of the diagonal 4 × 4 linear system (60)

T c = κ c-1 l c-1 T c-1 + κc l c+1 T c κ c-1 l c-1 + κc l c+1 .
Substituting the above result in the half-edge flux expression (57a) leads to

q c c = 1 l c-1 κ c-1 + l c+1 κc (T c -T c-1 ).
We have recovered the expression of the normal flux which corresponds to the classical five-point scheme characterized by a weighted harmonic averaging of the scalar conductivity at the cell interface.

Construction of the diffusion matrix

After having expressed the half-edge temperatures in terms of the mean cell temperatures, we are now in position to achieve the construction of the scheme by gathering the previous results.

Local diffusion matrix at a generic point

We start by deriving the local diffusion matrix at a generic point p. To this end, we first recall the semi-discrete scheme which has been obtained in Section 2.3

m c C vc d dt T c + p∈P(c) l - pc q - pc + l + pc q + pc = m c r c .
In this equation, the half-edge fluxes (q - pc , q + pc ) attached to sub-cell ω pc , express in terms of the halfedge temperatures by means of [START_REF] Potier | Schéma volumes finis monotone pour des opérateurs de diffusion fortement anisotropes sur des maillages de triangles non structurés[END_REF]. In addition, the half-edge temperatures express in terms of the mean cell temperatures surrounding point p through the use of the solution of the linear system (67). Therefore, to write the semi-discrete scheme in terms of the primary unknowns, that is, the mean cell temperatures, it remains to substitute the expression of the half-edge temperatures in terms of the mean cell temperatures into the half-edge normal fluxes. To perform this substitution, it is convenient to define the contribution of the sub-cell ω pc to the diffusion flux as

Q pc = l - pc q - pc + l + pc q + pc . (79)
Using the local notation at point p introduced in Section 3.3, it turns out that Q pc rewrites as

Q pc = l c q c c + l c+1 q c c+1 .
Using the expression of the half-edge fluxes in terms of the half-edge temperatures (57) yields

Q pc = -α c l c [l c (K c n c c • n c c ) + l c+1 (K c n c c • n c c+1 )](T c -T c ) (80) -α c l c+1 [l c (K c n c c+1 • n c c ) + l c+1 (K c n c c+1 • n c c+1 )](T c+1 -T c ).
Recalling the definition of the matrix S from (71)

Sc,c-1 = α c-1 l c [l c-1 (K c-1 n c-1 c • n c-1 c-1 ) + l c (K c-1 n c-1 c • n c-1 c )], Sc,c = α c l c [l c (K c n c c • n c c ) + l c+1 (K c n c c • n c c+1 )],
leads to recast (80) under the more compact form

Q pc = -Sc,c (T c -T c ) -Sc+1,c (T c+1 -T c ).
Due to the sparse structure of S, this last equation turns into

Q pc = - d∈C(p) St c,d (T d -T c ).
Finally, recalling that the vector of half-edge temperatures, T , is expressed in terms of the vector of the mean cell temperatures, T , through the use of the solution of the linear system (67), i.e., T = (M -1 S)T , allows to write Q pc as

Q pc = - d∈C(p) Γ p c,d (T d -T c ). (81)
Here, Γ p denotes the C p × C p matrix defined at point p by

Γ p = St M -1 S. (82)
In deriving (81), we have used that d∈C(p) Γ p c,d = d∈C(p) St c,d which follows from the fact that 1 C(p) is the eigenvector of M -1 S associated to the eigenvalue 1, refer to Comment 4. Let us note that the entries of Γ p have the physical dimension of a conductivity. Thus, Γ p can be viewed as the effective conductivity tensor at point p. More precisely, it follows from (81) that the entry Γ p c,d stands for the effective conductivity between cells c and d through the point p. We claim that the effective conductivity tensor at point p, Γ p , is symmetric positive definite provided that the physical conductivity tensor is symmetric positive definite. To prove this, observe that

Γ p T • T =( St M -1 S)T • T =M -1 (ST ) • ( ST ).
Since K c is symmetric positive definite, by virtue of (70) one deduces that S = S, in addition M is symmetric positive definite, which ends the proof.

In the general case, for which the physical conductivity tensor, K c , is only positive definite, we make the conjecture that Γ p is also positive definite since we are not able to prove this result directly.

Comment 6. We want to mention that in the case of a symmetric positive definite conductivity tensor, under the geometrical conditions (74) and (75), the entries of the matrix Γ p are nonnegative.

The global diffusion matrix

Gathering the previous results, leads to write the finite volume semi-discrete scheme over cell c as follows

m c C vc d dt T c - p∈P(c) d∈C(p) Γ p c,d (T d -T c ) = m c r c , (83)
where Γ p is the effective conductivity tensor defined at point p by (82). To put the above equation under a more compact form, let us introduce the following global notation. Let T be the vector of the cell averaged temperatures, that is, T = (T 1 , . . . T C D ) t , where C D denotes the total number of cells covering the domain D. Let us denote by M and C v the C D × C D diagonal matrices whose entries are given by m c δ c,d and C vc δ c,d . We also introduce R = (r 1 , . . . , r C D ) t as the source term vector. Finally, let A be the C D × C D matrix which stands for the global diffusion matrix associated to the above semi-discrete scheme. Observing (83), we readily deduce that its diagonal and off-diagonal entries write as

A c,c = p∈P(c) d∈C(p) Γ p c,d , (84a)
A c,d = - p∈P(c) Γ p c,d . (84b)
Bearing this notation in mind, our semi-discrete finite volume scheme reads

v dT dt + AT = MR. (85)
This results in a differential system satisfied by the vector of the cell averaged temperatures. Note that the above system has been derived without taking into account the boundary conditions. In writing the entries of the global diffusion matrix, we have to pay attention to the fact that indices in the left-hand side of (84) refer to the global numbering of the cells, whereas in the righthand side they refer to the local numbering of cells surrounding point p. In addition, index p refers to the local numbering of points belonging to cell c.

Let us point out that the global diffusion matrix results in assembling the small1 node-based C p × C p matrices Γ p . This node-based underlying data structure allows to handle easily general unstructured grids. However, the assembling of the global diffusion matrix requires the knowledge of the local matrix Γ p at each grid point p. This matrix is computed as Γ p = St M -1 S where matrices S, M and S are sparse C p × C p matrices explicitly given by formulas (71), ( 62) and (61). Within the framework of a time-marching algorithm, these matrices have to be stored at the beginning of each time step for each grid point. In addition, the computation of the inverse matrix M -1 is performed using an efficient algorithm well adapted to cyclic tridiagonal matrices [START_REF] Press | Numerical Recipes in Fortran 77[END_REF].

The stencil of the finite volume discretization (83) results directly from the structure of the above diffusion matrix. Being given a point p, its surrounding cells c and d are connected through point p by means of the diffusion exchange term Γ p c,d . Therefore, the stencil of cell corresponds to the set of neighboring cell d which shares a point with itself, refer to Figure 8. For a quadrangular grid, this results in a nine-point scheme.

Properties of the semi-discrete scheme

In this section, we state two fundamental properties which characterize our semi-discrete finite volume scheme. The first one concerns the positive semi-definiteness of matrix A and the second one consists of the L 2 -stability of the semi-discrete scheme.

Positive semi-definiteness of the global diffusion matrix

We claim that the global diffusion matrix is positive semi-definite, that is,

AT • T ≥ 0, ∀T ∈ IR C D . (86)
Let us emphasize that this statement will be of fundamental importance to ensure the solvability of the linear system associated to the time discretziation of (85). This result follows directly from the discrete inequality (42) satisfied by the numerical approximation of the half-edge fluxes [START_REF] Potier | Schéma volumes finis monotone pour des opérateurs de diffusion fortement anisotropes sur des maillages de triangles non structurés[END_REF].

To prove (86), let us remark that the cth entry of vector AT writes as

(AT ) c = p∈P(c) (l - pc q - pc + l + pc q + pc ).
Using the above equation allows to rewrite the left-hand side of (86) as

AT • T = C D c=1 p∈P(c) (l - pc q - pc + l + pc q + pc )T c .
Now, switching round the order of summation in the right-hand side leads to

AT • T = P D p c∈C(p) (l - pc q - pc + l + pc q + pc )T c , = P D p I p ,
where P D denotes the total number of points inside D and I p is given by [START_REF] Press | Numerical Recipes in Fortran 77[END_REF]. By virtue of ( 42), the nodal quantity I p is non-negative, thus the right-hand side of the above equation is also non-negative, which ends the proof. We conclude this paragraph by claiming that the global diffusion matrix inherits some properties of the small matrix Γ p . Namely, if Γ p is symmetric positive definite2 , is also symmetric positive but only semi-definite since by construction for a given row, the diagonal entry is equal to the sum of the off-diagonal entries, refer to (84). In addition, under the geometrical conditions ( 74) and (75), the diagonal entries of A are non-negative whereas the off-diagonal entries are non-positive, refer to Comment 6.

L 2 -stability of the semi-discrete scheme

In this section, we prove the stability of our semi-discrete finite volume scheme, in absence of source term, i.e., r = 0, in the sense of the discrete weighted L 2 norm defined as follows

T 2 w2 = C D c m c C vc T 2 c , (87)
where m c and C vc denote the piecewise constant approximation over cell c of the mass and the specific heat capacity, which are positive quantities. To study the L 2 stability of the semi-discrete scheme, we first recall that in absence of source term, vector T satisfies the following semi-discrete system

MC v dT dt + AT = 0,
where M and C v are the diagonal mass and heat capacity matrix. Dot-multiplying the above equation by T yields

MC v dT dt • T + AT • T = 0.
Knowing that A is positive semi-definite leads to

d dt (MC v T • T ) ≤ 0. ( 88
)
In writing this equation, we have supposed that M and C v does not depend on time. In addition, we are ignoring the contribution of the boundary terms, assuming for instance periodic boundary conditions. Using the definition (87) of the weighted L 2 norm, the above equations rewrites as

d dt ( 1 2 T 2 w2 ) ≤ 0. ( 89
)
This inequality ensures that the weighted L 2 norm of the semi-discrete solution remains bounded by the weighted L 2 norm of the initial data, which corresponds to the L 2 -stability of our semi-discrete finite volume scheme. This L 2 -stability is a direct consequence of the half-edge normal fluxes construction through the use of the sub-cell variational formulation. Once more, this shows the great importance of deriving the numerical approximation of the normal fluxes using a variational formulation.

p c T c+1 n c c+1 ω c T c c + 1 C p + 1 1 T 1 T C p +1 q ⋆ C p +1 T C p +1 n 1 1 T 1 q ⋆ 1 T c q c c q c c+1 n c c ω pc n C p C p +1
Figure 9: Fragment of an unstructured grid in the vicinity of a boundary point where are imposed Neumann boundary conditions. The prescribed fluxes, q ⋆ 1 and q ⋆ Cp+1 are displayed using green color.

Boundary conditions implementation

This section describes indications related to boundary conditions implementation. Let us point out that boundary conditions treatment relies on a straightforward extension of the half-edge temperatures elimination procedure which has been developed in Section 3.3 for internal nodes. The boundary conditions are prescribed at the two boundary half-edges connected to a point located on the boundary, they will be either of Neumann type (prescribed normal flux) or of Dirichlet type (prescribed temperature). Here, we are going to expose the main steps of the Neumann boundary conditions discretization without going into the detail of the computations. Concerning the Dirichlet boundary conditions discretization the interesting reader can refer to [START_REF] Breil | A cell-centered diffusion scheme on two-dimensional unstructured meshes[END_REF].

Let us consider a boundary node, p, which is surrounded by internal cells. As before, we denote by C p the number of cells surrounding point p. Note that the number of half-edges impinging on point p is equal to C p + 1. The edges impinging at point p are labelled using the subscript c ranging from 1 to C p + 1. The cell (sub-cell) numbering follows the edge numbering, that is, cell ω c (sub-cell ω pc ) is located between edges c and c + 1, refer to Figure 9. The unit outward normal to cell ω c at edge c is denoted by n c c whereas the unit outward normal to cell ω c at edge c + 1 is denoted by n c+1 c . The prescribed heat fluxes on the first and the last half-edge are denoted respectively by q ⋆ 1 and q ⋆ Cp+1 . Let us recall that the half-edge fluxes corresponding to the internal edges are expressed in terms of the half-edge and the mean cells temperatures using for c = 2, . . . ,

C p q c-1 c = -α c-1 [l c-1 (K c-1 n c-1 c-1 • n c-1 c )(T c-1 -T c-1 ) + l c (K c-1 n c-1 c • n c-1 c )(T c -T c-1 )], q c c = -α c [l c (K c n c c • n c c )(T c -T c ) + l c+1 (K c n c c+1 • n c c )(T c+1 -T c )]. (90)
Here, q c c (resp. q c-1 c ) denotes the half-edge normal flux at edge c viewed from cell (resp. c -1). In addition, T c is the temperature on the cth half-edge whereas T c corresponds to the mean temperature of cell c.

The elimination of the half-edge temperatures is obtained by writing the flux continuity conditions (59) for all internal half-edges. This system of C p -1 equations is completed by the two boundary conditions for the first and last half-edges. Finally, this results in the following system

l 1 q 1 1 = l 1 q ⋆ 1 , (91a)
l c q c-1 c + l c q c c = 0, for c = 2, . . . , C p , (91b)
l C+1 q C C+1 = l C+1 q ⋆ C+1 . (91c)
Substituting (90) into (91b) leads to the following (C p + 1) × (C p + 1) linear system satisfied by the

C p + 1 half-edges temperatures MT = ST -B. (92)
Here, T ∈ IR Cp+1 denotes the vector of half-edge temperatures whereas T ∈ IR Cp is the vector of cell centered temperatures. In addition, M and S are respectively a (C p + 1) × (C p + 1) tridiagonal matrix and a (C p + 1) × C p bidiagonal matrix. Their entries are computed developing (91). Boundary conditions are taken into account by means of B ∈ IR Cp+1 with B = (l 1 q ⋆ 1 , . . . , l Cp+1 q ⋆ Cp+1 ) t . Proceeding with the matrices M and S as before, one can show that (92) always admits a unique solution. This allows to compute the contribution of the sub-cell ω pc to the diffusion flux as

Q pc = - d∈C(p) Γ p c,d (T d T c ) + ( SM -1 B) c . (93)
Here, Γ p denotes the effective conductivity tensor at boundary point p, it is a C p ×C p matrix defined by Γ p = St M -1 S, where the matrix S is deduced from S by means of (70). From the above equation, we collect the contributions to the entries of the global diffusion matrix, A, following

A c,c = p∈P(c) d∈C(p) Γ p c,d , A c,d = - p∈P(c) Γ p c,d .
Having taken into account the boundary conditions, the semi-discrete finite volume scheme (85) turns into

MC v dT dt + AT = MR + S,
where S ∈ IR C D is the source term vector which represents the boundary conditions contribution.

According to (93), its cth entry reads as S c = -( SM -1 B) c . We conclude by stating that using similar arguments to those employed in Section 3.5, it is possible to show that matrix A is still positive semi-definite. Moreover, for homogeneous Neumann boundary conditions, i.e., q ⋆ 1 = q ⋆ Cp+1 = 0, the L 2 -stability of the semi-discrete scheme still holds.

Space discretization of the CCLADNS scheme

Here, we describe briefly the space discretization of the CCLADNS scheme which is obtained by employing a numerical expression of the half-edge normal fluxes resulting from a finite difference approximation. Our main motivation is to derive a finite difference approximation of the half-edge fluxes which preserves linear solutions regardless the shape of the cell. Doing so, we enforce the accuracy of our space discretization. However, since this approximation does not result from a variational formulation, in general, we cannot insure the transfer of the good properties of the conductivity tensor to the discrete approximation, i.e., its positive definiteness. For instance, it will not be possible to ensure that the inequality (42) holds for a finite difference approximation of the fluxes.

Half-edge fluxes expression by means of a finite difference approximation

Let ω c be an arbitrary cell and ω pc its sub-cell associated to point p. To define a finite difference approximation which is linear preserving, we reuse the piecewise linear approximation of the temperature over ω c introduced in the latter paragraph as

T h (x) = T c -K -1 c q c • (x -x c ), ∀x ∈ ω c .
Let us recall, that x c is the centroid of the cell, i.e., x c = 1 |ωc| ωc x dv, and T c is the cell averaged temperature. In writing the above equation we have used the piecewise constant definition of the heat flux, i.e., q c = -K c (∇T ) c . Let i ± denotes the midpoint of the segment [p, p ± ] and x i ± its corresponding position vector. Assuming that the half-edge temperature ± pc is given by

T ± pc = T h (x i ± ) leads to T - pc -T c = -K -1 c q c • (x - i -x c ), T + pc -T c = -K -1 c q c • (x + i -x c ). (94)
Setting µ ± pc = x ± ix c and using (25), we transform (94) by expressing the vectors q c and µ ± pc in terms of their normal components

T - pc -T c = -K -1 c J -t pc q - pc q + pc • J -t pc µ - pc • n - pc µ - pc • n + pc , T + pc -T c = -K -1 c J -t pc q - pc q + pc • J -t pc µ + pc • n - pc µ + pc • n + pc .
Introducing the sub-cell conductivity tensor K pc defined by [START_REF] Moulton | Workshop: Compatible Spatial Discretizations for Partial Differential Equations[END_REF] and after some manipulations, the above equation rewrites as

T - pc -T c T + pc -T c = - µ - pc • n - pc µ - pc • n + pc µ + pc • n - pc µ + pc • n + pc K -1 pc q - pc q + pc .
The rows of the first matrix in the right-hand side corresponds to the half-edge normal components of vectors µ - pc and µ + pc . This matrix is non-singular provided these two vectors are not colinear. Assuming this and solving the above linear system leads to the final expression of the half-edge normal fluxes

q - pc q + pc = - 1 ∆ pc K pc µ + pc • n + pc -µ - pc • n + pc -µ + pc • n - pc µ - pc • n - pc T - pc -T c T + pc -T c , (95)
where

∆ pc = (µ - pc • n - pc )(µ + pc • n + pc ) -(µ - pc • n + pc )(µ + pc • n - pc )
is the determinant of the matrix defined by the normal components of ± pc . Equation (95) constitutes a finite difference approximation of the half-edge normal fluxes which is linear preserving since it has been deduced from a piecewise linear approximation of the temperature field. However, the matrix form of this approximation reveals that the property transfer, which characterizes the flux approximation derived through the use of variational formulation, does not hold here. This is due to the fact that the matrix defined by the normal components of µ ± pc is not a symmetric positive definite matrix in general. Using such an approximation for the normal flux, leads to a finite volume discretization characterized in general by a non-symmetric discrete diffusion operator, which renders the resolution of the corresponding linear system less obvious. The main advantage in using this scheme lies in the fact that it preserves linear solution regardless the shape of the cells.

Comment 7. It is interesting to investigate further the finite difference approximation (4) in the case where ω pc is a parallelogram. In this particular case, we have µ - pc = -l + pc t + pc and µ + pc = l - pc t - pc where t - pc and t + pc denote the unit tangent vectors to the two half-edges impinging at point p. Bearing this in mind, we have µ - pc • n - pc = -l + pc sin θ pc and µ + pc • n + pc = l - pc sin θ pc . Using these results, the finite difference approximation (4) turns into

q - pc q + pc = - 1 l - pc l + pc sin θ pc K pc l - pc (T pc -T c ) l + pc (T + pc -T c ) . ( 96
)
We have recovered a form which coincides with the one resulting from the variational formulation [START_REF] Potier | Schéma volumes finis monotone pour des opérateurs de diffusion fortement anisotropes sur des maillages de triangles non structurés[END_REF]. In addition, we can identify the volume weight as being w pc = l - pc l + pc sin θ pc , which is exactly the weight (53) introduced in Section 3.2.3. This shows, that the flux approximation resulting from the sub-cell-based variational formulation, with the latter definition of the volume weight, is able to preserve linear solution on quadrangular grids made of parallelograms.

Construction of the CCLADNS scheme

Let us point out that the half-edge normal fluxes expression (95) can be rewritten under the following form

q - pc q + pc = - 1 w pc K pc l - pc (T - pc -T c ) l + pc (T + pc -T c ) . (97)
Here, w pc denotes the volume weight which is defined by w pc = ∆ pc , and the matrix K pc reads

K pc = K pc   µ + pc l - pc • n + pc - µ - pc l + pc • n + pc - µ + pc l - pc • n - pc µ - pc l + pc • n - pc   .
We note that this expression is very similar to the one derived using the local variational formulation, refer to [START_REF] Potier | Schéma volumes finis monotone pour des opérateurs de diffusion fortement anisotropes sur des maillages de triangles non structurés[END_REF]. The only difference lies in the expression and the property of the volume weight and the corner matrix K pc . Using (97), it is possible to perform a local elimination of the half-edge temperatures, around each node, invoking the continuity of the half-edge normal fluxes. This elimination procedure follows exactly the methodology described in Section 3.2. Namely, the halfedge fluxes continuity condition leads to a cyclic tridiagonal linear system. However, it will not be possible to prove that this linear system is invertible due to the lack of positive definiteness of matrix K pc . Finally, collecting the contribution of each node, we construct as previously the global diffusion matrix. Once more, we are not able to prove that this matrix is positive definite due the lack of positive definiteness of the corner matrix K pc .

Extension to cylindrical geometry

The purpose of this section is to present the extension of CCLAD scheme to cylindrical geometry. Note that the extension to cylindrical geometry of the mimetic finite difference method has been derived many years ago in [START_REF] Morel | A local support-operators diffusion discretization scheme for quadrilateral r-z meshes[END_REF]. We note that the case of Cartesian or cylindrical geometry can be combined by introducing the pseudo Cartesian frame (O, x, y), equipped with the orthonormal basis (e x , e y ), through the use of the pseudo radius

O x y e x e y a l(t) n R v = a R da s = l R dl
R(y) = 1 -α + αy,
where α = 1 for cylindrical geometry and α = 0 for Cartesian geometry. We remark that y corresponds to the radial coordinate in the cylindrical case. This means that we assume rotational symmetry about x axis, refer to Figure 10. We note that if we refer to standard cylindrical coordinates, (r, z), then x corresponds to z and y to r. In this framework, the volume v is obtained by rotating the area a about the x axis. Thus, the volume element, dv, writes dv = R da, where da = dxdy is the area element with respect to Cartesian coordinates (x, y). Note that we have omitted the factor 2π due to the integration in the azimuthal direction, namely we consider all integrated quantities to be defined per unit radian. The surface s, which bounds the volume v, is obtained by rotating, l, the boundary of the area a, about the x axis. Thus, the surface element, ds, writes ds = R dl, where dl the line element along the perimeter of a.

In view of subsequent spatial discretization, we shall express the volume integral associated to the divergence operator using the Green formula. We recall that, in the pseudo Cartesian frame, for an arbitrary vector φ ∈ IR 2 , the divergence operator writes

∇ • φ = ∂φ ∂x + 1 R ∂ ∂y (Rψ) = ∂φ ∂x + ∂ψ ∂y + α ψ R = 1 R ∂ ∂x (Rφ) + ∂ ∂y (Rψ) ,
where (φ, ψ) are the components of the vector φ. The gradient operator writes as usual

∇T = ∂T ∂x e x + ∂T ∂y e y .
Let us replace the volume integral form of the divergence operator by its surface integral form, employing the previous notation

v ∇ • φ dv = a 1 R ∂ ∂x (Rφ) + ∂ ∂y (Rψ) R da = a ∂ ∂x (Rφ) + ∂ ∂y (Rψ) da = l φ • nR dl,
where n is the unit outward normal to the contour l. Thus, the Green formula using pseudo Cartesian coordinates reads

v ∇ • φ dv = l φ • nR dl. ( 98
)
Applying this Green formula, the weak form of the heat conduction equation ( 18) rewrites as

d dt ωc ρC v T (x, t) dv + ∂ωc q • nR dl = ωc ρr(x, t)dv, (99)
where ω c denotes a generic cell of the computational domain, ∂ω its boundary and n the corresponding unit outward. Introducing as before a piecewise constant approximation of the physical variables allows to transform (99) into

m c C vc d dt T c + ∂ωc q • nR dl = m c r c , (100)
where the mass of the cell is given as usual by m c = ρ c | ω c |, whereas the cell volume is obtained rotating the polygonal cell ω c about x axis, i.e, | ω c |= ωc dv = ωc R da. In addition, T c = T c (t) denotes the mean cell temperature defined by

T c (t) = 1 | ω c | ωc T (x, t) dv.
To complete the space discretization it remains to discretize the surface integral in the above equation. To this end, we proceed as before introducing the piecewise constant approximation of the normal heat flux over each half-edge impinging on point

q ± pc = 1 l ± pc ∂ω ± pc q • n dl. (101)
The scalar q ± pc stands for the half-edge normal flux related to the half-edge ∂ω ± pc , refer to Figure 2. Note that the above piecewise constant approximation of the half-edge fluxes, coincides exactly with the one used in Cartesian geometry, refer to [START_REF] Hermeline | A finite volume method for approximating 3d diffusion operators on general meshes[END_REF]. Using the partition of the cell into sub-cells, i.e., ∂ω c = ∪ p∈P(c) ∂ω ± pc , the discretized heat conduction equation writes as

m c C vc d dt T c + p∈P(c) R p (l - pc q - pc + l + pc q + pc ) = m c r c , (102)
where R p = 1α + αy p and y p is the y-coordinate of point p. Note that the above equation has been derived by means of the following quadrature rule

∂ω - pc ∪∂ω + pc q • nR dl = R p (l - pc q - pc + l + pc q + pc ).
Now, the space discretization is achieved by using the half-edge flux approximations which have been previously constructed in Section 3.2 and Section 4. Then, the construction of the scheme is performed using the same steps than before. Let us point out that the above quadrature rule has been specifically chosen to ensure spherical symmetry preservation when solving heat conduction equation over equal angle polar grids. In addition, one has to use the volume weight defined for quadrangular cells by [START_REF] Yu | A generalization of the notion of Dirichlet cell for nonconvex domain[END_REF], as it has been already noticed in [START_REF] Morel | A local support-operators diffusion discretization scheme for quadrilateral r-z meshes[END_REF].

Time discretization

In this section, we describe the time discretization related to the finite volume scheme CCLAD. Let us recall that the semi-discrete scheme resulting from the space discretization writes under the form of the following system of differential equations

MC v dT dt + AT = MR + S, ( 103
)
where T is the cell centered temperatures vector, M and C v denote respectively the diagonal mass and heat capacity matrices, whereas A is the global diffusion matrix. In addition, R is the heat supply vector and S is the vector taking into account the prescribed boundary conditions. The above system is completed by prescribing the initial condition T (0) = T 0 . We solve the previous system over the time interval [0, T] using the subdivision 0 = t 1 < t 2 < . . . < t n < t n+1 < . . . < t N = T.

We denote by ∆t n generic time step, i.e., ∆t n = t n+1 -t n . The time approximation of a quantity at time n is denoted using the superscript n, for instance T n = T (t n ). Knowing that explicit time discretization of the diffusion operator requires a stability constraint on the time step which is quadratic with respect to the smallest cell size, we prefer to use an implicit time discretization for this term. Depending on whether the heat capacity and the conductivity tensor depend on the temperature or not, the above system might be linear or not with respect to temperature. This leads us to separate the description of the time discretization in two cases.

Linear case

In this case, we assume that the heat capacity and the conductivity tensor does not depend on temperature. Integrating (103) over the time interval [t n , t n+1 ] yields the first-order time implicit discrete scheme

MC v T n+1 -T n ∆t n + AT n+1 = MR n + S n . ( 104
)
The cell centered temperatures vector is updated by solving the following linear system

( MC v ∆t n + A)T n+1 = MC v ∆t n T n + MR n + S n .
For CCLADS, A is positive semi-definite. Knowing that the entries of the diagonal matrix MC v is always positive, we deduce that matrix MCv ∆t n + A is positive definite which implies that the above linear system always admits a unique solution. Note that, if A is symmetric, the matrix of the linear system is also symmetric. In this latter case, it also interesting to mention that in the absence of source term and for homogeneous boundary conditions, if the geometric conditions (74) and (75) hold, then MCv ∆t n + A is an M-matrix, which means that the positivity of the temperature field is preserved by the scheme.

Non-linear case

In this case, which frequently occurs in plasma physics, the material properties depend on temperature. It is more convenient to rewrite the heat conduction equation using the specific internal energy as

ρ ∂ε ∂t + ∇ • q = ρr.
Let us recall that the specific internal energy, ε, expresses in terms of the density and the temperature by means of an equation of state written under the form ε = ε(ρ, T ). Here, the specific heat capacity is computed as the partial derivative of the specific internal energy with respect to temperature, the density being fixed, i.e., C v = ( ∂ε ∂T ) ρ . The semi-discrete system corresponding to the above partial differential equation writes

M dE dt + AT = MR + S, (105)
where E denotes the specific internal energy vector which depends on ρ and T through the use of the equation of state. Assuming the absence of fluid motion, E exhibits only a non-linear dependency on T . Integrating (105) over the time interval [t n , t n+1 ] leads to the first-order time implicit discretization

M E n+1 -E n ∆t n + A n T n+1 = MR n + S n . (106)
In writing this equation, we have made the choice of an explicit treatment of the non-linear dependency on temperature of the diffusion matrix. Namely, A n corresponds to an evaluation of the diffusion matrix wherein the conductivity tensor has been computed using the temperature at the beginning of the time step. The above equation results in a non-linear system of differential equations, which requires an iterative method to be solved. This iterative method consists in defining a sequence of approximations of (106) by using a Newton-like method. To this end, let us denote by T q the sequence of temperatures vector, where q is natural integer. We initialize this sequence by setting for q = 0, T 0 = T n . Introducing the increment of the temperatures vector as ∆T = T q+1 -T q allows to define the specific internal energies vector increment as

E(T q+1 ) = E(T q ) + C v (T q )∆T . ( 107
)
This equation has been obtained through the use of the first-order Taylor expansion

ε[T q c + (∆T ) c ] = ε(T q c ) + ( ∂ε ∂T ) ρ (T q c )(∆T ) c ,
where (∆T ) c is the cth entry of vector ∆T . Substituting, the Taylor expansion (107) into (106) leads to the following linear system satisfied by ∆T

M ∆t n C v (T q ) + A n ∆T = - M ∆t n [E(T q ) -E n ] + A n T q -MR n -S n . ( 108
)
Let us point out that this linear system admits always a unique solution since the matrix between bracket in the left-hand side is positive definite. The updated value of the sequence, T q+1 , is computed as T q+1 = T q + ∆T by solving the above linear system in ∆T . This process is repeated until a sufficiently accurate solution is reached. More precisely, we assess the convergence of the iterative method by computing the value of a residual characterizing the accuracy at which the non-linear equation is solved. A relevant choice for the residual consists in defining it as being equal to a certain norm of the left-hand side. Setting

R q = M ∆t n [E(T q ) -E n ] + A n T q -MR n -S n
, we define the stopping criterion of the iterative procedure as

R q R 0 ≤ η, ( 109
)
where η is a fixed in advance positive real number. The usually employed norm is the maximum norm, i.e., R q ∞ = max c (R q c ). The main advantage in using this stopping criterion lies in the fact that it provides an indication on the accuracy at which the energy conservation is ensured.

Comment 8. In both linear and non-linear cases, the numerical solution requires to solve a sparse linear system. If the matrix related to the linear system is symmetric positive definite, we make use of an Incomplete Cholesky Conjugate Gradient (ICCG) [START_REF] Lascaux | Analyse Numérique matricielle appliquée à l'art de l'ingénieur[END_REF] to solve the corresponding linear system. On the other hand, if the matrix is non-symmetric, the linear system is solved by means of the Generalized Minimal Residual method (GMRES) available in the PETSc library [START_REF] Saad | GMRES: a Generalized Minimal RESidual algorithm for solving nonsymmetric linear systems[END_REF][START_REF]Summary of sparse linear solvers available from petsc[END_REF].

Numerical results

The aim of this section is to assess the robustness and the accuracy of CCLAD scheme against analytical test cases using various types of triangular, quadrangular and polygonal grids. We will present results obtained using two versions of CCLAD scheme, that is the symmetric version named CCLADS wherein the flux discretization is derived using the sub-cell-based variational formulation, refer to Section 3.2 and the non-symmetric version for which the flux discretization is obtained by means of a finite difference approximation, refer to Section 4. After evaluating the accuracy of CCLADS and CCLADNS using various test cases in Cartesian geometry, we show numerical results related to the extension of both schemes to cylindrical geometry. Finally, we conclude this section by presenting two tests which are not very far from the problems encountered in the context of the numerical simulation of Inertial Confinement Fusion. We also want to mention that CCLADS scheme has been successfully used within a set of two-dimensional hydrodynamic simulations, which were performed to reproduce experimental measurements resulting from laser plasma experiments, wherein electron heat transport was strongly modified by self-generated magnetic fields, refer to [START_REF] Schurtz | Revisiting nonlocal electron-energy transport in inertial-fusion conditions[END_REF].

Methodology for convergence analysis

Most of the test cases presented in this section are performed using a standard test problem which consists in solving the following diffusion equation over the domain D

= [0, 1] 2 ρC v ∂T ∂t -∇ • (K∇T ) = ρr, (110a)
T (x, 0) = T 0 (x), (110b)
where r = r(x) is a source term. Most of the analytical solutions being stationary, we compute them starting with the initial condition T 0 (x) = 0, and we run the numerical simulation until the steady state is reached. Density and the specific heat capacity are specified such that ρ = 1 and C v = 1. The boundary conditions, the source term and the heat conductivity tensor, K, will be prescribed for each test case. Bearing this in mind, let us describe the methodology used to perform the convergence analysis. Knowing that the computational domain is paved using C D cells, we define the mesh resolution as

h = | D | C D ,
where | D | denotes the domain volume. Let T = T (x) be the steady analytical solution of (110).

Being given a computational grid characterized by h, we denote by T h c the value of the analytical solution computed at the centroid of the cell ω c , i.e., T h c = T (x c ). If T h c denotes the cell averaged temperature corresponding to the numerical solution obtained by the finite volume scheme, we define the asymptotic numerical errors based on the maximum norm and the l 2 norm as

E h max = max c=1...C D | T h c -T h c |, (111a)
E h l 2 = C D c=1 (T h c -T h c ) 2 | ω c |. (111b)
asymptotic error for both norms is estimated by

h α = C α h qα + O(h qα+1 ), for α = max, l 2 ,
where q α denotes the order of truncation error and C α the convergence rate-constant which is independent of h. Having computed the asymptotic errors corresponding to two different grids characterized by mesh resolutions h 1 and h 2 < h 1 , we deduce an estimation of the order of truncation error as

q α = log( E h 1 α E h 2 α ) log( h 1 h 2 ) . (112)

Anisotropic linear problem with discontinuous conductivity tensor

This problem consists in finding the steady solution of (110) with r = 0 and an anisotropic discontinuous conductivity tensor given by transformation of Cartesian grids. Following the approach described in [START_REF] Shashkov | Solving Diffusion Equations with Rough Coefficients in Rough Grids[END_REF], we first introduce the smooth distorted grids resulting from the mapping defined on the unit square by

K(x, y) =                K xx l K xy l K yx l K yy l if 0 ≤ x ≤ 1 2 , K xx r K xy r K yx r K yy r if 1 2 ≤ x ≤ 1.
x(ξ, η) = ξ + a 0 sin(2πξ) sin(2πη), y(ξ, η) = η + a 0 sin(2πξ) sin(2πη). (113)
The three smooth grids resulting from this mapping with a 0 = 0.1 are displayed in Figure 12. We also define randomly distorted grids by means of the following mapping defined on the unit square by x(ξ, η) = ξ + a 0 h cos(2πθ),

y(ξ, η) = η + a 0 h sin(2πθ), (114)
where h corresponds to the mesh spacing of the initial Cartesian grid and θ is a random number chosen in [0, 1]. We have plotted the three random grids obtained with this mapping for a 0 = 0.2 in Figure 13. Note that we do not apply this mapping to the nodes located on the line ξ = 0.5 in order to preserve the interface. The convergence analysis for smooth grids is performed computing the asymptotic errors and the corresponding orders of truncation error using formulas (111) and (112).

The results displayed in Table 1(a) show that the convergence rate is almost second-order in the l 2 norm and a little bit less in the maximum norm. Proceeding with the convergence analysis for random grids as before, we have displayed the corresponding results in Table 1(b). We observe an erratic behavior regarding the asymptotic errors and the rate of convergence in both norms which clearly shows a lack of convergence for our scheme with this type of random grids. Note that this behavior has been already observed in the case of isotropic diffusion test cases [START_REF] Breil | A cell-centered diffusion scheme on two-dimensional unstructured meshes[END_REF].

Polygonal grids

Now, we study the convergence of CCLADS scheme using a sequence of three polygonal grids displayed in Figure 14. The volume weight used in the numerical approximation of the half-edge

h E h max q h max E h l 2 q h l 2
1.00D-1 1.2D-2 0.79 3.02D-3 0.73 5.00D-2 6.9D-3 0.61 1.82D-3 0.74

2.50D-2 4.5D-3 - 1.09D-3 - (c) Polygonal grids. h E h max q h max E h l 2 q h l 2 1.00D-1 1.75D-2 0.91 5.91D-3 1.19 5.00D-2 9.33D-3 1.00 2.58D-3 1.25 2.50D-2 4.65D-3 - 1.08D-3 - CCLADS CCLADNS h E max q max E h l 2 q l 2 E h max q max E h l 2 q l 2 6.27D-2 9.92D-

Triangular grids

We run this test problem on a sequence of three triangular grids which are displayed in Figure 15. The convergence analysis results corresponding to the numerical simulations using the three triangular grids and both schemes CCLADS and CCLADNS are displayed in Table 2. First, they show that CCLADS scheme has a second-order convergence rate in l 2 norm on triangular grids. Second, we observe that CCLADNS is not converging, we point out that its order of convergence in maximum norm becomes negative. This fact is probably a consequence of the oscillatory behavior of this scheme. We point out that CCLADS scheme satisfies a L 2 -stability property whereas CCLADNS does not. This fact is clearly illustrated by the numerical results displayed in Figure 16 where we can observe the spurious oscillations produced by CCLADNS scheme, refer to Figure 16(b).

Quadrangular grids

Concerning the quadrangular grids we perform the convergence analysis on three types of grids: rectangular, smooth and random. We start by giving in Table 3(a) the convergence analysis data for a sequence of three rectangular grids. These data demonstrate that both CCLADS and CCLADNS schemes exhibit a second-order rate of convergence on rectangular grids.

CCLADS CCLADNS

h E h max q max E h l 2 q l 2 E h max q max E h l 2 q l 2 1.00D-1 3.97D-2 2.

CCLADS CCLADNS

h E h max q max E h l 2 q l 2 E h max q max E h l 2 q l 2 1.00D-1 1.09D-

CCLADS CCLADNS

h E h max q max E h l 2 q l 2 E h max q max E h l 2 q l 2 1.00D-1 7.88D-2 1.12

CCLADS CCLADNS

h E h max q max E h l 2 q l 2 E h max q max E h l 2 q l 2 1.00D-1 4.87D-2 1.67 1.69D-2 1.91 3.83D-2 1.94 1.58D-2 1.95 5.00D-2 1.53D-2 1.03 4.5D-3 1.91 9.96D-3 1.99 4.1D-3 1.92 2.50D-2 7.51D-3 - 1.2D-3 - 2.51D-3 - 1.08D-3 -
Table 5: Anisotropic linear problem with a non-uniform non-symmetric conductivity tensor: asymptotic errors in both maximum and l 2 norms and corresponding truncation error orders for quadrangular grids.

(a) Smooth grids. The boundary conditions prescribed at the top and bottom boundaries of the computational domain are Dirichlet boundary conditions deduced from the above analytical solution. For numerical applications, we choose a = 1 and b = 2. This test problem is used to perform a convergence analysis for both CCLADS and CCLADNS schemes employing three different types of quadrangular grids. We reuse the smooth and the random grids which has been introduced previously, refer to Figures 12 and17. We also make use of the Kershaw grids [START_REF] Kershaw | Differencing of the Diffusion Equation in Lagrangian Hydrodynamic Codes[END_REF], which are highly skewed grids displayed in Figure 19. For each type of grids, we perform a sequence of three computations increasing the mesh refinement. The resulting asymptotic errors and rate of convergence in both maximum and l 2 norms are displayed respectively in Tables 5(a) and 5(b) for smooth and Kershaw grids. From these results it follows that both schemes exhibit a rate of convergence which is located between first-order and secondorder for these types of grids. Note that the rate of convergence of CCLADNS scheme is the greatest. The results associated to the sequence of random grids are displayed in Table 5(c). Once more, they reveal an erratic behavior of CCLADS scheme regarding its convergence on random grids whereas CCLADNS scheme is converging whith a rate which is quite close to second-order.

CCLADS CCLADNS

h E h max q max E h l 2 q l 2 E h max q max E h l 2 q l 2 1.00D-1 1.64D-2 1.

CCLADS CCLADNS

h E h max q max E h l 2 q l 2 E h max q max E h l 2 q l 2 1.67D-1 8.99D-2 1.44 3.63D-2 1.54 5.91D-2 1.35 2.1D-2 1.65 5.55D-2 1.83D-2 1.67 6.66D-3 1.78 1.34D-2 1.71 3.4D-3 1.89 2.78D-2 5.65D-3 - 1.95D-3 - 4.11D-3 - 9.23D-4 - (c) Random grids. CCLADS CCLADNS h E h max q max E h l 2 q l 2 E h max q max E h l 2 q l 2 1.00D-1 1.14D-2 0.

CCLADS CCLADNS

h E h max q max E h l 2 q l 2 E h max q max E h l 2
q l 2 7.77D-2 7.71D-2 0.94 3.18D-2 0. The isotropic heat conductivity, κ, is a non-linear function with respect to temperature. Here, we set κ(T ) = T 5 2 , this type of non-linearity corresponds to the so-called Spitzer-Härm conductivity which is frequently used in plasma physics to describe electron thermal heat flux in the local regime, refer to [START_REF] Ya | Physics of Shock Waves and High-Temperature Hydrodynamic Phenomena[END_REF].

The above diffusion equation is solved over a cylindrical domain D defined by r ∈ [0, 1] and θ ∈ [0, π 2 ], where (r, θ) denote the classical polar coordinates, which express in terms of the Cartesian coordinates (x, y) as r = x 2 + y 2 and θ = arctan y

x . The mass density and the heat capacity are given by ρ = 1 and C v = 1. The initial condition is defined by T 0 (x) = 1. We prescribe symmetry boundary conditions along axis x = 0 and y = 0 and we impose the normal flux q ⋆ = 1000 at the outer radius r = 1. The unsteady solution is computed until time T = 3. 10 -3 . At this time a non-linear heat wave has propagated into the cold medium. This wave is characterized by a sharp transition zone displaying a strong temperature gradient. Due to the boundary conditions and the geometry of the domain, the solution of the diffusion equation exhibits a cylindrical symmetry, namely T (x) ≡ T (r).

Unfortunately, the non-linear equation under consideration does not admit exact analytic solutions. That is why, we compute with our finite volume scheme a reference numerical solution using a 50 × 40 polar grid which is displayed in Figure 22 . Now, we compute the numerical solution of the non-linear test problem on the above distorted grid, using two different schemes. These are: the classical five-point scheme and the CCLADS scheme which reduces to a nine-point scheme on quadrangular grids. The numerical solution resulting from the five-point scheme is plotted in Figure 23 using blue dots. We have displayed the temperatures in all cells as function of the cell center radius versus the reference solution. In this figure, we observe the main flaw of the five-point scheme: in spite of its robustness, it produces a numerical solution wherein the temperature front is aligned with the grid distortion. The corresponding numerical solution is not able to preserve the cylindrical symmetry. In addition, the comparison to the reference solution, shows that the timing of the thermal wave is completely wrong. Let us emphasize that this test case is not a fake problem. It is representative of situations which frequently occur in the framework of plasma physics simulation wherein the heat conduction equation is coupled with a numerical method solving Lagrangian hydrodynamics equations. In this case, grid distortions are induced by fluid motion and thus the use of the five-point scheme to solve the heat conduction equation leads to a very bad result. This weakness of the five-point scheme follows from the fact that its construction is based on a two-point flux approximation 3 , which becomes inaccurate in the 0.00e+000 -3.33e+006 -6.67e+006 -1.00e+007 According to [START_REF] Dautray | La fusion thermonucléaire par laser[END_REF] the anisotropic Braginskii conductivity tensor, K e , writes

K e = κ ⊥ κ ∧ -κ ∧ κ ⊥ ,
where the Braginskii transport coefficients are expressed in terms of the Spitzer-Härm conductivity by means of functions describing the magnetization of the heat flux, refer to Figure 1 and [START_REF] Braginskii | Transport Process in a Plasma[END_REF]. To better understand the action of the above anisotropic conductivity tensor, observe that for any arbitrary vector φ ∈ IR 2 K e φ = κ ⊥ φκ ∧ R π 2 φ, where R π 2 is the counterclockwise rotation through the angle π 2 . This equation shows that the action of K e decomposes in an isotropic part, which corresponds to a multiplication by κ ⊥ , completed by a rotation. This last term follows directly from the magnetic field since κ ∧ is directly proportional to the magnitude of the magnetic field. This corresponds to the so-called Righi-Leduc effect.

Bearing this in mind, we proceed to compute numerical solutions of the above problem using our finite volume scheme. The computational domain has been paved using a 50×50 Cartesian grid and stopping time is T = 25 10 -12 s. We prescribe Neumann homogeneous boundary conditions at the boundaries of the computational domain. The first computation has been done suppressing the Righi-Leduc term, that is, setting κ ∧ = 0. In this case the conductivity tensor reduces to the isotropic form K e = κ ⊥ I d . The corresponding temperature contours are displayed in Figure 25(a) wherein we can observe a leftward thermal wave propagating from hot to cold region. Note that the propagation of the thermal wave is strongly inhibited by the presence of the magnetic field. This follows from the fact that the Braginskii coefficient κ ⊥ is a decreasing function with respect to the magnitude of magnetic field, refer to Figure 1. The second computation is performed using the complete form of the anisotropic conductivity tensor, that is, including the Righi-Leduc effect. This term, as it can be observed in Figure 25(b), rotates the temperature contours.

Conclusion

In this paper, we have presented a high-order cell-centered finite volume scheme called CCLAD, which aims at solving anisotropic diffusion equation on two-dimensional distorted grids. This scheme is characterized by a local stencil and cell-centered unknowns. Its main feature lies in the introduction of two half-edge normal fluxes and two half-edge temperatures per cell interface. Two versions of this scheme are available. In the first version, named CCLADS, the half-edge fluxes expression is derived using a local variational formulation whereas in the second one, called CCLADNS, it is obtained by means of a finite difference approximation. For both versions, the half-edge normal fluxes related to a cell corner are expressed in terms of the half-edge temperatures and the cell-centered temperature. In both cases, the final scheme results from the elimination of the half-edge temperatures in terms of the cell-centered temperatures by enforcing the continuity of the half-edge normal fluxes across the edge impinging at a given node. This elimination procedure amounts to solve a local linear system at each node. The main differences between these two versions are the following:

• CCLADS scheme in its semi-discrete version satisfies a L 2 -stability property. It is also characterized by a global diffusion matrix which is positive semi-definite. Moreover, this scheme preserves linear field uniquely for triangular grids. Its accuracy is almost second-order on smooth distorted unstructured grids.

• CCLADNS scheme does not share the good stability properties of CCLADS. However, the flux discretization has been derived in a such a fashion that this scheme preserves linear fields regardless the shape of the cell. Namely, it shows a second-order rate of convergence on general unstructured grids.

Figure 2 :

 2 Figure 2: Notation related to polygonal cell ωc and one of its sub-cell ωpc.

Figure 3 :

 3 Figure 3: Continuity conditions for the half-edges fluxes and temperatures at a half-edge shared by two sub-cells attached to the same point. Labels c and d denote the indices of two neighboring cells. Labels p and r denote the indices of the same point relatively to the local numbering of points in cell c and d. The neighboring sub-cells are denoted by ωpc and ω rd . The half-edge fluxes, q ±pc , q ± rd and temperatures, T ± cp , T ± rd are displayed using blue color.

Figure 4 :

 4 Figure 4: Fragment of a polygonal cell ωc. Notation for the sub-cell ωpc: The half-edge fluxes, q ±pc , and temperatures, T ± pc are displayed using blue color.

  to obtain a discretization of the heat flux. More precisely, it leads to a | P(c) | × | P(c) | linear system satisfied by the edge-based normal components of the heat flux. This results in a non explicit expression of the edge-based normal components of the flux with respect to the edge-based temperatures and the mean cell temperature, which leads to a finite volume discretization characterized by edge-based and cell-based unknowns. In contrast to this approach, the sub-cell based variational formulation (31) yields a finite-volume discretization with one unknown per cell.

Figure 5 :

 5 Figure 5: Notation for a triangular cell. Half-edge degrees of freedom are displayed in blue color.

Figure 6 :

 6 Figure 6: Notation for sub-cells surrounding point p.

Here, L

  is the diagonal matrix defined by L c,d = l c δ c,d , where δ c,d denotes the Kronecker symbol, i.e., δ c,d = 1 if c = d and δ c,d = 0 if c = d. Matrix N is also a tridiagonal cyclic matrix and its non-zero components read as

Figure 7 :

 7 Figure 7: Limit angle θ lim c versus sub-cell aspect ratio ζ = lc l c+1 . The permitted angular domain to satisfy the geometric condition (77) is the region located between the red and the blue curves.

Figure 8 :

 8 Figure 8: Stencil of the finite volume scheme for a given cell ωc in a fragment of an unstructured grid. degrees of freedom are displayed using red squares. The blue arrow shows the diffusion flux exchange through point p between cell c and d.

Figure 10 :

 10 Figure 10: Notation related to cylindrical geometry.

  Temperature contours of the numerical solution.

Figure 11 :

 11 Figure 11: Anisotropic linear problem with discontinuous conductivity tensor on triangular grid.

Figure 15 :

 15 Figure 15: Triangular grids.

Figure 16 :

 16 Figure 16: Anisotropic linear problem with a non-uniform symmetric positive definite conductivity tensor. Numerical results obtained using the triangular grid displayed in Figure 15(a).

Figure 18 :

 18 Figure 18: Polygonal grids.

  [START_REF] De Groot | Non-equilibrium thermodynamics[END_REF] 7.35D-3 1.51 1.19D-2 1.7 5.08D-3 1.75 5.00D-2 7.48D-3 1.6 2.58D-3 1.82 3.66D-3 1.9 1.51D-3 1

- 1 ,

 1 45 4.D-3 0.96 6.04D-3 1.59 2.28D-3 1.81 5.00D-2 8.34D-3 0.65 2.06D-3 0.58 2.01D-3 1.59 6.49D-4 1.95 2in plasma physics in the presence of a magnetic field, refer to Section 2.2. Setting r = 0, the one-dimensional analytical steady solution of (110) corresponding to the Dirichlet boundary conditions, T (0) = 0 and T (1) = 1, writes as T (x) = exp( b a x) -1 exp( b a ) ∀x ∈ [0, 1].

6 .

 6 91 4.34D-2 1.80 2.12D-2 1Isotropic non-linear problemIn this section, we investigate the numerical solution of the following non-linear heat conduction equationρC v ∂T ∂t -∇ • (κ(T )∇T ) = 0,T (x, 0) = T 0 (x).

  (a). Note that this numerical solution preserves perfectly the cylindrical symmetry. Next, we construct a distorted grid from the polar 50 × 40 distorted polar grid.

Figure 22 := r 1 +

 221 Figure 22: Polar grids for the isotropic non-linear problem.

Figure 24 :

 24 Figure 24: Anisotropic non-linear problem: contours of z-component of the prescribed magnetic field.

  Anisotropic conductivity tensor without Righi-Leduc term, i.e., κ∧ = 0.

  Anisotropic conductivity tensor with Righi-Leduc term.

Figure 25 :

 25 Figure 25: Anisotropic non-linear problem: temperature contours at time T = 25 10 -12 s.

  Figure 1: Normalized Braginskii transport coefficients κ ⊥ and κ∧ versus normalized parameter Ωeτei.

	1				κ ⊥ /κ e κ ∧ /κ e
	0.8
	0.6
	0.4
	0.2
	0
	0	2	4	6	8	10
			Ω e τ ei

Table 1 :

 1 Anisotropic linear problem with discontinuous conductivity tensor: asymptotic errors in both maximum and l 2 norms and corresponding truncation error orders for quadrangular and polygonal grids.

	(a) Smooth grids.
	h 1.00D-1 7.86D-3 1.83 2.38D-3 1.72 E h max q h max E h l 2 q h l 2
	5.00D-2 2.2D-3 1.53 7.2D-4 1.85
	2.50D-2 7.6D-4	-	2.D-4	-
	(b) Random grids.

Table 2 :

 2 Anisotropic linear problem with a non-uniform symmetric positive definite conductivity tensor: asymptotic errors in both maximum and l 2 norms and corresponding truncation error orders for triangular grids.

	2 2.73 2.61D-2 2.59	1.01	0.44	0.28	1.51
	3.18D-2 1.55D-2 1.31 4.48D-3 2.26	0.75 -0.63	0.1	0.79
	1.58D-2 6.19D-3	-	9.2D-4	-	1.17	-	5.74D-2	-

Table 3 :

 3 Anisotropic linear problem with a non-uniform symmetric positive definite conductivity tensor: asymptotic errors in both maximum and l 2 norms and corresponding truncation error orders for quadrangular grids.

	(a) Rectangular grids.

Table 6 :

 6 Isotropic linear problem in cylindrical geometry: asymptotic errors in both maximum and l 2 norms and corresponding truncation error orders for polar grids.(a) Regular polar grids.

			CCLADS			CCLADNS
	h 7.77D-2	E h max 3.D-2	q max 1.83 3.68D-2 1.82 3.85D-2 1.73 2.06D-2 1.71 E h l 2 q l 2 E h max q max E h q l 2 l 2
	3.88D-2	8.4D-3 1.94 1.04D-2 1.94 1.15D-2 1.93 6.27D-3 1.92
	1.94D-2 2.18D-3	-	2.7D-3	-	3.01D-3	-	1.66D-3	-
				(b) Random distorted polar grids.

For quadrangular grids, Cp = 4.

This is precisely the case when the physical conductivity tensor is symmetric positive definite.

Table4: Anisotropic linear problem with a non-uniform symmetric positive definite conductivity tensor: asymptotic errors in both maximum and l 2 norms and corresponding truncation error orders for polygonal grids.

The two-point flux approximation corresponds to a finite difference approximation of the normal flux at an edge,

We have assessed the accuracy and the robustness of both schemes on various representative numerical test cases.

In future, we intend to extend CCLAD scheme capability to cope with two-dimensional nonconformal grids. We also plan to investigate the extension of CCLAD scheme to three-dimensional unstructured grids.

The one-dimensional solution, i.e., T = T (x) which corresponds to Dirichlet boundary conditions: T (0) = 0 and T (1) = 1, writes as

This is a linear continuous solution for which the heat flux q = -K∇T writes as = -

The normal component of the heat flux is continuous at the interface = 1 2 whereas its tangential component undergone a jump discontinuity since in general K yx l K xx r = K xx l K yx r . The boundary conditions applied on the top and the bottom boundaries of the computational domain are Dirichlet boundary conditions deduced from the analytical solution. For the numerical applications we have defined the entries of the conductivity tensor as

This test case is characterized by a linear solution. Since CCLADNS scheme by construction preserves linear solution, we uniquely assess the ability of CCLADS scheme to preserve the above linear solution. Using this scheme, we run the present test problem on a triangular grid and on various quadrangular grids.

Triangular grid

We compute the steady numerical solution using a triangular grid which is made of 246 cells. The grid is displayed in Figure 11(a). Note that the unstructured grid has been constructed such that the interface x = 1 2 coincides with the cell interfaces. The temperature contours of the numerical solution are plotted in Figure 11(b). These are vertical straight lines which match perfectly with the analytical solution. In addition, we observe that the obtained asymptotic errors are equal to zero up to machine precision. As expected, CCLADS preserves linear solutions on triangular grids.

Quadrangular grids

The numerical solution computed using a 10 × 10 Cartesian grid also matches perfectly with the analytical solution and the temperature isolines plot is identical to the one obtained previously, refer to Figure 11(b). This result is coherent with the fact that the flux approximation used in our finite volume scheme preserves linear solutions on cells which are parallelograms. Next, we study the accuracy of our scheme for a sequence of distorted grids which result from an analytical normal fluxes is defined according to (54b). The polygonal grids have been obtained by means of a Voronoi tesselation using an incremental algorithm which is described in [START_REF] Yu | A generalization of the notion of Dirichlet cell for nonconvex domain[END_REF]. More precisely, these polygonal grids have been derived applying the Lloyd's algorithm [START_REF] Du | Centroidal voronoi tesselations: applications and algorithms[END_REF] in order to achieve a better Doing so, we obtain polygonal grids which are almost centroidal Voronoi tesselations, i.e., tesselations wherein the generator attached to a cell coincides with its centroid. Note that we do not apply this algorithm to the cells located around the line x = 0.5 in order to preserve the interface. The convergence analysis for polygonal grids is performed computing the asymptotic errors and the corresponding orders of truncation error using formulas (111) and (112). The results displayed in Table 1(c) show that CCLADS scheme is converging and the convergence rate is only first-order in both l 2 norm and maximum norms.

Anisotropic linear problem with a non-uniform symmetric positive definite conductivity tensor

This test problem has been initially presented in [START_REF] Potier | Schéma volumes finis monotone pour des opérateurs de diffusion fortement anisotropes sur des maillages de triangles non structurés[END_REF]. Once more, it consists in finding the steady solution of (110). However, it is characterized by an anisotropic non-uniform conductivity tensor which writes for all (x, y)

where η is a positive parameter characterizing the level of anisotropy. This tensor is symmetric positive definite. Its eigenvalues are λ + = x 2 + y 2 and λ -= η(x 2 + y 2 ). Thus, its condition number which is equal to 1 η characterizes the anisotropy ratio. The source term, r, is computed such that the analytical solution (110) is given by T (x, y) = sin 2 (πx) sin 2 (πy).

We apply a homogeneous Dirichlet boundary condition on the boundaries of the computational domain, i.e., T (x, t) = 0, ∀x ∈ ∂D. For numerical applications, we choose η = 10 -2 . We assess the accuracy of CCLADS and CCLADNS schemes by running this test problem on sequence of triangular and distorted quadrangular grids. Next, we pursue our investigations using the sequence of the three smooth distorted grids plotted in Figure 12. The convergence analysis results obtained with these three grids are presented in Table 3(b). Once more, we observe a second-order convergence rate in l 2 norm for both schemes, whereas the convergence rate in maximum norm is almost second-order for CCLADS whereas it is second-order for CCLADNS.

Finally, we achieve the convergence analysis of the present problem by performing computations on a sequence of three random grids which are displayed in Figure 17. Note that these grids differ from the previous ones plotted in Figure 13 since this time all the nodes have been displaced according to (114). The results of the convergence analysis corresponding to this sequence of random grids are given in Table3(c). In comparison to the previous results, these ones show that the accuracy of CCLADNS scheme is superior to the accuracy of CCLADS which exhibits an order of convergence between first and second-order.

Polygonal grids

Concerning the polygonal grids we perform the convergence analysis on a sequence of three polygonal grids displayed in Figure 18. As previously, these polygonal grids have been obtained using centroidal Voronoi tesselation by means of the Llyod's algorithm. The numerical results are presented in Table 4. These data demonstrate that CCLADNS scheme exhibits a second-order rate of convergence on polygonal grids whereas CCLADS has an order of convergence in between first and second-order.

Anisotropic linear problem with a non-uniform and non-symmetric conductivity tensor

Here, we aim at assessing the accuracy of our finite volume scheme on an analytical problem characterized by a non-symmetric positive definite conductivity tensor which writes for all (x, y) ∈

Isotropic linear problem in cylindrical geometry

Here, we assess the accuracy of CCLADS and CCLADNS schemes in cylindrical geometry. To this end, let us denote by T = T (r) the analytical solution of the following steady heat equation with one-dimensional spherical symmetry

Here, r 1 > r 2 > 0, κ > 0 is the heat conductivity and α > 0. We determine the analytical solution by means of the change of variable r = exp(z) for r > 0. With this new variable, T = T (z) satisfies the ordinary linear differential equation

Its solution writes as T (z) = A exp(q + z) + B exp(q -z),

where q -and q + are given by

Returning to the initial variable, the solution writes

The integration constants A and B are computed by means of the boundary conditions (115b) For numerical applications we will use the following values of the parameters: r 1 = 0.1, r 2 = 1, φ = 10, κ = 2 and α = 1. This test case is run using both CCLADS and CCLADNS schemes on a sequence of regular and random distorted polar grids displayed respectively in Figure 20 and Figure 21. Using the aforementioned analytical solution, we perform the convergence analysis for both CCLADS and CCLADNS using regular and random distorted polar grids displayed in Figure 20 and Figure 21. The asymptotic errors and the trunction error orders are gathered in Table 6(a) for the regular grids and in Table 6(b) for the random distorted grids. Concerning regular grids, we observe that the order of convergence for both schemes is quite close to secondorder. For random distorted grids, CCLADS scheme is not converging, whereas CCLADNS is still converging with a truncation error order quite close to second-order. The numerical solutions on the distorted polar grid are computed using the five-point scheme (blue dots) and the CCLADS scheme (black dots).

presence of strong grid distortions. The correction of this flaw requires the use of finite volume schemes based on more accurate flux approximations. As it can be observed in Figure 23, our method, which is a nine-point scheme on quadrangular grids, brings us the expected improvement since it matches very well to the one-dimensional reference solution.

Anisotropic non-linear problem

The goal of this paragraph consists in assessing our finite volume scheme CCLADS against a test case which is representative of electron heat transport in a magnetized plasma [START_REF] Breil | Modelling of the magnetic field effects in hydrodynamic codes using a second order tensorial diffusion scheme[END_REF]. In this situation, the classical isotropic Spitzer-Härm conductivity has to be replaced by the anisotropic Braginskii tensor conductivity, refer to Section 2.2. We aim at solving the electron heat conduction equation in this particular context to assess the ability of our finite volume scheme to handle such a physical phenomenon. To this end, let us consider a two-dimensional plasma at rest and confined in the domain D = [0, δ] × [0, δ], where δ = 20 µm. The initial electron temperature of the plasma is defined for all x ∈ D as T 0 e (x) = 10 3 K if 0 ≤ x ≤ 18 µm, 3 10 7 K if 18 µm ≤ x ≤ 20 µm.

The density and the specific heat capacity of the plasma are given by ρ = 0.025 g/cm 3 and C v = 5 10 7 erg/K/g. Let us point out that these values correspond roughly to a layer of a hot wherein the gradient is approximated through the use of the temperatures of the two cells sharing that edge.