Full Discretisations for Nonlinear Evolutionary Inequalities Based on Stiffly Accurate Runge–Kutta and hp-Finite Element Methods

The convergence of full discretisations by implicit Runge–Kutta and nonconforming Galerkin methods applied to nonlinear evolutionary inequalities is studied. The scope of applications includes differential inclusions governed by a nonlinear operator that is monotone and fulfills a certain growth condition. A basic assumption on the considered class of stiffly accurate Runge–Kutta time discretisations is a stability criterion which is in particular satisfied by the Radau IIA and Lobatto IIIC methods. In order to allow nonconforming hp-finite element approximations of unilateral constraints, set convergence of convex subsets in the sense of Glowinski–Mosco–Stummel is utilised. An appropriate formulation of the fully discrete variational inequality is deduced on the basis of a characteristic example of use, a Signorini-type initial-boundary value problem. Under hypotheses close to the existence theory of nonlinear first-order evolutionary equations and inequalities involving a monotone main part, a convergence result for the piecewise constant in time interpolant is established.


Introduction
Scope of Applications Nowadays, transient heat transfer between different materials under nonlinear unilateral/bilateral transmission conditions [13] and timedependent contact between deformable bodies with possibly nonlinear material behaviour [14,35] receive great attention because of relevant applications in engineering. Time-dependent free boundary problems of a similar character occur in elastoplasticity [23], often connected with viscosity, damage, fatigue and other effects [36].

Time and Space Discretisations
Discretisations for related classes of problems that are based on implicit time integration methods and a Galerkin approach in space arouse the interest of many researchers. Numerous contributions confirm that implicit Runge-Kutta and linear multistep methods are favourable for the time discretisation of dynamic processes described by nonlinear evolution equations; as a small selection of literature we refer to [7,15,22,37] and the references given therein. In particular, due to the stiff nature of the problems, implicit methods show superior stability properties compared to explicit methods. The behaviour of piecewise linear finite element approximations has been studied among others in [18,[24][25][26] for the restricted static and quasistatic regime of linear material behaviour in elasticity and plasticity. The convergence of the more general h-version finite element method (FEM) with respect to the mesh size h, where the order p of the finite elements is low and fixed to p = 2 or p = 3, respectively, is established in [17,19] for static unilateral/non-smooth contact problems. Contrary, the hp-version finite element method achieves convergence not only by refining the mesh size h, but also by increasing the polynomial order p of the finite elements. By the pioneering work of Babuška and co-workers, the exponentially fast convergence of the hp-finite element method is well-known for linear elliptic problems. More recently, the superior convergence properties of the p-and hp-compared to the standard h-finite element method have been proven for certain nonlinear elliptic problems and demonstrated in numerical experiments. In [1] the nonlinear Laplace equation is treated, and in [27] a domain obstacle problem with a general nonlinear monotone elliptic operator of quadratic growth is studied; for the investigation of unilateral contact problems within the range of linear elasticity without and with Tresca friction and related elliptic variational inequalities, see [11,12,16,21,28].
Objective In the present work, our concern is to combine a class of implicit Runge-Kutta methods for time discretisation and the hp-finite element method for space discretisation to arrive at full discretisations of nonlinear first-order in time evolutionary variational inequalities. Our analysis covers problems involving a main part that is monotone and satisfies a certain, not necessarily quadratic, growth condition. A basic assumption on the considered class of stiffly accurate Runge-Kutta time discretisations is a stability criterion which implies algebraic stability and is in particular fulfilled by the widely used Radau IIA method, including as special case the implicit Euler method, and by the Lobatto IIIC methods. In order to allow hp-finite element approximations that lead to nonconforming approximations of the given unilateral constraints, we utilise set convergence of convex subsets in the sense of Glowinski-Mosco-Stummel [17,29,38]. On the basis of a relatively simple model application, a Signorini-type initial-boundary value problem, we deduce an appropriate formulation of the fully discrete variational inequality. Under hypotheses close to the existence theory of nonlinear evolutionary inequalities, dispensing with unnatural regularity requirements on the exact solution to the problem, we are able to establish a convergence result for the piecewise constant in time interpolant.
The present work generalises [8], where full discretisations based on the implicit Euler method and low-order finite element approximations applied to nonlinear evolutionary inequalities are studied. Moreover, it extends the analysis of stiffly accurate Runge-Kutta time discretisations for nonlinear evolution equations given in [15] to nonlinear evolutionary inequalities.
To our knowledge, this work is the first contribution, where a class of implicit Runge-Kutta methods is investigated for nonlinear evolutionary inequalities. For this reason, we justify in more detail the employed novel formulation of the fully discrete variational inequality and deduce the needed auxiliary results for stiffly accurate Runge-Kutta and hp-finite element discretisations. However, to keep the work at a reasonable length, in the convergence proof we only indicate the standard arguments used and refer to the literature for further details.
Outline The work has the following structure. In Sects. 2 and 3 we state a characteristic model problem, a Signorini-type initial-boundary value problem, deduce its formulation as variational inequality, and specify the considered class of time and space discretisations. In Sects. 4 and 5, we introduce a general analytical framework of nonlinear evolutionary variational inequalities and its fully discrete counterpart; in particular, we deduce a result on the existence and uniqueness of the fully discrete solution as well as a priori bounds for the discrete approximation values. A convergence result for the piecewise constant in time interpolant related to the relaxed formulation of the evolutionary variational inequality is finally established in Sect. 6.
Notations To conclude the introduction we collect auxiliary notations that are used throughout. For further details on the considered function spaces, we refer to [9].
We employ standard vector and matrix notations such as x = (x 1 , . . . , x d ) ∈ R d and A = (a ij ) 1≤i,j ≤s ∈ R s×s . The Euclidian inner product and the associated norm are given by x · y = (x|y) = x 1 y 1 + · · · + x d y d and x = √ (x|x) for x, y ∈ R d . The partial derivative with respect to x i is denoted by ∂ x i , and, as usual, we employ the abbreviations ∇, div, for the Nabla operator, divergence and Laplacian.
Throughout, we consider a finite time interval I = (0, T ) ⊂ R, where T > 0. We suppose that the spatial domain Ω ⊂ R d with d = 2 or d = 3 is bounded and Lipschitz; to simplify matters, in connection with the finite element approximation we consider a polyhedral domain. The Dirichlet boundary condition and the unilateral constraint are prescribed on the sets Γ D and Γ S , respectively, where we suppose ∂Ω = Γ = Γ D ∪ Γ S and Γ D ∩ Γ S = ∅. We require the Dirichlet boundary to be of positive measure, i.e. meas Γ D > 0, and, for simplicity of the exposition, we do not prescribe Neumann boundary conditions on a boundary part, although this would pose no additional difficulty in the analysis.
We suppose that the exponent, which arises in the Signorini-type initial-boundary value problem or the hypotheses on the functional defining the variational inequality, respectively, satisfies r ∈ [2, ∞); the associated exponent is denoted by r * = r r−1 and evidently satisfies r * ∈ (1,2]. The space of continuous functions C (Ω) and the Lebesgue and Sobolev spaces such as L 2 (Ω), W 1,r (Ω) and W −1,r * (Ω) are equipped with the standard norms.
We recall that a Gelfand triple X ⊂ H ⊂ X * is formed by a real and separable Hilbert space (H, (·|·) H , · H ) and a real, separable and reflexive Banach space (X, · X ) such that X is dense and continuously embedded in H . As standard, the dual space (X * , · X * ) is equipped with the norm x * X * = sup{|x * (x)| : x ∈ X, x X = 1}, and for x ∈ X and x * ∈ X * the duality pairing is given by x * |x X * ×X = x * (x), which extends the inner product (·|·) H .
For exponents q ∈ [1, ∞) and a Banach space Z the space L q (I, Z) of Bochner integrable abstract functions z : I → Z is endowed with the norm and the space C (I , Z) of Z-valued continuous functions on I is equipped with the topology of uniform convergence. For elements z ∈ L q (I, Z) and z * ∈ (L q (I, Z)) * = L q * (I, Z * ) with q * = q q−1 the duality pairing is defined by As the distinction is clear from the context, we use the same letters for instance for elements in a Banach space Z and functions in L q (I, Z), respectively. For a Hilbert space (Z, (·|·) Z , · Z ) the product space (Z s , (·|·) Z s , · Z s ) is endowed with inner product and corresponding norm given by for z = (z 1 , . . . , z s ) ∈ Z s and z = ( z 1 , . . . , z s ) ∈ Z s . Besides, for product spaces the duality pairing is defined componentwise; for instance, we set for z * = (z * 1 , . . . , z * s ) T ∈ (L q * (I, Z * )) s and z = (z 1 , . . . , z s ) T ∈ (L q (I, Z)) s . Finally, we denote by C > 0 a generic constant, possibly with different values at different occasions.

A Signorini-Type Initial-Boundary Value Problem and Its Formulation as Evolutionary Variational Inequality
In this section, we introduce a simple though characteristic model problem, which shall serve as an illustration of the general framework of nonlinear evolutionary inequalities. We first state the nonlinear parabolic initial-boundary value problem involving the r-Laplacian and a free boundary condition of Signorini type and then indicate its reformulation as evolutionary variational inequality. Moreover, we specify the underlying function spaces and the basic properties of the governing operators. We recall that r ∈ [2, ∞) and r * = r r−1 .
A Signorini-Type Initial-Boundary Value Problem We consider the following initial-boundary value problem for a real-valued function u : Ω × I → R: involving the given right-hand side f : Ω × I → R, the (time-independent) function g : ∂Ω → R defining the boundary condition, as well as the given initial condition u 0 .

Formulation as Evolutionary Variational Inequality
In regard to the fully discrete analogue derived in Sect. 3, we indicate the reformulation of the initial-boundary value problem (2.1) as evolutionary variational inequality. We meanwhile assume the solution u : Ω × I → R to the Signorini-type initial-boundary value problem to be sufficiently regular; the regularity requirements on u will be specified below. Moreover, we denote by v : Ω × I → R a sufficiently regular function that fulfills the Dirichlet boundary condition v = g on Γ D × I and the Signorini boundary condition v ≥ g on Γ S × I . Pointwise multiplication of the differential equation in (2.1) with v − u, integration over the spatial domain, and an application of Green's identity yields where n denotes the outer normal to the boundary and ∂ n u = ∇u · n; to keep the formulae in a compact format, we occasionally do not indicate the dependence of the integrand on the spatial variable. Due to the fact that the contribution from the boundary term is non-negative, the evolutionary variational inequality results.
Underlying Function Spaces The model problem suggests to consider the function spaces as underlying Gelfand triple. The solution space is given by and thus functions f ∈ L r * (I, X * ) defining the right-hand side of the differential equation are admitted. For a function g ∈ X = W 1,r (Ω) the imposed boundary conditions are incorporated in the (nonvoid) closed convex set K = {w ∈ X : w = g a.e. on Γ D and w ≥ g a.e. on Γ S }.
Due to the continuous embedding X ⊂ C (I , H ) the initial condition u(0) = u 0 is well-defined. Furthermore, in order to ensure consistency with the boundary conditions, initial values u 0 ∈ K are considered.

Governing Operators and Functionals
In regard to abstract evolutionary inequalities treated in Sect. 4, the first term in the variational inequality (2.2) corresponds to the duality pairing of ∂ t u(·, t) − f (·, t) ∈ X * and v(·, t) − u(·, t) ∈ X. The decisive second term is captured by the nonlinear functional ϕ : For the subsequent considerations it is essential that the functional ϕ is monotoneconvex [32] and satisfies a growth condition, see follows.
Connection to Nonlinear Evolution Equations For the above model problem or, more generally, for problems governed by a monotone operator, the relation ϕ(u, v) = Au|v − u X * ×X establishes a connection between the nonlinear operator A : X → X * corresponding to the weak formulation of the involved differential operator and the monotone-convex functional ϕ : X × X → R, see also Sect. 4 for further details.

Full Discretisation of the Signorini-Type Initial-Boundary Value Problem and Associated Discrete Variational Inequality
In the following, we introduce the considered full discretisations for nonlinear evolutionary inequalities. In Sect. 3.1 we specify the basic assumptions on the implicit Runge-Kutta time discretisations, and in Sect. 3.2 we treat the nonconforming hpfinite element space discretisations utilising the concept of set convergence of convex subsets in the sense of Glowinski-Mosco-Stummel. In order to deduce an appropriate fully discrete analogue to the variational inequality (2.2), we follow the arguments indicated in Sect. 2. In Sects. 4 and 5 the approach will be extended to more general nonlinear evolutionary inequalities.

Time Discretisation by Stiffly Accurate Runge-Kutta Methods
Implicit Runge-Kutta methods such as Radau IIA methods are widely used in the context of stiff ordinary differential equations and parabolic evolution equations, primarily due to their favourable stability properties. However, to our knowledge, the present work is the first contribution where implicit Runge-Kutta time discretisations are studied for nonlinear evolutionary inequalities. In the following, we introduce the basic assumptions on the considered class of implicit Runge-Kutta methods and a fundamental auxiliary result. In order to state the defining relation for the timediscrete solution and its reformulation as time-discrete variational inequality, we find it useful to focus on the model problem and to utilise an approach via abstract evolutionary differential inclusions, as this resembles the approach for evolution equations exploited in [15]. In this section, for better readability and as the positive integer M ∈ N corresponding to the time discretisations are meanwhile fixed, we do not indicate the dependence of the quantities I, t m , τ and u m on M; the initial condition and its approximation are distinguished by the notations u(0) and u 0 , respectively.

Stiffly Accurate Runge-Kutta Methods
Stiffly Accurate Runge-Kutta Methods We study implicit Runge-Kutta methods involving s stages, defined by a real matrix A = (a ij ) 1≤i,j ≤s ∈ R s×s , a vector of positive weights b = (b i ) 1≤i≤s ∈ R s and a vector of associated nodes c = (c i ) 1≤i≤s ∈ R s with 0 < c i ≤ 1 for 1 ≤ i ≤ s. Throughout, we use the abbreviations α = diag(a s1 , . . . , a ss ) ∈ R s×s and γ = diag(c 1 , . . . , c s ) ∈ R s×s ; further, we set 1 = (1, . . . , 1) T ∈ R s and e s = (0, . . . , 0, 1) T ∈ R s . We employ the following assumptions. accurate and consistent, that is, we have b i = a si > 0 for 1 ≤ i ≤ s and further a i1 + · · · + a is = c i for 1 ≤ i ≤ s, in particular a s1 + · · · + a ss = c s = 1.
(ii) The coefficient matrix A is invertible and the matrix Example Methods The above hypotheses are in particular satisfied by the first-order implicit Euler method with s = a 11 = c 1 = 1 and by the third-order two-stage Radau IIA method with coefficients in both cases the matrix C is equal to zero, and thus Hypothesis 3.1 implying algebraic stability [22] is obviously satisfied. More generally, as shown in [15], the s-stage Radau IIA method of order 2s − 1 and the s-stage Lobatto IIIC method of order 2s − 2 fulfill Hypothesis 3.1.

Basic Auxiliary Result
The following auxiliary result ensures that the first term in the employed discrete analogue of the variational inequality (2.2) defines a monotoneconvex functional, see also [15,Lemma 3.4]. It generalises the relation known for the implicit Euler method. For the two-and three-stage Radau IIA methods analogous identities are established, whereas in the general case a lower bound is derived.
Lemma 3.1 Hypothesis 3.1 implies the following statements.
(ii) (a) For the two-stage Radau IIA method the following identity is valid for arbi- (b) For the three-stage Radau IIA method the following identity holds for any Proof For the convenience of the reader, we indicate the proof of the first result found in [15].
(i) A straightforward calculation implies the following identity: A first Gauss elimination step applied to the symmetric matrix M yields applying the transposed transformation matrix from the right further implies By Hypothesis 3.1 the matrix C is positive semi-definite, which implies positive semidefiniteness of M and thus yields the statement.
(ii) (a) For the special case of the two-stage Radau IIA method Gauss elimination applied to M leads to the matrix decomposition Employing the abbreviation the claimed result is obtained as follows: (b) Similar considerations yield the claimed relation for the three-stage Radau IIA method.
Remark We conjecture that analogous relations to Lemma 3.1(ii) are valid for Radau IIA methods of higher stage number. However, as in the present work our focus is on the derivation of a convergence result for stiffly implicit Runge-Kutta methods applied to evolutionary inequalities, but not on their rate of convergence, we do not further exploit this point. In view of Theorem 5.1 we note that under the above hypotheses the matrix A = αA −1 is positive definite; in particular, for the two-stage Radau IIA method this follows from the stated relation (with x 0 = 0) and elementary considerations.

Time Discretisation of Differential Inclusions
Differential Inclusions In order to account for the introduction of stiffly accurate Runge-Kutta time discretisations for evolutionary inequalities (2.2), we employ an abstract formulation of the Signorini-type initial-boundary value problem (2.1) as differential inclusion Here, the nonlinear operator A is related to the r-Laplacian and the function f corresponds to the time-dependent right-hand side; furthermore, the set K as given in (2.3) and the normal cone map [4] capture the prescribed Dirichlet boundary condition and the unilateral constraint of the Signorini boundary condition.

Time-Discrete Solution For an equidistant partition of the time interval
and a given initial approximation u 0 ≈ u(0), the above differential inclusion motivates to determine approximation values u m ≈ u(t m ) for 1 ≤ m ≤ M through a recurrence relation of the form In compact matrix and vector notation, the relation for the stages can be rewritten as which will be justified in Sect. 3.1.3; note that αN K (U m−1 ) = N K (U m−1 ), due to the positivity of the weights.

Time-Discrete Variational Inequality
Time-Discrete Solution to Model Problem The application of a stiffly accurate Runge-Kutta method to the Signorini-type initial-boundary value problem (2.1) yields where the stages are subject to the Dirichlet boundary conditions U m−1,i = g on Γ D and the unilateral constraints

Time-Discrete Variational Inequality
In order to deduce an appropriate timediscrete variational inequality, we imitate the procedure for the time-continuous prob-

i and applying integration by parts implies
we point out that it is essential to employ the reformulation (3.1b) instead of (3.1a), since it is then evident that the boundary term is non-negative, due to boundary conditions and the positivity of the weights α i > 0 for 1 ≤ i ≤ s. Finally, summation of the relations for 1 ≤ i ≤ s yields an appropriate time-discrete analogue of (2.2) for in compact notation, we thus obtain with the functional ϕ defined in (2.4).
Governing Functionals It is essential that the term which results from the employed reformulation (3.1b) and summation, defines a monotone-convex functional; this follows from the fact that the associated real func- Provided that the functional ϕ is monotone-convex, it is evident that the second integral involving the positive weights α i for 1 ≤ i ≤ s defines a monotoneconvex functional.

Space Discretisation by the hp-Finite Element Method
In the following, we introduce approximations based on the hp-finite element method and the concept of Glowinski-Mosco-Stummel set convergence of convex subsets. For our considerations in the subsequent sections, it is convenient to employ a general setting of nonconforming Galerkin methods. Finally, we state the fully discrete analogue of the evolutionary variational inequality, which gives rise to a finitedimensional and hence computable problem.

Approximations Based on the hp-Finite Element Method
Basic Assumptions For simplicity and as this is no restriction of generality, we consider a polygonal planar domain Ω ⊂ R 2 ; in fact, the p-and hp-finite element approximation on curvilinear domains is well-understood, see [5], and the analysis to follow can be extended to higher dimensional domains by tensor product approximation. Moreover, we assume that there is only a finite number of end points Γ D ∩ Γ S , which guarantees the density relation K ∩ C ∞ (Ω) = K, see [25]. For ease of notation we do not indicate the dependence of the quantities h, i, p, q, G, E on the positive integer N ∈ N.
Finite Element Subspaces We let (T (N ) ) N ∈N be a sequence of shape regular meshes, see [34] that covers the bounded polygonal domain Ω ⊂ R 2 and consists of affine quadrilaterals Q ∈ T (N ) with diameter h Q such that all corners of the boundary Γ and all end points Γ D ∩ Γ S are nodes of the mesh T (N ) . Obviously, for every edge E of T (N ) there exists a unique element Q E ∈ T (N ) such that E is an edge of Q E . For each affine quadrilateral Q ∈ T (N ) we denote by p Q ∈ N the associated polynomial degree, assuming that neighbouring elements have comparable polynomial degrees, that is, there exists a constant c > 0 such that the relation we denote by Π p (Q) the tensor product space of polynomials of degree p on Q and define the finite element subspaces (X (N ) ) N ∈N through Gauss-Lobatto Quadrature Similar to [11,12,28] we employ Gauss-Lobatto quadrature in the discretisation procedure. To this end, we introduce the Gauss-Lobatto quadrature nodes (ξ j ) 0≤j ≤q , given as zeros of the function where L q denotes the Legendre polynomial of degree q ≥ 1. Obviously, the nodes ξ 0 = −1 and ξ q = 1 coincide with the end points of the reference interval [−1, 1]. As well-known, there exist positive weights such that the associated quadrature formula is exact for all polynomials φ up to degree 2q − 1, Interpolation Operators For any edge E of a mesh T (N ) we introduce the quadrature order q E = p Q E , and by affine transformation F E : [−1, 1] → E we define the set G E of q E + 1 associated Gauss-Lobatto nodes. Local and global interpolation operators associated with the Gauss-Lobatto nodes are thus given by Convex Subsets In addition, we introduce the sets of edges on the Dirichlet and Signorini boundary, respectively, which provides the associated sets of Gauss-Lobatto nodes Choosing the Gauss-Lobatto nodes as control points of the boundary conditions, we define a family of closed convex subsets (K (N ) ) N ∈N by We note that the approximation is nonconforming, since the subset K (N ) is generally not contained in the set K, in particular, for polynomial degree ≥ 2 or a non-convex obstacle g. Instead we have and we are able to prove set convergence, see Lemma 3.3.

Interpolation Operators Likewise, local and global interpolation operators associated with pairs of Gauss-Lobatto nodes
Polynomial Interpolation Error For later use we recall the following result on the polynomial interpolation error in the reference interval E = (−1, 1) and the reference square Q = (−1, 1) 2 , respectively.

Set Convergence
Set Convergence We employ the notion of Glowinski-Mosco-Stummel set convergence of convex subsets, introduced in [29,38], further analysed in [3] and refined in [17,25]. Provided that the following conditions are satisfied, the sequence (K (N ) ) N ∈N is said to G-converge to the set K for N → ∞, that is, (i) For any subsequence (N ) ∈N such that N → ∞ for → ∞ and for any sequence (x N ) ∈N such that x N ∈ K (N ) and (x N ) ∈N converges to x for → ∞, weakly in X, it follows x ∈ K. (ii) There exist a dense subset K ⊂ K and mappings (N ) : K → K (N ) for N ∈ N such that for any element x ∈ K the sequence ( (N ) (x)) N ∈N converges to x for N → ∞, strongly in X, and further we have (N ) As shown in [3,20], e.g., G-convergence of sets gives rise to convergence (sometimes called epiconvergence) of convex lower semi-continuous functionals, when considering their epigraphs, which are convex closed sets. The following result ensures set convergence of the sequence (K (N ) ) N ∈N defined in (3.3a), (3.3b) towards K, see also (2.3) for the definition of K.
Proof Classical h-FEM convergence for a similar variational problem is already treated in [19], where Newton-Cotes formulae are used instead of Gauss-Lobatto quadrature. Inspecting the proof of [19,Theorem 4.1] shows that the norm convergence for a fixed quadrature order hinges on the positiveness of the quadrature weights, what is satisfied for all quadrature orders with Gauss-Lobatto quadrature. Therefore in the following we may focus on the case where h Q is fixed for all Q ∈ T (N ) and min{p Q : Q ∈ T (N ) } → ∞. Moreover, as every equality constraint w = g can be expressed through the two inequality constraints w ≤ g and w ≥ g, it suffices to treat the Signorini boundary condition. In order to verify the first requirement of G-convergence we have to show that for any λ ∈ C (Γ ) with λ| Γ S ≥ 0 it follows that using duality with respect to (L 1 , L ∞ ). Moreover, since the mesh T (N ) is supposed to be independent of N , we can simply consider the above integral on any fixed edge E. Thus, we fix λ ∈ C (E) with λ ≥ 0 and also the polynomial degree q = q E . Similarly as [28] we approximate the function λ by a combination of Bernstein polynomials B q , and with the local mapping F E : and then on E, Using on the other hand that λ q−1 (w (N ) − g (N ) )| E is a polynomial of degree 2q − 1 and that hence the above integral can be evaluated exactly by the Gauss-Lobatto quadrature we obtain for Due to the fact that the quadrature weights ω j > 0 are positive and further λ q−1 ≥ 0 as well as ((w (N ) − g (N ) ) • F E )(ξ j ) ≤ 0 as w (N ) ∈ K (N ) , we arrive at In view of (3.6) this proves our claim (3.4). It remains to prove the second requirement. By the finiteness assumption, due to [25] it follows that K ∩ C ∞ (Ω) is dense in K. Consequently, we may consider the dense subset K = K ∩ C ∞ (Ω) and may define (N ) as the Lagrange interpolation operator on K in X (N ) . Moreover, as w ∈ K satisfies the constraints in K pointwise, we have (N ) w ∈ K (N ) for all w ∈ K. By Theorem 3.2(ii) it follows (N ) w → w as N → ∞ in W 1,r (Ω). Altogether, this yields the stated result.

Fully Discrete Variational Inequality
Assumptions on Approximating Spaces In accordance with Sects. 3.2.1 and 3.2.2, we henceforth employ the following assumptions on the general Galerkin method that can be realized by the hp-finite element method.

Nonlinear Evolutionary Variational Inequalities
In Sects. 2 and 3 the focus is on a specific application, the Signorini-type initialboundary value problem (2.1), its formulation as evolutionary variational inequality (2.2) and the derivation of the fully discrete analogue (3.7). In this section, we introduce a general framework of nonlinear evolutionary variational inequalities involving monotone-convex functionals. Furthermore, following [8] we introduce an integrated and relaxed reformulation of the problem and state a uniqueness result.

Analytical Framework
Underlying Function Spaces and Governing Functional The underlying Gelfand triple X ⊂ H ⊂ X * , with continuous and dense embeddings, is formed by the Banach space X, its dual space X * and the Hilbert space H , as usual identified with its dual space. Basic assumptions on the governing nonlinear functional are collected in the following hypotheses. (i) For any x ∈ X the relation ϕ(x, x) = 0 holds. (ii) Convexity. For any z ∈ X the function ψ = ϕ(z, ·) : X → R is convex, that is, for all x, y ∈ X and σ ∈ [0, 1] we have (iii) Sequential lower semi-continuity. For any z ∈ X the function ψ = ϕ(z, ·) : X → R is sequentially lower semi-continuous, that is, for every x ∈ X and for any sequence (iv) Hemi-continuity. The functional ϕ : X × X → R is hemi-continuous, that is, for all x, x, y ∈ X the mapping is continuous.
(vi) Coercivity. The functional ϕ : X × X → R is coercive, that is, there exist constants C > 0 and C 0 ≥ 0 such that for all x ∈ X ϕ(x, 0) ≤ C 0 − C x r X .
(vii) Growth condition. There exists a constant C > 0 such that the operator ϕ : X × X → R satisfies the following growth condition for any x, y ∈ X:

Connection to Nonlinear Evolution Equations
The above hypotheses on the functional ϕ are in accordance with the model problem (2.1) governed by the r-Laplacian or, more generally, evolutionary problems governed by a nonlinear operator A : X → X * that is hemi-continuous, monotone, coercive and satisfies a certain growth condition, see for instance [15] and references given therein. In this case, the functional inherits the properties of A. Sequential lower semi-continuity and hemi-continuity are evident; moreover, monotonicity, coercitivity and a certain growth condition hold.
(a) Provided that the operator A is monotone, that is, for all x, y ∈ X the relation Ay − Ax|y − x X * ×X ≥ 0 holds, the associated functional ϕ is monotone, since ϕ(x, y) + ϕ(y, x) = − Ay − Ax|y − x X * ×X ≤ 0 for all x, y ∈ X. (b) If A is coercive with exponent r ∈ [2, ∞), that is, for any x ∈ X a bound of the form Ax|x X * ×X ≥ C x r X − C 0 holds with constants C > 0 and C 0 ≥ 0, the associated functional ϕ is coercive, since −ϕ(x, 0) = Ax|x X * ×X ≥ C x r X − C 0 and thus ϕ(x, 0) ≤ C 0 − C x r X for all x ∈ X. (c) Finally, provided that the operator A satisfies the growth condition Ax X * ≤ C(1 + x r−1 X ) for all x ∈ X, the associated functional ϕ fulfills the bound |ϕ(x, y)| ≤ Ax X * y − x X ≤ C(1 + x r−1 X ) y − x X for any x, y ∈ X. These conditions are also satisfied for functionals given by ϕ(x, y) = Ax|y − x X * ×X + f (y) − f (x), provided that additionally f : X → R is convex and lower semi-continuous.

Related Functionals and Operators
Under the assumptions of Hypothesis 4.1, the growth condition on ϕ and an application of Young's inequality yield the estimate A sufficient condition for measurability of the function I → R : t → ϕ(x(t), y(t)), where x, y : I → X are assumed to be Bochner integrable, is that ϕ is a Caratheodory function, which here amounts to continuity of (x, y) → ϕ(x, y) with respect to the norm topology, see also [2]. As a consequence, the Nemytskii operator associated with the nonlinear functional ϕ and a related functional, are well-defined, since the bound is ensured for all x, y ∈ L r (I, X). Moreover, the functional Φ inherits the property of hemi-continuity, see [8] for further details. For the derivation of the main result, an enhancement of sequential lower semi-continuity is needed.

Hypothesis 4.2 The functional Φ defined in (4.1) is continuous with respect to the second argument and further satisfies the LSC condition, that is, for all sequences (x j ) j ∈N and (y j ) j ∈N
such that x j converges to x for j → ∞, strongly in L r (I, X) and y j converges to y for j → ∞, weakly in L r (I, X), it follows that Solution Space In the following, the Banach space which is continuously embedded in the space of continuous functions C (I , H ), will be the solution space, and a nonvoid closed convex set K ⊂ X shall capture the imposed boundary conditions. Throughout, we employ the abbreviations L r (I, K) = w ∈ L r (I, X) : w(t) ∈ K almost everywhere in I , and utilise the continuous embedding X K ⊂ C (I , K). For the problem data, the function defining the right-hand side and the initial value, we suppose f ∈ L r * (I, X * ) and u 0 ∈ K. For the sake of brevity, we do not indicate the dependence of X and X K on the considered time interval I , the exponent r ∈ [2, ∞) given by Hypothesis 4.1, the associated exponent r * = r r−1 ∈ (1, 2] and the underlying Gelfand triple.

Evolutionary Inequality and Reformulations
Nonlinear Evolutionary Variational Inequality The formulation of the model problem (2.1) as evolutionary variational inequality (2.2) motivates the following general form of a nonlinear evolutionary variational inequality involving a nonlinear convexmonotone functional, see also Sect. 4.1 for the hypotheses on ϕ and the definition of the solution space.

Problem 4.1
For a given initial value u 0 ∈ K and a given function f ∈ L r * (I, X * ) find u ∈ X K such that u(0) = u 0 and for v(t) ∈ K and almost all t ∈ I .

Problem 4.2 For given
for all v ∈ L r (I, K), where the functional Φ is defined by (4.1).
Relaxed Formulation of the Evolutionary Inequality A reformulation of (4.2) which presupposes reduced regularity requirements on the solution is as follows.

Problem 4.3
For given u 0 ∈ K and f ∈ L r * (I, X * ) find u ∈ L r (I, K) ∩ C (I , H ) such that u(0) = u 0 and Any solution to the integrated variational inequality is a solution to the relaxed variational inequality. However, in order to show that a solution to (4.4) with d dt u ∈ L r * (I, X * ) is likewise a solution to (4.3), additional assumptions are needed. Namely, it is required that in a decomposition of the functional Φ, associated with the decomposition of the monotone-convex functional as ϕ(x, y) = ψ(y) − ψ(y) + ϕ(x, y), where ψ : X → R is assumed to be convex and lower semi-continuous, the functional Φ(u, ·) is continuous on L r (I, X); for further details, see [8, Lemmas 2.2 and 2.6].

A Uniqueness Result
A Uniqueness Result The following result ensures uniqueness of the solution to the integrated variational inequality and thus to the evolutionary variational inequality, see also [8,Lemma 2.7]. Under the above stated additional requirements guaranteeing the equivalence of the integrated and relaxed formulations, also uniqueness of the solution to the relaxed variational inequality follows.

Full Discretisations of Nonlinear Evolutionary Inequalities
In the following, we state the fully discrete counterparts of the considered nonlinear evolutionary variational inequality and of the integrated and relaxed formulations. Moreover, we deduce an existence and uniqueness result as well as a priori bounds for the discrete solution, which are an essential ingredient in the convergence analysis. Throughout, we assume that Hypothesis 4.1 on the governing nonlinear functional as well as Hypotheses 3.1 and 3.2 on the time and space discretisations are satisfied.
Notations For notational simplicity, as the integers M, N ∈ N corresponding to the time and space discretisations are meanwhile fixed, we do not indicate the dependence of the discrete solution values U m−1,i and u m−1 = U m−1,s on M, N ; as before, the initial condition and its approximation are distinguished by the notations u(0) and u 0 , respectively. We recall the compact vector notation U m−1 = (U m−1,1 , . . . , U m−1,s ) T and the relation u m = U m−1,s for any 1 ≤ m ≤ M. In view of the integrated formulation of the discrete variational inequality, we henceforth identify (in notation) the stage values and their piecewise constant interpolant U = (U ·,1 , . . . , U ·,s ) T , defined by

Fully Discrete Variational Inequality and Reformulations
Fully Discrete Nonlinear Evolutionary Variational Inequality In regard to (3.7) we consider the following discrete analogue of Problem 4.1.

Problem 5.1
For a given initial approximation u 0 ∈ K (N ) and a given approximation F ∈ (P 0 (I, X * )) s find U ∈ (P 0 (I, K (N ) )) s such that U(0) = u 0 1 and Integrated Discrete Nonlinear Evolutionary Variational Inequality For piecewise constant in time interpolants, the integral over the time domain simplifies to a Riemann sum. In particular, for any 1 ≤ i ≤ s we have see also (4.1). As a consequence, multiplication of (5.1) by the time increment and summation yields the following discrete analogue of Problem 4.2.
Problem 5.2 For u 0 ∈ K (N ) and F ∈ (P 0 (I, X * )) s find U ∈ (P 0 (I, K (N ) )) s such that U(0) = u 0 1 and for all V ∈ (P 0 (I, K (N ) )) s . Proof For the following considerations it is convenient to rewrite (5.1) as

Relaxed Discrete Nonlinear Evolutionary Variational Inequality
(a) Existence: The left-hand side in the above relation involves the positive definite matrix A, the positive weights α i for 1 ≤ i ≤ s and the monotone-convex functional ϕ; the right-hand side is a fixed linear functional applied to the argument V m−1 − U m−1 . Thus, at each time step, the existence of the stage value U m−1 is ensured by [32,39].
(b) Uniqueness: In order to prove uniqueness of the stage values, we assume the existence of two solutions U m−1 and U m−1 such that Inserting V m−1 = U m−1 into the first relation and V m−1 = U m−1 into the second identity, adding both relations and applying Lemma 3.1 yield

A Priori Estimates for the Discrete Solution
A Priori Estimates A priori estimates for the discrete approximation values and increments are basic tools in the proof of the convergence result.
depending in particular on the problem data. Furthermore, the estimate Proof (a) Inserting V m−1 = 0 into the discrete variational inequality (5.1) yields

by a slight reformulation the relation
The coercivity of ϕ with exponent r ∈ [2, ∞) implies see also Hypothesis 4.1. By means of Lemma 3.1 and an application of Young's inequality with additional (small) scaling ε > 0 we obtain where C(r, ε) = ε r r and C(r * , 1 ε ) = 1 r * 1 ε r * , and thus by absorption, for ε > 0 chosen sufficiently small, the relation

Multiplication with the time increment and summation further implies
and thus, by the a priori bound for the stage values and the positive definiteness of the matrix A the stated result follows.
Remark In order to utilise the a priori bounds provided by Theorem 5.2 in the proof of the convergence result, it is essential to ensure that the arising quantities

Convergence Result
In the following, our concern is to establish a convergence result for the piecewise constant in time interpolant, under hypotheses close to the existence theory of nonlinear evolutionary equations and inequalities involving a monotone main part. Auxiliary definitions and results are given in Sects. 6.1 and 6.2. In order to keep the work at a reasonable length, especially in the proofs of the density and feasibility results, we only indicate the standard arguments used and instead refer to the literature for further details.

Interpolation, Difference, and Restriction Operators
In this section, we collect several auxiliary results on the piecewise constant and linear in time interpolation operators as well as the difference and restriction operators associated with a stiffly accurate Runge-Kutta method.
In the following, we let Z denote an arbitrary reflexive Banach space or a closed subset thereof. For notational simplicity, we omit the dependences of certain quantities on the integer M ∈ N as long as M is fixed; further, with a slight abuse of notation, we denote the vector of stage values by Z m−1 for 1 ≤ m ≤ M and components of the associated piecewise constant interpolant by Z i for 1 ≤ i ≤ s. We recall that under Hypothesis 3.1 the nodes of the considered stiffly accurate Runge-Kutta method fulfill 0 < c i < 1 for 1 ≤ i ≤ s − 1 and c s = 1, see also Sect. 3.1.1.

Piecewise Constant Interpolation Operator
The space of Z -valued piecewise constant functions associated with the partition of the time interval is given by The piecewise constant interpolation operator associated with a stiffly accurate Runge-Kutta method I 0 : Analogously, for discrete values z 0 ∈ Z and Z m−1 ∈ Z s for 1 ≤ m ≤ M we consider the associated piecewise constant interpolant Z ∈ (P 0 (I, Z )) s and define D : (P 0 (I, Z )) s → (P 0 (I, Z )) s through for any 1 ≤ m ≤ M; as before, we denote z m = Z m−1,s for 1 ≤ m ≤ M.

A Density and Feasibility Result
The following density and feasibility results are utilised in the derivation of our convergence result for full discretisations of nonlinear evolutionary inequalities by stiffly accurate Runge-Kutta and hp-finite element methods.
A Density Result Following the proof of [8] for the implicit Euler method, an analogous statement for the difference operator associated with a stiffly accurate Runge-Kutta method is obtained. We recall the definition of the solution space X K , see Sect. 4.1.

Lemma 6.3
Suppose that v ∈ X K for some exponent 2 ≤ r < ∞ and that f ∈ L r * (I, X * ) for r * = r r−1 .
(iii) There exists a sequence (F (N ) ) N ∈N such that F (N ) ∈ P (M N ) 0 (I, (X (N ) ) * ) for all N ∈ N and Proof For a given function v ∈ X K ⊂ C (I , X) we consider the piecewise constant and linear in time interpolants I (M) v ∈ P (M) (I, K) for = 0, 1. An application of Lemma 6.2 ensures strong convergence of the interpolants as well as convergence of the finite differences The approach used in [8] is to construct piecewise constant and linear in time approximations to the discrete values (v(t m )) M m=1 . Differentiation of the piecewise linear approximation yields an approximation to the first time derivative d dt v. Furthermore, approximations of the initial value v(0) in H and of the function f in L r * (I, X * ) are constructed. For detailed explanations, we refer to [8,Lemma 3.7].

A Feasibility Result
The following result ensures that the considered time and space approximations admit limit points that belong to the convex set K. The proof of the lemma is found in [8, Proof of Lemma 3.13] and [20], see also Sect. 3.2.2 for the definition of G-convergence and Sect. 4.1 for the definition of L r (I, K).  (K (N ) ) N ∈N G-converges to K for N → ∞. Then, any limit point of U (M,N ) or U (M,N ) , respectively, with respect to weak convergence in (L r (I, X)) s belongs to (L r (I, K)) s .

Main Result
Convergence Result With the above stated auxiliary results at hand, we are able to establish the following convergence result for full discretisations of Problem 4.1 by stiffly accurate Runge-Kutta methods and hp-finite element approximations. In regard to the length of the work, we focus on the relaxed formulation of the problem; in this case, differentiability in time is not needed and thus it suffices to employ the first a priori estimate for the discrete solution. ; we recall that by the consistency condition of Hypothesis 3.1 we have A1 = c and consequently Aγ 1 = αA −1 c = α1 as well as α 1 + · · · + α s = 1.

Conclusions
In the present work, our main objective has been to provide a convergence analysis for full discretisations by stiffly accurate Runge-Kutta methods and hp-finite element approximations, applied to nonlinear evolutionary inequalities stemming for instance from a Signorini-type initial-boundary value problem. Following [8] it is straightforward to extend the analysis to time-dependent monotone-convex functionals involving an explicit time-dependency. Further related questions that remain open comprise the extension of the convergence result to nonuniform time grids and other time integration methods such as BDF-methods. In this work, we studied the relaxed formulation of the problem, and it remains open to extend the arguments to the (integrated) variational inequality. For future work, it is also intended to numerically exploit the attainable convergence rate and to deduce error estimates for problems, where the solution satisfies certain regularity requirements; in this context, we refer to Hölder regularity results in space and time for parabolic domain obstacle problems in [31] and to L 2 (I, L r (Ω)) regularity of ∂ t ∇u for the solution u to an evolutionary inequality involving the r-Laplacian in [30]. However, even for free problems on the boundary like a Signorini-type initial-boundary value problem with regular data it is supposed that due to the imposed constraints a severe order reduction compared to classic problems will occur, in general, and that a higher-order convergence rate is only obtained in the interior of the spatial domain.