domain: spi math

text:
Mathieu Fabre

Cédric Pozzolini

Yves Renard

Nitsche-based models for the unilateral contact of plates

This paper aims to present different Nitsche-based models for the unilateral contact of plate structures. Our analysis is based on the consideration of Nitsche's method on a 3D structure with kinematic assumptions of thin or thick plate theories. This approach is compared to that of Gustafsson, Stenberg and Videman which consists of Nitsche's method applied directly on a 2D plate model. To simplify the presentation, we focus on the contact of an elastic plate with a rigid obstacle. The different approaches are compared numerically in terms of reliability compared to the 3D elastic model.

Introduction and kinematics of plate models

The weak enforcement of unilateral contact conditions has already been treated by several techniques, the main ones being the penalty method, mixed/mortar methods (with standard or augmented Lagrangians) and, more recently, Nitsche's method. The penalty method is the simplest one but it leads to a supplementary approximation (it is not consistent in a strong sense) and quality of the approximation is sensitive to the choice of the penalty coefficient (see [START_REF] Chouly | An adaptation of Nitsche's method to the Tresca friction problem[END_REF][START_REF] Chouly | An overview of recent results on Nitsche's method for contact problems[END_REF][START_REF] Chouly | Symmetric and non-symmetric variants of Nitsche's method for contact problems in elasticity: Theory and numerical experiments[END_REF]). Mixed methods introduce an additional variable which, in that case, represents the contact stress. The corresponding weak form is consistent and leads numerically to an optimal convergence rate whenever an inf-sup condition is fulfilled (see [START_REF] Wohlmuth | A mortar finite element method using dual spaces for the lagrange multiplier[END_REF]), which can be a difficulty for instance with immersed boundaries or simply on the interface of several boundaries for multi-body contact. Augmented Lagrangian techniques can be used for numerical solution for instance with Uzawa techniques (see [START_REF] Renard | Generalized Newton's methods for the approximation and resolution of frictional contact problems in elasticity[END_REF]).

Nitsche's method [START_REF] Nitsche | Über ein variationsprinzip zur lösung von Dirichlet-problemen bei verwendung von teilräumen, die keinen randbedingungen unterworfen sind[END_REF] allows to treat boundary or interface conditions in a weak sense. The extension to the unilateral contact condition is proposed in [START_REF] Chouly | A Nitsche-based method for unilateral contact problems: Numerical analysis[END_REF][START_REF] Juntunen | Nitsche's method for general boundary conditions[END_REF][START_REF] Wriggers | A formulation for frictionless contact problems using a weak form introduced by Nitsche[END_REF]]. An advantage of Nitsche's method is to allow a consistent approximation without any additional unknown and without discrete inf-sup condition to satisfy. In an overview [START_REF] Chouly | An overview of recent results on Nitsche's method for contact problems[END_REF] it is presented a simple and systematic procedure to derive a family of Nitsche's formulations that have different symmetry properties and different degrees of dependency on the stabilization parameter and for various boundary and interface conditions in 3D contact problems. This family is indexed by the so-called Nitsche symmetric parameter. Let us mention that, in the context of the standard finite element method, the skew-symmetric formulation proves to be very robust with respect to the stabilization parameter (see [START_REF] Chouly | Symmetric and non-symmetric variants of Nitsche's method for contact problems in elasticity: Theory and numerical experiments[END_REF]).

We refer to [START_REF] Drouet | Optimal convergence for discrete variational inequalities modelling Signorini contact in 2D and 3D without additional assumptions on the unknown contact set[END_REF] for a more detailed review on a priori error estimates for contact problems in elasticity. The case of contact between two elastic bodies, still in small deformations framework, is addressed in [START_REF] Chouly | Symmetric and non-symmetric variants of Nitsche's method for contact problems in elasticity: Theory and numerical experiments[END_REF][START_REF] Chouly | An unbiased Nitsche's approximation of the frictional contact between two elastic structures[END_REF]. In [START_REF] Chouly | An unbiased Nitsche's formulation of large deformation frictional contact and self-contact[END_REF] a variant of an unbiased Nitsche's method is proposed for the contact between two elastic bodies without making any difference between master and slave contact surfaces, which is an advantage for treatment of self-contact or multi-body contact.

⇤ ESI Group, 70 Rue Robert, 69006 Lyon, France. email: Mathieu.Fabre@esi-group.com † ESI Group, 70 Rue Robert, 69006 Lyon, France. email: Cedric.Pozzolini@esi-group.com ‡ Université de Lyon, CNRS, INSA-Lyon, ICJ UMR5208, LaMCoS UMR5259, F-69621, Villeurbanne, France. email: Yves.Renard@insa-lyon.fr With regard to unilateral contact of thin structures, a Nitsche approximation of the unilateral contact between a Kirchhoff-Love plate and a rigid obstacle is presented in [START_REF] Gustafsson | A stabilised finite element method for the plate obstacle problem[END_REF][START_REF] Gustafsson | Finite Element Methods for Contact Problems[END_REF]. In these references, the origin of the proposed method is the Kirchhoff-Love plate model subject to a unilateral contact condition. It is shown that an approximation using Lagrange multipliers with an additional stabilization term leads to a Nitsche method by static condensation of the multiplier.

The approach we consider in this paper is different. We start from the three-dimensional problem on which we apply the unilateral contact condition using Nitsche's method, as in [START_REF] Chouly | Symmetric and non-symmetric variants of Nitsche's method for contact problems in elasticity: Theory and numerical experiments[END_REF]. Then, we apply the kinematics of different models of thin and thick plates, which allows us to write different approximated models of the unilateral contact of plate structures.

Let us now introduce the notations used and the plate models that we will consider. A thick or a thin elastic plate is a flat structure one dimension of which, called thickness, is small or very small compared to the others. For this kind of structure, starting from a priori hypotheses on the expression of the displacement fields, a two-dimensional problem is usually derived from the three-dimensional elasticity formulation by means of integration in the thickness direction (see for instance [START_REF] Ph | Mathematical Analysis of Thin Plate Models[END_REF]). The unknown variables are then described on the mid-plane of the plate.

In what follows, we use the canonical orthonormal basis (e 1 , e 2 , e 3 ) for the coordinate space R 3 . Let ⌦ be an open, bounded, connected subset of the plane Span(e 1 , e 2 ), with Lipschitz-continuous boundary which defines the mid-plane of the plate. Then, the plate in its stress-free reference configuration coincides with the domain

⌦ " := ⌦ ⇥ ( ", ") = (x 1 , x 2 , x 3 ) 2 R 3 | (x 1 , x 2 ) 2 ⌦ and x 3 2 ( ", ") ,
where 2" > 0 is called the thickness. Note that in reference configuration the mid-plane is in the plane x 3 = 0. We define on the boundary of ⌦ " an outward unit normal vector n and a unit tangential vector t such that the pair (n, t) defines a direct orthonormal coordinate system. We consider a plate made of an elastic, homogeneous and isotropic material, whose mechanical constants are its Young's modulus E and its Poisson's ratio ⌫, with, as usual, E > 0 and 0  ⌫ < 1 2 . In addition, ij is Kronecker's symbol and the summation convention over repeated indices is adopted, Greek indices varying in {1, 2} and Latin indices in {1, 2, 3}. In the following, @ i stands for the partial derivative with respect to x i and the second derivatives are @ 2 ij := @ 2 @x i @x j .

In general plate theory, it is assumed that a mid-surface plane can be used to represent the three-dimensional plate in two-dimensional form. We consider the following very common linear approximation of the three-dimensional displacements for x = (x 1 , x 2 , x 3 ) 2 ⌦ " : 8 < :

u 1 (x 1 , x 2 , x 3 ) = u 1 (x 1 , x 2 ) + x 3 ✓ 1 (x 1 , x 2 ), u 2 (x 1 , x 2 , x 3 ) = u 2 (x 1 , x 2 ) + x 3 ✓ 2 (x 1 , x 2 ), u 3 (x 1 , x 2 , x 3 ) = u 3 (x 1 , x 2 ) + x 3 ✓ 3 (x 1 , x 2 ). (1)
In these expressions, u 1 and u 2 are the membrane displacements of the mid-plane points, u 3 is the deflection, while ✓ 1 , ✓ 2 are the section rotations, and ✓ 3 is the pinching. In the case of a homogeneous isotropic material, the variational plate model often splits into two independent problems: the first, called the membrane problem, deals only with membrane displacements, while the second, called the bending problem, concerns deflection and rotations. Note also that some enhanced models have been proposed with a quadratic approximation in the variable x 3 for the transverse displacement u 3 instead of a linear one (see [START_REF] Büchter | Three-dimensional extension of non-linear shell formulation based on the enhanced assumed strain concept[END_REF][START_REF] Hauptmann | A systematic development of 'solid-shell' element formulations for linear and non-linear analyses employing only displacement degrees of freedom[END_REF] for instance) especially when a significant pinching of the plate occurs. However, we limit ourselves to a linear approximation along x 3 .

In this paper, we consider two types of plates models. On the one hand, the classical Kirchhoff-Love and Mindlin-Reissner models which are based on an additional plane stress approximation assumption, and on the other hand some simple solid plate models which do not include this additional assumption.

The Kirchhoff-Love plate model is based on the plane stress approximation and the following kinematic assumptions which can be seen as a special case of (1):

8 <

:

✓ 1 = @ 1 u 3 , ✓ 2 = @ 2 u 3 , ✓ 3 = 0. (2)
Therefore, deflection is the only unknown for the bending Kirchhoff-Love plate problem and this displacement is independent of membrane displacements (i.e in-plane displacements of the midsurface) in the homogeneous isotropic case.

In the Mindlin-Reissner theory, a plane stress approximation is also considered and the normal to the mid-surface remains straight but not necessarily perpendicular to the mid-surface. Shear stress is no longer neglected. The kinematics is also a special case of (1) obtained by neglecting the pinch, assuming

✓ 3 = 0. (3)
Finally, we consider two simple plate models which are not based on the plane stress approximation: a Mindlin-Reissner like plate model and the more general model (1) which will be called in the following the solid plate model, since it is the restriction to a plane geometry of a solid shell model (see [START_REF] Hauptmann | A systematic development of 'solid-shell' element formulations for linear and non-linear analyses employing only displacement degrees of freedom[END_REF] for instance) as long as the finite element approximation is not concerned.

The rest of the paper is outlined as follows. In Section 2 we introduce the classical formulation of the 3D frictional contact problem, with some specific notations from plate theory. We focus on a model problem which consists of a unilateral contact between an elastic body and a rigid support (Signorini's problem). Then, we introduce Nitsche's method for the contact condition. Section 3 is devoted to obtaining the Nitsche-based formulation for the classical Kirchhoff-Love and Mindlin-Reissner plate models. We obtain in Section 4 the Nitsche-based formulation for two additional solid plate models: the first one is based on Mindlin-Reissner kinematics but without the plane-stress assumption and the second one takes into account the enriched kinematics [START_REF] Alart | A generalized Newton method for contact problems with friction[END_REF]. In Section 5, the different approach of Nitsche's method proposed in [START_REF] Gustafsson | A stabilised finite element method for the plate obstacle problem[END_REF] for the Kirchhoff-Love plate model is recalled for comparison. Then, Section 6 deals with the finite element approximation of the proposed formulations and presents some numerical tests. The convergence in L 2 -norm and H 1 -semi-norm of the displacement toward a fully 3D-solution computed on a very fine mesh is analyzed. A brief study of the choice of Nitsche's parameter is given in Section 7 and the paper is ended by some concluding remarks in Section 8.

3D-Contact problems: weak formulations

Loading and boundary conditions

The boundary @⌦ " of the domain ⌦ " is assumed to be partitioned into three parts. Due to the particular shape of the domain, the partition is composed of the "lateral" one " 0 , the "upper" one " + and the "lower" one " , i.e.

@⌦ " = " 0 [ " + [ " , " 0 = @⌦ ⇥ ( ", "), " + = ⌦ ⇥ {"}, " = ⌦ ⇥ { "}.
We suppose that the lateral part " 0 consists in two non-overlapping sub-parts " D and " N . On " D (resp. " N ) displacements u (resp. tractions) are prescribed. We consider the following situation:

• For the sake of simplicity, the body is clamped on " D := D ⇥ ( ", ") which is assumed to be a non-zero Lebesgue measure part of the boundary " 0 ⇢ @⌦ " , that is to say u " = 0 on " D , and is traction free on the complementary part " N .

• In addition, the body can be subjected to a volume force f V 2 L 2 (⌦ " ; R 3 ) (such as gravity);

• As far as loading is concerned, the upper part, " + , is loaded by a surface force `2 L 2 ( " + ; R 3 )

• The actual surface on which the body comes into contact with the obstacle is not known in advance, but is assumed contained in " .

The space of admissible displacements for the three-dimensional problem is

V " := { v 2 H 1 (⌦ " ; R 3 ) | v = 0 on " D }.

Unilateral boundary conditions and weak form

On the potential contact boundary " , the vector n still denote the outward unit vector (n = e 3 , here). Moreover, we also consider the unit vector n, the contact normal. In common models, n can be chosen different from n (see Figure 1), for instance, to be the normal vector to the obstacle or, in a two-body contact, the normal vector to a second deformable body or even to an intermediate surface.

We decompose the displacement field u : ⌦ " ! R 3 and density of surface forces (u)n :

" ! R in normal and tangential parts as follows:

u = u n n + u t and (u)n = n (u)n + t (u).
O bs ta cl e Let us now introduce the static equation with Signorini's conditions along the plate. We assume that the motion of the plate is limited by a rigid obstacle, located below the plate, determined by

g : ⌦ " ! R [ {+1},
such that the displacement is constrained to belong to the convex set K " ⇢ V " given by

K " := {v 2 V " | v n  g on ⌦ " }.
The three-dimensional contact problem in linear elasticity consists in finding the displacement field u and the traction (u)n verifying the strong equations (4) and the contact conditions described hereafter:

8 > > < > > : div (u) = f V in ⌦ " , (u)n = `on " + , u = 0 on " D , (u)n = 0 on " N , (4)
where stands for the stress tensor field and div denotes the divergence operator of tensor valued functions. Denoting the linearized strain tensor field ij (u) := 1 2 (@ i u j + @ j u i ) ,

we consider Hooke's law for isotropic material:

ij = E⌫ (1 + ⌫)(1 2⌫) kk ij + E 2(1 + ⌫) ij .
The unilateral contact conditions (or Kuhn-Tucker conditions) are given by:

u n  g, n (u)  0 and (u n g) n (u) = 0 on " . (5)
In the frictionless contact case, this condition is simply completed by

t (u) = 0, on " , (6)
whereas the so-called Coulomb's static friction condition on " reads

( k t (u)k  F n (u) if u t = 0, t (u) = F n (u) u t ku t k otherwise, (7)
where F 0 is the Coulomb coefficient of friction and k.k stands for the euclidean norm in R 3 .

Remark 2.1. Recall that the Coulomb's static friction condition is rather artificial but that is leads to an interesting intermediate problem. The true Coulomb's friction law involves tangential contact velocity and not tangential displacement. However, a problem similar to the one discussed here is obtained by time discretization of the quasi-static frictional contact evolution problem. In this case u, f and `stand for u((i + 1) t), f ((i + 1) t)) and `((i + 1) t) respectively and u t has to be replaced by [u t ((i + 1) t) u t (i t)]/ t, where t denotes the time step. For simplicity and without any loss of generality only the static case described above will be considered in the following.

From Green's formula and equations (4), we get the weak variational formulation for all v 2 V " :

Z ⌦ " (u) : (v) d⌦ Z " (u)n • v d = Z ⌦ " f V • v d⌦ + Z " + `• v d . (8)
In the following, the membrane stress resultants (N ↵ ), the bending moments stress resultants (M ↵ ), Q , P the pinch stress and the shear stress (T ) are defined as (for more details see [START_REF] Ph | Mathematical Analysis of Thin Plate Models[END_REF] and [START_REF] Ciarlet | A justification of two-dimensional linear plate model[END_REF])

N ↵ := Z " " ↵ (u) dx 3 , M ↵ := Z " " x 3 ↵ (u) dx 3 , T := Z " " 3 (u) dx 3 , Q := Z " " x 3 3 (u) dx 3 , P := Z " " 33 (u) dx 3 .
Remark 2.2. The Korn inequality that is satisfied for this problem is simply adapted from the classical one, which allows to prove the existence and the uniqueness of the solution to problem (8) subject to (5), ( 6) for any fixed " > 0, thanks to Lax-Milgram theorem, see [START_REF] Fichera | Problemi elastostatici con vincoli unilaterali: il problema di Signorini con ambigue condizioni al contorno[END_REF]. For the general problem (8) subject to (5),( 7) the existence of a solution is only available for small Coulomb coefficient of friction F 0 (see [START_REF] Necas | On the solution of the variational inequality to the Signorini problem with small friction[END_REF][START_REF] Jarusek | Contact problems with bounded friction, Coercive case[END_REF][START_REF] Eck | Existence results for the static contact problem with Coulomb friction[END_REF]).

Weak formulation using Nitsche's method

A reformulation of the previous contact conditions comes from the augmented Lagrangian formulation of contact problems. Let r be a given real positive number. As in [START_REF] Alart | A generalized Newton method for contact problems with friction[END_REF][START_REF] Chouly | A Nitsche-based method for unilateral contact problems: Numerical analysis[END_REF][START_REF] Chouly | Symmetric and non-symmetric variants of Nitsche's method for contact problems in elasticity: Theory and numerical experiments[END_REF], the contact conditions ( 5) are rewritten as

n (u) = [ n (u) + r(g u n )] , (9)
where we denote [x] := (|x| x)/2 the negative part. We can reformulate as well Coulomb's friction condition using the projection P B(0,s) onto the ball B(0, s) defined by

P B(0,s) (q) := 8 < : ⇧ t (q) if k⇧ t (q)k  s, s ⇧ t (q) k⇧ t (q)k otherwise,
where ⇧ t is the tangential projection defined as ⇧ t := I n ⌦ n, with I the identity operator. Still for a given real positive number r, Coulomb's friction condition is equivalent to the non-smooth equation

t (u) = P B(0,F [ n(u)+r(g un)] ) ( t (u) ru t ). ( 10
)
Let now O 2 R be a fixed parameter that we use to recover different variants of the Nitsche's method, as in the 3D linear elastic setting (see, e.g., [START_REF] Chouly | Symmetric and non-symmetric variants of Nitsche's method for contact problems in elasticity: Theory and numerical experiments[END_REF]) for O = 1, 0, 1 corresponding to symmetric, non-symmetric and skew-symmetric contributions of the contact terms, respectively. We insert finally the expressions ( 9) and ( 10) into ( 8) and obtain, formally, our Nitsche-based formulation for frictional contact. In the decomposition

Z " (u)n • v d = Z " n (u)v n (v) d + Z " t (u) • v t d , (11)
the normal contact term (obtained from ( 9) and the classical decomposition introduced in [START_REF] Chouly | Symmetric and non-symmetric variants of Nitsche's method for contact problems in elasticity: Theory and numerical experiments[END_REF][START_REF] Fabre | A fictitious domain method for frictionless contact problems in elasticity using Nitsche's method[END_REF][START_REF] Chouly | An unbiased Nitsche's formulation of large deformation frictional contact and self-contact[END_REF]) can be written:

Z " n (u)v n d = 1 r Z " n (u) (rv n O n (v)) d + O r Z " n (u) n (v) d = 1 r Z " [ n (u) + r(g u n )] (rv n O n (v)) d + O r Z " n (u) n (v) d = 1 r Z ⌦ [ n (u) + r(g u n )] (rv n O n (v)) x 3 = " d⌦ + O r Z ⌦ n (u) n (v) x 3 = " d⌦, (12)
and with frictional contact term:

Z " t (u) • v t d = 1 r Z ⌦ P B(0, F [ n(u)+r(g un)] ) ( t (u) ru t ) • (rv t O t (v)) x 3 = " d⌦ + O r Z ⌦ t (u) • t (v) x 3 = " d⌦. (13)
Now, inserting the expressions ( 12) and ( 13) of ( 11) into the weak problem [START_REF] Chouly | An adaptation of Nitsche's method to the Tresca friction problem[END_REF], we obtain the so-called Nitsche-based method:

8 > > > > > > > < > > > > > > > : Z ⌦ " (u) : (v) d⌦ + 1 r Z ⌦ [ n (u) + r(g u n )] (rv n O n (v)) x 3 = " d⌦ 1 r Z ⌦ P B(0,F [ n(u)+r(g un)] ) ( t (u) ru t ) • (rv t O t (v)) x 3 = " d⌦ O r Z ⌦ ( (u)n) • ( (v)n) x 3 = " d⌦ = Z ⌦ " f V • v d⌦ + Z " + `• v d . (14)
Note that because of the terms in (u)n and (v)n, this formulation only makes sense for both (u)n and (v)n in L 2 ( " ; R 3 ) (except for F = O = 0). In our case, this regularity is automatically obtained for u, v 2 V " satisfiying additionnaly the kinematic assumptions (1) which are considered below. Note that with these assumptions, the decomposition into normal and tangential components is as follows:

v t := v v n n = v (v ↵ n↵ + x 3 ↵ n↵ + v 3 n3 + x 3 3 n3 ) n, u t := u u n n = u (u ↵ n↵ + x 3 ✓ ↵ n↵ + u 3 n3 + x 3 ✓ 3 n3 ) n.
3 Weak formulation for classical plate models under plane stress assumption

Weak formulation under Kirchhoff-Love assumptions

We now aim to make explicit each term of the formulation ( 14) under Kirchhoff-Love assumptions. The plane stress assumption consists in assuming that the components i3 (u) are negligible compared to ↵ (u) which, for an isotropic material under Hooke's law, leads to the relations

↵ (u) = 1 + ⌫ E ↵ (u) ⌫ E ⌘⌘ (u) ↵ , ↵ (u) = E 1 ⌫ 2 [(1 ⌫) ↵ (u) + ⌫ ⌘⌘ (u) ↵ ].
The linearized strain tensor field reads in that case

↵ (u) = ↵ (u) x 3 @ 2 ↵ u 3 ,
and the stress tensor

↵ (u) = E 1 ⌫ 2 (1 ⌫)[ ↵ (u) x 3 @ 2 ↵ u 3 ] + ⌫[ ⌘⌘ (u) x 3 @ 2 ⌘⌘ u 3 ] ↵ .
Then, we obtain

Z ⌦ " (u) : (v) d⌦ = Z ⌦ 2 6 6 6 4 ✓Z " " ↵ (u) dx 3 ◆ | {z } N ↵ ↵ (v) ✓Z " " x 3 ↵ (u) dx 3 ◆ | {z } M ↵ @ 2 ↵ v 3 3 7 7 7 5 d⌦.
Thanks to the model and according to Hooke's law, the membrane stress resultants (N ↵ ) and the bending moments stress resultants (M ↵ ) read

N ↵ = Z " " ↵ (u) dx 3 = 2"E 1 ⌫ 2 [(1 ⌫) ↵ (u) + ⌫ ⌘⌘ (u) ↵ ], M ↵ = Z " " x 3 ↵ (u) dx 3 = 2" 3 E 3(1 ⌫ 2 ) [(1 ⌫)@ 2 ↵ u 3 + ⌫@ 2 ⌘⌘ u 3 ↵ ].
We rewrite the right-hand side of ( 14) as follows:

Z ⌦ " f V .v d⌦ + Z " + `• v d = Z ⌦ F i v i G ↵ @ ↵ v 3 d⌦, with
F i := Z " " f V i dx 3 + `i, G ↵ := Z " " x 3 f V ↵ dx 3 + "`↵.
Neglecting the components i3 (u) according to the plane stress assumption, the terms involved in [START_REF] Chouly | An unbiased Nitsche's approximation of the frictional contact between two elastic structures[END_REF] can be expressed as follows:

[ n (u) + r(g u n )] = [ ↵ (u)n n↵ + r(g u ↵ n↵ "@ ↵ u 3 n↵ u 3 n3 )] (15)
rv n O n (v) = r (v ↵ n↵ + "@ ↵ v 3 n↵ + v 3 n3 ) O( ↵ (v)n n↵ ). (16)
Now, the term in the frictional part ( 13) can be made explicit as

t (u) ru t = ↵ (u)n e ↵ ( ↵ (u)n n↵ )n r(u u n n). (17)
In our situation, we have n = e 3 , then the three last expressions are then simplified into:

[ n (u) + r(g u n )] = r [g u ↵ n↵ "@ ↵ u 3 n↵ u 3 n3 ] , rv n O n (v) = r (v ↵ n↵ + "@ ↵ v 3 n↵ + v 3 n3 ) , t (u) ru t = ru t . Special case n = n = e 3
In the special, but common case n = n = e 3 , the terms (15), ( 16) and ( 17) reduce to:

[ n (u) + r(g u n )] = r [g + u 3 ] , rv n O n (v) = rv 3 , t (u) ru t = ru t = r(u ↵ e ↵ "@ ↵ u 3 e ↵ ).
Finally, denoting

V KL := { v 2 H 1 (⌦; R 3 ) | v 3 2 H 2 (⌦), v = 0 on D , @ n v 3 = 0 on D },
we end up with the following formulation considering the Kirchhoff-Love assumptions, still in the case n = n = e 3 :

8 > > > > > > > > > > > < > > > > > > > > > > > : Find u 2 V KL such that for all v 2 V KL Z ⌦ N ↵ ↵ (v) M ↵ @ 2 ↵ v 3 d⌦ r Z ⌦ [g + u 3 ] v 3 d⌦ +r Z ⌦ P B(0,rF [g+u 3 ] ) (u ↵ e ↵ "@ ↵ u 3 e ↵ ) • (v ↵ e ↵ "@ ↵ v 3 e ↵ ) d⌦ = Z ⌦ F i v i G ↵ @ ↵ v 3 d⌦.
We can now briefly review the remaining terms when Kirchhoff-Love assumptions are taken into account. The terms involved in unilateral contact and friction conditions are reduced to n (u) + r(g u n ) = r(g + u 3 ) and t (u) ru t = ru t (in the case n = n = e 3 ), which means that both contact and friction are approximated by a penalty method with r as the penalty coefficient. Indeed, with Kirchhoff-Love assumptions, no additional contribution than penalty terms can be obtained with the proposed Nitsche's method mainly because of the plane stress approximation.

Weak formulation under Mindlin-Reissner assumptions

Let us now perform the same analysis but considering Mindlin-Reissner assumptions. In the Mindlin-Reissner plate model, the pinch stress 33 is still neglected and Hooke's law still leads to the relations

↵ (u) = 1 + ⌫ E ↵ (u) ⌫ E ⌘⌘ (u) ↵ , ↵ (u) = E 1 ⌫ 2 [(1 ⌫) ↵ (u) + ⌫ ⌘⌘ (u) ↵ ].
However, the shear stress is no longer neglected and is approximated by

↵3 (u) =  E 1 + ⌫ ↵3 (u),
with  the shear correction factor [START_REF] Batoz | Modélisation des Structures par Éléments Finis[END_REF], which standard value is  = 5 6 . The linearized strain tensor field reads in that case

↵ (u) = 1 2 (@ ↵ u + @ u ↵ ) = ↵ (u) + x 3 ↵ (✓) and ↵3 (u) = 1 2 (@ ↵ u 3 + ✓ ↵ ),
and with

↵ (u) = E 1 ⌫ 2 ((1 ⌫)[ ↵ (u) x 3 ↵ (✓)] + ⌫[ ⌘⌘ (u) x 3 ⌘⌘ (✓)] ↵ ) , ↵3 (u) =  E 2(1 + ⌫) (@ ↵ u 3 + ✓ ↵ ).
The work of the internal efforts end up taking the form Thanks to the model and according to Hooke's law, the membrane stress resultants (N ↵ ), the bending moments stress resultants (M ↵ ) and the shear stress resultants (T ) are given by

Z ⌦ " (u) : (v) d⌦ = Z ⌦ 2 6 6 6 4 ✓Z " " ↵ (u) dx 3 ◆ | {z } N ↵ ↵ (v) + ✓Z " " x 3 ↵ (u) dx 3 ◆ | {z } M ↵ ↵ ( ) 3
N ↵ = 2"E 1 ⌫ 2 [(1 ⌫) ↵ (u) + ⌫ ⌘⌘ (u) ↵ ], M ↵ = 2" 3 E 3(1 ⌫ 2 ) [(1 ⌫) ↵ (✓) + ⌫ ⌘⌘ (✓) ↵ ], T = "E 1 + ⌫ [@ u 3 + ✓ ].
The right-hand side of ( 14) is now rewritten as

Z ⌦ " f V .v d⌦ + Z " + `• v d = Z ⌦ F i v i G ↵ ↵ d⌦,
with still

F i = Z " " f V i dx 3 + `i, G ↵ = Z " " x 3 f V ↵ dx 3 + "`↵.
We pass to the non-penetration condition terms from [START_REF] Chouly | An unbiased Nitsche's approximation of the frictional contact between two elastic structures[END_REF]. We obtain

[ n (u) + r(g u n )] = [ ij (u)n j ni + r(g u ↵ n↵ "✓ ↵ n↵ u 3 n3 )] (18)
rv n O n (v) = r (v ↵ n↵ + " ↵ n↵ + v 3 n3 ) O ij (v)n j ni . (19)
still assuming 33 = 0 and the corresponding expression for ↵ (v), 33 (v) and ↵3 (v). The term in the friction condition ( 13) reads now

t (u) ru t = ij (u)n j e i ( ij (u)n j ni )n r(u u n n), (20)
and the simplifications coming from n = e 3 leads to

[ n (u) + r(g u n )] = [ ↵3 (u)n ↵ + r(g u ↵ n↵ "✓ ↵ n↵ u 3 n3 )] , rv n O n (v) = r (v ↵ n↵ + " ↵ n↵ + v 3 n3 ) + O ↵3 (v)n ↵ , t (u) ru t = ↵3 (u)(n ↵ n e ↵ ) ru t . Special case n = n = e 3
In the special case n = n = e 3 , we obtain:

[ n (u) + r(g u n )] = r [(g + u 3 )] , rv n O n (v) = rv 3 , t (u) ru t = ↵3 (u)e ↵ ru t = E 2(1 + ⌫) (@ ↵ u 3 + ✓ ↵ )e ↵ r(u ↵ e ↵ + "✓ ↵ e ↵ ),
rv t O t (v) = rv t + O ↵3 (v)e ↵ = r(v ↵ e ↵ + " ↵ e ↵ ) + O E 2(1 + ⌫) (@ ↵ v 3 + ↵ )e ↵ .
Finally, denoting

V MR := { (v, ) 2 H 1 (⌦; R 3 ) ⇥ H 1 (⌦; R 2 ) | v = 0 on D , = 0 on D },
we end up with the following formulation under the Mindlin-Reissner assumptions, in the case

n = n = e 3 : 8 > > > > > > > > > > > < > > > > > > > > > > > : Find (u, ✓) 2 V MR such that for all (v, ) 2 V MR Z ⌦ N ↵ ↵ (v) + M ↵ ↵ ( ) + T (@ v 3 + ) d⌦ r Z ⌦ [g + u 3 ] v 3 d⌦ + 1 r Z ⌦ P B(0,rF [g+u 3 ] ) ( ↵3 (u)e ↵ + ru t ) • (rv t + O ↵3 (v)e ↵ ) d⌦ O r Z ⌦ ↵3 (u) ↵3 (v) d⌦ = Z ⌦ F i v i G ↵ ↵ d⌦.
Note in particular the presence of the shear stress in the friction term t (u) ru t which disappear in the Kirchhoff-Love plate model. This means that the friction condition is no longer reduced to a penalty approximation. The Mindlin-Reissner model is richer in this sense. However, the contact term remains a pure penalty term.

4 Weak formulation for some simple solid plate models

Weak formulation under Mindlin-Reissner kinematics without any plane stress approximation

To retrieve a non-trivial Nitsche contact term, it is necessary to abandon the plane stress approximation 33 = 0. This leads to plate models that are sometimes called solid plate models. The simplest one we will consider is built on the Mindlin-Reissner kinematic assumptions. According to the three-dimensional Hooke's law, we obtain for a compressible material

↵ (u) = E 1 + ⌫  ↵ (u) + x 3 ↵ (✓) + ⌫ 1 2⌫ [ ⌘⌘ (u) + x 3 ⌘⌘ (✓)] ↵ , ↵3 (u) = E 2(1 + ⌫) (@ ↵ u 3 + ✓ ↵ ), 33 (u) = E⌫ (1 + ⌫)(1 2⌫) [ ⌘⌘ (u) + x 3 ⌘⌘ (✓)] ,
and the linearized strain tensor field is still given by

↵ (u) = 1 2 (@ ↵ u + @ u ↵ ) = ↵ (u) + x 3 ↵ (✓) and ↵3 = 1 2 (@ ↵ u 3 + ✓ ↵ ).
The membrane stress resultants (N ↵ ), the bending moments stress resultants (M ↵ ) and the shear stress resultants (T ) have close but not identical expression compared to the Mindlin-Reissner model:

N ↵ = 2"E (1 + ⌫)(1 2⌫) [(1 2⌫) ↵ (u) + ⌫ ⌘⌘ (u) ↵ ], M ↵ = 2" 3 E 3(1 + ⌫)(1 2⌫) [(1 2⌫) ↵ (✓) + ⌫ ⌘⌘ (✓) ↵ ], T = "E 1 + ⌫ [@ u 3 + ✓ ].
The non-penetration condition terms from [START_REF] Chouly | An unbiased Nitsche's approximation of the frictional contact between two elastic structures[END_REF] becomes

[ n (u) + r(g u n )] = [ ij (u)n j ni + r(g u ↵ n↵ "✓ ↵ n↵ u 3 n3 )] (21)
rv n O n (v) = r (v ↵ n↵ + " ↵ n↵ + v 3 n3 ) O ij (v)n j ni . ( 22
)
The term in the friction condition ( 13) reads now

t (u) ru t = ij (u)n j e i ( ij (u)n j ni )n r(u u n n), (23)
and the simplifications coming from n = e 3 lead to

[ n (u) + r(g u n )] = [ i3 (u)n i + r(g u ↵ n↵ "✓ ↵ n↵ u 3 n3 )] , rv n O n (v) = r (v ↵ n↵ + " ↵ n↵ + v 3 n3 ) + O i3 (v)n i , t (u) ru t = i3 (u)(n i n e i ) ru t . Special case n = n = e 3
In the special case n = n = e 3 , we obtain:

[ n (u) + r(g u n )] = [ 33 (u) + r(g + u 3 )] =  E⌫ (1 + ⌫)(1 2⌫) [ ⌘⌘ (u) " ⌘⌘ (✓)] + r(g + u 3 ) , rv n O n (v) = rv 3 O 33 (v), = rv 3 OE⌫ (1 + ⌫)(1 2⌫) [ ⌘⌘ (v) " ⌘⌘ ( )] , t (u) ru t = ↵3 (u)e ↵ ru t = E 2(1 + ⌫) (@ ↵ u 3 + ✓ ↵ )e ↵ r(u ↵ e ↵ + "✓ ↵ e ↵ ),
rv t O t (v) = rv t + O ↵3 (v)e ↵ = r(v ↵ e ↵ + " ↵ e ↵ ) + O E 2(1 + ⌫) (@ ↵ v 3 + ↵ )e ↵ .
We then end up with the following formulation in this special case:

8 > > > > > > > > > > > > > > > < > > > > > > > > > > > > > > > : Find (u, ✓) 2 V MR such that for all (v, ) 2 V MR Z ⌦ N ↵ ↵ (v) + M ↵ ↵ ( ) + T (@ v 3 + ) d⌦ 1 r Z ⌦ [ 33 (u) + r(g + u 3 )] (O 33 (v) + rv 3 ) d⌦ + 1 r Z ⌦ P B(0,F [ 33 (u)+r(g+u 3 )] ) ( ↵3 (u)e ↵ + ru t ) • (rv t + O ↵3 (v)e ↵ ) d⌦ O r Z ⌦ 33 (u) 33 (v) + ↵3 (u) ↵3 (v) d⌦ = Z ⌦ F i v i G ↵ ↵ d⌦.
We note the presence of the pinch stress 33 (u) into the contact term, which means that the contact condition do not reduce to a pure penalty term. However, the pinch stress is proportional to ⌘⌘ (u), the divergence of u and this only represent a coupling with the Poisson effect. This is due to the fact that no pinch deformation is considered with the Mindlin-Reissner kinematics.

Weak formulation under solid plate displacements assumptions

Now we focus our attention on the richest kinematics (1) that we insert into [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF] and obtain, as previously, our Nitsche-based formulation for frictional contact. The work of elastic forces reads

Z ⌦ " (u) : (v) d⌦ = Z ⌦ 2 6 6 6 4 ✓Z " " ↵ (u) dx 3 ◆ | {z } N ↵ ↵ (v) + ✓Z " " x 3 ↵ (u) dx 3 ◆ | {z } M ↵ ↵ ( ) d⌦ + Z ⌦ 2 6 6 6 4 ✓Z " " 3 (u) dx 3 ◆ | {z } T (@ v 3 + )
Taking into account (1) the linear approximation of the three-dimensional displacements in ⌦ " , isotropic Hooke's law is given by

↵ (u) = E 1 + ⌫ [ ↵ (u) + ⌫ 1 2⌫ kk (u) ↵ ], ↵3 (u) = E 1 + ⌫ ↵3 (u), 33 (u) = E 1 + ⌫ [ 33 (u) + ⌫ 1 2⌫ kk (u)],
with the linearized strain tensor field

↵ (u) = ↵ (u) + x 3 ↵ (✓), ↵3 (u) = 1 2 (@ ↵ u 3 + x 3 @ ↵ ✓ 3 + ✓ ↵ ) , 33 (u) = ✓ 3 ,
and we get

↵ (u) = E 1 + ⌫  ↵ (u) + x 3 ↵ (✓) + ⌫ 1 2⌫ [ ⌘⌘ (u) + x 3 ⌘⌘ (✓) + ✓ 3 ] ↵ , ↵3 (u) = E 2(1 + ⌫) (@ ↵ u 3 + x 3 @ ↵ ✓ 3 + ✓ ↵ ), 33 (u) = E 1 + ⌫  ✓ 3 + ⌫ 1 2⌫ [ ⌘⌘ (u) + x 3 ⌘⌘ (✓) + ✓ 3 ] .
The expressions of N ↵ , M ↵ ,T ↵ and P are given by

N ↵ = 2"E 1 + ⌫  ↵ (u) + ⌫ 1 2⌫ ⇣ ⌘⌘ (u) + ✓ 3 ⌘ ↵ , M ↵ = 2" 3 E 3(1 + ⌫)  ↵ (✓) + ⌫ 1 2⌫ ⌘⌘ (✓) ↵ , T = "E 1 + ⌫ [@ u 3 + ✓ ] , Q = " 3 E 3(1 + ⌫) @ ✓ 3 , P = 2"E 1 + ⌫  ✓ 3 + ⌫ 1 2⌫ ⇣ ⌘⌘ (u) + ✓ 3 ⌘ .
The non-penetration condition terms from ( 12) reads

[ n (u) + r(g u n )] = [ ij (u)n j ni + r(g u i ni + "✓ i ni )] (25)
rv n O n (v) = r (v i ni " i ni ) O ij (v)n j ni . (26)
The term in the friction condition ( 13) is still

t (u) ru t = ij (u)n j e i ( ij (u)n j ni )n r(u u n n), (27)
and the simplifications coming form n = e 3 leads to

[ n (u) + r(g u n )] = [ i3 (u)n i + r(g u i ni + "✓ i ni )] , rv n O n (v) = r (v i ni " i ni ) + O i3 (v)n i , t (u) ru t = i3 (u)(n i n e i ) ru t .
Special case n = n = e 3

In the special case n = n = e 3 , we obtain:

[ n (u) + r(g u n )] = [ 33 (u) + r(g + u 3 "✓ 3 )] , =  E 1 + ⌫  ✓ 3 + ⌫ 1 2⌫ [ ⌘⌘ (u) " ⌘⌘ (✓) + ✓ 3 ] + r(g + u 3 "✓ 3 ) , rv n O n (v) = r(v 3 " 3 ) O 33 (v), = r(v 3 " 3 ) OE 1 + ⌫  3 + ⌫ 1 2⌫ [ ⌘⌘ (v) " ⌘⌘ ( ) + 3 ] , t (u) ru t = ↵3 (u)e ↵ ru t = E 2(1 + ⌫) (@ ↵ u 3 "@ ↵ ✓ 3 + ✓ ↵ )e ↵ r(u ↵ e ↵ + "✓ ↵ e ↵ ),
rv t O t (v) = rv t + O ↵3 (v)e ↵ = r(v ↵ e ↵ + " ↵ e ↵ ) + OE 2(1 + ⌫) (@ ↵ v 3 "@ ↵ 3 + ↵ )e ↵ .
Finally, denoting

V SP := { (v, ) 2 H 1 (⌦; R 3 ) ⇥ H 1 (⌦; R 3 ) | v = 0 on D , = 0 on D },
and defining additionnaly G 3 := Z " "

x 3 f V 3 dx 3 + "`3, the formulation of the solid plate model with Nitsche's contact reads in this case:

8 > > > > > > > > > > > > > > > < > > > > > > > > > > > > > > > : Find (u, ✓) 2 V SP such that for all (v, ) 2 V SP Z ⌦ N ↵ ↵ (v) + M ↵ ↵ ( ) + T (@ v 3 + ) + Q @ 3 + P 3 d⌦ 1 r Z ⌦ [ 33 (u) + r(g + u 3 "✓ 3 )] (O 33 (v) + rv 3 " 3 ) d⌦ + 1 r Z ⌦ P B(0,F [ 33 (u)+r(g+u 3 "✓ 3 )] ) ( ↵3 (u)e ↵ + ru t ) • (rv t + O ↵3 (v)e ↵ ) d⌦ O r Z ⌦ 33 (u) 33 (v) + ↵3 (u) ↵3 (v) d⌦ = Z ⌦ F i v i G i i d⌦.
As expected, we see that compared to the previous plate models, the solid plate model, due to its enriched transverse displacement, shows both pinch stress and pinch deformation in contact Nitsche terms.

Gustafsson, Stenberg and Videman discrete Nitsche's formulation of the Kirchhoff-Love plate model

In this section, we describe the Nitsche formulation of the Kirchhoff-Love plate model constrained by a rigid obstacle introduced in [START_REF] Gustafsson | A stabilised finite element method for the plate obstacle problem[END_REF]. For simplicity, only the transverse displacement variable (ũ := u 3 : ⌦ ! R) is considered. We address this formulation mainly to compare numerical results we obtain below with the models proposed previously. The variational formulation for a Kirchhoff-Love thin elastic clamped/free plate made of a homogeneous and isotropic material and lying over a rigid obstacle, consists in the following variational inequality:

8 > > > > > > < > > > > > > : Find (ũ, ˜ n ) 2 Ṽ ⇥ ⇤ such that for any ( w, µ) 2 Ṽ ⇥ ⇤ Z ⌦ D ⇥ (1 ⌫) @ 2 ↵ ũ + n ũ ↵ ⇤ @ 2 ↵ w d⌦ h˜ n , wi Ṽ0 , Ṽ = Z ⌦ f w d⌦, hµ ˜ n , wi Ṽ0 , Ṽ 0, (28)
with ˜ n the Lagrange multiplier for the contact condition (the contact reaction), f :=

Z " " f V 3 dx 3
the resulting transverse loading and the bending modulus

D := 2 E " 3 3 (1 ⌫ 2 )
the bending modulus, " > 0 being assumed to be constant all along the plate. The plate is assumed to be clamped on a part D of the boundary @⌦. Then, the space for the displacements is

Ṽ := { w 2 H 2 (⌦) | w(x) = @ n w(x) = 0, 8x 2 D },
where @ n w is the normal derivative along D , and the convex of admissible Lagrange multipliers is

⇤ := {µ 2 Ṽ0 | hµ, wi Ṽ0 , Ṽ 0, 8w 2 Ṽ, w  0}.
The well-posedness of ( 28) is guaranteed by the following result [START_REF] Pozzolini | Vibro-Impact of a plate on rigid obstacles: Existence theorem, convergence of a scheme and numerical simulations[END_REF].

Lemma 5.1. The bilinear form ã : Ṽ ⇥ Ṽ ! R defined by ã(ũ, ṽ) :=

Z ⌦ D ⇥ (1 ⌫) @ 2 ↵ ũ + ⌫ ũ ↵ ⇤ @ 2 ↵ ṽ d⌦,
is a scalar product on Ṽ which is equivalent to the canonical scalar product of H 2 (⌦) defined on Ṽ.

The contact reaction ˜ n (u) = D 2 ũ f is a positive measure, possibly singular, and its support included in the contact zone, denoted ⌦ C . For more details see [START_REF] Pozzolini | A Stability Result Concerning the Obstacle Problem for a Plate[END_REF].

The reaction force between the obstacle and the plate can be approximated by

˜ k n = 1 k [ũ g] + , (29)
where k > 0 is the penalty coefficient. Classically, in the limit k ! 0, ˜ k n tends to ˜ n (u), which means that the obstacle becomes rigid and the problem is reduced to that of constrained minimization [START_REF] Necas | On the solution of the variational inequality to the Signorini problem with small friction[END_REF]. Now, we summarize the technique proposed in [START_REF] Gustafsson | A stabilised finite element method for the plate obstacle problem[END_REF]. This approach is analogous to the derivation of Nitsche's method for 3D-contact problems. Let Th be a conforming shape regular triangulation of ⌦ which we assume to be polygonal, and the associated finite element subspace is Ṽh ⇢ Ṽ. Optimality is achieved when the Lagrange multiplier is approximated with four degrees smaller polynomials than that of the displacement variable, for example, when the Argyris element is coupled with a piecewise linear and discontinuous approximation of the Lagrange multiplier.

A reformulation of the normal contact conditions with a Nitsche-based method consists in giving an explicit expression of the discrete Lagrange multiplier from the previous formulation of contact problem [START_REF] Necas | On the solution of the variational inequality to the Signorini problem with small friction[END_REF]. As mentioned above, the contact reaction is a possibly singular measure, therefore the following L 2 -approximation of the reaction force is considered in [START_REF] Gustafsson | A stabilised finite element method for the plate obstacle problem[END_REF]:

˜ ⇡ h n (u) K := 1 k + rh 4 K h ⇡ h (g) K ⇡ h (ũ h ) K + rh 4 K (⇡ h ( 2 ũh ) ⇡ h ( f )) i + 8K 2 Th , ( 30
)
where ⇡ h is the L 2 -projection on Ṽh , k 0 is the penalty coefficient from ( 29) and r 0 is the stabilization or Nitsche's parameter. This hybrid formulation [START_REF] Poulios | An unconstrained integral approximation of large sliding frictional contact between deformable solids[END_REF] of the Lagrange multiplier corresponds to a pure penalty approach when r = 0 and a pure Nitsche's method if k = 0.

We introduce the function H : ⌦ ! R such that H K = h K , for all K 2 Th . A finite element approximation of (28) based on hybrid Nitsche's method ( 30) is:

8 < :

Find ũh 2 Ṽh such that ãh (ũ h , ṽh ; ũh ) = Lh (ṽ h ; ũh ) 8ṽ h 2 Ṽh , (31)
where

ãh (ũ h , ṽh ; ũh ) := ã(ũ h , ṽh ) + Z ⌦ C ( wh ) 1 k + rH 4 ũh ṽh d⌦ Z ⌦ C ( wh ) rH 4 k + rH 4 D 2 (ũ h ) ṽh d⌦ Z ⌦ C ( wh ) rH 4 k + rH 4 D 2 (ṽ h ) ũh d⌦ Z ⌦ C ( wh ) krH 4 k + rH 4 D 2 2 (ũ h ) 2 ṽh d⌦ Z ⌦\⌦ C ( wh ) rH 4 D 2 2 (ũ h ) 2 ṽh d⌦, Lh (ṽ h ; ũh ) := Z ⌦ f ṽh d⌦ + Z ⌦ C ( wh ) 1 k + rH 4 g ṽh d⌦ Z ⌦ C ( wh ) rH 4 k + rH 4 D 2 (ṽ h ) g d⌦ Z ⌦ C ( wh ) rH 4 k + rH 4 f ṽh d⌦ Z ⌦ C ( wh ) rkH 4 k + rH 4 D 2 (ṽ h ) f d⌦ Z ⌦\⌦ C ( wh ) rH 4 D f 2 ṽh d⌦.
The contact set ⌦ C ( wh ) above is given by

⌦ C ( wh ) := {(x, y) 2 ⌦ | ˜ h n ( wh ) > 0}.
The practical solution algorithm for Problem (31) is a fixed point process where at each step the contact set ⌦ C is approximated using the displacement field from the previous iteration so that System (31) becomes linear. The process is complete when the norm of the displacement field is below a predetermined tolerance. The stopping criterion is formulated with respect to the strain energy norm kwk ã := p ã(w, w) .

Numerical results

The numerical comparison of the different approaches is done on the frictionless contact between a rectangular elastic plate and two different rigid obstacles: a flat and a paraboloid one. The Kirchhoff-Love, Mindlin-Reissner and solid plate models are compared to a reference solution obtained directly from the 3D elasticity model and using an augmented Lagrangian multiplier technique to approximate the contact condition as in [START_REF] Alart | A generalized Newton method for contact problems with friction[END_REF][START_REF] Renard | Generalized Newton's methods for the approximation and resolution of frictional contact problems in elasticity[END_REF]. To avoid as much as possible the shear locking for the reference solution, quadratic elements (Q 2 -Lagrange elements) with two layers of elements in the thickness and a very fine mesh (100 ⇥ 100 ⇥ 2 elements) are used. The main goal is then to determinate which model will be the closest to the 3D reference solution.

Numerical tests are performed on a domain ⌦ " = ] 10, 10[ 2 ⇥ ] ", "[, for different values of the thickness " and subjected to different distributed transverse loads f . The contact surface C is still the lower part of ⌦ " (i.e. C := " = {x 2 @⌦ " : x 3 = "}). For the sake of simplicity, the plate is clamped along all its lateral sides. Because of the symmetry of this problem, only a quarter of the geometry is taken into account, where two of the lateral boundaries are clamped (see Figure 2) . We assume that the plate is made of an isotropic elastic material with a Young modulus of E = 200 GPa and a Poisson ratio ⌫ = 0.3. Following the justification given in Section 7, Proposition 7.1, Nitsche's parameter is chosen equal to r = E "

. Moreover, only the symmetric variant O = 1 is considered. The relative error is given by the ratio

u h u ref ku ref k ,
with L 2 norm or H 1 semi-norm on ⌦ " . This means in particular that the norms are applied on a 3D reconstruction of the 2D solutions obtained from the different plate models.

The discrete contact problem is solved using a generalized Newton method on our open-source finite element library GetFEM++ (see [START_REF] Renard | Automated FE modeling of multiphysics problems based on a generic weak form language[END_REF] and http://getfem.org).

Finite element approximation of the plate models

The numerical convergence tests will be performed on a family of quadrangulations T h of the midplane domain

⌦ such that ⌦ = [ K2T h K. We denote h K the diameter of K 2 T h and h = max K2T h h K .
For the Kirchhoff-Love and Mindlin-Reissner models, problems can be decoupled in a plane elasticity problem (membrane displacement) and a transverse displacement problem (bending problem). The membrane displacements are approximated with standard Q 1 -Lagrange quadrilateral elements (also called QUAD4 element).

The bending problem of the Kirchhoff-Love model corresponds to a fourth-order partial differential equation. Thus, a conformal finite element approximation space must be a subspace of H 2 (⌦) which requires the use of C 1 -finite elements (C 0 -interior penalty method such as the one presented in [START_REF] Brenner | A quadratic C 0 interior penalty method for the displacement obstacle problem of clamped Kirchhoff plates[END_REF] can be an interesting alternative). Among the available elements having this regularity (see [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF]), we chose the reduced FVS (Fraeijs de Veubeke-Sanders) quadrangles for our numerical tests. The FVS element basis functions are piecewise P 3 polynomials on a division into four sub-triangles of the quadrilateral, with C 1 -continuity across each internal edge. In addition, to decrease the number of degrees of freedom, the normal derivative is assumed to vary linearly along the external edges of the elements (this assumption does not hold on the internal edges). Consequently, this element has the same computational cost in terms of degrees of freedom as the standard first-order elements used for Mindlin-Reissner bending problem (3 degrees of freedom per vertex, the deflection and the two section rotations).

For Mindlin-Reissner and solid plate models, the transverse displacement, the sections rotations and the pinching are approximated with Q 1 -Lagrange quadrilateral elements. A very important issue is that this element is subject to the so-called shear locking so that a direct Galerkin procedure does not give a satisfactory approximation. There are several ways to circumvent the shear locking problem such as reduced integration and MITC projection. We use here the MITC projection (Mixed Interpolation of Tensorial Components) which was first developed to reduce shear locking for the 4-node continuum mechanics based MITC4 element [START_REF] Bathe | A four node plate bending element based on Mindlin-Reissner plate theory and mixed interpolation[END_REF]. The error analysis [START_REF] Bathe | On the convergence of a four-node plate bending element based on Mindlin-Reissner plate theory and mixed interpolation[END_REF] performed on this element showed that the MITC4 element is optimally convergent for deflections and rotations on regular meshes and static analysis.

The solid plate model is subjected to a supplementary locking between pinching and membrane deformation due to the term

2"⌫E (1 + ⌫)(1 2⌫) ⇣ ⌘⌘ (u) + ✓ 3 ⌘⇣ ⌘⌘ (v) + 3 ⌘ ,
in the weak formulation [START_REF] Gustafsson | A stabilised finite element method for the plate obstacle problem[END_REF]. This supplementary locking is treated with a reduced integration technique, in the sense that a single Gauss point is used for this term.

6.2 Numerical results for a very thin plate, a flat obstacle and a uniform volume load A very thin plate with " = 0.01m is first considered with a uniform volume load f V = 1Mpa and a flat obstacle represented by the initial gap g = 0.1 ". The reference 3D solution is represented in Figure 3 as an illustration.

Figure 4: Relative L 2 (⌦ " )-norm of the error in displacement with respect to the mesh size for the different plate models, comparatively to the 3D solution. Very thin plate on a flat obstacle with a uniform volume load.

Figure 5: Relative H 1 (⌦ " )-norm of the error in displacement with respect to the mesh size for the different plate models, comparatively to the 3D solution. Very thin plate on a flat obstacle with a uniform volume load.

The convergence test for the different models is represented on Figure 4 for the L 2 (⌦ " )-norm of the displacement and on Figure 5 for the H 1 (⌦ " )-semi-norm. In these figures, and in the rest of the paper, "GSV" refers to the method of Gustafsson, Stenberg and Videman (with the penalty coefficient k = 0), "KL" refers to the model with Kirchhoff-Love assumptions of Section 3.1, "MR" to the model with the Mindlin-Reissner assumptions of Section 3.2, "MR no ps" to the model of Section 4.1 and "SP" to the solid plate model of Section 4.2.

As the plate is very thin, the plate models give some good approximation of the 3D solution even for small mesh sizes. This results in a decrease in error for small mesh sizes. One can see that the curves start to slow down their decrease for h = 0.1, which is the smallest mesh size used.

Kirchhoff-Love models are a little more accurate for coarse meshes but are outpaced by Mindlin-Reissner and solid plate models for finer meshes. For very fine meshes, the most ac-curate models are those that do not take into account the plane stress assumption ("MR no ps" and "SP").

Numerical results for a thicker plate, a flat obstacle and a very large Gaussian volume load

We now consider a thicker plate with " = 0.05, still with a flat obstacle defined by the gap g = 0.15 and subjected to a very large Gaussian volume load f V (x) = 50000 exp( (x 2 1 + x 2 2 )) MPa. A large uniform load has not been considered because it concentrates the deformation close to the boundaries where the plate is clamped. The reference solution is plotted in Figure 6. The convergence performances are studied for L 2 (⌦ " )-norm in Figure 7 and for H 1 (⌦ " )-seminorm in Figure 8. Of course, in this experiment and due to the very large load, the pinch deformation is more important. It is again observed on these curves that the solid plate model "SP" is the most accurate on very fine meshes.

An interesting comparison here is the one between the two models based on Kirchhoff-Love assumptions. The main difference is that in the "GSV" model with the penalty parameter k = 0, the contact condition is prescribed by a Nitsche method which is equivalent to a stabilized Lagrange multiplier method with static condensation of the Lagrange multiplier. This of course guarantees a stricter contact condition than the "KL" model presented in Section 3.1 which corresponds to a penalized contact condition. Figure 8: Relative H 1 (⌦ " )-semi-norm of the error in displacement with respect to the mesh size for the different plate models, comparatively to the 3D solution. Thicker plate on a flat obstacle with a Gaussian volume load. One possible interpretation on the fact that "KL" model is far more accurate that "GSV" model here is that the penalized contact condition with a penalty coefficient E/" close to the stiffness of the elastic material gives a more realistic contact condition in the sense that it better reflects the contribution of the pinching of the plate. This interpretation is supported by the curves in Figures 9, 10 and 11. On these Figures, the transverse displacement is drawn close to the contact zone. The difference between the three Figures is that the transverse displacement is taken respectively on the bottom of the plate (contact boundary " ), on the mid-plane and on the top of the plate respectively. Of course, for the models based on the Mindlin-Reissner or Kirchhoff-Love assumptions, the three curves are identical since the transverse displacement is constant across the thickness. A difference is visible only for the "SP" model and the reference solution. We can see the very good behavior of the "SP" model which is close to the reference solution on the three Figures with practically no penetration on the contact boundary " (see Figure 9), conversely to the other models. We can see also that most of the models give a good approximation of the deformation on the mid-plane, except the "GSV" model since it corresponds to a quasi-strict non-interpenetration condition when k = 0. We can conclude that for an important pinching of the plate, the "SP" model leads to a more realistic approximation of the contact condition and conversely, the "GSV" is less accurate due to the more strict non-interpenetration condition.

Numerical results with a paraboloid obstacle and no volume load

For this last numerical case, we still consider a plate of medium thickness " = 0.5, on ⌦ = [ 1, 1] ⇥ [ 1, 1], without load (f V = 0), a paraboloid obstacle represented by the initial gap g = x 2 1 + x 2 2 and a nonzero transverse displacement u " 3 = 0.4 prescribed on the boundary " D .

Figure 12: Thicker square plate (" = 0.5m) on a paraboloid obstacle: 3D reference solution with the augmented Lagrangian method and Q 2 -Lagrange elements.

Figure 13: Relative L 2 (⌦ " )-norm of the error in displacement with respect to the mesh size for the different plate models, comparatively to the 3D solution. Thicker plate on a paraboloid obstacle with no volume load.

The reference solution is shown in Figure 12. In that case, one would expect to be closer to a plane stress situation since the only force acting on the plate is the contact force which acting primarily at the frontier of the effective contact area.

The convergence test is shown in Figures 13 and14. The numerical results are similar to those of the previous case with the flat obstacle. There is still a big difference between the "GSV" model and the "KL" ones although they are based on the same kinematic assumptions (but with two different approximations of the contact condition), and the "SP" model is still the one offering better accuracy for fine meshes although the plate pinch is a priori limited for this case. Figure 14: Relative H 1 (⌦ " )-semi-norm of the error in displacement with respect to the mesh size for the different plate models, comparatively to the 3D solution. Thicker plate on a paraboloid obstacle with no volume load.

Nitsche's parameter setting in the solid plate case

A very large value of Nitsche's parameter can deteriorate the solvability of the numerical approximation, especially when a Newton method is used. Conversely, at least for O 6 = 1, a too small value of Nitsche's parameter will not ensure the coercivity of formulation [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF] and may lead to non-existence or even spurious solutions (see [START_REF] Chouly | Symmetric and non-symmetric variants of Nitsche's method for contact problems in elasticity: Theory and numerical experiments[END_REF]). In this section, we present a mathematical justification of the choice of a Nitsche parameter r close to E ✏ as we used it in numerical tests.

For the sake of simplicity, this justification is produced only for the solid plate model (although it can be adapted to other models) and for Tresca's friction instead of Coulomb's friction. Tresca's friction is a simplification of Coulomb's friction in which the friction threshold is assumed to be known and does not depend on the normal stress.

The following results give an estimate of the optimal Nitsche's parameter that ensures the coercivity of formulation [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF] in the solid plate case. Proof. The monotonicity of ( 14) for Tresca friction with a friction threshold (see [START_REF] Chouly | An adaptation of Nitsche's method to the Tresca friction problem[END_REF] for Nitsche's method for Tresca friction) is ensured with kinematic assumptions (1) whenever for two arbitrary displacements u 1 , u 2 2 V " satisfying the kinematic assumptions (1) and w = u The boundary " D being assumed to be of non-zero measure and using Korn's inequality (see [START_REF] Chouly | Symmetric and non-symmetric variants of Nitsche's method for contact problems in elasticity: Theory and numerical experiments[END_REF]), the monotonicity is ensured for r > r 0 and r 0 satisfying + "E 1 + ⌫ (@ w 3 + ✓ )(@ w 3 + ✓ ) + " 3 E 3(1 + ⌫) (@ ✓ 3 )(@ ✓ 3 )

+ 2"E 1 + ⌫ ✓ 2 3 , | (w)n| 2 x 3 = " = E 2 (1 + ⌫) 2 ✓ ✓ 3 + ⌫ 1 2⌫
( ⌘⌘ (w) " ⌘⌘ (✓) + ✓ 3 )

◆ 2 + 1 4
E 2 (1 + ⌫) 2 (@ w 3 "@ ✓ 3 + ✓ ) (@ w 3 "@ ✓ 3 + ✓ ) Using Young's inequality, we write for 1 > 0 ✓

✓ 3 + ⌫ 1 2⌫
( ⌘⌘ (w) " ⌘⌘ (✓) + ✓ 3 )

◆ 2  (1 + 2 1 )✓ 2 3 + (2 + 1 1 )
⌫ 2 (1 2⌫) 2 (( ⌘⌘ (w) + ✓ 3 ) 2 + " 2 ( ⌘⌘ (✓)) 2 ), and for 2 > 0 (@ w 3 "@ ✓ 3 + ✓ ) (@ w 3 "@ ✓ 3 + ✓ )  (1 + 2 ) (@ w 3 + ✓ ) (@ w 3 + ✓ ) + (1 + 1

2

)" 2 @ ✓ 3 @ ✓ 3 .

Then, inequality [START_REF] Renard | Generalized Newton's methods for the approximation and resolution of frictional contact problems in elasticity[END_REF] is satisfied for r 0 verifying r 0

(1 + O) and 2 = 3. The existence and uniqueness of a solution can be deduced with the same arguments as in [START_REF] Chouly | Symmetric and non-symmetric variants of Nitsche's method for contact problems in elasticity: Theory and numerical experiments[END_REF]. ⇤

Note that the monotonicity of Tresca's friction law is crucial to prove the previous result. Its adaptation to Coulomb's friction remains an open problem.

Conclusions

We presented a systematic way to derive Nitsche's formulations for different kind of kinematics of plate models and studied these approaches in the context of finite element discretization. Several numerical studies were performed to illustrate the behavior of the proposed models. The numerical results allow us to draw the following main conclusions:

• For standard Kirchhoff-Love and Mindlin-Reissner plate models, which include a plane stress assumptions, the application of Nitsche's methods, as we propose in Section 3, simply leads to a penalized contact condition (although it is not the case for Gustafsson-Stenberg-Videman approach).

• The Gustafsson-Stenberg-Videman approach developed in [START_REF] Gustafsson | A stabilised finite element method for the plate obstacle problem[END_REF] for Kirchhoff-Love plate model leads to a stricter non-interpenetration condition when its penalty coefficient k is set to 0.

• The numerical comparison with a fully 3D solution give a clear advantage to the penalized contact condition compared to a strict non-interpenetration condition for standard plate models. The presented Nitsche's approach for the Kirchhoff-Love and Mindlin-Reissner plate models can then be considered as a kind of justification of the use of a penalized contact condition for plate models since it allows to recover a certain elasticity in the transverse direction for the use of a convenient penalty coefficient (r = E/"). Of course, The Gustafsson-Stenberg-Videman approach with a similar penalty coefficient would give some similar results.

• We proposed a solid plate model consisting in a kinematics enriched by a pinching component combined with Nitsche's method for the contact condition. It combines the advantages of a good approximation and a stricter non-interpenetration condition on the contact surface. This model seems numerically to be richer for capturing the effects of contact stresses in the thickness of plate and reproduce the real 3D behavior of the structure in contact with an obstacle.

Figure 2 :

 2 Figure 2: Square plate under transverse load, only a quarter of the plate (blue area) is modeled.

Figure 3 :

 3 Figure 3: Very thin square plate (" = 0.01m) under uniform volume load and flat obstacle: 3D reference solution with the augmented Lagrangian method and Q 2 -Lagrange elements.

Figure 6 :

 6 Figure 6: Thicker square plate (" = 0.05m) under very large Gaussian volume load and flat obstacle: 3D reference solution with the augmented Lagrangian method and Q 2 -Lagrange elements.

Figure 7 :

 7 Figure 7: Relative L 2 (⌦ " )-norm of the error in displacement with respect to the mesh size for the different plate models, comparatively to the 3D solution. Thicker plate on a flat obstacle with a Gaussian volume load.

Figure 9 :

 9 Figure 9: Lower-plane displacement on lateral edge: interpenetration with the rigid obstacle.

Figure 10 :

 10 Figure 10: Mid-plane displacement on lateral edge.

Figure 11 :

 11 Figure 11: Upper-plane displacement on lateral edge.

Proposition 7 . 1 .Remark 7 . 2 .

 7172 For the solid plate model of Section 4.2 and for special case n = n = e 3 and with Tresca's friction instead of Coulomb's friction, a Nitsche's parameter value r > ⇣(⌫, O) of the weak form (14) with kinematic assumptions (1) and consequently the existence and uniqueness of a solution. The coefficient ⇣(⌫, O) is close to unity for O = 1 and a Poisson ratio ⌫ = 1/3. It goes to infinity in the incompressible limit.

1 u 2 ⌦| (w)n| 2 x 3 = " d⌦ ↵ku 1 u 2 k 2 1 ,⌦ " . ( 32 )x 3 = " d⌦ ↵kwk 2 1

 2313232 one has 8 > > > > > > > > > > > < > > > > > > > > > > > : (0, ) ( t (u 1 ) ru 1 t ) P B(0, ) ( t (u 2 ) ru 2 t )) • (r(w t ) O t (w)) x 3 = " d⌦ O r ZUsing as in ([START_REF] Chouly | An adaptation of Nitsche's method to the Tresca friction problem[END_REF], Theorem 3.3) the monotonicity of the projections and Schwarz's inequality, we obtain the sufficient condition Z

Acknowledgments

The work presented in this paper has been carried out with the generous support of ESI Group.


J Guidez
email: joel.guidez@cea.fr

E Merle

D Heuer

S Bourg

G Campioni

M Allibert

S Delpech

P Gauthé

A Laureau

J Martinet

J Serp

Molten Salt Reactor to close the fuel cycle: example of MSFR multi-recycling applications

Keywords: Multi-recycling, waste reduction, fast reactors, Molten Salt Reactor, breeder and burner

On October 1, 2018, 452 nuclear water reactors were operating in the world mainly of second generation (PWR or BWR). The first third generation reactors, EPR and AP1000, have just started in China, and 2 two other ones should start next year in Europe. All these reactors produce spent fuel. Some countries have decided to manage this spent fuel as a final waste. In France, a reprocessing of this fuel is carried out in order to separate the uranium and the plutonium remaining in this spent fuel (96% of the total mass of the spent fuel) and reused as MOX fuel. Only 4% of the spent fuel is sent in to the final waste. Reprocessed uranium is stored as a strategic resource whereas plutonium is reused once in MOX fuel. Today, only a mono-recycling of the Pu plutonium is implemented, because of the difficulties to use it a second time in the REPPWR, mainly due to its isotopic composition. Only a reactor with a fast neutron spectrum, for example a Sodium Fast Reactor (SFR), makes the multi-recycling of the plutonium efficient by taking advantage of a higher ratio of fission vs over capture. This strategy was validated by reprocessing MOX assemblies from the French Phenix SFR, and by remanufacturing fresh fuel from the recycled Pu that was irradiated again in Phenix. This experiment demonstrates the feasibility of this actinides recycling by fast reactors. Molten Salt Reactors (MSR) with a fast neutron spectrum would also be potentially able to burn this waste while producing energy. The final performance obtained by burning plutonium or final waste is fairly comparable for any fast reactor. But on the other hand, the MSRs would allow this operation to be achieved in one step, through a continuous multi-recycling and without cumbersome operations to make new fuel assemblies, which entails considerable time and money savings. A concept of fast MSR of 3000 MWt, the MSFR (Molten Salt Fast Reactor) has been studied in France since 15 years by the CNRS and then further developed within European projects. The main technical choices for this reactor, used for this simulation, are explained. Several options are given on the way to use the MSFR as a burner of final products today available after reprocessing. In the first case, the reactor is used as burner of available depleted uranium and plutonium oxide recovered from the reprocessing. These products are fluorinated and introduced into the reactor to achieve criticality. The evolution of the composition of the salt during the operation is given as well as the periodic additions made, and the corresponding consumption by TWh produced. The evolution of the composition of other specific elements that could alter the physico-chemical properties of the salt or the neutronic performances of the reactor (such as the lanthanides) are also given since a specific treatment of the salt has to be developed to extract them. The same exercise is repeated in a second case with the assumption that a reprocessing of a spent MOX has been made that gives a new composition of the plutonium. In Finally in a third option, all the minor actinides (americium, curium and neptunium) are recovered and recycled together, with uranium and plutonium. The possible consumption of waste by TWh produced is determined. Changes in the composition of salt over time give also an idea of the quantities of products to be extracted by salt treatment. In conclusion, this paper gives quantitative estimations of the possibilities to operate a MSR in the U/Pu fuel cycle using only the available products obtained by reprocessing the spent nuclear fuel of a fleet of second-or third-generation water reactors.

1) Introduction

The aim of this paper is to study the neutronic operation of a molten salt reactor using the products available at the La Hague reprocessing plant, as part of the current operating cycle of the French fleet. This simulation is based on the concept of molten salt fast reactor, the MSFR, initially defined at the CNRS and subject of the EVOL and SAMOFAR European projects. The main parameters of the reactor are specified in this paper.

We are in the context of neutronic evolving calculations to verify the criticality of the concept, the consumptions made, and the necessary reloads. The evolutions of the main products are given, including the products to be extracted during operation by reprocessing.

2) Review of the French fuel cycle and related issues

The enrichment of the uranium needed to manufacture the fuel of the 58 PWRs of the French fleet leads to an annual production of about 7000 tons of depleted uranium (depU) per year. About 340 000 tons are now stocked, available and unemployed.

The reprocessing of the PWR fuel leads to an annual production of about 10 tons of Pu and 950 tons of reprocessed uranium (URT). These 10 tons of Pu, mixed with 110 tons of depU allow manufacturing annually 120 tons of MOX fuel used today in twenty-two 900MW PWR reactors. Until 2013, the URT was re-enriched to produce a fuel called URE, used in four 900MW reactors. The 800 tons of URT / year without employment were stored. The resumption of URE fuel fabrication is announced for 2023.

The reprocessing of used MOX is currently not carried out because the uranium and plutonium that would be obtained would not have the isotopic quality allowing their effective use in PWRs. They are currently stored, waiting for reprocessing for a fourth-generation reactor capable of burning them.

It is known that fast neutron reactors are able to burn these products after separation, and in particular fast sodium reactors whose industrial feasibility is demonstrated. In this paper, we look how a MSR could theoretically perform this same operation.

Several cases of calculations are looked at using the different available materials, i.e. the depleted uranium, the uranium and the plutonium obtained after the reprocessing of the PWR assemblies, the uranium and the plutonium obtained after reprocessing of the used MOX, as well as the actinides obtained in this reprocessing if one decided to separate them.

3) Methods and parameters used for these studies

Code used

Molten salt reactors have different ways of operation and management than solid fuel reactors. Tools have been specifically designed to simulate these reactors within the CNRS /LPSC. The code used for the studies presented in this paper is the coupling of a probabilistic neutronic code to a material evolution code developed to be adapted to the operating characteristics of a liquid fuel reactor (REM code [2]). This tool, developed over the past twenty years, has already been validated on the MSRE and MSBR historical concepts [START_REF] Mathieu | The Thorium Molten Salt Reactor: Moving on from the MSBR[END_REF] and then via a European benchmark on the MSFR reactor [START_REF] Brovchenko | Neutronic Benchmark of the Molten Salt Fast Reactor in the frame of EVOL and MARS Collaborative Projects[END_REF]. This code works by assuming an addition of plutonium and uranium during reactor operation, via a salt control unit, to maintain criticality.

In these neutronic calculations, an extraction of the fission products is simulated. The gaseous fission products (He, N, Ne, Ar, Kr and Xe), the noble metals (Nb, Mo, Tc, Ru, Rh, Pd, Ag) as well as Sb and Te are extracted in a very short time (about 30 s) either by bubbling or by deposition on a trap for metals. The lanthanides, Y, Zr and the alkaline and alkaline earth elements (Na, Mg, K, Ca, Rb, Sr, Cs, Ba, Fr, Ra) are extracted with an hypothesis of characteristic time of 1 year.

Parameters of the MSFR concept of studies

The power considered for the "reference MSFR" used here is 3000 MWth, which can produce about 1400 MWe [START_REF] Allibert | Chapter 7 -Molten Salt Fast Reactors[END_REF][START_REF] Heuer | Towards the Thorium Fuel Cycle with Molten Salt Fast Reactors[END_REF]. This is the power retained for most paper reactors under study (Russian concept MOSART [START_REF] Afonichkin | Molten Salt Actinide Recycler and Transforming System with and without Th-U Support: Fuel Cycle Flexibility and Key Materials Properties[END_REF], etc. ..) to allow comparisons on a typical case. The MSR concept used for the calculations presented in this article has been optimized as breeder in the thorium fuel cycle. The idea here was to check the ability of a fast spectrum MSR to both burn and breed in the U/Pu fuel cycle before further studies on a MSR dedicated to such an application. The points that will be optimized are thus mentioned all along the article as they appear in the studies.

-Choice of salt

The fuel salt is circulating in the MSFR and acts also as the coolant. This places specific constraints on the liquid used as fuel in such a reactor. There are two main salt families that can be used: chloride and fluoride.

Chloride salts have several qualities: they have a harder spectrum, they are able to operate at slightly lower temperatures and the Pu solubility is higher. For example, they were selected by two American start-ups who are aiming at the incineration of waste (Terrapower and Elysium). However, they have some disadvantages:

-If an enrichment of the natural chlorine in 37 Cl is not carried out, one obtains a significant production of 36 Cl, which is a problematic waste because of its period (300 000 years) and its difficult containment in final waste.

-Lack of experimental feedback on its functioning and on corrosion problems.

-The harder neutron spectrum induces a larger amount of irradiation damages (dpa) on the structures.

Nevertheless the choice of salt for such applications is under working and was not optimized in this study.

Fluoride salts have been used for the studies presented here, as for the Russian project MOSART [START_REF] Ignatiev | Molten-salt reactor as necessary element for the closure of the nuclear fuel cycle for all actinides[END_REF], with the corresponding experience feedback available from the MSRE. In this family, several salts are available. In particular the FLiBe salt with a certain percentage of beryllium allowing lower operating temperatures. In our calculations, a pure LiF (Lithium Fluoride) salt was considered as a solvent, corresponding to the reference MSFR concept,

-

Choice of temperatures and specific powers

For this simulation, the current data of the reference MSFR (333 MWth/m 3 of salt in the core) and a core input / output temperature rise of 100°C with a mean fuel temperature of 700°C were used.

ICAPP 2019 -International Congress on Advances in Nuclear Power Plants

France, Juan-les-pins -2019, May 12 │15

-Solubility U / Pu / Actinides

In the presented calculations, we imposed a proportion of 22.5 mol% of heavy nuclei to be on a eutectic which makes it possible to lower the temperatures. In this salt, the maximum U / Pu solubilities are known. The solubility limit of Pu is of the order of 6 to 8% at reactor operating temperatures [START_REF] Lizin | Joint Solubility of PuF3 and UF4 in Melt of Lithium, Sodium and Potassium Fluorides[END_REF].

-Reprocessing

During the reprocessing the extractions are driven by the redox potentials, which imposes an order of extraction which is not the one desired. Some extracted elements (uranium / actinides...) are then to be returned to the fuel salt. At the end of the extraction, elements such as lanthanides and zirconium will constitute the final waste. In our case, during the reprocessing the uranium will be extracted in the first steps (for return to the fuel salt), requiring a necessary increase of the temperature of the mixture to continue the reprocessing in liquid phase. Here again, an optimization of the salt is necessary, for a fuel reprocessing which is not yet finalized today.

Conclusion

The calculations were made on the basis of the initial parameters of the MSFR, but as mentioned a certain number of parameters could be optimized in the context of the use of the reactor considered here, both as regenerator in the U / Pu cycle and as incinerator. This optimization was not performed as part of this study.

4) Results of calculation case #1

We have chosen to use the widely available depleted uranium, and the plutonium from the reprocessing of spent UOX PWR fuel (also available but used today for the manufacture of PWR MOX).

Other choices, as using URT also available, or only plutonium in a burner option, will be investigated also in the near future. Indeed the presence of fertile uranium in the fuel salt brings a breeding aspect with creation of plutonium, which reduces its consumption and thus the capacity of incineration. The idea here is to study a mixed concept combining both burning and breeding capabilities. It should also be noted that the use of blankets would make it possible to imagine breeder options.

In the simulation code, fertile material (here depleted uranium) and fissile material (here Pu and minor actinides) are added periodically to keep the core critical while maintaining the proportion of heavy nuclei initially chosen (22.5%) which impacts the physicochemical properties of the salt which is also the coolant.

For the same reason, the code continuously extracts a number of neutron poisons (lanthanides, zirconium ... see section 3b), which will be the final waste.

The isotopic composition of depleted uranium is 0.25% of 235 U The isotopic composition of Pu is that obtained after reprocessing a spent UOX PWR fuel, after arbitrary but classical aging of 10 years, separation and then aging for 5 years. The isotopic composition obtained is as follows:

5) Results of calculation case #2

In this second study, we start from the assumption of a reprocessing of the used MOX available today where we do not separate the minor actinides from the fission products. So we recover only U and Pu, with a degraded isotopy but usable in the fast MSR. Here the reactor will be started up and fed with depleted uranium as fertile material and with Pu as fissile material with the isotopic vector of Pu from spent MOX, aged 10 years then separated and aged 5 years (leading to the production of 241 Here again the parameters of the data could be changed in future systematic studies, either on the choice at the level of uranium (one could have taken the uranium resulting from the reprocessing, or the URT, ...), or on the time of MOX aging. At startup, the fuel contains 30.2 tons of depleted uranium and 13 tons of Pu. This composition has been adjusted to have an exactly critical core. Figure 4 then gives the evolution of the heavy nuclei compositions in the fuel salt for this second case.

The evolution of plutonium inventories is presented in Figure 5. The conclusions are identical to that of the previous study.

6) Results of calculation case #3

For this third study, we consider the theoretical case where the reprocessing of used MOX would take place with separation of the minor actinides. In this case, the fissile material could be a mixture of Pu and MA coming from used MOX fuel. It should be noted that this option corresponds to the implementation of a minor actinide/fission products separation technique that is feasible from the research results, but is not in place today. The main interest of such a separation is to reduce the toxicity of final waste, down to a period of 300 years, after which time this toxicity would correspond to that of natural uranium. This could help revisiting some long-term storage options and also improving social acceptance.

We then look in our calculations how these actinides can be integrated into the MSFR cycle. Here the reactor will be fed with a fertile material still essentially of depleted uranium and as fissile material the fissile isotopes of Pu and minor actinides (isotopic vector from a used MOX fuel aged 15 years). The isotopic composition obtained for the fissile material is as follows (in atomic percentages): At start-up, the fuel contains 28.5 tons of depleted uranium, 13.4 tons of Pu, 100 kg of Np, 2.1 tons of Am and 380 kg of Cm. This composition was adjusted to have an exactly critical core by fixing the proportion of heavy nuclei (U, Pu and MA) to 22.5 mol% and respecting the isotopic compositions given above. Here again other calculation hypotheses will be considered in coming systematic studies, such as the choice of uranium or the quantities of minor actinides to be injected or the cooling time of spent MOX fuel. Figure 6 shows the evolution of the heavy nuclei compositions in the fuel salt of the fuel circuit. The evolution of the composition of the Pu isotopes is shown in Figure 7. Here again the conclusions are identical to those of the first 2 studies as to the equilibrium of the inventory of heavy nuclei in the fuel salt.

7) Assessment of U / Pu consumption in the three calculation cases presented

Figure 8 shows the depleted uranium consumption in the three cases presented. We see that this consumption is identical in all three cases. This is normal because it corresponds to the injection of the fertile matter that does not affect the criticality level required but is led by the power produced by the reactor. We therefore have an annual consumption of depleted uranium of about 1.15 tons for an MSR of 1400 MWe. Figure 9 gives the consumption of corresponding Pu in the three cases studied. It can be seen that the initial isotopic composition of Pu plays an important role in the establishment of the neutron equilibrium in the core and that in the case of Pu from a PWR fuel a higher Pu supply is required in the first years. The Pu supply is higher in the first years for the first case (Pu of the used UOX) compared to the other two cases using Pu from spent MOX because the isotopic vector of Pu is degraded in the first case given the high proportion of 239 Pu injected, element consumed quickly. In both cases based on Pu from MOX, the presence in the reactor does not degrade the isotopic vector similarly, on the contrary the presence of minor actinides leads to the production of fissile elements by neutron capture.

On the other hand, the feeding slopes are approximately identical at about 32 kg / year. This amount of external Pu used to feed the reactor, whose composition is different for the 3 cases considered but which is very small compared to the whole Pu consumption per year (around 1 ton per GWe mainly produced from 238 U) explains why the Pu isotopic compositions tend to the same values at equilibrium in the 3 calculations (see figures 3, 5 and 7).

8) Summary of extractions

To maintain criticality, the code assumes a continuous extraction of lanthanide or zirconium type poisons. It is a "theoretical" extraction that does not correspond to a reprocessing scheme defined and validated. In particular, the extraction coefficients are "ideal" at 100% which is never the case. Figure 12 gives the evolution of the compositions in the fuel salt, with the extraction hypotheses retained. It can be seen that the concentrations of the elements which are extracted (lanthanides, zirconium, alkalis, etc.) are stabilized quite rapidly. Figure 13 gives the quantities of products extracted from the mixture during the operation of the reactor.

We see an extraction roughly equivalent in mass to the reactor feed. But the elements obtained (metals, lanthanides, alkaline earths, zirconium...) are final wastes with a significantly reduced activity. One can notice the reduction of the toxicity duration of the final waste. Figure 14 gives a calculation of the evolution over time of the radiotoxicity of this residual waste. The calculation is carried out after 30 years of output from the reactor. This waste is at a level of toxicity corresponding to natural uranium, after about 300 years.

9) Conclusions

The neutronic calculations presented were carried out with the initial parameters of the existing "reference MSFR" without further optimization of them (choice of salt, specific power, etc.) vis-à-vis the mode of use of the reactor in the U / Pu cycle.

Three calculation cases are presented using depleted uranium with successively Pu from a reprocessed PWR used UOX fuel, then Pu from a used MOX, finally Pu and minor actinides from a used MOX fuel. The results show the great neutronic flexibility of use of this type of reactor with in all three cases a consumption of natural uranium of about 1.15 t / year and about 30 kg of plutonium per year. So we are almost in iso-breeder mode. Further studies are led to use plutonium without uranium to optimize a Pu burner version of the concept.

As Lavoisier said, "nothing is lost, nothing is created, everything is transformed". It is therefore necessary to extract quite the same quantities of final waste, i.e. about 1.15 tons per year of metals, alkaline earths and lanthanides.( in fact with E=m C2 , calculations show that only about one kilogram /year is converted in energy) It should be noted that this extraction is made theoretically by the code but does not correspond today to a reprocessing process defined and validated. This final waste then has activities and periods much lower. After about 300 years, this activity becomes equivalent to that of natural uranium, which has consequences for the storage options and would also allow better societal acceptability of this final waste.

Figure 1 :

 1 Figure 1: Cycle of materials in the French 400 TWh/year fleet of pressurized water reactorsIt may be noted that this actinides separation could allow a closure of the cycle with final waste whose activity would be strongly reduced.

ICAPP 2019 -

 2019 International Congress on Advances in Nuclear Power Plants France, Juan-les-pins -2019, May 12 │15 main balances are put in place before the first twenty years.

Figure 2 :Figure 3 :

 23 Figure 2: Evolution of the heavy nuclei inventory for the case using Pu from used UOX

Figure 4 :Figure 5 :

 45 Figure 4: Evolution of the heavy nuclei inventory for the case using Pu from used MOX fuel

Figure 6 :Figure 7 :

 67 Figure 6: Evolution of the heavy nuclei inventory for the case using Pu + MA from MOX used fuel

Figure 8 :Figure 9 :

 89 Figure 8: Feeding in fertile matter (depleted uranium) for the three studies presented

Figure 10 :

 10 Figure 10: Breeding ratio of the core for the 3 studied cases

Figure 11 :

 11 Figure 11: Burning rate of the reactor for the 3 cases studied

Figure 12 :Figure 13 :

 1213 Figure 12: Molar proportions of the different elements in the fuel salt

Figure 14 :

 14 Figure 14: Time evolution of the final waste radiotoxicity

Acknowledgment

The authors wish to thank the NEEDS French Interdisciplinary program , the IN2P3 department of the National Centre for Scientific Research (CNRS), Grenoble Institute of Technology, and the European program SAMOFAR (H2020) for their support.


A HIGH-ORDER MOC INCLUDING A SPATIAL POLYNOMIAL EXPANSION FOR CROSS SECTIONS

INTRODUCTION

The TDT (Two and Three Dimensional Transport) deterministic code of the APOLLO3 R calculation platform utilises the Method of Characteristics (MOC) to solve the Boltzmann equation for neutrons. In order to reduce the bias and the errors of the current schemes, in the recent years the TDT application has been extended to three-dimensional extruded geometries [START_REF] Santandrea | A Neutron Transport Characteristics Method for 3D Axially Extruded Geometries Coupled with a Fine Group Self-Shielding Environment[END_REF], to provide a better representation of 3D environment effects also in view of a future one-step core calculation.

Even more recently, a polynomial description for the neutron flux has been introduced [2]: along the axial (extruded) direction a per-region constant ("step-constant") flux representation has been substituted by a description assigning to each region a set of coefficients for the polynomial expansion. At the price of a greater complication in the equations, this has permitted to reduce the number of axial meshes by a factor higher than 10, the memory required by a factor ∼ 3 and the computation time by a factor ∼ 2.

The present work constitutes the natural continuation of the one above, as the same treatment is given to both the neutron flux and the cross sections. The aim is to apply the polynomial expansion in the case of evolution calculations: as a matter of fact, the combination of a polynomial flux expansion and step-constant cross sections over extended axial meshes is highly inaccurate if non-zero burnup cases are considered, since a refined axial mesh is required to describe the modification of cross-section profiles due to nuclide depletion. A polynomial expansion for cross sections is therefore needed, together with the one for the flux, to fully utilize the MOC 3D polynomial method.

THE POLYNOMIAL EXPANSION

If a polynomial profile is assumed for their spatial dependence, the angular mono-group neutron flux ψ( r, Ω) and the arbitrary mono-group macroscopic cross section Σ( r) can be written as

ψ( r, Ω) ∼ P ( r ) • ψ r ( Ω), Σ( r ) ∼ P ( r ) • Σ r P ( r ) = {P p ( r ), p = 0, N p }, ψ r ( Ω) = {ψ r,p ( Ω), p = 0, N p } Σ r = {Σ r,p , p = 0, N p }, (1)
where P is the set of expansion polynomials, and ψ r ( Ω) and Σ r the associated flux and cross-section coefficients, respectively. N p is the chosen degree for the expansion. At least for a first implementation, the choice has been made to use the same expansion base for the two physical quantities, which seems reasonable given their close interdependence. As in [2], the polynomial representation is adopted only along the axial direction, while step-constant values characterize the radial one; moreover, exactly the same polynomials are employed:

P (z r ) = {(z r ) p , p = 0, N p }, zr = z r -zr ∆z r /2 , ( 2
)
where r is the region index, zr the region mid-height axial coordinate, ∆z r the region height and z r the absolute axial coordinate ranging from the bottom to the top of region r. In this way, (z r ) p ∈ [-1, 1] for every region.

In order to highlight how this impacts the numerical problem, one may consider the mono-group integral transport equation evaluated along a characteristic line of direction Ω:

ψ( r 0 + l Ω, Ω) = ψ( r 0 , Ω) e -τ (l) + l 0 dt q( r 0 + t Ω, Ω) e τ (t)-τ (l) , (3)
l being the length of the considered chord and r 0 the chord entry endpoint. According to Eq. ( 1), the optical path length corresponding to t, τ (t), reads

τ (t) = t 0 dt Σ( r(t )) = t 0 dt P (z r (t )) • Σ r = Np p=0 t 0 dt z in r -zr + µt ∆z r /2 p Σ r,p = Np i=0 Λ i0 t + Np-1 i=0 Λ i1 t 2 + Np-2 i=0 Λ i2 t 3 + • • • + Λ 0Np t Np+1 , (4)
where

Λ ij = P i (z in r )Σ r,i+j 2µ ∆zr j 1 j+1 j+i i
, µ is the cosine of the polar component of Ω and z in r is the axial component of r 0 . The angular dependence of the emission source q( r, Ω) is classically expressed in terms of an expansion over real spherical harmonics, up to the chosen anisotropy order K: q( r, Ω) =

K l=0 m=l m=-l A m l ( Ω) q m l ( r) = Nm n=1 A n ( Ω) q n ( r)
, where a single index is used to indicate the angular moments, whose number is N m = (K + 1) 2 . The emission source for group g can therefore be written as

q g ( r, Ω) = = n A n ( Ω) g Σ g →g s,n ( r ) Φ g n ( r ) + 1 k ef f i χ g i g νΣ g f,i ( r ) Φ g 0 ( r ) = n A n ( Ω) g Σ g →g s,n,r P (z r ) ⊗ P (z r ) Φ g n,r + 1 k ef f i χ g i g -→ νΣ g f,i,r P (z r ) ⊗ P (z r ) Φ g 0,r = P 2Np (z r ) • Q 2Np,g r ( Ω), (5)
where

Q 2Np,g r ( Ω) = {Q g r,p ( Ω), p = 0, 2N p }, Q g r,p ( Ω) = n A n ( Ω) g k+l=p Σ g →g s,n,r,k Φ g n,r,l + 1 k ef f i χ g i g k+l=p νΣ g f,i,r,k Φ g 0,r,l . (6)
2/5

In Eq. ( 5) the spatial behaviour of the angular flux moments {Φ g n ( r )}, of the average numbers of neutrons emerging from fission per unit path for each fissile isotope {νΣ g f,i ( r )} and of the coefficients of the scattering cross-section expansion over Legendre polynomials {Σ g →g s,n ( r )} is expressed in terms of the polynomials of Eq. (2). Σ g and Σ i represent the sums over energy groups and fissile isotopes, respectively, and the superscript 2N p indicates that vectors contain the elements corresponding to polynomial orders from 0 to 2N p . Substituting the emission source into Eq. (3) (and omitting the energy index g), this latter results as

ψ( r 0 + l Ω, Ω) = ψ( r 0 , Ω) e -τ (l) + 2Np k=0 P k (z in r ) 2N p p=k Q 2Np r,p ( Ω) p k 2µ ∆z r p-k l 0 dt t p-k e τ (t)-τ (l) , (7)
where

P k (z in r ) = z in r -zr ∆zr/2 k
. The employment of the equation above requires knowing how to evaluate the integrals over the chord length, which, for how the optical path length is expressed (Eq. ( 4)), are of the kind

E n (a 1 , a 2 , . . . , a N p+1 , l) = l 0 dt t n e a 1 (t-l)+a 2 (t 2 -l 2 )+ . . . +a N p+1 (t N p+1 -l N p+1 ) . ( 8
)
Since no analytical solution is available for N p > 1 and since it would be desirable to treat higher-order profiles (at least parabolic), we have decided to adopt a Prony approximation to express e τ (t)-τ (l) as

e a 1 (t-l)+a 2 (t 2 -l 2 )+ . . . +a N p+1 (t N p+1 -l N p+1 ) ≈ N i /2 i=1 C i e b iR x cos(b iI x) sin(b iI x) , (9)
N i being the chosen number of interpolation points. The presence of a cosine or a sine function depends on the b i coefficients: for each pair of complex conjugates, there will be two subsequent terms having the same coefficient real parts (b iR and b i+1R ) and imaginary parts (b iI and b i+1I ), the former with the cosine function and the latter with the sine one.

FROM STEP-CONSTANT TO POLYNOMIAL CROSS SECTIONS

To obtain a first polynomial cross-section representation, we started translating a step-constant mesh into a polynomial one: for this purpose, we considered a set of macro-regions {r}, each one composed of a set of regions {i} having constant cross-section values. To determine the coefficients {Σ r,p , p = 0, N p } of the polynomial expansion within region r, one can first compute the cross-section moments with respect to the polynomials, which by definition read

Σ r = 1 V r r d r P (z r )Σ( r) = 1 ∆z r z max r z min r dz P (z r )Σ(z) = 1 2 1 -1 dz r P (z r )Σ(z r ), (10)
V r being the region volume. The integral over zr is equal to the sum of the contributions due to all the step-constant regions contained in r: therefore, the p th -order moment can be written as

Σ r,p = 1 2 i Σ i zup i zlow i dz r P p (z r ) = 1 2(p + 1) i Σ i (z up i ) p+1 -(z low i ) p+1 , (11)
where

zup i = zi -zr ∆z r /2 + ∆z i ∆z r , zlow i = zi -zr ∆z r /2 - ∆z i ∆z r . (12)
Once the moments have been the computed, the expansion coefficients can be retrieved as

Σ r = P-1 Σ r , (13)
where the elements of matrix P read

P ij = 1 V r r d r P i (z r )P j (z r ) = 1 ∆z r zr+∆zr/2 zr-∆zr/2 dz r z r -zr ∆z r /2 i+j = 1 i+j+1 f or i + j even 0 f or i + j odd . (14)
Applying this to a cell of the ASTRID reactor [3], for an expansion of order 2 we have obtained the results

shown in Fig. 1 for the fissile region (15 cm) and the fertile one above (10 cm). Figures refer to the energy group #1896 in a library of 1968 groups, in the case of burnup = 37500 MWd/t. The choice of this group is due both to the relatively strong variation it undergoes during evolution and to the relatively high values of the first-and second-order coefficients required for the polynomial expansion.

Figure 1: Total cross section plots for the fissile (left) and the fertile region (right).

PERSPECTIVES

The described method is currently being implemented in TDT, where Eq. ( 7) will be used as transmission equation for the trajectory sweep, at the inner iteration level of the power method. The problem will be closed by a suitable formulation of a balance equation, which has not been reported here because of lack of space and that will also be implemented. Once the code supports the polynomial cross-section environment, it will be necessary to accelerate the simulation.

[2] S. Santandrea, L. Graziano, D. Sciannandrone, "Accelerated polynomial axial expansions for full 3D neutron transport MOC in the APOLLO3 R code system as applied to the ASTRID fast breeder reactor," Annals of Nuclear Energy, 113, pp. 194-236 (2017).

[3] B. Fontaine, N. Devictor, P. Le Coz, A. Zaetta, D. Verwaerde, J.-M. Hamy, "The french R&D on SFR core design and ASTRID project," In: Global. Makuhari, Japan (2011).

5/5