A reduced basis method for frictional contact problems formulated with Nitsche’s method

We develop an eﬃcient reduced basis method for the frictional contact problem formulated using Nitsche’s method. We focus on the regime of small deformations and on Tresca friction. The key idea ensuring the computational eﬃciency of the method is to treat the nonlinearity resulting from the contact and friction conditions by means of the Empirical Interpolation Method. The proposed algorithm is applied to the Hertz contact problem between two half-disks with parameter-dependent radius. We also highlight the beneﬁts of the present approach with respect to the mixed (primal-dual) formulation.


Introduction
The reduced basis method (RBM) is a model reduction technique [29,5,30,19].The goal is to reduce the complexity of a parametrized model problem in computational studies where the parameters vary.The idea is to replace the high fidelity (HF) discretization space, which is supposed to be of very large dimension, by a small-dimensional subspace (called reduced space) which can be constructed by sampling the HF model.This allows one to organize the calculations in two phases.The first phase, called offline, is the construction phase of the reduced model.For this purpose, one considers a sample of the parameter space (assumed to be sufficiently representative) for which expensive calculations are performed by solving the HF problem for each parameter of the sample.The results of these calculations are then used to construct a small-dimensional subspace of the HF space, and the reduced model is built by replacing the HF space by the reduced subspace in a Galerkin-type approximation of the model.The second phase, called online, is a phase in which a large number of new values of the parameter are considered, for which accurate approximations of the HF solution are calculated by using the reduced model.The online phase is where substantial computational gains are achieved.
In this work, we are interested in the application of the RBM to the contact problem formulated with Nitsche's method.We focus on the regime of small deformations and on Tresca friction.The problem of mechanical contact [20,34] with or without friction is present in many structural problems encountered in several industrial fields.The variational formulation of this problem leads to a variational inequality of the first or second kind depending on whether there is friction or not [15,13].There are different approaches to solve these variational inequalities.We can mention mixed (primal-dual) methods [16,21,1] where Lagrange multipliers are introduced to enforce the contact and friction conditions.In this case, the problem to be solved is a saddle-point problem where one seeks a primal unknown (the displacement) and a dual unknown (the contact forces).One of the difficulties with these methods is that they require the contact operator to satisfy an inf-sup condition.In the literature, there is already some work on model reduction for the frictionless contact problem in the framework of a mixed formulation.For example, [18] derives model reduction methods in the general framework of variational inequalities including the unilateral contact problem.In [2], a projection-based method is proposed to reduce the contact problem under small deformations.In [14], an application of the hyper-reduction technique is presented for the contact problem under small deformations.We also mention [4] where a new dual basis construction is proposed for the RBM applied to the unilateral contact problem under large deformations.Finally, the recent work [27] proposes a stable and efficient model reduction method for the unilateral contact problem.Let us also mention [35] where the authors use the Progressive Generalized Method in order to build reduced-order models for problems with multiple contacts, and [24] where the authors use a hyper-reduction approach based on a reduced integration domain for the dual reduced basis.
In contrast to the mixed formulation approach, there are other methods to approximate the mechanical contact problem which are purely primal, i.e., they do not require the introduction of additional unknowns.These methods have the advantage of leading to unconstrained minimization problems (thus easier to solve) but do not guarantee that the contact and friction conditions are strictly satisfied.One example are penalty methods [32].Here, we focus on another primal approach based on Nitsche's method [28].This method was originally introduced for the reformulation of Dirichlet boundary conditions and extended in [9] to the frictionless contact problem in the framework of the finite element method.The main characteristic of Nitsche's method is that it is consistent, in contrast to classical penalty methods.In the last few years, many contributions have been made to this approach.For example, in [6], the method is extended to the case of contact with Tresca's friction; in [11], symmetric and nonsymmetric variants are presented; in [26], an extension to Coulomb's friction and large deformations is discussed; in [7], a nonconforming high-order discretization is considered; in [10], existence results for the contact problem with Coulomb friction are given in the context of static and dynamic finite element formulations.A state of the art on recent advances on Nitsche's method can be found in [8].To the best of our knowledge, there is no previous work on model reduction for the contact problem formulated with Nitsche's method.
In this paper, we propose to fill this gap for the frictionless contact problem and the contact problem with Tresca friction.The main challenge is that the classical RBM leads to an inefficient reduced model owing to the nonlinearity of Nitsche's formulation (even with small deformations).To overcome this problem, we propose a combination of the RBM with the Empirical Interpolation Method (EIM) [3,25].The realization of this idea is by no means straightforward since one needs to consider at the same time the parameter value and the iteration counter in the nonlinear iterative solver.The second important point addressed in this work is the comparison of the present approach with the inf-sup stable mixed formulation in terms of accuracy and efficiency.The two key advantages offered by Nitsche's method are the handling of unconstrained minimization problems and a higher effectiveness of the RBM since it is well-known that the dual basis is particularly hard to compress, as was highlighted in particular in [22].
The rest of this paper is organized as follows.In Section 2, we briefly recall the unilateral contact problem with friction and its variational formulation under the assumption of small deformations.In Section 3, we derive the formulation of this problem using Nitsche's method in a form suitable to the RBM.In Section 4, we present our main result, namely the procedure for building the reduced model with Nitsche's method using the RBM and the EIM.In Section 5, we provide numerical results showcasing the efficiency and the robustness of the proposed procedure and comparing it to the mixed formulation.We consider as test case the Hertz contact problem between two half-disks with parameter-dependent radius.The extension to Coulomb friction is briefly discussed at the end of Section 5.

Model problems
Let D ⊂ R m , m ∈ N * := N \ {0}, be the parameter set.For all µ ∈ D, we consider an elastic body whose reference configuration is the parameter-dependent domain Ω(µ The body is clamped at the boundary Γ D (µ), free of traction at the boundary Γ N (µ), and Γ c (µ) denotes the potential contact boundary with a given rigid support.We denote by n(µ) the unit outward normal on Γ(µ) and by τ (µ 1) an orthonormal basis of the hyperplane orthogonal to n(µ) in R d .For simplicity, we just write n and τ whenever there is no ambiguity.The body in its reference configuration is located at some distance from a rigid support and we denote by g(µ) ∈ L 2 (Γ c (µ); R + ) the corresponding gap function.An external load (µ) : Ω(µ) → R d is applied to the body, and we assume to be in the case of small deformations.For a generic R d -valued displacement field v, the R d×d -valued linearized strain tensor ε(v) and the R d×d -valued stress tensor σ(v) are given by with C the elastic coefficient tensor.At the boundary, we decompose the displacement field, v, and the normal component of the stress tensor, σ(v)n, in normal and tangential components as follows: with The frictionless contact problem (also called Signorini problem) consists in finding the displacement field u(µ) : Ω(µ) → R d satisfying, for all µ ∈ D, −div(σ(u(µ))) = (µ), in Ω(µ), (3a) In the case of contact problems with friction, the condition (3e) on Γ c (µ) has to be replaced by a condition depending on the considered friction law [21,26].Here, we focus on Tresca friction which leads to the following conditions: where • denotes the Euclidean norm in R d−1 and s > 0 is a given threshold, taken to be constant for simplicity (units in Pa).The model problem consisting of equations (3a)-(3b)-(3c)-(3d) and ( 4) is called Tresca frictional contact problem.We introduce the finite-dimensional space V(µ) and the admissible set K(µ) such that We notice that K(µ) is a non-empty convex set.The space V(µ) is typically built as a finite element space associated with a fine mesh of Ω(µ).The bilinear form a(µ; •, •) : V(µ) × V(µ) → R associated with the equilibrium equation (3a) in Ω(µ) is defined as and the linear form f (µ; •) : V(µ) → R associated with the external load (µ) as The weak formulation of the Signorini contact problem (3) consists of solving the following variational inequality of the first kind: For all µ ∈ D, find u(µ) ∈ K(µ) such that For all µ ∈ D, Stampacchia's theorem [33] ensures that there is a unique solution to (8) which is also the unique solution to the following constrained minimization problem: Find u(µ) ∈ K(µ) such that where the energy functional J (µ; •) : V(µ) → R is defined as follows: In the case of Tresca friction, we need to consider the friction functional F such that This leads to the following variational inequality: For all µ ∈ D, find u(µ) ∈ K(µ) such that

Nitsche's method
The main idea in the original Nitsche method [28] is the enforcement of Dirichlet boundary conditions by means of a consistent penalty method.As shown in [9,11,8], it is possible to generalize this idea to frictional contact problems.The main advantage of Nitsche's method is that the problem to be solved is unconstrained.Hence, in contrast to the mixed formulation, one does not need any additional unknowns such as Lagrange multipliers.Moreover, a higher effectivity of the RBM is expected since it is well-known that the dual basis is particularly hard to compress.The price to be paid, though, is that the constraint is not exactly enforced.
To derive Nitsche's method, we are going to assume that all the considered functions are smooth enough so that the associated normal stress tensor can be defined pointwise at the boundary.Recall that, for all µ ∈ D, the HF finite-dimensional space V(µ) in (5a) results from a finite element discretization of the Hilbert space H 1 (Ω(µ); R d ).Hence, the above assumption is indeed met.

Frictionless case
Following [12], the starting observation is that the Signorini conditions (3d) can be equivalently reformulated as follows: where z − := min z, 0 denotes the negative part of a generic real number z and where γ > 0 is a user-defined parameter (taken to be constant for simplicity).In practice, the parameter γ should be chosen large enough (see Section 5 for further discussion).For all µ ∈ D, one introduces the energy functional J Nitsche (µ; recalling that the energy functional J is defined in (10).Nitsche's method consists in finding u(µ) ∈ V(µ) solution to the following unconstrained minimization problem: Find u(µ) ∈ V(µ) such that The first-order optimality condition associated with (15) reads with the bilinear form a n γ (µ; •, •) : V(µ) × V(µ) → R defined as and the operators P n γ,g (µ; •), P n γ,0 (µ; •) : V(µ) → L 2 (Γ c (µ)) defined as With this notation, we can rewrite the energy functional J Nitsche (µ; •) as The problem ( 16) is nonlinear.To solve it, we use an iterative method.Given u 0 (µ) ∈ V(µ) and a tolerance δ u ∈ R + , we look, for all k ≥ 0, for the solution u k+1 (µ) in the form u k+1 (µ) = u k (µ)+δu k (µ).Ideally, we seek δu k (µ) ∈ V(µ) such that However, since the problem ( 20) is nonlinear, it is expensive to solve it directly.Therefore, we approximate the solution δu k (µ) by linearizing the problem.In order to do this, we observe that P n γ,g (µ; u k (µ) + δu k (µ)) = P n γ,g (µ; u k (µ)) + P n γ,0 (µ; δu k (µ)) and approximate the term P n γ,g (µ; u k (µ) + δu k (µ)) − as follows: where H(•) denotes the Heaviside function.Using this approximation in (20) (for simplicity, we keep the same notation for the unknown δu k (µ)), we consider the following sequence of problems: For all k ≥ 0, find δu k (µ) ∈ V(µ) such that where We iterate on k until the following convergence criterion is reached: In what follows, we denote by k cv (µ) ∈ N * the number of iterations required to solve (22), for all µ ∈ D.
We denote the converged solution to the sequence of problems (22) as u cv (µ) ∈ V(µ).
To introduce the algebraic formulation, we assume that for all µ ∈ D, the HF (finite-dimensional) space V(µ) is such that Notice that the dimension N HF of V(µ) is parameter-independent.We illustrate in Section 4.1 how to accomplish this property.Furthermore, we adopt the following decompositions for the sequence of solutions to (22): The algebraic formulation of the sequence of problems ( 22) then reads as follows: For all k ≥ 0, find where for all n, m ∈ {1:N HF }, The convergence criterion for (28) is still (25) using the reconstructed functions (see (27)).We denote the converged solution to the sequence of problems (28) as U cv (µ) ∈ R N HF .Furthermore, let us denote by Θ n γ (µ, w) (resp.F (µ)) the algebraic representation of the operator θ n γ (µ; w; •) (resp.f (µ; •)) such that for all m ∈ {1:N HF }, With this notation, we have the following decomposition:

Friction case
For the Tresca frictional contact problem, we use the following reformulation of the friction conditions given in [8]: where, for a positive real number r, the notation defines the projection of x ∈ R d−1 onto the ball centered at the origin and of radius r.Let us introduce the operator We start from the energy functional J Nitsche F (µ; •) : V(µ) → R defined as Nitsche's method consists in solving the following unconstrained minimization problem: Find u(µ) ∈ V(µ) such that Let us consider the functional J : R d−1 → R defined as follows: One easily verifies that this functional is Gâteaux-differentiable and that its differential is given by DJ(x) = x − x s for all x ∈ R d−1 .Therefore, using the same approach as for the frictionless contact problem, we obtain the following Nitsche's formulation for the Tresca frictional contact problem: Find where the bilinear form a nτ γ (µ; The linearization of the problem (38) leads to the following sequence of problems: For all k ≥ 0, find with where and G s (•), the differential of • s , is given by with I d−1 the identity matrix of order d − 1.The convergence criterion for (40) is still (25).We denote the converged solution to the sequence of problems (40) as u cv (µ) ∈ V(µ).
The algebraic formulation of the sequence of problems (40) reads as follows: For all k ≥ 0, find where for all n, m ∈ {1:N HF }, The convergence criterion for (44) is still (25) using the reconstructed functions.We denote the converged solution to the sequence of problems (44) as U cv (µ) ∈ R N HF .Furthermore, let us denote by With this notation, we have the following decompositions:

Reduced-basis formulation
In this section, we derive the RBM.We describe the problem in detail for the frictionless contact problem, and briefly highlight the (simple) adaptations needed to account for friction.

Geometric mapping
We recall that the dimension N HF of the HF (finite-dimensional) space V(µ) is parameter-independent.This is important since, in order to compress the space generated by the snapshots, it is necessary that all the snapshots live in the same space.For this purpose, since the geometry is parameter-dependent, we use a parameter-independent reference domain Ω.We assume that for all µ ∈ D, there is a smooth invertible geometric mapping h(µ) : Ω → Ω(µ).We denote by Γ := ∂ Ω the boundary of Ω and we assume that it can be partitioned as Γ = Γ D ∪ Γ N ∪ Γ c in such a way that, for all µ ∈ D, Then, the mesh of Ω(µ) is generated by generating a mesh of Ω matching the partition of the boundary Γ and applying the mapping h(µ) to the mesh of the reference domain Ω.The finite element basis functions are generated from the reference basis functions by using a pullback.

Naive approach
The goal of this section is to present the naive reduced model resulting from the application of a plain RBM to the contact problem formulated with Nitsche's method and highlight the computational inefficiency of such a formulation.This problem will be circumvented in the next section eventually leading to a computationally effective RBM.The major difficulty comes from the nonlinearity of Nitsche's formulation.
To build the reduced basis (RB), the starting point is to compute (in the offline phase) a family {U cv (µ p )} p∈{1:P } ⊂ R N HF of HF solutions to the frictionless contact problem (28) by using a training subset D train := {µ p } p∈{1:P } ⊂ D of cardinality P ∈ N * .Using the Proper Orthogonal Decomposition (POD) [17,23] based on the canonical inner product of H 1 ( Ω; R d ) and the geometric mapping h(µ), one can construct an orthonormal family {Ξ n } n∈{1:N } ⊂ R N HF of N ∈ N * (N ≤ P ) vectors.Let us denote by V N the reduced space generated by the family {Ξ n } n∈{1:N } , i.e., In algebraic form, the RB formulation of the sequence of HF problems (28) reads as follows: For all k ≥ 0, find where with The convergence criterion for (50) is still (25) using the reconstructed functions (see Remark 4.1 for more details).At this stage, the RBM consists of the following two stages: -Offline stage 1.Select a training subset D train := {µ p } p∈{1:P } ⊂ D.
3. Compute the reduced space V N by using POD on snapshots.
-Online stage: For any µ ∈ D \ D train , It is crucial to derive a reduced problem that is independent of the high-fidelity dimension N HF in order to obtain an inexpensive online stage.This condition is not yet satisfied with the current formalism.The main issue is the manipulation of large-dimensional arrays in (50).We propose in Section 4.3 a procedure to overcome this issue in order to construct a computationally efficient RBM.
With this notation, solving the RB problem (50) in algebraic form leads to the following reconstructed solutions:

Computationally efficient approach
In this section, we describe the strategy to avoid the manipulation of large-dimensional arrays in the problem (50).The idea consists in introducing appropriate affine parametric decompositions of the parameter-dependent and "parameter/iteration"-dependent operators involved in the problem by using the Empirical Interpolation Method (EIM) [3,25].Specifically, our goal is to separate the dependence on µ and k from the dependence on the indices in the large-dimensional arrays ).This operation is performed during the offline stage.For this purpose, using the EIM, we obtain the following approximations: where the large-dimensional arrays A n γ,s , B n γ,s , Θ n γ,s and F s are now independent of the parameter µ and the iteration counter k, whereas the functions α a n s and α f s (resp., α b n s and α θ n s ) only depend on µ (resp., (µ, k)).We obtain the following approximation of the residual: To build the large-dimensional arrays A n γ,s , B n γ,s , Θ n γ,s and F s , we respectively use the training sets E a n train , E b n train , E θ n train and E f train defined as follows: Notice that a different training set (possibly richer) than D train can be used instead.We introduce the index subsets {(n and ) and E f (µ) defined in (54).Notice that we use the HF solution u k (µ) instead of the RB solution u N,k (µ) to compute the functions α b n s and α θ n s .In the online phase, for every new value of the parameter pair (µ, k) ∈ D × N, the functions α a n s , α b n s , α θ n s and α f s are approximated by functions α a n N,s , α b n N,s , α θ n N,s and α f N,s which solve the following linear systems: where the vector-valued functions Notice that here, in the online phase, we use the RB solutions u N,k (µ).The parameter-independent interpolation matrices By construction, these matrices are lower-triangular with unit diagonal.Consequently, these matrices are invertible and their inverse can be easily computed once and for all during the offline phase.Combining (51) with (54), we obtain the following approximate decompositions: which lead to an efficient offline/online decomposition since the parameter-independent arrays A n γ,N,s , B n γ,N,s , Θ n γ,N,s and F N,s are small-dimensional and can be computed once and for all during the offline phase.Finally, using the approximations from (61) in (50) (for simplicity, we keep the same notation for the unknown ∆U N,k (µ)), we consider the following sequence of problems: For all (µ, k) where E r n N (µ, k) is given by The convergence criterion for (62) is still (25) using the reconstructed functions.We denote the converged solution to the sequence of problems (62) as U N,cv (µ) ∈ R N and the associated reconstructed solution as u N,cv (µ) ∈ V(µ).
To summarize, our RB procedure is organized as follows: -Offline stage 1.Select a training subset D train := {µ p } p∈{1:P } ⊂ D.
3. Compute the reduced space V N by using POD on the snapshots.

Compute the high-dimensional arrays
A n γ,s , B n γ,s , Θ n γ,s and F s by using the EIM.

Compute the small-dimensional arrays A n
γ,N,s , B n γ,N,s , Θ n γ,N,s and F N,s by using (61).

Evaluate E a n
N (µ) and E f N (µ) by using (61a) and (61d).For the contact problem with friction, the RB formulation is obtained in exactly the same way.We simply replace the forms a n γ (µ;

Numerical results
We consider the Hertz contact problem between the two half-disks as represented in Figure 1.The upper half-disk occupies the deformable domain Ω 1 (µ) ⊂ R 2 of parametric radius and the lower half-disk the rigid domain Ω 2 ⊂ R 2 of fixed radius R 2 := 1m.The initial gap between the two half-disks is equal to g 0 > 0. We impose a displacement of −d on Γ top 1 (µ) of Ω 1 (µ) with d ≥ g 0 .The initial gap g 0 and the imposed displacement d are, respectively, set to g 0 := 0.001m and d := 0.09m.This latter value, which is less than 10% of the maximum value of R 1 (µ) allows us to remain within the validity of the small deformation assumption.Notice that since Ω 2 is rigid and fixed, we only mesh the domain Ω 1 (µ) and set Ω(µ) := Ω 1 (µ) to build the HF space.The material parameters are E := 15Pa for the Young modulus and ν := 0.35 for the Poisson coefficient.A first training set is typically chosen as D train := 0.7 + 0.0075i, 0 ≤ i ≤ 60 (m) (altogether P = 61 points), and the validation set D valid is generated by choosing 30 elements in D randomly with a uniform distribution.A second richer training set (with P = 201 points) will also be considered.
We consider the reference domain Ω 1 := Ω 1 (1) and introduce the geometric mapping h 1 (µ) : Ω 1 → Ω 1 (µ) defined as h 1 (µ)(x) := µx, for all x ∈ Ω 1 , with the origin located at the center of Ω 1 .We use a mesh composed of 1956 nodes with 633 nodes on the potential contact manifold Γ c 1 which is the part of the half circle Γ c 1 of angle θ ∈ [− 5π 8 , − 3π 8 ] with respect to the horizontal axis.For all µ ∈ D, we equip the space V(µ) with the norm • V(µ) defined as follows: where the characteristic length := 1 is the radius of Ω 1 and is introduced for dimensional consistency.
The HF and RB computations use the python library of the finite element software getfem [31].

Frictionless case
We first consider the frictionless Hertz contact problem.
Figure 2 displays the deformed configurations resulting from the HF displacement fields u(µ) for µ = 0.7m (left panel) and for µ = 1.3m (right panel).We can see that we use a symmetric mesh.This is important because it guarantees the symmetry of the HF snapshots.Indeed, if the snapshots are not symmetric, the resulting POD modes will not be either.Consequently, the reduced model looses this symmetry property, leading to reduced solutions of poorer quality.Moreover, we have discretized a complete half-disk instead of a quarter-disk to avoid some difficulties when enforcing the symmetry condition at the lowest point of Ω 1 (µ) in the case of the frictional contact problems (see Section 5.2).Column 3).We see that the normal stress is of good quality (with almost no spurious oscillations) and matches very well with its counterpart P n γ,g (µ; u(µ)) − resulting from the Alart-Curnier reformulation.We also see that the negativity condition on the gap u n (µ) − g(µ) is satisfied on the whole potential contact manifold Γ c (µ).To have a better look at Signorini's contact conditions, we display in Table 1 the relative error on the Alart-Curnier reformulation of Signorini's contact conditions defined as follows: where the discrete 2 (Γ c (µ))-norms are sampled at the mesh nodes located on Γ c (µ).For the three values of the parameter µ ∈ {0.7, 1, 1.3}(m), we consider three values of the mesh size h, namely h = 5mm (coarse), h = 2.5mm (medium) and h = 1.25mm (fine).We notice that the relative error e n AC (µ) is smaller than 1.5% for the three parameters values and for the three mesh sizes.Moreover, we see that the error decreases when h decreases, thereby indicating the convergence of the approximation.More precisely, we observe a convergence of order 1 for the approximation of σ nn (u(µ)) (with a slight order reduction for the larger value of µ).Thus, we can say that the Signorini contact conditions are globally satisfied with a good accuracy although they are not strictly enforced.Relative error e n AC (µ) for µ ∈ {0.7, 1.0, 1.3}(m) and the mesh sizes h ∈ {5, 2.5, 1.25}(mm).
Let us consider the relative POD projection error defined as follows: where Π V N denotes the orthogonal projection onto V N ⊂ R N HF and W(µ p ) the Gram matrix of the inner product associated with H 1 (Ω(µ p ); R d ).We consider two different training sets, the first with cardinality 61 and the second with cardinality 201. Figure 4 shows the relative projection error e POD (N ) produced by the POD algorithm as a function of the number of vectors composing the reduced basis for both training sets.In all cases, we notice that the projection error decreases sufficiently fast so that indeed the linear spaces generated by the snapshots can be approximated by small-dimensional subspaces.We also observe a fast decrease of the POD error for the first 15 modes before a slower decrease occurs at error levels between 10 −5 and 10 −6 (resp.10 train and E θ n train introduced in (56) are then of cardinality 897.We fix a tolerance δ EIM := 10 −6 to bound the errors resulting from (54).With this choice, we obtain S b n = 619 N HF × N HF and S θ n = 281 N HF .Notice that in the present test case, we do not need to perform an EIM decomposition on A n γ (µ) since this matrix is already linearly dependent on µ; moreover, F (µ) vanishes since we only use an imposed displacement for the load.Figure 5 shows the relative EIM interpolation where for a generic matrix (resp.vector) V ∞ (j) := max and with D * either equal to D train or to D valid .Let us first consider the training set D train .We observe that both errors decrease fast enough to allow accurate approximations.For the tangent matrix B n γ (µ, u k (µ)), we notice a quasi-uniform decrease of the relative error e b n EIM (S b n , D train ) with an acceleration at error levels between 10 −3 and 10 −6 .For the residual vector Θ n γ (µ, u k (µ)), we observe a fast decrease of the relative error e θ n EIM (S θ n , D train ) for ranks S θ n between 1 and 40 yielding errors between 1 down to 10 −2 , and then a significant drop of the error at about S θ n = 40 before a slower decrease occurs at error levels between 10 −3 and 10 −6 .Considering the validation set D valid , we additionally plot the relative EIM approximation errors e b n ,cv EIM (S b n ) and e θ n ,cv EIM (S θ n ) defined as These errors correspond to the relative EIM approximation error at convergence of the iterative algorithm, i.e., when k = k cv (µ).For the tangent matrix we notice a quite modest decrease of the error for ranks S b n between 1 and 100 with errors values between 0.95 and 0.7, and then a stagnation of the error at around 0.65.Considering e b n ,cv EIM (S b n ), we observe instead a rather uniform decrease of the error, with values quite close to those of e b n EIM (S b n , D train ), before a stagnation occurs at a value of about 10 −3 .For the residual vector Θ n γ (µ, u k (µ)), considering first e θ n EIM (S θ n , D valid ), we observe a stagnation of the error for ranks S θ n between 1 and 80 with error values between 1 and 0.9, and then a slower decrease at error levels between 0.8 and 5 • 10 −2 with some stagnation phases.Considering e θ n ,cv EIM (S θ n ), we instead observe a stagnation at about 0.5 for ranks S θ n between 2 and 80, and then a slower decrease with error values between 0.3 down to 10 −2 .We conclude that the EIM approximation is not very accurate for µ ∈ D valid and small values of k, whereas the accuracy becomes more satisfactory as k → k cv (µ).Therefore, we may expect some difficulties in achieving convergence in the iterative solvers applied to reduced problems, but if convergence is indeed achieved, the accuracy should be reasonable.We denote by e u N (µ) (resp.e nn N (µ)) the relative RB approximation error on the displacement (resp.normal stress) defined as and introduce the relative error measures e u N,max and e nn N,max defined as Figure 6 displays e u N,max (left panel) and e nn N,max (right panel) as a function of the number of vectors composing the reduced basis.We only consider RB dimensions N larger than 10.Indeed, for smaller values, the iterative algorithm does not converge for some values of the parameter µ.This can be explained by the poor quality of the EIM approximation of the tangent matrix B n γ,N (µ, k) for small values of N due to the fact that the RB solution at the first iterations of the iterative algorithm on the reduced model is quite far from the HF solution on which the training of the EIM is performed.We also observe some convergence difficulties for values of the parameter µ larger than 1.2m.For this reason, we consider a validation set D valid restricted to the interval [0.7, 1.18](m).In Figure 6, we superpose the plain RBM approximation error (without any EIM, thus computationally inefficient) and the RBM-EIM approximation error.We observe that similar errors are obtained for plain RBM and RBM-EIM.This confirms the good quality of the EIM approximations at convergence as claimed above.With both approaches, we notice a stagnation of the relative error at about 10 −4 from about 40 modes, in agreement with the stagnation observed on the POD projection error.

Comparison with the mixed formulation
For the comparison, we consider the primal-dual formulation employed in [27] with P 2 finite elements for displacement and P 1 finite elements for the Lagrange multiplier.With this choice of the discretization, we can compare on a fair basis the HF displacements obtained with the mixed formulation and with Nitsche's method.We display in Figure 7 the HF energy J (µ; u(µ)) (see (10)) for all µ ∈ D train .We notice that we obtain (in the eyeball norm) the same values for the two methods for all µ ∈ D train .Thus, although the constraints are not strictly imposed with Nitsche's method, we obtain a satisfactory accuracy for the quality of the solution in comparison with the mixed formulation.We can also see that the energy decreases with the parameter radius µ of the half-disk Ω(µ).The complementary condition is not reported since it is on the order of the machine precision (10 −14 N/m) as expected with the primal-dual formulation since the constraints are exactly enforced.
In contrast to Nitsche's method, in the mixed formulation, it is necessary to stabilize the RBM in order to ensure inf-sup stability for the pair of primal/dual reduced spaces.For this purpose, we use the Projected Greedy Algorithm (PGA) algorithm from [27] with a tolerance δ PGA := inf µ∈Dtrain β HF (µ) c HF (µ) = 0.047 so that the stability condition established in [27,Prop 3.1] is fulfilled (the quantities β HF (µ) and c HF (µ) are defined therein).We denote by e u N,R (µ) (resp.e λ N,R (µ)) the primal (resp.dual) relative RB approximation error defined as follows: and introduce the relative primal (resp.dual) error measure e u N,R,max (resp.e λ N,R,max ) defined as  stagnation for all errors after a certain number of primal modes N .For the primal error, there is a clear decrease of the error as a function of the dimension of the dual basis R.However, for the dual error, although the error decreases, we observe some oscillations for certain values of the dimension of the primal basis N (compare the errors for R = 30 and R = 40).This can be explained by the fact that the dual basis obtained with the mCPG algorithm is not orthonormal (owing to the sign constraints on the Lagrange multiplier, since an orthonormalization process cannot be performed).Moreover, the higher the number of vectors in the dual basis, the more the orthogonality property is lost, and this fact introduces noise in the reduced model.In terms of accuracy, we observe that we have a better approximation for the primal variable (error of the order of 10 −4 m) than for the dual variable (error of the order of 10 −2 Pa) in the mixed formulation.These results illustrate the fact that it is very difficult to reduce the dual space, and therefore further motivate the use of purely primal methods like Nitsche's method in the RBM framework applied to contact problems.Finally, comparing the mixed and Nitsche formulations, we observe that the displacement error for the latter is better than for the former for the three values of the dimension of the dual reduced cone.Instead, comparing the error on the normal stress, we observe that the error for the former is slightly better when R = 40, but the dimension of the primal reduced space is much larger owing to the need to ensure inf-sup stability.

Friction case
Let us now consider the Tresca frictional Hertz contact problem.We choose a threshold s = 0.1.The parameter γ and the mesh size h are the same as for the frictionless case (see Section 5.1.1). Figure 9 displays the superposition of the tangential stress σ nτ (u(µ)) and its Alart-Curnier reformulation [P τ γ (µ; u(µ))] s as a function of the abscissa along Γ c (µ) for µ = 0.7m (Column 1), µ = 1m (Column 2), and µ = 1.3m (Column 3).We observe that the tangential stress matches very well with its counterpart [P τ γ (µ; u(µ))] s resulting from the Alart-Curnier reformulation.However, we observe some oscillations on σ nτ (u(µ)) at the end of the effective contact zone (contact/non-contact transition zone) and at the end of the potential contact zone.We also notice that the tangential stress is not zero outside the effective contact zone, which is consistent with a known (undesirable) feature of Tresca's model, namely to predict friction without contact.To have a better look at the friction  . (75) For the three values of the parameter µ ∈ {0.7, 1, 1.3}(m), we consider three values of the mesh size h, namely h = 5mm (coarse), h = 2.5mm (medium) and h = 1.25mm (fine).We notice that the relative error e nτ ,T AC (µ) is smaller than 6% for the three parameters values and for the three mesh sizes.Moreover, we see that the errors decrease when h decreases, thereby indicating the convergence of the approximation.More precisely, we observe a convergence of order 1  2 for the approximation of σ nτ (u(µ)).Notice that this order is different from the one observed in the frictionless case (order 1).Thus, we can say that the Tresca friction conditions are globally satisfied with a good accuracy as for the Signorini contact conditions although they are not strictly enforced.Relative errors e nτ ,T AC (µ) for µ ∈ {0.7, 1.0, 1.3}(m) and the mesh sizes h ∈ {5, 2.5, 1.25}(mm).
Figure 10 shows the relative projection error e POD (N ) produced by the POD algorithm (see (67)) as a function of the number of vectors composing the reduced basis.We notice that the projection error decreases sufficiently fast so that indeed the linear spaces generated by the snapshots can be approximated by small-dimensional subspaces.We also observe a fast decrease of the POD error for the first 15 modes before a slower decrease occurs at relative error levels between 10 −5 and 10 −6 .
Let us discuss the EIM approximation for which the training sets E b nτ train , E θ n train and E θ τ train are of cardinality 898.We fix a tolerance δ EIM := 10 −6 .With this choice, we obtain S b nτ = 630 N HF ×N HF , S θ n = 291 N HF and S θ τ = 3 N HF .For the same reason as for the frictionless case, we do not need to perform an EIM decomposition on A nτ γ (µ) and on F (µ).
and e θ n ,cv EIM (S θ n ) defined in (70b).Notice that for Tresca friction, we perform the EIM on the tangent matrix B nτ γ (µ; u k (µ)) instead of performing it separately on the tangent matrices B n γ (µ; u k (µ)) and B τ γ (µ; u k (µ)).We see that the result is close to the one obtained for the frictionless case because the contribution of B τ γ (µ; u k (µ)) is negligible.Indeed, as can be seen in Figure 9 (Row 1), the Alart-Curnier reformulation of the tangential stress is equal to ±s almost everywhere (except at 3 nodes) which results in the nullity of G s (P τ γ (µ; u k (µ))) almost everywhere.For the same reason, the dependence of the residual vector Θ τ γ (µ, u k (µ)) on (µ, k) is almost of rank one and therefore its EIM approximation is very simple (S θ τ = 3); hence, its relative EIM approximation error is not reported.Altogether, the behaviour of the residual vector Θ n γ (µ, u k (µ)) remains similar to that observed in the frictionless case.We denote by e nτ N (µ) the relative RB approximation error on the tangential stress defined as

Toward Coulomb friction
A more realistic model for friction is given by Coulomb conditions which read as follows: HF solutions for the Hertz test case obtained using the above nested loop can be obtained to generate snapshots for the parameter values from the training set.Three snapshots are illustrated in Figure 13 (compare with Figure 9).Moreover, the set of snapshots can be compressed by using POD, leading to similar results to those obtained for Tresca friction (see Figure 10).The main challenge within the current approach lies in the realization of the EIM since the nested iterative loop now requires to separate the dependencies on the triple (µ, n, k).To overcome this difficulty, one possibility is to consider only converged solutions for the inner iteration (index k) when computing the EIM decompositions of the tangent matrix and the residual.However, these decompositions turn out to be, so far, rather inaccurate at the early stages of the iterative procedure, thereby hampering convergence.This difficulty will be further investigated in future work.

Figure 3
Figure 3 displays the superposition of the normal stress σ nn (u(µ)) and its Alart-Curnier reformulation P n γ,g (µ; u(µ)) − (Row 1) and the gap on the deformed configuration u n (µ) − g(µ) (Row 2) as a function of the abscissa along Γ c (µ) for µ = 0.7m (Column 1), µ = 1m (Column 2), and µ = 1.3m (Column 3).We see that the normal stress is of good quality (with almost no spurious oscillations) and matches very well with its counterpart P n γ,g (µ; u(µ)) − resulting from the Alart-Curnier reformulation.We also see that the negativity condition on the gap u n (µ) − g(µ) is satisfied on the whole potential contact manifold Γ c (µ).To have a better look at Signorini's contact conditions, we display in Table1the relative error on the Alart-Curnier reformulation of Signorini's contact conditions defined as follows:

−5 and 10 − 7 )
for the first (resp.second) training set.For the EIM approximation, we use the training set D train of cardinality 61, and the training sets E b n

Figure 5 :
Figure 5: Frictionless Hertz test case, Nitsche's method: Relative EIM approximation errors as a function of the rank S b n or S θ n of the approximation.Left: e b n EIM (S b n , D train ), e b n EIM (S b n , D valid ) and e b n ,cv EIM (S b n ); Right: e θ n EIM (S θ n , D train ), e θ n EIM (S θ n , D valid ) and e θ n ,cv EIM (S θ n ).

Figure 7 :
Figure 7: Frictionless Hertz test case: Comparison of the HF energy J (µ; u(µ)) between the mixed formulation and Nitsche's method for all µ ∈ D train .

Figure 8
Figure 8 displays the quantity e u N,R,max (resp.e λ N,R,max ) on the left (resp.right) panel as a function of the dimension of the reduced primal basis (after stabilization) for three values of the dimension of the reduced dual basis R, namely R = 10, R = 30 and R = 40.In addition, we plot the relative error e u N,max (resp.e nn N,max ) in the left (resp.right) panel in order to compare the displacement (resp.normal stress) error between the mixed and Nitsche approaches.Considering the mixed approach, we notice a

Figure 8 :
Figure 8: Frictionless Hertz test case, mixed method: RBM approximation errors.Right: displacement errors e u N,R,max and e u N,max ; Left: normal stress errors e λ N,R,max and e nn N,max .

Figure 10 :
Figure 10: Frictional Hertz test case, Nitsche's method: Relative POD projection error e POD (N ) as function of N , for |D train | = 61.

Figure 11 :
Figure 11: Frictional Hertz test case, Nitsche's method: Relative EIM approximation errors as a function of the rank S b nτ or S θ n of the approximation.Left: e b nτ EIM (S b nτ , D train ), e b nτ EIM (S b nτ , D valid ) and e b nτ ,cv EIM (S b nτ ); Right: e θ n EIM (S θ n , D train ), e θ n EIM (S θ n , D valid ) and e θ n ,cv EIM (S θ n ).
and S f respectively, corresponding to the indices selected by the EIM for the approximation of A n n (µ)) n a n s m a n s