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

Idrissa Niakh, Guillaume Drouet, Virginie Ehrlacher, Alexandre Ern

## ► To cite this version:

Idrissa Niakh, Guillaume Drouet, Virginie Ehrlacher, Alexandre Ern. A reduced basis method for frictional contact problems formulated with Nitsche's method. 2024. hal-04168418v2

HAL Id: hal-04168418

https://inria.hal.science/hal-04168418v2

Preprint submitted on 16 Jan 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A reduced basis method for frictional contact problems formulated with Nitsche's method

IDRISSA NIAKH *†‡, GUILLAUME DROUET*, VIRGINIE EHRLACHER †‡, ALEXANDRE ERN †‡

January 16, 2024

## Abstract

We develop an efficient 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 efficiency 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 benefits of the present approach with respect to the mixed (primal-dual) formulation.

***Keywords*** — model reduction, variational inequalities, reduced basis method, contact problems, Nitsche's method, Tresca friction, Coulomb friction.

## 1 Introduction

The reduced basis method (RBM) is a model reduction technique [31, 5, 32, 21]. 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 [22, 36] 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

<footnote>*EDF R&D, 7 Boulevard Gaspard Monge, 91120 Palaiseau, France</footnote>

<footnote>†CERMICS, École des Ponts, 6-8 avenue Blaise Pascal, 77455 Marne-la-Vallée cedex 2, France</footnote>

<footnote>‡INRIA Paris, 2 Rue Simone Iff, 75012 Paris, France</footnote>
# A reduced basis method for frictional contact problems formulated with Nitsche's method

IDRISSA NIAKH *†‡, GUILLAUME DROUET, VIRGINIE EHRLACHER †‡, ALEXANDRE ERN †‡

January 16, 2024

## Abstract

We develop an efficient 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 efficiency 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 benefits of the present approach with respect to the mixed (primal-dual) formulation.

**Keywords** — model reduction, variational inequalities, reduced basis method, contact problems, Nitsche's method, Tresca friction, Coulomb friction.

## 1 Introduction

The reduced basis method (RBM) is a model reduction technique [31, 5, 32, 21]. 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 [22, 36] 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

<footnote>*EDF R&D, 7 Boulevard Gaspard Monge, 91120 Palaiseau, France</footnote>

<footnote>†CERMICS, École des Ponts, 6-8 avenue Blaise Pascal, 77455 Marne-la-Vallée cedex 2, France</footnote>

<footnote>‡INRIA Paris, 2 Rue Simone Iff, 75012 Paris, France</footnote>
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, 23, 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, [20] 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 [29] proposes a stable and efficient model reduction method for the unilateral contact problem. Let us also mention [37] where the authors use the Progressive Generalized Method in order to build reduced-order models for problems with multiple contacts, and [26] 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 [34]. Here, we focus on another primal approach
based on Nitsche's method [30]. 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 [28], 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 review on the state of the art for Nitsche's method
can be found in [8], while the two more recent contributions [18, 17] leverage on the viewpoint of
stabilized finite elements to lead the analysis, including a posteriori error estimates. 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, 27]. 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 [24].

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
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, 23, 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, [20] 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 [29] proposes a stable and efficient model reduction method for the unilateral contact problem. Let us also mention [37] where the authors use the Progressive Generalized Method in order to build reduced-order models for problems with multiple contacts, and [26] 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 [34]. Here, we focus on another primal approach
based on Nitsche's method [30]. 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 [28], 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 review on the state of the art for Nitsche's method
can be found in [8], while the two more recent contributions [18, 17] leverage on the viewpoint of
stabilized finite elements to lead the analysis, including a posteriori error estimates. 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, 27]. 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 [24].

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.

## 2 Model problems

Let $\mathcal{D} \subset \mathbb{R}^m, m \in \mathbb{N}^* := \mathbb{N} \setminus \{0\}$, be the parameter set. For all $\mu \in \mathcal{D}$, we consider an elastic body whose reference configuration is the parameter-dependent domain $\Omega(\mu) \subset \mathbb{R}^d$, with $d \in \{2, 3\}$. The boundary $\Gamma(\mu) := \partial\Omega(\mu)$ is partitioned as $\Gamma(\mu) = \Gamma^D(\mu) \cup \Gamma^N(\mu) \cup \Gamma^c(\mu)$. The body is clamped at the boundary $\Gamma^D(\mu)$, free of traction at the boundary $\Gamma^N(\mu)$, and $\Gamma^c(\mu)$ denotes the potential contact boundary with a given rigid support. We denote by $\mathbf{n}(\mu)$ the unit outward normal on $\Gamma(\mu)$ and by $\boldsymbol{\tau}(\mu) := [\boldsymbol{\tau}_1(\mu) \cdots \boldsymbol{\tau}_{d-1}(\mu)] \in \mathbb{R}^{d \times (d-1)}$ an orthonormal basis of the hyperplane orthogonal to $\mathbf{n}(\mu)$ in $\mathbb{R}^d$. For simplicity, we just write $\mathbf{n}$ and $\boldsymbol{\tau}$ 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(\mu) \in L^2(\Gamma^c(\mu); \mathbb{R}_+)$ the corresponding gap function. An external load $\ell(\mu) : \Omega(\mu) \to \mathbb{R}^d$ is applied to the body, and we assume to be in the case of small deformations. For a generic $\mathbb{R}^d$-valued displacement field $\mathbf{v}$, the $\mathbb{R}^{d \times d}$-valued linearized strain tensor $\varepsilon(\mathbf{v})$ and the $\mathbb{R}^{d \times d}$-valued stress tensor $\sigma(\mathbf{v})$ are given by

$$ \varepsilon(\mathbf{v}) := \frac{1}{2}(\nabla \mathbf{v} + \nabla \mathbf{v}^{\top}), \quad \sigma(\mathbf{v}) := \mathcal{C}\varepsilon(\mathbf{v}), \qquad (1) $$

with $\mathcal{C}$ the elastic coefficient tensor. At the boundary, we decompose the displacement field, $\mathbf{v}$, and the normal component of the stress tensor, $\sigma(\mathbf{v})\mathbf{n}$, in normal and tangential components as follows:

$$ \mathbf{v} = v_n \mathbf{n} + \boldsymbol{\tau} \mathbf{v}_{\boldsymbol{\tau}}, \quad \sigma(\mathbf{v}) \mathbf{n} = \sigma_{nn}(\mathbf{v}) \mathbf{n} + \boldsymbol{\tau} \sigma_{n\boldsymbol{\tau}}(\mathbf{v}), \qquad (2) $$

with $v_n \in \mathbb{R}$, $v_{\boldsymbol{\tau}} \in \mathbb{R}^{d-1}$, $\sigma_{nn}(\mathbf{v}) \in \mathbb{R}$ and $\sigma_{n\boldsymbol{\tau}}(\mathbf{v}) \in \mathbb{R}^{d-1}$.

The frictionless contact problem (also called Signorini problem) consists in finding the displacement field $\mathbf{u}(\mu) : \Omega(\mu) \to \mathbb{R}^d$ satisfying, for all $\mu \in \mathcal{D}$,

$$ -\mathrm{div}(\sigma(\mathbf{u}(\mu))) = \ell(\mu), \quad \text{in } \Omega(\mu), \qquad (3a) $$

$$ \mathbf{u}(\mu) = \mathbf{0}, \quad \text{on } \Gamma^D(\mu), \qquad (3b) $$

$$ \sigma(\mathbf{u}(\mu))\mathbf{n} = \mathbf{0}, \quad \text{on } \Gamma^N(\mu), \qquad (3c) $$

$$ u_n(\mu) \le g(\mu), \quad \sigma_{nn}(\mathbf{u}(\mu)) \le 0, \quad \sigma_{nn}(\mathbf{u}(\mu))(u_n(\mu) - g(\mu)) = 0, \quad \text{on } \Gamma^c(\mu), \qquad (3d) $$

$$ \sigma_{n\boldsymbol{\tau}}(\mathbf{u}(\mu)) = \mathbf{0}, \quad \text{on } \Gamma^c(\mu). \qquad (3e) $$

In the case of contact problems with friction, the condition (3e) on $\Gamma^c(\mu)$ has to be replaced by a condition depending on the considered friction law [23, 28]. Here, we focus on Tresca friction which leads to the following conditions:

$$ \left\{ \begin{array}{ll} \|\sigma_{n\boldsymbol{\tau}}(\mathbf{u}(\mu))\| \le s, & \text{if } \mathbf{u}_{\boldsymbol{\tau}}(\mu) = \mathbf{0}, \\ \sigma_{n\boldsymbol{\tau}}(\mathbf{u}(\mu)) = -s \frac{\mathbf{u}_{\boldsymbol{\tau}}(\mu)}{\|\mathbf{u}_{\boldsymbol{\tau}}(\mu)\|}, & \text{otherwise,} \end{array} \right. \qquad (4) $$

where $\| \cdot \|$ denotes the Euclidean norm in $\mathbb{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.
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.

## 2 Model problems

Let $\mathcal{D} \subset \mathbb{R}^m, m \in \mathbb{N}^* := \mathbb{N} \setminus \{0\}$, be the parameter set. For all $\mu \in \mathcal{D}$, we consider an elastic body whose reference configuration is the parameter-dependent domain $\Omega(\mu) \subset \mathbb{R}^d$, with $d \in \{2, 3\}$. The boundary $\Gamma(\mu) := \partial\Omega(\mu)$ is partitioned as $\Gamma(\mu) = \Gamma^D(\mu) \cup \Gamma^N(\mu) \cup \Gamma^c(\mu)$. The body is clamped at the boundary $\Gamma^D(\mu)$, free of traction at the boundary $\Gamma^N(\mu)$, and $\Gamma^c(\mu)$ denotes the potential contact boundary with a given rigid support. We denote by $\mathbf{n}(\mu)$ the unit outward normal on $\Gamma(\mu)$ and by $\boldsymbol{\tau}(\mu) := [\boldsymbol{\tau}_1(\mu) \cdots \boldsymbol{\tau}_{d-1}(\mu)] \in \mathbb{R}^{d \times (d-1)}$ an orthonormal basis of the hyperplane orthogonal to $\mathbf{n}(\mu)$ in $\mathbb{R}^d$. For simplicity, we just write $\mathbf{n}$ and $\boldsymbol{\tau}$ 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(\mu) \in L^2(\Gamma^c(\mu); \mathbb{R}_+)$ the corresponding gap function. An external load $\ell(\mu) : \Omega(\mu) \to \mathbb{R}^d$ is applied to the body, and we assume to be in the case of small deformations. For a generic $\mathbb{R}^d$-valued displacement field $\mathbf{v}$, the $\mathbb{R}^{d \times d}$-valued linearized strain tensor $\varepsilon(\mathbf{v})$ and the $\mathbb{R}^{d \times d}$-valued stress tensor $\sigma(\mathbf{v})$ are given by

$$ \varepsilon(\mathbf{v}) := \frac{1}{2}(\nabla \mathbf{v} + \nabla \mathbf{v}^{\top}), \quad \sigma(\mathbf{v}) := \mathcal{C}\varepsilon(\mathbf{v}), \qquad (1) $$

with $\mathcal{C}$ the elastic coefficient tensor. At the boundary, we decompose the displacement field, $\mathbf{v}$, and the normal component of the stress tensor, $\sigma(\mathbf{v})\mathbf{n}$, in normal and tangential components as follows:

$$ \mathbf{v} = v_n \mathbf{n} + \boldsymbol{\tau} \mathbf{v}_{\boldsymbol{\tau}}, \quad \sigma(\mathbf{v}) \mathbf{n} = \sigma_{nn}(\mathbf{v}) \mathbf{n} + \boldsymbol{\tau} \sigma_{n\boldsymbol{\tau}}(\mathbf{v}), \qquad (2) $$

with $v_n \in \mathbb{R}$, $v_{\boldsymbol{\tau}} \in \mathbb{R}^{d-1}$, $\sigma_{nn}(\mathbf{v}) \in \mathbb{R}$ and $\sigma_{n\boldsymbol{\tau}}(\mathbf{v}) \in \mathbb{R}^{d-1}$.

The frictionless contact problem (also called Signorini problem) consists in finding the displacement field $\mathbf{u}(\mu) : \Omega(\mu) \to \mathbb{R}^d$ satisfying, for all $\mu \in \mathcal{D}$,

$$ -\mathrm{div}(\sigma(\mathbf{u}(\mu))) = \ell(\mu), \quad \text{in } \Omega(\mu), \qquad (3a) $$

$$ \mathbf{u}(\mu) = \mathbf{0}, \quad \text{on } \Gamma^D(\mu), \qquad (3b) $$

$$ \sigma(\mathbf{u}(\mu))\mathbf{n} = \mathbf{0}, \quad \text{on } \Gamma^N(\mu), \qquad (3c) $$

$$ u_n(\mu) \le g(\mu), \quad \sigma_{nn}(\mathbf{u}(\mu)) \le 0, \quad \sigma_{nn}(\mathbf{u}(\mu))(u_n(\mu) - g(\mu)) = 0, \quad \text{on } \Gamma^c(\mu), \qquad (3d) $$

$$ \sigma_{n\boldsymbol{\tau}}(\mathbf{u}(\mu)) = \mathbf{0}, \quad \text{on } \Gamma^c(\mu). \qquad (3e) $$

In the case of contact problems with friction, the condition (3e) on $\Gamma^c(\mu)$ has to be replaced by a condition depending on the considered friction law [23, 28]. Here, we focus on Tresca friction which leads to the following conditions:

$$ \left\{ \begin{array}{ll} \|\sigma_{n\boldsymbol{\tau}}(\mathbf{u}(\mu))\| \le s, & \text{if } \mathbf{u}_{\boldsymbol{\tau}}(\mu) = \mathbf{0}, \\ \sigma_{n\boldsymbol{\tau}}(\mathbf{u}(\mu)) = -s \frac{\mathbf{u}_{\boldsymbol{\tau}}(\mu)}{\|\mathbf{u}_{\boldsymbol{\tau}}(\mu)\|}, & \text{otherwise,} \end{array} \right. \qquad (4) $$

where $\| \cdot \|$ denotes the Euclidean norm in $\mathbb{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 $\mathcal{V}(\mu)$ and the admissible set $\mathcal{K}(\mu)$ such that

$$ \mathcal{V}(\mu) \subset \{v \in H^1(\Omega(\mu); \mathbb{R}^d) \mid v = 0 \text{ on } \Gamma^\mathrm{D}(\mu)\}, \quad (5a) $$

$$ \mathcal{K}(\mu) := \{v \in \mathcal{V}(\mu) \mid v_n \le g(\mu) \text{ on } \Gamma^c(\mu)\}. \quad (5b) $$

We notice that $\mathcal{K}(\mu)$ is a non-empty convex set. The space $\mathcal{V}(\mu)$ is typically built as a finite element space associated with a fine mesh of $\Omega(\mu)$. The bilinear form $a(\mu; \cdot, \cdot) : \mathcal{V}(\mu) \times \mathcal{V}(\mu) \to \mathbb{R}$ associated with the equilibrium equation (3a) in $\Omega(\mu)$ is defined as

$$ a(\mu; \mathbf{u}, \mathbf{v}) := \int_{\Omega(\mu)} \sigma(\mathbf{u}) : \varepsilon(\mathbf{v}) \, d\Omega(\mu), \quad (6) $$

and the linear form $f(\mu; \cdot) : \mathcal{V}(\mu) \to \mathbb{R}$ associated with the external load $\ell(\mu)$ as

$$ f(\mu; \mathbf{v}) := \int_{\Omega(\mu)} \ell(\mu) \cdot \mathbf{v} \, d\Omega(\mu). \quad (7) $$

The weak formulation of the Signorini contact problem (3) consists of solving the following variational inequality of the first kind: For all $\mu \in \mathcal{D}$, find $\mathbf{u}(\mu) \in \mathcal{K}(\mu)$ such that

$$ a(\mu; \mathbf{u}(\mu), \mathbf{v} - \mathbf{u}(\mu)) \ge f(\mu; \mathbf{v} - \mathbf{u}(\mu)), \quad \forall \mathbf{v} \in \mathcal{K}(\mu). \quad (8) $$

For all $\mu \in \mathcal{D}$, Stampacchia's theorem [35] ensures that there is a unique solution to (8) which is also the unique solution to the following constrained minimization problem: Find $\mathbf{u}(\mu) \in \mathcal{K}(\mu)$ such that

$$ \mathbf{u}(\mu) = \underset{\mathbf{v} \in \mathcal{K}(\mu)}{\operatorname{argmin}} \mathcal{J}(\mu; \mathbf{v}), \quad (9) $$

where the energy functional $\mathcal{J}(\mu; \cdot) : \mathcal{V}(\mu) \to \mathbb{R}$ is defined as follows:

$$ \mathcal{J}(\mu; \mathbf{v}) := \frac{1}{2} a(\mu; \mathbf{v}, \mathbf{v}) - f(\mu; \mathbf{v}). \quad (10) $$

In the case of Tresca friction, we need to consider the friction functional $\mathcal{F}$ such that

$$ \mathcal{F}(\mu; \mathbf{v}) := \int_{\Gamma^c(\mu)} -s \|\mathbf{v}_{\tau}\| \, d\Gamma(\mu). \quad (11) $$

This leads to the following variational inequality: For all $\mu \in \mathcal{D}$, find $\mathbf{u}(\mu) \in \mathcal{K}(\mu)$ such that

$$ a(\mu; \mathbf{u}(\mu), \mathbf{v} - \mathbf{u}(\mu)) + \mathcal{F}(\mu; \mathbf{v}) - \mathcal{F}(\mu; \mathbf{u}(\mu)) \ge f(\mu; \mathbf{v} - \mathbf{u}(\mu)), \quad \forall \mathbf{v} \in \mathcal{K}(\mu). \quad (12) $$

# 3 Nitsche's method

The main idea in the original Nitsche method [30] 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 $\mu \in \mathcal{D}$, the HF finite-dimensional space $\mathcal{V}(\mu)$ in (5a) results from a finite element discretization of the Hilbert space $H^1(\Omega(\mu); \mathbb{R}^d)$. Hence, the above assumption is indeed met.
We introduce the finite-dimensional space $\mathcal{V}(\mu)$ and the admissible set $\mathcal{K}(\mu)$ such that

$$ \mathcal{V}(\mu) \subset \{\mathbf{v} \in H^1(\Omega(\mu); \mathbb{R}^d) \mid \mathbf{v} = \mathbf{0} \text{ on } \Gamma^\mathrm{D}(\mu)\}, \qquad (5a) $$

$$ \mathcal{K}(\mu) := \{\mathbf{v} \in \mathcal{V}(\mu) \mid v_n \le g(\mu) \text{ on } \Gamma^c(\mu)\}. \qquad (5b) $$

We notice that $\mathcal{K}(\mu)$ is a non-empty convex set. The space $\mathcal{V}(\mu)$ is typically built as a finite element space associated with a fine mesh of $\Omega(\mu)$. The bilinear form $a(\mu; \cdot, \cdot) : \mathcal{V}(\mu) \times \mathcal{V}(\mu) \to \mathbb{R}$ associated with the equilibrium equation (3a) in $\Omega(\mu)$ is defined as

$$ a(\mu; \mathbf{u}, \mathbf{v}) := \int_{\Omega(\mu)} \sigma(\mathbf{u}) : \varepsilon(\mathbf{v}) \, d\Omega(\mu), \qquad (6) $$

and the linear form $f(\mu; \cdot) : \mathcal{V}(\mu) \to \mathbb{R}$ associated with the external load $\ell(\mu)$ as

$$ f(\mu; \mathbf{v}) := \int_{\Omega(\mu)} \ell(\mu) \cdot \mathbf{v} \, d\Omega(\mu). \qquad (7) $$

The weak formulation of the Signorini contact problem (3) consists of solving the following variational inequality of the first kind: For all $\mu \in \mathcal{D}$, find $\mathbf{u}(\mu) \in \mathcal{K}(\mu)$ such that

$$ a(\mu; \mathbf{u}(\mu), \mathbf{v} - \mathbf{u}(\mu)) \ge f(\mu; \mathbf{v} - \mathbf{u}(\mu)), \quad \forall \mathbf{v} \in \mathcal{K}(\mu). \qquad (8) $$

For all $\mu \in \mathcal{D}$, Stampacchia's theorem [35] ensures that there is a unique solution to (8) which is also the unique solution to the following constrained minimization problem: Find $\mathbf{u}(\mu) \in \mathcal{K}(\mu)$ such that

$$ \mathbf{u}(\mu) = \underset{\mathbf{v} \in \mathcal{K}(\mu)}{\operatorname{argmin}} \mathcal{J}(\mu; \mathbf{v}), \qquad (9) $$

where the energy functional $\mathcal{J}(\mu; \cdot) : \mathcal{V}(\mu) \to \mathbb{R}$ is defined as follows:

$$ \mathcal{J}(\mu; \mathbf{v}) := \frac{1}{2} a(\mu; \mathbf{v}, \mathbf{v}) - f(\mu; \mathbf{v}). \qquad (10) $$

In the case of Tresca friction, we need to consider the friction functional $\mathcal{F}$ such that

$$ \mathcal{F}(\mu; \mathbf{v}) := \int_{\Gamma^c(\mu)} -s \|\mathbf{v}_{\tau}\| \, d\Gamma(\mu). \qquad (11) $$

This leads to the following variational inequality: For all $\mu \in \mathcal{D}$, find $\mathbf{u}(\mu) \in \mathcal{K}(\mu)$ such that

$$ a(\mu; \mathbf{u}(\mu), \mathbf{v} - \mathbf{u}(\mu)) + \mathcal{F}(\mu; \mathbf{v}) - \mathcal{F}(\mu; \mathbf{u}(\mu)) \ge f(\mu; \mathbf{v} - \mathbf{u}(\mu)), \quad \forall \mathbf{v} \in \mathcal{K}(\mu). \qquad (12) $$

# 3 Nitsche's method

The main idea in the original Nitsche method [30] 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 $\mu \in \mathcal{D}$, the HF finite-dimensional space $\mathcal{V}(\mu)$ in (5a) results from a finite element discretization of the Hilbert space $H^1(\Omega(\mu); \mathbb{R}^d)$. Hence, the above assumption is indeed met.
## 3.1 Frictionless case

Following [12], the starting observation is that the Signorini conditions (3d) can be equivalently reformulated as follows:

$$ \sigma_{nn}(\mathbf{u}(\mu)) = [\sigma_{nn}(\mathbf{u}(\mu)) - \gamma(u_n(\mu) - g(\mu))]_{-}, \quad (13) $$

where $[z]_{-} := \min(z, 0)$ denotes the negative part of a generic real number $z$ and where $\gamma > 0$ is a user-defined parameter (taken to be constant for simplicity). In practice, the parameter $\gamma$ should be chosen large enough (see Section 5 for further discussion).

For all $\mu \in \mathcal{D}$, one introduces the energy functional $J^{\text{Nitsche}}(\mu; \cdot) : \mathcal{V}(\mu) \to \mathbb{R}$ such that

$$ J^{\text{Nitsche}}(\mu; \mathbf{v}) := \mathcal{J}(\mu; \mathbf{v}) - \frac{1}{2} \int_{\Gamma^c(\mu)} \frac{1}{\gamma} |\sigma_{nn}(\mathbf{v})|^2 d\Gamma(\mu) \\ + \frac{1}{2} \int_{\Gamma^c(\mu)} \frac{1}{\gamma} [\sigma_{nn}(\mathbf{v}) - \gamma(v_n - g(\mu))]^2 d\Gamma(\mu), \quad (14) $$

recalling that the energy functional $\mathcal{J}$ is defined in (10). Nitsche's method consists in finding $\mathbf{u}(\mu) \in \mathcal{V}(\mu)$ solution to the following unconstrained minimization problem: Find $\mathbf{u}(\mu) \in \mathcal{V}(\mu)$ such that

$$ \mathbf{u}(\mu) = \underset{\mathbf{v} \in \mathcal{V}(\mu)}{\operatorname{argmin}} J^{\text{Nitsche}}(\mu; \mathbf{v}). \quad (15) $$

The first-order optimality condition associated with (15) reads

$$ a_{\gamma}^{n}(\mu; \mathbf{u}(\mu), \mathbf{v}) + \int_{\Gamma^{c}(\mu)} \frac{1}{\gamma} [P_{\gamma,g}^{n}(\mu; \mathbf{u}(\mu))]_{-} P_{\gamma,0}^{n}(\mu; \mathbf{v}) d\Gamma(\mu) = f(\mu; \mathbf{v}), \quad \forall \mathbf{v} \in \mathcal{V}(\mu), \quad (16) $$

with the bilinear form $a_{\gamma}^{n}(\mu; \cdot, \cdot) : \mathcal{V}(\mu) \times \mathcal{V}(\mu) \to \mathbb{R}$ defined as

$$ a_{\gamma}^{n}(\mu; \mathbf{u}, \mathbf{v}) := a(\mu; \mathbf{u}, \mathbf{v}) - \int_{\Gamma^{c}(\mu)} \frac{1}{\gamma} \sigma_{nn}(\mathbf{u}) \sigma_{nn}(\mathbf{v}) d\Gamma(\mu), \quad \forall \mathbf{u}, \mathbf{v} \in \mathcal{V}(\mu), \quad (17) $$

and the operators $P_{\gamma,g}^{n}(\mu; \cdot, \cdot)$, $P_{\gamma,0}^{n}(\mu; \cdot, \cdot) : \mathcal{V}(\mu) \to L^2(\Gamma^c(\mu))$ defined as

$$ P_{\gamma,g}^{n}(\mu; \mathbf{v}) := \sigma_{nn}(\mathbf{v}) - \gamma(v_n - g(\mu)), \quad (18a) $$

$$ P_{\gamma,0}^{n}(\mu; \mathbf{v}) := \sigma_{nn}(\mathbf{v}) - \gamma v_n. \quad (18b) $$

With this notation, we can rewrite the energy functional $J^{\text{Nitsche}}(\mu; \cdot)$ as

$$ J^{\text{Nitsche}}(\mu; \mathbf{v}) := \mathcal{J}(\mu; \mathbf{v}) - \frac{1}{2} \int_{\Gamma^{c}(\mu)} \frac{1}{\gamma} |\sigma_{nn}(\mathbf{v})|^2 d\Gamma(\mu) + \frac{1}{2} \int_{\Gamma^{c}(\mu)} \frac{1}{\gamma} [P_{\gamma,g}^{n}(\mathbf{v})]_{-}^2 d\Gamma(\mu). \quad (19) $$

The problem (16) is nonlinear. To solve it, we use an iterative method. Given $\mathbf{u}_0(\mu) \in \mathcal{V}(\mu)$ and a tolerance $\delta_{\mathbf{u}} \in \mathbb{R}_+$, we look, for all $k \ge 0$, for the solution $\mathbf{u}_{k+1}(\mu)$ in the form $\mathbf{u}_{k+1}(\mu) = \mathbf{u}_k(\mu) + \delta \mathbf{u}_k(\mu)$. Ideally, we seek $\delta \mathbf{u}_k(\mu) \in \mathcal{V}(\mu)$ such that

$$ \left\{ \begin{array}{l} a_{\gamma}^{n}(\mu; \mathbf{u}_k(\mu) + \delta \mathbf{u}_k(\mu), \mathbf{v}) \\ + \int_{\Gamma^c(\mu)} \frac{1}{\gamma} [P_{\gamma,g}^{n}(\mu; \mathbf{u}_k(\mu) + \delta \mathbf{u}_k(\mu))]_{-} P_{\gamma,0}^{n}(\mu; \mathbf{v}) d\Gamma(\mu) = f(\mu; \mathbf{v}), \quad \forall \mathbf{v} \in \mathcal{V}(\mu). \end{array} \right. \quad (20) $$

However, since the problem (20) is nonlinear, it is expensive to solve it directly. Therefore, we approximate the solution $\delta \mathbf{u}_k(\mu)$ by linearizing the problem. In order to do this, we observe that
## 3.1 Frictionless case

Following [12], the starting observation is that the Signorini conditions (3d) can be equivalently reformulated as follows:

$$ \sigma_{nn}(\mathbf{u}(\mu)) = [\sigma_{nn}(\mathbf{u}(\mu)) - \gamma(u_n(\mu) - g(\mu))]_{-}, \quad (13) $$

where $[z]_{-} := \min(z, 0)$ denotes the negative part of a generic real number $z$ and where $\gamma > 0$ is a user-defined parameter (taken to be constant for simplicity). In practice, the parameter $\gamma$ should be chosen large enough (see Section 5 for further discussion).

For all $\mu \in \mathcal{D}$, one introduces the energy functional $J^{\text{Nitsche}}(\mu; \cdot) : \mathcal{V}(\mu) \to \mathbb{R}$ such that

$$ J^{\text{Nitsche}}(\mu; \mathbf{v}) := \mathcal{J}(\mu; \mathbf{v}) - \frac{1}{2} \int_{\Gamma^c(\mu)} \frac{1}{\gamma} |\sigma_{nn}(\mathbf{v})|^2 d\Gamma(\mu) \\ + \frac{1}{2} \int_{\Gamma^c(\mu)} \frac{1}{\gamma} [\sigma_{nn}(\mathbf{v}) - \gamma(v_n - g(\mu))]^2 d\Gamma(\mu), \quad (14) $$

recalling that the energy functional $\mathcal{J}$ is defined in (10). Nitsche's method consists in finding $\mathbf{u}(\mu) \in \mathcal{V}(\mu)$ solution to the following unconstrained minimization problem: Find $\mathbf{u}(\mu) \in \mathcal{V}(\mu)$ such that

$$ \mathbf{u}(\mu) = \underset{\mathbf{v} \in \mathcal{V}(\mu)}{\operatorname{argmin}} J^{\text{Nitsche}}(\mu; \mathbf{v}). \quad (15) $$

The first-order optimality condition associated with (15) reads

$$ a_{\gamma}^{n}(\mu; \mathbf{u}(\mu), \mathbf{v}) + \int_{\Gamma^{c}(\mu)} \frac{1}{\gamma} [P_{\gamma,g}^{n}(\mu; \mathbf{u}(\mu))]_{-} P_{\gamma,0}^{n}(\mu; \mathbf{v}) d\Gamma(\mu) = f(\mu; \mathbf{v}), \quad \forall \mathbf{v} \in \mathcal{V}(\mu), \quad (16) $$

with the bilinear form $a_{\gamma}^{n}(\mu; \cdot, \cdot) : \mathcal{V}(\mu) \times \mathcal{V}(\mu) \to \mathbb{R}$ defined as

$$ a_{\gamma}^{n}(\mu; \mathbf{u}, \mathbf{v}) := a(\mu; \mathbf{u}, \mathbf{v}) - \int_{\Gamma^{c}(\mu)} \frac{1}{\gamma} \sigma_{nn}(\mathbf{u}) \sigma_{nn}(\mathbf{v}) d\Gamma(\mu), \quad \forall \mathbf{u}, \mathbf{v} \in \mathcal{V}(\mu), \quad (17) $$

and the operators $P_{\gamma,g}^{n}(\mu; \cdot, \cdot)$, $P_{\gamma,0}^{n}(\mu; \cdot, \cdot) : \mathcal{V}(\mu) \to L^2(\Gamma^c(\mu))$ defined as

$$ P_{\gamma,g}^{n}(\mu; \mathbf{v}) := \sigma_{nn}(\mathbf{v}) - \gamma(v_n - g(\mu)), \quad (18a) $$

$$ P_{\gamma,0}^{n}(\mu; \mathbf{v}) := \sigma_{nn}(\mathbf{v}) - \gamma v_n. \quad (18b) $$

With this notation, we can rewrite the energy functional $J^{\text{Nitsche}}(\mu; \cdot)$ as

$$ J^{\text{Nitsche}}(\mu; \mathbf{v}) := \mathcal{J}(\mu; \mathbf{v}) - \frac{1}{2} \int_{\Gamma^{c}(\mu)} \frac{1}{\gamma} |\sigma_{nn}(\mathbf{v})|^2 d\Gamma(\mu) + \frac{1}{2} \int_{\Gamma^{c}(\mu)} \frac{1}{\gamma} [P_{\gamma,g}^{n}(\mathbf{v})]_{-}^2 d\Gamma(\mu). \quad (19) $$

The problem (16) is nonlinear. To solve it, we use an iterative method. Given $\mathbf{u}_0(\mu) \in \mathcal{V}(\mu)$ and a tolerance $\delta_{\mathbf{u}} \in \mathbb{R}_+$, we look, for all $k \ge 0$, for the solution $\mathbf{u}_{k+1}(\mu)$ in the form $\mathbf{u}_{k+1}(\mu) = \mathbf{u}_k(\mu) + \delta \mathbf{u}_k(\mu)$. Ideally, we seek $\delta \mathbf{u}_k(\mu) \in \mathcal{V}(\mu)$ such that

$$ \left\{ \begin{array}{l} a_{\gamma}^{n}(\mu; \mathbf{u}_k(\mu) + \delta \mathbf{u}_k(\mu), \mathbf{v}) \\ + \int_{\Gamma^c(\mu)} \frac{1}{\gamma} [P_{\gamma,g}^{n}(\mu; \mathbf{u}_k(\mu) + \delta \mathbf{u}_k(\mu))]_{-} P_{\gamma,0}^{n}(\mu; \mathbf{v}) d\Gamma(\mu) = f(\mu; \mathbf{v}), \quad \forall \mathbf{v} \in \mathcal{V}(\mu). \end{array} \right. \quad (20) $$

However, since the problem (20) is nonlinear, it is expensive to solve it directly. Therefore, we approximate the solution $\delta \mathbf{u}_k(\mu)$ by linearizing the problem. In order to do this, we observe that
$P_{\gamma,g}^n(\mu; \mathbf{u}_k(\mu) + \delta\mathbf{u}_k(\mu)) = P_{\gamma,g}^n(\mu; \mathbf{u}_k(\mu)) + P_{\gamma,0}^n(\mu; \delta\mathbf{u}_k(\mu))$ and approximate the term $[P_{\gamma,g}^n(\mu; \mathbf{u}_k(\mu) + \delta\mathbf{u}_k(\mu))]_-$ as follows:

$$[P_{\gamma,g}^n(\mu; \mathbf{u}_k(\mu) + \delta\mathbf{u}_k(\mu))]_- \approx [P_{\gamma,g}^n(\mu; \mathbf{u}_k(\mu))]_- + H(-P_{\gamma,g}^n(\mu; \mathbf{u}_k(\mu)))P_{\gamma,0}^n(\mu; \delta\mathbf{u}_k(\mu)), \quad (21)$$

where $H(\cdot)$ denotes the Heaviside function. Using this approximation in (20) (for simplicity, we keep the same notation for the unknown $\delta\mathbf{u}_k(\mu)$), we consider the following sequence of problems: For all $k \ge 0$, find $\delta\mathbf{u}_k(\mu) \in \mathcal{V}(\mu)$ such that

$$a_{\gamma}^{n}(\mu; \delta \mathbf{u}_{k}(\mu), \mathbf{v}) + b_{\gamma}^{n}(\mu; \mathbf{u}_{k}(\mu); \delta \mathbf{u}_{k}(\mu), \mathbf{v}) = -r_{\gamma}^{n}(\mu; \mathbf{u}_{k}(\mu); \mathbf{v}), \quad \forall \mathbf{v} \in \mathcal{V}(\mu), \quad (22)$$

where

$$b_{\gamma}^{n}(\mu; \mathbf{w}; \mathbf{u}, \mathbf{v}) := \int_{\Gamma^{c}(\mu)} \frac{1}{\gamma} H(-P_{\gamma,g}^{n}(\mu; \mathbf{w})) P_{\gamma,0}^{n}(\mu; \mathbf{u}) P_{\gamma,0}^{n}(\mu; \mathbf{v}) d\Gamma(\mu), \quad (23a)$$

$$r_{\gamma}^{n}(\mu; \mathbf{w}; \mathbf{v}) := a_{\gamma}^{n}(\mu; \mathbf{w}, \mathbf{v}) + \theta_{\gamma}^{n}(\mu; \mathbf{w}; \mathbf{v}) - f(\mu; \mathbf{v}), \quad (23b)$$

with

$$\theta_{\gamma}^{n}(\mu; \mathbf{w}, \mathbf{v}) := \int_{\Gamma^{c}(\mu)} \frac{1}{\gamma} [P_{\gamma,g}^{n}(\mu; \mathbf{w})]_- P_{\gamma,0}^{n}(\mu; \mathbf{v}) d\Gamma(\mu). \quad (24)$$

We iterate on $k$ until the following convergence criterion is reached:

$$\frac{\|\mathbf{u}_{k+1}(\mu) - \mathbf{u}_k(\mu)\|_{\mathcal{V}(\mu)}}{\|\mathbf{u}_{k+1}(\mu)\|_{\mathcal{V}(\mu)}} \le \delta_{\mathbf{u}}. \quad (25)$$

In what follows, we denote by $k^{cv}(\mu) \in \mathbb{N}^*$ the number of iterations required to solve (22), for all $\mu \in \mathcal{D}$. We denote the converged solution to the sequence of problems (22) as $\mathbf{u}_{cv}(\mu) \in \mathcal{V}(\mu)$.

To introduce the algebraic formulation, we assume that for all $\mu \in \mathcal{D}$, the HF (finite-dimensional) space $\mathcal{V}(\mu)$ is such that

$$\mathcal{V}(\mu) := \operatorname{Span}(\{\varphi_n(\mu)\}_{n \in \{1:\mathcal{N}^{\mathrm{HF}}\}}). \quad (26)$$

Notice that the dimension $\mathcal{N}^{\mathrm{HF}}$ of $\mathcal{V}(\mu)$ 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):

$$\mathbf{u}_k(\mu) := \sum_{n \in \{1:\mathcal{N}^{\mathrm{HF}}\}} \mathbf{u}_k^n(\mu)\varphi_n(\mu), \quad \mathbf{U}_k(\mu) := (\mathbf{u}_k^n(\mu))_{n \in \{1:\mathcal{N}^{\mathrm{HF}}\}} \in \mathbb{R}^{\mathcal{N}^{\mathrm{HF}}}, \quad (27a)$$

$$\delta \mathbf{u}_k(\mu) := \sum_{n \in \{1:\mathcal{N}^{\mathrm{HF}}\}} \delta \mathbf{u}_k^n(\mu) \varphi_n(\mu), \quad \Delta \mathbf{U}_k(\mu) := (\delta \mathbf{u}_k^n(\mu))_{n \in \{1:\mathcal{N}^{\mathrm{HF}}\}} \in \mathbb{R}^{\mathcal{N}^{\mathrm{HF}}}. \quad (27b)$$

The algebraic formulation of the sequence of problems (22) then reads as follows: For all $k \ge 0$, find $\Delta \mathbf{U}_k(\mu) \in \mathbb{R}^{\mathcal{N}^{\mathrm{HF}}}$ such that

$$A_{\gamma}^{n}(\mu) \Delta \mathbf{U}_{k}(\mu) + B_{\gamma}^{n}(\mu, \mathbf{u}_{k}(\mu)) \Delta \mathbf{U}_{k}(\mu) = -R_{\gamma}^{n}(\mu, \mathbf{u}_{k}(\mu)), \quad (28)$$

where for all $n, m \in \{1:\mathcal{N}^{\mathrm{HF}}\}$,

$$(A_{\gamma}^{n}(\mu))_{mn} := a_{\gamma}^{n}(\mu; \varphi_n(\mu), \varphi_m(\mu)), \quad (29a)$$

$$(B_{\gamma}^{n}(\mu, \mathbf{u}_{k}(\mu)))_{mn} := b_{\gamma}^{n}(\mu; \mathbf{u}_{k}(\mu); \varphi_{n}(\mu), \varphi_{m}(\mu)), \quad (29b)$$

$$(R_{\gamma}^{n}(\mu, \mathbf{u}_{k}(\mu)))_{m} := r_{\gamma}^{n}(\mu; \mathbf{u}_{k}(\mu); \varphi_{m}(\mu)). \quad (29c)$$
$P_{\gamma,g}^n(\mu; \mathbf{u}_k(\mu) + \delta\mathbf{u}_k(\mu)) = P_{\gamma,g}^n(\mu; \mathbf{u}_k(\mu)) + P_{\gamma,0}^n(\mu; \delta\mathbf{u}_k(\mu))$ and approximate the term $[P_{\gamma,g}^n(\mu; \mathbf{u}_k(\mu) + \delta\mathbf{u}_k(\mu))]_-$ as follows:

$$[P_{\gamma,g}^n(\mu; \mathbf{u}_k(\mu) + \delta\mathbf{u}_k(\mu))]_- \approx [P_{\gamma,g}^n(\mu; \mathbf{u}_k(\mu))]_- + H(-P_{\gamma,g}^n(\mu; \mathbf{u}_k(\mu)))P_{\gamma,0}^n(\mu; \delta\mathbf{u}_k(\mu)), \quad (21)$$

where $H(\cdot)$ denotes the Heaviside function. Using this approximation in (20) (for simplicity, we keep the same notation for the unknown $\delta\mathbf{u}_k(\mu)$), we consider the following sequence of problems: For all $k \ge 0$, find $\delta\mathbf{u}_k(\mu) \in \mathcal{V}(\mu)$ such that

$$a_{\gamma}^{n}(\mu; \delta \mathbf{u}_{k}(\mu), \mathbf{v}) + b_{\gamma}^{n}(\mu; \mathbf{u}_{k}(\mu); \delta \mathbf{u}_{k}(\mu), \mathbf{v}) = -r_{\gamma}^{n}(\mu; \mathbf{u}_{k}(\mu); \mathbf{v}), \quad \forall \mathbf{v} \in \mathcal{V}(\mu), \quad (22)$$

where

$$b_{\gamma}^{n}(\mu; \mathbf{w}; \mathbf{u}, \mathbf{v}) := \int_{\Gamma^{c}(\mu)} \frac{1}{\gamma} H(-P_{\gamma,g}^{n}(\mu; \mathbf{w})) P_{\gamma,0}^{n}(\mu; \mathbf{u}) P_{\gamma,0}^{n}(\mu; \mathbf{v}) d\Gamma(\mu), \quad (23a)$$

$$r_{\gamma}^{n}(\mu; \mathbf{w}; \mathbf{v}) := a_{\gamma}^{n}(\mu; \mathbf{w}, \mathbf{v}) + \theta_{\gamma}^{n}(\mu; \mathbf{w}; \mathbf{v}) - f(\mu; \mathbf{v}), \quad (23b)$$

with

$$\theta_{\gamma}^{n}(\mu; \mathbf{w}, \mathbf{v}) := \int_{\Gamma^{c}(\mu)} \frac{1}{\gamma} [P_{\gamma,g}^{n}(\mu; \mathbf{w})]_- P_{\gamma,0}^{n}(\mu; \mathbf{v}) d\Gamma(\mu). \quad (24)$$

We iterate on $k$ until the following convergence criterion is reached:

$$\frac{\|\mathbf{u}_{k+1}(\mu) - \mathbf{u}_k(\mu)\|_{\mathcal{V}(\mu)}}{\|\mathbf{u}_{k+1}(\mu)\|_{\mathcal{V}(\mu)}} \le \delta_{\mathbf{u}}. \quad (25)$$

In what follows, we denote by $k^{cv}(\mu) \in \mathbb{N}^*$ the number of iterations required to solve (22), for all $\mu \in \mathcal{D}$. We denote the converged solution to the sequence of problems (22) as $\mathbf{u}_{cv}(\mu) \in \mathcal{V}(\mu)$.

To introduce the algebraic formulation, we assume that for all $\mu \in \mathcal{D}$, the HF (finite-dimensional) space $\mathcal{V}(\mu)$ is such that

$$\mathcal{V}(\mu) := \operatorname{Span}(\{\varphi_n(\mu)\}_{n \in \{1:\mathcal{N}^{\mathrm{HF}}\}}). \quad (26)$$

Notice that the dimension $\mathcal{N}^{\mathrm{HF}}$ of $\mathcal{V}(\mu)$ 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):

$$\mathbf{u}_k(\mu) := \sum_{n \in \{1:\mathcal{N}^{\mathrm{HF}}\}} \mathbf{u}_k^n(\mu)\varphi_n(\mu), \quad \mathbf{U}_k(\mu) := (\mathbf{u}_k^n(\mu))_{n \in \{1:\mathcal{N}^{\mathrm{HF}}\}} \in \mathbb{R}^{\mathcal{N}^{\mathrm{HF}}}, \quad (27a)$$

$$\delta \mathbf{u}_k(\mu) := \sum_{n \in \{1:\mathcal{N}^{\mathrm{HF}}\}} \delta \mathbf{u}_k^n(\mu) \varphi_n(\mu), \quad \Delta \mathbf{U}_k(\mu) := (\delta \mathbf{u}_k^n(\mu))_{n \in \{1:\mathcal{N}^{\mathrm{HF}}\}} \in \mathbb{R}^{\mathcal{N}^{\mathrm{HF}}}. \quad (27b)$$

The algebraic formulation of the sequence of problems (22) then reads as follows: For all $k \ge 0$, find $\Delta \mathbf{U}_k(\mu) \in \mathbb{R}^{\mathcal{N}^{\mathrm{HF}}}$ such that

$$A_{\gamma}^{n}(\mu) \Delta \mathbf{U}_{k}(\mu) + B_{\gamma}^{n}(\mu, \mathbf{u}_{k}(\mu)) \Delta \mathbf{U}_{k}(\mu) = -R_{\gamma}^{n}(\mu, \mathbf{u}_{k}(\mu)), \quad (28)$$

where for all $n, m \in \{1:\mathcal{N}^{\mathrm{HF}}\}$,

$$(A_{\gamma}^{n}(\mu))_{mn} := a_{\gamma}^{n}(\mu; \varphi_n(\mu), \varphi_m(\mu)), \quad (29a)$$

$$(B_{\gamma}^{n}(\mu, \mathbf{u}_{k}(\mu)))_{mn} := b_{\gamma}^{n}(\mu; \mathbf{u}_{k}(\mu); \varphi_{n}(\mu), \varphi_{m}(\mu)), \quad (29b)$$

$$(R_{\gamma}^{n}(\mu, \mathbf{u}_{k}(\mu)))_{m} := r_{\gamma}^{n}(\mu; \mathbf{u}_{k}(\mu); \varphi_{m}(\mu)). \quad (29c)$$
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}(\mu) \in \mathbb{R}^{\mathcal{N}^{\text{HF}}}$. Furthermore, let us denote by $\Theta_{\gamma}^{n}(\mu, \mathbf{w})$ (resp. $F(\mu)$) the algebraic representation of the operator $\theta_{\gamma}^{n}(\mu; \mathbf{w}; \cdot)$ (resp. $f(\mu; \cdot)$) such that for all $m \in \{1:\mathcal{N}^{\text{HF}}\}$,

$$ (\Theta_{\gamma}^{n}(\mu, \mathbf{w}))_{m} := \theta_{\gamma}^{n}(\mu; \mathbf{w}; \varphi_{m}(\mu)), \quad (F(\mu))_{m} := f(\mu; \varphi_{m}(\mu)). \qquad (30) $$

With this notation, we have the following decomposition:

$$ R_{\gamma}^{n}(\mu, \mathbf{u}_{k}(\mu)) := A_{\gamma}^{n}(\mu)U_{k}(\mu) + \Theta_{\gamma}^{n}(\mu, \mathbf{u}_{k}(\mu)) - F(\mu). \qquad (31) $$

## 3.2 Friction case

For the Tresca frictional contact problem, we use the following reformulation of the friction conditions given in [8]:

$$ \sigma_{n\tau}(\mathbf{u}(\mu)) = [\sigma_{n\tau}(\mathbf{u}(\mu)) - \gamma\mathbf{u}_{\tau}(\mu)]_{s}, \qquad (32) $$

where, for a positive real number $r$, the notation

$$ [\mathbf{x}]_r := \begin{cases} \mathbf{x}, & \text{if } \|\mathbf{x}\| \le r, \\ r \frac{\mathbf{x}}{\|\mathbf{x}\|}, & \text{otherwise,} \end{cases} \qquad (33) $$

defines the projection of $\mathbf{x} \in \mathbb{R}^{d-1}$ onto the ball centered at the origin and of radius $r$. Let us introduce the operator $P_{\gamma}^{\tau}(\mu; \cdot) : \mathcal{V}(\mu) \to L^{2}(\Gamma^{\text{c}}(\mu); \mathbb{R}^{d-1})$ such that

$$ P_{\gamma}^{\tau}(\mu; \mathbf{v}) := \sigma_{n\tau}(\mathbf{v}) - \gamma \mathbf{v}_{\tau}. \qquad (34) $$

We start from the energy functional $J_{\mathcal{F}}^{\text{Nitsche}}(\mu; \cdot) : \mathcal{V}(\mu) \to \mathbb{R}$ defined as

$$ J_{\mathcal{F}}^{\text{Nitsche}}(\mu; \mathbf{v}) := J^{\text{Nitsche}}(\mu; \mathbf{v}) - \frac{1}{2} \int_{\Gamma^{\text{c}}(\mu)} \frac{1}{\gamma} \|\boldsymbol{\sigma}_{n\tau}(\mathbf{v})\|^2 d\Gamma(\mu) \\ + \frac{1}{2} \int_{\Gamma^{\text{c}}(\mu)} \frac{1}{\gamma} \|P_{\gamma}^{\tau}(\mathbf{v})\|^2 d\Gamma(\mu) - \frac{1}{2} \int_{\Gamma^{\text{c}}(\mu)} \frac{1}{\gamma} \|P_{\gamma}^{\tau}(\mathbf{v}) - [P_{\gamma}^{\tau}(\mathbf{v})]_{s}\|^2 d\Gamma(\mu). \qquad (35) $$

Nitsche's method consists in solving the following unconstrained minimization problem: Find $\mathbf{u}(\mu) \in \mathcal{V}(\mu)$ such that

$$ \mathbf{u}(\mu) = \underset{\mathbf{v} \in \mathcal{V}(\mu)}{\operatorname{argmin}} J_{\mathcal{F}}^{\text{Nitsche}}(\mu; \mathbf{v}). \qquad (36) $$

Let us consider the functional $J: \mathbb{R}^{d-1} \to \mathbb{R}$ defined as follows:

$$ J(\mathbf{x}) = \frac{1}{2} \|\mathbf{x} - [\mathbf{x}]_s\|^2. \qquad (37) $$

One easily verifies that this functional is Gâteaux-differentiable and that its differential is given by $DJ(\mathbf{x}) = \mathbf{x} - [\mathbf{x}]_s$ for all $\mathbf{x} \in \mathbb{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 $\mathbf{u}(\mu) \in \mathcal{V}(\mu)$ such that

$$ \left\{ \begin{aligned} & a_{\gamma}^{n\tau}(\mu; \mathbf{u}(\mu), \mathbf{v}) + \int_{\Gamma^{\text{c}}(\mu)} \frac{1}{\gamma} [P_{\gamma,g}^{n}(\mu; \mathbf{u}(\mu))]_{s} P_{\gamma,0}^{n}(\mu; \mathbf{v}) d\Gamma(\mu) \\ & + \int_{\Gamma^{\text{c}}(\mu)} \frac{1}{\gamma} [P_{\gamma}^{\tau}(\mu; \mathbf{u}(\mu))]_{s} P_{\gamma}^{\tau}(\mu; \mathbf{v}) d\Gamma(\mu) = f(\mu; \mathbf{v}), \quad \forall \mathbf{v} \in \mathcal{V}(\mu), \end{aligned} \right. \qquad (38) $$
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}(\mu) \in \mathbb{R}^{\mathcal{N}^{\text{HF}}}$. Furthermore, let us denote by $\Theta_{\gamma}^{n}(\mu, \mathbf{w})$ (resp. $F(\mu)$) the algebraic representation of the operator $\theta_{\gamma}^{n}(\mu; \mathbf{w}; \cdot)$ (resp. $f(\mu; \cdot)$) such that for all $m \in \{1:\mathcal{N}^{\text{HF}}\}$,

$$ (\Theta_{\gamma}^{n}(\mu, \mathbf{w}))_{m} := \theta_{\gamma}^{n}(\mu; \mathbf{w}; \varphi_{m}(\mu)), \quad (F(\mu))_{m} := f(\mu; \varphi_{m}(\mu)). \quad (30) $$

With this notation, we have the following decomposition:

$$ R_{\gamma}^{n}(\mu, \mathbf{u}_{k}(\mu)) := A_{\gamma}^{n}(\mu)U_{k}(\mu) + \Theta_{\gamma}^{n}(\mu, \mathbf{u}_{k}(\mu)) - F(\mu). \quad (31) $$

## 3.2 Friction case

For the Tresca frictional contact problem, we use the following reformulation of the friction conditions given in [8]:

$$ \sigma_{n\tau}(\mathbf{u}(\mu)) = [\sigma_{n\tau}(\mathbf{u}(\mu)) - \gamma\mathbf{u}_{\tau}(\mu)]_{s}, \quad (32) $$

where, for a positive real number $r$, the notation

$$ [\mathbf{x}]_r := \begin{cases} \mathbf{x}, & \text{if } \|\mathbf{x}\| \le r, \\ r \frac{\mathbf{x}}{\|\mathbf{x}\|}, & \text{otherwise,} \end{cases} \quad (33) $$

defines the projection of $\mathbf{x} \in \mathbb{R}^{d-1}$ onto the ball centered at the origin and of radius $r$. Let us introduce the operator $P_{\gamma}^{\tau}(\mu; \cdot) : \mathcal{V}(\mu) \to L^{2}(\Gamma^{\text{c}}(\mu); \mathbb{R}^{d-1})$ such that

$$ P_{\gamma}^{\tau}(\mu; \mathbf{v}) := \sigma_{n\tau}(\mathbf{v}) - \gamma \mathbf{v}_{\tau}. \quad (34) $$

We start from the energy functional $J_{\mathcal{F}}^{\text{Nitsche}}(\mu; \cdot) : \mathcal{V}(\mu) \to \mathbb{R}$ defined as

$$ \begin{aligned} J_{\mathcal{F}}^{\text{Nitsche}}(\mu; \mathbf{v}) := & J^{\text{Nitsche}}(\mu; \mathbf{v}) - \frac{1}{2} \int_{\Gamma^{\text{c}}(\mu)} \frac{1}{\gamma} \|\boldsymbol{\sigma}_{n\tau}(\mathbf{v})\|^2 d\Gamma(\mu) \\ & + \frac{1}{2} \int_{\Gamma^{\text{c}}(\mu)} \frac{1}{\gamma} \|P_{\gamma}^{\tau}(\mathbf{v})\|^2 d\Gamma(\mu) - \frac{1}{2} \int_{\Gamma^{\text{c}}(\mu)} \frac{1}{\gamma} \|P_{\gamma}^{\tau}(\mathbf{v}) - [P_{\gamma}^{\tau}(\mathbf{v})]_{s}\|^2 d\Gamma(\mu). \end{aligned} \quad (35) $$

Nitsche's method consists in solving the following unconstrained minimization problem: Find $\mathbf{u}(\mu) \in \mathcal{V}(\mu)$ such that

$$ \mathbf{u}(\mu) = \underset{\mathbf{v} \in \mathcal{V}(\mu)}{\operatorname{argmin}} J_{\mathcal{F}}^{\text{Nitsche}}(\mu; \mathbf{v}). \quad (36) $$

Let us consider the functional $J: \mathbb{R}^{d-1} \to \mathbb{R}$ defined as follows:

$$ J(\mathbf{x}) = \frac{1}{2} \|\mathbf{x} - [\mathbf{x}]_s\|^2. \quad (37) $$

One easily verifies that this functional is Gâteaux-differentiable and that its differential is given by $DJ(\mathbf{x}) = \mathbf{x} - [\mathbf{x}]_s$ for all $\mathbf{x} \in \mathbb{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 $\mathbf{u}(\mu) \in \mathcal{V}(\mu)$ such that

$$ \left\{ \begin{aligned} & a_{\gamma}^{n\tau}(\mu; \mathbf{u}(\mu), \mathbf{v}) + \int_{\Gamma^{\text{c}}(\mu)} \frac{1}{\gamma} [P_{\gamma,g}^{n}(\mu; \mathbf{u}(\mu))]_{s} P_{\gamma,0}^{n}(\mu; \mathbf{v}) d\Gamma(\mu) \\ & + \int_{\Gamma^{\text{c}}(\mu)} \frac{1}{\gamma} [P_{\gamma}^{\tau}(\mu; \mathbf{u}(\mu))]_{s} P_{\gamma}^{\tau}(\mu; \mathbf{v}) d\Gamma(\mu) = f(\mu; \mathbf{v}), \quad \forall \mathbf{v} \in \mathcal{V}(\mu), \end{aligned} \right. \quad (38) $$
where the bilinear form $a_{\gamma}^{n\tau}(\mu; \cdot, \cdot) : \mathcal{V}(\mu) \times \mathcal{V}(\mu) \to \mathbb{R}$ is defined as

$$a_{\gamma}^{n\tau}(\mu; \mathbf{u}, \mathbf{v}) := a(\mu; \mathbf{u}, \mathbf{v}) - \int_{\Gamma^c(\mu)} \frac{1}{\gamma} \sigma(\mathbf{u}) \mathbf{n} \cdot \sigma(\mathbf{v}) \mathbf{n} \, d\Gamma(\mu). \quad (39)$$

The linearization of the problem (38) leads to the following sequence of problems: For all $k \ge 0$, find $\delta \mathbf{u}_k(\mu) \in \mathcal{V}(\mu)$ such that

$$a_{\gamma}^{n\tau}(\mu; \delta \mathbf{u}_k(\mu), \mathbf{v}) + b_{\gamma}^{n\tau}(\mu; \mathbf{u}_k(\mu); \delta \mathbf{u}_k(\mu), \mathbf{v}) = -r_{\gamma}^{n\tau}(\mu; \mathbf{u}_k(\mu); \mathbf{v}), \quad \forall \mathbf{v} \in \mathcal{V}(\mu), \quad (40)$$

with

$$b_{\gamma}^{n\tau}(\mu; \mathbf{w}; \mathbf{u}, \mathbf{v}) := b_{\gamma}^{n}(\mu; \mathbf{w}; \mathbf{u}, \mathbf{v}) + b_{\gamma}^{\tau}(\mu; \mathbf{w}; \mathbf{u}, \mathbf{v}), \quad (41a)$$

$$r_{\gamma}^{n\tau}(\mu; \mathbf{w}; \mathbf{v}) := a_{\gamma}^{n\tau}(\mu; \mathbf{w}, \mathbf{v}) + \theta_{\gamma}^{n}(\mu; \mathbf{w}, \mathbf{v}) + \theta_{\gamma}^{\tau}(\mu; \mathbf{w}, \mathbf{v}) - f(\mu; \mathbf{v}), \quad (41b)$$

where

$$b_{\gamma}^{\tau}(\mu; \mathbf{w}; \mathbf{u}, \mathbf{v}) := \int_{\Gamma^c(\mu)} \frac{1}{\gamma} (\mathcal{G}_s(P_{\gamma}^{\tau}(\mu; \mathbf{w})) \cdot P_{\gamma}^{\tau}(\mu; \mathbf{u})) P_{\gamma}^{\tau}(\mu; \mathbf{v}) \, d\Gamma(\mu), \quad (42a)$$

$$\theta_{\gamma}^{\tau}(\mu; \mathbf{w}, \mathbf{v}) := \int_{\Gamma^c(\mu)} \frac{1}{\gamma} [P_{\gamma}^{\tau}(\mu; \mathbf{w})]_s P_{\gamma}^{\tau}(\mu; \mathbf{v}) \, d\Gamma(\mu), \quad (42b)$$

and $\mathcal{G}_s(\cdot)$, the differential of $[\cdot]_s$, is given by

$$\mathcal{G}_s(\mathbf{x}) := \begin{cases} \mathbb{1}_{d-1}, & \text{if } \|\mathbf{x}\| \le s, \\ \frac{s}{\|\mathbf{x}\|} (\mathbb{1}_{d-1} - \frac{\mathbf{x} \otimes \mathbf{x}}{\|\mathbf{x}\|^2}), & \text{otherwise,} \end{cases} \quad (43)$$

with $\mathbb{1}_{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 $\mathbf{u}_{cv}(\mu) \in \mathcal{V}(\mu)$.

The algebraic formulation of the sequence of problems (40) reads as follows: For all $k \ge 0$, find $\mathbf{\Delta} \mathbf{U}_k(\mu) \in \mathbb{R}^{\mathcal{N}^{\text{HF}}}$ such that

$$A_{\gamma}^{n\tau}(\mu) \mathbf{\Delta} \mathbf{U}_k(\mu) + B_{\gamma}^{n\tau}(\mu, \mathbf{u}_k(\mu)) \mathbf{\Delta} \mathbf{U}_k(\mu) = -R_{\gamma}^{n\tau}(\mu, \mathbf{u}_k(\mu)), \quad (44)$$

where for all $n, m \in \{1:\mathcal{N}^{\text{HF}}\}$,

$$(A_{\gamma}^{n\tau}(\mu))_{mn} := a_{\gamma}^{n\tau}(\mu; \varphi_n(\mu), \varphi_m(\mu)), \quad (45a)$$

$$(B_{\gamma}^{n\tau}(\mu, \mathbf{u}_k(\mu)))_{mn} := b_{\gamma}^{n\tau}(\mu; \mathbf{u}_k(\mu); \varphi_n(\mu), \varphi_m(\mu)), \quad (45b)$$

$$(R_{\gamma}^{n\tau}(\mu, \mathbf{u}_k(\mu)))_m := r_{\gamma}^{n\tau}(\mu; \mathbf{u}_k(\mu); \varphi_m(\mu)). \quad (45c)$$

The convergence criterion for (44) is still (25) using the reconstructed functions. We denote the converged solution to the sequence of problems (44) as $\mathbf{U}_{cv}(\mu) \in \mathbb{R}^{\mathcal{N}^{\text{HF}}}$ . Furthermore, let us denote by $B_{\gamma}^{\tau}(\mu, \mathbf{w})$ (resp. $\Theta_{\gamma}^{\tau}(\mu, \mathbf{w})$) the algebraic representation of the operator $b_{\gamma}^{\tau}(\mu; \mathbf{w}; \cdot, \cdot)$ (resp. $\theta_{\gamma}^{\tau}(\mu; \mathbf{w}; \cdot)$) such that for all $n, m \in \{1:\mathcal{N}^{\text{HF}}\}$,

$$(B_{\gamma}^{\tau}(\mu, \mathbf{w}))_{mn} := b_{\gamma}^{\tau}(\mu; \mathbf{w}; \varphi_n(\mu), \varphi_m(\mu)), \quad (\Theta_{\gamma}^{\tau}(\mu, \mathbf{w}))_m := \theta_{\gamma}^{\tau}(\mu; \mathbf{w}; \varphi_m(\mu)). \quad (46)$$

With this notation, we have the following decompositions:

$$B_{\gamma}^{n\tau}(\mu, \mathbf{u}_k(\mu)) := B_{\gamma}^{n}(\mu, \mathbf{u}_k(\mu)) + B_{\gamma}^{\tau}(\mu, \mathbf{u}_k(\mu)), \quad (47a)$$

$$R_{\gamma}^{n\tau}(\mu, \mathbf{u}_k(\mu)) := A_{\gamma}^{n\tau}(\mu)\mathbf{U}_k(\mu) + \Theta_{\gamma}^{n}(\mu, \mathbf{u}_k(\mu)) + \Theta_{\gamma}^{\tau}(\mu, \mathbf{u}_k(\mu)) - F(\mu). \quad (47b)$$
where the bilinear form $a_{\gamma}^{n\tau}(\mu; \cdot, \cdot) : \mathcal{V}(\mu) \times \mathcal{V}(\mu) \to \mathbb{R}$ is defined as

$$a_{\gamma}^{n\tau}(\mu; \mathbf{u}, \mathbf{v}) := a(\mu; \mathbf{u}, \mathbf{v}) - \int_{\Gamma^c(\mu)} \frac{1}{\gamma} \sigma(\mathbf{u}) \mathbf{n} \cdot \sigma(\mathbf{v}) \mathbf{n} \, d\Gamma(\mu). \quad (39)$$

The linearization of the problem (38) leads to the following sequence of problems: For all $k \ge 0$, find $\delta \mathbf{u}_k(\mu) \in \mathcal{V}(\mu)$ such that

$$a_{\gamma}^{n\tau}(\mu; \delta \mathbf{u}_k(\mu), \mathbf{v}) + b_{\gamma}^{n\tau}(\mu; \mathbf{u}_k(\mu); \delta \mathbf{u}_k(\mu), \mathbf{v}) = -r_{\gamma}^{n\tau}(\mu; \mathbf{u}_k(\mu); \mathbf{v}), \quad \forall \mathbf{v} \in \mathcal{V}(\mu), \quad (40)$$

with

$$b_{\gamma}^{n\tau}(\mu; \mathbf{w}; \mathbf{u}, \mathbf{v}) := b_{\gamma}^{n}(\mu; \mathbf{w}; \mathbf{u}, \mathbf{v}) + b_{\gamma}^{\tau}(\mu; \mathbf{w}; \mathbf{u}, \mathbf{v}), \quad (41a)$$

$$r_{\gamma}^{n\tau}(\mu; \mathbf{w}; \mathbf{v}) := a_{\gamma}^{n\tau}(\mu; \mathbf{w}, \mathbf{v}) + \theta_{\gamma}^{n}(\mu; \mathbf{w}, \mathbf{v}) + \theta_{\gamma}^{\tau}(\mu; \mathbf{w}, \mathbf{v}) - f(\mu; \mathbf{v}), \quad (41b)$$

where

$$b_{\gamma}^{\tau}(\mu; \mathbf{w}; \mathbf{u}, \mathbf{v}) := \int_{\Gamma^c(\mu)} \frac{1}{\gamma} (\mathcal{G}_s(P_{\gamma}^{\tau}(\mu; \mathbf{w})) \cdot P_{\gamma}^{\tau}(\mu; \mathbf{u})) P_{\gamma}^{\tau}(\mu; \mathbf{v}) \, d\Gamma(\mu), \quad (42a)$$

$$\theta_{\gamma}^{\tau}(\mu; \mathbf{w}, \mathbf{v}) := \int_{\Gamma^c(\mu)} \frac{1}{\gamma} [P_{\gamma}^{\tau}(\mu; \mathbf{w})]_s P_{\gamma}^{\tau}(\mu; \mathbf{v}) \, d\Gamma(\mu), \quad (42b)$$

and $\mathcal{G}_s(\cdot)$, the differential of $[\cdot]_s$, is given by

$$\mathcal{G}_s(\mathbf{x}) := \begin{cases} \mathbb{1}_{d-1}, & \text{if } \|\mathbf{x}\| \le s, \\ \frac{s}{\|\mathbf{x}\|} (\mathbb{1}_{d-1} - \frac{\mathbf{x} \otimes \mathbf{x}}{\|\mathbf{x}\|^2}), & \text{otherwise,} \end{cases} \quad (43)$$

with $\mathbb{1}_{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 $\mathbf{u}_{cv}(\mu) \in \mathcal{V}(\mu)$.

The algebraic formulation of the sequence of problems (40) reads as follows: For all $k \ge 0$, find $\mathbf{\Delta} \mathbf{U}_k(\mu) \in \mathbb{R}^{\mathcal{N}^{\text{HF}}}$ such that

$$A_{\gamma}^{n\tau}(\mu) \mathbf{\Delta} \mathbf{U}_k(\mu) + B_{\gamma}^{n\tau}(\mu, \mathbf{u}_k(\mu)) \mathbf{\Delta} \mathbf{U}_k(\mu) = -R_{\gamma}^{n\tau}(\mu, \mathbf{u}_k(\mu)), \quad (44)$$

where for all $n, m \in \{1: \mathcal{N}^{\text{HF}}\}$,

$$(A_{\gamma}^{n\tau}(\mu))_{mn} := a_{\gamma}^{n\tau}(\mu; \varphi_n(\mu), \varphi_m(\mu)), \quad (45a)$$

$$(B_{\gamma}^{n\tau}(\mu, \mathbf{u}_k(\mu)))_{mn} := b_{\gamma}^{n\tau}(\mu; \mathbf{u}_k(\mu); \varphi_n(\mu), \varphi_m(\mu)), \quad (45b)$$

$$(R_{\gamma}^{n\tau}(\mu, \mathbf{u}_k(\mu)))_m := r_{\gamma}^{n\tau}(\mu; \mathbf{u}_k(\mu); \varphi_m(\mu)). \quad (45c)$$

The convergence criterion for (44) is still (25) using the reconstructed functions. We denote the converged solution to the sequence of problems (44) as $\mathbf{U}_{cv}(\mu) \in \mathbb{R}^{\mathcal{N}^{\text{HF}}}$. Furthermore, let us denote by $B_{\gamma}^{\tau}(\mu, \mathbf{w})$ (resp. $\Theta_{\gamma}^{\tau}(\mu, \mathbf{w})$) the algebraic representation of the operator $b_{\gamma}^{\tau}(\mu; \mathbf{w}; \cdot, \cdot)$ (resp. $\theta_{\gamma}^{\tau}(\mu; \mathbf{w}; \cdot)$) such that for all $n, m \in \{1: \mathcal{N}^{\text{HF}}\}$,

$$(B_{\gamma}^{\tau}(\mu, \mathbf{w}))_{mn} := b_{\gamma}^{\tau}(\mu; \mathbf{w}; \varphi_n(\mu), \varphi_m(\mu)), \quad (\Theta_{\gamma}^{\tau}(\mu, \mathbf{w}))_m := \theta_{\gamma}^{\tau}(\mu; \mathbf{w}; \varphi_m(\mu)). \quad (46)$$

With this notation, we have the following decompositions:

$$B_{\gamma}^{n\tau}(\mu, \mathbf{u}_k(\mu)) := B_{\gamma}^{n}(\mu, \mathbf{u}_k(\mu)) + B_{\gamma}^{\tau}(\mu, \mathbf{u}_k(\mu)), \quad (47a)$$

$$R_{\gamma}^{n\tau}(\mu, \mathbf{u}_k(\mu)) := A_{\gamma}^{n\tau}(\mu)\mathbf{U}_k(\mu) + \Theta_{\gamma}^{n}(\mu, \mathbf{u}_k(\mu)) + \Theta_{\gamma}^{\tau}(\mu, \mathbf{u}_k(\mu)) - F(\mu). \quad (47b)$$
# 4 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.

## 4.1 Geometric mapping

We recall that the dimension $\mathcal{N}^{\text{HF}}$ of the HF (finite-dimensional) space $\mathcal{V}(\mu)$ 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 $\hat{\Omega}$. We assume that for all $\mu \in \mathcal{D}$, there is a smooth invertible geometric mapping $h(\mu) : \hat{\Omega} \to \Omega(\mu)$. We denote by $\hat{\Gamma} := \partial\hat{\Omega}$ the boundary of $\hat{\Omega}$ and we assume that it can be partitioned as $\hat{\Gamma} = \hat{\Gamma}^{\text{D}} \cup \hat{\Gamma}^{\text{N}} \cup \hat{\Gamma}^{\text{c}}$ in such a way that, for all $\mu \in \mathcal{D}$,

$$ \hat{\Gamma}^{\text{D}} := h_d^{-1}(\mu)(\Gamma^{\text{D}}(\mu)), \quad \hat{\Gamma}^{\text{N}} := h_n^{-1}(\mu)(\Gamma^{\text{N}}(\mu)), \quad \hat{\Gamma}^{\text{c}} := h_c^{-1}(\mu)(\Gamma^{\text{c}}(\mu)), \qquad (48) $$

with $h_d(\mu) := h(\mu)|_{\hat{\Gamma}^{\text{D}}}$, $h_n(\mu) := h(\mu)|_{\hat{\Gamma}^{\text{N}}}$ and $h_c(\mu) := h(\mu)|_{\hat{\Gamma}^{\text{c}}}$. Then, the mesh of $\Omega(\mu)$ is generated by generating a mesh of $\hat{\Omega}$ matching the partition of the boundary $\hat{\Gamma}$ and applying the mapping $h(\mu)$ to the mesh of the reference domain $\hat{\Omega}$. The finite element basis functions are generated from the reference basis functions by using a pullback.

## 4.2 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}(\mu_p)\}_{p \in \{1:P\}} \subset \mathbb{R}^{\mathcal{N}^{\text{HF}}}$ of HF solutions to the frictionless contact problem (28) by using a training subset $\mathcal{D}_{\text{train}} := \{\mu_p\}_{p \in \{1:P\}} \subset \mathcal{D}$ of cardinality $P \in \mathbb{N}^*$. Using the Proper Orthogonal Decomposition (POD) [19, 25] based on the canonical inner product of $H^1(\hat{\Omega}; \mathbb{R}^d)$ and the geometric mapping $h(\mu)$, one can construct an orthonormal family $\{\Xi_n\}_{n \in \{1:N\}} \subset \mathbb{R}^{\mathcal{N}^{\text{HF}}}$ of $N \in \mathbb{N}^*$ ($N \le P$) vectors. Let us denote by $V_N$ the reduced space generated by the family $\{\Xi_n\}_{n \in \{1:N\}}$, i.e.,

$$ V_N := \operatorname{Span}(\{\Xi_n\}_{n \in \{1:N\}}) \subset \mathbb{R}^{\mathcal{N}^{\text{HF}}}. \qquad (49) $$

In algebraic form, the RB formulation of the sequence of HF problems (28) reads as follows: For all $k \ge 0$, find $\Delta U_{N,k}(\mu) \in \mathbb{R}^N$ such that

$$ A_{\gamma,N}^{n}(\mu) \Delta U_{N,k}(\mu) + B_{\gamma,N}^{n}(\mu, k) \Delta U_{N,k}(\mu) = -R_{\gamma,N}^{n}(\mu, k), \qquad (50) $$

where

$$ A_{\gamma,N}^{n}(\mu) := Z^{\top} A_{\gamma}^{n}(\mu) Z \in \mathbb{R}^{N \times N}, \qquad (51a) $$

$$ B_{\gamma,N}^{n}(\mu, k) := Z^{\top} B_{\gamma}^{n}(\mu, u_{N,k}(\mu)) Z \in \mathbb{R}^{N \times N}, \quad R_{\gamma,N}^{n}(\mu, k) := Z^{\top} R_{\gamma}^{n}(\mu, u_{N,k}(\mu)) \in \mathbb{R}^{N}, \qquad (51b) $$

with $Z := [\Xi_1 \cdots \Xi_N] \in \mathbb{R}^{\mathcal{N}^{\text{HF}} \times N}$. 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:
# 4 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.

## 4.1 Geometric mapping

We recall that the dimension $\mathcal{N}^{\text{HF}}$ of the HF (finite-dimensional) space $\mathcal{V}(\mu)$ 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 $\hat{\Omega}$. We assume that for all $\mu \in \mathcal{D}$, there is a smooth invertible geometric mapping $h(\mu) : \hat{\Omega} \to \Omega(\mu)$. We denote by $\hat{\Gamma} := \partial\hat{\Omega}$ the boundary of $\hat{\Omega}$ and we assume that it can be partitioned as $\hat{\Gamma} = \hat{\Gamma}^{\text{D}} \cup \hat{\Gamma}^{\text{N}} \cup \hat{\Gamma}^{\text{c}}$ in such a way that, for all $\mu \in \mathcal{D}$,

$$ \hat{\Gamma}^{\text{D}} := h_d^{-1}(\mu)(\Gamma^{\text{D}}(\mu)), \quad \hat{\Gamma}^{\text{N}} := h_n^{-1}(\mu)(\Gamma^{\text{N}}(\mu)), \quad \hat{\Gamma}^{\text{c}} := h_c^{-1}(\mu)(\Gamma^{\text{c}}(\mu)), \qquad (48) $$

with $h_d(\mu) := h(\mu)|_{\hat{\Gamma}^{\text{D}}}$, $h_n(\mu) := h(\mu)|_{\hat{\Gamma}^{\text{N}}}$ and $h_c(\mu) := h(\mu)|_{\hat{\Gamma}^{\text{c}}}$. Then, the mesh of $\Omega(\mu)$ is generated by generating a mesh of $\hat{\Omega}$ matching the partition of the boundary $\hat{\Gamma}$ and applying the mapping $h(\mu)$ to the mesh of the reference domain $\hat{\Omega}$. The finite element basis functions are generated from the reference basis functions by using a pullback.

## 4.2 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}(\mu_p)\}_{p \in \{1:P\}} \subset \mathbb{R}^{\mathcal{N}^{\text{HF}}}$ of HF solutions to the frictionless contact problem (28) by using a training subset $\mathcal{D}_{\text{train}} := \{\mu_p\}_{p \in \{1:P\}} \subset \mathcal{D}$ of cardinality $P \in \mathbb{N}^*$. Using the Proper Orthogonal Decomposition (POD) [19, 25] based on the canonical inner product of $H^1(\hat{\Omega}; \mathbb{R}^d)$ and the geometric mapping $h(\mu)$, one can construct an orthonormal family $\{\Xi_n\}_{n \in \{1:N\}} \subset \mathbb{R}^{\mathcal{N}^{\text{HF}}}$ of $N \in \mathbb{N}^*$ ($N \le P$) vectors. Let us denote by $V_N$ the reduced space generated by the family $\{\Xi_n\}_{n \in \{1:N\}}$, i.e.,

$$ V_N := \text{Span}(\{\Xi_n\}_{n \in \{1:N\}}) \subset \mathbb{R}^{\mathcal{N}^{\text{HF}}}. \qquad (49) $$

In algebraic form, the RB formulation of the sequence of HF problems (28) reads as follows: For all $k \ge 0$, find $\Delta U_{N,k}(\mu) \in \mathbb{R}^N$ such that

$$ A_{\gamma,N}^{n}(\mu) \Delta U_{N,k}(\mu) + B_{\gamma,N}^{n}(\mu, k) \Delta U_{N,k}(\mu) = -R_{\gamma,N}^{n}(\mu, k), \qquad (50) $$

where

$$ A_{\gamma,N}^{n}(\mu) := Z^{\top} A_{\gamma}^{n}(\mu) Z \in \mathbb{R}^{N \times N}, \qquad (51a) $$

$$ B_{\gamma,N}^{n}(\mu, k) := Z^{\top} B_{\gamma}^{n}(\mu, u_{N,k}(\mu)) Z \in \mathbb{R}^{N \times N}, \quad R_{\gamma,N}^{n}(\mu, k) := Z^{\top} R_{\gamma}^{n}(\mu, u_{N,k}(\mu)) \in \mathbb{R}^{N}, \qquad (51b) $$

with $Z := [\Xi_1 \cdots \Xi_N] \in \mathbb{R}^{\mathcal{N}^{\text{HF}} \times N}$. 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 $\mathcal{D}_{\text{train}} := \{\mu_p\}_{p \in \{1:P\}} \subset \mathcal{D}$.

2. Compute the snapshots $\{U_{cv}(\mu_p)\}_{p \in \{1:P\}} \subset \mathbb{R}^{\mathcal{N}^{\text{HF}}}$ by solving (28).

3. Compute the reduced space $V_N$ by using POD on snapshots.

— Online stage: For any $\mu \in \mathcal{D} \setminus \mathcal{D}_{\text{train}}$,

1. Compute $A_{\gamma,N}^n(\mu)$ using (51a).

2. Loop on $k \ge 0$,

(a) Compute $B_{\gamma,N}^n(\mu, k)$ and $R_{\gamma,N}^n(\mu, k)$ using (51b).

(b) Solve (50).

(c) Check convergence; if not, set $k \leftarrow k + 1$ and go back to (2a).

It is crucial to derive a reduced problem that is independent of the high-fidelity dimension $\mathcal{N}^{\text{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.

*Remark 4.1* (Reconstructed functions). Let us introduce the following reconstructed functions:

$$ \xi_n(\mu) := \sum_{i \in \{1:N\}} \Xi_n^i \varphi_i(\mu) \in \mathcal{V}(\mu), \quad \forall n \in \{1:N\}. \qquad (52) $$

With this notation, solving the RB problem (50) in algebraic form leads to the following reconstructed solutions:

$$ \mathbf{u}_{N,k}(\mu) := \sum_{n \in \{1:N\}} \mathbf{u}_{N,k}^n(\mu) \xi_n(\mu) \in \mathcal{V}(\mu), \quad \mathbf{U}_{N,k}(\mu) := (\mathbf{u}_{N,k}^n(\mu))_{n \in \{1:N\}} \in \mathbb{R}^N, \qquad (53a) $$

$$ \delta \mathbf{u}_{N,k}(\mu) := \sum_{n \in \{1:N\}} \delta \mathbf{u}_{N,k}^n(\mu) \xi_n(\mu) \in \mathcal{V}(\mu), \quad \mathbf{\Delta} \mathbf{U}_{N,k}(\mu) := (\delta \mathbf{u}_{N,k}^n(\mu))_{n \in \{1:N\}} \in \mathbb{R}^N. \qquad (53b) $$

## 4.3 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, 27].

Specifically, our goal is to separate the dependence on $\mu$ and $k$ from the dependence on the indices in the large-dimensional arrays $A_{\gamma}^n(\mu)$, $B_{\gamma}^n(\mu, \mathbf{u}_{N,k}(\mu))$ and $R_{\gamma}^n(\mu, \mathbf{u}_{N,k}(\mu))$. This operation is performed during the offline stage. For this purpose, using the EIM, we obtain the following approximations:

$$ A_{\gamma}^n(\mu) \approx E^{a^n}(\mu) := \sum_{s \in \{1:S^{a^n}\}} \alpha_s^{a^n}(\mu) A_{\gamma,s}^n, \quad A_{\gamma,s}^n \in \mathbb{R}^{\mathcal{N}^{\text{HF}} \times \mathcal{N}^{\text{HF}}}, \alpha_s^{a^n}(\mu) \in \mathbb{R}, \qquad (54a) $$

$$ B_{\gamma}^n(\mu, \mathbf{u}_k(\mu)) \approx E^{b^n}(\mu, k) := \sum_{s \in \{1:S^{b^n}\}} \alpha_s^{b^n}(\mu, k) B_{\gamma,s}^n, \quad B_{\gamma,s}^n \in \mathbb{R}^{\mathcal{N}^{\text{HF}} \times \mathcal{N}^{\text{HF}}}, \alpha_s^{b^n}(\mu, k) \in \mathbb{R}, \qquad (54b) $$

$$ \Theta_{\gamma}^n(\mu, \mathbf{u}_k(\mu)) \approx E^{\theta^n}(\mu, k) := \sum_{s \in \{1:S^{\theta^n}\}} \alpha_s^{\theta^n}(\mu, k) \Theta_{\gamma,s}^n, \quad \Theta_{\gamma,s}^n \in \mathbb{R}^{\mathcal{N}^{\text{HF}}}, \alpha_s^{\theta^n}(\mu, k) \in \mathbb{R}, \qquad (54c) $$

$$ F(\mu) \approx E^f(\mu) := \sum_{s \in \{1:S^f\}} \alpha_s^f(\mu) F_s, \quad F_s \in \mathbb{R}^{\mathcal{N}^{\text{HF}}}, \alpha_s^f(\mu) \in \mathbb{R}, \qquad (54d) $$
— Offline stage

1. Select a training subset $\mathcal{D}_{\text{train}} := \{\mu_p\}_{p \in \{1:P\}} \subset \mathcal{D}$.

2. Compute the snapshots $\{U_{cv}(\mu_p)\}_{p \in \{1:P\}} \subset \mathbb{R}^{\mathcal{N}^{\text{HF}}}$ by solving (28).

3. Compute the reduced space $\mathbf{V}_N$ by using POD on snapshots.

— Online stage: For any $\mu \in \mathcal{D} \setminus \mathcal{D}_{\text{train}}$,

1. Compute $A_{\gamma,N}^n(\mu)$ using (51a).

2. Loop on $k \ge 0$,

(a) Compute $B_{\gamma,N}^n(\mu, k)$ and $R_{\gamma,N}^n(\mu, k)$ using (51b).

(b) Solve (50).

(c) Check convergence; if not, set $k \leftarrow k + 1$ and go back to (2a).

It is crucial to derive a reduced problem that is independent of the high-fidelity dimension $\mathcal{N}^{\text{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.

*Remark 4.1 (Reconstructed functions).* Let us introduce the following reconstructed functions:

$$ \xi_n(\mu) := \sum_{i \in \{1:N\}} \Xi_n^i \varphi_i(\mu) \in \mathcal{V}(\mu), \quad \forall n \in \{1:N\}. \qquad (52) $$

With this notation, solving the RB problem (50) in algebraic form leads to the following reconstructed solutions:

$$ \mathbf{u}_{N,k}(\mu) := \sum_{n \in \{1:N\}} \mathbf{u}_{N,k}^n(\mu) \xi_n(\mu) \in \mathcal{V}(\mu), \quad \mathbf{U}_{N,k}(\mu) := (\mathbf{u}_{N,k}^n(\mu))_{n \in \{1:N\}} \in \mathbb{R}^N, \qquad (53a) $$

$$ \delta \mathbf{u}_{N,k}(\mu) := \sum_{n \in \{1:N\}} \delta \mathbf{u}_{N,k}^n(\mu) \xi_n(\mu) \in \mathcal{V}(\mu), \quad \Delta \mathbf{U}_{N,k}(\mu) := (\delta \mathbf{u}_{N,k}^n(\mu))_{n \in \{1:N\}} \in \mathbb{R}^N. \qquad (53b) $$

## 4.3 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, 27].

Specifically, our goal is to separate the dependence on $\mu$ and $k$ from the dependence on the indices in the large-dimensional arrays $A_{\gamma}^n(\mu)$, $B_{\gamma}^n(\mu, \mathbf{u}_{N,k}(\mu))$ and $R_{\gamma}^n(\mu, \mathbf{u}_{N,k}(\mu))$. This operation is performed during the offline stage. For this purpose, using the EIM, we obtain the following approximations:

$$ A_{\gamma}^n(\mu) \approx E^{a^n}(\mu) := \sum_{s \in \{1:S^{a^n}\}} \alpha_s^{a^n}(\mu) A_{\gamma,s}^n, \quad A_{\gamma,s}^n \in \mathbb{R}^{\mathcal{N}^{\text{HF}} \times \mathcal{N}^{\text{HF}}}, \alpha_s^{a^n}(\mu) \in \mathbb{R}, \qquad (54a) $$

$$ B_{\gamma}^n(\mu, \mathbf{u}_k(\mu)) \approx E^{b^n}(\mu, k) := \sum_{s \in \{1:S^{b^n}\}} \alpha_s^{b^n}(\mu, k) B_{\gamma,s}^n, \quad B_{\gamma,s}^n \in \mathbb{R}^{\mathcal{N}^{\text{HF}} \times \mathcal{N}^{\text{HF}}}, \alpha_s^{b^n}(\mu, k) \in \mathbb{R}, \qquad (54b) $$

$$ \Theta_{\gamma}^n(\mu, \mathbf{u}_k(\mu)) \approx E^{\theta^n}(\mu, k) := \sum_{s \in \{1:S^{\theta^n}\}} \alpha_s^{\theta^n}(\mu, k) \Theta_{\gamma,s}^n, \quad \Theta_{\gamma,s}^n \in \mathbb{R}^{\mathcal{N}^{\text{HF}}}, \alpha_s^{\theta^n}(\mu, k) \in \mathbb{R}, \qquad (54c) $$

$$ F(\mu) \approx E^f(\mu) := \sum_{s \in \{1:S^f\}} \alpha_s^f(\mu) F_s, \quad F_s \in \mathbb{R}^{\mathcal{N}^{\text{HF}}}, \alpha_s^f(\mu) \in \mathbb{R}, \qquad (54d) $$
where the large-dimensional arrays $A_{\gamma,s}^n$, $B_{\gamma,s}^n$, $\Theta_{\gamma,s}^n$ and $F_s$ are now independent of the parameter $\mu$ and the iteration counter $k$, whereas the functions $\alpha_s^{an}$ and $\alpha_s^f$ (resp., $\alpha_s^{bn}$ and $\alpha_s^{\theta n}$) only depend on $\mu$ (resp., $(\mu, k)$). We obtain the following approximation of the residual:

$$R_{\gamma}^{n}(\mu, \mathbf{u}_{k}(\mu)) \approx E^{r^{n}}(\mu, k) := E^{a^{n}}(\mu)\mathbf{U}_{k}(\mu) + E^{\theta^{n}}(\mu, k) - E^{f}(\mu). \quad (55)$$

To build the large-dimensional arrays $A_{\gamma,s}^n$, $B_{\gamma,s}^n$, $\Theta_{\gamma,s}^n$ and $F_s$, we respectively use the training sets $\mathcal{E}_{\text{train}}^{an}$, $\mathcal{E}_{\text{train}}^{bn}$, $\mathcal{E}_{\text{train}}^{\theta n}$ and $\mathcal{E}_{\text{train}}^{f}$ defined as follows:

$$\mathcal{E}_{\text{train}}^{a^n} = \mathcal{E}_{\text{train}}^{f} := \mathcal{D}_{\text{train}}, \quad \mathcal{E}_{\text{train}}^{b^n} = \mathcal{E}_{\text{train}}^{\theta^n} := \{(\mu, k) \mid \mu \in \mathcal{D}_{\text{train}}, k \in \{1:k^{\text{cv}}(\mu)\}\}. \quad (56)$$

Notice that a different training set (possibly richer) than $\mathcal{D}_{\text{train}}$ can be used instead. We introduce the index subsets $\{(n_s^{an}, m_s^{an})\}_{s \in \{1:S^{an}\}}$, $\{(n_s^{bn}, m_s^{bn})\}_{s \in \{1:S^{bn}\}} \subset \{1:\mathcal{N}^{\text{HF}}\} \times \{1:\mathcal{N}^{\text{HF}}\}$, $\{m_s^{\theta n}\}_{s \in \{1:S^{\theta n}\}} \subset \{1:\mathcal{N}^{\text{HF}}\}$ and $\{m_s^f\}_{s \in \{1:S^f\}} \subset \{1:\mathcal{N}^{\text{HF}}\}$ of cardinality $S^{an}$, $S^{bn}$, $S^{\theta n}$, and $S^f$ respectively, corresponding to the indices selected by the EIM for the approximation of $A_{\gamma}^n$, $B_{\gamma}^n$, $\Theta_{\gamma}^n$ and $F$, respectively. Then, the functions $\alpha_s^{an}$, $\alpha_s^{bn}$, $\alpha_s^{\theta n}$ and $\alpha_s^f$ are such that

$$\forall \mu \in \mathcal{E}_{\text{train}}^{a^n}, \quad (E^{a^n}(\mu))_{n_s^{an} m_s^{an}} = (A_{\gamma}^{n}(\mu))_{n_s^{an} m_s^{an}}, \quad \forall s \in \{1:S^{a^n}\}, \quad (57a)$$

$$\forall (\mu, k) \in \mathcal{E}_{\text{train}}^{b^n}, \quad (E^{b^n}(\mu, k))_{n_s^{bn} m_s^{bn}} = (B_{\gamma}^{n}(\mu, \mathbf{u}_k(\mu)))_{n_s^{bn} m_s^{bn}}, \quad \forall s \in \{1:S^{b^n}\}, \quad (57b)$$

$$\forall (\mu, k) \in \mathcal{E}_{\text{train}}^{\theta^n}, \quad (E^{\theta^n}(\mu, k))_{m_s^{\theta^n}} = (\Theta_{\gamma}^{n}(\mu, \mathbf{u}_k(\mu)))_{m_s^{\theta^n}}, \quad \forall s \in \{1:S^{\theta^n}\}, \quad (57c)$$

$$\forall \mu \in \mathcal{E}_{\text{train}}^{f}, \quad (E^{f}(\mu))_{m_s^{f}} = (F(\mu))_{m_s^{f}}, \quad \forall s \in \{1:S^{f}\}, \quad (57d)$$

and $E^{a^n}(\mu)$, $E^{b^n}(\mu, k)$, $E^{\theta^n}(\mu, k)$ and $E^f(\mu)$ defined in (54). Notice that we use the HF solution $\mathbf{u}_k(\mu)$ instead of the RB solution $\mathbf{u}_{N,k}(\mu)$ to compute the functions $\alpha_s^{bn}$ and $\alpha_s^{\theta n}$.

In the online phase, for every new value of the parameter pair $(\mu, k) \in \mathcal{D} \times \mathbb{N}$, the functions $\alpha_s^{an}$, $\alpha_s^{bn}$, $\alpha_s^{\theta n}$ and $\alpha_s^f$ are approximated by functions $\alpha_{N,s}^{an}$, $\alpha_{N,s}^{bn}$, $\alpha_{N,s}^{\theta n}$ and $\alpha_{N,s}^f$ which solve the following linear systems:

$$Q^{a^n} \alpha_N^{a^n}(\mu) = T^{a^n}(\mu), \quad T^{a^n}(\mu) := ((A_{\gamma}^{n}(\mu))_{n_s^{a^n} m_s^{a^n}})_{s \in \{1:S^{a^n}\}} \in \mathbb{R}^{S^{a^n}}, \quad (58a)$$

$$Q^{b^n} \alpha_N^{b^n}(\mu, k) = T^{b^n}(\mu, k), \quad T^{b^n}(\mu, k) := ((B_{\gamma}^{n}(\mu, \mathbf{u}_{N,k}(\mu)))_{n_s^{b^n} m_s^{b^n}})_{s \in \{1:S^{b^n}\}} \in \mathbb{R}^{S^{b^n}}, \quad (58b)$$

$$Q^{\theta^n} \alpha_N^{\theta^n}(\mu, k) = T^{\theta^n}(\mu, k), \quad T^{\theta^n}(\mu, k) := ((\Theta_{\gamma}^{n}(\mu, \mathbf{u}_{N,k}(\mu)))_{m_s^{\theta^n}})_{s \in \{1:S^{\theta^n}\}} \in \mathbb{R}^{S^{\theta^n}}, \quad (58c)$$

$$Q^f \alpha_N^f(\mu) = T^f(\mu), \quad T^f(\mu) := ((F(\mu))_{m_s^f})_{s \in \{1:S^f\}} \in \mathbb{R}^{S^f}, \quad (58d)$$

where the vector-valued functions $\alpha_N^{an}$, $\alpha_N^{bn}$, $\alpha_N^{\theta n}$ and $\alpha_N^f$ are such that

$$\alpha_N^{a^n}(\mu) := (\alpha_{N,s}^{a^n}(\mu))_{s \in \{1:S^{a^n}\}} \in \mathbb{R}^{S^{a^n}}, \quad (59a)$$

$$\alpha_N^{b^n}(\mu, k) := (\alpha_{N,s}^{b^n}(\mu))_{s \in \{1:S^{b^n}\}} \in \mathbb{R}^{S^{b^n}}, \quad (59b)$$

$$\alpha_N^{\theta^n}(\mu, k) := (\alpha_{N,s}^{\theta^n}(\mu))_{s \in \{1:S^{\theta^n}\}} \in \mathbb{R}^{S^{\theta^n}}, \quad (59c)$$

$$\alpha_N^f(\mu) := (\alpha_{N,s}^f(\mu))_{s \in \{1:S^f\}} \in \mathbb{R}^{S^f}. \quad (59d)$$

Notice that here, in the online phase, we use the RB solutions $\mathbf{u}_{N,k}(\mu)$. The parameter-independent interpolation matrices $Q^{a^n} \in \mathbb{R}^{S^{a^n} \times S^{a^n}}$, $Q^{b^n} \in \mathbb{R}^{S^{b^n} \times S^{b^n}}$, $Q^{\theta^n} \in \mathbb{R}^{S^{\theta^n} \times S^{\theta^n}}$ and $Q^f \in \mathbb{R}^{S^f \times S^f}$ are such that

$$(Q^{a^n})_{is} := (A_{\gamma,s}^n)_{n_i^{a^n} m_i^{a^n}}, \quad \forall i, s \in \{1:S^{a^n}\}, \quad (60a)$$

$$(Q^{b^n})_{is} := (B_{\gamma,s}^n)_{n_i^{b^n} m_i^{b^n}}, \quad \forall i, s \in \{1:S^{b^n}\}, \quad (60b)$$

$$(Q^{\theta^n})_{is} := (\Theta_{\gamma,s}^n)_{m_i^{\theta^n}}, \quad \forall i, s \in \{1:S^{\theta^n}\}, \quad (60c)$$

$$(Q^f)_{is} := (F_s)_{m_i^f}, \quad \forall i, s \in \{1:S^f\}. \quad (60d)$$
where the large-dimensional arrays $A_{\gamma,s}^n, B_{\gamma,s}^n, \Theta_{\gamma,s}^n$ and $F_s$ are now independent of the parameter $\mu$ and the iteration counter $k$, whereas the functions $\alpha_s^{an}$ and $\alpha_s^f$ (resp., $\alpha_s^{bn}$ and $\alpha_s^{\theta n}$) only depend on $\mu$ (resp., $(\mu, k)$). We obtain the following approximation of the residual:

$$R_{\gamma}^{n}(\mu, \mathbf{u}_{k}(\mu)) \approx E^{r^{n}}(\mu, k) := E^{a^{n}}(\mu)\mathbf{U}_{k}(\mu) + E^{\theta^{n}}(\mu, k) - E^{f}(\mu). \quad (55)$$

To build the large-dimensional arrays $A_{\gamma,s}^n, B_{\gamma,s}^n, \Theta_{\gamma,s}^n$ and $F_s$, we respectively use the training sets $\mathcal{E}_{\text{train}}^{an}, \mathcal{E}_{\text{train}}^{bn}, \mathcal{E}_{\text{train}}^{\theta n}$ and $\mathcal{E}_{\text{train}}^{f}$ defined as follows:

$$\mathcal{E}_{\text{train}}^{an} = \mathcal{E}_{\text{train}}^{f} := \mathcal{D}_{\text{train}}, \quad \mathcal{E}_{\text{train}}^{bn} = \mathcal{E}_{\text{train}}^{\theta n} := \{(\mu, k) \mid \mu \in \mathcal{D}_{\text{train}}, k \in \{1:k^{\text{cv}}(\mu)\}\}. \quad (56)$$

Notice that a different training set (possibly richer) than $\mathcal{D}_{\text{train}}$ can be used instead. We introduce the index subsets $\{(n_s^{an}, m_s^{an})\}_{s \in \{1:S^{an}\}}, \{(n_s^{bn}, m_s^{bn})\}_{s \in \{1:S^{bn}\}} \subset \{1:\mathcal{N}^{\text{HF}}\} \times \{1:\mathcal{N}^{\text{HF}}\}, \{m_s^{\theta n}\}_{s \in \{1:S^{\theta n}\}} \subset \{1:\mathcal{N}^{\text{HF}}\} \text{ and } \{m_s^f\}_{s \in \{1:S^f\}} \subset \{1:\mathcal{N}^{\text{HF}}\} \text{ of cardinality } S^{an}, S^{bn}, S^{\theta n}, \text{ and } S^f \text{ respectively, corresponding to the indices selected by the EIM for the approximation of } A_{\gamma}^n, B_{\gamma}^n, \Theta_{\gamma}^n \text{ and } F, \text{ respectively. Then, the functions } \alpha_s^{an}, \alpha_s^{bn}, \alpha_s^{\theta n} \text{ and } \alpha_s^f \text{ are such that}$

$$\forall \mu \in \mathcal{E}_{\text{train}}^{an}, \quad (E^{a^{n}}(\mu))_{n_s^{an} m_s^{an}} = (A_{\gamma}^{n}(\mu))_{n_s^{an} m_s^{an}}, \quad \forall s \in \{1:S^{a^{n}}\}, \quad (57a)$$

$$\forall (\mu, k) \in \mathcal{E}_{\text{train}}^{bn}, \quad (E^{b^{n}}(\mu, k))_{n_s^{bn} m_s^{bn}} = (B_{\gamma}^{n}(\mu, \mathbf{u}_k(\mu)))_{n_s^{bn} m_s^{bn}}, \quad \forall s \in \{1:S^{b^{n}}\}, \quad (57b)$$

$$\forall (\mu, k) \in \mathcal{E}_{\text{train}}^{\theta n}, \quad (E^{\theta n}(\mu, k))_{m_s^{\theta n}} = (\Theta_{\gamma}^{n}(\mu, \mathbf{u}_k(\mu)))_{m_s^{\theta n}}, \quad \forall s \in \{1:S^{\theta n}\}, \quad (57c)$$

$$\forall \mu \in \mathcal{E}_{\text{train}}^{f}, \quad (E^{f}(\mu))_{m_s^{f}} = (F(\mu))_{m_s^{f}}, \quad \forall s \in \{1:S^{f}\}, \quad (57d)$$

and $E^{a^n}(\mu), E^{b^n}(\mu, k), E^{\theta^n}(\mu, k)$ and $E^f(\mu)$ defined in (54). Notice that we use the HF solution $\mathbf{u}_k(\mu)$ instead of the RB solution $\mathbf{u}_{N,k}(\mu)$ to compute the functions $\alpha_s^{bn}$ and $\alpha_s^{\theta n}$.

In the online phase, for every new value of the parameter pair $(\mu, k) \in \mathcal{D} \times \mathbb{N}$, the functions $\alpha_s^{an}$, $\alpha_s^{bn}$, $\alpha_s^{\theta n}$ and $\alpha_s^f$ are approximated by functions $\alpha_{N,s}^{an}$, $\alpha_{N,s}^{bn}$, $\alpha_{N,s}^{\theta n}$ and $\alpha_{N,s}^f$ which solve the following linear systems:

$$Q^{a^n} \alpha_N^{a^n}(\mu) = T^{a^n}(\mu), \quad T^{a^n}(\mu) := ((A_{\gamma}^{n}(\mu))_{n_s^{a^n} m_s^{a^n}})_{s \in \{1:S^{a^n}\}} \in \mathbb{R}^{S^{a^n}}, \quad (58a)$$

$$Q^{b^n} \alpha_N^{b^n}(\mu, k) = T^{b^n}(\mu, k), \quad T^{b^n}(\mu, k) := ((B_{\gamma}^{n}(\mu, \mathbf{u}_{N,k}(\mu)))_{n_s^{b^n} m_s^{b^n}})_{s \in \{1:S^{b^n}\}} \in \mathbb{R}^{S^{b^n}}, \quad (58b)$$

$$Q^{\theta^n} \alpha_N^{\theta^n}(\mu, k) = T^{\theta^n}(\mu, k), \quad T^{\theta^n}(\mu, k) := ((\Theta_{\gamma}^{n}(\mu, \mathbf{u}_{N,k}(\mu)))_{m_s^{\theta^n}})_{s \in \{1:S^{\theta^n}\}} \in \mathbb{R}^{S^{\theta^n}}, \quad (58c)$$

$$Q^f \alpha_N^f(\mu) = T^f(\mu), \quad T^f(\mu) := ((F(\mu))_{m_s^f})_{s \in \{1:S^f\}} \in \mathbb{R}^{S^f}, \quad (58d)$$

where the vector-valued functions $\alpha_N^{an}, \alpha_N^{bn}, \alpha_N^{\theta n}$ and $\alpha_N^f$ are such that

$$\alpha_N^{a^n}(\mu) := (\alpha_{N,s}^{a^n}(\mu))_{s \in \{1:S^{a^n}\}} \in \mathbb{R}^{S^{a^n}}, \quad (59a)$$

$$\alpha_N^{bn}(\mu, k) := (\alpha_{N,s}^{bn}(\mu))_{s \in \{1:S^{bn}\}} \in \mathbb{R}^{S^{bn}}, \quad (59b)$$

$$\alpha_N^{\theta^n}(\mu, k) := (\alpha_{N,s}^{\theta^n}(\mu))_{s \in \{1:S^{\theta^n}\}} \in \mathbb{R}^{S^{\theta^n}}, \quad (59c)$$

$$\alpha_N^f(\mu) := (\alpha_{N,s}^f(\mu))_{s \in \{1:S^f\}} \in \mathbb{R}^{S^f}. \quad (59d)$$

Notice that here, in the online phase, we use the RB solutions $\mathbf{u}_{N,k}(\mu)$. The parameter-independent interpolation matrices $Q^{a^n} \in \mathbb{R}^{S^{a^n} \times S^{a^n}}, Q^{bn} \in \mathbb{R}^{S^{bn} \times S^{bn}}, Q^{\theta n} \in \mathbb{R}^{S^{\theta n} \times S^{\theta n}}$ and $Q^f \in \mathbb{R}^{S^f \times S^f}$ are such that

$$(Q^{a^n})_{is} := (A_{\gamma,s}^n)_{n_i^{a^n} m_i^{a^n}}, \quad \forall i, s \in \{1:S^{a^n}\}, \quad (60a)$$

$$(Q^{b^n})_{is} := (B_{\gamma,s}^n)_{n_i^{b^n} m_i^{b^n}}, \quad \forall i, s \in \{1:S^{b^n}\}, \quad (60b)$$

$$(Q^{\theta^n})_{is} := (\Theta_{\gamma,s}^n)_{m_i^{\theta^n}}, \quad \forall i, s \in \{1:S^{\theta^n}\}, \quad (60c)$$

$$(Q^f)_{is} := (F_s)_{m_i^f}, \quad \forall i, s \in \{1:S^f\}. \quad (60d)$$
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:

$$A_{\gamma,N}^{n}(\mu) \approx E_{N}^{a^{n}}(\mu) := \sum_{s \in \{1:S^{a^{n}}\}} \alpha_{N,s}^{a^{n}}(\mu) A_{\gamma,N,s}^{n}, \quad A_{\gamma,N,s}^{n} := Z^{\top} A_{\gamma,s}^{n} Z \in \mathbb{R}^{N \times N}, \quad (61a)$$

$$B_{\gamma,N}^{n}(\mu, k) \approx E_{N}^{b^{n}}(\mu, k) := \sum_{s \in \{1:S^{b^{n}}\}} \alpha_{N,s}^{b^{n}}(\mu, k) B_{\gamma,N,s}^{n}, \quad B_{\gamma,N,s}^{n} := Z^{\top} B_{\gamma,s}^{n} Z \in \mathbb{R}^{N \times N}, \quad (61b)$$

$$\Theta_{\gamma,N}^{n}(\mu, k) \approx E_{N}^{\theta^{n}}(\mu, k) := \sum_{s \in \{1:S^{\theta^{n}}\}} \alpha_{N,s}^{\theta^{n}}(\mu, k) \Theta_{\gamma,N,s}^{n}, \quad \Theta_{\gamma,N,s}^{n} := Z^{\top} \Theta_{\gamma,s}^{n} \in \mathbb{R}^{N}, \quad (61c)$$

$$F_N(\mu) \approx E_N^f(\mu) := \sum_{s \in \{1:S^f\}} \alpha_{N,s}^f(\mu) F_{N,s}, \quad F_{N,s} := Z^\top F_s \in \mathbb{R}^N, \quad (61d)$$

which lead to an efficient offline/online decomposition since the parameter-independent arrays $A_{\gamma,N,s}^n$, $B_{\gamma,N,s}^n$, $\Theta_{\gamma,N,s}^n$ 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 $\Delta U_{N,k}(\mu)$), we consider the following sequence of problems: For all $(\mu, k) \in \mathcal{D} \times \mathbb{N}$, find $\Delta U_{N,k}(\mu) \in \mathbb{R}^N$ such that

$$E_{N}^{a^{n}}(\mu) \Delta U_{N,k}(\mu) + E_{N}^{b^{n}}(\mu, k) \Delta U_{N,k}(\mu) = -E_{N}^{r^{n}}(\mu, k), \quad (62)$$

where $E_N^{r^n}(\mu, k)$ is given by

$$E_{N}^{r^{n}}(\mu, k) := E_{N}^{a^{n}}(\mu) U_{N,k}(\mu) + E_{N}^{\theta^{n}}(\mu, k) - E_{N}^{f}(\mu). \quad (63)$$

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,\text{cv}}(\mu) \in \mathbb{R}^N$ and the associated reconstructed solution as $u_{N,\text{cv}}(\mu) \in \mathcal{V}(\mu)$.

To summarize, our RB procedure is organized as follows:

## Offline stage

1. Select a training subset $\mathcal{D}_{\text{train}} := \{\mu_p\}_{p \in \{1:P\}} \subset \mathcal{D}$.

2. Compute the HF snapshots $\{U_{cv}(\mu_p)\}_{p \in \{1:P\}} \subset \mathbb{R}^{N^{\text{HF}}}$ by solving (28) until convergence on $k$.

3. Compute the reduced space $V_N$ by using POD on the snapshots.

4. Compute the high-dimensional arrays $A_{\gamma,s}^n$, $B_{\gamma,s}^n$, $\Theta_{\gamma,s}^n$ and $F_s$ by using the EIM.

5. Invert the interpolation matrices $Q^{a^n}$, $Q^{b^n}$, $Q^{\theta^n}$ and $Q^f$.

6. Compute the small-dimensional arrays $A_{\gamma,N,s}^n$, $B_{\gamma,N,s}^n$, $\Theta_{\gamma,N,s}^n$ and $F_{N,s}$ by using (61).

## Online stage: For any $\mu \in \mathcal{D} \setminus \mathcal{D}_{\text{train}}$,

1. Evaluate $T^{a^n}(\mu)$ and $T^f(\mu)$ and compute $\alpha_N^{a^n}(\mu)$ and $\alpha_N^f(\mu)$ by using (58a) and (58d).

2. Evaluate $E_N^{a^n}(\mu)$ and $E_N^f(\mu)$ by using (61a) and (61d).

3. Loop on $k \ge 0$,

(a) Evaluate $T^{b^n}(\mu, k)$ and $T^{\theta^n}(\mu, k)$ and compute $\alpha_N^{b^n}(\mu, k)$ and $\alpha_N^{\theta^n}(\mu, k)$ by using (58b) and (58c).

(b) Evaluate $E_N^{b^n}(\mu, k)$ and $E_N^{\theta^n}(\mu, k)$ by using (61b) and (61c).
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:

$$A_{\gamma,N}^{n}(\mu) \approx E_{N}^{a^{n}}(\mu) := \sum_{s \in \{1:S^{a^{n}}\}} \alpha_{N,s}^{a^{n}}(\mu) A_{\gamma,N,s}^{n}, \quad A_{\gamma,N,s}^{n} := Z^{\top} A_{\gamma,s}^{n} Z \in \mathbb{R}^{N \times N}, \quad (61a)$$

$$B_{\gamma,N}^{n}(\mu, k) \approx E_{N}^{b^{n}}(\mu, k) := \sum_{s \in \{1:S^{b^{n}}\}} \alpha_{N,s}^{b^{n}}(\mu, k) B_{\gamma,N,s}^{n}, \quad B_{\gamma,N,s}^{n} := Z^{\top} B_{\gamma,s}^{n} Z \in \mathbb{R}^{N \times N}, \quad (61b)$$

$$\Theta_{\gamma,N}^{n}(\mu, k) \approx E_{N}^{\theta^{n}}(\mu, k) := \sum_{s \in \{1:S^{\theta^{n}}\}} \alpha_{N,s}^{\theta^{n}}(\mu, k) \Theta_{\gamma,N,s}^{n}, \quad \Theta_{\gamma,N,s}^{n} := Z^{\top} \Theta_{\gamma,s}^{n} \in \mathbb{R}^{N}, \quad (61c)$$

$$F_N(\mu) \approx E_N^f(\mu) := \sum_{s \in \{1:S^f\}} \alpha_{N,s}^f(\mu) F_{N,s}, \quad F_{N,s} := Z^\top F_s \in \mathbb{R}^N, \quad (61d)$$

which lead to an efficient offline/online decomposition since the parameter-independent arrays $A_{\gamma,N,s}^n$, $B_{\gamma,N,s}^n$, $\Theta_{\gamma,N,s}^n$ 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 $\Delta U_{N,k}(\mu)$), we consider the following sequence of problems: For all $(\mu, k) \in \mathcal{D} \times \mathbb{N}$, find $\Delta U_{N,k}(\mu) \in \mathbb{R}^N$ such that

$$E_{N}^{a^{n}}(\mu)\Delta U_{N,k}(\mu) + E_{N}^{b^{n}}(\mu, k)\Delta U_{N,k}(\mu) = -E_{N}^{r^{n}}(\mu, k), \quad (62)$$

where $E_N^{r^n}(\mu, k)$ is given by

$$E_{N}^{r^{n}}(\mu, k) := E_{N}^{a^{n}}(\mu)U_{N,k}(\mu) + E_{N}^{\theta^{n}}(\mu, k) - E_{N}^{f}(\mu). \quad (63)$$

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,\text{cv}}(\mu) \in \mathbb{R}^N$ and the associated reconstructed solution as $u_{N,\text{cv}}(\mu) \in \mathcal{V}(\mu)$.

To summarize, our RB procedure is organized as follows:

## Offline stage

1. Select a training subset $\mathcal{D}_{\text{train}} := \{\mu_p\}_{p \in \{1:P\}} \subset \mathcal{D}$.

2. Compute the HF snapshots $\{U_{cv}(\mu_p)\}_{p \in \{1:P\}} \subset \mathbb{R}^{N^{\text{HF}}}$ by solving (28) until convergence on $k$.

3. Compute the reduced space $V_N$ by using POD on the snapshots.

4. Compute the high-dimensional arrays $A_{\gamma,s}^n$, $B_{\gamma,s}^n$, $\Theta_{\gamma,s}^n$ and $F_s$ by using the EIM.

5. Invert the interpolation matrices $Q^{a^n}$, $Q^{b^n}$, $Q^{\theta^n}$ and $Q^f$.

6. Compute the small-dimensional arrays $A_{\gamma,N,s}^n$, $B_{\gamma,N,s}^n$, $\Theta_{\gamma,N,s}^n$ and $F_{N,s}$ by using (61).

## Online stage: For any $\mu \in \mathcal{D} \setminus \mathcal{D}_{\text{train}}$,

1. Evaluate $T^{a^n}(\mu)$ and $T^f(\mu)$ and compute $\alpha_N^{a^n}(\mu)$ and $\alpha_N^f(\mu)$ by using (58a) and (58d).

2. Evaluate $E_N^{a^n}(\mu)$ and $E_N^f(\mu)$ by using (61a) and (61d).

3. Loop on $k \ge 0$,

(a) Evaluate $T^{b^n}(\mu, k)$ and $T^{\theta^n}(\mu, k)$ and compute $\alpha_N^{b^n}(\mu, k)$ and $\alpha_N^{\theta^n}(\mu, k)$ by using (58b) and (58c).

(b) Evaluate $E_N^{b^n}(\mu, k)$ and $E_N^{\theta^n}(\mu, k)$ by using (61b) and (61c).
(c) Solve (62).

(d) Check convergence; if not, set $k ← k + 1$ and go back to (3a).

For the contact problem with friction, the RB formulation is obtained in exactly the same way. We simply replace the forms $a_{\gamma}^{n}(\mu; \cdot, \cdot)$, $b_{\gamma}^{n}(\mu; \cdot, \cdot, \cdot)$ and $r_{\gamma}^{n}(\mu; \cdot, \cdot)$ by the forms $a_{\gamma}^{n\tau}(\mu; \cdot, \cdot)$, $b_{\gamma}^{n\tau}(\mu; \cdot, \cdot, \cdot)$ and $r_{\gamma}^{n\tau}(\mu; \cdot, \cdot)$. Notice that for the EIM approximation, additional affine parametric decompositions are performed on the large-dimensional arrays $B_{\gamma}^{\tau}(\mu, k)$ and $\Theta_{\gamma}^{\tau}(\mu, k)$.

# 5 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 $\Omega_1(\mu) \subset \mathbb{R}^2$ of parametric radius

$$R_1(\mu) := \mu \quad \text{with } \mu \in \mathcal{D} := [0.7, 1.3](\text{m}), \qquad (64)$$

and the lower half-disk the rigid domain $\Omega_2 \subset \mathbb{R}^2$ of fixed radius $R_2 := 1\text{m}$. The initial gap between the two half-disks is equal to $g_0 > 0$. We impose a displacement of $-d$ on $\Gamma_1^{\text{top}}(\mu)$ of $\Omega_1(\mu)$ with $d \ge g_0$. The initial gap $g_0$ and the imposed displacement $d$ are, respectively, set to $g_0 := 0.001\text{m}$ and $d := 0.09\text{m}$. This latter value, which is less than 10% of the maximum value of $R_1(\mu)$ allows us to remain within the validity of the small deformation assumption. Notice that since $\Omega_2$ is rigid and fixed, we only mesh the domain $\Omega_1(\mu)$ and set $\Omega(\mu) := \Omega_1(\mu)$ to build the HF space. The material parameters are $E := 15\text{Pa}$ for the Young modulus and $\nu := 0.35$ for the Poisson coefficient. A first

<figure><img src="image_1.png" /><figcaption>Figure 1: Frictionless Hertz test case: Configuration.</figcaption></figure>

training set is typically chosen as $\mathcal{D}_{\text{train}} := \{0.7 + 0.0075i, 0 \le i \le 60\}$ (m) (altogether $P = 61$ points), and the validation set $\mathcal{D}_{\text{valid}}$ is generated by choosing 30 elements in $\mathcal{D}$ randomly with a uniform distribution. A second richer training set (with $P = 201$ points) will also be considered.

We consider the reference domain $\hat{\Omega}_1 := \Omega_1(1)$ and introduce the geometric mapping $h_1(\mu) : \hat{\Omega}_1 \to \Omega_1(\mu)$ defined as $h_1(\mu)(\boldsymbol{x}) := \mu\boldsymbol{x}$, for all $\boldsymbol{x} \in \hat{\Omega}_1$, with the origin located at the center of $\hat{\Omega}_1$. We use a mesh composed of 1956 nodes with 633 nodes on the potential contact manifold $\hat{\Gamma}_1^c$ which is the part of the half circle $\hat{\Gamma}_1^c$ of angle $\theta \in [-\frac{5\pi}{8}, -\frac{3\pi}{8}]$ with respect to the horizontal axis. For all $\mu \in \mathcal{D}$, we equip the space $\mathcal{V}(\mu)$ with the norm $\| \cdot \|_{\mathcal{V}(\mu)}$ defined as follows:

$$\|\boldsymbol{v}\|_{\mathcal{V}(\mu)} := \left( \|\boldsymbol{v}\|_{L^2(\Omega(\mu))}^2 + \hat{l}^2 \|\nabla \boldsymbol{v}\|_{L^2(\Omega(\mu))}^2 \right)^{\frac{1}{2}}, \quad \forall \boldsymbol{v} \in \mathcal{V}(\mu), \qquad (65)$$
(c) Solve (62).

(d) Check convergence; if not, set $k ← k + 1$ and go back to (3a).

For the contact problem with friction, the RB formulation is obtained in exactly the same way. We simply replace the forms $a_{\gamma}^{n}(\mu; \cdot, \cdot)$, $b_{\gamma}^{n}(\mu; \cdot, \cdot, \cdot)$ and $r_{\gamma}^{n}(\mu; \cdot, \cdot)$ by the forms $a_{\gamma}^{n\tau}(\mu; \cdot, \cdot)$, $b_{\gamma}^{n\tau}(\mu; \cdot, \cdot, \cdot)$ and $r_{\gamma}^{n\tau}(\mu; \cdot, \cdot)$. Notice that for the EIM approximation, additional affine parametric decompositions are performed on the large-dimensional arrays $B_{\gamma}^{\tau}(\mu, k)$ and $\Theta_{\gamma}^{\tau}(\mu, k)$.

# 5 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 $\Omega_1(\mu) \subset \mathbb{R}^2$ of parametric radius

$$R_1(\mu) := \mu \quad \text{with } \mu \in \mathcal{D} := [0.7, 1.3](\text{m}), \qquad (64)$$

and the lower half-disk the rigid domain $\Omega_2 \subset \mathbb{R}^2$ of fixed radius $R_2 := 1\text{m}$. The initial gap between the two half-disks is equal to $g_0 > 0$. We impose a displacement of $-d$ on $\Gamma_1^{\text{top}}(\mu)$ of $\Omega_1(\mu)$ with $d \ge g_0$. The initial gap $g_0$ and the imposed displacement $d$ are, respectively, set to $g_0 := 0.001\text{m}$ and $d := 0.09\text{m}$. This latter value, which is less than 10% of the maximum value of $R_1(\mu)$ allows us to remain within the validity of the small deformation assumption. Notice that since $\Omega_2$ is rigid and fixed, we only mesh the domain $\Omega_1(\mu)$ and set $\Omega(\mu) := \Omega_1(\mu)$ to build the HF space. The material parameters are $E := 15\text{Pa}$ for the Young modulus and $\nu := 0.35$ for the Poisson coefficient. A first

<figure><img src="image_2.png" /><figcaption>Figure 1: Frictionless Hertz test case: Configuration.</figcaption></figure>

training set is typically chosen as $\mathcal{D}_{\text{train}} := \{0.7 + 0.0075i, 0 \le i \le 60\}$ (m) (altogether $P = 61$ points), and the validation set $\mathcal{D}_{\text{valid}}$ is generated by choosing 30 elements in $\mathcal{D}$ randomly with a uniform distribution. A second richer training set (with $P = 201$ points) will also be considered.

We consider the reference domain $\hat{\Omega}_1 := \Omega_1(1)$ and introduce the geometric mapping $h_1(\mu) : \hat{\Omega}_1 \to \Omega_1(\mu)$ defined as $h_1(\mu)(\boldsymbol{x}) := \mu\boldsymbol{x}$, for all $\boldsymbol{x} \in \hat{\Omega}_1$, with the origin located at the center of $\hat{\Omega}_1$. We use a mesh composed of 1956 nodes with 633 nodes on the potential contact manifold $\hat{\Gamma}_1^c$ which is the part of the half circle $\hat{\Gamma}_1^c$ of angle $\theta \in [-\frac{5\pi}{8}, -\frac{3\pi}{8}]$ with respect to the horizontal axis. For all $\mu \in \mathcal{D}$, we equip the space $\mathcal{V}(\mu)$ with the norm $\| \cdot \|_{\mathcal{V}(\mu)}$ defined as follows:

$$\|\boldsymbol{v}\|_{\mathcal{V}(\mu)} := \left( \|\boldsymbol{v}\|_{L^2(\Omega(\mu))}^2 + \hat{l}^2 \|\nabla \boldsymbol{v}\|_{L^2(\Omega(\mu))}^2 \right)^{\frac{1}{2}}, \quad \forall \boldsymbol{v} \in \mathcal{V}(\mu), \qquad (65)$$
where the characteristic length $\hat{l} := 1$ is the radius of $\hat{\Omega}_1$ and is introduced for dimensional consistency. The HF and RB computations use the python library of the finite element software `getfem` [33].

# 5.1 Frictionless case

We first consider the frictionless Hertz contact problem.

## 5.1.1 Results using Nitsche's method

For the discretization, we use $\mathbb{P}_2$ Lagrange finite elements leading to $\mathcal{N}^{\text{HF}} := 14918$ degrees of freedom. We choose $\gamma := \frac{\gamma_0}{h}$ with $\gamma_0 := 50\mu^{\text{Lamé}}$, where $h := 2.5\text{mm}$ refers to the mesh size and $\mu^{\text{Lamé}} := \frac{E}{2(1+\nu)}\text{Pa}$ refers to the second Lamé parameter.

Figure 2 displays the deformed configurations resulting from the HF displacement fields $\mathbf{u}(\mu)$ for $\mu = 0.7\text{m}$ (left panel) and for $\mu = 1.3\text{m}$ (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 loses 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 $\hat{\Omega}_1(\mu)$ in the case of the frictional contact problems (see Section 5.2).

<figure><img src="image_2.png" /><figcaption>Figure 2: Frictionless Hertz test case, Nitsche's method (HF): $\mathbf{u}(\mu)$. Left: $\mu = 0.7\text{m}$; Right: $\mu = 1.3\text{m}$.</figcaption></figure>

Figure 3 displays the superposition of the normal stress $\sigma_{nn}(\mathbf{u}(\mu))$ and its Alart-Curnier reformulation $[P_{\gamma,g}^n(\mu; \mathbf{u}(\mu))]_{-}$ (Row 1) and the gap on the deformed configuration $u_n(\mu) - g(\mu)$ (Row 2) as a function of the abscissa along $\Gamma^c(\mu)$ for $\mu = 0.7\text{m}$ (Column 1), $\mu = 1\text{m}$ (Column 2), and $\mu = 1.3\text{m}$ (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_{\gamma,g}^n(\mu; \mathbf{u}(\mu))]_{-}$ resulting from the Alart-Curnier reformulation. We also see that the negativity condition on the gap $u_n(\mu) - g(\mu)$ is satisfied on the whole potential contact manifold $\Gamma^c(\mu)$. 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:

$$e_{\text{AC}}^{n}(\mu) := \frac{\|\sigma_{nn}(\mathbf{u}(\mu)) - [P_{\gamma,g}^{n}(\mu; \mathbf{u}(\mu))]_{-}\|_{\ell^2(\Gamma^c(\mu))}}{\|\sigma_{nn}(\mathbf{u}(\mu))\|_{\ell^2(\Gamma^c(\mu))}}, \quad (66)$$

where the discrete $\ell^2(\Gamma^c(\mu))$-norms are sampled at the mesh nodes located on $\Gamma^c(\mu)$. For the three values of the parameter $\mu \in \{0.7, 1, 1.3\}(\text{m})$, we consider three values of the mesh size $h$, namely $h = 5\text{mm}$ (coarse), $h = 2.5\text{mm}$ (medium) and $h = 1.25\text{mm}$ (fine). We notice that the relative error $e_{\text{AC}}^n(\mu)$ 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 $\sigma_{nn}(\mathbf{u}(\mu))$ (with a slight order reduction for the larger value of $\mu$). Thus, we can say that the Signorini contact conditions are globally satisfied with a good accuracy although they are not strictly enforced.
where the characteristic length $\hat{l} := 1$ is the radius of $\hat{\Omega}_1$ and is introduced for dimensional consistency. The HF and RB computations use the python library of the finite element software `getfem` [33].

# 5.1 Frictionless case

We first consider the frictionless Hertz contact problem.

## 5.1.1 Results using Nitsche's method

For the discretization, we use $\mathbb{P}_2$ Lagrange finite elements leading to $\mathcal{N}^{\text{HF}} := 14918$ degrees of freedom. We choose $\gamma := \frac{\gamma_0}{h}$ with $\gamma_0 := 50\mu^{\text{Lamé}}$, where $h := 2.5\text{mm}$ refers to the mesh size and $\mu^{\text{Lamé}} := \frac{E}{2(1+\nu)}\text{Pa}$ refers to the second Lamé parameter.

Figure 2 displays the deformed configurations resulting from the HF displacement fields $u(\mu)$ for $\mu = 0.7\text{m}$ (left panel) and for $\mu = 1.3\text{m}$ (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 loses 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 $\Omega_1(\mu)$ in the case of the frictional contact problems (see Section 5.2).

<figure><img src="image_3.png" /><figcaption>Figure 2: Frictionless Hertz test case, Nitsche's method (HF): $u(\mu)$. Left: $\mu = 0.7\text{m}$; Right: $\mu = 1.3\text{m}$.</figcaption></figure>

Figure 3 displays the superposition of the normal stress $\sigma_{nn}(u(\mu))$ and its Alart-Curnier reformulation $[P_{\gamma,g}^n(\mu; u(\mu))]_-$ (Row 1) and the gap on the deformed configuration $u_n(\mu) - g(\mu)$ (Row 2) as a function of the abscissa along $\Gamma^c(\mu)$ for $\mu = 0.7\text{m}$ (Column 1), $\mu = 1\text{m}$ (Column 2), and $\mu = 1.3\text{m}$ (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_{\gamma,g}^n(\mu; u(\mu))]_-$ resulting from the Alart-Curnier reformulation. We also see that the negativity condition on the gap $u_n(\mu) - g(\mu)$ is satisfied on the whole potential contact manifold $\Gamma^c(\mu)$. 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:

$$e_{\text{AC}}^{n}(\mu) := \frac{\|\sigma_{nn}(\mathbf{u}(\mu)) - [P_{\gamma,g}^{n}(\mu; \mathbf{u}(\mu))]_{-} \|_{\ell^2(\Gamma^c(\mu))}}{\|\sigma_{nn}(\mathbf{u}(\mu))\|_{\ell^2(\Gamma^c(\mu))}}, \quad (66)$$

where the discrete $\ell^2(\Gamma^c(\mu))$-norms are sampled at the mesh nodes located on $\Gamma^c(\mu)$. For the three values of the parameter $\mu \in \{0.7, 1, 1.3\}(\text{m})$, we consider three values of the mesh size $h$, namely $h = 5\text{mm}$ (coarse), $h = 2.5\text{mm}$ (medium) and $h = 1.25\text{mm}$ (fine). We notice that the relative error $e_{\text{AC}}^n(\mu)$ 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 $\sigma_{nn}(u(\mu))$ (with a slight order reduction for the larger value of $\mu$). Thus, we can say that the Signorini contact conditions are globally satisfied with a good accuracy although they are not strictly enforced.
<figure><img src="image_3.png" /><figcaption>Figure 3: Frictionless Hertz test case, Nitsche’s method (HF). Row 1: $\sigma_{nn}(\mathbf{u}(\mu))$ and $[P_{\gamma,g}^n(\mu; \mathbf{u}(\mu))]$. Row 2: $u_n(\mu) - g(\mu)$. Column 1: $\mu = 0.7\text{m}$. Column 2: $\mu = 1.0\text{m}$. Column 3: $\mu = 1.3\text{m}$.</figcaption></figure>

<table><caption>Table 1: Frictionless Hertz test case, Nitsche’s method (HF): Relative error $e_{AC}^n(\mu)$ for $\mu \in \{0.7, 1.0, 1.3\}(\text{m})$ and the mesh sizes $h \in \{5, 2.5, 1.25\}(\text{mm})$.</caption><thead><tr><th>μ(m)</th><th colspan="3">0.7</th><th colspan="3">1</th><th colspan="3">1.3</th></tr></thead><tbody><tr><th>h(mm)</th><td>5</td><td>2.5</td><td>1.25</td><td>5</td><td>2.5</td><td>1.25</td><td>5</td><td>2.5</td><td>1.25</td></tr><tr><th>e<sup>n</sup><sub>AC</sub>(%)</th><td>1</td><td>0.52</td><td>0.28</td><td>1.45</td><td>0.72</td><td>0.33</td><td>1.17</td><td>1.1</td><td>0.43</td></tr></tbody></table>

Let us consider the relative POD projection error defined as follows:

$$e_{\text{POD}}(N) := \frac{\left( \sum_{p \in \{1:P\}} \| (\mathbb{I}_{\mathbb{R}^{\text{HF}}} - \Pi_{V_N})(\mathbf{U}_{\text{cv}}(\mu_p)) \|_{\mathcal{W}(\mu_p)}^2 \right)^{\frac{1}{2}}}{\left( \sum_{p \in \{1:P\}} \| \mathbf{U}_{\text{cv}}(\mu_p) \|_{\mathcal{W}(\mu_p)}^2 \right)^{\frac{1}{2}}}, \quad (67)$$

where $\Pi_{V_N}$ denotes the orthogonal projection onto $V_N \subset \mathbb{R}^{\mathcal{N}^{\text{HF}}}$ and $\mathcal{W}(\mu_p)$ the Gram matrix of the inner product associated with $H^1(\Omega(\mu_p); \mathbb{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_{\text{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^{-5}$ and $10^{-7}$) for the first (resp. second) training set.

For the EIM approximation, we use the training set $\mathcal{D}_{\text{train}}$ of cardinality 61, and the training sets $\mathcal{E}_{\text{train}}^{b^n}$ and $\mathcal{E}_{\text{train}}^{\theta^n}$ introduced in (56) are then of cardinality 897. We fix a tolerance $\delta_{\text{EIM}} := 10^{-6}$ to bound the errors resulting from (54). With this choice, we obtain $S^{b^n} = 619 \ll \mathcal{N}^{\text{HF}} \times \mathcal{N}^{\text{HF}}$ and $S^{\theta^n} = 281 \ll \mathcal{N}^{\text{HF}}$. Notice that in the present test case, we do not need to perform an EIM decomposition on $A_{\gamma}^n(\mu)$ since this matrix is already linearly dependent on $\mu$; moreover, $F(\mu)$ vanishes since we only use an imposed displacement for the load. Figure 5 shows the relative EIM interpolation
<figure><img src="image_4.png" /><figcaption>Figure 3: Frictionless Hertz test case, Nitsche’s method (HF). Row 1: $\sigma_{nn}(u(\mu))$ and $[P_{\gamma,g}^n(\mu; u(\mu))]$. Row 2: $u_n(\mu) - g(\mu)$. Column 1: $\mu = 0.7\text{m}$. Column 2: $\mu = 1.0\text{m}$. Column 3: $\mu = 1.3\text{m}$.</figcaption></figure>

<table><caption>Table 1: Frictionless Hertz test case, Nitsche’s method (HF): Relative error $e_{AC}^n(\mu)$ for $\mu \in \{0.7, 1.0, 1.3\}(\text{m})$ and the mesh sizes $h \in \{5, 2.5, 1.25\}(\text{mm})$.</caption><thead><tr><th>μ(m)</th><th colspan="3">0.7</th><th colspan="3">1</th><th colspan="3">1.3</th></tr></thead><tbody><tr><th>h(mm)</th><td>5</td><td>2.5</td><td>1.25</td><td>5</td><td>2.5</td><td>1.25</td><td>5</td><td>2.5</td><td>1.25</td></tr><tr><th>e<sup>n</sup><sub>AC</sub>(%)</th><td>1</td><td>0.52</td><td>0.28</td><td>1.45</td><td>0.72</td><td>0.33</td><td>1.17</td><td>1.1</td><td>0.43</td></tr></tbody></table>

Let us consider the relative POD projection error defined as follows:

$$e_{\text{POD}}(N) := \frac{\left( \sum_{p \in \{1:P\}} \| (\mathbb{I}_{\mathbb{R}^{\text{HF}}} - \Pi_{V_N})(\mathbf{U}_{\text{cv}}(\mu_p)) \|_{\mathcal{W}(\mu_p)}^2 \right)^{\frac{1}{2}}}{\left( \sum_{p \in \{1:P\}} \| \mathbf{U}_{\text{cv}}(\mu_p) \|_{\mathcal{W}(\mu_p)}^2 \right)^{\frac{1}{2}}}, \quad (67)$$

where $\Pi_{V_N}$ denotes the orthogonal projection onto $V_N \subset \mathbb{R}^{\mathcal{N}^{\text{HF}}}$ and $\mathcal{W}(\mu_p)$ the Gram matrix of the inner product associated with $H^1(\Omega(\mu_p); \mathbb{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_{\text{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^{-5}$ and $10^{-7}$) for the first (resp. second) training set.

For the EIM approximation, we use the training set $\mathcal{D}_{\text{train}}$ of cardinality 61, and the training sets $\mathcal{E}_{\text{train}}^{b^n}$ and $\mathcal{E}_{\text{train}}^{\theta^n}$ introduced in (56) are then of cardinality 897. We fix a tolerance $\delta_{\text{EIM}} := 10^{-6}$ to bound the errors resulting from (54). With this choice, we obtain $S^{b^n} = 619 \ll \mathcal{N}^{\text{HF}} \times \mathcal{N}^{\text{HF}}$ and $S^{\theta^n} = 281 \ll \mathcal{N}^{\text{HF}}$. Notice that in the present test case, we do not need to perform an EIM decomposition on $A_{\gamma}^n(\mu)$ since this matrix is already linearly dependent on $\mu$; moreover, $F(\mu)$ vanishes since we only use an imposed displacement for the load. Figure 5 shows the relative EIM interpolation
<figure><img src="image_4.png" /><figcaption>Figure 4: Frictionless Hertz test case, Nitsche’s method: Relative POD projection error $e_{\text{POD}}(N)$ as function of $N$, for $|\mathcal{D}_{\text{train}}| = 61$ and $|\mathcal{D}_{\text{train}}| = 201$.</figcaption></figure>

errors for the tangent matrix $B_{\gamma}^{n}(\mu, \mathbf{u}_{k}(\mu))$ (left panel) and the residual vector $\Theta_{\gamma}^{n}(\mu, \mathbf{u}_{k}(\mu))$ (right panel) as a function of the rank $S^{b^n}$ or $S^{\theta^n}$, i.e., we plot

$$e_{\text{EIM}}^{b^n}(S^{b^n}, \mathcal{D}_*) := \frac{\max_{\mu \in \mathcal{D}_*} \max_{k \in \{1:k^{\text{cv}}(\mu)\}} \|B_{\gamma}^{n}(\mu, \mathbf{u}_k(\mu)) - E^{b^n}(\mu, k)\|_{\ell^{\infty}(ij)} \max_{\mu \in \mathcal{D}_*} \max_{k \in \{1:k^{\text{cv}}(\mu)\}} \|B_{\gamma}^{n}(\mu, \mathbf{u}_k(\mu))\|_{\ell^{\infty}(ij)}}{\max_{\mu \in \mathcal{D}_*} \max_{k \in \{1:k^{\text{cv}}(\mu)\}} \|\Theta_{\gamma}^{n}(\mu, \mathbf{u}_k(\mu)) - E^{\theta^n}(\mu, k)\|_{\ell^{\infty}(j)} \max_{\mu \in \mathcal{D}_*} \max_{k \in \{1:k^{\text{cv}}(\mu)\}} \|\Theta_{\gamma}^{n}(\mu, \mathbf{u}_k(\mu))\|_{\ell^{\infty}(j)}}, \quad (68a)$$

$$e_{\text{EIM}}^{\theta^n}(S^{\theta^n}, \mathcal{D}_*) := \frac{\max_{\mu \in \mathcal{D}_*} \max_{k \in \{1:k^{\text{cv}}(\mu)\}} \|\Theta_{\gamma}^{n}(\mu, \mathbf{u}_k(\mu)) - E^{\theta^n}(\mu, k)\|_{\ell^{\infty}(j)} \max_{\mu \in \mathcal{D}_*} \max_{k \in \{1:k^{\text{cv}}(\mu)\}} \|\Theta_{\gamma}^{n}(\mu, \mathbf{u}_k(\mu))\|_{\ell^{\infty}(j)}}{\max_{\mu \in \mathcal{D}_*} \max_{k \in \{1:k^{\text{cv}}(\mu)\}} \|\Theta_{\gamma}^{n}(\mu, \mathbf{u}_k(\mu))\|_{\ell^{\infty}(j)}}, \quad (68b)$$

where for a generic matrix (resp. vector) $M \in \mathbb{R}^{\mathcal{N}^{\text{HF}} \times \mathcal{N}^{\text{HF}}}$ (resp. $V \in \mathbb{R}^{\mathcal{N}^{\text{HF}}}$),

$$\|M\|_{\ell^{\infty}(ij)} := \max_{(i,j) \in \{1:\mathcal{N}^{\text{HF}}\} \times \{1:\mathcal{N}^{\text{HF}}\}} |M_{ij}|, \quad (69a)$$

$$\|V\|_{\ell^{\infty}(j)} := \max_{j \in \{1:\mathcal{N}^{\text{HF}}\}} |V_j|, \quad (69b)$$

and with $\mathcal{D}_*$ either equal to $\mathcal{D}_{\text{train}}$ or to $\mathcal{D}_{\text{valid}}$. Let us first consider the training set $\mathcal{D}_{\text{train}}$. We observe that both errors decrease fast enough to allow accurate approximations. For the tangent matrix $B_{\gamma}^{n}(\mu, \mathbf{u}_{k}(\mu))$, we notice a quasi-uniform decrease of the relative error $e_{\text{EIM}}^{b^n}(S^{b^n}, \mathcal{D}_{\text{train}})$ with an acceleration at error levels between $10^{-3}$ and $10^{-6}$. For the residual vector $\Theta_{\gamma}^{n}(\mu, \mathbf{u}_{k}(\mu))$, we observe a fast decrease of the relative error $e_{\text{EIM}}^{\theta^n}(S^{\theta^n}, \mathcal{D}_{\text{train}})$ for ranks $S^{\theta^n}$ between 1 and 40 yielding errors between 1 down to $10^{-2}$, and then a significant drop of the error at about $S^{\theta^n} = 40$ before a slower decrease occurs at error levels between $10^{-3}$ and $10^{-6}$. Considering the validation set $\mathcal{D}_{\text{valid}}$, we additionally plot the relative EIM approximation errors $e_{\text{EIM}}^{b^n, \text{cv}}(S^{b^n})$ and $e_{\text{EIM}}^{\theta^n, \text{cv}}(S^{\theta^n})$ defined as

$$e_{\text{EIM}}^{b^n, \text{cv}}(S^{b^n}) := \frac{\max_{\mu \in \mathcal{D}_{\text{valid}}} \|B_{\gamma}^{n}(\mu, \mathbf{u}_{\text{cv}}(\mu)) - E^{b^n}(\mu, k^{\text{cv}}(\mu))\|_{\ell^{\infty}(ij)}}{\max_{\mu \in \mathcal{D}_{\text{valid}}} \|B_{\gamma}^{n}(\mu, \mathbf{u}_{\text{cv}}(\mu))\|_{\ell^{\infty}(ij)}}, \quad (70a)$$

$$e_{\text{EIM}}^{\theta^n, \text{cv}}(S^{\theta^n}) := \frac{\max_{\mu \in \mathcal{D}_{\text{valid}}} \|\Theta_{\gamma}^{n}(\mu, \mathbf{u}_{\text{cv}}(\mu)) - E^{\theta^n}(\mu, k^{\text{cv}}(\mu))\|_{\ell^{\infty}(j)}}{\max_{\mu \in \mathcal{D}_{\text{valid}}} \|\Theta_{\gamma}^{n}(\mu, \mathbf{u}_{\text{cv}}(\mu))\|_{\ell^{\infty}(j)}}, \quad (70b)$$

These errors correspond to the relative EIM approximation error at convergence of the iterative algorithm, i.e., when $k = k^{\text{cv}}(\mu)$. For the tangent matrix $B_{\gamma}^{n}(\mu, \mathbf{u}_{k}(\mu))$, considering first $e_{\text{EIM}}^{b^n}(S^{b^n}, \mathcal{D}_{\text{valid}})$,
<figure><img src="image_5.png" /><figcaption>Figure 4: Frictionless Hertz test case, Nitsche’s method: Relative POD projection error $e_{\text{POD}}(N)$ as function of $N$, for $|\mathcal{D}_{\text{train}}| = 61$ and $|\mathcal{D}_{\text{train}}| = 201$.</figcaption></figure>

errors for the tangent matrix $B_{\gamma}^{n}(\mu, \mathbf{u}_{k}(\mu))$ (left panel) and the residual vector $\Theta_{\gamma}^{n}(\mu, \mathbf{u}_{k}(\mu))$ (right panel) as a function of the rank $S^{b^n}$ or $S^{\theta^n}$, i.e., we plot

$$e_{\text{EIM}}^{b^n}(S^{b^n}, \mathcal{D}_*) := \frac{\max_{\mu \in \mathcal{D}_*} \max_{k \in \{1:k^{\text{cv}}(\mu)\}} \|B_{\gamma}^{n}(\mu, \mathbf{u}_k(\mu)) - E^{b^n}(\mu, k)\|_{\ell^{\infty}(ij)}}{\max_{\mu \in \mathcal{D}_*} \max_{k \in \{1:k^{\text{cv}}(\mu)\}} \|B_{\gamma}^{n}(\mu, \mathbf{u}_k(\mu))\|_{\ell^{\infty}(ij)}}, \quad (68a)$$

$$e_{\text{EIM}}^{\theta^n}(S^{\theta^n}, \mathcal{D}_*) := \frac{\max_{\mu \in \mathcal{D}_*} \max_{k \in \{1:k^{\text{cv}}(\mu)\}} \|\Theta_{\gamma}^{n}(\mu, \mathbf{u}_k(\mu)) - E^{\theta^n}(\mu, k)\|_{\ell^{\infty}(j)}}{\max_{\mu \in \mathcal{D}_*} \max_{k \in \{1:k^{\text{cv}}(\mu)\}} \|\Theta_{\gamma}^{n}(\mu, \mathbf{u}_k(\mu))\|_{\ell^{\infty}(j)}}, \quad (68b)$$

where for a generic matrix (resp. vector) $M \in \mathbb{R}^{\mathcal{N}^{\text{HF}} \times \mathcal{N}^{\text{HF}}}$ (resp. $V \in \mathbb{R}^{\mathcal{N}^{\text{HF}}}$),

$$\|M\|_{\ell^{\infty}(ij)} := \max_{(i,j) \in \{1:\mathcal{N}^{\text{HF}}\} \times \{1:\mathcal{N}^{\text{HF}}\}} |M_{ij}|, \quad (69a)$$

$$\|V\|_{\ell^{\infty}(j)} := \max_{j \in \{1:\mathcal{N}^{\text{HF}}\}} |V_j|, \quad (69b)$$

and with $\mathcal{D}_*$ either equal to $\mathcal{D}_{\text{train}}$ or to $\mathcal{D}_{\text{valid}}$. Let us first consider the training set $\mathcal{D}_{\text{train}}$. We observe that both errors decrease fast enough to allow accurate approximations. For the tangent matrix $B_{\gamma}^{n}(\mu, \mathbf{u}_{k}(\mu))$, we notice a quasi-uniform decrease of the relative error $e_{\text{EIM}}^{b^n}(S^{b^n}, \mathcal{D}_{\text{train}})$ with an acceleration at error levels between $10^{-3}$ and $10^{-6}$. For the residual vector $\Theta_{\gamma}^{n}(\mu, \mathbf{u}_{k}(\mu))$, we observe a fast decrease of the relative error $e_{\text{EIM}}^{\theta^n}(S^{\theta^n}, \mathcal{D}_{\text{train}})$ for ranks $S^{\theta^n}$ between 1 and 40 yielding errors between 1 down to $10^{-2}$, and then a significant drop of the error at about $S^{\theta^n} = 40$ before a slower decrease occurs at error levels between $10^{-3}$ and $10^{-6}$. Considering the validation set $\mathcal{D}_{\text{valid}}$, we additionally plot the relative EIM approximation errors $e_{\text{EIM}}^{b^n, \text{cv}}(S^{b^n})$ and $e_{\text{EIM}}^{\theta^n, \text{cv}}(S^{\theta^n})$ defined as

$$e_{\text{EIM}}^{b^n, \text{cv}}(S^{b^n}) := \frac{\max_{\mu \in \mathcal{D}_{\text{valid}}} \|B_{\gamma}^{n}(\mu, \mathbf{u}_{\text{cv}}(\mu)) - E^{b^n}(\mu, k^{\text{cv}}(\mu))\|_{\ell^{\infty}(ij)}}{\max_{\mu \in \mathcal{D}_{\text{valid}}} \|B_{\gamma}^{n}(\mu, \mathbf{u}_{\text{cv}}(\mu))\|_{\ell^{\infty}(ij)}}, \quad (70a)$$

$$e_{\text{EIM}}^{\theta^n, \text{cv}}(S^{\theta^n}) := \frac{\max_{\mu \in \mathcal{D}_{\text{valid}}} \|\Theta_{\gamma}^{n}(\mu, \mathbf{u}_{\text{cv}}(\mu)) - E^{\theta^n}(\mu, k^{\text{cv}}(\mu))\|_{\ell^{\infty}(j)}}{\max_{\mu \in \mathcal{D}_{\text{valid}}} \|\Theta_{\gamma}^{n}(\mu, \mathbf{u}_{\text{cv}}(\mu))\|_{\ell^{\infty}(j)}}, \quad (70b)$$

These errors correspond to the relative EIM approximation error at convergence of the iterative algorithm, i.e., when $k = k^{\text{cv}}(\mu)$. For the tangent matrix $B_{\gamma}^{n}(\mu, \mathbf{u}_{k}(\mu))$, considering first $e_{\text{EIM}}^{b^n}(S^{b^n}, \mathcal{D}_{\text{valid}})$,
we notice a quite modest decrease of the error for ranks $S^{bn}$ 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_{EIM}^{b^n, cv}(S^{bn})$, we observe instead a rather uniform decrease of the error, with values quite close to those of $e_{EIM}^{b^n}(S^{bn}, D_{train})$, before a stagnation occurs at a value of about $10^{-3}$. For the residual vector $\Theta_{\gamma}^n(\mu, u_k(\mu))$, considering first $e_{EIM}^{\theta^n}(S^{\theta^n}, D_{valid})$, we observe a stagnation of the error for ranks $S^{\theta^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 \cdot 10^{-2}$ with some stagnation phases. Considering $e_{EIM}^{\theta^n, cv}(S^{\theta^n})$, we instead observe a stagnation at about 0.5 for ranks $S^{\theta^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 $\mu \in D_{valid}$ and small values of $k$, whereas the accuracy becomes more satisfactory as $k \to k^{cv}(\mu)$. 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.

<figure><img src="image_5.png" /><figcaption>Figure 5: Frictionless Hertz test case, Nitsche's method: Relative EIM approximation errors as a function of the rank $S^{bn}$ or $S^{\theta^n}$ of the approximation. Left: $e_{EIM}^{b^n}(S^{bn}, D_{train})$, $e_{EIM}^{b^n}(S^{bn}, D_{valid})$ and $e_{EIM}^{b^n, cv}(S^{bn})$; Right: $e_{EIM}^{\theta^n}(S^{\theta^n}, D_{train})$, $e_{EIM}^{\theta^n}(S^{\theta^n}, D_{valid})$ and $e_{EIM}^{\theta^n, cv}(S^{\theta^n})$.</figcaption></figure>

We denote by $e_N^u(\mu)$ (resp. $e_N^{nn}(\mu)$) the relative RB approximation error on the displacement (resp. normal stress) defined as

$$e_N^u(\mu) := \frac{\|u_{cv}(\mu) - u_{N, cv}(\mu)\|_{V(\mu)}}{\|u_{cv}(\mu)\|_{V(\mu)}}, \qquad (71a)$$

$$e_N^{nn}(\mu) := \frac{\|\sigma_{nn}(u_{cv}(\mu)) - \sigma_{nn}(u_{N, cv}(\mu))\|_{\ell^2(\Gamma^c(\mu))}}{\|\sigma_{nn}(u_{cv}(\mu))\|_{\ell^2(\Gamma^c(\mu))}}, \qquad (71b)$$

and introduce the relative error measures $e_{N, \max}^u$ and $e_{N, \max}^{nn}$ defined as

$$e_{N, \max}^u := \max_{\mu \in D_{\text{valid}}} e_N^u(\mu), \qquad e_{N, \max}^{nn} := \max_{\mu \in D_{\text{valid}}} e_N^{nn}(\mu). \qquad (72)$$

Figure 6 displays $e_{N, \max}^u$ (left panel) and $e_{N, \max}^{nn}$ (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 $\mu$. This can be explained by the poor quality of the EIM approximation of the tangent matrix $B_{\gamma, N}^n(\mu, 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 $\mu$ larger than 1.2m. For this reason, we consider a validation set $D_{\text{valid}}$ restricted to the interval $[0.7, 1.18](\text{m})$. In Figure 6, we superpose the plain RBM approximation error (without any EIM, thus computationally inefficient)
we notice a quite modest decrease of the error for ranks $S^{bn}$ 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_{EIM}^{b^n, cv}(S^{bn})$, we observe instead a rather uniform decrease of the error, with values quite close to those of $e_{EIM}^{b^n}(S^{bn}, D_{train})$, before a stagnation occurs at a value of about $10^{-3}$. For the residual vector $\Theta_{\gamma}^n(\mu, \mathbf{u}_k(\mu))$, considering first $e_{EIM}^{\theta^n}(S^{\theta^n}, D_{valid})$, we observe a stagnation of the error for ranks $S^{\theta^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 \cdot 10^{-2}$ with some stagnation phases. Considering $e_{EIM}^{\theta^n, cv}(S^{\theta^n})$, we instead observe a stagnation at about 0.5 for ranks $S^{\theta^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 $\mu \in D_{valid}$ and small values of $k$, whereas the accuracy becomes more satisfactory as $k \to k^{cv}(\mu)$. 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.

<figure><img src="image_6.png" /><figcaption>Figure 5: Frictionless Hertz test case, Nitsche’s method: Relative EIM approximation errors as a function of the rank $S^{bn}$ or $S^{\theta^n}$ of the approximation. Left: $e_{EIM}^{b^n}(S^{bn}, D_{train})$, $e_{EIM}^{b^n}(S^{bn}, D_{valid})$ and $e_{EIM}^{b^n, cv}(S^{bn})$; Right: $e_{EIM}^{\theta^n}(S^{\theta^n}, D_{train})$, $e_{EIM}^{\theta^n}(S^{\theta^n}, D_{valid})$ and $e_{EIM}^{\theta^n, cv}(S^{\theta^n})$.</figcaption></figure>

We denote by $e_N^u(\mu)$ (resp. $e_N^{nn}(\mu)$) the relative RB approximation error on the displacement (resp. normal stress) defined as

$$e_N^u(\mu) := \frac{\|\mathbf{u}_{cv}(\mu) - \mathbf{u}_{N, cv}(\mu)\|_{\mathcal{V}(\mu)}}{\|\mathbf{u}_{cv}(\mu)\|_{\mathcal{V}(\mu)}}, \qquad (71a)$$

$$e_N^{nn}(\mu) := \frac{\|\sigma_{nn}(\mathbf{u}_{cv}(\mu)) - \sigma_{nn}(\mathbf{u}_{N, cv}(\mu))\|_{\ell^2(\Gamma^c(\mu))}}{\|\sigma_{nn}(\mathbf{u}_{cv}(\mu))\|_{\ell^2(\Gamma^c(\mu))}}, \qquad (71b)$$

and introduce the relative error measures $e_{N, \max}^u$ and $e_{N, \max}^{nn}$ defined as

$$e_{N, \max}^u := \max_{\mu \in D_{\text{valid}}} e_N^u(\mu), \quad e_{N, \max}^{nn} := \max_{\mu \in D_{\text{valid}}} e_N^{nn}(\mu). \qquad (72)$$

Figure 6 displays $e_{N, \max}^u$ (left panel) and $e_{N, \max}^{nn}$ (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 $\mu$. This can be explained by the poor quality of the EIM approximation of the tangent matrix $B_{\gamma, N}^n(\mu, 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 $\mu$ larger than 1.2m. For this reason, we consider a validation set $D_{valid}$ restricted to the interval $[0.7, 1.18](\text{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.

<figure><img src="image_6.png" /><figcaption>Figure 6: Frictionless Hertz test case, Nitsche's method: RBM approximation errors for $|D_{\text{train}}| = 61$. Left: $e_{N,\text{max}}^u$; Right: $e_{N,\text{max}}^{nn}$.</figcaption></figure>

### 5.1.2 Comparison with the mixed formulation

For the comparison, we consider the primal-dual formulation employed in [29] with $\mathbb{P}_2$ finite elements for displacement and $\mathbb{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 $\mathcal{J}(\mu; \mathbf{u}(\mu))$ (see (10)) for all $\mu \in D_{\text{train}}$. We notice that we obtain (in the eyeball norm) the same values for the two methods for all $\mu \in D_{\text{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 $\mu$ of the half-disk $\Omega(\mu)$. The complementary

<figure><img src="image_7.png" /><figcaption>Figure 7: Frictionless Hertz test case: Comparison of the HF energy $\mathcal{J}(\mu; \mathbf{u}(\mu))$ between the mixed formulation and Nitsche's method for all $\mu \in D_{\text{train}}$.</figcaption></figure>

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.
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.

<figure><img src="image_7.png" /><figcaption>Figure 6: Frictionless Hertz test case, Nitsche's method: RBM approximation errors for $|\mathcal{D}_{\text{train}}| = 61$. Left: $e_{N,\text{max}}^u$; Right: $e_{N,\text{max}}^{nn}$.</figcaption></figure>

### 5.1.2 Comparison with the mixed formulation

For the comparison, we consider the primal-dual formulation employed in [29] with $\mathbb{P}_2$ finite elements for displacement and $\mathbb{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 $\mathcal{J}(\mu; u(\mu))$ (see (10)) for all $\mu \in \mathcal{D}_{\text{train}}$. We notice that we obtain (in the eyeball norm) the same values for the two methods for all $\mu \in \mathcal{D}_{\text{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 $\mu$ of the half-disk $\Omega(\mu)$. The complementary

<figure><img src="image_8.png" /><figcaption>Figure 7: Frictionless Hertz test case: Comparison of the HF energy $\mathcal{J}(\mu; u(\mu))$ between the mixed formulation and Nitsche's method for all $\mu \in \mathcal{D}_{\text{train}}$.</figcaption></figure>

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 [29] with a tolerance $\delta_{PGA} := \inf_{\mu \in D_{train}} \frac{\beta_{HF}(\mu)}{c_{HF}(\mu)} = 0.047$ so that the stability condition established in [29, Prop 3.1] is fulfilled (the quantities $\beta_{HF}(\mu)$ and $c_{HF}(\mu)$ are defined therein). We denote by $e_{N,R}^{u}(\mu)$ (resp. $e_{N,R}^{\lambda}(\mu)$) the primal (resp. dual) relative RB approximation error defined as follows:

$$e_{N,R}^{u}(\mu) := \frac{\|\mathbf{u}(\mu) - \mathbf{u}_{N}(\mu)\|_{V(\mu)}}{\|\mathbf{u}(\mu)\|_{V(\mu)}}, \quad e_{N,R}^{\lambda}(\mu) := \frac{\|\lambda(\mu) - \lambda_{R}(\mu)\|_{\ell^{2}(\Gamma^{c}(\mu))}}{\|\lambda(\mu)\|_{\ell^{2}(\Gamma^{c}(\mu))}}, \quad (73)$$

and introduce the relative primal (resp. dual) error measure $e_{N,R,\max}^{u}$ (resp. $e_{N,R,\max}^{\lambda}$) defined as

$$e_{N,R,\max}^{u} := \max_{\mu \in D_{\text{valid}}} e_{N,R}^{u}(\mu), \quad e_{N,R,\max}^{\lambda} := \max_{\mu \in D_{\text{valid}}} e_{N,R}^{\lambda}(\mu). \quad (74)$$

Figure 8 displays the quantity $e_{N,R,\max}^{u}$ (resp. $e_{N,R,\max}^{\lambda}$) 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_{N,\max}^{u}$ (resp. $e_{N,\max}^{nn}$) 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><img src="image_8.png" /><figcaption>Figure 8: Frictionless Hertz test case, mixed method: RBM approximation errors. Right: displacement errors $e_{N,R,\max}^{u}$ and $e_{N,R,\max}^{\lambda}$; Left: normal stress errors $e_{N,R,\max}^{\lambda}$ and $e_{N,R,\max}^{nn}$.</figcaption></figure>

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.
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 [29] with a tolerance $\delta_{PGA} := \inf_{\mu \in D_{train}} \frac{\beta_{HF}(\mu)}{c_{HF}(\mu)} = 0.047$ so that the stability condition established in [29, Prop 3.1] is fulfilled (the quantities $\beta_{HF}(\mu)$ and $c_{HF}(\mu)$ are defined therein). We denote by $e_{N,R}^{u}(\mu)$ (resp. $e_{N,R}^{\lambda}(\mu)$) the primal (resp. dual) relative RB approximation error defined as follows:

$$e_{N,R}^{u}(\mu) := \frac{\|\mathbf{u}(\mu) - \mathbf{u}_{N}(\mu)\|_{\mathcal{V}(\mu)}}{\|\mathbf{u}(\mu)\|_{\mathcal{V}(\mu)}}, \quad e_{N,R}^{\lambda}(\mu) := \frac{\|\lambda(\mu) - \lambda_{R}(\mu)\|_{\ell^{2}(\Gamma^{c}(\mu))}}{\|\lambda(\mu)\|_{\ell^{2}(\Gamma^{c}(\mu))}}, \quad (73)$$

and introduce the relative primal (resp. dual) error measure $e_{N,R,\max}^{u}$ (resp. $e_{N,R,\max}^{\lambda}$) defined as

$$e_{N,R,\max}^{u} := \max_{\mu \in D_{\text{valid}}} e_{N,R}^{u}(\mu), \quad e_{N,R,\max}^{\lambda} := \max_{\mu \in D_{\text{valid}}} e_{N,R}^{\lambda}(\mu). \quad (74)$$

Figure 8 displays the quantity $e_{N,R,\max}^{u}$ (resp. $e_{N,R,\max}^{\lambda}$) 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_{N,\max}^{u}$ (resp. $e_{N,\max}^{nn}$) 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><img src="image_9.png" /><figcaption>Figure 8: Frictionless Hertz test case, mixed method: RBM approximation errors. Right: displacement errors $e_{N,R,\max}^{u}$ and $e_{N,R,\max}^{\lambda}$; Left: normal stress errors $e_{N,R,\max}^{u}$ and $e_{N,R,\max}^{nn}$.</figcaption></figure>

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.
## 5.2 Friction case

Let us now consider the Tresca frictional Hertz contact problem. We choose a threshold $s = 0.1$. The parameter $\gamma$ 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 $\sigma_{n\tau}(\mathbf{u}(\mu))$ and its Alart-Curnier reformulation $[P_{\gamma}^{\tau}(\mu; \mathbf{u}(\mu))]_s$ as a function of the abscissa along $\Gamma^c(\mu)$ for $\mu = 0.7\text{m}$ (Column 1), $\mu = 1\text{m}$ (Column 2), and $\mu = 1.3\text{m}$ (Column 3). We observe that the tangential stress matches very well with its counterpart $[P_{\gamma}^{\tau}(\mu; \mathbf{u}(\mu))]_s$ resulting from the Alart-Curnier reformulation. However, we observe some oscillations on $\sigma_{n\tau}(\mathbf{u}(\mu))$ 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

<figure><img src="image_9.png" /><figcaption>Figure 9: Frictional Hertz test case; Nitsche's method (HF). $\sigma_{n\tau}(\mathbf{u}(\mu))$ and $[P_{\gamma}^{\tau}(\mu; \mathbf{u}(\mu))]_s$. Column 1: $\mu = 0.7\text{m}$. Column 2: $\mu = 1.0\text{m}$. Column 3: $\mu = 1.3\text{m}$.</figcaption></figure>

conditions, we display in Table 2 the relative errors on the Alart-Curnier reformulation of the Tresca friction conditions defined as follows:

$$e_{AC}^{n\tau,\mathrm{T}}(\mu) := \frac{\|\sigma_{n\tau}(\mathbf{u}(\mu)) - [P_{\gamma}^{\tau}(\mu; \mathbf{u}(\mu))]_s\|_{\ell^2(\Gamma^c(\mu))}}{\|\sigma_{nn}(\mathbf{u}(\mu))\|_{\ell^2(\Gamma^c(\mu))}} \quad (75)$$

For the three values of the parameter $\mu \in \{0.7, 1, 1.3\}(\text{m})$, we consider three values of the mesh size $h$, namely $h = 5\text{mm}$ (coarse), $h = 2.5\text{mm}$ (medium) and $h = 1.25\text{mm}$ (fine). We notice that the relative error $e_{AC}^{n\tau,\mathrm{T}}(\mu)$ 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 $\frac{1}{2}$ for the approximation of $\sigma_{n\tau}(\mathbf{u}(\mu))$. 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.

<table><caption>Table 2: Frictional Hertz test case, Nitsche's method (HF): Relative errors $e_{AC}^{n\tau,\mathrm{T}}(\mu)$ for $\mu \in \{0.7, 1.0, 1.3\}(\text{m})$ and the mesh sizes $h \in \{5, 2.5, 1.25\}(\text{mm})$.</caption><thead><tr><th>μ(m)</th><th colspan="3">0.7</th><th colspan="3">1</th><th colspan="3">1.3</th></tr></thead><tbody><tr><th>h(mm)</th><td>5</td><td>2.5</td><td>1.25</td><td>5</td><td>2.5</td><td>1.25</td><td>5</td><td>2.5</td><td>1.25</td></tr><tr><th>e<sup>nτ, T</sup><sub>AC</sub>(%)</th><td>5.49</td><td>3.33</td><td>2.36</td><td>5.5</td><td>3.35</td><td>2.33</td><td>5.76</td><td>3.4</td><td>2.37</td></tr></tbody></table>

Figure 10 shows the relative projection error $e_{\mathrm{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
## 5.2 Friction case

Let us now consider the Tresca frictional Hertz contact problem. We choose a threshold $s = 0.1$. The parameter $\gamma$ 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 $\sigma_{n\tau}(\mathbf{u}(\mu))$ and its Alart-Curnier reformulation $[P_{\gamma}^{\tau}(\mu; \mathbf{u}(\mu))]_s$ as a function of the abscissa along $\Gamma^c(\mu)$ for $\mu = 0.7\text{m}$ (Column 1), $\mu = 1\text{m}$ (Column 2), and $\mu = 1.3\text{m}$ (Column 3). We observe that the tangential stress matches very well with its counterpart $[P_{\gamma}^{\tau}(\mu; \mathbf{u}(\mu))]_s$ resulting from the Alart-Curnier reformulation. However, we observe some oscillations on $\sigma_{n\tau}(\mathbf{u}(\mu))$ 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

<figure><img src="image_10.png" /><figcaption>Figure 9: Frictional Hertz test case; Nitsche's method (HF). $\sigma_{n\tau}(\mathbf{u}(\mu))$ and $[P_{\gamma}^{\tau}(\mu; \mathbf{u}(\mu))]_s$. Column 1: $\mu = 0.7\text{m}$. Column 2: $\mu = 1.0\text{m}$. Column 3: $\mu = 1.3\text{m}$.</figcaption></figure>

conditions, we display in Table 2 the relative errors on the Alart-Curnier reformulation of the Tresca friction conditions defined as follows:

$$e_{AC}^{n\tau,\text{T}}(\mu) := \frac{\|\sigma_{n\tau}(\mathbf{u}(\mu)) - [P_{\gamma}^{\tau}(\mu; \mathbf{u}(\mu))]_s\|_{\ell^2(\Gamma^c(\mu))}}{\|\sigma_{nn}(\mathbf{u}(\mu))\|_{\ell^2(\Gamma^c(\mu))}} \quad (75)$$

For the three values of the parameter $\mu \in \{0.7, 1, 1.3\}(\text{m})$, we consider three values of the mesh size $h$, namely $h = 5\text{mm}$ (coarse), $h = 2.5\text{mm}$ (medium) and $h = 1.25\text{mm}$ (fine). We notice that the relative error $e_{AC}^{n\tau,\text{T}}(\mu)$ 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 $\frac{1}{2}$ for the approximation of $\sigma_{n\tau}(\mathbf{u}(\mu))$. 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.

<table><caption>Table 2: Frictional Hertz test case, Nitsche's method (HF): Relative errors $e_{AC}^{n\tau,\text{T}}(\mu)$ for $\mu \in \{0.7, 1.0, 1.3\}(\text{m})$ and the mesh sizes $h \in \{5, 2.5, 1.25\}(\text{mm})$.</caption><thead><tr><th>μ(m)</th><th colspan="3">0.7</th><th colspan="3">1</th><th colspan="3">1.3</th></tr></thead><tbody><tr><th>h(mm)</th><td>5</td><td>2.5</td><td>1.25</td><td>5</td><td>2.5</td><td>1.25</td><td>5</td><td>2.5</td><td>1.25</td></tr><tr><th>e<sup>nτ, T</sup><sub>AC</sub>(%)</th><td>5.49</td><td>3.33</td><td>2.36</td><td>5.5</td><td>3.35</td><td>2.33</td><td>5.76</td><td>3.4</td><td>2.37</td></tr></tbody></table>

Figure 10 shows the relative projection error $e_{\text{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
<figure><img src="image_10.png" /><figcaption>Figure 10: Frictional Hertz test case, Nitsche’s method: Relative POD projection error $e_{\text{POD}}(N)$ as function of $N$, for $|\mathcal{D}_{\text{train}}| = 61$.</figcaption></figure>

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 $\mathcal{E}_{\text{train}}^{b^{n\tau}}$, $\mathcal{E}_{\text{train}}^{\theta^n}$ and $\mathcal{E}_{\text{train}}^{\theta^\tau}$ are of cardinality 898. We fix a tolerance $\delta_{\text{EIM}} := 10^{-6}$. With this choice, we obtain $S^{b^{n\tau}} = 630 \ll \mathcal{N}^{\text{HF}} \times \mathcal{N}^{\text{HF}}$, $S^{\theta^n} = 291 \ll \mathcal{N}^{\text{HF}}$ and $S^{\theta^\tau} = 3 \ll \mathcal{N}^{\text{HF}}$. For the same reason as for the frictionless case, we do not need to perform an EIM decomposition on $A_{\gamma}^{n\tau}(\mu)$ and on $F(\mu)$. Figure 11 shows the relative EIM interpolation errors for the tangent matrix $B_{\gamma}^{n\tau}(\mu, \mathbf{u}_k(\mu))$ (left panel) and the residual vector $\Theta_{\gamma}^{n}(\mu, \mathbf{u}_k(\mu))$ (right panel) as a function of the rank $S^{b^{n\tau}}$ or $S^{\theta^n}$, i.e., we plot $e_{\text{EIM}}^{\theta^n}(S^{\theta^n}, \mathcal{D}_*)$ defined in (68b) and

$$e_{\text{EIM}}^{b^{n\tau}}(S^{b^{n\tau}}, \mathcal{D}_*) := \frac{\max_{\mu \in \mathcal{D}_*} \max_{k \in \{1:k^{\text{cv}}(\mu)\}} \|B_{\gamma}^{n\tau}(\mu, \mathbf{u}_k(\mu)) - E^{b^{n\tau}}(\mu, k)\|_{\ell^\infty(ij)}}{\max_{\mu \in \mathcal{D}_*} \max_{k \in \{1:k^{\text{cv}}(\mu)\}} \|B_{\gamma}^{n\tau}(\mu, \mathbf{u}_k(\mu))\|_{\ell^\infty(ij)}}, \quad (76)$$

with $\mathcal{D}_*$ either equal to $\mathcal{D}_{\text{train}}$ or to $\mathcal{D}_{\text{valid}}$. For the validation set $\mathcal{D}_{\text{valid}}$, we additionally plot the relative EIM approximation errors

$$e_{\text{EIM}}^{b^{n\tau}, \text{cv}}(S^{b^{n\tau}}) := \frac{\max_{\mu \in \mathcal{D}_{\text{valid}}} \|B_{\gamma}^{n\tau}(\mu, \mathbf{u}_{\text{cv}}(\mu)) - E^{b^{n\tau}}(\mu, k^{\text{cv}}(\mu))\|_{\ell^\infty(ij)}}{\max_{\mu \in \mathcal{D}_{\text{valid}}} \|B_{\gamma}^{n\tau}(\mu, \mathbf{u}_{\text{cv}}(\mu))\|_{\ell^\infty(ij)}}, \quad (77)$$

and $e_{\text{EIM}}^{\theta^n, \text{cv}}(S^{\theta^n})$ defined in (70b). Notice that for Tresca friction, we perform the EIM on the tangent matrix $B_{\gamma}^{n\tau}(\mu; \mathbf{u}_k(\mu))$ instead of performing it separately on the tangent matrices $B_{\gamma}^{n}(\mu; \mathbf{u}_k(\mu))$ and $B_{\gamma}^{\tau}(\mu; \mathbf{u}_k(\mu))$. We see that the result is close to the one obtained for the frictionless case because the contribution of $B_{\gamma}^{\tau}(\mu; \mathbf{u}_k(\mu))$ is negligible. Indeed, as can be seen in Figure 9 (Row 1), the Alart-Curnier reformulation of the tangential stress is equal to $\pm s$ almost everywhere (except at 3 nodes) which results in the nullity of $\mathcal{G}_s(P_{\gamma}^{\tau}(\mu; \mathbf{u}_k(\mu)))$ almost everywhere. For the same reason, the dependence of the residual vector $\Theta_{\gamma}^{\tau}(\mu, \mathbf{u}_k(\mu))$ on $(\mu, k)$ is almost of rank one and therefore its EIM approximation is very simple ($S^{\theta^\tau} = 3$); hence, its relative EIM approximation error is not reported. Altogether, the behaviour of the residual vector $\Theta_{\gamma}^{n}(\mu, \mathbf{u}_k(\mu))$ remains similar to that observed in the frictionless case.

We denote by $e_{N}^{n\tau}(\mu)$ the relative RB approximation error on the tangential stress defined as

$$e_{N}^{n\tau}(\mu) := \frac{\|\sigma_{n\tau}(\mathbf{u}_{\text{cv}}(\mu)) - \sigma_{n\tau}(\mathbf{u}_{N,\text{cv}}(\mu))\|_{\ell^2(\Gamma^c(\mu))}}{\|\sigma_{n\tau}(\mathbf{u}_{\text{cv}}(\mu))\|_{\ell^2(\Gamma^c(\mu))}}, \quad (78)$$
<figure><img src="image_11.png" /><figcaption>Figure 10: Frictional Hertz test case, Nitsche’s method: Relative POD projection error $e_{\text{POD}}(N)$ as function of $N$, for $|\mathcal{D}_{\text{train}}| = 61$.</figcaption></figure>

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 $\mathcal{E}_{\text{train}}^{b^{n\tau}}$, $\mathcal{E}_{\text{train}}^{\theta^n}$ and $\mathcal{E}_{\text{train}}^{\theta^\tau}$ are of cardinality 898. We fix a tolerance $\delta_{\text{EIM}} := 10^{-6}$. With this choice, we obtain $S^{b^{n\tau}} = 630 \ll \mathcal{N}^{\text{HF}} \times \mathcal{N}^{\text{HF}}$, $S^{\theta^n} = 291 \ll \mathcal{N}^{\text{HF}}$ and $S^{\theta^\tau} = 3 \ll \mathcal{N}^{\text{HF}}$. For the same reason as for the frictionless case, we do not need to perform an EIM decomposition on $A_{\gamma}^{n\tau}(\mu)$ and on $F(\mu)$. Figure 11 shows the relative EIM interpolation errors for the tangent matrix $B_{\gamma}^{n\tau}(\mu, \mathbf{u}_k(\mu))$ (left panel) and the residual vector $\Theta_{\gamma}^{n}(\mu, \mathbf{u}_k(\mu))$ (right panel) as a function of the rank $S^{b^{n\tau}}$ or $S^{\theta^n}$, i.e., we plot $e_{\text{EIM}}^{\theta^n}(S^{\theta^n}, \mathcal{D}_*)$ defined in (68b) and

$$e_{\text{EIM}}^{b^{n\tau}}(S^{b^{n\tau}}, \mathcal{D}_*) := \frac{\max_{\mu \in \mathcal{D}_*} \max_{k \in \{1:k^{\text{cv}}(\mu)\}} \|B_{\gamma}^{n\tau}(\mu, \mathbf{u}_k(\mu)) - E^{b^{n\tau}}(\mu, k)\|_{\ell^\infty(ij)}}{\max_{\mu \in \mathcal{D}_*} \max_{k \in \{1:k^{\text{cv}}(\mu)\}} \|B_{\gamma}^{n\tau}(\mu, \mathbf{u}_k(\mu))\|_{\ell^\infty(ij)}}, \quad (76)$$

with $\mathcal{D}_*$ either equal to $\mathcal{D}_{\text{train}}$ or to $\mathcal{D}_{\text{valid}}$. For the validation set $\mathcal{D}_{\text{valid}}$, we additionally plot the relative EIM approximation errors

$$e_{\text{EIM}}^{b^{n\tau}, \text{cv}}(S^{b^{n\tau}}) := \frac{\max_{\mu \in \mathcal{D}_{\text{valid}}} \|B_{\gamma}^{n\tau}(\mu, \mathbf{u}_{\text{cv}}(\mu)) - E^{b^{n\tau}}(\mu, k^{\text{cv}}(\mu))\|_{\ell^\infty(ij)}}{\max_{\mu \in \mathcal{D}_{\text{valid}}} \|B_{\gamma}^{n\tau}(\mu, \mathbf{u}_{\text{cv}}(\mu))\|_{\ell^\infty(ij)}}, \quad (77)$$

and $e_{\text{EIM}}^{\theta^n, \text{cv}}(S^{\theta^n})$ defined in (70b). Notice that for Tresca friction, we perform the EIM on the tangent matrix $B_{\gamma}^{n\tau}(\mu; \mathbf{u}_k(\mu))$ instead of performing it separately on the tangent matrices $B_{\gamma}^{n}(\mu; \mathbf{u}_k(\mu))$ and $B_{\gamma}^{\tau}(\mu; \mathbf{u}_k(\mu))$. We see that the result is close to the one obtained for the frictionless case because the contribution of $B_{\gamma}^{\tau}(\mu; \mathbf{u}_k(\mu))$ is negligible. Indeed, as can be seen in Figure 9 (Row 1), the Alart-Curnier reformulation of the tangential stress is equal to $\pm s$ almost everywhere (except at 3 nodes) which results in the nullity of $\mathcal{G}_s(P_{\gamma}^{\tau}(\mu; \mathbf{u}_k(\mu)))$ almost everywhere. For the same reason, the dependence of the residual vector $\Theta_{\gamma}^{\tau}(\mu, \mathbf{u}_k(\mu))$ on $(\mu, k)$ is almost of rank one and therefore its EIM approximation is very simple ($S^{\theta^\tau} = 3$); hence, its relative EIM approximation error is not reported. Altogether, the behaviour of the residual vector $\Theta_{\gamma}^{n}(\mu, \mathbf{u}_k(\mu))$ remains similar to that observed in the frictionless case.

We denote by $e_{N}^{n\tau}(\mu)$ the relative RB approximation error on the tangential stress defined as

$$e_{N}^{n\tau}(\mu) := \frac{\|\sigma_{n\tau}(\mathbf{u}_{\text{cv}}(\mu)) - \sigma_{n\tau}(\mathbf{u}_{N,\text{cv}}(\mu))\|_{\ell^2(\Gamma^c(\mu))}}{\|\sigma_{n\tau}(\mathbf{u}_{\text{cv}}(\mu))\|_{\ell^2(\Gamma^c(\mu))}}, \quad (78)$$
<figure><img src="image_11.png" /><figcaption>Figure 11: Frictional Hertz test case, Nitsche’s method: Relative EIM approximation errors as a function of the rank $S^{b^{n\tau}}$ or $S^{\theta^n}$ of the approximation. Left: $e_{EIM}^{b^{n\tau}}(S^{b^{n\tau}}, \mathcal{D}_{train})$, $e_{EIM}^{b^{n\tau}}(S^{b^{n\tau}}, \mathcal{D}_{valid})$ and $e_{EIM}^{b^{n\tau},cv}(S^{b^{n\tau}})$; Right: $e_{EIM}^{\theta^n}(S^{\theta^n}, \mathcal{D}_{train})$, $e_{EIM}^{\theta^n}(S^{\theta^n}, \mathcal{D}_{valid})$ and $e_{EIM}^{\theta^n, cv}(S^{\theta^n})$.</figcaption></figure>

and introduce the relative error measure $e_{N,\max}^{n\tau}$ defined as

$$e_{N,\max}^{n\tau} := \max_{\mu \in \mathcal{D}_{\text{valid}}} e_N^{n\tau}(\mu). \qquad (79)$$

Figure 12 displays for Tresca’s friction the relative errors $e_{N,\max}^u$ (left panel) and $e_{N,\max}^{n\tau}$ (right panel) as a function of the number of vectors composing the reduced basis. As for the frictionless case, we only consider RB dimensions $N$ larger than 10 and the validation set $\mathcal{D}_{valid}$ is restricted to the interval $[0.7, 1.18](m)$. In Figure 12, we superpose the plain RBM approximation error (without any EIM) and the RBM-EIM approximation error. As for the frictionless case, we observe that similar errors are obtained for plain RBM and RBM-EIM. This again confirms the satisfactory quality of the EIM approximations at convergence.

<figure><img src="image_12.png" /><figcaption>Figure 12: Frictional Hertz test case, Nitsche’s method: RBM approximation errors for $|\mathcal{D}_{train}| = 61$. Left: $e_{N,\max}^u$; Right: $e_{N,\max}^{n\tau}$.</figcaption></figure>

## 5.3 Toward Coulomb friction

A more realistic model for friction is given by Coulomb conditions which read as follows:

$$\left\{ \begin{array}{ll} \|\sigma_{n\tau}(\mathbf{u}(\mu))\| \le \nu_F |\sigma_{nn}(\mathbf{u}(\mu))|, & \text{if } \mathbf{u}_{\tau}(\mu) = \mathbf{0}, \\ \sigma_{n\tau}(\mathbf{u}(\mu)) = -\nu_F |\sigma_{nn}(\mathbf{u}(\mu))| \frac{\mathbf{u}_{\tau}(\mu)}{\|\mathbf{u}_{\tau}(\mu)\|}, & \text{otherwise,} \end{array} \right. \qquad (80)$$
<figure><img src="image_12.png" /><figcaption>Figure 11: Frictional Hertz test case, Nitsche’s method: Relative EIM approximation errors as a function of the rank $S^{b^{n\tau}}$ or $S^{\theta^n}$ of the approximation. Left: $e_{EIM}^{b^{n\tau}}(S^{b^{n\tau}}, \mathcal{D}_{train})$, $e_{EIM}^{b^{n\tau}}(S^{b^{n\tau}}, \mathcal{D}_{valid})$ and $e_{EIM}^{b^{n\tau},cv}(S^{b^{n\tau}})$; Right: $e_{EIM}^{\theta^n}(S^{\theta^n}, \mathcal{D}_{train})$, $e_{EIM}^{\theta^n}(S^{\theta^n}, \mathcal{D}_{valid})$ and $e_{EIM}^{\theta^n, cv}(S^{\theta^n})$.</figcaption></figure>

and introduce the relative error measure $e_{N,\max}^{n\tau}$ defined as

$$e_{N,\max}^{n\tau} := \max_{\mu \in \mathcal{D}_{\text{valid}}} e_N^{n\tau}(\mu). \qquad (79)$$

Figure 12 displays for Tresca’s friction the relative errors $e_{N,\max}^u$ (left panel) and $e_{N,\max}^{n\tau}$ (right panel) as a function of the number of vectors composing the reduced basis. As for the frictionless case, we only consider RB dimensions $N$ larger than 10 and the validation set $\mathcal{D}_{valid}$ is restricted to the interval $[0.7, 1.18](m)$. In Figure 12, we superpose the plain RBM approximation error (without any EIM) and the RBM-EIM approximation error. As for the frictionless case, we observe that similar errors are obtained for plain RBM and RBM-EIM. This again confirms the satisfactory quality of the EIM approximations at convergence.

<figure><img src="image_13.png" /><figcaption>Figure 12: Frictional Hertz test case, Nitsche’s method: RBM approximation errors for $|\mathcal{D}_{train}| = 61$. Left: $e_{N,\max}^u$; Right: $e_{N,\max}^{n\tau}$.</figcaption></figure>

## 5.3 Toward Coulomb friction

A more realistic model for friction is given by Coulomb conditions which read as follows:

$$ \left\{ \begin{array}{ll} \|\sigma_{n\tau}(\mathbf{u}(\mu))\| \le \nu_F |\sigma_{nn}(\mathbf{u}(\mu))|, & \text{if } \mathbf{u}_{\tau}(\mu) = \mathbf{0}, \\ \sigma_{n\tau}(\mathbf{u}(\mu)) = -\nu_F |\sigma_{nn}(\mathbf{u}(\mu))| \frac{\mathbf{u}_{\tau}(\mu)}{\|\mathbf{u}_{\tau}(\mu)\|}, & \text{otherwise,} \end{array} \right. \qquad (80) $$
where $\nu_F > 0$ is a given nondimensional coefficient (which can be taken to be constant for simplicity). A rather well-established approach to solve the Coulomb frictional problem is to use a fixed-point method on the Tresca frictional problem. Specifically, one introduces the mapping $\phi(\mu; \cdot) : \mathbf{V}(\mu) \to \mathbf{V}(\mu)$ defined for all $\mathbf{v} \in \mathbf{V}(\mu)$ by requiring that $\phi(\mu; \mathbf{v})$ solves (38) with the threshold $s := s^c(\mu; \mathbf{v})$ with $s^c(\mu; \mathbf{v}) := \nu_F |[P_{\gamma,g}^n(\mu; \mathbf{v})]|$. It is shown in [10, Thm 4.3] that, under certain conditions, the mapping $\phi(\mu; \cdot)$ is contractive. Thus, the Coulomb frictional problem can be solved by means of a nested loop, where the outer iteration index, say $n$, refers to the fixed-point iteration on $\phi(\mu; \cdot)$, and the inner iteration index, say $k$, refers as in (38) to the iterative solution of the Tresca frictional problem.

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,

<figure><img src="image_13.png" /><figcaption>Figure 13: Frictional Hertz test case with Coulomb friction law, Nitsche's method (HF). $\sigma_{n\tau}(\mathbf{u}(\mu))$ and $[P_{\gamma}^{\tau}(\mu; \mathbf{u}(\mu))]_{s^c(\mu;\mathbf{u}(\mu))}$. Column 1: $\mu = 0.7\text{m}$. Column 2: $\mu = 1.0\text{m}$. Column 3: $\mu = 1.3\text{m}$.</figcaption></figure>

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 $(\mu, 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.

## References

[1] L. Baillet and T. Sassi. Mixed finite element methods for the Signorini problem with friction. *Numer. Methods Partial Differential Equations*, 22(6):1489–1508, 2006.

[2] M. Balajewicz, D. Amsallem, and C. Farhat. Projection-based model reduction for contact problems. *Internat. J. Numer. Methods Engrg.*, 106(8):644–663, 2016.

[3] M. Barrault, Y. Maday, N. C. Nguyen, and A. T. Patera. An 'empirical interpolation' method: application to efficient reduced-basis discretization of partial differential equations. *C. R. Math. Acad. Sci. Paris*, 339(9):667–672, 2004.

[4] A. Benaceur, A. Ern, and V. Ehrlacher. A reduced basis method for parametrized variational inequalities applied to contact mechanics. *Internat. J. Numer. Methods Engrg.*, 121(6):1170–1197, 2020.

[5] A. Buffa, Y. Maday, A. T. Patera, C. Prud'homme, and G. Turinici. A priori convergence of the greedy algorithm for the parametrized reduced basis method. *ESAIM Math. Model. Numer. Anal.*, 46(3):595–603, 2012.
where $\nu_F > 0$ is a given nondimensional coefficient (which can be taken to be constant for simplicity). A rather well-established approach to solve the Coulomb frictional problem is to use a fixed-point method on the Tresca frictional problem. Specifically, one introduces the mapping $\phi(\mu; \cdot) : \mathbf{V}(\mu) \to \mathbf{V}(\mu)$ defined for all $\mathbf{v} \in \mathbf{V}(\mu)$ by requiring that $\phi(\mu; \mathbf{v})$ solves (38) with the threshold $s := s^c(\mu; \mathbf{v})$ with $s^c(\mu; \mathbf{v}) := \nu_F |[P_{\gamma,g}^n(\mu; \mathbf{v})]|$. It is shown in [10, Thm 4.3] that, under certain conditions, the mapping $\phi(\mu; \cdot)$ is contractive. Thus, the Coulomb frictional problem can be solved by means of a nested loop, where the outer iteration index, say $n$, refers to the fixed-point iteration on $\phi(\mu; \cdot)$, and the inner iteration index, say $k$, refers as in (38) to the iterative solution of the Tresca frictional problem.

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,

<figure><img src="image_14.png" /><figcaption>Figure 13: Frictional Hertz test case with Coulomb friction law, Nitsche's method (HF). $\sigma_{n\tau}(\mathbf{u}(\mu))$ and $[P_{\gamma}^{\tau}(\mu; \mathbf{u}(\mu))]_{s^c(\mu;\mathbf{u}(\mu))}$. Column 1: $\mu = 0.7\text{m}$. Column 2: $\mu = 1.0\text{m}$. Column 3: $\mu = 1.3\text{m}$.</figcaption></figure>

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 $(\mu, 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.

## References

[1] L. Baillet and T. Sassi. Mixed finite element methods for the Signorini problem with friction. *Numer. Methods Partial Differential Equations*, 22(6):1489–1508, 2006.

[2] M. Balajewicz, D. Amsallem, and C. Farhat. Projection-based model reduction for contact problems. *Internat. J. Numer. Methods Engrg.*, 106(8):644–663, 2016.

[3] M. Barrault, Y. Maday, N. C. Nguyen, and A. T. Patera. An 'empirical interpolation' method: application to efficient reduced-basis discretization of partial differential equations. *C. R. Math. Acad. Sci. Paris*, 339(9):667–672, 2004.

[4] A. Benaceur, A. Ern, and V. Ehrlacher. A reduced basis method for parametrized variational inequalities applied to contact mechanics. *Internat. J. Numer. Methods Engrg.*, 121(6):1170–1197, 2020.

[5] A. Buffa, Y. Maday, A. T. Patera, C. Prud'homme, and G. Turinici. A priori convergence of the greedy algorithm for the parametrized reduced basis method. *ESAIM Math. Model. Numer. Anal.*, 46(3):595–603, 2012.
[6] F. Chouly. An adaptation of Nitsche's method to the Tresca friction problem. *J. Math. Anal. Appl.*, 411(1):329-339, 2014.

[7] F. Chouly, A. Ern, and N. Pignet. A hybrid high-order discretization combined with Nitsche's method for contact and Tresca friction in small strain elasticity. *SIAM J. Sci. Comput.*, 42(4):A2300-A2324, 2020.

[8] F. Chouly, M. Fabre, P. Hild, R. Mlika, J. Pousin, and Y. Renard. An overview of recent results on Nitsche's method for contact problems. In *Geometrically unfitted finite element methods and applications*, volume 121 of *Lect. Notes Comput. Sci. Eng.*, pages 93-141. Springer, Cham, 2017.

[9] F. Chouly and P. Hild. A Nitsche-based method for unilateral contact problems: numerical analysis. *SIAM J. Numer. Anal.*, 51(2):1295-1307, 2013.

[10] F. Chouly, P. Hild, V. Lleras, and Y. Renard. Nitsche method for contact with Coulomb friction: existence results for the static and dynamic finite element formulations. *J. Comput. Appl. Math.*, 416, 2022.

[11] F. Chouly, P. Hild, and Y. Renard. Symmetric and non-symmetric variants of Nitsche's method for contact problems in elasticity: theory and numerical experiments. *Math. Comp.*, 84(293):1089-1112, 2015.

[12] A. Curnier and P. Alart. A generalized Newton method for contact problems with friction. *J. Méc. Théor. Appl.*, 7(suppl. 1):67-82, 1988.

[13] G. Duvaut and J.-L. Lions. *Les inéquations en mécanique et en physique*. Travaux et Recherches Mathématiques, No. 21. Dunod, Paris, 1972.

[14] J. Fauque, I. Ramière, and D. Ryckelynck. Hybrid hyper-reduced modeling for contact mechanics problems. *Internat. J. Numer. Methods Engrg.*, 115(1):117-139, 2018.

[15] G. Fichera. Problemi elastostatici con vincoli unilaterali: Il problema di Signorini con ambigue condizioni al contorno. *Atti Accad. Naz. Lincei Mem. Cl. Sci. Fis. Mat. Natur. Sez. Ia (8)*, 7:91-140, 1963/64.

[16] M. Fortin and R. Glowinski. *Augmented Lagrangian methods*, volume 15 of *Studies in Mathematics and its Applications*. North-Holland Publishing Co., Amsterdam, 1983. Applications to the numerical solution of boundary value problems, Translated from the French by B. Hunt and D. C. Spicer.

[17] T. Gustafsson, R. Stenberg, and J. Videman. On Nitsche's method for elastic contact problems. *SIAM J. Sci. Comput.*, 42(2):B425-B446, 2020.

[18] T. Gustafsson and J. Videman. Stabilized finite elements for tresca friction problem. *ESAIM Math. Model. Numer. Anal.*, 56(4):1307-1326, 2022.

[19] B. Haasdonk. Convergence rates of the POD-greedy method. *ESAIM Math. Model. Numer. Anal.*, 47(3):859-873, 2013.

[20] B. Haasdonk, J. Salomon, and B. Wohlmuth. A reduced basis method for parametrized variational inequalities. *SIAM J. Numer. Anal.*, 50(5):2656-2676, 2012.

[21] J. S. Hesthaven, G. Rozza, and B. Stamm. *Certified reduced basis methods for parametrized partial differential equations*, volume 590. Springer, 2016.
[6] F. Chouly. An adaptation of Nitsche's method to the Tresca friction problem. *J. Math. Anal. Appl.*, 411(1):329-339, 2014.

[7] F. Chouly, A. Ern, and N. Pignet. A hybrid high-order discretization combined with Nitsche's method for contact and Tresca friction in small strain elasticity. *SIAM J. Sci. Comput.*, 42(4):A2300-A2324, 2020.

[8] F. Chouly, M. Fabre, P. Hild, R. Mlika, J. Pousin, and Y. Renard. An overview of recent results on Nitsche's method for contact problems. In *Geometrically unfitted finite element methods and applications*, volume 121 of *Lect. Notes Comput. Sci. Eng.*, pages 93-141. Springer, Cham, 2017.

[9] F. Chouly and P. Hild. A Nitsche-based method for unilateral contact problems: numerical analysis. *SIAM J. Numer. Anal.*, 51(2):1295-1307, 2013.

[10] F. Chouly, P. Hild, V. Lleras, and Y. Renard. Nitsche method for contact with Coulomb friction: existence results for the static and dynamic finite element formulations. *J. Comput. Appl. Math.*, 416, 2022.

[11] F. Chouly, P. Hild, and Y. Renard. Symmetric and non-symmetric variants of Nitsche's method for contact problems in elasticity: theory and numerical experiments. *Math. Comp.*, 84(293):1089-1112, 2015.

[12] A. Curnier and P. Alart. A generalized Newton method for contact problems with friction. *J. Méc. Théor. Appl.*, 7(suppl. 1):67-82, 1988.

[13] G. Duvaut and J.-L. Lions. *Les inéquations en mécanique et en physique*. Travaux et Recherches Mathématiques, No. 21. Dunod, Paris, 1972.

[14] J. Fauque, I. Ramière, and D. Ryckelynck. Hybrid hyper-reduced modeling for contact mechanics problems. *Internat. J. Numer. Methods Engrg.*, 115(1):117-139, 2018.

[15] G. Fichera. Problemi elastostatici con vincoli unilaterali: Il problema di Signorini con ambigue condizioni al contorno. *Atti Accad. Naz. Lincei Mem. Cl. Sci. Fis. Mat. Natur. Sez. Ia (8)*, 7:91-140, 1963/64.

[16] M. Fortin and R. Glowinski. *Augmented Lagrangian methods*, volume 15 of *Studies in Mathematics and its Applications*. North-Holland Publishing Co., Amsterdam, 1983. Applications to the numerical solution of boundary value problems, Translated from the French by B. Hunt and D. C. Spicer.

[17] T. Gustafsson, R. Stenberg, and J. Videman. On Nitsche's method for elastic contact problems. *SIAM J. Sci. Comput.*, 42(2):B425-B446, 2020.

[18] T. Gustafsson and J. Videman. Stabilized finite elements for tresca friction problem. *ESAIM Math. Model. Numer. Anal.*, 56(4):1307-1326, 2022.

[19] B. Haasdonk. Convergence rates of the POD-greedy method. *ESAIM Math. Model. Numer. Anal.*, 47(3):859-873, 2013.

[20] B. Haasdonk, J. Salomon, and B. Wohlmuth. A reduced basis method for parametrized variational inequalities. *SIAM J. Numer. Anal.*, 50(5):2656-2676, 2012.

[21] J. S. Hesthaven, G. Rozza, and B. Stamm. *Certified reduced basis methods for parametrized partial differential equations*, volume 590. Springer, 2016.
[22] K. L. Johnson. *Contact mechanics*. Cambridge University Press, 1987.

[23] N. Kikuchi and J. T. Oden. *Contact Problems in Elasticity: A Study of Variational Inequalities and Finite Element Methods*, volume 8 of SIAM Studies in Applied Mathematics. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 1988.

[24] K. S. Kollepara, J. M. Navarro-Jiménez, Y. Le Guennec, L. Silva, and J. V. Aguado. On the limitations of low-rank approximations in contact mechanics problems. *Internat. J. Numer. Methods Engrg.*, 124(1):217-234, 2023.

[25] K. Kunisch and S. Volkwein. Galerkin proper orthogonal decomposition methods for parabolic problems. *Numer. Math.*, 90(1):117-148, 2001.

[26] S. Le Berre, I. Ramière, J. Fauque, and D. Ryckelynck. Condition number and clustering-based efficiency improvement of reduced-order solvers for contact problems using Lagrange multipliers. *Mathematics*, 10:1495-1520, 2022.

[27] Y. Maday, N. C. Nguyen, A. T. Patera, and G. S. H. Pau. A general multipurpose interpolation procedure: the magic points. *Commun. Pure Appl. Anal.*, 8(1):383-404, 2009.

[28] R. Mlika, Y. Renard, and F. Chouly. An unbiased Nitsche's formulation of large deformation frictional contact and self-contact. *Comput. Methods Appl. Mech. Engrg.*, 325:265-288, 2017.

[29] I. Niakh, G. Drouet, V. Ehrlacher, and A. Ern. Stable model reduction for linear variational inequalities with parameter-dependent constraints. M2AN, to appear, https://hal.archives-ouvertes.fr/hal-03611982, March 2022.

[30] J. Nitsche. Über ein Variationsprinzip zur Lösung von Dirichlet-Problemen bei Verwendung von Teilräumen, die keinen Randbedingungen unterworfen sind. *Abh. Math. Sem. Univ. Hamburg*, 36:9-15, 1971.

[31] C. Prud'Homme, D. V. Rovas, K. Veroy, L. Machiels, Y. Maday, A. T. Patera, and G. Turinici. Reliable real-time solution of parametrized partial differential equations: Reduced-basis output bound methods. *J. Fluids Eng.*, 124(1):70-80, 2002.

[32] A. Quarteroni, A. Manzoni, and F. Negri. *Reduced basis methods for partial differential equations: an introduction*, volume 92. Springer, 2015.

[33] Y. Renard and K. Poulios. GetFEM: Automated FE modeling of multiphysics problems based on a generic weak form language. https://hal.archives-ouvertes.fr/hal-02532422, 2020.

[34] M. Sofonea and A. Matei. *Mathematical models in contact mechanics*, volume 398. Cambridge University Press, 2012.

[35] G. Stampacchia. Formes bilinéaires coercitives sur les ensembles convexes. *C. R. Acad. Sci. Paris*, 258:4413-4416, 1964.

[36] P. Wriggers. *Computational Contact Mechanics*, volume 2. Springer Berlin, Heidelberg, 2006.

[37] D. Zeka, P.-A. Guidault, D. Néron, M. Guiton, and G. Enchéry. Preliminary study for the simulation of wire ropes using a model reduction approach suitable for multiple contacts. In *25ème Congrès Français de Mécanique*, 2022.
[22] K. L. Johnson. *Contact mechanics*. Cambridge University Press, 1987.

[23] N. Kikuchi and J. T. Oden. *Contact Problems in Elasticity: A Study of Variational Inequalities and Finite Element Methods*, volume 8 of SIAM Studies in Applied Mathematics. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 1988.

[24] K. S. Kollepara, J. M. Navarro-Jiménez, Y. Le Guennec, L. Silva, and J. V. Aguado. On the limitations of low-rank approximations in contact mechanics problems. *Internat. J. Numer. Methods Engrg.*, 124(1):217-234, 2023.

[25] K. Kunisch and S. Volkwein. Galerkin proper orthogonal decomposition methods for parabolic problems. *Numer. Math.*, 90(1):117-148, 2001.

[26] S. Le Berre, I. Ramière, J. Fauque, and D. Ryckelynck. Condition number and clustering-based efficiency improvement of reduced-order solvers for contact problems using Lagrange multipliers. *Mathematics*, 10:1495-1520, 2022.

[27] Y. Maday, N. C. Nguyen, A. T. Patera, and G. S. H. Pau. A general multipurpose interpolation procedure: the magic points. *Commun. Pure Appl. Anal.*, 8(1):383-404, 2009.

[28] R. Mlika, Y. Renard, and F. Chouly. An unbiased Nitsche's formulation of large deformation frictional contact and self-contact. *Comput. Methods Appl. Mech. Engrg.*, 325:265-288, 2017.

[29] I. Niakh, G. Drouet, V. Ehrlacher, and A. Ern. Stable model reduction for linear variational inequalities with parameter-dependent constraints. M2AN, to appear, https://hal.archives-ouvertes.fr/hal-03611982, March 2022.

[30] J. Nitsche. Über ein Variationsprinzip zur Lösung von Dirichlet-Problemen bei Verwendung von Teilräumen, die keinen Randbedingungen unterworfen sind. *Abh. Math. Sem. Univ. Hamburg*, 36:9-15, 1971.

[31] C. Prud'Homme, D. V. Rovas, K. Veroy, L. Machiels, Y. Maday, A. T. Patera, and G. Turinici. Reliable real-time solution of parametrized partial differential equations: Reduced-basis output bound methods. *J. Fluids Eng.*, 124(1):70-80, 2002.

[32] A. Quarteroni, A. Manzoni, and F. Negri. *Reduced basis methods for partial differential equations: an introduction*, volume 92. Springer, 2015.

[33] Y. Renard and K. Poulios. GetFEM: Automated FE modeling of multiphysics problems based on a generic weak form language. https://hal.archives-ouvertes.fr/hal-02532422, 2020.

[34] M. Sofonea and A. Matei. *Mathematical models in contact mechanics*, volume 398. Cambridge University Press, 2012.

[35] G. Stampacchia. Formes bilinéaires coercitives sur les ensembles convexes. *C. R. Acad. Sci. Paris*, 258:4413-4416, 1964.

[36] P. Wriggers. *Computational Contact Mechanics*, volume 2. Springer Berlin, Heidelberg, 2006.

[37] D. Zeka, P.-A. Guidault, D. Néron, M. Guiton, and G. Enchéry. Preliminary study for the simulation of wire ropes using a model reduction approach suitable for multiple contacts. In *25ème Congrès Français de Mécanique*, 2022.