A NITSCHE METHOD FOR THE ELASTOPLASTIC TORSION PROBLEM

. This study is concerned with the elastoplastic torsion problem, in dimension 𝑛 ≥ 1, and in a polytopal, convex or not, domain. In the physically relevant case where the source term is a constant, this problem can be reformulated using the distance function to the boundary. We combine the aforementioned reformulation with a Nitsche-type discretization as in Burman et al . [ Comput. Methods Appl. Mech. Eng. 313 (2017) 362–374]. This has two advantages: (1) it leads to optimal error bounds in the natural norm, even for nonconvex domains; (2) it is easy to implement within most of finite element libraries. We establish the well-posedness and convergence properties of the method, and illustrate its behavior with numerical experiments.


Introduction
Problems written with weak formulations involving variational inequalities represent various nonlinear phenomena which occur in mechanics and physics [14,24].We focus on the elastoplastic torsion problem, as presented in, e.g., [18] (see also [7,19]).In the aforementioned reference [18], a direct piecewise affine Lagrange finite element approximation of the variational inequality is also presented, as well as a convergence result (Thm.3.3), and two error estimates in the  1 -norm, in dimension one (Thm.3.4) and in dimension two (Thm.3.5).The error estimate in one dimension is optimal ((ℎ)), whereas it remained suboptimal in dimension two, as it is of order (ℎ problem where the constraint involves the distance to the boundary, so the obstacle is nonsmooth and the usual techniques from the obstacle problem cannot be directly applied: for instance in Theorem 5.1.2 of [11], the obstacle is supposed of Sobolev regularity  2 .In a previous paper [8], a direct finite element approximation of the variational inequality has been proposed, that makes use of piecewise affine, continuous, Lagrange finite elements, and in which the constraint involving the distance function is imposed at each node.When the domain is convex, error estimates have been established in any dimension  = 1, 2, 3, with an optimal error bound of (ℎ), for a regular enough continuous solution.In the case of a nonconvex domain, an error bound of (ℎ 3 4 ) has been proven for a solution of Sobolev regularity   ,  ≥ 7/4.
In the present paper, we propose a new method that combines both the reformulation with the distance function, as in [8], and a Nitsche term that allows to incorporate weakly the inequality constraint, following [6] and related works on Nitsche's method for variational inequalities, see, e.g., [10] and references therein.For this discretization, we manage to derive optimal error estimates, for linear and quadratic finite elements, and even in the nonconvex situation, which improves the result of [8].Moreover, this method is easy to implement into modern finite element librairies, and we provide also some numerical experiments, that allow to confirm the expected theoretical convergence rates.
As usual, we denote by   (•),  ∈ R, the Sobolev spaces.The usual norm of   () is denoted by ‖ • ‖ , , and the corresponding semi-norm is denoted by | • | , .The space  1 0 () is the subspace of functions in  1 () with vanishing trace on .The letter  stands for a generic constant, independent of the mesh size, the value of which can change at different occurences.

The elastoplastic torsion problem
Let Ω ⊂ R  ,  ≥ 1, be an open bounded polytope, connected and with Lipschitz boundary.We consider the variational inequality which, for  = 2, models the torsion of an infinitely long elastoplastic cylinder of cross section Ω and plasticity yield  > 0. To simplify we assume that  = 1.The problem is to find the stress potential  such that where  :  1 0 (Ω) ×  1 0 (Ω) → R is the bilinear form given by: and with  ∈  2 (Ω).The notation  1 represents the nonempty closed convex set of admissible stress potentials: where | • | denotes the euclidian norm in R  .From Stampacchia's theorem we deduce that problem (1) admits a unique solution (see also, e.g., [14,18,19,24]).
Next we suppose that  =  is a constant function (when  is not a constant, there is no clear physical meaning associated to problem (1)).In this case and according to [4] (see also [23]) the problem (1) can be rewritten as follows: find the stress potential  such that with and  Ω denotes the (interior) distance function with respect to the boundary Ω: Note that (2) still admits a unique solution from Stampacchia's theorem.To lighten the discussion we can suppose without loss of generality that  > 0 (see [8], Rem.2.2), so problem (2) can be rewritten as follows: find the stress potential  such that We reformulate (4) using a Lagrange multiplier , and get: We introduce the following notation for the positive part: [] + := max(0, ), for  ∈ R, and recall the relationship for ,  ∈ R.
Following [6], the Kuhn-Tucker condition (5) 3−5 can equivalently be reformulated as with  an arbitrary positive function on the domain Ω.

A Nitsche finite element method
Let   ℎ be a family of Lagrange finite element spaces of degree  ≥ 1 indexed by ℎ, and coming from a family  ℎ of simplicial meshes of the domain Ω (ℎ := max  ∈ ℎ ℎ  where ℎ  is the diameter of  ∈  ℎ ).The family of meshes is assumed regular.More precisely we have: Each simplex  of the mesh  ℎ is supposed to be closed, and we denote by T the interior of  .We define a piecewise polynomial discrete Laplacian as follows, for every  ℎ in   ℎ , and every simplex  ∈  ℎ : The value of ∆ ℎ  ℎ on the facets of the mesh is of no importance, and can be set in practice to 0, for instance.We define also: Observe that, for  = 1, ∆ ℎ  ℎ = 0 and  ℎ ( ℎ ) = .The Nitsche-type method proposed for the discretization of the elastoplastic torsion problem (5) reads: find for all  ℎ ∈   ℎ .Above the notation (•, •) stands for the  2 (Ω)-scalar product, and the function  ℎ is defined cell-wise as follows: where  0 > 0 is the Nitsche parameter.Again, the value of  ℎ on the facets of the mesh is of no importance, and can be set in practice to 0, for instance.
Remark 3.1.As in [9,22], for any parameter  ∈ R, we can write a whole family of methods: Method ( 8) corresponds to  = 0 and can be called an incomplete method, using the terminology widespread for discontinuous Galerkin methods.This method involves less terms and is the easiest to extend to more complex problems [27].A symmetric method is recovered when  = 1, that corresponds to the Galerkin Least Squares technique of [6] and the Nitsche method of [21]: this symmetric method can be recovered thanks to a minimization argument, and the tangent system has a symmetric Jacobian.Provided that the Nitsche parameter  0 is large enough, the analysis below, for  = 0, can be extended without difficulty to other values of .
Remark 3.2.Note that, for the distance function, there holds  Ω ∈  1 (Ω) ∩ C 0,1 (Ω), see [8] and references therein.In [6] the assumption made on the obstacle function is stronger and this function is supposed to be C 1,1 (Ω).In fact, such Nitsche or Galerkin Least Squares formulations do not require so much regularity on the obstacle function.
The following local inverse inequality will be helpful in the sequel, that holds for an arbitrary  ℎ ∈   ℎ and every  ∈  ℎ : where  I > 0 is a constant that depends on the shape regularity of the mesh and of the polynomial order , but not of ℎ and  ∈  ℎ .See, e.g., [2,15] for the proof.

Numerical analysis: well-posedness and error estimate
We first state a preliminary consistency result: The above result is a direct consequence of ( 5)-( 7) and the inclusion   ℎ ⊂  1 0 (Ω).

Well-posedness
We make use of the results from Brezis (see, e.g., [3]) for M-type and pseudo-monotone operators in vector spaces.It consists in showing that the operator associated to problem (8) is bijective.
) admits one unique solution  ℎ ∈   ℎ .Proof.We introduce the nonlinear operator  ℎ :   ℎ →   ℎ defined as follows: for all  ℎ ∈   ℎ and where (•, •) 1,Ω is the scalar product in  1 (Ω).Note that the right-hand side of ( 12) is linear with respect to  ℎ and, hence,  ℎ is well defined thanks to Riesz' Theorem.The existence and uniqueness of solution for problem ( 8) is equivalent to the property of  ℎ to be bijective.To this purpose, according to Corollary 15, p.126 of [3], it suffices to show that by  ℎ is monotone and hemicontinuous.
For the monotonicity we first note that, owing to (12), we have for all  ℎ ,  ℎ ∈   ℎ .On the other hand, by adding and subtracting suitable terms into  1 , we get Hence, using the inequality (6), Cauchy-Schwarz, Young inequalities and the inverse inequality (10), it follows that The monotonicity of  ℎ then follows by inserting this estimate into (13) under the condition  0 ≥  2  2 I .For the hemicontinuity, we must show that the real function  : [0, 1] → R defined as which means that  is Lipschitz and, thus,  ℎ is hemicontinuous.

A priori error estimate
We provide first an abstract estimate, as in [6].
Theorem 4.2.Assume that the solution (, ) ( = ∆ + ) to the elastoplastic torsion problem (5) belongs to  ×  2 (Ω).For  0 > 0 large enough, the approximation  ℎ provided by (8) satisfies the following error estimate: Proof.Let  ℎ ∈   ℎ .We first use the  -ellipticity and the continuity of (•, •), as well as Young's inequality, to obtain: with  > 0 the ellipticity constant.We can transform the last term using the consistency property (11) for  and the finite element formulation (8) for  ℎ .This yields We now transform the expression  ℎ −  ℎ : Hence, by inserting this expression into (16), we get The first term is estimated by using the inequality (6), which yields For the term  2 we use Cauchy-Schwarz and Young's inequality, to obtain There remains to bound the last term above.Let  ∈  ℎ be a mesh cell.Using a triangular inequality, we bound first With the inverse inequality (10) and triangular inequality, we bound Therefore the last term in (19) can be bounded as We combine estimates ( 15)-( 17)-( 18)-( 19)-( 20), which yields Choosing  0 large enough, we obtain the desired bound on the first two terms in (14).The bound on the error on the multiplier  comes from combination of ( 20) and (21).
The optimal convergence of the method for P 1 and P 2 Lagrange finite elements is stated below, for a Sobolev regularity  ≤ 2. Theorem 4.3.Suppose that  = 1, 2, and that the solution  belongs to  1 0 (Ω) ∩   (Ω), with the Sobolev regularity  that satisfies max(1, /2) <  ≤ 2. Suppose that  belongs to  2 (Ω).Suppose finally that the Nitsche parameter  0 is large enough.The solution  ℎ to problem (8) satisfies the following error estimate: with  > 0 a constant, independent of ℎ and , but not of  0 .
Proof.We consider first the case  = 1.We start from the estimate (14).We take  ℎ = ℐ 1 ℎ , the Lagrange interpolant of  onto  1 ℎ .So, first, from standard interpolation estimates, we obtain And then, for  = 1, we can simply proceed as follows: Figure 4.The error between  ℎ and a reference solution as a function of the mesh parameter ℎ for the convex example.The reference solution is obtained using a quadratic finite element method on a suitably refined mesh.
For  = 2, we proceed exactly as above.We start from the estimate (14).Since  1 ℎ ⊂  2 ℎ , we take , the Lagrange interpolant of  onto  1 ℎ , and get the same estimates as above for the three terms.This ends the proof.
The next statement is for P 2 Lagrange finite elements, where a better convergence rate than (ℎ) can be expected if the solution  is regular enough.Theorem 4.4.Suppose that  = 2, and that the solution  belongs to  1 0 (Ω)∩  (Ω), with the Sobolev regularity  that satisfies max(2, /2) <  ≤ 3. Suppose that the Nitsche parameter  0 is large enough.The solution  ℎ to problem (8) satisfies the following error estimate: with  > 0 a constant, independent of ℎ and , but not of  0 .
Proof.We proceed as previously, in Theorem 4.3, but with the following modifications.We start from the estimate (14).We take  ℎ = ℐ 2 ℎ , the Lagrange interpolant of  onto  2 ℎ .Still from standard interpolation estimates, we obtain For the last term, we need the following local error bound, for each  ∈  ℎ : The error between  ℎ and a reference solution as a function of the mesh parameter ℎ for the nonconvex example.The reference solution is obtained using a quadratic finite element method on a suitably refined mesh.
which is possible since  > 2. Then we use standard (local) interpolation error estimates to get: By summation on the simplices  of the mesh  ℎ , we get finally This ends the proof.

Numerical experiments
The following results are computed with the help of scikit-fem [20] for finite element assembly and autograd [26] for automatic differentiation.In particular, the nonlinear problem ( 8) is solved using the semismooth/generalized Newton's method (see, e.g., [13,25] and references therein).The generalized Jacobian is calculated automatically by autograd given the source code for the finite element assembly of the linear form (where  ℎ and  ℎ are in   ℎ )

Numerical convergence rates
We consider two experiments introduced in [19] with convex and nonconvex domains.The parameters are chosen as  =  0 = 10.First, a convergence study is performed using linear elements and the mesh sequences are depicted in Figure 1.Error in the stress potential  is computed against a reference solution which is obtained using quadratic finite elements and a sufficiently refined mesh.The resulting discrete solutions  ℎ and the magnitudes of the gradient |∇ ℎ | are depicted in Figures 2 and 3, and the error in Figures 4 and 5 for the convex and nonconvex examples, respectively.The convergence rate of the  1 error is (ℎ) which is also an expected consequence of Theorem 4.3.Moreover, the convergence rate of the  2 -error is (ℎ 2 ), a rate we are not able to prove.It is well-known that  2 -error estimates are difficult to prove for problems modelled by variational inequalities.
We continue by solving the same experiments using quadratic elements.For this purpose, a reference solution is computed using cubic elements and a suitably refined mesh.The resulting discrete solutions are given in Figure 6 with the convergence rates in Figure 7.This time we observe that the  1 error is approximately (ℎ 1.5 ).This is in agreement with Theorem 4.4.However this is lower than the rate (ℎ 2 ) for quadratic elements and a smooth solution.The reduction in the convergence rate is explained by the Sobolev regularity of the exact solution as  ∈  2 (Ω) while  ̸ ∈  3 (Ω).In one dimension, it is easy to check that the second derivative is expected to have a step discontinuity at the free boundary between elastic and plastic zones which suggests  ∈   (Ω),  < 5/2 (see also Rem. 2.1).This also means that adaptive techniques could be used to improve the convergence rate with respect to the number of degrees-of-freedom.

Comparison with solving the discrete variational inequality
For comparison purposes, we present here results obtained with the discrete variational inequality, introduced in [8].An Uzawa solver is used to handle the constraint  ℎ ≤  Ω on each Lagrange node.First is observed Figure 7.The error between  ℎ and a reference solution as a function of the mesh parameter ℎ for the convex example using quadratic elements.The reference solution is obtained using a cubic finite element method on a suitably refined mesh.that the Uzawa iteration takes really long to converge.Of course this issue can be improved through a better, optimized choice, of the Uzawa parameter, but is not present at all when solving problem (8). Figure 8 depicts the final solutions, that are very similar to those obtained by solving problem (8) (Fig. 3).
Figure 9 presents the convergence rates with the method of [8], still for the nonconvex domain, and for P 1 Lagrange finite elements (the method of [8] is for piecewise affine Lagrange finite elements).The optimal rate of (ℎ) is observed, and is better in this case than the best theoretical rate of (ℎ 3/4 ), expected from Theorem 4.1 in [8].This means that the convergence analysis in [8] may be improved.

Conclusion
The Nitsche method (8) for the elastoplastic torsion problem with a positive constant source term can be implemented easily within many finite element libraries, since it introduces no additional unknown and only requires a standard solver for the nonlinear term that comes from the constraint.When a solver such as semismooth Newton is used, the solution can be obtained in a few iterations and without having to tune an extra parameter.Moreover we managed to prove in Theorems 4.3 and 4.4 optimal convergence rates, for piecewise affine and piecewise quadratic Lagrange finite elements, and for possibly nonconvex domains.Numerical experiments allowed to confirm in practice the predicted rates.
An important idea consists in reformulating the initial problem, with a constraint on the gradient of the solution, as an obstacle problem that involves the distance function to the boundary.Conversely to the most standard cases for the obstacle problem, this distance function is nonsmooth.However, some recent studies for the obstacle can circumvent this issue, as, for instance, Cicuttin et al. [12] for a discrete variational inequality obtained with the Hybrid High Order method.Some perspectives of this work may concern Hybrid High Order methods or adaptive finite element techniques.Also, though it has no clear physical motivation, it might be interesting, at least from the mathematical viewpoint, to study the case with a nonconstant source term.The error between  ℎ and a reference solution as a function of the mesh parameter ℎ for the nonconvex example using the method of [8].Check with your library that it subscribes to the journal, or consider making a personal donation to the S2O programme by contacting subscribers@edpsciences.org.
More information, including a list of supporters and financial transparency reports, is available at https://edpsciences.org/en/subscribe-to-open-s2o.

Figure 1 .
Figure 1.The first three meshes for the convex (top) and nonconvex (bottom) examples from the uniform mesh sequences.The side length of the larger square is 1 while the smaller square in the nonconvex example has a side length of 0.2.

Figure 6 .
Figure 6.The first three discrete solutions  ℎ (top) and gradient magnitudes |∇ ℎ | (bottom) for the convex example using quadratic elements.

Figure 8 .
Figure 8.The first three discrete solutions  ℎ (top) and gradient magnitudes |∇ ℎ | (bottom) for the nonconvex example using the method of [8].

Figure 9 .
Figure 9.The error between  ℎ and a reference solution as a function of the mesh parameter ℎ for the nonconvex example using the method of[8].

[
28] K. Mouallif, Approximation du probléme de la torsionélasto-plastique d'une barre cylindrique par régularisation et discrétisation d'un probléme inf-sup sur  1 0 (Ω) ×  ∞ + (Ω).Travaux Sém.Anal.Convexe 12 (1982) 24.Please help to maintain this journal in open access!This journal is currently published in open access under the Subscribe to Open model (S2O).We are thankful to our subscribers and supporters for making it possible to publish this journal in open access in the current year, free of charge for authors and readers.