Aquib Qazi

Honoré Yin

Michel Sebès

Hugues Chollet

Cédric Pozzolini

A semi-analytical numerical method

Keywords: wheel-rail contact, Hertz, non-elliptic contact, vehicle dynamics, contact mechanics, boundary element method

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

Introduction

Safety, track fatigue analyses, and maintenance of railway vehicles are only some of the applications that highlight the importance of wheel-rail contact models. The complexity of different operations necessitates very different levels of modelling: a simplified and coarse model may often be sufficient. For hardware-in-the-loop simulations, where the calculations must be carried out in real time over railway tracks which are several kilometres long, very approximate wheel-rail contact models are generally used. Conversely, other applications require more sophisticated models, such as in the estimation of wear in urban rail networks, or the assessment of rolling contact fatigue. The shape of the contact area and the stress distribution within it categorise the different normal contact models available in the literature [START_REF] Kalker | Survey of wheel-rail rolling contact theory[END_REF][START_REF] Piotrowski | Wheel-rail contact models for vehicle system dynamics including multi-point contact[END_REF][START_REF] Knothe | History of wheel/rail contact mechanics: from Redtenbacher to Kalker[END_REF][START_REF] Meymand | A survey of wheel-rail contact models for rail vehicles[END_REF].

The case of normal contact between two continuous and non-conforming bodies was first treated by Hertz in 1882 [START_REF] Hertz | Über die Berührung fester elastischer Körper[END_REF], assuming frictionless contact and each solid body as an elastic half-space. The Hertzian solution is characterised by an elliptical contact CONTACT Aquib Qazi Email: aquib.qazi@univ-eiffel.fr patch, with a semi-ellipsoidal normal pressure distribution acting within the contact area. This approach is limited to profiles with constant curvature, and therefore may lead to approximate results in railway applications where the curvature of the profiles changes along the width of the rail section.

The use of discrete methods offers a more realistic approximation of the contact conditions in real-life applications. The finite element method (FEM) has already been used in various studies for elasto-plastic stress analyses and creep force characteristics, employing both implicit and explicit approaches [START_REF] Zhao | The solution of frictional wheel-rail rolling contact with a 3d transient finite element model: Validation and error analysis[END_REF][START_REF] Toumi | Finite element analysis of the frictional wheel-rail rolling contact using explicit and implicit methods[END_REF]. Certain recent works have also used Nitsche's method for frictional contact problems as opposed to the classical penalty methods or mixed methods, which offers good convergence and robustness [START_REF] Chouly | An overview of recent results on Nitsche's method for contact problems[END_REF]. Variational methods as the one described by Kalker's complete theory [START_REF] Kalker | Rolling contact phenomena: Linear elasticity[END_REF] and implemented in the program CONTACT [START_REF] Vollebregt | User guide for CONTACT, rolling and sliding contact with friction[END_REF], the boundary element method (BEM) like the one presented by Knothe and Le The [START_REF] Knothe | A contribution to the calculation of the contact stress distribution between two elastic bodies of revolution with non-elliptical contact area[END_REF], or the matrix inversion method (MIM) by Johnson [START_REF] Johnson | Contact mechanics[END_REF] also enable a good estimation of the contact properties. Despite the advances in computational prowess in recent years, discrete methods continue to struggle in terms of their capability to be used online in vehicle dynamics simulations. A very fine discretisation is usually required with FEM approaches to satisfactorily represent the contact patch boundaries, while classical boundary element methods require iterative techniques over large potential contact grids to verify the stress constraints.

In most software applications, a compromise is reached between the level of accuracy that is deemed acceptable and the time it takes to solve the contact problem. This has given rise to an intermediate level of approximate methods based on the theory of virtual penetration of the contacting surfaces such as the Linder model [START_REF] Linder | Verschleiss von eisenbahnrädern mit unrundheiten[END_REF], the Kik-Piotrowski model [START_REF] Piotrowski | A simplified model of wheel/rail contact mechanics for non-Hertzian problems and its application in rail vehicle dynamic simulations[END_REF], the semi-Hertzian model (STRIPES) [START_REF] Ayasse | Determination of the wheel rail contact patch in semi-Hertzian conditions[END_REF], and the extended Kik-Piotrowski (EKP) model [START_REF] Liu | A non-Hertzian method for solving wheel-rail normal contact problem taking into account the effect of yaw[END_REF] amongst others. The contact area in these approaches is approximated using the area in which the surfaces would overlap if there is no deformation. This overlapping area is bigger than the actual contact patch, and a scaling factor is used to prescribe the virtual penetration such that the interpenetration zone approaches the real contact area. The technique used for determining the scaling parameter is what broadly differentiates these approaches [START_REF] Sichani | Comparison of non-elliptic contact models: Towards fast and accurate modelling of wheel-rail contact[END_REF]. There exists in parallel another family of models wherein the non-elliptical part of the contact area is approximated using a series of individual ellipses. This so-called multi-Hertzian method, proposed by Pascal and Sauvage [START_REF] Pascal | The available methods to calculate the wheel/rail forces in non Hertzian contact patches and rail damaging[END_REF], is further simplified by Ayasse et al. [START_REF] Ayasse | Paramètres caractéristiques du contact roue-rail[END_REF] to an analytical approach. Some of the questions raised regarding these models when compared against virtual penetration methods in [START_REF] Piotrowski | Wheel-rail contact models for vehicle system dynamics including multi-point contact[END_REF] have been addressed recently by Pascal [START_REF] Pascal | Solving conformal contacts using multi-Hertzian techniques[END_REF].

A proposed improvement to these models is through the introduction of an analytical approximation for the surface deformation in the method ANALYN [START_REF] Sichani | A novel method to model wheel-rail normal contact in vehicle dynamics simulation[END_REF]. Although the results are promising in comparison to other simplified methods, this method remains sensitive to the negative-curvature correction, a characteristic of all point-based contact models based on Hertz' theory. This makes it challenging to implement the method as a true black-box solution in dynamic vehicle simulation. It is generally noted that this approach gives a good approximation of the contact area when the pressure distribution estimated is also accurate. A better approximation of the pressure distribution using a more thorough approach in the above family of methods should lead to better estimation of the contact patch, which forms the basis of the new semi-analytical approach proposed in this article.

In the first part, a reduced formulation of the contact between two elastic half-spaces is used together with a Hertzian model to solve the normal contact, following a semi-analytical methodology. The novel method is then validated using the results from STRIPES [START_REF] Ayasse | Determination of the wheel rail contact patch in semi-Hertzian conditions[END_REF], ANALYN [START_REF] Sichani | A novel method to model wheel-rail normal contact in vehicle dynamics simulation[END_REF], and CONTACT [START_REF] Vollebregt | User guide for CONTACT, rolling and sliding contact with friction[END_REF] with theoretical as well as wheelrail profiles. A comparison of the required computational resources is also carried out using different numerical approaches, before the final conclusions.

Methodology

Theoretical background

The problem describing the stress field in a semi-infinite half-space subjected to a concentrated normal surface force was studied in detail by Boussinesq [START_REF] Boussinesq | Application des potentiels à l'étude de l'équilibre et du mouvement des solides élastiques[END_REF]. For the half-space consideration to be valid locally, the contact area dimensions must be significantly smaller as compared to the principal radii of the contacting surfaces. This is generally true when considering several common contact scenarios such as the tread contact between wheels and rails. The following assumptions are also made:

• The bodies are elastic, homogeneous and isotropic.

• The contacting surfaces have a continuous profile.

• The normal and tangential problems can be solved sequentially.

One possible solution is using the potential theory of Boussinesq [START_REF] Boussinesq | Application des potentiels à l'étude de l'équilibre et du mouvement des solides élastiques[END_REF], described exhaustively in the book by Johnson [START_REF] Johnson | Contact mechanics[END_REF]. The normal elastic displacement ū(x, y) is related to the normal pressure distribution p(x, y) by the integral equation [START_REF] Cesbron | Numerical and experimental study of multicontact on an elastic half-space[END_REF]

∀(x, y) ∈ Σ, ū(x, y) = Σc p(ξ, η)T (x, y; ξ, η)dξdη, (1)
where Σ is the half-space domain under consideration, Σ C is the zone of the contact, and T is Boussinesq's influence function defined as

∀(x, y; ξ, η) ∈ Σ 2 , T (x, y; ξ, η) = 1 πE 1 (x -ξ) 2 + (y -η) 2 . ( 2
)
E is the effective Young's modulus of elasticity for the two contacting bodies Ω 1 and Ω 2 such that

1 E = 1 -ν 2 1 E 1 + 1 -ν 2 2 E 2 , (3)
where E 1 and E 2 are the modulus of elasticity, while ν 1 and ν 2 are the Poisson's ratios for the two solids. With a defined rigid approach δ, the only geometrical input required to solve the contact problem is the separation between the two undeformed surfaces h(x, y). The normal surface deformation ū(x, y) = ū1 (x, y)+ ū2 (x, y) must be evaluated at each point within the contact zone. This contact problem can be summarised in Figure 1. Considering Σ C as the area of contact, the conditions of Signorini give the cases of separation and of contact:

∀(x, y) / ∈ Σ C , δ -ū(x, y) -h(x, y) < 0 and p(x, y) = 0, (4)
∀(x, y) ∈ Σ C , δ -ū(x, y) -h(x, y) = 0 and p(x, y) > 0.

(5)

Figure 1. Two elastic bodies in contact

The boundary conditions defined in Equations ( 4) and ( 5) ensure that the normal pressure can only act within the contact zone, vanishing at the boundary of the contact area. The two bodies are also prevented from interpenetrating. If the rigid approach δ is given, the problem to be solved can be completely described by Equations ( 1) and ( 5). If the total resultant force N is known in advance, the following equilibrium condition is added:

∀(x, y) ∈ Σ C , N = Σc p(ξ, η)dξdη. (6)

Towards a reduced approach

The classical direct method to solve the contact problem is the matrix inversion method (MIM) found in [START_REF] Johnson | Contact mechanics[END_REF]. The potential contact area in the boundary element problem is overestimated in both the lateral and the rolling directions and divided into a rectangular grid, wherein the contact constraints of Equations ( 4) and ( 5) are then evaluated. The elements in which the pressure distribution has negative values are discarded, and the procedure is repeated until all the remaining elements satisfy the boundary conditions. The elastic deformation is thus calculated at each iteration for each element within the potential contact area, resulting in a large expenditure of computational resources. In the proposed new approach, the pressure distribution is assumed to be symmetric and elliptic about the x = 0 plane and the potential area of contact is discretised only in the lateral y direction. A similar strategy for discretisation is also used by Reusner for the treatment of roller bearings in [START_REF] Reusner | Druckflächenbelastung und oberflächen-verschiebung im wälzkontakt von rotationskörpern[END_REF], and by Knothe and Le-The for more arbitrary elastic bodies in [START_REF] Knothe | A contribution to the calculation of the contact stress distribution between two elastic bodies of revolution with non-elliptical contact area[END_REF]. The novel approach in the present paper is to consider the contact patch boundaries as a quasi-known quantity dependent on the form of the normal stress distribution. The potential contact area is first divided into thin strips with the larger dimension in the rolling direction, as shown in Figure 2, from where the name MIM-1D is chosen for the new method. The unknown in this case reduces to the maximum pressure values p 0i at the centre of each strip i. At each iteration, the half-length of the contact strip a i in the rolling direction is then computed as a function of the maximum pressure distribution p 0i , using an approximate analytical formulation based on Hertz' theory.

Discrete problem

The half-length and the half-width of each strip element are given as a i and b i respectively. The pressure distribution over each strip is assumed to be semi-elliptical in the x direction, and constant in the y direction. Thus, the expression for the pressure distribution over each strip may be written as

p(x, y i ) = p 0i 1 - x a i 2 , (7)
where p 0i is the maximum pressure at the centre (x i , y i ) of the strip i. The deformation and the separation at the centre are denoted as u i and h i respectively. Equations ( 1) and ( 5) in discrete form are:

∀i ∈ [1, n], u i = n j=1 C ij p 0j , ( 8
) ∀i ∈ [1, n] (x i , y i ) ∈ Σ C , δ -u i -h i = 0 and p 0i > 0, (9)
where n denotes the total number of elements in the potential contact area. The terms C ij from Equation ( 8) are called the coefficients of influence, and are defined as

C ij = 2 πE yj+bj yj-bj aj 0 1 -( ξ aj ) 2 (x i -ξ) 2 + (y i -η) 2 dξdη. (10)
The expression for C ij describes the influence of normal stress distribution in the jth element to induce elastic displacement in the ith element. The x i terms are ignored as they always remain zero in the case of MIM-1D. Equation ( 10) is written as:

C ij = 2 πE   yj+bj yj-bj aj 0 1 -( ξ aj ) 2 -1 ξ 2 + (y i -η) 2 dξdη + yj+bj yj-bj aj 0 dξdη ξ 2 + (y i -η) 2   . (11)
The first integral is regular and can be evaluated numerically using Gaussian quadrature. A higher number of integration points may be chosen closer to the diagonal terms. The second integral is singular when the denominator approaches zero. This expression represents the case of uniform normal pressure acting on a rectangular area of 2a i × 2b i , and an analytical solution of this problem has been presented in detail by Love [START_REF] Johnson | Contact mechanics[END_REF][START_REF] Love | The stress produced in a semi-infinite solid by pressure on part of the boundary[END_REF]. We have

2 yj+bj yj-bj aj 0 dξdη ξ 2 + (y i -η) 2 = (y + b j ) log   a j + (y + b j ) 2 + a 2 j -a j + (y + b j ) 2 + a 2 j   + (y -b j ) log   -a j + (y -b j ) 2 + a 2 j a j + (y -b j ) 2 + a 2 j   + 2a j log   (y + b j ) + (y + b j ) 2 + a 2 j (y -b j ) + (y -b j ) 2 + a 2 j   , (12)
where

y = y i -y j .
An additional advantage of treating the original integral expression in such a manner is the introduction of an analytical solution into the numerical results, which may consequently reduce the integration error linked to the quadrature method. This also permits the use of lesser number of integration points, and consequently faster calculation times. Equation ( 8) in matrix form is written as

u = C p, (13)
where the vectors p and u are given as {p 01 , ..., p 0n } T and {δ -h 1 , ..., δ -h n } T respectively. The matrix of influence coefficients C is

C =    C 11 . . . C 1n . . . . . . . . . C n1 . . . C nn    . (14)
This is the method of resolution generally followed in the case where the rigid approach δ is known in advance. Unlike in the classical method, C needs to be evaluated at the beginning of each iteration as the size of the elements in the x direction does not remain the same.

In the case where the normal force is prescribed instead of rigid body approach δ, the normal contact problem is solved with an additional iteration for δ. The resultant force at the end of each iteration is calculated using Equation ( 6):

N = π n i=1 (a i b i )p 0i . (15)
The initial value of δ can be taken as the Hertzian rigid body approach. The subsequent values of δ may be evaluated using a dichotomy or an iterative scheme based on Hertz' relations [START_REF] Knothe | A contribution to the calculation of the contact stress distribution between two elastic bodies of revolution with non-elliptical contact area[END_REF]:

δ (m+1) = δ (m) N Ñ (m) 2 3 , (16)
where Ñ (m) represents the resultant normal force at the end of mth iteration.

For numerical simulations, it is also possible to modify Equation ( 13) using the normal force over each strip N i (see Appendix A).

Calculation of the contact patch half-length a i

The half-length of the contact patch a i is evaluated analytically, using the expressions from Hertz' theory locally. a i is updated at each iteration to define the new potential contact zone, which is then used to construct the matrix C. The equation for the Hertzian contact ellipse is

x a 0 2 + y b 0 2 = 1, ( 17
)
where a 0 and b 0 are the semi-axes of the ellipse [START_REF] Ayasse | Wheel-rail contact[END_REF]:

a 0 = m 3 2 N 1 2E 1 (A + B) 1 3 , (18)
b 0 = n 3 2 N 1 2E 1 (A + B) 1 3 . ( 19
)
Here, A and B are the relative curvatures in the longitudinal and lateral directions respectively, while m and n are the non-dimensional Hertzian coefficients assessed using the curvatures (see Appendix in [START_REF] Ayasse | Determination of the wheel rail contact patch in semi-Hertzian conditions[END_REF], or [START_REF] Ayasse | Wheel-rail contact[END_REF]). The normal pressure distribution p(x, y) over the contact area is semi-ellipsoidal, and given by

p(x, y) = p 0H 1 - x a 0 2 - y b 0 2 , (20)
with the maximum normal pressure at the centre of ellipse,

p 0H = 3 2 N πa 0 b 0 . (21)
Equation ( 17) over each strip can be written as

a i a 0 2 + y b 0 2 = 1. ( 22
)
Combining Equations ( 20) and ( 22), the normal pressure distribution in the longitudinal x direction for the ith strip is

p(x, y i ) = p 0H a i a 0 1 - x a i 2 . ( 23
)
Comparing Equations ( 7) and ( 23), it is possible to deduce

p 0i = p 0H a i a 0 = 3 2 N πa 0 b 0 a i a 0 . ( 24
)
Using the expressions for the semi-axes a 0 and b 0 defined previously, the contact patch half-length can be found using

a i = π 2E * m 2 n A + B p 0i . (25)
This Hertzian expression for the contact patch boundary depends on the geometric properties of the profiles in contact, and the normal pressure at the centre of the strip under consideration. In non-elliptic cases, the curvatures and the Hertzian coefficients are replaced by their local values. Thus,

a i = π 2E * m 2 i n i A i + B i p 0i . (26)
The same expression for a i can also be obtained using the theory of ANALYN (see Appendix B). Hertz' solution remains valid only for positive values of B i . If the lateral curvature is negative at a given point, a correction must be carried out. Moreover, if the curvature is discontinuous, a smoothing is also applied. This correction and smoothing is done using the procedure described in [START_REF] Ayasse | Determination of the wheel rail contact patch in semi-Hertzian conditions[END_REF].

Iterative resolution

Equation ( 26) is incorporated into the iterative algorithm presented in Figure 3 to solve the contact problem described by Equation [START_REF] Linder | Verschleiss von eisenbahnrädern mit unrundheiten[END_REF]. In order to speed up the computation, the matrix of influence coefficients C is constructed only using the elements i where the separation h i is less than a predefined maximum value h max . When the resultant normal force N is known, the same algorithm is repeated for each value of δ (m) , evaluated using Equation ( 16). The total normal force at the mth iteration Ñ (m) is computed using Equation [START_REF] Ayasse | Determination of the wheel rail contact patch in semi-Hertzian conditions[END_REF]. The algorithm in this case converges when a user defined tolerance value = | Ñ (m) -N | is attained.

Results and discussion

The proposed new approach MIM-1D is implemented as a Matlab function. The approximate surface deformation method ANALYN [START_REF] Sichani | A novel method to model wheel-rail normal contact in vehicle dynamics simulation[END_REF], and the virtual penetration method STRIPES [START_REF] Ayasse | Determination of the wheel rail contact patch in semi-Hertzian conditions[END_REF] are also programmed using Matlab. The results from the commercial version of the program CONTACT (v20.1) [START_REF] Vollebregt | User guide for CONTACT, rolling and sliding contact with friction[END_REF] are used as the reference. It should be noted that the ANALYN results are sensitive to the method used for the negative-curvature correction. In the original publication, the negative lateral curvature values are replaced using a fifth-degree polynomial. However, this correction strategy has not been explained in further detail, and here this is done heuristically to obtain results as close as possible to the ones presented in [START_REF] Sichani | A novel method to model wheel-rail normal contact in vehicle dynamics simulation[END_REF].

Theoretical profiles

Hertzian

To validate the new approach, a Hertzian case is considered with the contact between a sphere (R = 40 mm) and a flat surface. The separation curve between the two surfaces The rigid approach between the two bodies is taken as 1 mm. The contact patches obtained using different approaches and the maximum pressure distribution p 0 (y) are shown in Figure 5. All the methods can be observed to be in good agreement with each other. The relative error in the contact area for MIM-1D is found to be within 1% of Hertz's analytical solution, which can be attributed to the accuracy of numerical procedure used.

Non-Hertzian

The contact between a flat surface, and a body of revolution (rolling radius R n = 400 mm) with a non-Hertzian theoretical profile developed using two different radii R 1 = 40 mm and R 2 = 500 mm on either side of the point of first contact is considered next. The separation curve between the two surfaces is shown in Figure 4b. The material properties are the same as in the Hertzian case, with the bodies pushed 1 mm towards each other.

The results for the contact area and the maximum pressure distribution are presented in Figure 6. From Figure 6a, it can be seen that MIM-1D and ANALYN correspond reasonably well with the reference results from CONTACT. The relative error in the contact area for MIM-1D is within 1% of the reference results. STRIPES notably underestimates the width of the contact patch: this is expected, as neglecting the surface deformation should lead to a smaller contact zone. Although the pressure distribution curve for ANALYN in Figure 6b follows the same trend as that of the reference, the peak of the pressure curve remains significantly higher. The results using MIM-1D can be observed to be in a better agreement with the reference method.

Wheel-rail contact

The case of wheel-rail contact is presented using the standard wheel profile S1002 over the rail profile UIC60, with an inclination of 1:40. The material for both the wheel and the rail is steel. The resultant normal contact force is taken as 78500 N. The results are presented for various positions of the wheel, displaced from its centre position over the rail, denoted by ∆y. The sign convention is taken the same as in [START_REF] Sichani | A novel method to model wheel-rail normal contact in vehicle dynamics simulation[END_REF], where a positive ∆y signifies the movement of the contact point towards the wheel flange. The discretisation size is taken as approximately 0.2 mm for all the considered cases.

The test cases are chosen to remain in the tread region with low contact angles, where the half-space assumption is not violated. From Figure 7, it can be seen the MIM- 1D results match closely with the reference CONTACT results in all of the presented cases. This is highlighted more prominently in the maximum pressure distribution over the contact patch length. It can also be remarked that MIM-1D manages to accurately capture the characteristic slight variations in the pressure distribution, such as those presented in the tail end of the case ∆y = -1, a trait missing in the other simplified methods. ANALYN remains more precise compared to the virtual penetration methods, as neglecting the surface deformation in STRIPES again leads to an underestimated contact zone. The negative curvature compensation procedure used in STRIPES may have an effect as well [START_REF] Quost | Assessment of a semi-Hertzian method for determination of wheel-rail contact patch[END_REF]. This sensibility of STRIPES related to the processing of the curvature is found to be true in the case of ANALYN as well, with some fine tuning required to obtain the desired results. A correction strategy dependent on the applied contact force and the separation may possibly improve the results [START_REF] Sichani | A novel method to model wheel-rail normal contact in vehicle dynamics simulation[END_REF]. Figure 8 presents the relative comparison of the contact area for a range of ∆y values to the CONTACT results, emphasising the improvement using the new method as compared to existing fast approaches.

Computational cost

To make a representative comparison of the computational cost, MIM-1D is tested for different mesh sizes using the Hertzian profile against the other complete numerical methods i.e. Kalker's CONTACT, and the classical matrix inversion method (MIM) implemented as Matlab function. These comparisons are made using a 64 bits 2.70 GHz Intel processor. Only the computation time for the normal contact problem is considered. The time for CONTACT is taken from the generated output file, while the other two methods are measured using the elapsed CPU time averaged over a finite number of runs. These results are tabulated in Table 1. The gain in the computation speed using the new method with respect to the classical MIM is clear, with up to a 1000-fold speed up in cases where a comparison is possible. The conventional desktop employed is in fact unable to supply enough memory to carry out the simulations when using the fine meshes with MIM. CONTACT uses a boundconstrained conjugate gradient (BCCG) method with a fast Fourier transform (FFT) pre-conditioner which permits quicker resolution, even with a very fine discretisation [START_REF] Vollebregt | 100-fold speed-up of the normal contact problem and other recent developments in "CONTACT[END_REF]. The MIM-1D implementation here simply employs Matlab's inbuilt direct solver yet permits a significant gain in runtime, in no small part due to the reduced semianalytical formulation of the method.

Conclusion

A simplified boundary element formulation is presented in this paper and tested against two existing approximate methods, and a complete numerical method used as the reference. In essence, this approach further develops the discretisation strategy presented in [START_REF] Knothe | A contribution to the calculation of the contact stress distribution between two elastic bodies of revolution with non-elliptical contact area[END_REF], by using a semi-analytical methodology to determine the contact patch dimension in the rolling direction. The novel method, implemented in the algorithm MIM-1D, provides a precise approach comparable to more rigorous complete methods such as CONTACT, with lesser computing effort. The results from the theoretical and wheel-rail test cases presented allow the following conclusions to be drawn:

• The proposed approach MIM-1D enables an improved approximation of the pressure distribution and the contact area as compared to the other existing virtual penetration [START_REF] Ayasse | Determination of the wheel rail contact patch in semi-Hertzian conditions[END_REF] and approximate surface deformation [START_REF] Sichani | A novel method to model wheel-rail normal contact in vehicle dynamics simulation[END_REF] methods. • Comparison with a complete numerical method for wheel-rail contact shows close agreement, with the relative error in the contact area as compared to the reference results being less than 3% in all the tested cases using theoretical and wheel-rail profiles. • The reduction in the number of system unknowns as compared to CONTACT due to a semi-analytical approach provides up to a 10-fold speed up with the current Matlab implementation. This should ideally improve further by using a programming language closer to machine language such as Fortran. Similarly, optimisation with respect to the iterative solver used should make the algorithm run faster.

It is important to keep in context that the latest versions of CONTACT incorporate advanced numerical techniques to enable faster solving of the normal contact problem [START_REF] Vollebregt | 100-fold speed-up of the normal contact problem and other recent developments in "CONTACT[END_REF]. Incorporating similar numerical optimisation strategies may permit further improvement in the performance of MIM-1D as well. Further developments must also be implemented in order to accurately treat the flange contact, where the contact angle varies significantly in a small zone, thus violating the half-space assumption. These improvements are currently under study, within the scope of the implementation of MIM-1D in the dynamic multi-body simulation code VOCO (VOitures en COurbe). The eventual implementation should account for the coupling with the tangential contact problem, while taking the effect of friction into account.

With proper optimisation, solving the classical normal contact problem using only strip elements should enable MIM-1D to be used as a good reference for other coarse models commonly employed in dynamic vehicle simulations.

Combining Equations (B5) and (B6), the half-length of the contact patch a i can finally be written as

a i = π 2E * m 2 i n i A i + B i p 0i . (B7)

Figure 2 .

 2 Figure2. Contact area Σ C divided into strips, and the normal pressure distribution[START_REF] Knothe | A contribution to the calculation of the contact stress distribution between two elastic bodies of revolution with non-elliptical contact area[END_REF]

Figure 3 .

 3 Figure 3. MIM-1D algorithm with a given rigid body approach δ

Figure 4 .Figure 5 .

 45 Figure 4. Separation profiles of the two contacting surfaces

Figure 6 .

 6 Figure 6. Results for a non-Hertzian profile using existing methods and MIM-1D

Figure 7 .

 7 Figure 7. Contact patch (left), and the maximum pressure distribution p 0 (y) (right) for wheel-rail contact cases (from top to bottom): (a) ∆y = -1 mm, (b) ∆y = 0 mm, (c) ∆y = 1 mm, (d) ∆y = 2 mm, and (e) ∆y = 5 mm.

Figure 8 .

 8 Figure 8. Comparison of the relative error in the contact area with reference results from CONTACT

Table 1 .

 1 Comparison of CPU time with existing numerical approaches

		n	Time (s)
	Mesh size (mm) Number of elements a	MIM	CONTACT MIM-1D
	0.2	51	1.002	0.1	0.010
	0.1	101	40.28	0.2	0.037
	0.02	501		5.2	0.342
	0.01	1001		14.2	1.111

a The total number of elements in the case of MIM and CONTACT is n × n.

Funding

The authors would like to thank Association Nationale de la Recherche et de la Technologie (ANRT), and ESI Group (France) for their financial support.

Appendix A. Solving the contact problem for the normal force N i

The normal force acting on each strip N i can be evaluated using Equation [START_REF] Ayasse | Determination of the wheel rail contact patch in semi-Hertzian conditions[END_REF]:

(A1)

The system of equations in matrix form is

where N = {N 1 , ..., N n } T , and the matrix of influence coefficients C is given by

The half-length of the contact patch at the end of each iteration can be computed using

Equations ( A2) and (A4) can subsequently be employed in the algorithm presented in Figure 3 to solve the contact problem.

Appendix B. Contact patch half-length a i using ANALYN

The expression for the contact patch boundaries in ANALYN [START_REF] Sichani | A novel method to model wheel-rail normal contact in vehicle dynamics simulation[END_REF] is

where

Here, d i is the penetration and the term β i h i takes the surface deformation into account analytically as opposed to methods based on the virtual penetration where the deformation is neglected. The coefficients α i and β i are defined as:

where A i and B i are the relative longitudinal and lateral curvatures respectively, while m i , n i , and r i are non-dimensional Hertzian coefficients calculated using the local curvatures. The maximum pressure value p 0i is

Squaring Equation (B1) and taking into account the expression for α i , we have


Pierre Anguill
email: pierre.anguill@cea.fr

Patricia Cargo
email: patricia.cargo@cea.fr

Cedric Énaux
email: cedric.enaux@cea.fr

Philippe Hoch
email: philippe.hoch@cea.fr

Emmanuel Labourasse
email: emmanuel.labourasse@cea.fr

Gerald Samba
email: gerald.samba@cea.fr

An asymptotic preserving method for the linear transport equation on general meshes

Keywords: Finite Volume, radiative transfer, asymptotic analysis, computational transport, monotone anistotropic diffusion, unstructured meshes

While many numerical methods for the linear transport equation are available in the literature in 1D or on Cartesian meshes, fewer works are dedicated to the resolution of this model on unstructured meshes. In the context of radiative hydrodynamics, we need a method capable to handle a wide range of radiation regimes going from freestreaming to diffusion and to be coupled with a Lagrangian hydrodynamics solver. In this paper we design a method based on the micro-macro paradigm and to the Discrete Ordinates (S N ) angular discretization, which fulfills these requirements. It allows to choose the limit transport scheme and the limit diffusion scheme. It is compared on challenging test problems to a Discontinuous Finite Element (DFE) method.

Introduction

In this paper, we are interested in the transport of photons in astrophysical phenomena or in Inertial Confinement Fusion experiment (refer to [START_REF] Castor | Radiation Hydrodynamics[END_REF][START_REF] Mihalas | Foundations of radiation hydrodynamics[END_REF][START_REF] Zel'dovich | Physics of Shock Waves and High-Temperature Hydrodynamic Phenomena[END_REF] for an overview of radiation driven physics). To this end, we consider the grey radiative transfer equations [START_REF] Mihalas | Foundations of radiation hydrodynamics[END_REF]

1 c ∂ t I + ω • ∇I + σI = σB(T M ), ∂ t E M = cσ (E r -aT 4 M ) ,
where c is the speed of light, I denotes the grey radiative intensity, σ the cross section, B is the Planck function at the matter temperature T M , ω the direction, E M the matter energy density, and E r the radiative energy density. One of the problem encountered when dealing with these equations is that they are non-linear and coupled. One solution to deal with this difficulty is the linearization used by Fleck and Cummings in [START_REF] Fleck | An implicit Monte Carlo scheme for calculating time and frequency dependent nonlinear radiation transport[END_REF]. This brings us, at each time iteration n, to the following linear problem on the time interval [t n , t n + ∆t]

1 c ∂ t I + ω • ∇I + σ n I = f n σ n B n + (1 -f n )σ n < I >, ∂ t E M = f n cσ n (E r -Φ n ) , (1)
where < • > accounts for the angular mean value, f = 1/ (1 + βcσ∆t) denotes the factor of Fleck, β = dΦ dE M , Φ = aT 4 M , and the superscript n denotes quantities taken at time t n .

The System (1) involves very different length and time scales. These scales are related to the scaled mean free path of the photons, called the Knudsen number. For Knudsen numbers close to unity, the transport of photons can be modeled by a hyperbolic kinetic equation. However, as the Knudsen number tends to zero (for instance in opaque media), it is well-known (refer for instance to [START_REF] Bardos | The nonaccretive radiative transfer equations: Existence of solutions and Rosseland approximation[END_REF]) that this equation tends to a diffusion equation (that is of parabolic type). That poses major difficulties for the numerical simulation, since a standard discretization constrains the time and space step proportionally to the Knudsen number. One way to circumvent this problem, is to decompose the calculation domain into a diffusive domain (small Knudsen number) and a transport domain (large Knudsen number), as explained for instance in [START_REF] Bal | Coupling of transport and diffusion models in linear transport theory[END_REF][START_REF] Degond | A smooth transition model between kinetic and hydrodynamic equations[END_REF][START_REF] Degond | A multiscale kinetic-fluid solver with dynamic localization of kinetic effects[END_REF][START_REF] Golse | A domain decomposition analysis for a two-scale linear transport problem[END_REF][START_REF] Hauck | A collision-based hybrid method for time-dependent, linear, kinetic transport equations[END_REF][START_REF] Klar | Asymptotic-induced domain decomposition methods for kinetic and drift diffusion semiconductor equations[END_REF][START_REF] Larsen | Asymptotic solution of neutron transport problems for small mean free paths[END_REF]. Unfortunately, this solution is complicated to apply in our case, because the boundaries between diffusion and transport are fuzzy and moving, due to ablation processes. The alternative is to design numerical schemes able to handle all the radiative regimes without loss of accuracy or too drastic time-step restriction. These type of scheme is called Asymptotic Preserving (AP), and has been the subject of a considerable amount of work since the seminal work of Larsen et al [START_REF] Larsen | Asymptotic solutions of numerical transport problems in optically thick, diffusive regimes II[END_REF][START_REF] Larsen | Asymptotic solutions of numerical transport problems in optically thick, diffusive regimes[END_REF] (refer also to Jin and Levermore [START_REF] Levermore | Numerical schemes for hyperbolic conservation laws with stiff relaxation terms[END_REF], who introduced the expression Asymptotic Preserving). A comprehensive description of the work that has since been done in this area is beyond the scope of this paper, and we invite Since we want to couple our radiation solver with a Lagrangian hydrodynamic scheme for the Euler system (as [START_REF] Caramana | The construction of compatible hydrodynamics algorithms utilizing conservation of total energy[END_REF][START_REF] Carré | A cell-centered Lagrangian hydrodynamics scheme on general unstructured meshes in arbitrary dimension[END_REF][START_REF] Maire | A cell-centered Lagrangian scheme for two-dimensional compressible flow problems[END_REF][START_REF] Von Neumann | A method for the numerical calculation of hydrodynamic shocks[END_REF]), we desire our solver to fulfill the following properties:

P1 to be consistent on general meshes (meaning any polytopal valid mesh), P2 to enforce conservation of radiative energy, P3 to be able to handle all the radiation regimes from free-streaming to diffusion (Asymptotic Preserving),

P4 to have the radiative energy degrees of freedom located at the centers of the elements to ensure the compatibility with the hydrodynamic scheme.

Moreover, one of the major drawback of the Discrete Finite Element method on general meshes, is their algorithmic cost in the diffusion regime. While they are indeed Asymptotic Preserving, in the sense that the accuracy of the result and the time-step are independent of the Knudsen number, the CPU consumption strongly depends on acceleration technics as DSA (Diffusion Synthetic Acceleration) [START_REF] Alcouffe | Diffusion synthetic acceleration methods for the diamond-differenced discrete-ordinates equations[END_REF] or Transport Synthetic Acceleration (TSA) [START_REF] Ramone | A Transport Synthetic Acceleration Method for Transport Iterations[END_REF]. The development of these methods is still a field of research. This is why we want to free ourselves from these algorithms, and even be able P5 to choose the diffusion scheme towards which our method tends when the Knudsen number becomes low.

To our knowledge, no existing scheme fulfills all these requirements. That is why we were interested in designing a new scheme. To achieve this, we were inspired by a very popular method in the AP community, called micro-macro formulation. This method was introduced in [START_REF] Klar | Numerical passage from radiative heat transfer to nonlinear diffusion models[END_REF] and popularized by Lemou and co-authors [START_REF] Bennoune | Uniformly stable numerical schemes for the Boltzmann equation preserving the compressible Navier-Stokes asymptotics[END_REF][START_REF] Lemou | A New Asymptotic Preserving Scheme Based on Micro-Macro Formulation for Linear Kinetic Equations in the Diffusion Limit[END_REF]. It consists in decomposing the unknown (the radiative intensity in our case) into an equilibrium part (the angular mean value) and a deviation. Since the radiation is isotropic (does not depend on angle) in the diffusion limit, the deviation part tends to zero with the Knudsen number. This principle has been widely used for designing numerical AP schemes during the past two decades ( [START_REF] Crestetto | Asymptotically complexity diminishing schemes (ACDS) for kinetic equations in the diffusive scaling[END_REF][START_REF] Crouseilles | An asymptotic preserving scheme based on a micro-macro decomposition for collisional Vlasov equations: diffusion and high-field scaling limits[END_REF][START_REF] Einkemmer | An asymptotic-preserving dynamical low-rank method for the multi-scale multi-dimensional linear transport equation[END_REF][START_REF] Shi | A micro-macro decomposition-based asymptoticpreserving scheme for the multispecies boltzmann equation[END_REF][START_REF] Laiu | A Positive Asymptotic-Preserving Scheme for Linear Kinetic Transport Equations[END_REF][START_REF] Roger | A hybrid transportdiffusion model for radiative transfer in absorbing and scattering media[END_REF] among others), and has similarities with some earlier methods called Variable Eddington Factor (described for instance in [START_REF] Mihalas | Foundations of radiation hydrodynamics[END_REF] p503). However, to our knowledge, no extension to general meshes has been achieved. It may be because the limit diffusion scheme arising from the discretization proposed in [START_REF] Lemou | A New Asymptotic Preserving Scheme Based on Micro-Macro Formulation for Linear Kinetic Equations in the Diffusion Limit[END_REF] is a two-point flux approximation, which is known to be inconsistent on most of the meshes (refer for instance to [START_REF] Eymard | TP or not TP, that is the question[END_REF]). In this paper, we propose a micro-macro scheme for the linear transport equation which can be proven to fulfill the propositions P1 to P5. To achieve this, we first perform the angular semi-discretization of the transport equation, using a Discrete Ordinate method. Then we perform the micro-macro decomposition, and obtain a transport equation for the angular dependant variables (deviation), and a diffusion equation for the angular main value (Section 2). The spatial and semi-implicit temporal discretization are described in Section 3. It is shown that our temporal discretization permits to choose any diffusion scheme for the main angular value. It allows us to choose a diffusion scheme which preserves the positivity of the energy e, which is an important property, especially in this regime (in which radiation is strongly coupled with the matter temperature). Finally, we assess the new method on several test problems in a wide range of regime, and compare it either to analytical solutions (if they exist) or to Discrete Finite Element method results in Section 5.

Problem and properties

Continuous transport equation

We focus on the following linear integro-differential transport equation, which is indeed of the form of Eq. ( 1)

1 v ∂ t u(t, x, ω)+ω • ∇u(t, x, ω)+σ t u(t, x, ω) = q(t, x, ω)+σ s S 2 u(t, x, ω ) dω 4π , (2)
where

• x ∈ Ω is the space coordinates in the open domain Ω ⊂ R d , with d the dimension of the problem;

• t ∈ R + is the time;

• u(t, x, ω) ∈ R the distribution function;

• v ∈ R + is the norm of the velocity of the particles (equal to the speed of light in the case of photons);

• ω ∈ S 2 accounts for the direction of the velocity of the particles, S 2 being the unit sphere;

• σ t = σ a +σ s , with σ t , σ a and σ s respectively the total, absorption, and scattering cross sections;

• q(t, x, ω) is a given source term.

We introduce the two first moments of u with respect to ω e =< u >=

S 2 u(ω ) dω 4π , f = S 2 ω u(ω ) dω 4π . (3)
By averaging Eq. ( 2) over all directions ω, we obtain a balance equation satisfied by e

1 v ∂ t e + ∇ • f + σ a e = s, (4)
with the isotropic part of the source

s(t, x) = S 2 q(ω ) dω 4π . (5)
It is shown in [START_REF] Case | Linear Transport Theory[END_REF] that if q > 0, σ a > 0 and σ s > 0, Eq. (2) satisfies positivity. Moreover, in [START_REF] Dautray | Analyse mathématique et calcul numérique pour les sciences et les techniques[END_REF], we have that if σ a = 0, and q = 0, u and e satisfies a conservation law and verify a global maximum principle (u min ≤ u ≤ u max and e min ≤ e ≤ e max , where u min , u max , e min , e max depend only on the initial state and the boundary conditions). It is well known (refer for instance to [START_REF] Castor | Radiation Hydrodynamics[END_REF][START_REF] Mihalas | Foundations of radiation hydrodynamics[END_REF]), that Eq. (2) satisfies a diffusion limit, in thick media, in the sense that when σ s → ∞, σ a → 0 and q → 0,

up to O( 1 σ 2 s ): u = e - 1 σ s ω • ∇e, f = - 1 3σ s ∇e and 1 v ∂ t e = ∇• 1 3σ s ∇e .

Semi-discrete angular equation

We introduce K > 0 angular directions, noted ω k . Each direction is weighted by a positive weight w k . Following the recommendations of [START_REF] Carlson | Transport theory: discrete ordinates quadrature over the unit sphere[END_REF][START_REF] Carlson | Computing Methods in Reactor Physics[END_REF], we impose on the couples (w k , ω k ) the following constraints:

                     K k=1 w k = 1, K k=1 w k ω k = 0, K k=1 w k ω k ⊗ ω k = 1 3 I K . ( 6
)
120

This leads to an even number of directions, opposite in pairs and of equal weights 1/K. This quadrature is called quadrature Equal Weights [START_REF] Carlson | Transport theory: discrete ordinates quadrature over the unit sphere[END_REF]. It is then sufficient to discretize the sphere uniformly and then apply a renormalization to all the discrete directions to restore the coefficient 1/3 of the moment of order 2.In this work, we have chosen to use a simple uniform angular grid which is controlled by its number 125 of degrees of freedom N ≥ 1, and yields a number of discrete ordinates K = 4N 2 (in 2D). Thus, with u k ≡ u(t, x, ω k ), we are able to write the K semi-discrete equations associated to the transport equation ( 2):

1 v ∂ t u k + ω k • ∇u k + σ t u k = q k + σ s K k =1 u k w k . (7)
It yields the semi-discrete formula for e, f and s:

e(t, x) ≈ K k=1 w k u k , f(t, x) ≈ K k=1 w k u k ω k and s(t, x) ≈ K k=1 w k q k .
Likewise, the semi-discrete diffusion limit comes down to:

∀t, when σ s → ∞, σ a → 0 and q → 0, and up to O 1

σ 2 s , u k = e -1 σs ω k • ∇e (∀k), f = -1 3σs ∇e, 1 v ∂ t e = ∇ • 1 3σs ∇e . 6

Semi-discrete micro-macro formulation

In order to have an Asymptotic-Preserving scheme, and to clearly distinguish between transport terms and diffusion terms, we introduce δu k ≡ u k -e, the deviation to isotropy.

By substracting (4) to [START_REF] Bal | Coupling of transport and diffusion models in linear transport theory[END_REF], we obtain:

1 v ∂ t δu k + ω k • ∇δu k + σ t δu k = ∇ • f -ω k • ∇e + δq k , (8)
135

where δq k ≡ q k -s. Thus, the first two angular moments of δu k are:

K k=1 w k δu k = 0, K k=1 w k δu k ω k = f . (9)
It leaves unchanged the balance law for e (4) except for the definition (9) of f as a function of the δu k . Since no further approximation has been made, we still have the semi-discrete diffusion limit, hence:

when σ s → ∞, σ a → 0 and q k → 0, and up to O 1

σ 2 s , δu k = -1 σs ω k • ∇e (∀k), f = -1 3σs ∇e, 1 v ∂ t e = ∇ • 1 3σs ∇e .

2D planar spatio-temporal discretization on unstructured mesh

We aim at solving the system (4)- [START_REF] Bardos | The nonaccretive radiative transfer equations: Existence of solutions and Rosseland approximation[END_REF] with the closure relations [START_REF] Bennoune | Uniformly stable numerical schemes for the Boltzmann equation preserving the compressible Navier-Stokes asymptotics[END_REF], on a 2D unstructured mesh. Since our ambition is to couple it with Lagrangian hydro schemes as [START_REF] Carré | A cell-centered Lagrangian hydrodynamics scheme on general unstructured meshes in arbitrary dimension[END_REF][START_REF] Maire | A cell-centered Lagrangian scheme for two-dimensional compressible flow problems[END_REF][START_REF] Von Neumann | A method for the numerical calculation of hydrodynamic shocks[END_REF] (even if it is out of the scope of this paper), the geometric quality of the mesh is not controlled. This is why our discretization choices are guided by the concern of robustness. The conservation property of e is very important for the targeted applications and has to be preserved. Moreover, we want our numerical method to be able to deal with all range of radiation regime, going from pure transport to diffusion. We give ourselves a regular mesh M of Ω. We consider here that M is defined by a finite collection of cells j that partition Ω(t). Specifically, we have

1.∀j ∈ M, j ⊂ Ω, 2.∀x ∈ Ω, ∃j ∈ M s.t. x ∈ j, 3.∀j 1 , j 2 ∈ M, j 1 ∩ j 2 = ∅ ⇐⇒ j 1 = j 2 , or ∂j 1 ∩ ∂j 2 = ∅ and j 1 ∩ j 2 = ∅.
The last relationship indicates that if two distinct cells intersect, it can only be through a piece of their edge (possibly a vertex). We define the following notations depicted on Fig. 1. Notations

• j : cell index, and by extension the cell,

• r : node index, and by extension the dual cell,

• n : time iteration index,

• p : optional sub-iteration index,

• V j : volume of cell j,

• V r : volume of dual cell r,

• x j = 1 V j j x : center of cell j,
• x r : position of node r,

• x j 1 j 2 : vector equals to x j 2 -x j 1 ,

• l : edge of the primal cell j,

• l : edge of the dual cell associated to a vertex r,

• r(j) : designates a node r of the cell j,

• j(r) : designates a cell j which owns the node r,

• r (r) : designates a node r which is adjacent to the node r,

• l(j) : designates an edge l which is contained by the cell j,

• C jr := ∇ xr V j : outgoing unit normal at node r with respect to cell j,

• C rj := ∇ x j V r : outgoing unit normal at the center of the cell j with respect to the dual cell r,

• n jl : outgoing normal at edge l with respect to the cell j,

• n r l : outgoing normal at edge l with respect to the dual cell r. We also need to define a dual mesh. It is built by joining the center of the cells j, and is depicted on Fig. 2.

r • j • l n jl C jr Figure 1: Notations (primal mesh) j • r • r • l n r l C rj
We use a finite volume method for the spatial discretization. We integrate the Eq. ( 4) on the primal mesh, which gives

1 v ∂ t j e + j ∇ • f + j σ a e = j s, (10)
We use the Green formula, a second-order accurate quadrature at nodes to approximate the boundary integrals, and define 1 V j j ϕ := ϕ j , which yields

V j 1 v ∂ t e j + r(j) f r • C jr + V j σ a,j e j = V j s j , (11)
With this quadrature formula, the flux f r has to be evaluated at the nodes. As we saw previously, the angular discretization of f r , writes:

f r = K k=1 w k δu k,r ω k . ( 12
)
This is why we choice to perform a staggered discretization, and to writes the discrete balance of δu k,r at the nodes. Eq. ( 8) is then approximated with:

1 v ∂ t δu k,r + [ω k • ∇δu k ] r + σ t,r δu k,r = K k =1 w k [ω k • ∇δu k ] r -ω k • [∇e] r + δq k,r . (13)
The notation [ϕ] r accounts for a spatial discretization of ϕ at the nodes, which is explained in the following.

Remark 3.1. Another possible choice of staggered discretization could be performed, in locating the δu k at the center of the edges instead at the nodes. In this case, the advection step would be performed on the diamond dual mesh, and the final global scheme would have the same properties than the one presented in this paper.

Temporal discretization

For this scheme, we use an implicit formulation. But, in order to avoid to solve a large linear system, we perform an hybrid time discretization, with sub-iterations. The discrete version of Eq. ( 13) with this time discretization is:

1 v δu n+1,p+1 k,r -δu n k,r ∆t + [ω k • ∇δu k ] n+1,p+1/2 r + σ n t,r δu n+1,p+1 k,r = K k =1 w k [ω k • ∇δu k ] n+1,p+1/2 r -ω k • [∇e] n+1,p+1 r + δq n k,r , (14)
where n is the time index and p the sub-iteration index. The index p + 1/2 means that a part of discretization is performed at the sub-iteration p, while the other is performed at the iteration p + 1 as explained in the following. The term [∇e] n+1,p+1 r is taken implicitly, which relies on an implicit diffusion equation on e, also explained in the following.

A very important property for this formulation to be meaningful is the following. It means that e indeed remains the discrete angular mean of u k with time.

The advection fluxes [ω k •∇δu k ] r are discretized on the dual mesh, with an upwind implicit iterative scheme

[ω k • ∇δu k ] n+1,p+1/2 r = 1 V r l, n r l•ω k >0 | l| n r l • ω k δu n+1,p+1 k,r + 1 V r l, n r l•ω k <0 | l| n r l • ω k δu n+1,p k,r (r) , (15)
where l is the edge intersecting the segment [rr ] joining the nodes r and r (two centers of the dual cells), see Fig. 2.

Proposition 3.2 (stability of the advection scheme). The advection scheme corresponding to [START_REF] Blanc | A positive scheme for diffusion problems on deformed meshes[END_REF] with the flux [START_REF] Brunner | Two-dimensional time dependent Riemann solvers for neutron transport[END_REF], and without source terms, is unconditionally stable (see proof in Appendix E).

It remains to discretize the term

ω k • [∇e] n+1,p+1 r
, to this end, we use the same nodal quadrature rule as before:

ω k • [∇e] n+1,p+1 r ≡ 1 V r j(r) e n+1,p+1 j C rj • ω k .
Thanks to our time discretization, we can provide an analytical expression for δu n+1,p+1 k,r . Indeed, defining δu r as the vector whose kth component is δu k,r , the Eq. ( 14) can be put in the form of a K × K system, local to each node:

I K -a n r ⊗ b n r δu n+1,p+1 r = A n r + B n+1,p r + C n+1,p+1 r , (16)
where I K is the identity matrix of rank K, a n r and b n r are the vectors of size K of components

a n k,r = 1/   1 + v∆tσ n t,r + v∆t Vr l, n r l•ω k >0 | l| n r l • ω k   , b n k,r = w k v∆t Vr l, ,n r l•ω k >0 | l| n r l • ω k . ( 17
)
and A n r , B n+1,p r and C n+1,p+1 r are the vectors of size K defined by:

A n k,r = a n k,r δu n k,r + v∆tδq n k,r , B n+1,p k,r = a n k,r v∆t( K k =1 w k 1 V r l, n r l•ω k <0 | l| n r l • ω k δu n+1,p k ,r -1 Vr l, n r l•ω k <0 | l| n r l • ω k δu n+1,p k,r
),

C n+1,p+1 k,r = -a n k,r v∆t ω k • [∇e] n+1,p+1 r . ( 18
)
215

Here, to prove the invertibility of the linear system ( 16), we use a corollary of the Sherman-Morrison lemma whose proof is provided in Appendix C. The statement of the corollary is: Let I K be the identity matrix of size K × K, and let a and b be two vectors (columns) of size K. The matrix

I K -a ⊗ b is invertible if a • b = 1,
and its inverse is:

I K -a ⊗ b -1 = I K + a ⊗ b 1 -a • b .
Using that k w k = 1, it is easy to show that a n r • b n r < 1, and consequently that the system ( 16) is invertible. Moreover, we are able to express its solution δu n+1,p+1 r

δu n+1,p+1 r = α n r + β n+1,p r + γ n+1,p+1 r , (19)
where α n r , β n+1,p r and γ n+1,p+1 r are the vectors of size K defined by:

α n r = I K + a n r ⊗b n r 1-a n r •b n r A n r , β n+1,p r = I K + a n r ⊗b n r 1-a n r •b n r B n+1,p r , γ n+1,p+1 r = I K + a n r ⊗b n r 1-a n r •b n r C n+1,p+1 r . ( 20
)
220

The expression of γ n+1,p+1 r can be simplified thanks to the following argument. Recalling that the directions are opposite in pairs and of equal weights, meaning that ∀k ∈

[1; K], ∃ k ∈ [1; K] s.t. k = k, wk = w k , ωk = -ω k , we infer from Eq. (17) that a i k,n = a i k,n and b i k,n = b i k,n , s.t. ∀k: γ n+1,p+1 k,r = -v∆t      a n k,r ω k + a n k,r K k =1 a n k ,r b n k ,r 1 -a n r • b n r ω k ≡0      • [∇e] n+1,p+1 r ,
which means that γ n+1,p+1 r reduces to:

γ n+1,p+1 r = C n+1,p+1 r .
It gives the discrete approximate solution of Eq. ( 8).

It remains now to solve the Eq. ( 4). To avoid a drastic parabolic CFL, we perform an implicit backward Euler time discretization of (11):

1 v e n+1 j -e n j ∆t + 1 V j r(j) f n+1 r • C jr + σ n a,j e n+1 j = s n j . (21)
Since

f n+1,p+1 r = K k=1 w k δu n+1,p+1 k,r ω k , the expression (19) of δu n+1,p+1 k,r yields 225 f n+1,p+1 r = N k=1 w k ω k α n k,r + β n+1,p k,r + γ n+1,p+1 k,r , = K k=1 w k ω k α n k,r + K k=1 w k ω k β n+1,p k,r -D n r [∇e] n+1,p+1 r , (22)
where D n r is an anisotropic tensor diffusion coefficient defined by Inserting the discrete flux [START_REF] Carré | A cell-centered Lagrangian hydrodynamics scheme on general unstructured meshes in arbitrary dimension[END_REF] into Eq. ( 21), we found the following discrete balance equation for e j :

D n r = v∆t K k=1 w k a n k,r ω k ⊗ ω k . ( 23
1 + v∆tσ n a,j e n+1,p+1 j - v∆t V j r(j) D n r [∇e] n+1,p+1 r • C jr = e n j + v∆t s n j - v∆t V j r(j) K k=1 w k ω k α n k,r + β n+1,p k,r • C jr , (24)
We have the important property for the tensor D • C jr in Eq. ( 24) is a discrete diffusion operator.

However, it would lead to a matrix which is not unconditionally invertible, and a diffusion scheme which is not positive [START_REF] Franck | Design and numerical analysis of asymptotic preserving schemes on unstructured meshes[END_REF]. This is why we reformulate Eq. ( 24), in order to use the diffusion scheme of our choice, and for which the outgoing fluxes are taken at the edges and not at the nodes:

r(j) D n r [∇e] n+1,p+1 r • C jr ≈ l(j) |l| D n l [∇e] n+1,p+1 l • n jl ,
where [∇e] l is the gradient of e evaluated at the center of the edge l. This is a valid approximation since In order to obtain such a discretization, we define, for each edge l (defined by the nodes r and r )

D n l = 1 2 D n r + 1 2 D n r
which is also a second-order accurate approximation. Thereby, we obtain a linear system on e:

1 + v∆tσ n a,j e n+1,p+1 j -v∆t V j l(j) |l| D n l [∇e] n+1,p+1 l • n jl = e n j + v∆t s n j -v∆t V j r(j) K k=1 w k ω k α n k,r + β n+1,p k,r • C jr , (25)
245

Finally, to complete the scheme, the last ingredient we need is the calculation of the diffusion operator:

v∆t V j l(j) |l| D n l [∇e] n+1,p+1 l • n jl .
To achieve that, we use an implicit solver for diffusion which is an extension of [START_REF] Blanc | A positive scheme for diffusion problems on deformed meshes[END_REF] to tensorial diffusion coefficient, and is presented in the following subsection.

Diffusion

Consider the following tensorial anisotropic diffusion problem :

∂ t e -∇ • (D∇e) = h in (0, T ) × Ω e = g on (0, T ) × ∂Ω • Ω bounded open set of R 2 , • h ∈ L 2 (Ω), 250 • g ∈ H 1/2 (∂Ω), • D : (0, T ) × Ω → M 2,2 is a bounded measurable function, s.t. D is a definite positive symmetric tensor, • e(t = 0, x) = e 0 ∈ H 1 (Ω) initial condition.
By integration of the equation in a cell j, we obtain

j ∂e ∂t + l∈∂j - l ∇e • ñjl = j h,
where ñjl = D l n jl , with n jl the outgoing normal (defined on previous section, refer 255 to Fig. 1) to the edge l of the cell j. We consider the flux (notations are given on Fig. 3) Let us define x jr 1 and x jr 2 s.t. ñjl has non-negative coefficients when decomposed into the basis (x jr 1 , x jr 2 ) (refer to Fig. 3 for different examples). Then, we express ñjl in the basis (x jr 1 , x jr 2 )

F j,l = - l ∇e • ñjl . ( 26
) × j × m • • r 1 r 3 r 2 r 4 ñjl ñml × × j m • • • ñjl ñml r 1 r 3 r 2 r 4
ñjl = λ j x jr 1 ||x jr 1 || + η j x jr 2 ||x jr 2 ||
, with λ j ∈ R + and η j ∈ R + . Inserting this decomposition into (26) yields

F j,l = - l λ j ∇e • x jr 1 ||x jr 1 || + η j ∇e • x jr 2 ||x jr 2 || .
Thanks to a Taylor expansion, we have

∇e • x jr i ||x jr i || = e(x r i ) -e(x j ) ||x jr i || + O(∆x),
and then,

F j,l = -|l| λ j e r 1 -e j ||x jr 1 || + η j e r 2 -e j ||x jr 2 || + O(∆x 2 ).
This last formulation involves e r 1 and e r 2 which are nodal evaluations of e at the nodes r i . To achieve the computation of e at these nodes, we use a least-squares procedure as in [START_REF] Coudière | Convergence rate of a finite volume scheme for a two dimensional convection-diffusion problem[END_REF]. We perform the same procedure for approximating the flux F m,l .

We can now define

F 1 = -|l| λ j e r 1 -
F j,l = µ 1 (e)F 1 -µ 2 (e)F 2 , (= -F m,l ),
with µ 1 + µ 2 = 1. It leads to the following expression for our discrete flux F j,l

F j,l = µ 1 |l| λ j ||x jr 1 || + η j ||x jr 2 || e j -µ 2 |l| λ m ||x mr 3 || + η m ||x mr 4 || e m + µ 1 R 1 -µ 2 R 2 ,
where R 1 and R 2 depend respectively on e r 1 , e r 2 and e r 3 , e r 4 . Choosing µ 1 and µ 2 as the solution of the following system

µ 1 + µ 2 = 1, µ 1 R 1 -µ 2 R 2 = 0,
we end up with a two-point but non-linear approximation of the flux F j,l

F j,l = µ 1 |l| λ j ||x jr 1 || + η j ||x jr 2 || e j -µ 2 |l| λ m ||x mr 3 || + η m ||x mr 4 || e m ,
since µ 1 and µ 2 now depend on e. We note that µ 1 and µ 2 are non-negative, provided that the nodal interpolations of e are non-negative. To this end, we apply the same trick as in [START_REF] Blanc | A positive scheme for diffusion problems on deformed meshes[END_REF]. As we explain in the next subsection, this procedure enforces the positivity of the diffusion scheme. However, a fixed point algorithm is required to solve the non-linearity of the fluxes. It remains to express the boundary conditions, for which we use the following approximation ∇e • x jt = e j -e t + O(∆x), using λ j x jt = D l n j,l = ñj,l . With e t = g(x t ), we have Finally, the implicit diffusion scheme is

F j,l = |l|λ j (e j -g(x t )). j ∂Ω × ñjl t •
1 ∆t (e n+1 -e n ) + M (e n+1 )e n+1 = h + g e n+1 ≥ 0 , with Φ j,l = µ 1 |l| λ j ||x jr 1 || + η j ||x jr 2 || , Φ m,l = µ 2 |l| λ m ||x mr 3 || + η m ||x mr 4 || , and        [M (e)] jj = 1 V j   l∈M\∂Ω Φ j,l + l∈M∩∂Ω |l|λ j   [M (e)] jm = - Φ m,l V j if j = m .

Properties of the diffusion scheme and of the global scheme

The diffusion scheme fulfills the following properties.

Proposition 3.4 (Positivity (monotony) of the diffusion scheme). The diffu-270 sion scheme is positive and is unconditionnaly stable.

This property is also called monotony in other publications [START_REF] Lipnikov | A monotone finite volume method for advection-diffusion equations on unstructured polygon meshes[END_REF][START_REF] Nikitin | A monotone nonlinear finite volume method for advection-diffusion equations on unstructured polyhedral meshes in 3D. Russian[END_REF][START_REF] Sheng | A finite volume scheme for diffusion equations on distorted quadrilateral meshes[END_REF]. It means that if e n ≥ 0, h ≥ 0 and g ≥ 0, then e n+1 ≥ 0. The proof is displayed in Appendix G.

Proposition 3.5 (Convergence).

Providing that e n ≥ 0, h ≥ 0 and g ≥ 0, the scheme has a unique solution, and the fixed-point strategy can be proven to converge under a parabolic CFL restriction.

The proof can be found in [START_REF] Blanc | A positive scheme for diffusion problems on deformed meshes[END_REF].

Proposition 3.6 (Conservation of the diffusion scheme). The diffusion scheme is conservative.

See the proof in Appendix F.

We are now able to demonstrate the following property on the global balance of energy [START_REF] Chaland | Discrete Ordinates Method for the Transport Equation Preserving One-Dimensional Spherical Symmetry in Two-Dimensional Cylindrical Geometry[END_REF]. Proposition 3.7 (Conservation of the global scheme). The global scheme is conservative.

See the proof in Appendix H. We finally prove that the global scheme is Asymptotic Preserving. Proposition 3.8 (Asymptotic Preserving behaviour of the scheme). In the limit v∆tσ n t,r >> 1, v∆tσ n a,j << 1, v∆t|s n j | << e n j , and v∆t|q n k,r | << |δ n k,r |, the system ( 14)-( 25) tends to the following discrete implicit diffusion operator e n+1 -e n = ∆tM e n+1 , where e ∈ R #j is the vector of components e j and of size the number of cells #j, and M ∈ R #j × R #j is the diffusion matrix corresponding to the discretization of

l(j) |l| 1 3σ s,j
[∇e] n+1 l

• n jl with the algorithm described in this section.

See the proof in Appendix J.

Remark 3.2 (About the Asymptotic Preserving property). Despite that the transport scheme used to solve Eq. ( 8) is a simple upwind discretization, the global scheme is Asymptotic Preserving. This is due to the fact that all the δu k tends to zero in the diffusive regime, and then only the equation on e matters.

Summary of the scheme

Time iterations (n):

calculation of a n r , b n r , and α n r , with Eq. ( 17), ( 18) and ( 20), Do (sub-iterations p):

• construction of the diffusion matrix (see section (3.2)),

• initialization of δu n+1,p=0 k,r = δu n k,r , • calculation of β n+1,p r with ( 18) and ( 20),

• construction of the system on e (25),

• calculation of e n+1,p+1 j by solving the linear system (25),

• deduction of γ n+1,p+1 r with ( 18) and ( 20 It is important to notice that there is only one loop on p for the fixed point algorithm for e and the implicit solve of δu.

Numerical experiments

In this last section, we are first interested in four 1D tests, which aim to validate the numerical scheme. We perform them on a Cartesian mesh and on an unstructured mesh, and compare the results with those of a Discontinuous Finite Element (DFE) code on a Cartesian mesh, and to an analytical solution if it exists. In order to test the scheme efficiency on a fully 2D configuration, we perform a fifth and final test: the lattice problem. This problem has been studied in numerous publications in the past decades [START_REF] Buet | Asymptotic preserving schemes on distorted meshes for Friedrichs systems with stiff relaxation: application to angular models in linear transport[END_REF][START_REF] Crockatt | An arbitraryorder, fully implicit, hybrid kinetic solver for linear radiative transport using integral deferred correction[END_REF][START_REF] Hermeline | A discretization of the multigroup PN radiative transfer equation on general meshes[END_REF][START_REF] Himpe | Hierarchical approximate proper orthogonal decomposition[END_REF][START_REF] Mcclarren | Robust and accurate filtered spherical harmonics expansions for radiative transfer[END_REF][START_REF] Radice | A new spherical harmonics scheme for multi-dimensional radiation transport I. Static matter configurations[END_REF][START_REF] Schäfer | Diffusive corrections to $p_n$ approximations[END_REF][START_REF] Seibold | StaRMAP-A Second Order Staggered Grid Method for Spherical Harmonics Moment Equations of Radiative Transfer[END_REF], and the reader is invited to refer to them for comparison. Moreover, we also provide the results obtained with the DFE method.

Before detailing the test problems and the results, we give a short description of the DFE method we use for the sake of comparison.

Short description of the Discontinuous Finite Element (DFE) method

The P1 Discontinuous Finite Element (P1 DFE) is widely used to solve the transport equation especially for photons. Its ability to give good results in transparent and opaque media and in the last case, without meshing at the mean free path length, explains this popularity. Several extensions to general polygonal (resp. polyhedral) meshes in 2D (resp. 3D) which degenerate in 1D into P1 DFE and still possess the diffusion limit have been proposed: the Upstream Corner Balance, Simple Corner

Balance [START_REF] Adams | Subcell Balance Methods for Radiative Transfer on Arbitrary Grids[END_REF], Bailey methods [START_REF] Bailey | The piecewise linear discontinuous finite element method applied to the RZ and XYZ transport equations[END_REF]. The unknowns of these three methods are the values of the radiative intensity at the vertices of the cell. The main features of the DFE method are:

• Accuracy The P1 DFE is second order accurate in space, at least on Cartesian mesh. As a consequence, positivity is not guaranteed if no limiter is used.

• Diffusion limit 1D P1 DFE method has the diffusion limit [START_REF] Larsen | Asymptotic solutions of numerical transport problems in optically thick, diffusive regimes II[END_REF]. The limit scheme is a valid discretization of the diffusion equation by P1 continuous finite elements. 1D P1 DFE also captures the boundary layers accurately even using a crude mesh [START_REF] Larsen | Asymptotic solutions of numerical transport problems in optically thick, diffusive regimes II[END_REF]. Note that the P0 Discontinuous Finite Element (the well known donor cell scheme), does not recover the diffusion limit.

• Implicitness In order not to constraint the time step, the scheme is made implicit in time. The source term which couples all directions is usually lagged. This allows to solve the angular discretized problem by direction. This technique is named the source iteration method. In opaque media, the convergence of the iterated process can be arbitrary slow. Solving between iterations some low order problem (DSA [START_REF] Alcouffe | Diffusion synthetic acceleration methods for the diamond-differenced discrete-ordinates equations[END_REF], TSA [START_REF] Ramone | A Transport Synthetic Acceleration Method for Transport Iterations[END_REF]) is often used to accelerate the convergence.

The discretization of the transport term leads to solve by direction a large linear system whose the unknowns are the values of the radiative intensity at the vertices of the cells. The properties of the transport operator allow to make this system block triangular by ordering the cells. Solving the whole linear system directly block by block is much cheaper than solving it with an iterative method. By cell, each block is a N (j) × N (j) matrix with N (j) the number of the vertices of the cell. Finally, compute the unknowns needs to solve a N (j) × N (j) linear system per cell. Such ordering exists in 2D if only convex cells are present. In 3D, this is no longer the case. Graph algorithms allow to recover such an ordering at the expense of lagging some intensities at the previous iteration in the source iteration method [START_REF] Plimpton | Parallel Sn sweeps on unstructured grids: Algorithms for prioritization, grid partitioning, and cycle detection[END_REF].

In this article, the Bailey method for the spatial discretization and an acceleration technique TSA [START_REF] Ramone | A Transport Synthetic Acceleration Method for Transport Iterations[END_REF]) is used.

1D tests

Regarding the first four tests:

• The first test free streaming limit aims to check the consistency and the good behavior (stability, precision) of the new scheme within the free streaming limit.

• The second test diffusion in homogeneous media aims to check the Asymptotic-Preserving property of the new scheme, as well as its ability to support arbitrarily large time steps as would do a classical solver for implicit diffusion.

• The third test transport/diffusion interfaces aims to evaluate the behavior (robustness, precision) of the new scheme in the presence of an heterogeneous media, made of a succession of transport zones and diffusion zones.

• The fourth test stationary boundary layer aims to check the ability of the new scheme to solve a non-trivial stationary problem, and to evaluate the accuracy of the results on a coarse mesh in the presence of a boundary layer.

In the figures, one can observe two or three curves: one curve corresponding to our finite volume scheme (named FV), one curve corresponding to the DFE method, and for the tests one and two, a last curve corresponds to the analytical solution.

To compare Cartesian mesh to unstructured mesh results in 1D, we use the spatial discretization of the boundary obtained on the Cartesian mesh as the basis for the unstructured discretization. For instance, when we say that the discretization corresponds to 100 cells along the x-axis, it means that if we use 100 edges of squares (resp. 100 edges of triangles) discretizing the y = y min boundary for the Cartesian mesh (resp. for the unstructured mesh).

First test: free streaming limit

The analytical solution is obtained using the characteristic method. We work on the domain Ω = [-1, +1] × [-0.08, 0.08]. The type of unstructured mesh we use for this test is shown on Fig. 5. Initial condition is

u(t = 0, x, µ) = 1 1 x∈[-0.5;+0.5] , q = σ a = σ s = 0.
Zero incoming flux conditions are applied on boundaries where n x = (1, 0) T is the unit normal along the x axis. Final time of the calculation is t = 0.2, and we use a time-step ∆t = 0.02. The number of cells along the x-axis is 100 for the Cartesian cells. We compare the plot e(x) obtained with our new method (FV) to the analytical solution, and the DFE result, for N = 1 (K = 4 directions, Fig. 6 and7), and N = 6 (K = 144 directions Fig. 8 and9). Results are very close to each other and in good agreement with the analytical solution, whatever the kind of mesh (Cartesian or unstructured). Our scheme is more diffusive than DFE, but DFE exhibits small overshoots and undershoots, which makes the solution goes below zero, which relies on the second-order accuracy of these methods.

u(t, x = -1, y, ω • n x > 0) = u(t, x = 1, y, ω • n x < 0) = 0,
The angular convergence is illustrated for N = 6 (K = 144), on Fig. 8 and9.

We have done tests with very large CFL numbers, to verify that the new method indeed restores the purely transport regime in an unconditionally stable way, whereas the algorithmic complexity of the solver is comparable to the one of an implicit diffusion type solver.

Second test: diffusion in homogeneous media

In this test, initial and boundary conditions are unchanged, except for the scattering cross section, which is taken as σ s = 10 5 . In these conditions, the problem is steep enough for the diffusion limit to be valid. This problem is intended to verify that our new scheme is Asymptotic Preserving. In addition, we do not expect any dependence of the result w.r.t. angular discretization. The analytical solution of the test is the solution of the diffusion equation ∂ t e = D∂ 2

x 2 e with D = v/(3σ s ) associated with the initial condition e(t = 0, x) = 1 1 x∈[-0,5;+0,5] . The time-step ∆t = 1 is used for all the results in this section. The type of unstructured mesh used for this test is the same than for the previous one (see Fig. 5, with The solution does not depend on the angular discretization, and is superimposed to the analytical solution, even on these coarse meshes. The unstructured mesh result is as satisfactory as that of the Cartesian mesh, and respects the maximum principle on e. Moreover, it behaves like a classical implicit diffusion solver, whereas the DFE 415 method requires the use of an acceleration technique (DSA or TSA).

Third test: transport/diffusion interfaces

This problem requires a sufficient number of angular directions to be interesting, and we take N = 11 (K = 484 directions) for the finite volume method and K = 496 directions for the in the DFE method in the following1 . For this test, the scattering cross section is discontinuous: The initial condition is

σ s = 10
u(t = 0, x, µ) = 1 1 x∈[-0.1;+0.1] .
The unstructured mesh is depicted on Fig. 14. We draw the profile e(t, x) at time t = 1, 10 and 100, and compare our new scheme with the DFE method. Results for t=1 are given on Fig. 15 and16.

420

All the curves agree well for |x| > 0, 4. The precision with which these waves are captured in the x < 0.3 zone obviously does not influence the precision for x ≥ 0.3, which is reassuring. In particular, the height of the flat area for x ∈]0.4; 0.5[ fits perfectly between methods, whatever the mesh (Cartesian or unstructured).

Results for t=10 and for t=100 are displayed on Fig. 17 and 20, respectively. The results of our method fit very well the results of the DFE method. In addition, the positions of the discontinuities of the derivative of e match those of the cross section, as expected. We observe that our method fulfills the maximum principle on e, while small oscillations occur in the DFE calculations, again because of the second order accuracy of the method.

Fourth test: stationary boundary layer

This test is freely inspired from the analysis of Larsen & Morel [START_REF] Larsen | Asymptotic solutions of numerical transport problems in optically thick, diffusive regimes II[END_REF]. It aims at verifying the ability of our method to handle unresolved boundary layers due to the steep variation of the cross sections. The objective is to be able to calculate the correct slope of e on coarse meshes. As for the previous problem, this problem requires a sufficient number of angular directions to be interesting. We take N = 4 (K = 64 directions in 2D) in the following. We are then interested in the profile of the moment of order 0, e, in particular close to the interface located at x = 0.8 where a boundary layer develops.

For this test, we have : In order to simplify the initialization (of δu and e), we divide the calculation domain in three different areas.

The results on Cartesian and unstructured meshes are displayed on Fig. 22. For the comparison purpose, the DFE method is used with a very fine mesh, to obtain a DFE reference solution. Our results, on both structured and unstructured grids, agree very well with the reference DFE calculation. We now focus on the behaviour of e at the boundary layer, close to x = 0.8. Results are given on Fig. 24.

On this coarse mesh, the discontinuity of e at x = 0.8 is poorly restored in the FV code, but the good agreement of the slope in the diffusion zone (x > 0.8) testifies that the net flux at the interface is correct. With the DFE code, we can observe a discontinuity of e at x = 0.8. As DFE methods are known to very well capture these boundary layers [START_REF] Larsen | Asymptotic solutions of numerical transport problems in optically thick, diffusive regimes II[END_REF], the very good agreement with our method is very satisfactory.

It means that the net flux at the interface is correctly calculated by the new method.

Then we perform a spatial convergence study on the test 4. Results are displayed on Fig. [START_REF] Chaland | Discrete Ordinates Method for the Transport Equation Preserving One-Dimensional Spherical Symmetry in Two-Dimensional Cylindrical Geometry[END_REF]. It shows that the slope is correctly calculated even on the 24 cells mesh. In addition, it shows the quick convergence of our new algorithm to the reference DFE solution obtained on a very fine mesh (20000 degrees of freedom).

Lattice problem

We consider the lattice problem as described in [START_REF] Brunner | Two-dimensional time dependent Riemann solvers for neutron transport[END_REF]. The domain is

Ω = [0, 7] × [0, 7].
We divide the domain as in Fig. 26. The white area, bulk of the lattice, corresponds to a scattering material and eleven blue area to absorbing regions. The black region in the center is a scattering region, where a source of photons is emitted. Let Ω A be the union of the eleven striped squares (absorption region) and let Ω C = [START_REF] Adams | Asymptotic analysis of a computational method for time-and frequency-dependent radiative transfer[END_REF][START_REF] Alcouffe | Diffusion synthetic acceleration methods for the diamond-differenced discrete-ordinates equations[END_REF] The source is concentrated in the center (in the black square), such as

s(x) = 1 if x ∈ Ω C s(x) = 0 else
We used a 140 × 140 mesh, with ∆t = 10 -2 . The test is performed on Cartesian and unstructured meshes for our method, and on a Cartesian mesh for the DFE method. The kind of unstructured mesh used for this first test is given on Fig. 27.

We plot the map of e(x, t) at time t = 3.2 on Fig. 28 to Fig. 33. Our results are 475 comparable to those obtained in other publications [START_REF] Buet | Asymptotic preserving schemes on distorted meshes for Friedrichs systems with stiff relaxation: application to angular models in linear transport[END_REF][START_REF] Crockatt | An arbitraryorder, fully implicit, hybrid kinetic solver for linear radiative transport using integral deferred correction[END_REF][START_REF] Hermeline | A discretization of the multigroup PN radiative transfer equation on general meshes[END_REF][START_REF] Himpe | Hierarchical approximate proper orthogonal decomposition[END_REF][START_REF] Mcclarren | Robust and accurate filtered spherical harmonics expansions for radiative transfer[END_REF][START_REF] Radice | A new spherical harmonics scheme for multi-dimensional radiation transport I. Static matter configurations[END_REF][START_REF] Schäfer | Diffusive corrections to $p_n$ approximations[END_REF][START_REF] Seibold | StaRMAP-A Second Order Staggered Grid Method for Spherical Harmonics Moment Equations of Radiative Transfer[END_REF].

We compare our results in 2D on Cartesian mesh, unstructured mesh and a test carried out with DFE on Fig. 34 and Fig. 35. The left figure represents the results obtained on Cartesian mesh with the new scheme, the right figure represents the results obtained with the new scheme on unstructured mesh, and finally, the middle

Conclusion

We have presented in this article a finite-volume scheme for a linear transport equation derived from a linearization of the radiative transfer equations. This scheme as been implemented on 2D unstructured meshes, and satisfies the following proper- P1 to be consistent on general meshes, P2 to enforce the conservation of radiative energy, P3 to be able to handle all the radiation regimes from free-streaming to diffusion (Asymptotic Preserving), P4 to have the radiative energy degrees of freedom located at the centers of the elements to ensure the compatibility with the hydrodynamic scheme.

P5 to allow us to use the limit diffusion scheme of our choice.

As our scheme is implicit, we used a system of sub-iterations to avoid to solve a global linear system (of size equal to the cells number × directions number), while remaining 495 stable. In addition, this sub-iterations system is acting as a fixed point loop, in order to get rid of the non-linearity of the diffusion scheme. Thanks to the properties of our diffusion scheme, we enforce the energy to remain positive in this regime, for which radiation and matter temperature are strongly coupled. We have conducted numerical 1D tests on structured and unstructured 2D meshes, which assess that the 500 expected properties are respected. Finally, we carried out the Lattice problem test to compare our results with existing methods, in particular the Discontinuous Finite Element (DFE) discretization. In the future, we aim at designing a scheme which enforces the positivity of the energy in all regimes. In the following we prove that

         1 |∂j| ∂j ϕn = 1 |∂j| l(j) |l|ϕ l n jl + O(h 2 ), 1 |∂j| ∂j ϕn = 1 |∂j| r(j) ϕ r C jr + O(h 2 ), (A.1)
ϕ l n jl + (∇ϕ) l (x -x l ) + O(h 2 ) , Taylor expansion, = l(j) |l|ϕ l n jl + O(h 3 ), l(j) (x -x l ) = 0.
Since |∂j| = O(h), it proves the first line of (A.1). Since there is at least d independent vectors ω k , the set of ω k spans R d , meaning that x is orthogonal to all vectors of R d . So x = 0. Then, we suppose that a • b = 1 and we develop the product of matrices:

I K -a ⊗ b -1 = I K + a ⊗ b 1 -a • b .
I K -a ⊗ b I K + a ⊗ b 1 -a • b = I K -a ⊗ b + a ⊗ b 1 -a • b - (a ⊗ b) (a ⊗ b) 1 -a • b = I K ,
that proves the corollary Appendix F. Proof that the diffusion scheme is conservative 535

We consider our diffusion scheme: We denote by S the set of the index of the cells. Then, for all j ∈ S, we have w k δu m,q k,r = 0 ∀m and ∀q.

because of the constraints on the angular quadrature [START_REF] Bailey | A Piecewise Bi-Linear Discontinuous Finite Element Spatial Discretization of the Sn Transport Equation[END_REF]. The reformulation [START_REF] Chaland | Discrete Ordinates Method for the Transport Equation Preserving One-Dimensional Spherical Symmetry in Two-Dimensional Cylindrical Geometry[END_REF] does not change this analysis, which implies the result.

Figure 2 :

 2 Figure 2: Notations (dual mesh)

Proposition 3 . 1 (

 31 Preservation of the meaning of e during time). With such a temporal discretization, K k=1 w k δu n k,r = 0, ∀n, ∀r, whatever the nodal approximation of [ω k • ∇δu k ] r and [∇e] r (see proof in Appendix I).

  ) and α n k,r , β n+1,p k,r and γ n+1,p+1 k,r are the kth components of vectors α n r , β n+1,p r and γ n+1,p+1 r .

  Proposition 3.3 (definite positiveness of the diffusion tensor). D n r is definite positive (see proof in Appendix B). The term r(j) D n r [∇e] n+1,p+1 r

  r(j) ϕ r C jr and l(j) |l|ϕ l n jl are two second-order accurate approximations of ∂j ϕn (we recall the proof in Appendix A). It is the subject of the following section to explain how we approximate [∇e] n+1,p+1 l . Since, the diffusion operator comes from the expression of γ n+1,p+1 r , and more particularly from D n r [∇e] n+1,p+1 r , it remains to define D n l from D n r .

Figure 3 :

 3 Figure 3: Expression of ñjl and ñml

Figure 4 :

 4 Figure 4: Boundary conditions notations

Figure 5 :

 5 Figure 5: Unstructured mesh used for the free streaming problem.

Figure 6 :

 6 Figure 6: Comparison of e between finite volume and DFE, on Cartesian mesh, for the free streaming problem, with N = 1 (K = 4) and ∆t = 0.02.

Figure 7 :

 7 Figure 7: Comparison of e between finite volume and DFE, on unstructured mesh, for the free streaming problem, with N = 1 (K = 4) and ∆t = 0.02.

2 1 1

 1 |x|∈[0.3;0.4] + 10 3 1 1 |x|∈[0.5;0.6] + 10 4 1 1 |x|∈[0.7;0.8] 1 1 |x|≥0.9 , and we have q = 0, σ a = 0.

Figure 8 :

 8 Figure 8: Comparison of the angular moment e between finite volume and DFE, on Cartesian mesh, for the free streaming problem, with N = 12 (K = 576 for FV and K = 544 for DFE) and ∆t = 0.02.

  and 18 and Fig.19

  425

Figure 9 :

 9 Figure 9: Comparison of the angular moment e between finite volume and DFE, on unstructured mesh, for the free streaming problem, with N = 12 (K = 576 for FV and K = 544 for DFE) and ∆t = 0.02.

Figure 10 :

 10 Figure 10: Comparison of the angular moment e between finite volume and DFE, on Cartesian mesh, for the diffusion in homogeneous media problem, with N = 1 (K = 4) and ∆t = 10.

Figure 11 :

 11 Figure 11: Comparison of the angular moment e between finite volume and DFE, on unstructured mesh, for the diffusion in homogeneous media problem, with N = 1 (K = 4), ∆t = 10.

Figure 12 :

 12 Figure 12: Comparison of the angular moment e between finite volume and DFE, on Cartesian mesh, for the diffusion in homogeneous media problem, with ∆t = 10.

Figure 13 :

 13 Figure 13: Comparison of the angular moment e between finite volume and DFE, on unstructured mesh, for the diffusion in homogeneous media problem, with ∆t = 10.

2

 2

Figure 14 :

 14 Figure 14: Unstructured mesh used for the transport/diffusion interfaces problem.

Figure 15 :

 15 Figure 15: Comparison of the angular moment e between finite volume and DFE, on Cartesian mesh, for the transport/diffusion interfaces problem, with N = 11 (K = 484) and ∆t = 1, plot at time t = 1.

Figure 16 :

 16 Figure 16: Comparison of the angular moment e between finite volume and DFE, on unstructured mesh, for the transport/diffusion interfaces problem, with N = 11 (K = 484) and ∆t = 1, plot at time t = 1.

Figure 17 :

 17 Figure 17: Comparison of the angular moment e between finite volume and DFE, on Cartesian mesh, for the transport/diffusion interfaces problem, with N = 11 (K = 484) and ∆t = 1, plot at time t = 10.

Figure 18 :

 18 Figure 18: Comparison of the angular moment e between finite volume and DFE, on unstructured mesh, for the transport/diffusion interfaces problem, with N = 11 (K = 484) and ∆t = 1, plot at time t = 10.

Figure 19 :

 19 Figure 19: Comparison of the angular moment e between finite volume and DFE, on Cartesian mesh, for the transport/diffusion interfaces problem, with N = 11 (K = 484) and ∆t = 1, plot at time t = 100.

Figure 20 : 2 r 2 r

 2022 Figure 20: Comparison of the angular moment e between finite volume and DFE, on unstructured mesh, for the transport/diffusion interfaces problem, with N = 11 (K = 484) and ∆t = 1, plot at time t = 100.

Figure 21 :

 21 Figure 21: Unstructured mesh used for test four: stationary boundary layer.

Figure 22 :

 22 Figure 22: Comparison of the angular moment e between the DFE method and the FV method in Cartesian and unstructured meshes, for the stationary boundary layer problem, with N = 4 (K = 64) and ∆t = 10 9 . The DFE reference uses 10000 cells (20000 degrees of freedom).

Figure 23 :

 23 Figure 23: Comparison of the angular moment e between finite volume and DFE, on Cartesian mesh, for the stationary boundary layer problem, with N = 4 (K = 64) and ∆t = 10 9 . Zoom around x = 0.8.

  Proof:We start by noting that(a ⊗ b) (a ⊗ b) = (a • b) (a ⊗ b) ; indeed: For i, j from 1 to K: [(a ⊗ b) (a ⊗ b)] i,j = K k=1 a i b k a k b j = (a • b) [(a ⊗ b)] i,j .

Figure 24 :

 24 Figure 24: Comparison of the angular moment e between finite volume and DFE, on unstructured mesh, for the stationary boundary layer problem, with N = 4 (K = 64) and ∆t = 10 9 . Zoom around x = 0.8.

Figure 25 :

 25 Figure 25: Comparison of the angular moment e of the FV code, for different numbers of cells, for test 4: stationary boundary layer, with N = 4 (K = 64) and ∆t = 10 9 . The DFE reference use 10000 cells (20000 degrees of freedom). Zoom around x = 0.8.

Figure 26 :

 26 Figure 26: Lattice calculation domain.

Figure 28 :

 28 Figure 28: Cartesian mesh 140 × 140, lattice problem, map of e, N = 1 (K = 4), time t = 3.2.

Figure 30 :Figure 31 :

 3031 Figure 30: Cartesian mesh 140 × 140, lattice problem, map of e, N = 3 (K = 36), time t = 3.2.

Figure 33 :

 33 Figure 33: Cartesian mesh 140 × 140, lattice problem, map of e, N = 8 (K = 256), time t = 3.2.

Figure 34 :

 34 Figure 34: Comparison of results on lattice problem. Left: FV method, Cartesian mesh 140×140. Middle: DFE code, Cartesian mesh 140 × 140. Right: FV method unstructured mesh. K = 4, time t = 3.2.

  e j ||x jr 1 || + η j e r 2 -e j ||x jr 2 || and F 2 = -|l| λ m e r 3 -e m ||x mr 3 || + η m e r 4 -e m ||x mr 4 || , which gives us two consistent approximations of the flux. It yields our expression of the flux from the cell j to the cell m as a convex combination of F 1 and -F 2

  1 ∆t (e n+1 -e n ) + M (e n+1 )e n+1 = h + g e n+1 ≥ 0 With Φ j,l = µ 1 |l| λ j ||x jr 1 || + η j ||x jr 2 || , Φ m,l = µ 2 |l| λ m ||x mr 3 || + η m ||x mr 4 ||.

The angular quadrature formula implemented being different for the two methods, we cannot use exactly the same number of directions. However, we checked that the angular convergence is achieved.

figure is that obtained with the DFE.

The initial condition is :

u(t = 0, x, µ) = 0.

The code used being unsteady, we reach the solution of this stationary problem by taking in practice:

u (t = 0, x, µ) = 0, t = 10 9 , Appendix D. Proof that the advection scheme is conservative At convergence of the sub-iterations, we have an advection scheme of the form:

where l is the edge intersecting segment [rr ].

equal to zero. It remains to show that the value V r δu k,r is conserved.

In the following, we consider periodic boundary conditions. Let r and r be two adjacent nodes. We write the flux from node r to node r , we get

525 Appendix E. Proof that the advection scheme is unconditionally stable

The advection scheme has the form

We can rewrite it with the expression of

, which gives

Then, factorising by δu n+1,p+1 k,r

, we get ⇒ δu n+1,p+1 k,r

.

Given that

depending on δu n k,r and δu n+1,p k,r

, ∀ν, which means ∀∆t.

We consider periodic boundary conditions, so we have [M (e)] jj = 1

We need to prove that if h = 0, then ∀n ≥ 0,

, where V j corresponds to the volume of cell j.

The continuous version of it is d

dt Ω e = 0.

Appendix G. Proof that the diffusion scheme respect positivity and is 545 unconditionally stable

Positivity. Setting h = g = 0, our diffusion scheme writes: e n+1 = e n -∆tM (e n+1 )e n+1 , with Thus, the diffusion scheme is positive.

Stability. Concerning the stability of the scheme, we refer to [START_REF] Blanc | A positive scheme for diffusion problems on deformed meshes[END_REF].

550

In this article, it is shown that for h ≥ 0, g ≥ 0 and e n ≥ 0, it exists a solution e n+1 , and the fixed point strategy converges to it.

Appendix H. Proof that the global scheme is conservative

Summing over all the cells Eq. ( 25), considering the absorption term σ a = 0 and 555 the source s = 0, one gets

The term

• n jl = 0 because of property 3.6.

C jr = 0, refer for instance to [START_REF] Carré | A cell-centered Lagrangian hydrodynamics scheme on general unstructured meshes in arbitrary dimension[END_REF].

We prove it by induction. Initially, K k=1 δu 0 k,r = 0. Let us assume that

We consider the scheme:

Then, we multiply the equation by w k and we sum over k.

Using that k w k = 1 , we have

Then, we have

Then, we have that

Appendix J. Proof that the global scheme is Asymptotic Preserving

From the hypothesis σ s → ∞, σ a → 0 and q → 0, up to O( 1

), we define the discrete Knudsen number ε = 1/(v∆tσ n t,r ). In addition, we assume that the initial 575 condition is well-prepared, meaning that δu n k,r = O(ε). Then from definition [START_REF] Buet | Proof of uniform convergence for a cell-centered AP discretization of the hyperbolic heat equation on general meshes[END_REF], we infer that: