Rémi Feuillet
email: remi.feuillet@inria.fr

Olivier Coulaud
email: olivier.coulaud@cenaero.be

Adrien Loseille
email: adrien.loseille@inria.fr

ANISOTROPIC ERROR ESTIMATE FOR HIGH-ORDER PARAMETRIC SURFACE MESH GENERATION

Keywords: Surface mesh generation, High-order error estimates, High-order mesh generation

Parametric surface mesh generation is one of the crucial step of the computational pipeline. Standard techniques, that are now mature, control the deviation to the tangent plane by using intrinsic quantities as the minimum and maximum curvatures. However, for high-order meshes, deriving intrinsic quantities that have the ability to control the mesh generation process is much more challenging. Indeed, those provided by the first and second fundamental forms of a surface are not sufficient when high order curved meshes are employed. In this paper, we introduce a new set of error estimates for high-order surface mesh generation. It is based on performing a Taylor expansion of the underlying surface in the tangent plane. The independence to the parametric space is obtained by using an inversion formula. High-order terms of this expansion are then used to derive an optimal metric by using the log-simplex approach. Examples are shown to prove the efficiency of the method.

INTRODUCTION

In industrial applications, the definition of the computational domain (or of a design) is provided by a continuous description composed by a collection of patches using a CAD (Computer Aided Design) system. If several continuous representations of a patch exist via an implicit equation or a solid model, we focus on the boundary representation (BREP). In this description, the topology and the geometry are defined conjointly. For the topological part, a hierarchical description is used from top level topological objects to lower level objects, we have: model → bodies → faces → loops → edges → nodes Each entity of upper level is described by a list of entities of lower level. This is represented in Figure 1 for an Onera M6 model, where a face, a loop and corresponding edges are depicted. Note that most of the time, only the topology of a face is provided, the topol- ogy between all the faces (patches) needs to be recovered. This piece of information is needed to have a watertight valid surface mesh on output for the whole computational domain. This step makes the surface mesh generation of equal difficulty as volume mesh generation and have been shown to be not trivial [START_REF] Alleaume | Automatic Non-manifold Topology Recovery and Geometry Noise Removal[END_REF]. For node, edge, and face, a geometry representation is also associated to the entity. For node, it is generally the position in space, while for edge and face a parametric representation is used. It consists in defining a mapping from a bounded domain of R 2 onto R 3 such that (x, y, z) = σ(u, v) where (u, v) are the parameters (Figure 2). Generally, σ is a NURBS function (Non-uniform rational B-spline) as it is a common tool in geometry modeling and CAD systems [START_REF] Piegl | The NURBS Book[END_REF]. From a conceptual point of view, meshing a parametric surface consists in meshing a 2D domain in the parametric space. The linear case has been widely studied for years [START_REF] De Cougny | Surface meshing using vertex insertion[END_REF][START_REF] Tristano | Advancing front surface mesh generation in parametric space using a riemannian surface definition[END_REF][START_REF] Miranda | Mesh generation on high-curvature surfaces based on a background quadtree structure[END_REF][START_REF] Wang | EQSM: An efficient high quality surface grid generation method based on remeshing[END_REF][START_REF] Laug | Some aspects of parametric surface meshing[END_REF][START_REF] Siqueira | A Hierarchical Adaptive Mesh Generation Strategy for Parametric Surfaces Based on Tree Structures[END_REF][START_REF] De Siqueira | An Adaptive Parametric Surface Mesh Generation Method Guided by Curvatures[END_REF][START_REF] Dey | A robust conforming NURBS tessellation for industrial applications based on a mesh generation approach[END_REF] and some approaches consist of meshing the 2D domain according to a curvature-based metric [START_REF] Borouchaki | Maillage de surfaces paramétriques. Partie I: Aspects théoriques[END_REF]. The use of this metric enables an independence to the used parameters space as the curvature is an intrinsic data. The generation of highorder meshes is on the contrary relatively new. The most common idea is to generate a linear mesh and then to project the high-order nodes on the geometry. However, the position of the nodes may not be suitable for a high-order representation of the boundary. For this reason, optimization procedures are applied to the mesh to improve its shape [START_REF] Toulorge | Optimizing the geometrical accuracy of curvilinear meshes[END_REF][START_REF] Ruiz-Gironès | Defining an L 2 -disparity Measure to Check and Improve the Geometric Accuracy of Non-interpolating Curved High-order Meshes[END_REF][START_REF] Ruiz-Gironès | High-order mesh curving by distortion minimization with boundary nodes free to slide on a 3D CAD representation[END_REF][START_REF] Ruiz-Gironès | Generation of curved high-order meshes with optimal quality and geometric accuracy[END_REF][START_REF] Turner | High-order mesh generation for CFD solvers[END_REF]. The procedure can be done by solving an optimization problem or performing a spring analogy. This can also be performed directly in the parameters space [START_REF] Gargallo-Peiró | A distortion measure to validate and generate curved high-order meshes on CAD surfaces with independence of parameterization[END_REF] by minimizing a distortion measure. In all these cases, this is r-adaptation which is performed. The intent of the presented work is to provide tools for high-order parametric surface mesh generation. The use of the high-order estimates developed in [START_REF] Coulaud | Very High Order Anisotropic Metric-Based Mesh Adaptation in 3D[END_REF] will provide a node distribution which will be specifically tailored for the high-order with a given threshold. Note that high order surface meshes have a large set of applications. It can be naturally used as an input for the generation of 3D curved meshes. In our case, high-order surface mesh is used advantageously as a surrogate CAD (geometry) model. Indeed, it provides fast forward and inverse evaluation as required in classic linear mesh adaptation.

The first section briefly recalls the framework used for the high-order error estimates and how a metric can be deduced from it. The second section shows how the standard parametric surface mesh generation can be extended to the high-order using this framework.

HIGH ORDER METRIC BASED

MESH ADAPTATION

Continuous mesh framework

In this section, some basic tools concerning the continuous mesh framework are recalled. For further details about these techniques, one refers to [START_REF] Loseille | Continuous mesh framework part I: well-posed continuous interpolation error[END_REF][START_REF] Loseille | Continuous mesh framework part II: validations and applications[END_REF] and the references therein. The main idea underlying the metric based mesh adaptation is to change the way distances are computed, via a continuous metric field M = (M(x)) x∈Ω , where for all x ∈ Ω, the matrix M(x) is a metric of R n , namely a symmetric definite positive matrix of R n . As it is explained below, the notion of unit mesh with respect to a metric field establishes a link between a metric field and a mesh, which reduces the problem of finding an optimal mesh to the problem of finding an optimal metric field. First of all, the scalar product induced by a constant metric M is given by

x, y M = t x M y, for all x, y ∈ R n . (1)
The following distance results from this scalar product.

x M = √ t x M x, for all x ∈ R n . (2)
Likewise, in the Banach space R n , . M , the distance between two points a, b ∈ R n is given by

M (ab) = √ t ab M ab, (3)
where ab stands for the vector linking a to b. In the same way, if M = (M(x)) x∈Ω is a smooth metric field on Ω, the length of the segment between a and b can be computed by using the integral formula

M (ab) = 1 0 t ab M (γ(t)) ab dt, (4)
where γ(t) = (1 -t)a + tb.

As the next definition shows, there is a strong correspondence between continuous metric spaces and meshes, through the notion of unit mesh with respect to a metric field.

Definition 1.1 A mesh H of a domain Ω ⊂ R n is said to be unit with respect to a continuous metric field M = (M(x)) x∈Ω if the edges lengths of its elements equal 1. If {e1, ...en+1} are the edges of an element K of H, then M (ei) = 1, for all i ∈ {1, ..., n + 1} .

In most of the metric based mesh adaptation methods, the problem is reduced to an optimization problem of finding the optimal metric field for which unit meshes minimize the interpolation error of the solution. In this article, the problem of surface approximation is seen as a 2D high order metric based mesh adaptation in the parameters space.

Error model

The examined error model is the on studied in [START_REF] Coulaud | Very High Order Anisotropic Metric-Based Mesh Adaptation in 3D[END_REF]. Let u be a smooth solution on the domain Ω, H be a mesh of Ω and k be an arbitrary positive integer.

In what follows, Π k u denotes the projection of u onto the finite elements space P k (H), whose functions are polynomials of degree k on each element K of H. For all x0 ∈ Ω, it is well known that there exists a positive constant C such that, for all x0 ∈ Ω and

x ∈ R n , |u(x) -Π k u(x)| ≤ C d (k+1) u(x0)(x -x0) +O x -x0 k+2 , (5)
where d k+1 u(x0) is the differential form of u of order k + 1 at x0, |.| is the absolute value function and . denotes the Euclidean norm of R n . In the high order case, the main idea is to replace the right hand side of ( 5) by a term governed by a metric field Q = (Q(x)) x∈Ω , which approximates the k + 1 differential form of u. More precisely, we are looking for Q such that for all x0 ∈ Ω,

x ∈ R n d k+1 u(x0)(x -x0) ≤ t (x -x0) Q(x0) (x -x0) k+1 2
(6) The main issue is to find the metric field Q such that the inequality (6) is as optimal as possible. From a geometrical point of view, the local problem is to find the largest ellipse in 2D (or the largest ellipsoid in 3D) included into the domain surrounded by the level set of level 1 of d (k+1) u(x0). Indeed, let BQ = x ∈ R n : t x Q x ≤ 1 be the unit ball of a metric Q, which is an ellipse in 2D (an ellipsoid in 3D). Now assume that, for all x ∈ R n such that d

(k+1) u(x0)(x) = 1, one has t x Q x ≥ 1. Let x ∈ R n and y = x |d (k+1) u(x0)(x)| 1 k+1
. In particular d (k+1) u(x0)(y) = 1, and:

d (k+1) u(x0)(y) ≤ t y Q y. Since d (k+1) u(x0) is a homogeneous polynomial of de- gree k + 1, it comes 1 ≤ t x Q x |d (k+1) u(x0)(x)| 2 k+1
, and consequently

d (k+1) u(x0)(x) ≤ t x Q x k+1 2
, for all x ∈ R n .

The purpose of the next section is to solve this minimization problem.

Log-simplex method

The high order mesh adaptation method which we use in this paper is the log-simplex method, introduced in [START_REF] Coulaud | Very High Order Anisotropic Metric-Based Mesh Adaptation in 3D[END_REF]. It can be seen as an extension of the method described in [START_REF] Loseille | Continuous mesh framework part I: well-posed continuous interpolation error[END_REF][START_REF] Loseille | Continuous mesh framework part II: validations and applications[END_REF] to the high order case. The main difference between the P 1 and the P k adaptation methods relies on the fact that Q is directly given by the Hessian matrix of u when dealing with P 1 adaptation, whereas it is mandatory to find a suitable metric field satisfying (6) for the P k adaptation. The log-simplex algorithm is a way to compute such a metric field. It is based on a sequence of linear problems written in terms of the logarithm matrix L = log(Q). In this section, the highlights of this method are recalled. For a full description of the algorithm, see [START_REF] Coulaud | Very High Order Anisotropic Metric-Based Mesh Adaptation in 3D[END_REF]. Given a homogeneous polynomial p of degree k + 1 on R n which stands for d (k+1) u(x0), a set of points {x1, ..., xm} of R n such that p(xi) = 1 for all i ∈ {1, ..., m} is considered. The optimization problem that the log-simplex method solves is the following.

Find a metric Q such that det(Q) is minimal, t xi Q xi ≥ 1, for all i ∈ {1, ..., m} . (7)
The first line of ( 7) translates the fact that we are looking for the metric with the largest area (or volume in 3D). Since the cost function of this problem is nonlinear, one rewrites it as a problem in L = log(Q). Notice that L is not a metric but only a symmetric matrix. This formulation also allows the discrete counterpart of the problem to be well-posed. Indeed, in [START_REF] Coulaud | Very High Order Anisotropic Metric-Based Mesh Adaptation in 3D[END_REF], it is shown that the discrete form of ( 7) is ill-posed. For det(Q) = exp(trace(L)), a linear cost function is recovered by replacing Q by L in [START_REF] Laug | Some aspects of parametric surface meshing[END_REF]. On the contrary, the constraints which are linear on Q become nonlinear when writing them in terms of L. This can lead to really expensive computations.

To avoid this problem, the convexity property of the exponential is used and replaces these constraints by approximated linear ones. More precisely, through the classic convexity inequality, if

x ∈ R n satisfies t x L x ≥ -x 2 log x 2 , it ensures that t x Q x ≥ 1.
By this way, the following linear optimization problem is obtained.

Find a symmetric matrix L such that trace(L) is minimal, t xi L xi ≥ -xi 2 log xi 2 , ∀i ∈ {1, ..., m} . (8)
Since this problem is linear in L, it can be solved by a simplex method (see for instance [START_REF] Dantzig | Linear Programming 2: Theory and Extensions[END_REF]). Unfortunately, in most of the cases, solving (8) once does not provide accurate metrics, in the sense that the unit ball of the obtain metric Q = exp(L) can be far from the level set of p (see Figure 3). This issue is dealt by an iterative process. More precisely, once we have computed the solution of ( 8) and recovered Q = exp(L), we apply the mapping

x → Q 1 2 x by replacing p by q = p • Q -1 2 .
Then, we take a new set of points {x1, ..., xm} such that q(xi) = 1, for all i ∈ {1, ..., m} and solve again [START_REF] Siqueira | A Hierarchical Adaptive Mesh Generation Strategy for Parametric Surfaces Based on Tree Structures[END_REF]. Finally, the log-simplex algorithm is the following.

input : A mesh H of Ω d (k+1) u(x), for all x ∈ H output: Q = (Q(x)) x∈H foreach x ∈ H do repeat
choose a set of points {x1, ..., xn} on the level set of p of level 1 perform the log-simplex algorithm and obtain a metric

Q replace p by p • Q -1 2 until convergence; end
In order to implement numerically the log-simplex method, notice that this algorithm must contain a polynomial reduction so that the possible infinite branches in the level set of d (k+1) u disappear. All the theoretical and numerical issues of the log-simplex method are described in details in [START_REF] Coulaud | Very High Order Anisotropic Metric-Based Mesh Adaptation in 3D[END_REF]. The objective of the following section is to find the function u on which the high-order error model should be applied and then to deduce a metric on which mesh adaptation will be performed in the parameters space.

HIGH-ORDER SURFACE MESH GENERATION 2.1 Metrics for linear surface mesh generation

In the case of parametric surface meshing, the whole problem is to find a suitable metric thanks to which a mesh adaptation process in the parameters space will be performed. First let us have a look on the case of curve meshing.

Metrics for curve meshing

When dealing with meshing of parametric curves, it is frequent to perform a local analysis on it. To do so, let us have a look on a Taylor expansion of a parametric curve t → γ(t) ∈ R 3 , that we will assume smooth enough, in the vicinity of t01

γ(t) = γ(t0)+γ (t0)(t-t0)+ γ (t0) 2 (t-t0) 2 +O((t-t0) 3 ).
Now if we consider a change of variable with s being the curvilinear abscissa such that s(t0) = 0 and ds dt = ||γ (t)|| then the Taylor expansion becomes:

γ(s) = γ(0) + sT + 1 2 κ(t0)s 2 N + O(s 3 )
where T = γ (t) ||γ (t)|| and (κ, N) are such that dT ds = κN. (κ, T, N) are intrinsic data [START_REF] Carmo | Differential Geometry of Curves and Surfaces: Revised and Updated Second Edition[END_REF]. κ is called the curvature and can be computed with:

κ(t) = ||γ (t) × γ (t)|| ||γ (t)|| 3 .
By setting B = T×N, (T, N, B) defines an orthornormal basis, of so-called Frénet frame.

Note that if we denote (x, y, z) the components of γ(t) in the Frénet frame defined in t0 with γ(t0) = (x0, y0, z0) in this frame, we have:

   y = y0 + 1 2 κ(t0)(x -x0) 2 + O(|x -x0| 3 ) z = z0 + O(|x -x0| 3 ) (9)
Based on the Frénet frame, and on the curvature, a 3D metric tensor can be deduced via:

M1 = t T t N t B    1 (2 √ (2-)ρ(t)) 2 0 0 0 λ 0 0 0 λ      T N B   where λ ∈ R is an arbitrary constant 2 , ρ(t) = 1 κ(t)
is the radius of curvature, and 2 (2 -) is a scaling coefficient which guarantees, for a second order approximation of the curve, to maintain a deviation gap between the mesh elements and the curve geometry of [START_REF] Frey | About surface remeshing[END_REF]. As the metric relies on only intrinsic data, it is independent of the parameterization. The metric can be mapped back to the parameter space via the following formula:

M1 = t γ (t)M1γ (t).
In this case, the formula simplifies to M1 = 1

h 2 1 = ||γ (t)|| 2
(2

√ (2-)ρ(t)) 2 .
Once the metrics for curves are set, the next step is to define metrics for the surfaces.

Metrics for surface meshing

The meshing process of parametric surfaces is a bit more complex. It relies on some differential geometry notions [START_REF] Carmo | Differential Geometry of Curves and Surfaces: Revised and Updated Second Edition[END_REF]. For this purpose, let us consider a parametric surface (u, v) → σ(u, v) ∈ R 3 that we will assume smooth enough. In this case, the first fundamental form I(du, dv) is defined as follows:

I(du, dv) = (du dv) ||σu|| 2 (σu, σv) (σu, σv) ||σv|| 2 du dv ,
where (du, dv) is an elementary displacement, and σu (resp. σv) the partial derivative of σ w.r.t. u (resp. v). The first fundamental form explains how the three dimensional distances are perceived in the two dimensional space. In particular, it provides a two-dimensional riemannian structure to the surface with a metric tensor defined as:

MI = ||σu|| 2 (σu, σv) (σu, σv) ||σv|| 2 .
In the same framework, comes also the second fundamental form II(du, dv) that is defined as follows:

II(du, dv) = (du dv) (σuu, N) (σuv, N) (σuv, N) (σvv, N) du dv ,
where (σuu, σuv, σvv) are the second derivatives of σ w.r.t (u, v) and N = σu×σv ||σu×σv || is the normal vector to the surface. The second fundamental form expresses the gap of a surface to its tangent plane at the order two.

Based on these two quadratic forms and their matrices, we are able, for a given point of the surface, to define the principal curvatures (κi)i=1,2 and principal directions (Vi)i=1,2 (in 3D) as solution of the generalized eigenvalue problem:

   MII vi = κiMI vi Vi = (σu σv) vi || (σu σv) vi|| i = 1, 2
with MII being the symmetric matrix associated to the second fundamental form. These quantities are independent of the parameterization and when κ1 = κ2, (V1, V2) forms an orthonormal basis of the tangent plane. If we complete the basis with N, they form a local basis (V1, V2, N) of R 3 . Note that if we denote (x, y, z) the components of σ(u, v) in this local basis defined in (u0, v0) with σ(u0, v0) = (x0, y0, z0) in this basis, we have:

z = z0 + 1 2 (κ1(u0, v0)(x -x0) 2 + κ2(u0, v0)(y -y0) 2 )
+ O( (x -x0, y -y0) 3 ). ( 10) Now, thanks to this basis and on the curvatures, we can define the following 3D metric tensor:

M2 = t V1 t V2 t N ×   1 (c 1 ρ 1 (u,v)) 2 0 0 0 1 (c 2 ρ 2 (u,v)) 2 0 0 0 λ     V1 V2 N   where λ ∈ R is an arbitrary constant, ρi(u, v) = 1 κ i (u,v)
for i = 1, 2 are the radii of curvature, with the convention ρ1(u, v) ≤ ρ2(u, v) and c1and c2 are scaling coefficients. For the direction of greater curvature (e.g. the direction given by V1), we want to control the deviation under a threshold which comes down to set c1 to the value of 2 (2 -). Now, as we want the same threshold in all the directions in the tangent plane, the coefficient c2 is set to 2 [START_REF] George | Meshing[END_REF]. Similarly to curves, the metric relies only on intrinsic data and is therefore independent of the parameterization Now, the metric M2 can be mapped back to the parameters space by applying the first fundamental form:

ρ 1 ρ 2 (2 -ρ 1 ρ 2 )
M2 = t σu t σv M2 σu σv .
This is the metric that will be used as an anisotropic metric for the mesh adaptation in the parameters space.

Computation of higher-order metrics

The object of this section is to deal with the extension of the previous framework for higher-order elements.

In particular, we seek for parameterization independent Taylor expansion similar to [START_REF] De Siqueira | An Adaptive Parametric Surface Mesh Generation Method Guided by Curvatures[END_REF] and to [START_REF] Dey | A robust conforming NURBS tessellation for industrial applications based on a mesh generation approach[END_REF] with terms of degree greater than 2. First let us have a look on the case of the curve.

Case of the curve

As seen previously, the metric should rely on intrinsic data to be independent of the parameterization. A way to do so is to have a look at the formula [START_REF] De Siqueira | An Adaptive Parametric Surface Mesh Generation Method Guided by Curvatures[END_REF]. This formula gives a Taylor expansion of the gap of a curve to the straight edge at the order two and shows that it is driven by the curvature. Moreover, this expansion is done with the physical coordinates which naturally gives an independence with respect to any parameterization. A natural idea is therefore to extend the previous Taylor expansion to get higher-order terms and deduce metrics that will be fitted to high-order approximation.

To do so, let us write γ(t) in the Frénet frame (T0, N0, B0) associated to t0, a regular point of γ. If we note X = x -x0, Y = y -y0, Z = z -z0, we have:

     X = (γ(t) -γ(t0), T0), Y = (γ(t) -γ(t0), N0), Z = (γ(t) -γ(t0), B0).
Let us note φ(t) = (γ(t) -γ(t0), T0). If t0 is a regular point of γ then γ (t0) = 0 and therefore φ (t0) = 0. The inversion function theorem can thus be applied and there exists a function ψ such that

ψ(X) = ψ(φ(t)) = t -t0 in the vicinity of t0. More- over, if φ is C k+1 then ψ is C k+1 and ψ (X) = 1 φ (t)
with X = φ(t).

Based on this statement, it is thus possible to get a Taylor expansion of t -t0 with respect to X up to order k. To do so, let us compute the higher-order derivatives of ψ in t0. As γ ( and consequently φ) is an analytical function issued from CAD model, all its derivatives can be computed using the implementation details of [START_REF] Piegl | The NURBS Book[END_REF]. The derivatives of ψ are then deduced using the following result [START_REF] Di | Note sur une nouvelle formule de calcul différentiel[END_REF]:

Theorem 2.1 (Faà di Bruno's Formula) Let us consider f, g : R → R of class C k+1 with k + 1 ≥ n, then

d n dx n (g(f (x))) = E n! m1!..mn! g (m 1 +...+mn) (f (x)) n j=1 f (j) (x) j! m j
,

where E = {(m1, .., mn) ∈ N n | n k=1 k.m k = n} 3 .
In our case, we set g = ψ and f = φ and it comes that

d n dx n (g(f (x))) = 0 for k + 1 ≥ n ≥ 2.
This brings us to that for k + 1 ≥ n ≥ 2:

(φ (t)) n ψ (n) (x) = F (ψ (x), .., ψ (n-1) (x), φ (t), .., φ (n) (t)),
where F is a function that can be directly deduced from Theorem 2.1. This result means that as long as φ (t) = 0, then ψ (n) (x) can be recursively computed given its previous derivatives and the derivatives of φ.

Thanks to this, we can now write:

t = t0 + k+1 n=1 anX n + O(|X| k+2 ),
where an have been computed with the derivatives of ψ. Now if we recall that γ is C k+1 , we also have:

             Y = k+1 n=1 bn(t -t0) n + O(|t -t0| k+2 ), Z = k+1 n=1 cn(t -t0) n + O(|t -t0| k+2 ).
By composition of both Taylor expansion, we then obtain a Taylor expansion of Y and Z in X at the order k + 1 ≥ 2, which is independent of the parameterization and a generalization of [START_REF] De Siqueira | An Adaptive Parametric Surface Mesh Generation Method Guided by Curvatures[END_REF]. From this, we have an intrinsic information of the gap of a curve to the straight edge up to the order k + 1. We can then write:

Y Z = F k (X) + A k+1 X k+1 + O(|X| k+2 ),
where F k (X) is a polynomial of degree k in X and

A k+1 ∈ R 2 .
For an approximation of a curve at the degree k, the leading term of the error is therefore A k+1 X k+1 [START_REF] Loseille | Continuous mesh framework part I: well-posed continuous interpolation error[END_REF][START_REF] Loseille | Continuous mesh framework part II: validations and applications[END_REF][START_REF] Coulaud | Very High Order Anisotropic Metric-Based Mesh Adaptation in 3D[END_REF]. Thus, if we control ||A k+1 X k+1 ||, the error of approximation will be controlled. Now, if we note that:

||A k+1 X k+1 || = (||A k+1 || 2 k+1 X 2 ) k+1 2 ,
we can then set κ k+1 = 2||A k+1 || 2 k+1 and reuse the metrics used for the linear meshing with ρ = 1 κ k+1 for radii of curvature and k+1 = 2 k+1 for threshold. This way, the classic formula is found for k = 1 and a generalization is proposed for k ≥ 2. Note that in some configurations, the found size for the order k can be significantly lower than the found size for the order k + 1. In this case the size given by order k + 1 is preferred. Now , let us interest to the more complex case of the surfaces.

Case of the surface

Like for curves, the idea is to start from an intrinsic representation of the surface. For this purpose, let us have a look at formula [START_REF] Dey | A robust conforming NURBS tessellation for industrial applications based on a mesh generation approach[END_REF]. This formula gives a Taylor expansion of the gap of a surface to its tangent plane in physical coordinates which is an intrinsic representation. As previously, our idea is to extend this Taylor expansion to higher-order terms. To do so, let us note (V1,0, V2,0, N0) the local basis defined for a regular point (u0, v0) of the surface. If we note X = x -x0, Y = y -y0, Z = z -z0, we have:

     X = (σ(u, v) -σ(u0, v0), V1,0), Y = (σ(u, v) -σ(u0, v0), V2,0), Z = (σ(u, v) -σ(u0, v0), N0). Now, if we define: Φ(u, v) = (σ(u, v) -σ(u0, v0), V1,0) (σ(u, v) -σ(u0, v0), V2,0) ,
and if (u0, v0) is a regular point of σ, then its jacobian matrix JΦ is invertible in (u0, v0). Consequently, the inverse function theorem can be applied and there ex-

ists a function Ψ such that Ψ(X, Y ) = Ψ(Φ(u, v)) = (u -u0, v -v0) in the vicinity of (u0, v0). Moreover, if Φ is C k+1 then Ψ is C k+1 and JΨ(X, Y ) = JΦ(u, v) -1 with (X, Y ) = Φ(u, v).
With this statement, we know that we can have a Taylor expansion of (u -u0, v -v0) with respect to (X, Y ) up to order k + 1. For this purpose, let us compute the higher-order derivatives of Ψ. As σ (and therefor Φ) is an analytical function issued from a CAD model, all its derivatives can be computed using the recipes in [START_REF] Piegl | The NURBS Book[END_REF]. The derivatives of Ψ can be then deduced using the following result [START_REF] Encinas | A short proof of the generalized Faà di Bruno's formula[END_REF]:

Theorem 2.2 (2D Faà di Bruno's Formula) Let us consider f, g : R 2 → R 2 of class C k+1 with k + 1 ≥ |α|, then ∂ |α| ∂x α (g(f (x))) = n |σ|=1 α! ∂ |σ| ∂x σ (g(f (x)))× Eσ 2 i=1 Aα 1 e iα i ! 1 α i ! ∂ |α i | ∂x α i (fi(x)) e iα i , where f (x) = (f1(x), f2(x)), α = (α1, α2) ∈ N 2 ,n = |α| = α1 + α2, σ = (σ1, σ2) ∈ N 2 Eσ =    (e 1α 1 , e 2α 2 ) ∈ N 2 , 1 ≤ |α i | ≤ n,   n |α i |=1 e iα i = σi,   i=1,2    and 4 Aα =    (α 1 , α 2 ) : 1 ≤ |α i | ≤ n, i = 1, 2, 2 i=1 n |α i |=1 e iα i .α i = α    .
In our case, we set g = Ψ and f = Φ and it comes that

∂ |α| ∂x α (g(f (x))) = 0 for k + 1 ≥ |α| ≥ 2.
4 Note that in practice, the sets Eσ and Aα can be precomputed once for all for the range of values of α and σ that are used.

If we consider all the n + 1 α such that |α| = n, then we have a system of equations of the kind:

A ∂ |α| Φ ∂x α |α|=n × ∂ |α| Ψ ∂x α |α|=n = F ∂ |α| Ψ ∂x α |α|<n , ∂ |α| Φ ∂x α |α|≤n ,
where A is (n + 1) × (n + 1) matrix, ∂ |α| Ψ ∂x α |α|=n is a vector of size n + 1 containing the n + 1 derivatives of Ψ of order n and F is vector function of size n + 1 that can be deduced from theorem 2.2. Moreover, it is shown in [START_REF] Encinas | A short proof of the generalized Faà di Bruno's formula[END_REF] that |A| = |JΦ| n , which proves that the system has always a solution if the inverse function theorem is successfully applied. This way, a recursive method to compute all the derivative of Ψ in (u0, v0) is set and the computation of the Taylor expansion is therefore possible:

u -u0 v -v0 = k+1 n=1 i+j=n A n ij X i Y j + O(||(X, Y )|| k+2 ),
where A n ij ∈ R 2 and is defined thanks to the partial derivatives of Ψ. But, as σ is C k+1 , we also have:

Z = k+1 n=1 i+j=n c n ij (u -u0) i (v -v0) j + O(||((u -u0), (v -v0))|| k+2 ).
By composition of both Taylor expansions, we then obtain a Taylor expansion of Z in (X, Y ) at the order k+1 ≥ 2 which is independent of the parameterization and a generalization of the formula [START_REF] Dey | A robust conforming NURBS tessellation for industrial applications based on a mesh generation approach[END_REF].

The gap to the tangent plane is thus expressed up to the order k + 1:

Z = F k (X, Y ) + R k+1 (X, Y ) + O(||(X, Y )|| k+2 ), (11)
where F k is a polynomial of degree k and R k+1 is an homogeneous polynomial of degree k + 1. For an approximation of the surface at the degree k, the leading term of the error is R k+1 (X, Y ) [START_REF] Loseille | Continuous mesh framework part I: well-posed continuous interpolation error[END_REF][START_REF] Loseille | Continuous mesh framework part II: validations and applications[END_REF][START_REF] Coulaud | Very High Order Anisotropic Metric-Based Mesh Adaptation in 3D[END_REF]. So, if we want to control the P k approximation, we need to control |R k+1 (X, Y )|. By applying the log-simplex algorithm explained in the previous section, we are able to find a metric that satisfies an inequality like [START_REF] Wang | EQSM: An efficient high quality surface grid generation method based on remeshing[END_REF], that is to say, we are able to compute a matrix Q k+1 such that:

|R k+1 (X, Y )| ≤ 1 2 (X Y )Q k+1 X Y k+1 2
, where Q k+1 is the optimal symmetric matrix (in a sense explained in the first section) that verifies this inequality.

If we note (κ i,k+1 , v i,k+1 )i=1,2, the eigenvalues and eigenvectors of Q k+1 , we can then reuse the metrics used for the linear meshing with ρi = 1 κ i,k+1 the radii of curvature, V i,k+1 = (V1,0V2,0) v i,k+1 the principal directions in the tangent plane and k+1 = 2 k+1 the threshold. This way, the classic formula is found for k = 1 and a generalisation is proposed for k ≥ 2. Like for curves, in some configurations, found sizes for the order k can be significantly lower than found sizes for the order k + 1. In this case sizes given by order k + 1 are preferred.

Meshing process

The mesh generation process is based on the classical unit-mesh concept, where a metric field, as M2 or M2, is used to drive the orientation and sizing of the elements. In the context of parametric surface meshing, several approaches are typically devised to generate a final 3D surface mesh. Full 2D methods are a convenient way to avoid 3D surface meshing and inverse projection to the geometry. However, a special care is needed to handle degenerated points, periodicity, highly non uniform (even discontinuous) parameterization or degenerated edges. Approaches that mix 2D and 3D methods tend to reduce the impact of the parametric space to the final mesh.

In the paper, we consider a rather classical approach. The core steps of the procedure are decomposed as follows:

1. For each Edge The EGADS API [START_REF] Haimes | On The Construction of Aircraft Conceptual Geometry for High-Fidelity Analysis and Design[END_REF] is used to perform the CAD linking and the planar mesh generation process is performed using a Delaunay triangulation-based algorithm [START_REF] George | Delaunay triangulation and meshing[END_REF].

EXAMPLES

We illustrate our approach on several examples. To compute error indicates, the computation of an exact Haussdorf distance would be too expansive. Therefore, a normalized and an absolute error in distance are computed. For each triangle of the generated surface mesh, a sampling of the underlying parameters space is performed and the distance between a point of the curved mesh and its surface counterpart is performed:

δ(ũ, ṽ, w) = ||σ( n i=0 φi(ũ, ṽ, w)(ui, vi)) - n i=0
φi(ũ, ṽ, w)σ(ui, vi)||, where (ui, vi) is the set of points of the geometry used to define the high-order triangle φi the associated shape functions and (ũ, ṽ, w) are the barycentric coordinates of the triangle.

To get the normalized error in distance, the absolute distance is normalized by multiplying it by the local curvature of greater value.

The worst of all these errors then gives the overall errors of the meshing process.

Note that the relative error can be higher than the absolute one. This process can be explained for several reasons. Firstly, the absolute error, is dependent on the system of units. In the case of the shuttle, the average order of magnitude for the bounding box is 10 2 -10 3 while it is of 10 0 -10 1 in the case of the torus. Secondly, the point of worst value that determines the error is not necessarily the same for both errors. The worst value of the absolute error can occur in high stretching zones where the curvature is high, which emphasizes the value. Thus, it provides levels of error of greater magnitude. Note that all the figures are obtained using Vizir [START_REF] Loseille | Vizir: High-order mesh and solution visualization using OpenGL 4.0 graphic pipeline[END_REF].

Sphere example. Despite its simplicity, the unit sphere example is used to illustrate that our error estimate is independent to the parameterization of the model. A classic boundary representation of a sphere is to consider a surface of revolution, where the two pole maps to two degenerated edges. In the vicinity of these points, depending on the underlying CAD kernel, the definition of normals or more generally the definition of the derivatives of u, v are either undefined or unstable.

We illustrate the obtained Taylor expansion on the point of the sphere using the inversion formula given by [START_REF] Borouchaki | Maillage de surfaces paramétriques. Partie I: Aspects théoriques[END_REF]. Note that the sphere is parameterized as a circle of revolution. far from the pole, we have :

Z = -0.5Y 2 -0.5X 2 -0.125Y 4 -0.25X 2 Y 2 -0.125X 4 -0.0625Y 6 -0.1875X 2 Y 4 -0.1875X 4 Y 2 -0.0625X 6 +O(||(X, Y )|| 7 ).
In the vicinity of the two poles, the expansion is

Z = -0.5Y 2 -0.5X 2 -0.125Y 4 -0.25X 2 Y 2 -0.125X 4 -1.80444e -9 XY 4 + 1.9886e -9 X 3 Y 2 -0.0625Y 6 -0.187501X 2 Y 4 -0.187499X 4 X 2 -0.0625X 6 +O(||(X, Y )|| 7 ).
We observe that a numerical noise appears while an almost perfect expansion is obtained as in regular points.

As expected, the second order terms reveal half the principal curvatures. This numerical noise is a consequence of the manipulation of the huge values at stake in the vicinity of the apex. As the value of the derivative of the parameters is going to 0, the normalization can lead to arbitrary huge values. In the end, everything is simplified but the numerical computations are a bit impacted by this.

Torus. The torus example is based on a NURBS representation is based on two NURBS of degree 5 both defined by 6 control points and 12 knots. P 1 , P 2 and P 3 meshes are generated for this case. The point-wise normalized and absolute error in distance are reported in Table 1. The meshes and point-wise errors are reported in Figures 4 and5 Also, to show the efficiency of the method on this case, a study has been made between a method where the P 1 metric will be used and ours where a P k metric is used to generate a meh of degree k. For various absolute error levels, P 1 -metric based and P k -metric based meshes are generated and comparison are done with regards to the number of needed degrees of freedom.

In Table 2 (resp 3), a comparison between a P 1 -metric and a P 2 -metric (resp. P 3 -metric) is done. It shows in particular, that for a given error, the number of requested degrees of freedom is higher using the P 1error than the P 2 (resp. P 3 ) one. If we have a look at the meshes with the point-wise error in Figure 6 (resp 7), we explain it by the fact that the high-order metric is able to detect high-order features that the P 1 -metric does not detect. In particular, we see that the error due to the presence of an annulus on the torus is mostly canceled. And from a more global point a view, the high-order metric provides a mesh with a more equitable distribution of the error with a lower number of degrees of freedom.

Error

CONCLUSION

Generating high-order curved surface meshes from a geometry require the derivation of intrinsic quantities of the surface in order to guarantee the approximation. In the paper, we have used a Taylor expansion coupled with an inversion formula to derive a local approximation of the underlying surface in the Frénet frame. To extend the notion of principal curvature, a log-simplex approach is used to approximate optimally the variation of the polynomial by a quadratic function. The eigenvalues and eigenvectors can be viewed as "high-order" curvatures. As a metric field is naturally derived, these estimates can be used directly in any adaptive anisotropic mesh generation process.

Natural perspectives are possible for this work. Notably, this work could be extended to the implicit surface meshing. Indeed, once the implicit function theorem is applied, it is possible to apply the previous framework in order to get a Taylor expansion in any point of the surface and consequently deduce a metric. The only limitation to this meshing method would be to find the set of points where the surface is evaluated, which is the same problem as in linear meshing of implicit surfaces.

Also, the parametric mesh generation process could be improved by performing metric based curvilinear mesh adaptation to optimize the position of the high-order nodes. This problem is currently under investigation.

Figure 1 :

 1 Figure 1: Topology hierarchy (Face, Loop, Edge) of the continuous representation of model using the Boundary REPresentation (BREP).

Figure 2 :

 2 Figure 2: Mapping between the parametric and physical space.

Finally

  examples are shown to highlight the process.

Figure 3 :

 3 Figure 3: Illustration of the log approximation for the constraints for an error level 1 (in grey). The optimal metric (in red) is far from the boundary of the error due to the convexity approximation.

Figure 4 :

 4 Figure 4: P 1 (top), P 2 (middle) and P 3 (bottom) meshes for the torus geometry.

Figure 5 :

 5 Figure 5: P 1 (top), P 2 (middle) and P 3 (bottom) point-wise distance to the torus geometry.

Figure 6 :

 6 Figure 6: P 2 meshes generated with a P 1 -metric (left) and P 2 -metric (right) for various error levels: 0.5 (top), 0.2 (middle), 0.03 (bottom)

Figure 7 :

 7 Figure 7: P 3 meshes generated with a P 1 -metric (left) and P 3 -metric (right) for various error levels: 0.1 (top), 0.002 (middle), 0.0003 (bottom)

Figure 8 :

 8 Figure 8: P 1 (top), P 2 (middle) and P 3 (bottom) meshes for the shuttle geometry.

Figure 9 :

 9 Figure 9: P 1 (top), P 2 (middle) and P 3 (bottom) point-wise distance to the shuttle geometry.

Table 1 :

 1 . Deviation to the geometry for torus geometry for P 1 to P 3 meshes.

	Order DOF Normalized	Absolute
	P 1	2787	0.0371489	0.275027
	P 2	666	0.0287625	0.18465
	P 3	252	0.0137922	0.0883003

Table 2 :

 2 Comparison between the number of needed degrees of freedom to get a given error between a P 1 -metric and a P 2 -metric.

		P 1 -metric #DOF P 2 -metric #DOF
	0.5	347	156
	0.2	687	242
	0.03	1556	1036
	Error	P 1 -metric #DOF P 3 -metric #DOF
	0.1	255	238
	0.002	6171	1848
	0.0002	11329	5362

Table 3 :

 3 Comparison between the number of needed degrees of freedom to get a given error between a P 1 -metric and a P 3 -metric.Shuttle. The shuttle geometry is based on two NURBS of degree 3 defined by 8 (resp. 13) control points and 12 (resp. 17) knots with strong variation in the parametric space. The P 1 , P 2 and P 3 meshes are depicted in Figures8 and 9. The error to the geometry are reported in Table4.

	Order DOF Normalized Absolute
	P 1	1647	254.551	197.814
	P 2	1690	61.6727	58.3507
	P 3	1791	7.02513	22.217

Table 4 :

 4 Deviation to the geometry for the shuttle geometry.

In the case of NURBS, it is easy to split the NURBS according the smooth parts[START_REF] Piegl | The NURBS Book[END_REF].

Note that in practice, the set E can be precomputed once for all for the range of values of n that are used.

ACKNOWLEDGMENTS

This work was supported by the following public grants : Investissement d'avenir project, reference ANR-11-LABX-0056-LMH, LabEx LMH and ANR IMPACTS, reference ANR-18-CE46-0003. The authors are also grateful to Bob Haimes (MIT) for giving some advice about the use of the EGADS API for high-order differentiation.


Michael A Sutton

François Hild

Helena Jin

Xia Li

Michel Grediac

published or not. The documents may come   L'archive ouverte pluridisciplinaire


L Martinelli

D J Young

S Gossé

S Bosonnet

Corrosion of 316L in Liquid Tellurium at 551°C

Keywords: Tellurium corrosion, telluride scale, dissolution, diffusion control

Grade 316L stainless steel was exposed to liquid tellurium at 551 o C. Corrosion was rapid, leading to more than 100 m loss of steel section in 30 min. The corrosion product was a mixed telluride scale, which thickened according to parabolic kinetics whilst simultaneously dissolving at its outer surface into the liquid tellurium. A mathematical model based on diffusion-controlled scale growth coupled with dissolution at a rate controlled by liquid phase diffusion is shown to describe scaling and dissolution kinetics successfully.

Introduction

Stainless steel grade 316L is often considered as fuel cladding in fast nuclear reactors, and its resistance to corrosion by fission products can be important. One such product is tellurium [1], which melts at 450 o C, and can be liquid at operating temperatures at the cladding surface. As caesium is also a fission product [1], it will reduce any oxide scale on the cladding internal surface. Thus direct reaction between metal and Te(l) can be of practical concern.

The corrosion process is also of fundamental interest, as it involves formation of solid metal tellurides [1]- [6] and their dissolution into the melt. This paper reports the results of a preliminary laboratory investigation of the kinetics of 316L corrosion in liquid tellurium, and the morphological evolution of the reaction products.

Thermodynamic considerations

The phase diagrams and thermodynamic interactions between the major alloying elements (Cr, Fe, Ni) and Te were studied by numerous authors: Cr-Te [7] [8] [9], Fe-Te [10]- [START_REF] Arvhult | Part II: Thermodynamic modelling[END_REF] and Ni-Te [START_REF] Ettenberg | [END_REF] [15]. Thermodynamic informationmetal solubilities in liquid Te, phases at equilibrium, standard Gibbs free energies of formation for the intermetallic telluridesis of major importance when dealing with corrosion issues. All these data made it possible to perform a thermodynamic assessment of these binary metal-chalcogen systems using the Calphad method [16] based on these thermodynamic and phase diagram data. These metal-tellurium binary and ternary systems are under development in order to obtain a thermodynamic database dedicated to fuelcladding interaction applications [START_REF] Arvhult | Thermodynamic Assessment of the Fe-Te System[END_REF] [START_REF] Arvhult | Part II: Thermodynamic modelling[END_REF].

From these thermodynamic modelling operations, the liquidus line of each binary system was calculated. In Figure 1, the metal solubility limits were calculated as a function of temperature inverse. At 823 K (10000/T=12.15 K -1 ), liquid Te has solubility limits of 5.33 at. % for iron [START_REF] Arvhult | Part II: Thermodynamic modelling[END_REF], 4.55 at. % for nickel and 2.54 at. % for chromium. From Calphad modelling, as shown in Figure 1, the solubility limit of iron, chromium and nickel in liquid Te is given by the following relations:

for 720<T(K)<919 for 919<T(K)<1039 for 1039<T(K)<1088 (2)
for 720<T(K)<1160

(3) These high solubilities highlight the possible extensive loss of steel depending on the amount of liquid tellurium, and the dissolution rate. The high temperature loops correspond to the congruent melting of the Te-rich binary intermetallics.

Experimental

Corrosion tests were performed using 316L stainless steel (composition in Table 1) ground with SiC paper to a 1200 grit finish before the experiments. Samples of dimensions 30x10x0.5 mm, were immersed one by one in a crucible containing 107 g of liquid Te at 824 K. The volume of liquid Te was 18.8 cm 3 in a pyrex crucible of 5.3 cm diameter, producing a (1) liquid height of 0.85 cm. A diagram of the facility is presented in Figure 2. Samples were fixed by a Mo wire to a pyrex sample holder and immersed in the liquid Te, above and approximately parallel to the bottom of the crucible. The crucible was contained in a glass reactor which was swept continuously with Argon 6.0 to maintain an inert environment. The reactor was contained in a glove box to avoid contamination from the atmosphere. Five samples were immersed, one after the other in the same melt, for various durations. In the event of a dissolution process, the liquid Te is expected to become continuously enriched in corrosion products. Thus the concentrations of Fe, Cr, Ni increase as a function of accumulated steel immersion time in the liquid. The corrosion rate is then expected to decrease as the liquid is enriching in dissolved Fe, Cr, Ni. The first sample was immersed for 1 min and the second for 10 min, leading in the latter case to total dissolution of that sample. The third sample was immersed also for 10 min, the fourth 30 min and the fifth 5 min. The sequence of immersion experiments and their duration are shown in Figure 3.

Results

Successive steel specimens were immersed in and withdrawn from the same melt in the sequence shown in Figure 3. Apart from the second specimen, which had completely dissolved, all specimens emerged from the melt covered with a thin film of liquid. Crosssections of the remaining steel together with the adhering solidified melt are shown in Figure 4. Solidification has led to precipitation from the parent melt of particles, the volume fraction of which increases with accumulated immersion time. Simultaneously, an adherent reaction product scale grows on the steel surface. Compositional analysis of the various phases by EDX produced the qualitative maps of Figure 5 and Figure 6, and the quantitative profiles of Figure 7. The particles precipitated during solidification of the melt are iron or chromium tellurides, of approximate composition (Fe,Ni,Cr)Te 2 . X-ray diffraction analysis of the solidified melt on a stainless steel specimen (shown in Figure 8) confirmed the presence of the phase FeTe 2 . From the Calphad modelling (Figure 9), the calculated Te rich domain of the Cr-Fe-Te isothermal section at 823 K confirms the presence of a three phase domain where liquid Te equilibrates with two different intermetallics: Te liq + Cr 3 Te 4 + FeTe 2 . In the present study the phase Cr 3 Te 4 could not be detected, neither with EDX (Figure 7) nor with XRD (Figure 8). Indeed, examination of the reaction product layer at the steel surface shows it to be continuous and apparently dense (Figure 5, Figure 6, Figure 10). Its composition is seen from the profiles in Figure 7 to correspond, on average, to metal telluride, MTe 2 . Within the layer are discrete regions of chromium enrichment, indicating the presence of a multi-phase telluride mixture. No depletion profiles are evident in the steel beneath its corroded surface, indicating that no preferential reaction of any element occurred. This is confirmed by a comparison of steel and scale compositions: the former has an Fe/Ni ratio of 7.5 and an Fe/Cr ratio of 3.7, whilst the scale has corresponding ratios measured as 8.4 and 3.7 (see Figure 7).

The scale continues to thicken with time. As seen in Figure 11, scale growth kinetics follow approximately parabolic kinetics (4) where X is the scale thickness grown in time, t, and k p the parabolic rate constant. It is noted that each point in Figure 11 represents a different specimen, each inserted with a clean, bare surface into the melt. Consequently the specimen immersed for 5 minutes (Experiment E, Figure 3) is the second point of the graph: although the specimen was immersed for only 5 minutes, the liquid Te was extensively enriched in dissolved Fe as a result of previous experiments (Figure 3).

Overall corrosion of the steel was assessed by measuring the thickness of the remaining steel section and calculating section loss for each specimen after a single immersion. The kinetics of this process are seen in Figure 12 to be neither linear nor parabolic. It is also clear that the rate of steel consumption is far higher than the rate of scale growth, consistent with the appearance of increasing amounts of Fe, Cr, Ni and Mo in the melt, as evidenced by the growing volume fraction of MTe 2 precipitates in the solidified liquid (Figure 5, Figure 6, Figure 10).

Figure 12: Steel corrosion kinetics

The amount of metal dissolved into the melt as steel is consumed is readily calculated from simple mass balance, using the measured loss of steel thickness. The resulting average concentrations are then found by taking into account the relative volumes of melt and steel consumed according to the following relation:

(5)

Here

is the average dissolved iron concentration in the liquid Te, µ 316L and µ Te are respectively the density of solid iron and liquid tellurium (g cm -3 ), M Fe is the molar mass of Fe (g mol -1 ), m Te is the mass of Te in the crucible (g), S is the sample surface area (cm 2 ), X 316 is the loss of thickness of 316L of one face of the sample (cm), and is the weight fraction of iron in 316. Edge effects have been ignored.

Increments in dissolved iron concentration can be calculated in this way from the specimen consumption observed in each single exposure experiment. In the first such exposure, the initial solute iron concentration was zero, and the calculated increment is therefore equal to the final concentration. In subsequent exposures, of course, it is not. For this reason, and because of the unusual sequence of exposures (see Figure 3), the results are shown in the bar chart of Figure 13 as concentration increments for the succession of corrosion processes specified in Figure 3.

Figure 13: Iron dissolution kinetics for exposures A to E (defined in Figure 3) as calculated from mass balance, and predicted from interface reaction model (Equation 11) and from diffusion model (Equation 12).

Discussion

Corrosion of 316L in Te(l) produces a thin surface layer of metal tellurides and extensive dissolution into the melt. During this process, the telluride-steel interface remains flat over the period observed (Figure 10).

On this basis, it is concluded that liquid metal is not present at the steel surface, and that dissolution takes place only at the outer surface of the telluride layer. If, conversely, the melt penetrated a porous telluride layer, the rapid dissolution of steel at penetration locations would lead to an irregular steel-telluride interface, an outcome not arrived at under the reaction conditions employed.

A solid reaction product layer develops on the steel surface, and thickens with time (Figure 4 to Figure 7, Figure 11). To investigate whether this is a steady-state situation, it is desirable to map the locus of composition observed across the steeltelluride layer system onto the appropriate phase diagram for the reaction system, to see if it is consistent with local equilibrium.

Unfortunately, no ternary phase diagrams are available and the thermodynamic data is limited in the Fe-Ni-Te ternary system. In the binaries, both Fe and Ni form the nonstoichiometric phases FeTe 2 and NiTe 2 at 823 K. The situation for the Cr-Te system is more complex, with a series of tellurides thought to be stable at 823 K, with the highest being Cr 3 Te 4 . As Ni is a minority element in both the steel and in the analysed telluride reaction product layer (Figure 7), it is ignored here, and the chemical system is approximated as Cr-Fe-Te.

The isothermal section shown in Figure 9 has been constructed on the basis of the binary diagrams. For simplicity, the Te rich side of the Cr-Te system is approximated as forming a single intermetallic Cr 3 Te 4 phase. The reciprocal Fe and Cr solubilities in Cr 3 Te 4 and FeTe 2 are unknown and not considered in the phase diagram. Nevertheless, due to the different structures of these intermetallics (Cr 3 S 4 for Cr 3 Te 4 [9] and FeS 2 pyrite prototype for FeTe 2 [10]), small reciprocal Cr and Fe solubilities may exist, but full intersolubility is regarded as unlikely.

The locus of composition in the solid part of the reacting system taken from the profiles in Figure 7 is mapped onto the approximate Fe-Cr-Te ternary diagram of Figure 9. It is seen to be consistent with local equilibrium, corresponding to a diffusion path. In the development which follows, mass transport in the solid telluride layer is taken to be diffusion controlled. The processes of telluride dissolution into the melt and transport of solute away from the telluride-melt interface, which represent the major part of the corrosion reaction are considered first.

Dissolution Reaction

The telluride product layer is closely similar in Fe-Ni-Cr content to the substrate steel. Thus the Fe/Ni ratio in the layer is found to be 8.4 compared with 7.5 in the steel, and the Fe/Cr ratio is 3.7 in both (see Figure 7). Since, moreover, little depletion is apparent in the subsurface alloy, it is concluded that the proportions of Fe, Ni and Cr solutes in the melt are the same as in the steel. Dissolution of the majority component, Fe, is considered and the flux of iron across the solid-liquid interface, , is evaluated.

In the simplest case, the iron content of the melt has no influence on dissolution rate, and = k (6) with k a constant. Then if mass transfer within the liquid phase is rapid in comparison to the dissolution rate (7) with the average concentration of Fe in the melt, and its original value. The prediction of linear kinetics for iron dissolution is clearly incorrect (Figure 12), and a more realistic description is sought.

Assuming now that the driving force for iron dissolution is the difference between solute iron activity, a Fe , and its value at equilibrium, , then

with k r the rate constant for the phase boundary mass transfer process. If it is also assumed that transport within the melt is rapid, and the value of a Fe is essentially uniform in the melt, then Eqn (8) yields (9) where is the activity coefficient for Fe in the melt, the iron concentration at equilibrium in liquid Te (assumed equal to iron solubility limit given in Eqn ( 4))and G a geometric factor relating the surface area of the specimen to the volume of liquid:

Here S is the sample surface area (cm 2 ), R is the crucible radius (cm) and X Te is the depth (cm) of liquid Te in the crucible.

Integration of Eqn (9) using the initial condition at t 0 (10) leads to the result (11) where .

This equation for interfacial reaction controlled dissolution can be tested by applying it to each of the successive stages of the experiment charted in Figure 3, using the indicated, successively higher values of as initial values in calculating incremental increases in dissolved iron concentration. The value of s -1 was arrived at in optimising the fit of Eqn (11) to experimental data. Concentrations calculated in this way are compared with those found from simple mass balance in Figure 13. It is seen that the interface reaction rate control model of Eqn (11) provides a qualitatively reasonable description. In particular, it successfully predicts the complete dissolution of the steel in experiment B, the first 10 min exposure. In addition, it succeeds in predicting that less dissolution occurs in the second of the two 10 min experiments, B and C, a consequence of the higher value of present in experiment C. However, it significantly underestimates the rate of dissolution for exposures A and E (Figure 3, Figure 13). It is evident that mass transfer within the liquid slows dissolution rates, and this process must be incorporated into the description.

Diffusion Model

Assuming now that the controlling factor for iron dissolution is the diffusion of iron in the finite volume of liquid Te, and that the rate of dissolution is much more rapid, the concentration of dissolved iron is constant at the solid/liquid interface and equals . If edge effects are ignored, and convection in the isothermal melt pool can also be neglected, the situation can be modelled as shown in Figure 14. If the steel specimen is located symmetrically within the melt, then uptake of iron by the two halves of the Te(l) reservoir can be modelled as sorption of a solute by a sheet of finite thickness, X Te , with surfaces at a constant concentration, . The average concentration of iron, , within the liquid increases as a function of time according to the following equation given by Crank [START_REF] Crank | The mathematics of diffusion[END_REF]: [START_REF] Arvhult | Thermodynamic Assessment of the Fe-Te System[END_REF] where D is the diffusion coefficient (cm 2 s -1 ) of dissolved iron in the liquid tellurium and X Te is the depth (cm) of liquid Te in the crucible. The value of is set at the solubility limit at this temperature, 2.5x10 -3 mol/cm -3 .

Just as for Eqn (11), Eqn [START_REF] Arvhult | Thermodynamic Assessment of the Fe-Te System[END_REF] for diffusion controlled dissolution can be tested by applying it to each of the successive stages of the experiment charted in Figure 3, using the indicated, successively higher values of as initial values in calculating incremental increases in dissolved iron concentration. Concentration increments calculated in this way are compared with those found from simple mass balance in the bar chart of Figure 13. The only unknown value of Eqn [START_REF] Arvhult | Thermodynamic Assessment of the Fe-Te System[END_REF] is the diffusion coefficient of Fe in liquid Te. This diffusion coefficient is adjusted to reach the best fit for Eqn [START_REF] Arvhult | Thermodynamic Assessment of the Fe-Te System[END_REF] with the experimental data for dissolved iron concentration. The best fit is obtained for a diffusion coefficient equal to 4x 10 -5 cm 2 s -1 , which agrees with the usual values of diffusion coefficients in liquids. Moreover, it is seen in Figure 13 that the volume mass transfer rate control model of Eqn [START_REF] Arvhult | Thermodynamic Assessment of the Fe-Te System[END_REF] provides a relatively improved description of the iron dissolution kinetics.

On this basis, it is concluded that the rate controlling step for the dissolution process is the diffusion of dissolved Fe in the Te melt. In addition, a description of the thickening kinetics of the adherent layer, or scale, is required.

Growth of the scale

The scale is simultaneously growing at the expense of the underlying steel, and dissolving at its outer surface. The dissolution part could, in principle, be either reaction rate controlled (Eqn 11) or diffusion rate controlled (Eqn 12). Considering these two dissolution mechanisms, the growth of the adherent layer becomes Eqn [START_REF] Arvhult | Part II: Thermodynamic modelling[END_REF] for parabolic scale growth coupled with an interface reaction dissolution process, and Eqn ( 14) for parabolic scale growth coupled with a diffusion controlled dissolution process:

(13) (14)
Here is the thickness of the corrosion layer, D is the diffusion coefficient of dissolved iron in the liquid tellurium and X Te is the depth of liquid Te in the crucible, is the concentration of Fe in the corrosion layer, approximated as FeTe 2 and  the apparent dissolution constant as defined in Eqn (11).

In order to fit Eqn [START_REF] Ettenberg | [END_REF] to experimental results, a simplified equation is used, considering only the first term in the summation (only n=0). This simplification underestimates the dissolution part and leads to:

Eqn ( 13) and ( 15) are numerically solved using Mathematica software. Assuming an apparent dissolution constant  equal to s -1 and a diffusion coefficient D equal to 4x 10 -5 cm 2 s -1 (as obtained fitting respectively Eqn (11) and [START_REF] Arvhult | Thermodynamic Assessment of the Fe-Te System[END_REF] on experimental data of Figure 13), the best fits for Eqn [START_REF] Arvhult | Part II: Thermodynamic modelling[END_REF] and Eqn (15) on experimental scale thickness are given in Figure 15 together with experimental data.

The fits are seen to be reasonably good, with the exception of the 5 min exposure. Optimisation of the fit leads to k p =5 x 10 -9 cm 2 s -1 for diffusion model (Eqn (15)) and to k p =6.7 x 10 -9 cm 2 s -1 for interface reaction model (Eqn [START_REF] Arvhult | Part II: Thermodynamic modelling[END_REF]). These values of k p can be compared to experimental results obtained in Te vapour by Lobb and Robins [4] and by Magara et al [START_REF] Magara | [END_REF]. In their case, as the Te was gaseous, there was no dissolution of the corrosion product telluride layer. For Lobb and Robins, at 550°C, for 2.6 minutes and 7 hours of corrosion they obtained, respectively, an 11 and 17 µm thick scale. It can firstly be noted that this pair of values does not correspond to parabolic growth. However, using the k p values (5-6.7 x 10 -9 cm 2 s -1 obtained by fitting Eqn ( 15)-( 13)) to the present experimental data, the scale thicknesses are calculated to be equal to 12-14 µm and 159-183 µm for, respectively, 2.6 minutes and 7 hours of corrosion. The scale thicknesses obtained for 2.6 minutes of corrosion using the k p of this study and the one measured directly by Lobb and Robins [4] in Te vapour are very close, providing further support for the proposed model of Eqn (15). However, the approximation to Eqn [START_REF] Ettenberg | [END_REF] leading to Eqn (15) is valid only at small values of t, and Eqn (15) fails at long times. In Magara's work, at 650°C [START_REF] Magara | [END_REF], an Fe-11.5Cr steel in contact with gaseous Te formed a parabolic growing telluride layer assumed to be (Cr,Fe) 2 Te 3 (on the basis of unreliable EPMA analyses) and for which k p was evaluated as 10 -7 g 2 cm -4 s -1 . Adjusting for the fact that the corrosion scale is FeTe 2 , the k p value becomes 2.4 x 10 -9 cm 2 s -1 , which is very close to the values obtained in this study. The limiting step for this layer growth was deduced, from the parabolic constant activation energy and from Fe diffusion measurement in various iron tellurides [START_REF] Magara | [END_REF]- [22], to be grain boundary diffusion of iron in '-iron telluride ('-Fe 0.66 Te). In the present study, the analysed scale phase is more probably -Fe 0.95 Te 2 , for which no Fe diffusion coefficient is available. Figure 15 shows that the models (Eqn ( 13) and ( 15)) describe the experimental results relatively well, except for the thickness value after 5 min of immersion (experiment E, Figure 3), which is overestimated. Indeed, there is an important amount of dissolved Fe in the melt, as a result of the preceding experiments, and this limits the scale dissolution. For scale growth (which does not depend on the dissolved iron level), the rate is high as it corresponds to an early stage of reaction, when the scale thickness is small. It is noted that both Eqn ( 13) and ( 15) are subject to error in evaluation, as they each calculate a small difference between two large terms (scale growth and scale dissolution). An additional complication is that the solute concentration in the melt is not directly measured: it is deduced from the observed steel section loss in a calculation which neglects scale formation.

The overprediction of scale thickness from Eqns [START_REF] Arvhult | Part II: Thermodynamic modelling[END_REF] or (15) for the 5 min specimen corresponds to the underprediction of dissolution in this case (Figure 15, second point, and Figure 13, last exposure). Both models (Eqn ( 13) and ( 15)) fail for this one point, but remain reasonable overall for the dataset.

The proposed values of D, k p and  given in this study are only order of magnitude estimates and should not be considered as reliable data a.

It is concluded on the basis of Figure 13 and Figure 15, that the dissolution rate is likely controlled by diffusion of dissolved Fe in the liquid Te melt, and the growth of the adherent telluride layer is controlled by the competition between its parabolic thickening and its dissolution, the latter process being supported by metal diffusion in the melt.

Conclusions

Using the Calphad method, a thermodynamic approach is used to study the interaction between Stainless steel 316L and liquid tellurium. The developed thermodynamic database makes it possible to calculate the solubilities of the major alloying elements (Cr-Fe-Ni) in liquid Te. However, in this preliminary work, the calculated value for iron seems to be underestimated. The metal tellurides formed at 823 K by interaction with liquid Te are predicted to be Cr 3 Te 4 and FeTe 2 Corrosion tests have now been carried out, and the observed intermetallics in the corrosion layers are MTe 2 type which is partly in accord with thermodynamic predictions.

From the kinetic point of view, steel consumption is faster than scale growth, indicating that simultaneous formation and dissolution of tellurides into the melt occurs.

A mathematical model based on simultaneous scale thickening according to parabolic kinetics and scale dissolution controlled by liquid phase diffusion describes reasonably well both scaling kinetics and solute accumulation in the melt.

In order to reinforce this model, more study is required: (i) further investigation should be undertaken for thermodynamics data acquisition in the quaternary system Fe-Ni-Cr-Te, in particular to determine the effects of Ni on telluride stabilities; (ii) corrosion experiments should be carried out in liquid Te flowing at controlled rates in order to discriminate between the interface reaction and diffusion controlled mechanisms. Finally, in this study, only Fe loss has been modelled, whereas the global 316L corrosion process has to be considered.

Figure 1 :

 1 Figure 1: Binary Cr, Fe, Ni solubilities in liquid Te (at. fr.) as a function of inverse temperature (K -1 )

Figure 2 :

 2 Figure 2: Schematic diagram of the facility (a) and enlargement of the sample holder (b).

Figure 3 :

 3 Figure 3: Sequence of exposures, showing their duration and corresponding accumulated concentration of dissolved iron calculated from the loss of thickness of each sample.

Figure 4 :

 4 Figure 4: Cross-sections of 316L showing adherent solidified liquid after immersion in melt for indicated times.

Figure 5 :

 5 Figure 5: SEM image (a) and EDX maps (b), (c) of corrosion products) and solidified liquid on 316L after immersion in the melt for 1 min.

Figure 6 :

 6 Figure 6: SEM image (a) and EDX maps (b), (c) of corrosion products and solidified liquid on 316L after immersion in the melt for 30 min.

Figure 7 :

 7 Figure 7: EDX line scans across reaction zone and into solidified melt after 10 min immersion (Exposure C).

Figure 8 :

 8 Figure 8: X-Ray diffraction patterns obtained from the sample surface (covered by solidified melt) after 10 min (red) and 30 min (blue) immersion. Identification of peaks confirms the presence of Te, MoTe 2 and FeTe 2 . Presence of MoTe 2 is due to the use of Mo suspension wire for the samples.

Figure 9 :

 9 Figure 9: Te rich domain of the Cr-Fe-Te isothermal section calculated at 823 K, extrapolation from the binary systems.

Figure 10 :

 10 Figure 10: SEM-QBSD image of cross section of 316L after immersion in the melt for 30 min (Exposure D).

Figure 11 :

 11 Figure 11: Scale growth kinetics

Figure 14 :

 14 Figure 14: Diagram of the diffusion model for solute entering the melt.

Figure 15 :

 15 Figure 15: Scale growth kinetics: measured from SEM observation, ( ), and predicted from Equation 13 ( ) and Equation 15 ( ).

Table 1 :

 1 Chemical composition of 316L in wt%

	Fe	Cr	Ni	Mo Mn Si	P	N	C	S
	Bal. 16.82 10.17 2.08 1.84 0.64 0.026 0.025 0.016 0.003