# The augmented Lagrangian method as a framework for stabilised methods in
computational mechanics

Erik Burman  Peter Hansbo  Mats G. Larson

###### Abstract

In this paper we will review recent advances in the application of the
augmented Lagrange multiplier method as a general approach for generating
multiplier–free stabilised methods. We first show how the method generates
Galerkin/Least Squares type schemes for equality constraints and then how it
can be extended to develop new stabilised methods for inequality constraints.
Application to several different problems in computational mechanics is given.

##  1 Introduction

The Augmented Lagrangian Method (ALM) has a long history in optimisation. In
its standard form it can be seen as augmenting standard Lagrange multiplier
methods with a penalty term, penalising the constraint equations. It was
introduced in order to combine the advantages of the penalty method and the
multiplier method in the context on constrained optimisation independently by
Hestenes and Powell in [72, 89]. It was then extended to the case of
optimization with inequality constraints by Rockafellar in [90, 91]. Soon
afterwards the potential of ALM for the numerical approximation of partial
differential equations (pde) and computational mechanics was explored in
Glowinski and Morocco [57] and by Fortin in [51]. For overviews of the early
results on augmented Lagrangian methods for approximation of pde we refer to
the monographs by Glowinski and coworkers [52, 55].

In computational mechanics, Lagrangian methods have the drawback of having to
fulfil an inf–sup condition to ensure stability of the discrete scheme such
that the balance between the discretisation of the primal variable and the
multiplier variable must be chosen carefully. Adding a penalty term does not
change this situation, and in computational mechanics ALM has therefore been
used mostly in an iterative approach (improving the conditioning of the
discrete system) [52, 55, 94, 104, 82, 105, 102], or as a way of strengthening
control of the constraints in cases where the discretisation is under-
constrained. It was also shown to improve convergence in some cases by making
the penalty parameter mesh dependent in [14]. Recently similar ideas have been
applied in the context of preconditioning solution methods for discretisations
of incompressible flows [50, 86]. The ideas of extending the ALM to
variational inequalities of [90, 91] were introduced in the context of contact
mechanics by Alart and Curnier in [1].

An early approach to weak boundary conditions for finite element methods was
introduced orignally by Nitsche in [85], using a method that is related to
ALM, but without any multiplier. Indeed here the multiplier has been replaced
by its physical representation, the normal boundary flux. Only recently this
possibility of substituting the multiplier by its physical interpretation in
the discrete augmented Lagrangian formulation has been explored in its
generality. This approach gives rise to schemes that are formally equivalent
to stabilised Lagrange multiplier methods, where the stabilisation is of
Galerkin/Least Squares (GLS) type [97].

There is, however a crucial difference between the ALM and GLS stabilisation
method, and that is the treatment of variational inequalities. The classical
GLS formulation for variational inequalities of Barbosa and Hughes [5] is very
close to standard multiplier schemes, whereas the ALM supplies an alternative
way to define the stabilisation mechanism which transforms the variational
inequalities to nonlinear equalities to which iterative schemes can be readily
applied.

There is a very large literature on on variational inequalities in pde and we
can not survey the whole field herein. Below we will focus on works on finite
element method formulation and error analysis. For theoretical background
material relevant to the material herein we refer to [43, 81, 44] and for a
review of computational aspects including design of special finite element
spaces, adaptive method and solvers we refer to [103] and references therein.

The theoretical foundation for finite element approximation of variational
inequalities was laid in the seminal works by Falk [48, 49], by Brezzi et al.
[17, 18] and Haslinger [69]. For early overviews on computational aspects we
refer to the monographs by Glowinski and co-workers [56, 54] and Kikuchi and
Oden [80]. More recent studies of the numerical analysis of finite element
methods for variational inequalities include [73, 8, 30, 9, 74, 7, 76]. For
further work on mixed finite element methods we refer to [70, 92, 40, 10, 95,
11, 93]. For stabilised finite element methods in the context of variational
inequalities see [5, 71, 75, 66, 59, 60, 62]. More recently discontinuous
Galerkin methods and other non-conforming methods allowing for polygonal
elements have been developed for different types of contact problems [100, 99,
28, 106, 107, 53, 101, 39]. Another recent development is the application of
isogeometric analysis to contact problems [98, 41, 77, 2]. Some results on
fourth order problems have been reported in [68, 92, 16, 15, 58, 62]. Some
early error analyses for augmented Lagrangian finite element methods applied
to variational inequalities have been proposed in [30, 79].

Optimal error estimates for the unilateral contact problem however remained
elusive and typically required some additional assumptions on the interface
between the zones of contact and no contact. The Nitsche ALM, where the
multiplier is replaced by its physical interpretation, was first introduced
and analysed for variational inequalities by Chouly and Hild [32] in the
setting of friction free small deformation elastic contact (without explicit
reference to augmented Lagrangians). In this context they also showed optimal
error estimates without additional a priori assumptions on the contact set. A
similar result for the Signorini problem using a Lagrange multiplier approach
(without ALM) was derived in [42]. The idea of using ALM with eliminated
multiplier for contact problems was then extended to various other models in
[31, 36, 34, 35, 38, 33]; for an overview, cf. [37]. Finite element methods
using ALM in the form of a nonlinear equality without eliminating the
multiplier was analysed in [25]. In the context of non-conforming
approximation the approach has been applied in [29] and using IGA in [46, 77].
It has been explored for CutFEM applications in [47, 23], for obstacle
problems in [27, 24], and for Signorini boundaries in the plate model in [26].
Typically the analysis of Nitsche’s method requires some additional regularity
assumptions in order to make sense of the non-conforming terms and we will
consider this case below. An analysis for low regularity solutions for Nitsche
type methods applied to contact problems was proposed in [62, 63]. The
reformulation of the variational inequality as a nonlinear equality with
elimination of the multiplier is also advantageous in multi physics
applications as illustrated in [20, 21] and to impose positivity in flow
problems [19].

Our main objective in this paper is to introduce the ALM in a model context,
starting with the original formulation for optimization under constraints and
then presenting the extension to pde approximation in an abstract framework.
Particular focus will be given to variational inequalities that are rewritten
as nonlinear equalities in the ALM framework. Here we prove existence and best
approximation estimates for the multiplier method under the assumption of
sufficient smoothness of the multiplier. We discuss stabilised methods and
sketch how these results generalize to the case where the multiplier is
eliminated. The versatility of the approach is then shown by applying it in
some different settings.

In section 2 we start by recalling the augmented Lagrangian method in the
finite dimensional setting both for equality and inequality constraints and
derive the augmented Lagrangian formulation for inequality constraints using
the equality constraint formulation and slack variables. In section 4 we then
discuss the use of the augmented Lagrangian in the context of partial
differential equations and present the properties of the formulation in an
abstract framework. We show how the necessary a priori bounds for existence of
discrete solutions are obtained and we derive best approximation estimates for
the augmented Lagrangian finite element method. In section 5 we proceed and
give a number of different applications drawing from fluid and solid
mechanics. The paper finishes with some numerical experiments in section 6
showing the versatility of the proposed framework.

##  2 The finite dimensional setting

We begin by recalling the ALM for finite dimensional optimisation problems and
by giving an informal introduction to some key ideas to be used in the
following. Below we will frequently use the notation $a\lesssim b$ for $a\leq
Cb$.

###  2.1 Optimisation with equality constraints

We consider the quadratic optimisation problem:

$\displaystyle\min_{{\boldsymbol{x}}\in\mathbb{R}^{n}}f({\boldsymbol{x}})\quad\text{subject
to}\quad g_{i}({\boldsymbol{x}})=0,\quad i=1,\dots,m$ (2.1)

This problem can be solved by the Lagrange multiplier method, seeking
stationary points to the function

$\mathcal{L}({\boldsymbol{x}},\lambda_{1},\ldots,\lambda_{m})=f({\boldsymbol{x}})+\sum_{i}\lambda_{i}g_{i}({\boldsymbol{x}})$
(2.2)

solving the system of equations

$\displaystyle\nabla f-\sum_{i}\lambda_{i}\nabla g_{i}={}$ $\displaystyle{\bf
0}$ (2.3) $\displaystyle g_{i}={}$ $\displaystyle 0,\;i=1,\ldots,m$ (2.4)

It can also be solved approximately by the penalty method, seeking the minimum
to the function

$\mathcal{L}_{\gamma}({\boldsymbol{x}})=f({\boldsymbol{x}})+\frac{\gamma}{2}\sum_{i}g_{i}({\boldsymbol{x}})^{2}$
(2.5)

where $\gamma\in\mathbb{R}^{+}$ is a given (large) penalty parameter. We note
that the penalty method has a strong regularising effect on the problem in the
sense that if some of the side conditions are (close to being) linear
combinations of each other, this does not matter; indeed even if
$g_{j}({\boldsymbol{x}})=g_{1}({\boldsymbol{x}})$ for all $j$ we simply solve

$\mathcal{L}_{\gamma}({\boldsymbol{x}})=f({\boldsymbol{x}})+m\frac{\gamma}{2}g_{1}({\boldsymbol{x}})^{2}$
(2.6)

which is a well posed problem. This is not the case in the multiplier method,
where the system (2.3)–(2.4) would then be ill posed. The key point is that
the side conditions do not come into play explicitly in the penalty method. On
the other hand, in general the minimiser of (2.5) coincides with that of (2.1)
only in the limit as $\gamma\rightarrow\infty$. The ALM is a combination of
the penalty method and the multiplier method: seek the stationary point to

$\displaystyle\mathcal{L}_{\gamma}({\boldsymbol{x}},{\boldsymbol{\lambda}})$
$\displaystyle=f({\boldsymbol{x}})-\sum_{i}\lambda_{i}g_{i}({\boldsymbol{x}})+\frac{\gamma}{2}\sum_{i}g_{i}({\boldsymbol{x}})^{2}$
(2.7)

This problem has the same stationary point as (2.2) and the same stability
problem in case of linearly independent side conditions. We note, however,
that the multiplier can be eliminated by first solving (2.3), which we
symbolically denote by

$\lambda_{i}=\frac{df}{dg_{i}}({\boldsymbol{x}})$ (2.8)

(the multipliers can be interpreted as the change in objective with respect to
change in the corresponding side condition), and seek the minimum to the
reduced Lagrangian

$\mathcal{L}_{A}({\boldsymbol{x}})=f({\boldsymbol{x}})+\sum_{i}\frac{\gamma}{2}g^{2}_{i}({\boldsymbol{x}})-\frac{df}{dg_{i}}({\boldsymbol{x}})g_{i}({\boldsymbol{x}})$
(2.9)

Like in the penalty method, the side conditions are then no longer explicit;
however, in case of linear dependence we still have an ill posed problem in
solving (2.3) and we cannot obtain the representation (2.8). But say that we
had an alternative way of computing the multiplier so that symbolically we had

$\lambda_{i}^{*}({\boldsymbol{x}})\approx\frac{df}{dg_{i}}({\boldsymbol{x}}),\quad\lambda_{i}^{*}({\boldsymbol{x}})\;\text{computable}$
(2.10)

Then we could consider the problem of minimising

$\mathcal{L}_{A}^{*}({\boldsymbol{x}})=f({\boldsymbol{x}})+\sum_{i}\left(\frac{\gamma}{2}g_{i}({\boldsymbol{x}})^{2}-\lambda_{i}^{*}({\boldsymbol{x}})g_{i}({\boldsymbol{x}})\right)$
(2.11)

The accuracy of this method would then depend on the accuracy of the
approximation (2.10) and the stability of the formulation. A typical situation
is that there is a constant such that

$\sum_{i}|\lambda_{i}^{*}({\boldsymbol{x}})|^{2}\leq Cf({\boldsymbol{x}})$
(2.12)

which gives

$\displaystyle\mathcal{L}_{A}^{*}({\boldsymbol{x}})$
$\displaystyle=f({\boldsymbol{x}})+\sum_{i}\left(\frac{\gamma}{2}g_{i}({\boldsymbol{x}})^{2}-\lambda_{i}^{*}({\boldsymbol{x}})g_{i}({\boldsymbol{x}})\right)$
(2.13)
$\displaystyle=f({\boldsymbol{x}})+\sum_{i}\left(\frac{\gamma}{2}g_{i}({\boldsymbol{x}})^{2}-\delta|\lambda_{i}^{*}({\boldsymbol{x}})|^{2}-\frac{1}{4\delta}g_{i}^{2}({\boldsymbol{x}})\right)$
(2.14) $\displaystyle\geq
f({\boldsymbol{x}})-\delta\sum_{i}|\lambda_{i}^{*}({\boldsymbol{x}})|^{2}+\sum_{i}\left(\frac{\gamma}{2}-\frac{1}{4\delta}\right)g_{i}^{2}({\boldsymbol{x}})$
(2.15) $\displaystyle\geq(1-\delta
C)f({\boldsymbol{x}})+\sum_{i}\left(\frac{\gamma}{2}-\frac{1}{4\delta}\right)g_{i}^{2}({\boldsymbol{x}})$
(2.16) $\displaystyle\gtrsim
f({\boldsymbol{x}})+\sum_{i}g_{i}^{2}({\boldsymbol{x}})$ (2.17)

where we obtained the last estimate by taking $\delta$ sufficiently small and
$\gamma$ sufficiently large. We conclude that the minimization problem for
$\mathcal{L}_{A}^{*}({\boldsymbol{x}})$ is well posed if $\gamma>\gamma_{C}$.
This is the basic idea that underlies the application of the ALM as a
stabilisation method, in cases where the multiplier can be eliminated.

###  2.2 Optimisation with inequality constraints

We consider next a quadratic optimisation problems of the type:

$\displaystyle\min_{x\in\mathbb{R}^{n}}f({\boldsymbol{x}})\quad\text{subject
to}\quad g_{i}({\boldsymbol{x}})\leq 0,\quad i=1,\dots,m$ (2.18)

The augmented Lagrangian for this problem proposed by Rockafellar [90,
Equation (7)] (here with $\gamma=2r$, and with the multiplier chosen negative)
takes the form for $\gamma\in\mathbb{R}^{+}$,

$\displaystyle\mathcal{L}_{A}({\boldsymbol{x}},{\boldsymbol{\lambda}})=$
$\displaystyle f({\boldsymbol{x}})+\frac{1}{2\gamma}\sum_{i}\left([\gamma
g_{i}({\boldsymbol{x}})-\lambda_{i}]_{+}^{2}-\lambda_{i}^{2}\right)$ (2.19)

where $[x]_{+}=\max(x,0)$.

Observe that another equivalent reformulation is given by

$\displaystyle\mathcal{L}_{A}({\boldsymbol{x}},{\boldsymbol{\lambda}})=$
$\displaystyle
f({\boldsymbol{x}})-\sum_{i}\lambda_{i}g_{i}({\boldsymbol{x}})+\frac{\gamma}{2}\sum_{i}g_{i}({\boldsymbol{x}})^{2}-\frac{1}{2\gamma}\sum_{i}[\gamma
g_{i}({\boldsymbol{x}})-\lambda_{i}]_{-}^{2}$ (2.20)

where $[x]_{-}=\min(x,0)$. This is easily seen by using that
$x=[x]_{+}+[x]_{-}$ and hence

$[x]_{+}^{2}=([x]_{+}+[x]_{-})^{2}-[x]_{-}^{2}-\underbrace{2[x]_{+}[x]_{-}}_{=0}=x^{2}-[x]_{-}^{2}$

Applying this in (2.19) with $x=\gamma g_{i}({\boldsymbol{x}})-\lambda_{i}$
leads to (2.20). In (2.20) we recognise the augmented Lagrangian for the
equality constraint (2.7) in the first three terms and the last term is the
non-linear switch that introduces the inequality constraint.

To see that (2.19) is indeed the natural formulation we introduce slack
variables $z_{i}\in\mathbb{R}_{+}$ and rewrite (2.18) in the form

$\displaystyle\min_{({\boldsymbol{x}},{\boldsymbol{z}})\in\mathbb{R}^{n}\times\mathbb{R}^{m}_{+}}f({\boldsymbol{x}})\quad\text{subject
to}\quad g_{i}({\boldsymbol{x}})+z_{i}=0,\quad i=1,\dots,m$ (2.21)

with corresponding augmented Lagrangian

$\displaystyle\mathcal{L}_{A}({\boldsymbol{x}},{\boldsymbol{z}},{\boldsymbol{\lambda}})=f({\boldsymbol{x}})-\sum_{i}\left\\{(g_{i}({\boldsymbol{x}})+z_{i})\lambda_{i}+\frac{\gamma}{2}(g_{i}({\boldsymbol{x}})+z_{i})^{2}\right\\}$
(2.22)

for which we seek stationary points, minimizing in
$({\boldsymbol{x}},{\boldsymbol{z}})$. Here we may now perform the
optimization over ${\boldsymbol{z}}\in\mathbb{R}_{+}^{m}$ explicitly by noting
that for each ${\boldsymbol{x}}$ and ${\boldsymbol{\lambda}}$ we obtain a sum
of quadratic polynomials in $z_{i}$ of the form

$\displaystyle-(g_{i}({\boldsymbol{x}})+z_{i})\lambda_{i}+\frac{\gamma}{2}(g_{i}({\boldsymbol{x}})+z_{i})^{2}=\frac{1}{2\gamma}\Big{(}(\gamma(g_{i}({\boldsymbol{x}})+z_{i})-\lambda_{i})^{2}-\lambda_{i}^{2}\Big{)}$
(2.23)

and therefore the minimum is attained at $\gamma z_{i}=-(\gamma
g_{i}-\lambda_{i})$ and taking the constraint $z_{i}\in\mathbb{R}_{+}$ into
account we find that $\gamma z_{i}=[-(\gamma
g_{i}({\boldsymbol{x}})-\lambda_{i})]_{+}$. Inserting this expression for
$\gamma z_{i}$ into (2.23) and using the identity $a+[-a]_{+}=[a]_{+}$ we
arrive at

$\displaystyle\mathcal{L}_{A}({\boldsymbol{x}},{\boldsymbol{z}},{\boldsymbol{\lambda}})=f({\boldsymbol{x}})+\frac{1}{2\gamma}\sum_{i}([\gamma
g_{i}({\boldsymbol{x}})-\lambda_{i}]_{+}^{2}-\lambda_{i}^{2})$ (2.24)

Alternatively we may seek stationary points to the standard Lagrangian

$\mathcal{L}({\boldsymbol{x}},\lambda)=f({\boldsymbol{x}})-\sum_{i}\lambda_{i}g_{i}({\boldsymbol{x}})$
(2.25)

under the Karush–Kuhn–Tucker (KKT) conditions

$\displaystyle g_{i}\leq{}$ $\displaystyle 0,\;i=1,\ldots,m$ (2.26)
$\displaystyle\lambda_{i}\leq{}$ $\displaystyle 0,\;i=1,\ldots,m$ (2.27)
$\displaystyle\lambda_{i}g_{i}={}$ $\displaystyle 0,\;i=1,\ldots,m$ (2.28)

Noting that the KKT conditions (2.26)–(2.28) are equivalent to the single
statement

$\lambda_{i}=-[\gamma g_{i}-\lambda_{i}]_{+}$ (2.29)

where $\gamma\in\mathbb{R}^{+}$ is an arbitrary positive number. We may then
rewrite the Lagrangian in the form

$\displaystyle f({\boldsymbol{x}})-\sum_{i}\lambda_{i}g_{i}({\boldsymbol{x}})$
$\displaystyle=f({\boldsymbol{x}})-\sum_{i}\lambda_{i}\Big{(}g_{i}({\boldsymbol{x}})-\frac{1}{\gamma}\lambda_{i}\Big{)}-\frac{1}{\gamma}\lambda_{i}^{2}$
(2.30)
$\displaystyle=f({\boldsymbol{x}})+\frac{1}{\gamma}\sum_{i}\underbrace{[\gamma
g_{i}({\boldsymbol{x}})-\lambda_{i}]_{+}(\gamma
g_{i}({\boldsymbol{x}})+\lambda_{i})}_{[\gamma
g_{i}({\boldsymbol{x}})-\lambda_{i}]_{+}^{2}}-\frac{1}{\gamma}\lambda_{i}^{2}$
(2.31)

where we used (2.29) and the fact that $[a]_{+}a=[a]_{+}^{2}$. The
substitutions $\lambda_{i}\mapsto\lambda_{i}/2$ and $\gamma\mapsto 2\gamma$
manufactures the Lagrangian (2.24).

Writing the optimality system of (2.24) results in the system of equations

$\displaystyle\nabla f+\sum_{i}[\gamma g_{i}-\lambda_{i}]_{+}\nabla g_{i}={}$
$\displaystyle{\bf 0}$ (2.32) $\displaystyle[\gamma g_{i}-\lambda_{i}]_{+}={}$
$\displaystyle-\lambda_{i},\;i=1,\ldots,m$ (2.33)

which is a nonlinear equality problem which explicitly includes the KKT
conditions.

Again, if we can use (2.10) we may instead seek the minima to

$\mathcal{L}_{A}^{*}({\boldsymbol{x}}):=f({\boldsymbol{x}})+\frac{1}{2\gamma}\sum_{i}[\gamma
g_{i}({\boldsymbol{x}})-\lambda_{i}^{*}({\boldsymbol{x}})]_{+}^{2}-(\lambda_{i}^{*}({\boldsymbol{x}}))^{2}$
(2.34)

##  3 Iterative solution using the augmented Lagrangian

The augmented Lagrangian is possibly most well known as the basis for an
iterative algorithm for constrained optimization problems. The stationary
points of the functional (2.7) can be approximated using the following
classical algorithm attributed to Usawa, with the application to augmented
Lagrangian methods developed in the works of Glowinski and co-workers [3, 56,
52, 55]. Following [52] we consider the situation where the model problem is
to minimize

$J({\boldsymbol{x}}):=\frac{1}{2}{\boldsymbol{x}}^{T}{\boldsymbol{A}}{\boldsymbol{x}}-{\boldsymbol{b}}^{T}{\boldsymbol{x}}$

over ${\boldsymbol{x}}\in\mathbb{R}^{n}$ under the constraint
${\boldsymbol{B}}{\boldsymbol{x}}={\boldsymbol{c}}\in\mathbb{R}^{m}$. Here
${\boldsymbol{A}}\in\mathbb{R}^{n\times n}$ is symmetric positive definite,
${\boldsymbol{b}}\in\mathbb{R}^{n}$ and
${\boldsymbol{B}}\in\mathbb{R}^{m\times n}$. The augmented Lagrangian (2.7)
then takes the form,

$\mathcal{L}_{A}({\boldsymbol{x}},{\boldsymbol{\lambda}}):=\frac{1}{2}{\boldsymbol{x}}^{T}{\boldsymbol{A}}{\boldsymbol{x}}-{\boldsymbol{b}}^{T}{\boldsymbol{x}}+{\boldsymbol{\lambda}}^{T}({\boldsymbol{B}}{\boldsymbol{x}}-{\boldsymbol{c}})+\frac{\gamma}{2}|{\boldsymbol{B}}\bf
x-{\boldsymbol{c}}|^{2}$ (3.1)

(Uzawa’s algorithm)

  1. 1.

Let ${\boldsymbol{\lambda}}^{0}\in\mathbb{R}^{m}$

  2. 2.

Find ${\boldsymbol{x}}^{n}\in\mathbb{R}^{n}$ such that
$\mathcal{L}_{A}({\boldsymbol{x}}^{n},{\boldsymbol{\lambda}}^{n})\leq\mathcal{L}_{A}({\boldsymbol{y}}^{n},{\boldsymbol{\lambda}}^{n})$
for all ${\boldsymbol{y}}^{n}\in\mathbb{R}^{n}$

  3. 3.

Update the multiplier:
${\boldsymbol{\lambda}}^{n+1}={\boldsymbol{\lambda}}^{n}+\rho_{n}({\boldsymbol{B}}{\boldsymbol{x}}^{n}-{\boldsymbol{c}}),\quad\rho_{n}>0$.

We note that step 2 is equivalent to solving the linear system, find
${\boldsymbol{x}}^{n}\in\mathbb{R}^{n}$ such that

$({\boldsymbol{A}}+\gamma{\boldsymbol{B}}^{T}{\boldsymbol{B}}){\boldsymbol{x}}^{n}=-{\boldsymbol{B}}^{T}{\boldsymbol{\lambda}}^{n}+{\boldsymbol{b}}+\gamma{\boldsymbol{B}}^{T}{\boldsymbol{c}}$

The iterates ${\boldsymbol{x}}^{n},{\boldsymbol{\lambda}}^{n}$ of the
iterative method converges to the saddle point of (3.1) provided the
steplength $\rho_{n}$ satisfies

$0<\alpha_{0}\leq\rho_{n}\leq\alpha_{1}<2\left(\gamma+\frac{1}{\beta^{2}}\right)$

where $\beta^{2}$ is the largest eigenvalue of the matrix
${\boldsymbol{A}}^{-1}{\boldsymbol{B}}^{T}{\boldsymbol{B}}$ defined by

$\beta^{2}=\max_{{\boldsymbol{v}}\neq
0}\frac{|{\boldsymbol{B}}{\boldsymbol{v}}|^{2}}{{\boldsymbol{v}}^{T}{\boldsymbol{A}}{\boldsymbol{v}}}$

For a proof of the convergence result we refer to [56, Chapter 2, Section 4]
or [52, Chapter 1, Section 2].

##  4 Augmented Lagrangian methods and Galerkin/Least squares

We now turn to the case where the Lagrangian is a functional taking values in
some Sobolev space and the numerical method is obtained by finding the
stationary points in a finite dimensional approximation space. Typically we
are interested in the discretisation of a problem where some energy is
minimised under a constraint. To illustrate this we consider the case with
equality constraints. Let $V$ and $H$ denote two Hilbert spaces, with dual
spaces $V^{\prime}$ and $H^{\prime}$, respectively. Let
$F:V\rightarrow\mathbb{R}$ denote a strictly convex $C^{2}$-functional and
$B:V\rightarrow H$ a linear operator. We are interested in minimising $F$
under a constraint defined by $B$. Given the data $f\in V^{\prime}$ and $g\in
H$ We consider the optimization problem

$u=\mbox{arginf}_{v\in V}F(v)-\left<f,v\right>_{V^{\prime},V}\mbox{ such that
}Bu=g.$ (4.1)

The Lagrangian takes the form

$\mathcal{L}(v,\mu):=F(v)-\left<f,v\right>_{V^{\prime},V}-\left<\mu,Bv-g\right>_{H^{\prime},H}.$
(4.2)

This problem can be shown to have unique solution under suitable hypothesis on
the spaces $V$ and $H$ and the operators $F$, $B$, $f$ and $g$ (see for
instance [83, Chapter 1, Section 2.1, Theorems 2.1 and 2.2]). Augmenting the
Lagrangian has no effect on the continuous level, but formally an augmented
version of (4.2), in the spirit of (2.7) can be written

$\mathcal{L}_{\text{A}}(v,\mu):=F(v)-\left<f,v\right>_{V^{\prime},V}-\left<\mu,Bv-g\right>_{H^{\prime},H}+\frac{\gamma}{2}\|Bv-g\|_{H}^{2}.$
(4.3)

The discrete version of the ALM based on (4.3), would then be obtained by
restricting $\mathcal{L}_{\text{A}}$ to finite dimensional spaces. As we saw
in the previous section the ALM on the discrete level combines the control of
the constraint given by the Lagrange multiplier and of the penalty. It also
gives us an iterative procedure to find the minimiser. When using the ALM in
the context of pde problems the ALM also gives enhanced control of the side
condition in the sense of a GaLS method, or a variational multiscale method.
To see this we assume that $H=H^{\prime}=L^{2}$ and that
$H^{\prime}_{h}\subset H^{\prime}$, $V_{h}\subset V$ are some finite
dimensional approximation spaces. Here $h$ denotes the characteristic
lengthscale (or mesh parameter) of the discrete space. We let
$\pi_{H}:H\mapsto H^{\prime}_{h}$ denote the $L^{2}$-orthogonal projection
onto $H^{\prime}_{h}$. Since

$\left<\mu_{h},Bv-g\right>_{L^{2}}=\left<\mu_{h},\pi_{H}(Bv-g)\right>_{L^{2}}$

we see that the Lagrange multiplier only gives control of the projection of
$Bv-g$ on the finite dimensional subspace $H_{h}^{\prime}$. This may be
insufficient for the stability of the method, in particular since
$H_{h}^{\prime}$ may need to be chosen small compared to $V_{h}$ for stability
reasons, i.e. to satisfy the inf-sup stability condition that we will discuss
below. A classical example is the stability of the incompressibility
constraint (in which case $B$ is the divergence operator) of the Brinkman
problem when the viscosity becomes negligible. Adding the term
$\|Bv-g\|_{H}^{2}$ enhances the stability, by adding control of
$(I-\pi_{H})(Bv-g)$ compared to the pure Lagrange multiplier method. This also
shows that a sufficient stabilization can be achieved by augmenting with
$\|(I-\pi_{H})(Bv-g)\|_{H}^{2}$. This we recognise as a stabilization of the
orthogonal subscales, which is a member of the family of variational
multiscale methods. Of course in the associated Euler-Lagrange equations these
terms take the form of GLS stabilizations of some residual quantities. Indeed
a number of ideas from the field of stabilized methods can be made to bear to
the ALM, but we will not explore this further herein. Instead we will show in
the examples below how the design of finite element methods using the ALM
allows us to recover some well known GLS methods from computational mechanics.

We can discern two different situations for the continuous problem (4.3):

  1. A.

The multiplier has enough regularity to define a scalar product with the side
condition.

  2. B.

The multiplier has only regularity enough to support a duality pairing with
the side condition.

In the first case we can use an analogue to the reformulation (2.19) which is
convenient for the treatment of inequality conditions, and formulate the
problem on the continuous level; in the second case this is not formally
correct. Indeed if the multiplier does not have sufficient regularity the
augmented continuous formulation does not lead to a well-defined problem,
unless the augmentation is taken in the continuous $H$-norm, which may be
inconvenient from computational standpoint. In this case the reformulation
(2.19) is not available. We emphasize that this is not a problem in the
discrete setting since we can use norm equivalence of discrete spaces to
obtain an ALM that has the right asymptotic scaling. However in order to carry
out a rigorous numerical analysis of the resulting finite element method the
assumption of additional regularity of the exact solution must be justifiable.
This is often, but not always the case. In that sense ALM methods in the
situation B can be seen as a non-conforming method.

For the discrete as well as the continuous problem we have two further cases:

  1. C.

The multiplier has a physical interpretation in terms of the primal variable.

  2. D.

The multiplier cannot be interpreted (or be easily interpreted) in terms of
the primal variable.

For the discrete case, we also have the problem of finding suitable
approximations to fulfil a discrete inf–sup condition. In case C we can use a
trick analogous to that of (2.10), which gives a class of problems where the
multiplier has been eliminated beforehand; alternatively, the multiplier can
be retained and stabilised by the addition of a GLS term, in the spirit of [4,
5]. These approaches give stability without balancing the discretisation of
the multiplier space and the space for the primal variable. In case D the
multiplier has to be retained, but the inequality case can still be handled in
the same way as above and stabilisation is still possible, for instance using
interior penalty stabilization where the stabilization acts on the multiplier
alone [22, 25].

###  4.1 Abstract framework

Since the rationale of the method is from numerical approximation we will only
consider formulations that work in the finite dimensional setting, then A and
B above are treated similarly. However it is only in case A that the
discussion holds also for the continuous case. The resulting numerical methods
can be shown to be optimally converging for sufficiently smooth exact
solutions, but the problem of convergence is not established for exact
solutions that has no additional regularity. The question of how to design
methods that are valid formulations also for the original pde problem is
subtle and requires the design of sophisticated stabilization operators, for
an interesting work in this direction we refer to [12]. Below we let
$\left<\cdot,\cdot\right>$ denote the $L^{2}$ scalar product over the the
domain of definition of functions in $H$ and we denote the associated norm
$\|v\|:=\left<v,v\right>^{\frac{1}{2}}$.

We are interested in minimising $F$ under a constraint defined by $B$, either
as an equality or an inequality constraint. We will now introduce some
sufficient conditions for the abstract analysis below to hold. We will then in
the examples show that the assumptions are verified.

  1. 1.

We assume that the operator $B$ is bounded and surjective from $V$ to $H$, so
that for every $\zeta\in H$ there exists $\xi\in V$ such that $B\xi=\zeta$ and
$\|\xi\|_{V}\leq C\|\zeta\|_{H}$. It follows that there exists $\alpha>0$ such
that for every $\mu\in H^{\prime}$ there holds

$\alpha\|\mu\|_{H^{\prime}}\leq\sup_{v\in
V}\frac{\left<Bv,\mu\right>_{H,H^{\prime}}}{\|v\|_{V}}$ (4.4)

  2. 2.

We also assume that $V_{h}$ and $H_{h}^{\prime}$ are chosen in such a way that
this property carries over to the finite dimensional setting, in the sense
that a so called Fortin interpolant exists, for all $v\in V$ such that $Bv\in
H$, there exists $\exists i_{F}v\in V_{h}$ such that for all $q_{h}\in
H_{h}^{\prime}$,

$\left<B(v-i_{F}v),q_{h}\right>=0,\quad\|i_{F}v\|_{V}+\|Bi_{F}v\|_{H_{h}}\lesssim\|v\|_{V}+\|Bv\|_{H_{h}}$
(4.5)

Note that for $v\in V_{h}$ there holds $i_{F}v=v$.

  3. 3.

We assume that the surjectivity also holds for the discrete spaces on the
following form: for all $\mu_{h}\in H_{h}^{\prime}$ there exists $v_{h}\in
V_{h}$ such that for all $q_{h}\in H_{h}^{\prime}$,

$\left<\mu_{h}-Bv_{h},q_{h}\right>=0,\quad\|v_{h}\|_{V}+\|Bv_{h}\|_{H_{h}}\lesssim\|\mu_{h}\|_{H_{h}}$
(4.6)

Discrete surjectivity is a consequence of the discrete inf-sup condition which
typically is equivalent with the existence of the Fortin interpolant [45,
Lemma 26.9]. We state both (4.5) and (4.6) separately here for future
reference and to highlight the difference of the norms required in the right
hand side. If we are in a non-conforming situation it is not immediately clear
that equivalence holds. Note however that if the spaces are such that
$\|Bv-\pi_{H_{h}}Bv\|_{H_{h}}\leq\|v\|_{V}$ then (4.6) implies (4.5).

The form of the stabilities in (4.5) and (4.6) appear a bit ad hoc here, but
as we shall see below this is the natural stability to require for the
analysis. Here the norm $\|\cdot\|_{H_{h}}$ is an $h$-weighted $L^{2}$-norm
and will be discussed below.

###  4.2 Equality constraints

We wish to solve the optimization problem (4.1) and recall the formal
augmented Lagrangian similar to (2.7) given by

$\mathcal{L}_{\text{A}}(v,\mu):=F(v)-\left<f,v\right>_{V^{\prime},V}-\left<\mu,Bv-g\right>_{H^{\prime},H}+\frac{\gamma}{2}\|Bv-g\|_{H}^{2}$
(4.7)

For later use with inequality constraints, we would now like to use the
analogy to (2.19). However, this is not possible unless $H^{\prime}=H:=L_{2}$,
where $L_{2}$ denotes the space of square integrable functions over the
pertinent domain, which is case A above. In this particular case, completing
the square, $-2ab+b^{2}=(a-b)^{2}-a^{2}$, results in the following equivalent
formulation

$\mathcal{L}_{\text{A}}(v,\mu):=F(v)+\left<f,v\right>_{V^{\prime},V}+\frac{\gamma}{2}\|Bv-g-\frac{1}{2\gamma}\mu\|^{2}-\frac{1}{2\gamma}\|\mu\|^{2}$
(4.8)

analogous to (2.19). We let the semi-linear form $a:V\times
V\rightarrow\mathbb{R}$ be defined by the Gateaux derivative of $F(v)$,

$a(u;v):=\left<\frac{\partial F}{\partial u}(u),v\right>_{V^{\prime},V}$ (4.9)

and we assume that the form $a$ satisfies the positivity, monotonicity and
continuity conditions

$\displaystyle a(v;v)\geq\alpha\|v\|_{V}^{2},\qquad\alpha>0$ (4.10)
$\displaystyle
a(w_{1};w_{1}-w_{2})-a(w_{2};w_{1}-w_{2})\geq\alpha\|w_{1}-w_{2}\|_{V}^{2}$
(4.11) $\displaystyle|a(w_{1};v)-a(w_{2};v)|\leq
C\|w_{1}-w_{2}\|_{V}\|v\|_{V}$ (4.12)

The optimality system obtained by differentiating (4.8) then reads: find
$(u,\lambda)\in V\times H^{\prime}$such that

$\displaystyle
a(u;v)-\left<\lambda,Bv\right>_{H^{\prime},H}-\left<\mu,Bu\right>_{H^{\prime},H}+{\gamma}\left<Bu,Bv\right>_{H}=\left<f,v_{h}\right>_{V^{\prime},V}+\left<g,\mu+{\gamma}Bv\right>_{H}$
(4.13)

for all $(v,\mu)\in V\times H^{\prime}$. Here we simply replace $V$ and
$H^{\prime}$ by $V_{h}$ and $H_{h}^{\prime}$ to obtain the discrete method.

We also want to handle case B. Then typically $Bv\in H:=H^{r}$ where $H^{r}$
denotes a (potentially fractional) Hilbert space with $r>0$, and consequently
$\mu\in H^{\prime}:=H^{-r}$, the dual to $H^{r}$. Since $H^{-r}\not\subset
H^{r}$ the formulation (4.8) no longer makes sense. Instead in the spirit of
discretize first then optimize we move to the discrete counterpart of (4.1)
and introduce discrete spaces $V_{h}\subset V$ and $H^{\prime}_{h}\subset
H^{\prime}$. The finite element method then amounts to seek stationary points
in $V_{h}$ and $H^{\prime}_{h}$ to the augmented Lagrangian (4.7). On the
finite dimensional finite element spaces we can approximate the continuous
norms $\|\cdot\|_{H}$ and $\|\cdot\|_{H^{\prime}}$ by discrete counterparts

$\|Bv\|_{H}^{2}\approx\|Bv\|_{H_{h}}^{2}:=\|h^{-r}Bv\|_{L_{2}}^{2}$ (4.14)

and

$\|\mu\|_{H^{\prime}}^{2}\approx\|\mu\|_{H^{\prime}_{h}}^{2}:=\|h^{r}\mu\|_{L_{2}}^{2}$
(4.15)

where $h$ is the local meshsize (assumed constant in the following for
simplicity) and $r\geq 0$ depends on the space $H$; loosely speaking $r$
corresponds to the number of derivatives present in the norm $\|\cdot\|_{H}$.
It is also immediate by the Cauchy-Schwarz inequality that the following
discrete duality property holds

$\left<v,\mu\right>_{H_{h},H_{h}^{\prime}}:=\left<v,\mu\right>\leq\|v\|_{H_{h}}\|\mu\|_{H^{\prime}_{h}}.$

This is done for two reasons

  1. 1.

To obtain a well conditioned method, we wish to have the same condition number
emanating from the penalty term as from the form $a(\cdot,\cdot)$.

  2. 2.

The analysis of the resulting methods requires that the discrete norms can be
bounded in terms of the form $a(\cdot,\cdot)$ which is only possible if they
scale the same way.

Now we can use the arbitrariness of $\gamma$ to set

$\gamma=\gamma_{0}/h^{2r}$ (4.16)

where $\gamma_{0}$ is a problem– and discretization–dependent constant.
Proceeding as above we find that on discrete spaces

$\mathcal{L}^{h}_{\text{A}}(v,\mu):=F(v)+\left<f,v\right>_{V^{\prime},V}+\frac{\gamma_{0}}{2h^{2r}}\|Bv-g-\frac{h^{2r}}{2\gamma_{0}}\mu\|^{2}-\frac{h^{2r}}{2\gamma_{0}}\|\mu\|^{2}$

and the discrete optimality system reads: find $(u_{h},\lambda_{h})\in
V_{h}\times H^{\prime}_{h}$ such that

$a(u_{h};v)-\left<\lambda_{h},Bv\right>-\left<\mu,Bu_{h}\right>+\frac{\gamma_{0}}{h^{2r}}\left<Bu_{h},Bv\right>=\left<f,v_{h}\right>_{V^{\prime},V}+\left<g,\mu+\frac{\gamma_{0}}{h^{2r}}Bv\right>$
(4.17)

for all $(v,\mu)\in V_{h}\times H^{\prime}_{h}$, where
$\langle\cdot,\cdot\rangle$ denotes the standard $L_{2}$ scalar product.
Introducing the global form

$A[(w,\eta);(v,\mu)]:=a(w;v)-\left<\eta,Bv\right>-\left<\mu,Bw\right>+\frac{\gamma_{0}}{h^{2r}}\left<Bw,Bv\right>,$

we can cast the optimality system on the compact form: find
$(u_{h},\lambda_{h})\in V_{h}\times H^{\prime}_{h}$ such that

$A[(u_{h},\lambda_{h});(v,\mu)]=\left<f,v\right>_{V^{\prime},V}+\left<g,\mu+\frac{\gamma_{0}}{h^{2r}}Bv\right>.$
(4.18)

for all $(v,\mu)\in V_{h}\times H^{\prime}_{h}$.

It follows by inspection that any solution to (4.1) that is sufficiently
smooth, i.e. $(u,\lambda)\in V\times H^{\prime}\cap L^{2}$ is a solution to
(4.17) and hence the formulation is consistent. Indeed the stationary point of
(4.2) is given by the solution to

$a(u;v)-\left<\lambda,Bv\right>_{H^{\prime},H}=\left<f,v\right>_{V^{\prime},V},\quad\forall
v\in V$

and

$\left<Bu,\mu\right>_{H,H^{\prime}}=\left<g,\mu\right>_{H,H^{\prime}}.$

If the solution is sufficiently regular these equalities hold with $H$ and
$H^{\prime}$ replaced by the $L^{2}$ norm and we see that in that case the
exact solution satisfies the finite element formulation,

$\underbrace{a(u;v)-\left<\lambda,Bv\right>}_{=\left<f,v\right>_{V^{\prime},V}}-\underbrace{\left<\mu,Bu\right>+\frac{\gamma_{0}}{h^{2r}}\left<Bu,Bv\right>}_{=\left<g,\mu+\frac{\gamma_{0}}{h^{2r}}Bv\right>}=\left<f,v_{h}\right>_{V^{\prime},V}+\left<g,\mu+\frac{\gamma_{0}}{h^{2r}}Bv\right>.$

We do not give a full analysis of the linear problem herein, but focus on the
nonlinear case in the next section. The analysis immediately also applies to
the linear case.

###  4.3 Inequality constraints

For the subsequent analysis, we will consider the discrete case and hence we
use the space $V_{h}$ for the primal variable and $H_{h}^{\prime}$ for the
dual variable. For simplicity we do not use the subscript $h$ on all variables
below. We wish to solve the continuous optimization problem

$u=\mbox{arginf}_{v\in V}F(v)-\left<f,v\right>_{V^{\prime},V}\mbox{ such that
}Bu\leq 0$ (4.19)

Where the inequality constraint must be interpreted in the sense of
distributions on $H$ and we will denote the continuous multiplier appearing in
the constrained optimization $\lambda\in H^{\prime}$. The weak formulation
characterizing the solution to the continuous problem is as follows. Find
$(u,\lambda)\in V\times K$ (where $K:=\\{\mu\in H^{\prime}:\mu\leq 0\\}$) such
that

$a(u;v)-\left<\lambda,Bv\right>_{H^{\prime},H}=\left<f,v\right>_{V^{\prime},V},\quad\forall
v\in V$ (4.20) $\left<Bu,\lambda-\mu\right>_{H,H^{\prime}}\leq
0,\quad\forall\mu\in K$ (4.21)

It follows by choosing $\lambda-\mu>0$ in (4.21) that $Bu\leq 0$. By taking
$\mu=0$ it follows that $\left<Bu,\lambda\right>_{H,H^{\prime}}\leq 0$ and
since both $Bu$ and $\lambda$ are negative it follows that $\lambda Bu=0$.

We have arrived at the following Kuhn–Tucker conditions on the multiplier and
side condition:

$Bu\leq 0,\quad\lambda\leq 0,\quad\lambda Bu=0.$ (4.22)

We now use the analogue to (2.29), to show that (4.22) formally is equivalent
to

$\lambda=-{\gamma}\,[Bu-\gamma^{-1}\,\lambda]_{+}$ (4.23)

To derive the finite element formulation we also proceed formally following
the discussion of section 2.1 applied to the problem (4.19) with the min taken
over the finite dimensional space $V_{h}$ and write the augmented Lagrangian,
for $\gamma\in\mathbb{R}^{+}$, $(v,\mu)\in V_{h}\times H_{h}^{\prime}$,

$\mathcal{L}_{\text{A}}(v,\mu):=F(v)-\left<f,v\right>_{V^{\prime},V}+\frac{\gamma}{2}\|[Bv-\mu/\gamma]_{+}\|^{2}-\frac{1}{2\gamma}\|\mu\|^{2}$
(4.24)

we note that if $\gamma$ is chosen as in (4.16) we may use (4.14) and (4.15)
to write

$\mathcal{L}_{\text{A}}(v,\mu):=F(v)-\left<f,v\right>_{V^{\prime},V}+\frac{\gamma_{0}}{2}\|[Bv-\mu/\gamma]_{+}\|_{H_{h}}^{2}-\frac{1}{2\gamma_{0}}\|\mu\|_{H_{h}^{\prime}}^{2}$
(4.25)

The finite element optimality system reads: find $(u_{h},\lambda_{h})\in
V_{h}\times H_{h}^{\prime}$ such that

$A[(u_{h},\lambda_{h});(v,\mu)]=\left<f,v\right>_{V^{\prime},V}$ (4.26)

for all $(v,\mu)\in V_{h}\times H_{h}^{\prime}$, where

$A[(w,\eta);(v,\mu)]:=a(w;v)+\left<\gamma[Bw-\zeta/\gamma]_{+},Bv-\mu/\gamma\right>-\left<\gamma^{-1}\zeta,\mu\right>.$
(4.27)

Note that in general ($H\neq H^{\prime}$) and it is not possible to prove
well-posedness of (4.26) in continuous spaces. Nevertheless also in this case
a sufficiently smooth solution of the original continuous problem will also be
solution to the formulation (4.26), showing that the formulation remains
consistent.

First we note that for smooth solutions $\lambda\in K$ and (4.21) are
equivalent to (4.23). Then evaluating (4.26) at a sufficiently smooth exact
solution $(u,\lambda)$ we see that for all $(v,\mu)\in V_{h}\times
H_{h}^{\prime}$

$A[(u,\lambda);(v,\mu)]:=a(u;v)+\underbrace{\left<\gamma[Bu-\lambda/\gamma]_{+},Bv-\mu/\gamma\right>}_{=-\left<\lambda,Bv-\mu/\gamma\right>\,\mbox{by
\eqref{eq:lambda0}}}-\left<\gamma^{-1}\lambda,\mu\right>=a(u;v)-\left<\lambda,Bv\right>_{H^{\prime},H}$
(4.28)

and hence by (4.20) the formulation (4.26) is consistent for exact solutions
$(u,\lambda)\in V\times H^{\prime}\cap L^{2}$.

To see the effect of the nonlinear formulation for active and non-active
constraints, first assume $[Bw-\zeta/\gamma]_{+}>0$ in (4.26). The constraint
is active and we see that the equation becomes

$a(w;v)-\left<\mu,Bw\right>-\left<\eta,Bv\right>+\left<\gamma Bw,Bv\right>=0$

which we recognise as the augmented Lagrangian form from (4.17) imposing the
equality constraint $Bw=0$. If on the other hand $[Bw-\eta/\gamma]_{+}=0$ then
the constraint is not active and the equation (4.26) takes the form

$a(w;v)-\left<\gamma^{-1}\eta,\mu\right>=0$

and we see that $Bw$ is free and $\eta=0$ is imposed. As expected the
formulation expresses the conditions of (4.22) and acts as a nonlinear switch
between imposing either $Bu=0$ and $\lambda=0$.

Using the parameter $\gamma$ introduced in (4.16) and the h-weighted norms
introduced in (4.14) and (4.15) together with the inequality
$|[a]_{+}-[b]_{+}|\leq|a-b|$ [32] we see that the following continuity holds

$\left<\gamma([Bw_{1}-\eta_{1}/\gamma]_{+}-[Bw_{2}-\eta_{2}/\gamma]_{+},Bv+\mu/\gamma\right>\\\
\lesssim(\|B(w_{1}-w_{2})\|_{H_{h}}+\|\eta_{1}-\eta_{2}\|_{H_{h}^{\prime}})(\|Bv\|_{H_{h}}+\|\mu\|_{H_{h}^{\prime}})$
(4.29)

Together with (4.12) this shows that the form $A$ is continuous. If $H\equiv
L^{2}$, the formulation (4.26) and (4.29) makes sense on the continuous level.
Observe that unless $r=0$ the norms are $h$ dependent and hence the bound
degenerates for decreasing $h$.

####  4.3.1 Stability, existence and uniqueness of solutions

We will now show that thanks to the properties (4.10) - (4.12) we can derive a
priori bounds on $(w,\eta)$ that allows us to prove existence of a solution in
the spaces $V_{h}\times H_{h}^{\prime}$, using fixed point arguments.

###### Proposition 1

Assume that (4.4)-(4.6) and (4.10)-(4.12) hold. Then for every fixed $h$ the
formulation (4.26) admits a unique solution $(u_{h},\lambda_{h})\in
V_{h}\times H_{h}^{\prime}$. The solution satisfies the a priori bound

$\boxed{\|u_{h}\|_{V}+\gamma_{0}^{\frac{1}{2}}\|[Bu_{h}-\gamma^{-1}\lambda_{h}]_{+}+\gamma^{-1}\lambda_{h}\|_{H_{h}}+\gamma_{0}^{-\frac{1}{2}}\|\lambda_{h}\|_{H_{h}^{\prime}}\lesssim\|f\|_{V^{\prime}}}$
(4.30)

Proof. If we can show that the operator $A$ is continuous and satisfies a
stability condition then existence follows using Brouwer’s fixed point theorem
and the arguments of [83, Chapter 2, Theorem 4.3] (see also [25, Proposition
4.3] for a discussion of finite element methods and augmented Lagrangian
methods). First note that continuity of $A$ follows by (4.29) and (4.12).
Since $h$ is fixed there is no need for the constant of the continuity to be
independent of $h$. Existence of discrete solutions follow from the stability
estimate, for all $w,\eta\in V_{h}\times H_{h}^{\prime}$,

$\boxed{A[(w,\eta),(w+\alpha_{\xi}\xi,-\eta)]\geq\frac{1}{2}\alpha\|w\|^{2}_{V}+\frac{1}{2}\gamma_{0}\|[Bw-\gamma^{-1}\eta]_{+}+\gamma^{-1}\eta\|_{H_{h}}^{2}+\frac{1}{2}\gamma_{0}^{-1}\alpha_{\xi}\|\eta\|_{H_{h}^{\prime}}^{2}}$
(4.31)

where $\xi\in V_{h}$ is a function such that

$\left<B\xi(\eta),q_{h}\right>=-\left<\eta/\gamma,q_{h}\right>,\mbox{ for all
}q_{h}\in H_{h}^{\prime}\mbox{ and
}\|\xi\|_{V}\lesssim\gamma_{0}^{-1}\|\eta\|_{H_{h}^{\prime}}$ (4.32)

(c.f (4.6)), $\gamma_{0}\geq 1$ and
$\alpha_{\xi}=1/2\min(C_{4.6}^{-2},C_{4.12}^{-2}\min(1,\gamma_{0}\alpha)$
where $C_{4.6}$ and $C_{4.12}$ are the constants in the bounds (4.6) and
(4.12) respectively. The bound (4.30) follows from (4.31) since for a solution
$(u_{h},\lambda_{h})$ there holds

$A[(u_{h},\lambda_{h}),(u_{h}+\alpha_{\xi}\xi(\lambda_{h}),-\lambda_{h})]=\left<f,u_{h}+\alpha_{\xi}\xi(\lambda_{h})\right>_{V,V^{\prime}}$

Using the duality pairing we see that

$\left<f,u_{h}+\alpha_{\xi}\xi(\lambda_{h})\right>_{V,V^{\prime}}\leq\|f\|_{V^{\prime}}(\|u_{h}\|_{V}+\alpha_{\xi}\|\xi(\lambda_{h})\|_{V})\lesssim\|f\|_{V^{\prime}}(\|u_{h}\|_{V}+\alpha_{\xi}\gamma_{0}^{-1}\|\lambda_{h}\|_{H_{h}^{\prime}})$

and the claim follows.

To show (4.31) observe that by testing with $(v,\mu)=(w,-\eta)$ we have

$A[(w,\eta),(w,-\eta)]=a(w;w)+\gamma_{0}^{-1}\|\eta\|_{H_{h}^{\prime}}^{2}+\left<\gamma[Bw-\eta/\gamma]_{+},Bw+\eta/\gamma\right>$
(4.33)

By completing the square we see that

$\gamma^{-1}\|\eta\|_{L^{2}}^{2}+\left<\gamma[Bw-\eta/\gamma]_{+},Bw+\eta/\gamma\right>=\gamma_{0}\|[Bw-\eta/\gamma]_{+}+\gamma^{-1}\eta\|_{H_{h}}^{2}$
(4.34)

We conclude that $A$ satisfies the following positivity property, for all
$(w,\eta)\in V_{h}\times H_{h}^{\prime}$,

$A[(w,\eta),(w,-\eta)]=a(w;w)+\gamma_{0}\|[Bw-\eta/(2\gamma)]_{+}+\gamma^{-1}\eta\|_{H_{h}}^{2}$
(4.35)

Then, since $\eta\in H_{h}^{\prime}$ we can use (4.6) to choose $\xi(\eta)\in
V_{h}$ satisfying (4.32), and test with $v=\xi$ and $\mu=0$ to obtain

$\displaystyle A[(w,\eta),(\xi,0)]$
$\displaystyle=a(w;\xi)+\gamma\left<[Bw-\eta/\gamma]_{+},B\xi(\eta)\right>$
(4.36)

Now observe that

$\gamma\left<[Bw-\eta/\gamma]_{+},B\xi(\eta)\right>=\gamma\left<[Bw-\eta/\gamma]_{+}+\gamma^{-1}\eta,B\xi(\eta)\right>-\underbrace{\left<\eta,B\xi(\eta)\right>}_{=-\gamma_{0}^{-1}\|\eta\|_{H_{h}^{\prime}}^{2}}\\\
\geq\gamma_{0}^{-1}\|\eta\|_{H_{h}^{\prime}}^{2}-\frac{1}{2}C_{4.6}^{2}\gamma_{0}\|[Bw+\eta/\gamma]_{+}+\gamma^{-1}\eta\|_{H_{h}}^{2}-\frac{1}{2}\gamma_{0}C_{4.6}^{-2}\|B\xi(\eta)\|_{H_{h}}^{2}$
(4.37)

and since $\gamma_{0}^{\frac{1}{2}}\|B\xi(\eta)\|_{H_{h}}\leq
C_{4.6}\gamma_{0}^{\frac{1}{2}}\|\eta/\gamma\|_{H_{h}}=C_{4.6}\gamma_{0}^{-\frac{1}{2}}\|\eta\|_{H_{h}^{\prime}}$
we see that

$\gamma\left<[Bw-\eta/\gamma]_{+},B\xi(\eta)\right>\geq\frac{1}{2}\gamma_{0}^{-1}\|\eta\|_{H_{h}^{\prime}}^{2}-\frac{1}{2}C_{4.6}^{2}\gamma_{0}\|[Bw+\eta/\gamma]_{+}-\gamma^{-1}\eta\|_{H_{h}}^{2}$

Combining (4.6) with (4.12) we see that using the boundedness $a(w;\xi)\leq
C_{4.12}\|w\|_{V}\|\eta/\gamma\|_{H_{h}}\leq
C_{4.12}\gamma_{0}^{-1}\|w\|_{V}\|\eta\|_{H_{h}^{\prime}}$

$a(w;\alpha_{\xi}\xi)+\gamma\left<[Bw-\eta/\gamma]_{+}+\gamma^{-1}\eta,\alpha_{\xi}B\xi(\eta)\right>\geq-\gamma_{0}^{-1}\alpha_{\xi}C_{4.12}^{2}\|w\|_{V}^{2}\\\
-\alpha_{\xi}C_{4.6}^{2}\gamma_{0}\|[Bw+\eta/\gamma]_{+}+\gamma^{-1}\eta\|_{H_{h}}^{2}+\frac{1}{2}\gamma_{0}^{-1}\alpha_{\xi}\|\eta\|_{H_{h}^{\prime}}^{2}$
(4.38)

The desired inequality then follow by adding (4.35) and (4.38) for
$\gamma_{0}\geq 1$ and

$\alpha_{\xi}=1/2\min(C_{4.6}^{-2},C_{4.12}^{-2})\min(1,\gamma_{0}\alpha).$

If $H\equiv L^{2}$ the analysis can be extended to the continuous case, for
details see [83, Chapter 1, Lemma 4.3].

Uniqueness follows in principle from [83, Chapter 2, Theorem 2.2], but for
completeness we give a simple proof below. Considering the nonlinearity
expressing the constraint we have using the monotonicity
$([a]_{+}-[b]_{+})(a-b)\geq([a]_{+}-[b]_{+})^{2}$, and setting,
$e=w_{1}-w_{2}$ and $\zeta=\eta_{1}-\eta_{2}$,

$\displaystyle\left<\gamma([Bw_{1}-\eta_{1}/\gamma]_{+}-[Bw_{2}-\eta_{2}/\gamma]_{+},Be+\zeta\gamma\right>+\left<\gamma^{-1}(\eta_{1}-\eta_{2}),\zeta\right>$
$\displaystyle\qquad=\left<\gamma([Bw_{1}-\eta_{1}/\gamma]_{+}-[Bw_{2}-\eta_{2}/\gamma]_{+},Be-\zeta/\gamma\right>$
$\displaystyle\qquad\qquad+2\left<\gamma([Bw_{1}-\eta_{1}/\gamma]_{+}-[Bw_{2}-\eta_{2}/\gamma]_{+},\zeta/\gamma\right>+\gamma_{0}^{-1}\|\zeta\|^{2}_{H_{h}^{\prime}}$
$\displaystyle\qquad\geq\gamma_{0}\|[Bw_{1}-\eta_{1}/\gamma]_{+}-[Bw_{2}-\eta_{2}/\gamma]_{+}+\gamma^{-1}\zeta\|_{H_{h}}^{2}.$
(4.39)

It follows from (4.11) and (4.39) that

$\boxed{E_{C}[(w_{1},\eta_{1}),(w_{2},\eta_{2})]^{2}+\alpha\|e\|_{V}^{2}\leq
A[(w_{1},\eta_{1}),(e,-\zeta)]-A[(w_{2},\eta_{2}),(e,-\zeta)]}$ (4.40)

where $E_{C}$ is the error in the approximation of the contact zone defined by

$E_{C}[(w_{1},\eta_{1}),(w_{2},\eta_{2})]:=\gamma_{0}^{\frac{1}{2}}\|[Bw_{1}+\eta_{1}/\gamma]_{+}-[Bw_{2}+\eta_{2}/\gamma]_{+}+\gamma^{-1}\zeta\|_{H_{h}}$

If we assume that both $\\{w_{1},\eta_{1}\\}$ and $\\{w_{2},\eta_{2}\\}$ are
solutions to (4.26) it follows that the right hand side of (4.40) is zero and

$E_{C}[(w_{1},\eta_{1}),(w_{2},\eta_{2})]^{2}+\alpha\|e\|_{V}^{2}=0.$

It follows that $e=0$ and the primal solution is unique. To see that also the
multiplier is unique once again choose $\xi(\eta)$ such that
$B\xi(\eta)=-\eta/\gamma$, in the sense that
$\left<B\xi(\zeta),q_{h}\right>=-\left<\zeta/\gamma,q_{h}\right>$, for all
$q_{h}\in H_{h}^{\prime}$, and test with $v=\xi$ and $\mu=0$, and use
arguments similar as those leading to (4.31) to see that

$\displaystyle 0=$ $\displaystyle
A[(w_{1},\eta_{1}),(\xi(\zeta),0)]-A[(w_{2},\eta_{1}),(\xi(\zeta),0)]$
$\displaystyle\qquad\geq-C^{2}\gamma_{0}^{-1}\underbrace{\|e\|_{V}^{2}}_{I_{1}}+\frac{1}{2}\gamma_{0}^{-1}\|\zeta\|^{2}_{H_{h}^{\prime}}-\gamma_{0}C^{2}\underbrace{E_{C}[(w_{1},\eta_{1}),(w_{2},\eta_{2})]^{2}}_{I_{2}}$
(4.41)

We have already shown in (4.40) that $I_{1}=I_{2}=0$ if both
$\\{w_{1},\eta_{1}\\}$ and $\\{w_{2},\eta_{2}\\}$ are solutions, hence we
conclude that $\|\zeta\|_{H_{h}^{\prime}}=0$ which finishes the discussion of
(discrete) well-posedness.

####  4.3.2 Best approximation results

In this section we will derive a best approximation result for the solution of
(4.26). Due to the nonconforming character of the ALM we need to assume that
the multiplier is in $H^{\prime}\cap L^{2}$. By specifying the approximation
properties of our finite element spaces optimal a priori error estimates can
be deduced.

###### Proposition 2

Assume that (4.4)-(4.6) and (4.10)-(4.12) hold. Let $(u,\lambda)\in
V\times(H^{\prime}\cap L^{2})$ be the solution to (4.19)-(4.22) and
$(u_{h},\lambda_{h})\in V_{h}\times H_{h}^{\prime}$ be the solution of (4.26).
Then if
$\Phi[(u,\lambda),(u_{h},\lambda_{h})]:=E_{C}[(u,\lambda),(u_{h},\lambda_{h})]+\|u-u_{h}\|_{V}+\gamma_{0}^{-\frac{1}{2}}\|\lambda-\lambda_{h}\|_{H_{h}^{\prime}}$
then there holds

$\boxed{\Phi[(u,\lambda),(u_{h},\lambda_{h})]\lesssim\inf_{(v_{h},\mu_{h})\in
V_{h}\times
H_{h}^{\prime}}(\|u-v_{h}\|_{V}+\gamma_{0}^{\frac{1}{2}}\|B(u-v_{h})\|_{H_{h}}+\gamma_{0}^{-\frac{1}{2}}\|\lambda-\mu_{h}\|_{H_{h}^{\prime}})}$
(4.42)

Proof. Since (4.40) holds for all $w_{1},w_{2}\in V$ and $\zeta\in
H^{\prime}\cap L^{2}$, if the exact solution $u,\lambda\in V\times
H^{\prime}\cap L^{2}$ we may apply it with $w_{1}=u$, $\eta_{1}=\lambda$ and
$w_{2}=u_{h}$, $\eta_{2}=\lambda_{h}$ to obtain, with $e=u-u_{h}$ and
$\zeta=\lambda-\lambda_{h}$,

$E_{C}[(u,\lambda),(u_{h},\lambda_{h})]^{2}+\alpha\|e\|_{V}^{2}\leq
A[(u,\lambda),(e,-\zeta)]-A[(u_{h},\lambda_{h}),(e,-\zeta)]$ (4.43)

Using the consistency of the method we have

$A[(u,\lambda),(e,-\zeta)]-A[(u_{h},\lambda_{h}),(e,-\zeta)]=A[(u,\lambda),(u-i_{F}u,\pi_{H}\lambda-\lambda)]-A[(u_{h},\lambda_{h}),(u-i_{F}u,\pi_{H}\lambda-\lambda)]$

By the continuity of $a$ we have

$a(u;u-i_{F}u)-a(u_{h},u-i_{F}u)\leq C\|e\|_{V}\|u-i_{F}u\|_{V}$

For the nonlinearity imposing the constraint we notice that by the
$L^{2}$-orthogonality of $\pi_{H}$,

$\gamma^{-1}\left<\zeta,\lambda-\pi_{H}\lambda\right>=\gamma_{0}^{-1}\|\lambda-\pi_{H}\lambda\|_{H_{h}^{\prime}}^{2}$

and using in addition the properties of $i_{F}u$ we have
$\left<\pi_{H}\zeta,B(u-i_{F}u)+(\lambda-\pi_{H}\lambda)/\gamma\right>=0$ and
hence using that
$\pi_{H}\zeta=\zeta+\pi_{H}\zeta-\zeta=\zeta-(\lambda-\pi_{H}\lambda)$,

$\left<\gamma([Bu+\lambda/\gamma]_{+}-[Bu_{h}+\lambda_{h}/\gamma]_{+}),B(u-i_{F}u)+(\lambda-\pi_{H}\lambda)/\gamma\right>\\\
=\left<\gamma([Bu+\lambda/\gamma]_{+}-[Bu_{h}+\lambda_{h}/\gamma]_{+})+\zeta,B(u-i_{F}u)+(\lambda-\pi_{H}\lambda)/\gamma\right>\\\
-\left<\lambda-\pi_{H}\lambda,B(u-i_{F}u)+(\lambda-\pi_{H}\lambda)/\gamma\right>$
(4.44)

Collecting the above inequalities we obtain using the Cauchy-Schwarz
inequality and the arithmetic-geometric inequality in each right hand side,

$A[(u,\lambda),(e,-\zeta)]-A[(u_{h},\lambda_{h}),(e,-\zeta)]\leq\frac{1}{2}(E_{C}[(u,\lambda),(u_{h},\lambda_{h})]^{2}+\alpha\|e\|_{V}^{2})\\\
+\frac{C}{\alpha}(\|u-i_{F}u\|_{V}^{2}+\gamma_{0}\|B(u-i_{F}u)\|_{H_{h}}^{2}+\gamma_{0}^{-1}\|\lambda-\pi_{H}\lambda\|_{H_{h}^{\prime}}^{2})$

It follows that the following error bound holds,

$E_{C}[(u,\lambda),(u_{h},\lambda_{h})]^{2}+\alpha\|e\|_{V}^{2}\lesssim\alpha\|u-i_{F}u\|_{V}^{2}+\gamma_{0}\|B(u-i_{F}u)\|^{2}_{H_{h}}+\gamma_{0}^{-1}\|\lambda-\pi_{H}\lambda\|_{H_{h}^{\prime}}^{2}$

By adding and subtracting $v_{h}$, applying the triangle inequality followed
by the stability of the Fortin operator (right inequality of (4.5)) there
holds

$\|u-i_{F}u\|_{V}+\gamma_{0}^{\frac{1}{2}}\|B(u-i_{F}u)\|_{H_{h}}\lesssim\|u-v_{h}\|_{V}+\gamma_{0}^{\frac{1}{2}}\|B(u-v_{h})\|_{H_{h}}$

and we conclude using also the definition of the $L^{2}$-projection $\pi_{H}$,
that

$\boxed{E_{C}[(u,\lambda),(u_{h},\lambda_{h})]+\|e\|_{V}\lesssim\inf_{(v_{h},\mu_{h})\in
V_{h}\times
H_{h}^{\prime}}(\|u-v_{h}\|_{V}+\gamma_{0}^{\frac{1}{2}}\|B(u-v_{h})\|_{H_{h}}+\gamma_{0}^{-1}\|\lambda-\mu_{h}\|_{H_{h}^{\prime}})}$
(4.45)

Turning to the error in the multiplier we have using (4.6)

$\gamma_{0}^{-1}\|\pi_{H}\zeta\|_{H_{h}^{\prime}}^{2}=\left<\pi_{H}\zeta,B\xi(\gamma^{-1}\pi_{H}\zeta)\right>$

where $\xi(\gamma^{-1}\pi_{H}\zeta)$ is defined by (4.6) with
$z_{h}=\gamma^{-1}\pi_{H}\zeta$ using the equation we see that

$\gamma_{0}^{-1}\|\pi_{H}\zeta\|_{H_{h}^{\prime}}^{2}=\left<\pi_{H}\lambda-\lambda,B\xi(\pi_{H}\zeta)\right>+\left<\gamma([Bu+\lambda/\gamma]_{+}-[Bu_{h}+\lambda`-h/\gamma]_{+})+\zeta,B\xi(\pi_{H}\zeta)\right>\\\
+a(u;\xi(\pi_{H}\zeta))-a(u_{h};\xi(\pi_{H}\zeta))$

Applying the bound (4.12) to the last two terms of the right hand side and the
Cauchy-Schwarz inequality to the others and applying the stability of (4.6) we
see that

$a(u;\xi(\pi_{H}\eta))-a(u_{h};\xi(\pi_{H}\eta))\leq
C_{4.12}\|e\|_{V}\|\xi(\pi_{H}\eta)\|_{V}\leq
C_{4.12}\|e\|_{V}\gamma_{0}^{-\frac{1}{2}}\|\pi_{H}\zeta\|_{H_{h}^{\prime}},$
$\left<\pi_{H}\lambda-\lambda,B\xi(\pi_{H}\zeta)\right>\leq\gamma_{0}^{-\frac{1}{2}}\|\pi_{H}\lambda-\lambda\|_{H_{h}^{\prime}}\gamma_{0}^{\frac{1}{2}}\|B\xi(\pi_{H}\zeta)\|_{H_{h}}\leq\gamma_{0}^{-\frac{1}{2}}\|\pi_{H}\lambda-\lambda\|_{H_{h}^{\prime}}\gamma_{0}^{-\frac{1}{2}}\|\pi_{H}\zeta\|_{H_{h}^{\prime}}$

and

$\left<\gamma([Bu+\lambda/\gamma]_{+}-[Bu_{h}+\lambda_{h}/\gamma]_{+})+\eta,B\xi(\pi_{H}\eta)\right>\leq
E_{C}[(u,\lambda),(u_{h},\lambda_{h})]\gamma_{0}^{-\frac{1}{2}}\|\pi_{H}\zeta\|_{H_{h}^{\prime}}.$

Collecting terms and dividing through by
$\gamma_{0}^{-\frac{1}{2}}\|\pi_{H}\zeta\|_{H_{h}^{\prime}}$ we have

$\gamma_{0}^{-\frac{1}{2}}\|\pi_{H}\zeta\|_{H_{h}^{\prime}}\lesssim
E_{C}[(u,\lambda),(u_{h},\lambda_{h})]+\|e\|_{V}+\gamma_{0}^{-\frac{1}{2}}\|\pi_{H}\lambda-\lambda\|_{H_{h}^{\prime}}$

We conclude by applying (4.42) to the right hand side and the triangle
inequality $\|\zeta\|\leq\|\lambda-\pi_{H}\lambda\|+\|\pi_{H}\zeta\|$ to
obtain,

$\boxed{\gamma_{0}^{-\frac{1}{2}}\|\zeta\|_{H_{h}^{\prime}}\lesssim\inf_{(v_{h},\mu_{h})\in
V_{h}\times
H_{h}^{\prime}}(\|u-v_{h}\|_{V}+\gamma_{0}^{\frac{1}{2}}\|B(u-v_{h})\|_{H_{h}})+\gamma_{0}^{-\frac{1}{2}}\|\lambda-\mu_{h}\|_{H_{h}^{\prime}})}$
(4.46)

The claim now follows by combining (4.42) and (4.45).
We observe that the natural norm for $\lambda$ here would be $H^{\prime}$, but
that we here consider the corresponding weighted $L^{2}$-norm $H_{h}^{\prime}$
instead. Since this is an $h$-weighted norm, the resulting $L^{2}$ error
estimate is subotimal compared to approximation. Recovering control of the
error in the $H^{\prime}$ norm would require an additional duality argument
that is beyond the scope of this work.

####  4.3.3 Remark on stabilized methods

If the discrete spaces $V_{h}$, $H_{h}^{\prime}$ do not satisfy the infsup
condition (4.5), one can introduce a stabilization operator $s(\cdot,\cdot)$
which is designed to control the unstable modes. If a stable pair $V_{h}$,
$\tilde{H}_{h}^{\prime}$, where $\tilde{H}_{h}^{\prime}$ has the same
approximation properties as $H_{h}^{\prime}$ up to a constant factor, is
known, i.e. (4.5) and (4.6) are satisfied for these spaces, then a convenient
way of choosing $s$ is by using the following design criteria

  1. 1.

Control of unstable modes:

$\gamma_{0}^{-1/2}\|\mu-\tilde{\pi}_{H}\mu\|_{H_{h}^{\prime}}\lesssim
s(\mu,\mu)^{\frac{1}{2}},\quad\forall\mu_{h}\in H_{h}^{\prime}+L^{2}$ (4.47)

where $\tilde{\pi}_{H}$ denotes the $L^{2}$ projection on $\tilde{H}_{h}$.

  2. 2.

Weak consistency:

$s(\mu-\tilde{\pi}_{H}\mu,\mu-\tilde{\pi}_{H}\mu)\sim\gamma_{0}^{-1}\|\mu-\tilde{\pi}_{H}\mu\|_{H_{h}^{\prime}}^{2},\quad\forall\mu\in
L^{2}$ (4.48)

Here the $\sim$ notation means that the two quantities have the same
asymptotics in $h$ for smooth enough $\mu$.

The simplest choice of $s$ is

$s(\eta,\mu)=\gamma^{-1}\left<(\pi_{H}-\tilde{\pi}_{H})\eta,\mu\right>$

The optimality system of the finite element formulation then reads: find
$(u_{h},\lambda_{h})\in V_{h}\times H_{h}^{\prime}$ such that

$A[(u_{h},\lambda_{h});(v,\mu)]-s(\lambda_{h},\mu)=\left<f,v\right>_{V^{\prime},V}$
(4.49)

for all $(v,\mu)\in V_{h}\times H_{h}^{\prime}$, with $A$ defined in (4.28).

It is then possible to use the monotonicity, the inf-sup stability (4.5)
together with (4.47) and (4.48) to obtain bounds similar to (4.42) for the
error of the stabilized Galerkin approximation. We only sketch the arguments.
The only modification of the stability is that the stabilization operator
appears in the left hand side. If $e=u-u_{h}$ and $\zeta=\lambda-\lambda_{h}$
then

$E_{C}[(u,\lambda),(u_{h},\lambda_{h})]^{2}+\alpha\|e\|_{V}^{2}+s(\zeta,\zeta)\leq
A[(u,\lambda),(e,-\zeta)]-A[(u_{h},\lambda_{h}),(e,-\zeta)]-s(\zeta,-\zeta)$
(4.50)

The key observation to obtain optimal approximation is to use Galerkin
orthogonality using $u_{h}-i_{F}u$ and $\lambda_{h}-\tilde{\pi}_{H}\lambda$
and then apply a modified continuity estimate. Indeed by the assumptions we
have $\left<\tilde{\pi}_{H}\zeta,B(u-i_{F}u)\right>=0$ and hence we can modify
the continuity (4.44) the following way,

$\left<\gamma([Bu+\lambda/\gamma]_{+}-[Bu_{h}+\lambda_{h}/\gamma]_{+}),B(u-i_{F}u)+(\lambda-\tilde{\pi}_{H}\lambda)/\gamma\right>\\\
=\left<\gamma([Bu+\lambda/\gamma]_{+}-[Bu_{h}+\lambda_{h}/\gamma]_{+})+\tilde{\pi}_{H}\zeta,B(u-i_{F}u)+(\lambda-\tilde{\pi}_{H}\lambda)/\gamma\right>\\\
=\left<\gamma([Bu+\lambda/\gamma]_{+}-[Bu_{h}+\lambda_{h}/\gamma]_{+})+\zeta,B(u-i_{F}u)+(\lambda-\tilde{\pi}_{H}\lambda)/\gamma\right>\\\
-\left<\zeta-\tilde{\pi}_{H}\zeta,B(u-i_{F}u)+(\lambda-\tilde{\pi}_{H}\lambda)/\gamma\right>$

where we used that
$\left<\tilde{\pi}_{H}\zeta,B(u-i_{F}u)+(\lambda-\tilde{\pi}_{H}\lambda)/\gamma\right>=0$.
In this expression all but the last term can be bounded in the same fashion as
before. For the last term we apply the Cauchy-Schwarz inequality and then
(4.47) to see that

$\left<\zeta-\tilde{\pi}_{H}\zeta,B(u-i_{F}u)+(\lambda-\tilde{\pi}_{H}\lambda)/\gamma\right>\leq
s(\zeta,\zeta)^{\frac{1}{2}}(\gamma_{0}^{\frac{1}{2}}\|B(u-i_{F}u)\|_{H_{h}}+\gamma_{0}^{-\frac{1}{2}}\|\lambda-\tilde{\pi}_{H}\lambda\|_{H_{h}^{\prime}})$

where now the right hand side is controlled by stability and approximation
respectively. This leads to an error estimate for $u-u_{h}$. The error in the
multiplier can also be estimated using that

$\|\zeta\|_{H^{\prime}_{h}}\leq\|\tilde{\pi}_{h}(\zeta-\zeta_{h})\|_{H^{\prime}_{h}}+s(\zeta,\zeta)^{\frac{1}{2}}$

and noting that the first term of the right hand side can be controlled as in
the infsup stable case and the second is bounded by (4.47).

###  4.4 Eliminating the multiplier

Now we assume that the multiplier can be expressed in the primal variable
through a linear operator $T$ on the continuous level, i.e. $\lambda=Tu$, such
that for $v_{h}\in V_{h}$, the following inequality that typically is of
inverse type, holds

$\|Tu_{h}\|_{H_{h}^{\prime}}^{2}\leq C_{\text{I}}\|u_{h}\|^{2}_{V}.$ (4.51)

where $C_{\text{I}}$ is a constant that may depend on the mesh geometry, but
not on the mesh size. We may then write the Nitsche type form of the equation
(4.26): find $u_{h}\in V_{h}$ such that

$A[(u_{h},Tu_{h});(v,Tv)]=\left<f,v\right>_{V^{\prime},V}$ (4.52)

for all $v\in V_{h}$, where $A_{h}$ was defined in (4.28). This formulation,
where the multiplier is eliminated is identified as a nonlinear GLS method.
For this GLS formulation existence and uniqueness is ensured without any inf-
sup condition [32, Theorem 3.3]. Stability is obtained thanks to the
continuity of the $T$ operator, (4.51).

We now revisit the analysis of the previous section and show that the same
results hold for the case when the multiplier has been eliminated.

####  4.4.1 Continuity and stability

We only need to verify (4.29) for the method (4.52). We immediately have for
$w_{1},w_{2},v\in V_{h}$,

$\left<\gamma([Bw_{1}-Tw_{1}/\gamma]_{+}-[Bw_{2}-Tw_{2}/\gamma]_{+},Bv+Tv/\gamma\right>\\\
\lesssim(\|B(w_{1}-w_{2})\|_{H_{h}}+\|T(w_{1}-w_{2})\|_{H_{h}^{\prime}})(\|Bv\|_{H_{h}}+\|Tv\|_{H_{h}^{\prime}})\\\
\leq
C(\|B(w_{1}-w_{2})\|_{H_{h}}+\|w_{1}-w_{2}\|_{V})(\|Bv\|_{H_{h}}+\|v\|_{V})$
(4.53)

where we used (4.51) for the second inequality. To prove the a priori estimate
that together with the continuity allows for the fixed point analysis we test
with $v=u_{h}$ in (4.52) to obtain using (4.10)

$\alpha\|u_{h}\|_{V}^{2}+\|[Bu_{h}-Tu_{h}/\gamma]_{+}\|^{2}_{H_{h}}-\gamma_{0}^{-1}\|Tu_{h}\|_{H_{h}^{\prime}}^{2}\leq
A[(u_{h},Tu_{h});(u_{h},Tu_{h})]$

Applying (4.51) to the last term of the right hand side we see that

$\boxed{(\alpha-
C_{I}/\gamma_{0})\|u_{h}\|_{V}^{2}+\gamma_{0}\|[Bu_{h}-Tu_{h}/\gamma]_{+}\|^{2}_{H_{h}}\leq
A[(u_{h},Tu_{h});(u_{h},Tu_{h})]}$

We conclude that the stability holds for $\gamma_{0}>C_{I}/\alpha$. Hence
under this condition there exists a discrete solution to (4.52)

####  4.4.2 Uniqueness and best approximation estimates

Uniqueness and best approximation follows using similar arguments, we only
detail the best approximation case. We assume that the exact solution $u$ to
(4.19) is sufficiently smooth that

$A_{h}[(u,Tu);(v,Tv)]=\left<f,v\right>_{V^{\prime},V},\forall v\in V_{h}$
(4.54)

Then we may write, $e=u-u_{h}$ and using the monotonicity of $a$ (4.11) and of
$[\cdot]_{+}$ we see that, using the notation
$E_{C}(u,u_{h}):=\gamma_{0}\|[Bu-
Tu\gamma]_{+}-[Bu_{h}-Tu_{h}/\gamma]_{+}\|_{H_{h}}^{2}$,

$A_{h}[(u,Tu);(e,Te)]-A_{h}[(u_{h},Tu_{h});(e,Te)]\\\
\geq\alpha\|e\|_{V}^{2}+E_{C}(u,u_{h})-\gamma_{0}^{-1}\|Te\|_{H_{h}^{\prime}}^{2}$

For the last term of the right hand side observe that

$\|Te\|_{H_{h}^{\prime}}\leq\|T(u-v_{h})\|_{H_{h}^{\prime}}+\|T(u_{h}-v_{h})\|_{H_{h}^{\prime}}\\\
\leq\|T(u-v_{h})\|_{H_{h}^{\prime}}+C_{I}^{1/2}\|u_{h}-v_{h}\|_{V}\\\
\leq\|T(u-v_{h})\|_{H_{h}^{\prime}}+C_{I}^{1/2}(\|e\|_{V}+\|u-v_{h}\|_{V})$

Hence

$A_{h}[(u,Tu);(e,Te)]-A_{h}[(u_{h},Tu_{h});(e,Te)]\\\
\geq(\alpha-3C_{I}/\gamma_{0})\|e\|_{V}^{2}+E_{C}(u,u_{h})-3\|T(u-v_{h})\|_{H_{h}^{\prime}}^{2}-3C_{I}\|u-v_{h}\|_{V}^{2}$
(4.55)

Fix $\gamma_{0}=6C_{I}/\alpha$ so that $\alpha-3C_{I}/\gamma_{0}=\alpha/2$.
Considering the left hand side we have using (4.54), for all $v_{h}\in V_{h}$

$A_{h}[(u,Tu);(e,Te)]-A_{h}[(u_{h},Tu_{h});(e,Te)]\\\
=A_{h}[(u,Tu);(u-v_{h},T(u-v_{h}))]-A_{h}[(u_{h},Tu_{h});(u-v_{h},T(u-v_{h}))]$

To conclude we use the continuity (4.12) and the arithmetic-geometric
inequality,

$a(u;u-v_{h})-a(u_{h};u-v_{h})\leq
C\|e\|_{V}\|u-v_{h}\|_{V}\leq\frac{\alpha}{4}\|e\|_{V}^{2}+C^{2}\|u-v_{h}\|_{V}^{2}$

together with the Cauchy-Schwarz inequality and the arithmetic-geometric
inequality,

$\left<\gamma([Bu-
Tu\gamma]_{+}-[Bu_{h}-Tu_{h}/\gamma]_{+},B(u-v_{h})+T(u-v_{h})/\gamma\right>\\\
\leq\frac{1}{2}E_{C}(u,u_{h})+\gamma_{0}\|B(u-v_{h})\|_{H_{h}}^{2}+\gamma_{0}^{-1}\|T(u-v_{h})\|_{H_{h}^{\prime}}^{2}.$

Applying these inequalities in (4.55) we see that for all $v_{h}\in V_{h}$

$\alpha\|e\|_{V}^{2}+E_{C}(u,u_{h})\\\
\lesssim\|u-v_{h}\|_{V}^{2}+\|B(u-v_{h})\|_{H_{h}}^{2}+\|T(u-v_{h})\|_{H_{h}^{\prime}}^{2}$

Taking square roots of both sides and the infimum over $v_{h}\in V_{h}$ in the
right hand side we conclude

$\boxed{E_{C}(u,u_{h})+\|e\|_{V}\lesssim\inf_{v_{h}\in
V_{h}}(\|u-v_{h}\|_{V}+\|B(u-v_{h})\|_{H_{h}}+\|T(u-v_{h})\|_{H_{h}^{\prime}})}$
(4.56)

We have sketched a best approximation result for the formulation (4.52).
Observe that no condition needs to be imposed on the finite element space in
this case. Instead stability is ensured by the inverse inequality (4.51) that
bounds the $H_{h}^{\prime}$-norm of the multiplier expressed in the primal
variable by the $V$-norm of the primal variable. By equivalence of norms on
finite dimensional spaces this bound is always true. The key to optimality of
the estimate is the proper $h$-scaling of the discrete norms given in (4.14)
and (4.15).

We now turn to specific examples.

##  5 Applications

###  5.1 The Stokes problem with cavitation

Consider a domain $\Omega$ in ${\mathbb{R}}^{n}$, $n=2$ or $n=3$ with boundary
$\partial\Omega$ that is composed of the two subsets $\Gamma_{D}$ and
$\Gamma_{N}$ such that $\partial\Omega=\bar{\Gamma}_{D}\cup\bar{\Gamma}_{N}$.
We consider a lubricant with viscosity $\mu$. The Stokes equation can then be
written

$-\mu\Delta\boldsymbol{u}+\nabla
p={\boldsymbol{f}}\;\text{and}\;\nabla\cdot\boldsymbol{u}=0\quad\text{in}\;\Omega,$
(5.1)

with $\boldsymbol{u}=0$ on $\Gamma_{D}$ and
$(-p\boldsymbol{I}+\mu\nabla\boldsymbol{u})\cdot\boldsymbol{n}=\boldsymbol{0}$
on $\Gamma_{N}$. Here, $\boldsymbol{u}$ is the velocity of the lubricant, $p$
is the pressure, and ${\boldsymbol{f}}$ is a force term. The lubricant cannot
support subatmospheric pressure, so an additional condition is $p\geq 0$ in
$\Omega$. In order to incorporate this condition into the model, it can be
written as a variational inequality as follows. Let

$a(\boldsymbol{u},\boldsymbol{v}):=\int_{\Omega}\mu\nabla\boldsymbol{u}:\nabla\boldsymbol{v}\,{\rm
d}\Omega,\quad
L(\boldsymbol{v}):=\int_{\Omega}\boldsymbol{f}\cdot\boldsymbol{v}\,{\rm
d}\Omega$

and

$K=\\{p\in L_{2}(\Omega):\quad p\geq 0\\}$

Seek $\boldsymbol{u}\in[H_{0}^{1}(\Omega)]^{n}$ and $p\in K$ such that

$a(\boldsymbol{u},\boldsymbol{v})-\int_{\Omega}p\,\nabla\cdot{\boldsymbol{v}}\,d\Omega=L(\boldsymbol{v}),$
(5.2)

for all ${\boldsymbol{v}}\in[H^{1}(\Omega)]^{n}$, and

$-\int_{\Omega}\nabla\cdot\boldsymbol{u}\,(q-p)\,d\Omega\leq 0,\quad\forall
q\in K$ (5.3)

To rewrite this problem as a variational equality, we use the Kuhn-Tucker
conditions

$p\geq 0,\quad\nabla\cdot\boldsymbol{u}\geq 0,\quad
p\,\nabla\cdot\boldsymbol{u}=0$ (5.4)

and again replace conditions (5.4) by the equivalent statement

$p=\gamma_{0}[\gamma_{0}^{-1}p-\nabla\cdot\boldsymbol{u}]_{+}$ (5.5)

with $\gamma_{0}$ a positive number. We note here that we can identify the
abstract spaces $H$ and $H^{\prime}$ with $L_{2}(\Omega)$ and that here the
pressure cannot easily be interpreted as coming from a linear operator on the
velocity, so we are in cases A and D from Section 4; the pressure has to be
retained but $r=0$ in the discrete norms.

Defining function spaces

$V=\\{\boldsymbol{v}\in[H^{1}(\Omega)]^{n}:\;\boldsymbol{v}=\boldsymbol{0}\;\text{on
$\Gamma^{\text{D}}$}\\},\quad Q=L_{2}(\Omega)$ (5.6)

and seeking $(\boldsymbol{u},p)\in V\times Q$ we seek stationary points to the
functional

$\mathcal{L}_{\text{A}}(\boldsymbol{u},p):=\frac{1}{2}a(\boldsymbol{u},\boldsymbol{u})-L(\boldsymbol{u})+\int_{\Omega}\frac{\gamma_{0}}{2}\left[\gamma_{0}^{-1}p-\nabla\cdot\boldsymbol{u}\right]_{+}^{2}d\Omega-\int_{\Omega}\frac{1}{2\gamma_{0}}p^{2}d\Omega$
(5.7)

analogously to (4.24).

For the discrete problem, we will use the inf–sup stable Taylor-Hood
approximation which utilises the finite element space

$\vec{V}^{h}=\\{\boldsymbol{v}:\boldsymbol{v}\in\left[C^{0}(\Omega)\right]^{d},\;\boldsymbol{v}|_{K}\in[P^{2}(K)]^{d},\
\forall K\in{\mathcal{T}}^{h},\;\text{$\boldsymbol{v}=\boldsymbol{0}$ on
$\Gamma^{\text{D}}$}\\}$

for the velocity, where $P^{2}(K)$ denotes the space of piecewise quadratic
polynomials on $K$, and the space $Q_{h}$ of piecewise linears for the
pressure:

$Q^{h}=\\{p\in C^{0}(\Omega):\;p|_{K}\in P^{1}(K),\,\forall
K\in{\mathcal{T}}^{h}\\}.$ (5.8)

The finite element method based on (5.7) is to find
$(\boldsymbol{u}^{h},p^{h})\in\vec{V}^{h}\times Q^{h}$ such that

$a(\boldsymbol{u}^{h},\boldsymbol{v})-\int_{\Omega}{\gamma_{0}}[\gamma_{0}^{-1}p^{h}-\nabla\cdot\boldsymbol{u}^{h}]_{+}\nabla\cdot\boldsymbol{v}d\Omega=(\boldsymbol{f},\boldsymbol{v})\quad\forall\boldsymbol{v}\in\vec{V}^{h},$
(5.9)

and

$\int_{\Omega}\left({\gamma_{0}}[\gamma_{0}^{-1}p^{h}-\nabla\cdot\boldsymbol{u}^{h}]_{+}-p^{h}\right)qd\Omega=0,\quad\forall
q\in Q^{h}.$ (5.10)

####  5.1.1 Satisfaction of assumptions for the abstract analysis

For the present problem we have $V=[H^{1}(\Omega)]^{n}$,
$H=H_{h}=H_{h}^{\prime}=L^{2}(\Omega)$. The constraint operator $B$ is the
divergence operator. It if well known that the Taylor-Hood element admits a
Fortin interpolant satisfying

$\|\nabla\pi_{F}\boldsymbol{v}\|_{\Omega}\lesssim\|\nabla\boldsymbol{v}\|_{\Omega},\forall\boldsymbol{v}\in[H^{1}(\Omega)]^{d}.$

Since
$\|\nabla\pi_{F}\boldsymbol{v}\|_{\Omega}\lesssim\|\nabla\pi_{F}\boldsymbol{v}\|_{\Omega}$
the relation (4.5) holds. This means that for all $\mu_{h}\in Q_{h}$ there
exists $\boldsymbol{v}_{h}\in\vec{V}^{h}$ such that
$(\nabla\cdot\boldsymbol{v}_{h},q_{h})_{\Omega}=(\mu_{h},q_{h})_{\Omega}$ and
$\|\boldsymbol{v}_{h}\|_{V}\lesssim\|\mu_{h}\|_{\Omega}$. Hence (4.6) is also
satisfied.

Since $a(\cdot,\cdot)$ is a linear operator in this case we see that
(4.10)-(4.12) are satisfied using standard arguments. Hence the assumptions of
section 4 are satisfied in this case and hence we conclude that the best
approximation estimate (4.42) holds.

###  5.2 Weak imposition of Dirichlet boundary conditions

####  5.2.1 Model problem

Let us first consider the Poisson model problem: find
$u:\Omega\rightarrow\mathbb{R}$ such that

$-\Delta u=f~{}\text{in}~{}\Omega,\quad
u=g~{}\text{on}~{}\Gamma:=\partial\Omega$ (5.11)

where $\Omega$ is a bounded domain in two or three space dimensions, with
outward pointing normal $\boldsymbol{n}$, and $f$ and $g$ are given functions.
For simplicity, we shall assume that $\Omega$ is polyhedral (polygonal). A
classical way of prescribing $u=g$ on the boundary is to pose the problem
(5.11) as a minimisation problem with side conditions and seek stationary
points to the functional

$\mathcal{L}(v,\mu):=\frac{1}{2}a(v,v)-\left<\mu,v-g\right>_{H^{1/2}(\Gamma),H^{-1/2}(\Gamma)}-(f,v)_{\Omega}$
(5.12)

where

$(f,v)_{\Omega}:=\int_{\Omega}fv\,d\Omega,\quad a(u,v):=\int_{\Omega}\nabla
u\cdot\nabla v\,d\Omega$ (5.13)

and $\left<\mu,v-g\right>_{H^{1/2}(\Gamma),H^{-1/2}(\Gamma)}$ is interpreted
as a duality pairing on $H^{-1/2}(\Gamma)\times H^{1/2}(\Gamma)$. We are thus
in case B of Sec. 4, and the method proposed will only make sense on discrete
spaces.

The stationary points to (5.12) are given by finding $(u,\lambda)\in
H^{1}(\Omega)\times H^{-1/2}(\Gamma)$ such that

$a(u,v)-\left<\lambda,v\right>_{H^{1/2}(\Gamma),H^{-1/2}(\Gamma)}=(f,v)\quad\forall
v\in H^{1}(\Omega)$ (5.14)
$\left<\mu,u\right>_{H^{1/2}(\Gamma),H^{-1/2}(\Gamma)}=\left<\mu,g\right>_{H^{1/2}(\Gamma),H^{-1/2}(\Gamma)}\quad\forall\mu\in
H^{-1/2}(\Gamma)$ (5.15)

As mentioned above, the discretisation of this problem requires balancing of
the discrete spaces for the multiplier $\lambda$ and the primal solution $u$
in order for the method to be stable.

####  5.2.2 The augmented Lagrangian method for boundary conditions

The Lagrangian in (5.12) is augmented by a penalty term scaled by a parameter
$\gamma\in\mathbb{R}^{+}$ so that we seek stationary points to

$\mathcal{L}_{A}(v,\mu):=\frac{1}{2}a(v,v)-\left<\mu,v-g\right>_{H^{1/2}(\Gamma),H^{-1/2}(\Gamma)}+\frac{\gamma}{2}\|(v-g)\|^{2}_{H^{1/2}(\Gamma)}-(f,v)_{\Omega}$
(5.16)

We note that the continuous norms imply $r=1/2$ in the discrete norms. To find
the stationary points we seek $(u,\lambda)$ such that

$\displaystyle
a(u,v)-\left<\lambda,v\right>_{H^{1/2}(\Gamma),H^{-1/2}(\Gamma)}+\gamma(u-g,v)_{H^{1/2}(\Gamma)}+\left<u,\mu\right>_{H^{1/2}(\Gamma),H^{-1/2}(\Gamma)}={}$
$\displaystyle(f,v)_{\Omega}$ $\displaystyle+\gamma\langle
g,v\rangle_{H^{1/2}(\Gamma)}$
$\displaystyle+\left<g,\mu\right>_{H^{1/2}(\Gamma),H^{-1/2}(\Gamma)}$

To determine the Lagrange multiplier $\lambda$ we set $\mu=0$, and integrate
by parts which gives

$\displaystyle(-\Delta
v+f,v)_{\Omega}+\langle\nabla_{n}u-\lambda,v\rangle_{H^{1/2}(\Gamma),H^{-1/2}(\Gamma)}=0$
(5.17)

For the exact solution the first term vanish and we conclude that
$\lambda=\nabla_{n}v$.

We now wish to find a stable discrete counterpart to this optimisation
problem. To this end, let $\mathcal{T}_{h}$ be a family of quasi–uniform
partitions, with mesh parameter $h$, of $\Omega$ into shape regular triangles
or tetrahedra $T$ and the discrete space

$V_{h}:=\\{v_{h}\in H^{1}(\Omega):v_{h}|_{T}\in\mathbb{P}_{k}(T),\,\forall
T\in{\mathcal{T}_{h}}\\},\quad\mbox{ for }k\geq 1$ (5.18)

and some discrete space $Q_{h}$ (not explicitly defined) for the approximation
of the Lagrange multiplier.

We first follow the idea of (4.14) and replace the $H^{1/2}$–norm by the
discrete counterpart $h^{-1/2}\|\cdot\|_{L_{2}(\Gamma)}$, which by an inverse
estimate dominates the $H^{1/2}(\Gamma)$ norm,

$\|v\|^{2}_{H^{1/2}(\Gamma)}\lesssim
h^{-1}\|v\|^{2}_{L_{2}(\Gamma)}\quad\text{$v\in V_{h}$}$ (5.19)

and introduce the problem of finding the stationary point in $V_{h}\times
Q_{h}$ of the discrete Lagrangian

$\mathcal{L}_{A}^{h}(v,\mu):=\frac{1}{2}a(v,v)-({\mu},{v-g})_{\Gamma}+\frac{\gamma_{0}}{2h}\|v-g\|^{2}_{L_{2}(\Gamma)}-(f,v)_{\Omega}$
(5.20)

Recalling next that formally the Lagrange multiplier in (5.14) is given by
$\mu=\nabla_{n}v$, which provides a direct way of computing the Lagrange
multiplier from the primal solution, we obtain

$\mathcal{L}_{A}^{h}(v):=\frac{1}{2}a(v,v)-(\nabla_{n}{v},{v-g})_{\Gamma}+\frac{\gamma_{0}}{2h}\|v-g\|^{2}_{L_{2}(\Gamma)}-(f,v)_{\Omega}$
(5.21)

This is our stabilised ALM, the minimiser to which solves the problem of
finding $u_{h}\in V_{h}$ such that

$a(u_{h},v)-({\nabla_{n}u_{h}},{v})_{\Gamma}-({\nabla_{n}v},{u_{h}})_{\Gamma}+\gamma_{0}h^{-1}({u_{h}},{v})_{\Gamma}=l(v)\quad\forall
v\in V_{h}$ (5.22)

where

$l(v):=(f,v)+(\gamma_{0}h^{-1}{v}-{\nabla_{n}v},{g})_{\Gamma}$ (5.23)

We identify the classical method of Nitsche [85], stable if $\gamma_{0}$ is
chosen so that $\gamma_{0}>\gamma_{C}$, where $\gamma_{C}$ is the constant in
the inverse inequality

$h\|\nabla_{n}v\|_{L_{2}(\Gamma)}^{2}\leq\gamma_{C}\|\nabla
v\|_{L_{2}(\Omega)}^{2}$ (5.24)

###### Remark 1

As shown by Stenberg [97] (and discussed in Sec. 4.4), Nitsche’s method can be
viewed as a particular instance of the GLS stabilisation method of
Barbosa–Hughes [4]; in this sense the ALM is a variant of GLS, with the
multiplier eliminated.

###### Remark 2

We note that the ALM leads to the symmetric form of Nitsche’s method. The
corresponding unsymmetric forms, as discussed, e.g., in [36], are derived
using different arguments.

###  5.3 Inequality boundary conditions

An important feature of the augmented Lagrangian approach is that it can be
extended to the case of inequality constraints, as first shown by Chouly and
Hild in the context of elastic contact [32]. We consider the problem: find
$u:\Omega\rightarrow\mathbb{R}$ such that

$-\Delta u=f~{}\text{in}~{}\Omega,\quad u-g\leq 0~{}\text{on}~{}\Gamma$ (5.25)

We have the following Kuhn–Tucker conditions on the multiplier and side
condition:

$u-g\leq 0,\quad\lambda\leq 0,\quad\lambda(u-g)=0.$ (5.26)

We now use the analogue to (2.29), that (5.26) is equivalent to

$\lambda=-{\gamma}\,[u-g-\gamma^{-1}\,\lambda]_{+}$ (5.27)

first used in this context by Alart and Curnier [1]. Now we can take another
route to the augmented Lagrangian method. Taking the discrete counterpart to
the standard multiplier equilibrium equation (5.14) we find

$(f,v)=a(u_{h},v)-({\lambda_{h}},{v})_{\Gamma}=a(u_{h},v)-({\lambda_{h}},{v-\gamma^{-1}\mu})_{\Gamma}-({\gamma^{-1}\lambda_{h}},{\mu})_{\Gamma}$
(5.28)

for all $v\in V_{h}$ and $\mu\in Q_{h}$ arbitrary. Using now (5.27) we find

$(f,v)=a(u_{h},v)+({{\gamma}\,[u_{h}-g-\gamma^{-1}\,\lambda_{h}]_{+}},{v-\gamma^{-1}\mu})_{\Gamma}-({\gamma^{-1}\lambda_{h}},{\mu})_{\Gamma}\quad\forall(v,\mu)\in
V_{h}\times Q_{h}.$ (5.29)

This is the optimality system for the Lagrangian

$\mathcal{L}_{A}^{h}(v,\mu):=\frac{1}{2}a(v,v)+\frac{1}{2}\|\gamma^{1/2}[v-g-\gamma^{-1}\mu]_{+}\|_{L_{2}(\Gamma)}^{2}-\|\gamma^{-1/2}\mu\|^{2}_{L_{2}(\Gamma)}-(f,v)_{\Omega}$
(5.30)

cf. [1]. Approximating $\lambda_{h}\approx\partial_{n}u_{h}$ and setting
$\mu=\partial_{n}v$, we seek $u_{h}\in V_{h}$ such that

$a(u_{h},v)+({{\gamma}\,[u_{h}-g-\gamma^{-1}\,\partial_{n}u_{h}]_{+}},{v-\gamma^{-1}\partial_{n}v})_{\Gamma}-({\gamma^{-1}\partial_{n}u_{h}},{\partial_{n}v})_{\Gamma}=(f,v)_{\Omega}\quad\forall\in
V_{h}$ (5.31)

The solution to this problem is the minimiser of the nonlinear augmented
Lagrangian

$\mathcal{L}_{A}^{h}(v):=\frac{1}{2}a(v,v)+\frac{1}{2}\|\gamma^{1/2}[v_{h}-g-\gamma^{-1}\partial_{n}v]_{+}\|_{L_{2}(\Gamma)}^{2}-\|\gamma^{-1/2}\partial_{n}v\|^{2}_{L_{2}(\Gamma)}-(f,v_{h})_{\Omega}$
(5.32)

Again, we choose $\gamma=\gamma_{0}/h$. Variants and several extensions of
(5.31) can be found in [23]. We remark here that (5.31) coincides with (5.22)
in case of contact and gives a penalty on $\partial_{n}u=0$ on $\Gamma$ in
case of no contact. This penalty does not destroy the coercivity of the
problem if (5.19) is satisfied.

###### Remark 3

In the GLS stabilisation for variational inequalities proposed by Barbosa and
Hughes [5], no penalty is added to the Lagrangian; the multiplier is not
eliminated, and their approach is a stabilised Lagrange multiplier method
which requires the solution of an inequality problem. It is also possible to
retain the multiplier in the ALM and add GLS stabilisation to the augmented
Lagrangian. This approach, which also leads to a nonlinear equality problem,
was explored in [67].

####  5.3.1 Satisfaction of assumptions for the abstract analysis

In this case $V=H^{1}(\Omega)$ and $H=H^{\frac{1}{2}}(\partial\Omega)$,
$H^{\prime}=H^{-\frac{1}{2}}(\partial\Omega)$. However since the solution to
(5.25) is known to have the additional regularity $u\in
H^{\frac{3}{2}+\epsilon}(\Omega)$, $\epsilon>0$ it follows that
$\partial_{n}u\in L^{2}(\Omega)$ and the discrete norms $H_{h}$ and
$H_{h}^{\prime}$ defined by (4.14) and (4.15) are well defined on the exact
solution. While (5.24) then is enough to make the formulation (5.31) satisfy
the assumptions necessary for the analysis of section 4.4, the formulation
(5.30) still requires the satisfaction of (4.5) and (4.6). For a
charaterisation of spaces satisfying these conditions (in the $h$-weighted
$L^{2}$-norm) we refer to [87]. An example of a construction is two space
dimension is to take element wise constant approximation for $Q_{h}$ and let
$V_{h}$ consist of piecewise quadratic continuous approximation, or piecewise
affine approximation enriched with a quadratic bubble added to elements
adjacent to the boundary on each boundary face. The Fortin interpolant can
then be constructed by first defining the nodal degrees of freedom using any
$H^{1}$-stable interpolant and then fixing the degree of freedom associated to
the bubble on each boundary faces so that (4.5) and (4.6) are satisfied.
Indeed here they are equivalent. The same construction may be used for the
forthcoming sections.

###  5.4 A model for elastic contact

####  5.4.1 Treatment of Robin boundary conditions

To show the versatility of the ALM we shall consider the equations of linear
elasticity in contact with a springy substrate. We start with the linear case
of a Robin boundary condition: Find the displacement
${\boldsymbol{u}}=\left[u_{i}\right]_{i=1}^{n}$ and the symmetric stress
tensor ${\boldsymbol{\sigma}}=\left[\sigma_{ij}\right]_{i,j=1}^{n}$ such that

$\displaystyle{\boldsymbol{\sigma}}={}$ $\displaystyle\frac{\nu
E}{(1+\nu)(1-2\nu)}~{}\text{tr}\,{\boldsymbol{\varepsilon}}({\boldsymbol{u}})\,{\boldsymbol{I}}+\frac{E}{(1+\nu)}{\boldsymbol{\varepsilon}}({\boldsymbol{u}})\quad\text{in}\quad\Omega,$
(5.33) $\displaystyle-\nabla\cdot{\boldsymbol{\sigma}}={}$
$\displaystyle{\boldsymbol{f}}\quad\text{in}\quad\Omega,$ (5.34)
$\displaystyle{\boldsymbol{S}}{\boldsymbol{u}}={}$
$\displaystyle-{\boldsymbol{\sigma}}\cdot{\boldsymbol{n}}\quad\text{on}\quad\partial\Omega_{\text{S}},$
(5.35) $\displaystyle{\boldsymbol{\sigma}}\cdot{\boldsymbol{n}}={}$
$\displaystyle{\bf
0}\quad\text{on}\quad\partial\Omega\setminus\partial\Omega_{\text{S}}.$ (5.36)

Here $\Omega$ is a closed subset of $\mathbb{R}^{n}$, $n=2$ or $n=3$, $E$ is
Young’s modulus and $\nu$ is Poisson’s ratio.
${\boldsymbol{\varepsilon}}\left({\boldsymbol{u}}\right)=\left[\varepsilon_{ij}({\boldsymbol{u}})\right]_{i,j=1}^{n}$
is the strain tensor with components

$\varepsilon_{ij}({\boldsymbol{u}})=\frac{1}{2}\left(\frac{\partial
u_{i}}{\partial x_{j}}+\frac{\partial u_{j}}{\partial x_{i}}\right),$

and trace

$\text{tr}\,{\boldsymbol{\varepsilon}}({\boldsymbol{u}})=\sum_{i}\varepsilon_{ii}({\boldsymbol{u}})=\nabla\cdot{\boldsymbol{u}}.$

Furthermore,
$\nabla\cdot{\boldsymbol{\sigma}}=\left[\sum_{j=1}^{n}\partial\sigma_{ij}/\partial
x_{j}\right]_{i=1}^{n}$,
${\boldsymbol{I}}=\left[\delta_{ij}\right]_{i,j=1}^{n}$ with $\delta_{ij}=1$
if $i=j$ and $\delta_{ij}=0$ if $i\neq j$, and ${\boldsymbol{f}}$ is a given
load. Finally, we assume that the boundary stiffness ${\boldsymbol{S}}$ is of
the form

${\boldsymbol{S}}=\alpha^{-1}{\boldsymbol{n}}\otimes{\boldsymbol{n}}+\beta^{-1}{\boldsymbol{P}},\quad{\boldsymbol{P}}:=({\boldsymbol{I}}-{\boldsymbol{n}}\otimes{\boldsymbol{n}})$

where $\alpha$ and $\beta$ are flexibility parameters in the normal and
tangential direction, respectively. The solution to (5.33)–(5.36) minimises
the functional

$\mathcal{L}_{S}({\boldsymbol{u}}):=\frac{1}{2}a({\boldsymbol{u}},{\boldsymbol{u}})-({\boldsymbol{f}},{\boldsymbol{u}})_{\Omega}+\langle{\boldsymbol{S}}{\boldsymbol{u}},{\boldsymbol{u}}\rangle_{\partial\Omega_{\text{S}}}$
(5.37)

where

$a({\boldsymbol{u}},{\boldsymbol{v}}):=({\boldsymbol{\sigma}}({\boldsymbol{u}}),{\boldsymbol{\varepsilon}}({\boldsymbol{v}}))_{\Omega}=\int_{\Omega}{\boldsymbol{\sigma}}({\boldsymbol{u}}):{\boldsymbol{\varepsilon}}({\boldsymbol{v}})\,d\Omega$

which is the usual foundation for a discrete method. However, to obtain a
robust method for the case of $\alpha\rightarrow 0$ or $\beta\rightarrow 0$,
we can introduce a new variable
${\boldsymbol{\lambda}}\in[L_{2}(\partial\Omega_{\text{S}})]^{n}$ and seek
stationary points to

$\mathcal{L}({\boldsymbol{u}},{\boldsymbol{\lambda}}):=\frac{1}{2}a({\boldsymbol{u}},{\boldsymbol{u}})-({\boldsymbol{f}},{\boldsymbol{u}})_{\Omega}-\frac{1}{2}\langle{\boldsymbol{K}}{\boldsymbol{\lambda}},{\boldsymbol{\lambda}}\rangle_{\partial\Omega_{\text{S}}}-\langle{\boldsymbol{\lambda}},{\boldsymbol{u}}\rangle_{\partial\Omega_{\text{S}}}$
(5.38)

where ${\boldsymbol{K}}:={\boldsymbol{S}}^{-1}$ is a flexibility matrix which
simply tends to the zero matrix if $\alpha,\beta\rightarrow 0$, and the Robin
condition becomes a Dirichlet condition. The stationary point to (5.38)
fulfils the variational equations of finding
$({\boldsymbol{u}},{\boldsymbol{\lambda}})\in[H^{1}(\Omega)]^{n}\times[L_{2}(\partial\Omega_{\text{S}})]^{n}$
such that

$\displaystyle
a({\boldsymbol{u}},{\boldsymbol{v}})-\langle{\boldsymbol{\lambda}},{\boldsymbol{v}}\rangle_{\partial\Omega_{\text{S}}}={}$
$\displaystyle({\boldsymbol{f}},{\boldsymbol{v}})_{\Omega}\quad\forall{\boldsymbol{v}}\in[H^{1}(\Omega)]^{n}$
(5.39)
$\displaystyle\langle{\boldsymbol{K}}{\boldsymbol{\lambda}}+{\boldsymbol{u}},{\boldsymbol{\mu}}\rangle_{\partial\Omega_{\text{S}}}={}$
$\displaystyle{0}\quad\forall{\boldsymbol{\mu}}\in[L_{2}(\partial\Omega_{\text{S}})]^{n}$
(5.40)

and we note that, formally,

${\boldsymbol{\lambda}}={\boldsymbol{\sigma}}({\boldsymbol{u}})\cdot{\boldsymbol{n}}$
(5.41)

In the discrete case, we can now formulate an ALM by adding a penalty term and
replacing ${\boldsymbol{\lambda}}$ using (5.41), looking for the minimiser of

$\displaystyle\mathcal{L}_{A}^{h}({\boldsymbol{u}}):={}$
$\displaystyle\frac{1}{2}a({\boldsymbol{u}},{\boldsymbol{u}})-({\boldsymbol{f}},{\boldsymbol{u}})_{\Omega}-\langle{\boldsymbol{\sigma}}({\boldsymbol{u}})\cdot{\boldsymbol{n}},{\boldsymbol{u}}\rangle_{\partial\Omega_{\text{S}}}$
$\displaystyle-\frac{1}{2}\langle{\boldsymbol{K}}{\boldsymbol{\sigma}}({\boldsymbol{u}})\cdot{\boldsymbol{n}},{\boldsymbol{\sigma}}({\boldsymbol{u}})\cdot{\boldsymbol{n}}\rangle_{\partial\Omega_{\text{S}}}+\frac{1}{2}\langle{\boldsymbol{S}}_{h}({\boldsymbol{K}}{\boldsymbol{\sigma}}({\boldsymbol{u}})\cdot{\boldsymbol{n}}+{\boldsymbol{u}}),{\boldsymbol{K}}{\boldsymbol{\sigma}}({\boldsymbol{u}})\cdot{\boldsymbol{n}}+{\boldsymbol{u}}\rangle_{\partial\Omega_{\text{S}}}$
(5.42)

where ${\boldsymbol{S}}_{h}$ is a discrete stiffness matrix, to be chosen. The
minimiser to (5.42) satisfies the variational equation of finding
${\boldsymbol{u}}_{h}\in[V_{h}]^{n}$ such that

$a_{{\boldsymbol{S}}_{h}}({\boldsymbol{u}}_{h},{\boldsymbol{v}})=({\boldsymbol{f}},{\boldsymbol{v}})_{\Omega}\quad\forall{\boldsymbol{v}}\in
V:=[V_{h}]^{n}$ (5.43)

where

$\displaystyle
a_{{\boldsymbol{S}}_{h}}({\boldsymbol{u}},{\boldsymbol{v}}):={}$
$\displaystyle
a({\boldsymbol{u}},{\boldsymbol{v}})-\langle{\boldsymbol{u}}+{\boldsymbol{K}}{\boldsymbol{\sigma}}({\boldsymbol{u}})\cdot{\boldsymbol{n}},{\boldsymbol{\sigma}}({\boldsymbol{v}})\cdot{\boldsymbol{n}}\rangle_{\partial\Omega_{\text{S}}}-\langle{\boldsymbol{\sigma}}({\boldsymbol{u}})\cdot{\boldsymbol{n}},{\boldsymbol{v}}+{\boldsymbol{K}}{\boldsymbol{\sigma}}({\boldsymbol{v}})\cdot{\boldsymbol{n}}\rangle_{\partial\Omega_{\text{S}}}$
(5.44)
$\displaystyle+\langle{\boldsymbol{K}}{\boldsymbol{\sigma}}({\boldsymbol{u}})\cdot{\boldsymbol{n}},{\boldsymbol{\sigma}}({\boldsymbol{v}})\cdot{\boldsymbol{n}}\rangle_{\partial\Omega_{\text{S}}}+\langle{\boldsymbol{S}}_{h}({\boldsymbol{u}}+{\boldsymbol{K}}{\boldsymbol{\sigma}}({\boldsymbol{u}})\cdot{\boldsymbol{n}}),{\boldsymbol{v}}+{\boldsymbol{K}}{\boldsymbol{\sigma}}({\boldsymbol{v}})\cdot{\boldsymbol{n}}\rangle_{\partial\Omega_{\text{S}}}$
(5.45)

which is related to the Nitsche method for interfaces in [64, 65], and a
variant of the method of Juntunen and Stenberg [78] for Poisson’s problem with
Robin boundary conditions. With the particular choice

${\boldsymbol{S}}_{h}=\left((h/\gamma_{0}){\boldsymbol{I}}+{\boldsymbol{K}}\right)^{-1}$
(5.46)

we regain the standard Nitsche method for the Dirichet problem if
${\boldsymbol{K}}$ is the zero matrix, and if ${\boldsymbol{K}}$ is nonzero we
approach the minimiser of (5.37) as $h\rightarrow 0$. Thus the method is
robust also in the limit of zero flexibility.

####  5.4.2 One–sided conditions in contact

We now wish to activate the Robin boundary only if
${\boldsymbol{u}}\cdot{\boldsymbol{n}}-g>0$, corresponding to contact with a
springy foundation at a distance $g$ from the elastic body. Since this
condition is only on the normal part of the displacement, we consider the case
of slip, i.e., we choose

${\boldsymbol{K}}=\alpha{\boldsymbol{n}}\otimes{\boldsymbol{n}}$

Setting
$\sigma_{n}:={\boldsymbol{n}}\cdot{\boldsymbol{\sigma}}\cdot{\boldsymbol{n}}$
and $u_{n}={\boldsymbol{u}}\cdot{\boldsymbol{n}}$, the linear case is then to
find stationary points to (5.38) simplified as

$\mathcal{L}({\boldsymbol{u}},\lambda_{n}):=\frac{1}{2}a({\boldsymbol{u}},{\boldsymbol{u}})-({\boldsymbol{f}},{\boldsymbol{u}})_{\Omega}-\frac{1}{2}\langle\alpha\lambda_{n},\lambda_{n}\rangle_{\partial\Omega_{\text{S}}}-\langle\lambda_{n},u_{n}-g\rangle_{\partial\Omega_{\text{S}}}$
(5.47)

where formally $\lambda_{n}=\sigma_{n}({\boldsymbol{u}})$. In the case of
contact we now have the KKT condition

$\displaystyle u_{n}-g+\alpha\lambda_{n}\leq{}$ $\displaystyle 0$ (5.48)
$\displaystyle\lambda_{n}\leq{}$ $\displaystyle 0$ (5.49)
$\displaystyle\lambda_{n}\left(u_{n}-g+\alpha\lambda_{n}\right)={}$
$\displaystyle 0$ (5.50)

which we can formally rewrite as

$\lambda_{n}=-\gamma[(u_{n}-g+\alpha\lambda_{n})-\gamma^{-1}\lambda_{n}]_{+}$
(5.51)

Proceeding as in (5.28), the equilibrium equation resulting from (5.47) is

$({\boldsymbol{f}},{\boldsymbol{v}})_{\Omega}=a({\boldsymbol{u}},{\boldsymbol{v}})-\langle\lambda_{n},v_{n}\rangle_{\partial\Omega_{\text{S}}}$
(5.52)

and seeing as

$\displaystyle-\langle\lambda_{n},v_{n}\rangle_{\partial\Omega_{\text{S}}}={}$
$\displaystyle-\langle\lambda_{n},v_{n}+\alpha\mu_{n}\rangle_{\partial\Omega_{\text{S}}}+\langle\alpha\lambda_{n},\mu_{n}\rangle_{\partial\Omega_{\text{S}}}$
$\displaystyle={}$
$\displaystyle-\langle\lambda_{n},v_{n}+(\alpha-\gamma^{-1})\mu_{n}\rangle_{\partial\Omega_{\text{S}}}+\langle(\alpha-\gamma^{-1})\lambda_{n},\mu_{n}\rangle_{\partial\Omega_{\text{S}}}$
(5.53)

with $\mu_{n}$ arbitrary, we find that the discrete augmented Lagrangian can
be written

$\displaystyle\mathcal{L}_{A}^{h}({\boldsymbol{u}},\lambda_{n}):={}$
$\displaystyle\frac{1}{2}a({\boldsymbol{u}},{\boldsymbol{u}})-({\boldsymbol{f}},{\boldsymbol{u}})_{\Omega}+\frac{1}{2}\|\gamma^{1/2}[(u_{n}-g+(\alpha-\gamma^{-1})\lambda_{n})]_{+}\|^{2}_{\partial\Omega_{\text{S}}}$
$\displaystyle+\frac{1}{2}\langle(\alpha-\gamma^{-1})\lambda_{n},\lambda_{n}\rangle_{\partial\Omega_{\text{S}}}$
(5.54)

and with $\lambda_{n}\approx\sigma_{n}({\boldsymbol{u}})$,

$\displaystyle\mathcal{L}_{A}^{h}({\boldsymbol{u}}):={}$
$\displaystyle\frac{1}{2}a({\boldsymbol{u}},{\boldsymbol{u}})-({\boldsymbol{f}},{\boldsymbol{u}})_{\Omega}+\frac{1}{2}\|\gamma^{1/2}[u_{n}-g+(\alpha-\gamma^{-1})\sigma_{n}({\boldsymbol{u}})]_{+}\|^{2}_{\partial\Omega_{\text{S}}}$
$\displaystyle+\frac{1}{2}\langle(\alpha-\gamma^{-1})\sigma_{n}({\boldsymbol{u}}),\sigma_{n}({\boldsymbol{u}})\rangle_{\partial\Omega_{\text{S}}}$
(5.55)

the minimiser of which is ${\boldsymbol{u}}\in V$ satisfying

$\displaystyle({\boldsymbol{f}},{\boldsymbol{v}})_{\Omega}={}$ $\displaystyle
a({\boldsymbol{u}},{\boldsymbol{v}})+\langle\gamma[u_{n}-g+(\alpha-\gamma^{-1})\sigma_{n}({\boldsymbol{u}})]_{+},v_{n}+(\alpha-\gamma^{-1})\sigma_{n}({\boldsymbol{v}})\rangle$
$\displaystyle+\langle(\alpha-\gamma^{-1})\sigma_{n}({\boldsymbol{u}}),\sigma_{n}({\boldsymbol{v}})\rangle_{\partial\Omega_{\text{S}}}\quad\forall{\boldsymbol{v}}\in
V$ (5.56)

which coincides with (5.43) in contact, and gives an additional penalty on the
condition $\sigma_{n}({\boldsymbol{u}})=0$ if there is no contact. Choosing
now

$\gamma=(h/\gamma_{0}+\alpha)^{-1}\;\Rightarrow\;\alpha-\gamma^{-1}=-\frac{h}{\gamma_{0}}$
(5.57)

we obtain the same penalty on the normal stress as in [32], which does not
destroy the positive definite nature of the problem if we take
$\gamma_{0}>\gamma_{C}$ where $\gamma_{C}$ is the (stiffness dependent)
constant in the inverse inequality

$\|h^{1/2}{\boldsymbol{\sigma}}({\boldsymbol{v}})\cdot{\boldsymbol{n}}\|_{\partial\Omega_{\text{S}}}^{2}\leq\gamma_{C}a({\boldsymbol{v}},{\boldsymbol{v}})\,\quad\forall{\boldsymbol{v}}\in
V$ (5.58)

###  5.5 Stabilising the Kirchhoff plate model

####  5.5.1 Approximation with independent rotations and displacement

In the Kirchhoff plate model, posed on a domain $\Omega\subset\mathbb{R}^{2}$
with boundary $\partial\Omega$, we seek an out–of–plane (scalar) displacement
$u$ to which we associate the strain (curvature) tensor

${\boldsymbol{\varepsilon}}(\nabla u):=\frac{1}{2}\left(\nabla\otimes(\nabla
u)+(\nabla u)\otimes\nabla\right)=\nabla\otimes\nabla u=\nabla^{2}u$ (5.59)

and the plate stress (moment) tensor

$\displaystyle{\boldsymbol{\sigma}}_{P}(\nabla u)$
$\displaystyle:=\text{D}\left({\boldsymbol{\varepsilon}}(\nabla
u)+\nu(1-{\nu})^{-1}\nabla\cdot\nabla u\,{\boldsymbol{I}}\right)$ (5.60)
$\displaystyle=\text{D}\left(\nabla^{2}u+\nu(1-\nu)^{-1}\Delta
u{\boldsymbol{I}}\right)$ (5.61)

where

$\text{D}=\frac{Et^{3}}{12(1+\nu)}$ (5.62)

where $t$ denotes the plate thickness.

The Kirchhoff clamped problem then takes the form: given the out–of–plane
(scaled) load $t^{3}f$, find the displacement $u$ such that

$\displaystyle\nabla\cdot\left(\nabla\cdot{\boldsymbol{\sigma}}_{P}(\nabla
u)\right)=t^{3}f$ in $\Omega$ (5.63) $\displaystyle u=0$ on $\partial\Omega$
(5.64) $\displaystyle{\boldsymbol{n}}\cdot\nabla u=0$ on $\partial\Omega$
(5.65)

The corresponding variational problem takes the form: Find the displacement
$u\in H^{2}_{0}(\Omega)$ such that

$a_{P}(\nabla u,\nabla v)=(f,v)_{\Omega}\qquad\forall v\in H^{2}_{0}(\Omega$
(5.66)

where

$\displaystyle a_{P}(\nabla v,\nabla w)$
$\displaystyle:=(t^{-3}{\boldsymbol{\sigma}}_{P}(\nabla
v),{\boldsymbol{\varepsilon}}(\nabla w))_{\Omega}$ (5.67)

From a computational point of view (5.66) is cumbersome since it requires
$C^{1}$–conforming elements or carefully constructed nonconforming
approximations. It is therefore common to use instead the Mindlin–Reissner
model which is described by the following partial differential equations:

$\begin{split}-t^{-3}\nabla\cdot{\boldsymbol{\sigma}}_{P}({\boldsymbol{\theta}})-\kappa\,t^{-2}\left(\nabla
u-{\boldsymbol{\theta}}\right)&=0,\quad\text{in
}\Omega\subset\mathbb{R}^{2},\\\ -\kappa\,t^{-2}\,\nabla\cdot\left(\nabla
u-{\boldsymbol{\theta}}\right)&=f,\quad\text{in }\Omega,\end{split}$ (5.68)

where ${\boldsymbol{\theta}}$ is the rotation of the median surface and
$\kappa$ is a shear correction factor. We note that this relaxes the
continuity requirement on $u$ and that, as $t\rightarrow 0$, tends to the
Kirchhoff model. However, the requirement on the approximation to allow
$|\nabla u-{\boldsymbol{\theta}}|\rightarrow 0$ is difficult to realise in the
discrete setting and if this condition cannot be met, shear locking occurs,
destroying the approximation properties of the discrete model. The ALM can
offer an alternative approach in which we enforce the requirement $\nabla
u={\boldsymbol{\theta}}$ by a Lagrange multiplier. To this end we consider the
Lagrangian

$\mathcal{L}(u,{\boldsymbol{\theta}},{\boldsymbol{\lambda}}):=\frac{1}{2}a_{P}({\boldsymbol{\theta}},{\boldsymbol{\theta}})+({\boldsymbol{\lambda}},\nabla
u-{\boldsymbol{\theta}})_{\Omega}-(f,u)_{\Omega}$ (5.69)

The Euler stationary points of (5.69) satisfy the weak system

$\displaystyle
a_{P}({\boldsymbol{\theta}},{\boldsymbol{\vartheta}})+({\boldsymbol{\lambda}},\nabla
v-{\boldsymbol{\vartheta}})_{\Omega}={}$
$\displaystyle(f,v)_{\Omega}\quad\forall(v,{\boldsymbol{\vartheta}})\in
H_{0}^{1}(\Omega)\times[H_{0}^{1}(\Omega)]^{2},$ (5.70) $\displaystyle(\nabla
u-{\boldsymbol{\theta}},{\boldsymbol{\mu}})_{\Omega}={}$ $\displaystyle
0\quad\forall{\boldsymbol{\mu}}\in[L_{2}(\Omega)]^{2},$ (5.71)

corresponding to the strong form

$\displaystyle-\nabla\cdot{\boldsymbol{\sigma}}_{P}({\boldsymbol{\theta}})$
$\displaystyle=t^{3}{\boldsymbol{\lambda}}\quad\text{in }\Omega$ (5.72)
$\displaystyle-\nabla\cdot{\boldsymbol{\lambda}}$
$\displaystyle=f\quad\text{in }\Omega$ (5.73) $\displaystyle\nabla
u-{\boldsymbol{\theta}}$ $\displaystyle=0\quad\text{in }\Omega$ (5.74)

We now wish to stabilise (5.69) using the ALM. To this end, we use (5.72) to
eliminate ${\boldsymbol{\lambda}}$ and add a penalty term on the side
condition to obtain the augmented discrete functional

$\mathcal{L}_{A}^{h}(u_{h},{\boldsymbol{\theta}}_{h}):=\frac{1}{2}a_{P}({\boldsymbol{\theta}}_{h},{\boldsymbol{\theta}}_{h})-(t^{-3}\nabla\cdot{\boldsymbol{\sigma}}_{P}({\boldsymbol{\theta}}_{h}),\nabla
u_{h}-{\boldsymbol{\theta}}_{h})_{h}+\frac{\gamma}{2}\|\nabla
u_{h}-{\boldsymbol{\theta}}_{h}\|^{2}_{\Omega}-(f,u_{h})_{\Omega}$ (5.75)

where $u_{h}\in V^{h}_{1}$ and ${\boldsymbol{\theta}}_{h}\in[V^{h}_{2}]^{2}$
for some discrete spaces $V^{h}_{1}$ and $V^{h}_{2}$. Here we use the notation

$({\boldsymbol{u}},{\boldsymbol{v}})_{h}:=\sum_{T\in\mathcal{T}_{h}}\int_{T}{\boldsymbol{u}}\cdot{\boldsymbol{v}}\,dxdy$
(5.76)

The Euler equations corresponding to the augmented system are

$A_{h}(({\boldsymbol{\theta}}_{h},u_{h}),({\boldsymbol{\vartheta}},v))=(f,v)\
$ (5.77)

for all $(v,{\boldsymbol{\vartheta}})\in V^{h}_{1}\times[V^{h}_{2}]^{2}$,
where

$A_{h}(({\boldsymbol{\theta}},u),({\boldsymbol{\vartheta}},v)):=a_{P}({\boldsymbol{\theta}},{\boldsymbol{\vartheta}})-(t^{-3}\nabla\cdot{\boldsymbol{\sigma}}_{P}({\boldsymbol{\theta}}),\nabla
v-{\boldsymbol{\vartheta}})_{h}-(\nabla
u-{\boldsymbol{\theta}},t^{-3}\nabla\cdot{\boldsymbol{\sigma}}_{P}({\boldsymbol{\vartheta}}))_{h}+\gamma(\nabla
u-{\boldsymbol{\theta}},\nabla v-{\boldsymbol{\vartheta}})$ (5.78)

Now, if $V_{2}^{h}$ is the space of piecewise linears, the terms
$(\cdot,\cdot)_{h}$ vanish and, seeing as $\theta\in H^{1}(\Omega)$ and thus
$\lambda\in H^{-1}(\Omega)$, we choose $r=1$ in (4.15) and
$\gamma=\gamma_{0}/h^{2}$ to obtain a scheme proposed by Pitkäranta [88]; for
higher order polynomial approximations we recover a GLS stabilisation method
due to Stenberg [96, 6].

####  5.5.2 The plate obstacle problem

We next consider applying the model from the previous Section to a regularised
plate obstacle problem. The continuous model is

$\displaystyle-\nabla\cdot{\boldsymbol{\sigma}}_{P}({\boldsymbol{\theta}})$
$\displaystyle=t^{3}{\boldsymbol{\lambda}}\quad\text{in }\Omega$ (5.79)
$\displaystyle-\nabla\cdot{\boldsymbol{\lambda}}+p$
$\displaystyle=f\quad\text{in }\Omega$ (5.80) $\displaystyle\nabla
u-{\boldsymbol{\theta}}$ $\displaystyle=0\quad\text{in }\Omega$ (5.81)
$\displaystyle p\geq 0,\;u-g+\beta p\geq 0,\;p(u-g+\beta p)$
$\displaystyle=0\quad\text{in }\Omega$ (5.82)

Here, $\beta$ is a given compliance which regularises the problem, in the
limit case of $\beta=0$ (rigid obstacle) we instead have the KKT conditions
$p\geq 0$, $u-g\geq 0$, and $p(u-g)=0$. Note that the regularity in the limit
case is insufficient for the analysis above. Indeed it is well known that
$u\not\in H^{4}(\Omega)$, which is insufficient for the multiplier to be in
$L_{2}$. It is however known that for $\beta>0$, $u\in H^{4}(\Omega)$ if the
interior angles of the domain are smaller than $126^{\circ}$ (see [13]).
Therefore the analysis is valid for all $\beta>0$, since we have

$p\in Q=\left\\{\begin{array}[]{>{\displaystyle}l}L_{2}(\Omega)\quad\text{if
$\beta>0$}\\\ H^{-2}(\Omega)\quad\text{if $\beta=0$}\end{array}\right.$ (5.83)

We see that, again, formally
${\boldsymbol{\lambda}}=-t^{-3}\nabla\cdot{\boldsymbol{\sigma}}_{P}({\boldsymbol{\theta}})$
and that
$p=f-t^{-3}\nabla\cdot\left(\nabla\cdot{\boldsymbol{\sigma}}_{P}({\boldsymbol{\theta}})\right)$.
Following the strategy from Sec. 5.4.2 we write

$p=\epsilon[(u_{n}-g+\beta p)-\epsilon^{-1}p]_{+}$ (5.84)

We need to also stabilise the rotations, and to this end we consider the
discrete Lagrangian

$\displaystyle\mathcal{L}_{A}^{h}({\boldsymbol{\vartheta}},v):=$
$\displaystyle{}\frac{1}{2}A_{h}(({\boldsymbol{\vartheta}},v),({\boldsymbol{\vartheta}},v))+\frac{1}{2}\|\epsilon^{1/2}[u-g-(\epsilon^{-1}-\beta)(f-t^{-3}\nabla\cdot\left(\nabla\cdot{\boldsymbol{\sigma}}_{P}({\boldsymbol{\vartheta}}))\right)]_{+}\|_{h}^{2}$
$\displaystyle-\frac{1}{2}\|(\epsilon^{-1}-\beta)^{1/2}(f-t^{-3}\nabla\cdot(\nabla\cdot{\boldsymbol{\sigma}}_{P}({\boldsymbol{\vartheta}})))\|_{h}^{2}-(f,v)_{\Omega}$
(5.85)

where, considering the limit case $p\in H^{-2}(\Omega)$, we choose $r=2$ and
thus

$\epsilon=(h^{4}/\gamma_{1}+\beta)^{-1}$ (5.86)

with $\gamma_{1}$ a sufficiently large constant. A similar approach has been
suggested by Gustafsson et al. [61, 62] in the context of $C^{1}$
approximations of the clamped Kirchhoff plate with GLS stabilisation, without
specific reference to augmented Lagrangian methods.

##  6 Numerical examples

###  6.1 Cavitation

The problem formulation is that of (5.9)–(5.10). Our numerical experience is
that for the chosen discretization $\gamma_{0}$ should not be chosen too
large; in our example we chose $\gamma_{0}=1/100$.

We consider a domain with an elliptically shaped pocket, with mesh shown in
Fig. 1. The boundary conditions are natural boundary conditions
$(-p\boldsymbol{I}+\mu\nabla\boldsymbol{u})\cdot\boldsymbol{n}=\boldsymbol{0}$
at the left- and right-hand sides. The velocity is set to zero along the floor
of the channel and pocket boundary, and the flow is driven by setting
$\boldsymbol{u}=(1,0)$ at the ceiling. The viscosity is $\mu=1$. We compare
the pressure solution with and without cavitation in Figs. 2–3 and note that
there is a pressure resultant in the cavitation case, creating a lifting
resultant force, cf. [84].

###  6.2 Elastic contact with flexible plane

In this example, we consider an elastic sphere of radius 1 under the load
${\boldsymbol{f}}=(0,0,-50)$ in contact with a flexible plane. The contact is
assumed friction–free, in accordance with the form (5.56). The moduli of
elasticity were chosen as $E=200$ and $\nu=0.33$ and the stabilisation
parameter was taken as $\gamma=100E$. In Figs. 4–6 we show the deformation and
contact pressure for increasing flexibilities of the contact plane.

###  6.3 Plate obstacle problem

The considered example, from [62], concerns a clamped square plate
$\Omega=(0,1)\times(0,1)$ in contact with a rigid obstacle ($\beta=0)$ in the
center of the plate, $g=100((x-1/2)^{2}+(y-1/2)^{2})$. Here $E=1$, $\nu=0$,
$t=1$, and we chose $\gamma_{1}=10E$ and $\gamma_{2}=E/10$. We present a
sample computation using continuous, piecewise $P^{2}$ approximations for both
displacement and rotations on triangular meshes, based on the variational
equations resulting from minimization of the Lagrangian (5.85). The mesh is
shown in Fig. 7 (left), and the corresponding soultion is given in Figs. 7
(right, with obstacle indicated) and 8. The computational solution agrees well
with that of [62].

## References

  * [1] P. Alart and A. Curnier.  A mixed formulation for frictional contact problems prone to Newton like solution methods.  Comput. Methods Appl. Mech. Engrg., 92(3):353–375, 1991.
  * [2] P. Antolin, A. Buffa, and M. Fabre.  A priori error for unilateral contact problems with Lagrange multipliers and isogeometric analysis.  IMA J. Numer. Anal., 39(4):1627–1651, 2019.
  * [3] K. J. Arrow, L. Hurwicz, and H. Uzawa.  Studies in linear and non-linear programming.  Stanford Mathematical Studies in the Social Sciences, II. Stanford University Press, Stanford, Calif., 1958.  With contributions by H. B. Chenery, S. M. Johnson, S. Karlin, T. Marschak, R. M. Solow.
  * [4] H. J. C. Barbosa and T. J. R. Hughes.  The finite element method with Lagrange multipliers on the boundary: circumventing the Babuška-Brezzi condition.  Comput. Methods Appl. Mech. Engrg., 85(1):109–128, 1991.
  * [5] H. J. C. Barbosa and T. J. R. Hughes.  Circumventing the Babuška-Brezzi condition in mixed finite element approximations of elliptic variational inequalities.  Comput. Methods Appl. Mech. Engrg., 97(2):193–210, 1992.
  * [6] L. Beirão da Veiga, J. Niiranen, and R. Stenberg.  A family of $C^{0}$ finite elements for Kirchhoff plates. I. Error analysis.  SIAM J. Numer. Anal., 45(5):2047–2071, 2007.
  * [7] Z. Belhachmi and F. B. Belgacem.  Quadratic finite element approximation of the Signorini problem.  Math. Comp., 72(241):83–104, 2003.
  * [8] F. Ben Belgacem.  Numerical simulation of some variational inequalities arisen from unilateral contact problems by the finite element methods.  SIAM J. Numer. Anal., 37(4):1198–1216, 2000.
  * [9] F. Ben Belgacem and S. C. Brenner.  Some nonstandard finite element estimates with applications to 3D Poisson and Signorini problems.  Electron. Trans. Numer. Anal., 12:134–148, 2001.
  * [10] F. Ben Belgacem and Y. Renard.  Hybrid finite element methods for the Signorini problem.  Math. Comp., 72(243):1117–1145, 2003.
  * [11] F. Ben Belgacem, Y. Renard, and L. Slimane.  A mixed formulation for the Signorini problem in nearly incompressible elasticity.  Appl. Numer. Math., 54(1):1–22, 2005.
  * [12] S. Bertoluzza.  Wavelet stabilization of the Lagrange multiplier method.  Numer. Math., 86(1):1–28, 2000.
  * [13] H. Blum and R. Rannacher.  On the boundary value problem of the biharmonic operator on domains with angular corners.  Math. Methods Appl. Sci., 2(4):556–581, 1980.
  * [14] D. Boffi and C. Lovadina.  Analysis of new augmented Lagrangian formulations for mixed finite element schemes.  Numer. Math., 75(4):405–419, 1997.
  * [15] S. C. Brenner, L.-Y. Sung, H. Zhang, and Y. Zhang.  A quadratic $C^{0}$ interior penalty method for the displacement obstacle problem of clamped Kirchhoff plates.  SIAM J. Numer. Anal., 50(6):3329–3350, 2012.
  * [16] S. C. Brenner, L.-y. Sung, and Y. Zhang.  Finite element methods for the displacement obstacle problem of clamped plates.  Math. Comp., 81(279):1247–1262, 2012.
  * [17] F. Brezzi, W. W. Hager, and P.-A. Raviart.  Error estimates for the finite element solution of variational inequalities.  Numer. Math., 28(4):431–443, 1977.
  * [18] F. Brezzi, W. W. Hager, and P.-A. Raviart.  Error estimates for the finite element solution of variational inequalities. II. Mixed methods.  Numer. Math., 31(1):1–16, 1978/79.
  * [19] E. Burman and A. Ern.  A nonlinear consistent penalty method weakly enforcing positivity in the finite element approximation of the transport equation.  Comput. Methods Appl. Mech. Engrg., 320:122–132, 2017.
  * [20] E. Burman, M. A. Fernández, and S. Frei.  A Nitsche-based formulation for fluid-structure interactions with contact.  ESAIM Math. Model. Numer. Anal., 54(2):531–564, 2020.
  * [21] E. Burman, M. A. Fernández, S. Frei, and F. M. Gerosa.  A mechanically consistent model for fluid-structure interactions with contact including seepage.  Comput. Methods Appl. Mech. Engrg., 392:Paper No. 114637, 28, 2022\.
  * [22] E. Burman and P. Hansbo.  Fictitious domain finite element methods using cut elements: I. A stabilized Lagrange multiplier method.  Comput. Methods Appl. Mech. Engrg., 199(41-44):2680–2686, 2010\.
  * [23] E. Burman and P. Hansbo.  Deriving robust unfitted finite element methods from augmented Lagrangian formulations.  In Geometrically unfitted finite element methods and applications, volume 121 of Lect. Notes Comput. Sci. Eng., pages 1–24. Springer, Cham, 2017.
  * [24] E. Burman, P. Hansbo, and M. G. Larson.  Augmented Lagrangian and Galerkin least-squares methods for membrane contact.  Internat. J. Numer. Methods Engrg., 114(11):1179–1191, 2018.
  * [25] E. Burman, P. Hansbo, and M. G. Larson.  Augmented Lagrangian finite element methods for contact problems.  ESAIM Math. Model. Numer. Anal., 53(1):173–195, 2019.
  * [26] E. Burman, P. Hansbo, and M. G. Larson.  Augmented Lagrangian method for thin plates with Signorini boundaries.  In Numerical mathematics and advanced applications—ENUMATH 2019, volume 139 of Lect. Notes Comput. Sci. Eng., pages 509–519. Springer, Cham, 2021.
  * [27] E. Burman, P. Hansbo, M. G. Larson, and R. Stenberg.  Galerkin least squares finite element method for the obstacle problem.  Comput. Methods Appl. Mech. Engrg., 313:362–374, 2017.
  * [28] R. Bustinza and F.-J. Sayas.  Error estimates for an LDG method applied to Signorini type problems.  J. Sci. Comput., 52(2):322–339, 2012.
  * [29] K. L. Cascavita, F. Chouly, and A. Ern.  Hybrid high-order discretizations combined with Nitsche’s method for Dirichlet and Signorini boundary conditions.  IMA J. Numer. Anal., 40(4):2189–2226, 2020.
  * [30] Z. Chen.  On the augmented Lagrangian approach to Signorini elastic contact problem.  Numer. Math., 88(4):641–659, 2001.
  * [31] F. Chouly.  An adaptation of Nitsche’s method to the Tresca friction problem.  J. Math. Anal. Appl., 411(1):329–339, 2014.
  * [32] F. Chouly and P. Hild.  A Nitsche-based method for unilateral contact problems: numerical analysis.  SIAM J. Numer. Anal., 51(2):1295–1307, 2013.
  * [33] F. Chouly, P. Hild, V. Lleras, and Y. Renard.  Nitsche-based finite element method for contact with Coulomb friction.  In Numerical mathematics and advanced applications—ENUMATH 2017, volume 126 of Lect. Notes Comput. Sci. Eng., pages 839–847. Springer, Cham, 2019.
  * [34] F. Chouly, P. Hild, and Y. Renard.  A Nitsche finite element method for dynamic contact: 1. Space semi-discretization and time-marching schemes.  ESAIM Math. Model. Numer. Anal., 49(2):481–502, 2015.
  * [35] F. Chouly, P. Hild, and Y. Renard.  A Nitsche finite element method for dynamic contact: 2. Stability of the schemes and numerical experiments.  ESAIM Math. Model. Numer. Anal., 49(2):503–528, 2015.
  * [36] 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.
  * [37] F. Chouly, F. Mathieu, P. Hild, R. Mlika, J. Pousin, and Y. Renard.  An overview of recent results on Nitsche’s method for contact problems.  In Geometrically Unfitted FEM, theory and applications. Proceedings from the 2016 UCL workshop. Springer Verlag, 2018.
  * [38] F. Chouly, R. Mlika, and Y. Renard.  An unbiased Nitsche’s approximation of the frictional contact between two elastic structures.  Numer. Math., 139(3):593–631, 2018.
  * [39] M. Cicuttin, A. Ern, and T. Gudi.  Hybrid high-order methods for the elliptic obstacle problem.  J. Sci. Comput., 83(1):Paper No. 8, 18, 2020.
  * [40] P. Coorevits, P. Hild, K. Lhalouani, and T. Sassi.  Mixed finite element methods for unilateral problems: convergence analysis and numerical studies.  Math. Comp., 71(237):1–25, 2002.
  * [41] L. De Lorenzis, J. A. Evans, T. J. R. Hughes, and A. Reali.  Isogeometric collocation: Neumann boundary conditions and contact.  Comput. Methods Appl. Mech. Engrg., 284:21–54, 2015.
  * [42] G. Drouet and P. Hild.  Optimal convergence for discrete variational inequalities modelling Signorini contact in 2D and 3D without additional assumptions on the unknown contact set.  SIAM J. Numer. Anal., 53(3):1488–1507, 2015.
  * [43] G. Duvaut and J.-L. Lions.  Inequalities in mechanics and physics, volume 219 of Grundlehren der Mathematischen Wissenschaften.  Springer-Verlag, Berlin-New York, 1976.  Translated from the French by C. W. John.
  * [44] C. Eck, J. Jarušek, and M. Krbec.  Unilateral contact problems, volume 270 of Pure and Applied Mathematics (Boca Raton).  Chapman & Hall/CRC, Boca Raton, FL, 2005.  Variational methods and existence theorems.
  * [45] A. Ern and J.-L. Guermond.  Finite elements II—Galerkin approximation, elliptic and mixed PDEs, volume 73 of Texts in Applied Mathematics.  Springer, Cham, [2021] ©2021.
  * [46] M. Fabre.  A priori error for unilateral contact problems with augmented lagrange multipliers and isogeometric analysis.  Technical report, 2018.
  * [47] M. Fabre, J. Pousin, and Y. Renard.  A fictitious domain method for frictionless contact problems in elasticity using Nitsche’s method.  SMAI J. Comput. Math., 2:19–50, 2016.
  * [48] R. S. Falk.  Error estimates for the approximation of a class of variational inequalities.  Math. Comput., 28:963–971, 1974.
  * [49] R. S. Falk.  Approximation of an elliptic boundary value problem with unilateral constraints.  Rev. Française Automat. Informat. Recherche Opérationnelle, 9(R-2):5–12, 1975.
  * [50] P. E. Farrell, L. Mitchell, and F. Wechsung.  An augmented Lagrangian preconditioner for the 3D stationary incompressible Navier-Stokes equations at high Reynolds number.  SIAM J. Sci. Comput., 41(5):A3073–A3096, 2019.
  * [51] M. Fortin.  An analysis of the convergence of mixed finite element methods.  RAIRO Anal. Numér., 11(4):341–354, iii, 1977.
  * [52] 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.
  * [53] T. Führer, N. Heuer, and E. P. Stephan.  On the DPG method for Signorini problems.  IMA J. Numer. Anal., 38(4):1893–1926, 2018.
  * [54] R. Glowinski.  Numerical methods for nonlinear variational problems.  Springer Series in Computational Physics. Springer-Verlag, New York, 1984\.
  * [55] R. Glowinski and P. Le Tallec.  Augmented Lagrangian and operator-splitting methods in nonlinear mechanics, volume 9 of SIAM Studies in Applied Mathematics.  Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 1989.
  * [56] R. Glowinski, J.-L. Lions, and R. Trémolières.  Numerical analysis of variational inequalities, volume 8 of Studies in Mathematics and its Applications.  North-Holland Publishing Co., Amsterdam-New York, 1981.  Translated from the French.
  * [57] R. Glowinski and A. Marrocco.  Sur l’approximation, par éléments finis d’ordre un, et la résolution, par pénalisation-dualité, d’une classe de problèmes de Dirichlet non linéaires.  Rev. Française Automat. Informat. Recherche Opérationnelle Sér. Rouge Anal. Numér., 9(R-2):41–76, 1975\.
  * [58] T. Gudi and K. Porwal.  A $C^{0}$ interior penalty method for a fourth-order variational inequality of the second kind.  Numer. Methods Partial Differential Equations, 32(1):36–59, 2016\.
  * [59] T. Gustafsson, R. Stenberg, and J. Videman.  Mixed and stabilized finite element methods for the obstacle problem.  SIAM J. Numer. Anal., 55(6):2718–2744, 2017.
  * [60] T. Gustafsson, R. Stenberg, and J. Videman.  On finite element formulations for the obstacle problem—mixed and stabilised methods.  Comput. Methods Appl. Math., 17(3):413–429, 2017.
  * [61] T. Gustafsson, R. Stenberg, and J. Videman.  Nitsche’s method for the obstacle problem of clamped Kirchhoff plates.  In Numerical mathematics and advanced applications—ENUMATH 2017, volume 126 of Lect. Notes Comput. Sci. Eng., pages 407–415. Springer, Cham, 2019.
  * [62] T. Gustafsson, R. Stenberg, and J. Videman.  A stabilised finite element method for the plate obstacle problem.  BIT, 59(1):97–124, 2019.
  * [63] T. Gustafsson, R. Stenberg, and J. Videman.  On Nitsche’s method for elastic contact problems.  SIAM J. Sci. Comput., 42(2):B425–B446, 2020.
  * [64] A. Hansbo and P. Hansbo.  A finite element method for the simulation of strong and weak discontinuities in solid mechanics.  Comput. Methods Appl. Mech. Engrg., 193(33-35):3523–3540, 2004\.
  * [65] P. Hansbo.  Nitsche’s method for interface problems in computational mechanics.  GAMM-Mitt., 28(2):183–206, 2005.
  * [66] P. Hansbo, A. Rashid, and K. Salomonsson.  Least-squares stabilized augmented Lagrangian multiplier method for elastic contact.  Finite Elem. Anal. Des., 116:32–37, 2016.
  * [67] P. Hansbo, A. Rashid, and K. Salomonsson.  Least-squares stabilized augmented Lagrangian multiplier method for elastic contact.  Finite Elem. Anal. Des., 116:32–37, 2016.
  * [68] J. Haslinger.  On numerical solution of a variational inequality of the 4th order by finite element method.  Apl. Mat., 23(5):334–345, 1978.  With a loose Russian summary.
  * [69] J. Haslinger.  Finite element analysis of the Signorini problem.  Comment. Math. Univ. Carolin., 20(1):1–17, 1979.
  * [70] J. Haslinger and I. Hlaváček.  Approximation of the Signorini problem with friction by a mixed finite element method.  J. Math. Anal. Appl., 86(1):99–122, 1982.
  * [71] P. Heintz and P. Hansbo.  Stabilized Lagrange multiplier methods for bilateral elastic contact with friction.  Comput. Methods Appl. Mech. Engrg., 195(33-36):4323–4333, 2006\.
  * [72] M. R. Hestenes.  Multiplier and gradient methods.  J. Optim. Theory Appl., 4:303–320, 1969.
  * [73] P. Hild.  Numerical implementation of two nonconforming finite element methods for unilateral contact.  Comput. Methods Appl. Mech. Engrg., 184(1):99–123, 2000.
  * [74] P. Hild and P. Laborde.  Quadratic finite element methods for unilateral contact problems.  Appl. Numer. Math., 41(3):401–421, 2002.
  * [75] P. Hild and Y. Renard.  A stabilized Lagrange multiplier method for the finite element approximation of contact problems in elastostatics.  Numer. Math., 115(1):101–129, 2010.
  * [76] P. Hild and Y. Renard.  An improved a priori error analysis for finite element approximations of Signorini’s problem.  SIAM J. Numer. Anal., 50(5):2400–2419, 2012.
  * [77] Q. Hu, F. Chouly, P. Hu, G. Cheng, and S. P. A. Bordas.  Skew-symmetric Nitsche’s formulation in isogeometric analysis: Dirichlet and symmetry conditions, patch coupling and frictionless contact.  Comput. Methods Appl. Mech. Engrg., 341:188–220, 2018.
  * [78] M. Juntunen and R. Stenberg.  Nitsche’s method for general boundary conditions.  Math. Comp., 78(267):1353–1374, 2009.
  * [79] T. Kärkkäinen, K. Kunisch, and P. Tarvainen.  Augmented Lagrangian active set methods for obstacle problems.  J. Optim. Theory Appl., 119(3):499–533, 2003.
  * [80] 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.
  * [81] D. Kinderlehrer and G. Stampacchia.  An introduction to variational inequalities and their applications, volume 88 of Pure and Applied Mathematics.  Academic Press, Inc. [Harcourt Brace Jovanovich, Publishers], New York-London, 1980.
  * [82] T. A. Laursen and V. G. Oancea.  Automation and assessment of augmented Lagrangian algorithms for frictional contact problems.  J. Appl. Mech.-Trans. ASME, 61(4):956–963, 1994.
  * [83] J.-L. Lions.  Quelques méthodes de résolution des problèmes aux limites non linéaires.  Dunod; Gauthier-Villars, Paris, 1969.
  * [84] B. Nilsson and P. Hansbo.  A Stokes model with cavitation for the numerical simulation of hydrodynamic lubrication.  Internat. J. Numer. Methods Fluids, 67(12):2015–2025, 2011.
  * [85] J. A. Nitsche.  Über ein Variationsprinzip zur Lösung von Dirichlet-Problemen bei Verwendung von Teilräumen, die keinen Randbedingungen unterworfen sind.  Abh. Math. Univ. Hamburg, 36:9–15, 1971.
  * [86] M. A. Olshanskii and A. Zhiliakov.  Recycling augmented Lagrangian preconditioner in an incompressible fluid solver.  Numer. Linear Algebra Appl., 29(2):Paper No. e2415, 15, 2022.
  * [87] J. Pitkäranta.  Local stability conditions for the Babuška method of Lagrange multipliers.  Math. Comp., 35(152):1113–1129, 1980.
  * [88] J. Pitkäranta.  Analysis of some low-order finite element schemes for Mindlin-Reissner and Kirchhoff plates.  Numer. Math., 53(1-2):237–254, 1988.
  * [89] M. J. D. Powell.  A method for nonlinear constraints in minimization problems.  In Optimization (Sympos., Univ. Keele, Keele, 1968), pages 283–298. Academic Press, London, 1969.
  * [90] R. T. Rockafellar.  A dual approach to solving nonlinear programming problems by unconstrained optimization.  Math. Programming, 5:354–373, 1973.
  * [91] R. T. Rockafellar.  The multiplier method of Hestenes and Powell applied to convex programming.  J. Optim. Theory Appl., 12:555–562, 1973.
  * [92] R. Scholz.  Mixed finite element approximation of a fourth order variational inequality by the penalty method.  Numer. Funct. Anal. Optim., 9(3-4):233–247, 1987.
  * [93] A. Schröder.  Mixed finite element methods of higher-order for model contact problems.  SIAM J. Numer. Anal., 49(6):2323–2339, 2011.
  * [94] J. C. Simo and T. A. Laursen.  An augmented Lagrangian treatment of contact problems involving friction.  Comput. & Structures, 42(1):97–116, 1992.
  * [95] L. Slimane, A. Bendali, and P. Laborde.  Mixed formulations for a class of variational inequalities.  M2AN Math. Model. Numer. Anal., 38(1):177–201, 2004.
  * [96] R. Stenberg.  A new finite element formulation for the plate bending problem.  In Asymptotic methods for elastic structures (Lisbon, 1993), pages 209–221. de Gruyter, Berlin, 1995.
  * [97] R. Stenberg.  On some techniques for approximating boundary conditions in the finite element method.  J. Comput. Appl. Math., 63(1-3):139–148, 1995.
  * [98] I. Temizer, P. Wriggers, and T. J. R. Hughes.  Contact treatment in isogeometric analysis with NURBS.  Comput. Methods Appl. Mech. Engrg., 200(9-12):1100–1112, 2011.
  * [99] F. Wang, W. Han, and X. Cheng.  Discontinuous Galerkin methods for solving the Signorini problem.  IMA J. Numer. Anal., 31(4):1754–1772, 2011.
  * [100] F. Wang, W. Han, and X.-L. Cheng.  Discontinuous Galerkin methods for solving elliptic variational inequalities.  SIAM J. Numer. Anal., 48(2):708–733, 2010.
  * [101] F. Wang and H. Wei.  Virtual element methods for the obstacle problem.  IMA J. Numer. Anal., 40(1):708–728, 2020.
  * [102] M. F. Wheeler, T. Wick, and W. Wollner.  An augmented-Lagrangian method for the phase-field approach for pressurized fractures.  Comput. Methods Appl. Mech. Engrg., 271:69–85, 2014.
  * [103] B. Wohlmuth.  Variationally consistent discretization schemes and numerical algorithms for contact problems.  Acta Numer., 20:569–734, 2011.
  * [104] P. Wriggers and G. Zavarise.  Application of augmented Lagrangian techniques for non-linear constitutive laws in contact interfaces.  Comm. Numer. Methods Engrg., 9(10):815–824, 1993.
  * [105] G. Zavarise, P. Wriggers, and B. A. Schrefler.  On augmented Lagrangian algorithms for thermomechanical contact problems with friction.  Internat. J. Numer. Methods Engrg., 38(17):2924–2949, 1995.
  * [106] Y. Zeng, J. Chen, and F. Wang.  Error estimates of the weakly over-penalized symmetric interior penalty method for two variational inequalities.  Comput. Math. Appl., 69(8):760–770, 2015.
  * [107] Y. Zeng, J. Chen, and F. Wang.  Convergence analysis of a modified weak Galerkin finite element method for Signorini and obstacle problems.  Numer. Methods Partial Differential Equations, 33(5):1459–1474, 2017.

Figure 1: Mesh for cavitation computations.

Figure 2: Pressure isolines without (left) and with (right) cavitation.

Figure 3: Pressure elevation without (left) and with (right) cavitation.

Figure 4: Deformations for $\alpha=0$ and associated contact pressure.

Figure 5: Deformations for $\alpha=10^{-3}$ and associated contact pressure.

Figure 6: Deformations for $\alpha=10^{-2}$ and associated contact pressure.

Figure 7: Computational mesh and elevation of displacements with obstacle
indicated. Figure 8: Deformation isoplot.