Two discretisations of the time-dependent Bingham problem

This paper introduces two methods for the fully discrete time-dependent Bingham problem in a three-dimensional domain and for the flow in a pipe also named after Mosolov. The first time discretisation is a generalised midpoint rule and the second time discretisation is a discontinuous Galerkin scheme. The space discretisation in both cases employs the non-conforming first-order finite elements of Crouzeix and Raviart. The a priori error analyses for both schemes yield certain convergence rates in time and optimal convergence rates in space. It guarantees convergence of the fully-discrete scheme with a discontinuous Galerkin time-discretisation for consistent initial conditions and a source term f∈H1(0,T;L2(Ω))\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f\in H^1(0,T;L^2(\Omega ))$$\end{document}.


Introduction
The mathematical model The parabolic variational inequality of this paper describes the flow of a Bingham fluid, which behaves as a rigid material, if the stress is below a certain threshold. If the flow is sufficiently slow (i.e., if the convective term can be omitted for small Reynolds numbers), the inertial terms can be neglected, and, in the case of a uni-directional flow, i.e., the flow in a pipe of cross-section ⊂ R 2 , the prob-The work of the authors was supported by the Deutsche Forschungsgemeinschaft(DFG) in the Priority Program 1748 "Reliable simulation techniques in solid mechanics: Development of non-standard discretization methods, mechanical and mathematical analysis," under the projects CA 151/22 and SCHE1885/1-1.
Discretisations in space Conforming finite element methods (FEM) for the discretisation in space suffer from the non-differentiable L 1 norm, which results in a suboptimal predicted convergence O(h 1/2 ) for general smooth solutions and in a convergence O(h log(1/h)) for a particular prototype example of a circular domain [16]. The traditional analysis for linear problems via the Strang-Fix lemma for non-conforming schemes relies qualitatively on the smoothness of the exact solution. This has been overcome recently by a medius analysis, which employs techniques from a posteriori analysis to prove a priori results [7,9,12,18]. Quasi-optimal error estimates for the stationary Bingham problem were proved for a non-conforming FEM [11] based on the finite element space of Crouzeix and Raviart [10] whenever a stress-like variable belongs to H 1 ( ). This results in the optimal convergence O(h), whenever the solution and a corresponding stress-like variable are smooth. The only other discretisations in space with proved optimal convergence rates were those of [14] (which is equivalent to the discretisation of [11]) and the mixed discretisations of the very recent paper [15]. The generalisation to the three-dimensional case is still open [15] and the combination of it with a time-discretisation is left for future research.
Time discretisation This paper combines the discretisation of [11] in space with two time-discretisations of (tB): A generalised midpoint rule and a discontinuous Galerkin discretisation.
The mathematical analysis of this paper is twofold: The dG(0) scheme is studied under the relatively weak assumption of consistent initial data and a source term f ∈ H 1 (0, T ; L 2 ( )) that guarantees a weak solution u ∈ H 1 (0, T ; L 2 ( )) ∩ L 2 (0, T ; H 1 0 ( )). This results in a discrete approximation u DG ∈ P 0 (K; CR 1 0 (T )) that is piecewise constant in time with minimal time-step k and maximal time-step k and a Crouzeix-Raviart function in space based on a shape-regular triangulation T with maximal mesh-size h with the plain convergence under the (mild) constraint h k. In fact, according to the best knowledge of the authors (and seemingly that of our anonymous referees) this is the first paper that proves convergence under realistic regularity assumptions on the weak solution. In the absence of further regularity of the weak solution u, any convergence rate in (1.1) remains unclear. This paper identifies certain best-approximation terms of ∇u and a stress-type variable σ ∈ ∂ W (∇u) (derived in Theorem 2.3) plus data oscillation terms of f as an upper bound for the errors in (1.1). The backward Euler scheme is a modified dG(0) scheme with a slightly different treatment of the source term. In this low regularity regime, there are no other convergence results for the generalised midpoint scheme. In all numerical examples below f is piecewise constant and so u DG = u h (for θ = 1 in the generalised mid-point rule with τ = t for the backward Euler scheme). The second look at the class of time-discretisations hence solely concerns the generalised mid-point rules under the (possibly) unrealistic regularity assumptions u ∈ W m,1 (0, T ; L 2 ( )) ∩ C(0, T ; H 1 0 ( )) and σ ∈ C(0, T ; L 2 ( ; R 2 )) that allow for a pointwise evaluation of the variational inequality at the prescribed times τ for m = 2, 3. In the second analysis for smooth solutions, the convergence rates are O(h + k) and for a uniform mid-point rule O(h + k 2 ). This improves the known result for a conforming FEM of [21].
Although discontinuous Galerkin methods for parabolic equations are already included in textbooks [19], to the best of the authors knowledge only [2,8] define a discontinuous Galerkin FEM for a time-dependent variational inequality.
An overview over some time-stepping schemes for problem (tB) can be found in [13,17]; see also the references therein. However, proofs of convergence rates are rare and only concern the combination of continuous linear finite elements in space with backward differences in time, which result in a convergence O(h 1/2 + k) in [21] under the much stronger regularity assumptions u ∈ L ∞ (0, T ; H 2 ( )),u ∈ L 2 (0, T ; H 1 ( )) ∩ L ∞ (0, T ; L 2 ( )), andü ∈ L 2 (0, T ; H −1 ( )). This suboptimal convergence in space should be contrasted with the optimal convergence rate obtained in this paper, which is a consequence of the non-conforming discretisation in space.
Structure of the paper The paper is organised as follows. Section 2 specifies notation and states the main results in Sect. 2.4. The proofs of the error estimates for the generalised midpoint rule follow in Sect. 3, while the error estimate for the discontinuous Galerkin FEM is proved in Sect. 4. Section 5 generalises the results to the Bingham flow in a three-dimensional domain. Section 6 concludes the paper with numerical experiments.
General notation Standard notation on Lebesgue and Sobolev spaces applies throughout the paper and the space of L 2 functions with values in R m reads L 2 ( ; R m ) and • dx denotes an integral over a domain ⊂ R 2 in space, while´ Q • d Q denotes an integral over a domain Q = (t 1 , t 2 ) × in space and time and • L 2 ( ) and • L 2 ( Q) denote the respective L 2 norms. The notation • abbreviates the identity mapping.
The space L 2 (t 0 , t 1 ; X ) denotes the space of L 2 functions v : [t 0 , t 1 ] → X for any linear space X , the space W m, p (t 0 , t 1 ; X ) denotes the space of m-times (weakly) differentiable L p (t 0 , t 1 ; X ) functions, whose derivatives belong to L p (t 0 , t 1 ; X ), while C(t 0 , t 1 ; X ) denotes the space of continuous functions from [t 0 , t 1 ] to X . We abbreviate • L p (L q ) := • L p (0,T ;L q ( )) . The formula A B abbreviates that there exists a positive generic mesh-size independent constant C > 0 such that A ≤ C B, A ≈ B abbreviates A B A.

Notation and main results
This section defines some notation in Sect. 2.1 and recalls two operators for the P 1 nonconforming finite element space in Sect. 2.2, before it defines the discrete problems in Sect. 2.3 and states the main results in Sect. 2.4.

Notation
A shape-regular triangulation T of a polygonal bounded Lipschitz domain ⊂ R n is a set of simplices such that = T and any two distinct simplices are either disjoint or share exactly one common face, edge or vertex. Let E denote the set of sides in T and N the set of vertices. Let P p (K ; R m ) denote the set of polynomials on K of total degree ≤ p and P p (T ; R m ) the set of piecewise polynomials, i.e., and let P p (T ) := P p (T ; R). Define the space of continuous piecewise affine functions with homogeneous boundary conditions by where C 0 ( ) denotes the space of continuous functions with homogeneous boundary conditions on ∂ .
For an edge E ∈ E, mid(E) denotes the midpoint of E. The integral mean over a domain ω ⊆ or ω ⊆ [0, T ] is defined as ffl ω • dx :=´ω • dx/|ω|, where |ω| denotes the volume of ω or the length of ω, if it is an edge or an interval. The L 2 -projection onto T -piecewise constant functions or vectors 0 : Let h T ∈ P 0 (T ) denote the piecewise constant mesh-size with h T | K := diam(K ) for all K ∈ T . The square of the oscillations of a function f ∈ L 2 ( ) is defined by The P 1 non-conforming finite element space [10], sometimes named after Crouzeix and Raviart, reads CR 1 0 (T ) := {v CR ∈ P 1 (T ) | v CR is continuous at midpoints of interior edges and vanishes at midpoints of boundary edges}.

Interpolation and companion operators
This subsection recalls the non-conforming interpolation operator and a companion operator for Crouzeix-Raviart functions in 2D.
Non-conforming interpolation operator The non-conforming interpolation operator and satisfies the integral mean property and the approximation and stability estimates Companion operators The design of three conforming companions to any v h ∈ CR 1 0 (T ) begins with the map J 1 : CR 1 0 (T ) → S 1 0 (T ). Given a vertex z ∈ N , let T (z) := {T ∈ T | z ∈ T } denote the set of triangles that contain z. Then J 1 is defined by with the conforming nodal basis function ϕ z ∈ S 1 0 (T ). This operator is also known as enriching operator in the context of fast solvers [3]. For a given edge E := conv{a, b} ∈ E let b E := 6ϕ a ϕ b denote the edge bubble function. Then the operator J 2 : CR 1 0 (T ) → P 2 (T ) ∩ C 0 ( ) is given by For any triangle T ∈ T with T = conv{a, b, c} define the element bubble function b T := 60ϕ a ϕ b ϕ c . The operator J 3 : CR 1 0 (T ) → P 3 (T ) ∩ C 0 ( ) is given by The operators J m : CR 1 0 (T ) → (P m (T ) ∩ C 0 ( )), m = 1, 2, 3, defined above satisfy the approximation and stability properties and the conservation propertieŝ
The following problem admits a unique solution [17, Chapter III, Theorem 6.1].
Definition 2.2 (Problem (tB)) Given initial data u 0 ∈ L 2 ( ), the time-dependent Mosolov problem in two dimensions seeks u ∈ L 2 (0, T ; H 1 0 ( ))∩H 1 (0, T ; H −1 ( )) such that u(0) = u 0 and, for a.e. t ∈ (0, T ) and all v ∈ H 1 0 ( ), it holds that The following theorem proves the existence of some stress-like variable σ . Furthermore, it proves regularity of the solution u of (tB) under smoothness assumptions on f and u 0 . In the stationary case, the existence of the stress-like variable σ is proved in [17, Chapter II, Theorem 6.3].
The generalised midpoint rule requires the evaluation of the right-hand side F(τ , •) at τ and, hence, it requires f ∈ C(0, T ; L 2 ( )), while the subsequent discontinuous Galerkin discretisation merely involves an integrability f ∈ L 2 (0, T ; L 2 ( )). For a piecewise smooth (in time) function v DG , define the left and right limit and the jump at t by Furthermore, define the time-space cylinders Q := I × , where I := (t −1 , t ). Definition 2.5 (dG(0) FEM) The dG(0) FEM seeks u DG ∈ P 0 (K; CR 1 0 (T )) such that u DG (0) = U 0 and, for all = 1, . . . , N and for all v DG ∈ P 0 (K;

Proposition 2.6 There exist unique solutions to problems (tB h ) and (dG(0)).
Proof The existence of solutions to problems (tB h ) and (dG(0)) follows in every time step from the equivalence to the minimization problems The uniqueness follows from the strong convexity of the energies.

Main results
This section states the error estimates for the generalised midpoint rule and the discontinuous Galerkin scheme. The proofs follow in Sects. 3-4.
The first two theorems state error estimates for the generalised midpoint rule (tB h ). The first theorem proves an error estimate for the case that θ > 1/2 for all = 1, . . . , N . Let e := u(t )−u h (t ) and e(τ ) := u(τ )−u h (τ ) and recall the definition of σ from Theorem 2.3.
The following theorem asserts quadratic convergence in time for the uniform midpoint rule.
The following theorem states an error estimate for the discontinuous Galerkin scheme. Let e := u − u DG ∈ L 2 (0, T ; L 2 ( )) and let e(t − ) := lim s t e(s) denote the right limit of e in t and u(t + −1 ) and u(t − ) denote the left and right limits of u in t −1 and t . Recall the non-conforming interpolation operator I NC from Sect. 2.2. Theorem 2.9 (Error estimate for dG(0)) Assume that u ∈ C(0, T ; L 2 ( )) and let f ∈ L 2 (0, T ; L 2 ( )) and u ∈ L 2 (0, T ; L 2 ( )) be defined by

Remark 2.10
A piecewise Poincaré-Friedrichs inequality [4] proves, that ∇ NC • L 2 (Q ) defines in fact a norm on L 2 ([t −1 , t ]; H 1 0 ( ) + CR 1 0 (T )). The following lemma proves a trace inequality with respect to a time interval. This is employed in Corollary 2.12 below, which states a convergence rate for dG(0).
Proof The fundamental theorem of calculus implies for a.e. s ∈ [a, b] An integration over s ∈ [a, b], a triangle inequality and two Hölder inequalities lead to This eventually proves

Proofs of theorems 2.7-2.8
This section starts with a stability result for the generalised midpoint rule for the two cases θ > 1/2 as well as θ = 1/2 for all = 1, . . . , n. Sect. 3.2 proves the error estimates from Theorems 2.7 and 2.8.

Stability
The stability of the generalised midpoint rule will be employed in the error estimate of Theorem 2.7. The stability of the uniform midpoint rule with θ = θ = 1/2, will not enter the proof of the error estimate of Theorem 2.8 below, and is merely given for completeness.

Error estimates
The following theorem controls the error in space for a fixed time τ . The proof follows the arguments of [11] and considers σ := σ (τ ) with σ from (2.6). Recall the constants C comp and C NC from Sect. 2.2.

Theorem 3.2 (Error estimate in space)
Assume that f ∈ C(0, T ; L 2 ( )), σ ∈ C(0, T ; L 2 ( ; R 2 )) andu ∈ C(0, T ; L 2 ( )). The error e(τ ) Proof Theorem 3.2 of [11] proves the existence of σ h, ∈ P 0 (T ; The definition of the subgradients σ h, ∈ ∂ W (∇ NC u h (τ )) and σ ∈ ∂ W (∇u(τ )) results inˆ Since ∇ NC I NC u(τ ) = 0 ∇u(τ ), Jensen's inequality proveŝ This, F(τ , •) −u(τ ) + div σ = 0 in H −1 ( ), and its discrete analogue (3.5) lead to The integral mean property of I NC implies that the last term vanishes. The conservation property (2.4) and the stability property (2.3) provê The conservation property (2.5) and the approximation property (2.3) imply The integral mean property and the approximation property of I NC yield The combination of the previously displayed formulas with two applications of Young's inequality with weight 2/μ lead to The approximation properties of J 3 and I NC , two applications of the Cauchy inequality, and two applications of a weighted Young's inequality with weight 4/μ lead tô The combination of the two previously displayed inequalities concludes the proof.
The following theorem proves a first error estimate for the discretisation (tB h ).

Theorem 3.3 (A priori error estimate for the full discretisation)
If k = k = T /N , θ = θ = 1/2, and, in addition to the above regularity assumptions u ∈ W 3,1 (0, T ; L 2 ( )), then Proof The error analysis with respect to the time discretisation follows as in [1, Proposition 5.1] for elastoplasticity with ||| • ||| replaced by the L 2 norm • L 2 ( ) and the scalar product a(•, •) replaced by the L 2 scalar product. Those arguments lead to max =0,...,N for the generalised midpoint rule with θ > 1/2 and for the uniform midpoint rule. The combination with the error estimate in space of ė(τ )e(τ ) dx from Theorem 3.2 then proves the assertion.

Proof of Theorem 2.7
The combination of Theorem 3.3 with the estimation of h Tuh (τ ) 2 L 2 ( ) from Theorem 3.1 leads to Theorem 2.7.

Proof of Theorem 2.8 For θ = 1/2, a Young inequality implies
with e := u − u h for the nodal interpolation u = I h u ∈ S 1 (K; ffl I u dt, the second term on the right-hand side in (3.7) is bounded by Let H := h T L ∞ ( ) abbreviate the maximal (spatial) mesh-size. Since u(t ) = u(t ) by the definition of u, it holds e(t ) = e . For the first term on the right-hand side of (3.7), the representatioṅ e(τ ) = ( e(t ) − e(t −1 ))/k = (e − e −1 )/k and a triangle inequality lead to Eventually, the right-hand side of (3.8) can be absorbed on the right-hand side of (3.6).
The combination of this with Theorem 3.3 concludes the proof of Theorem 2.8.

Proof of theorem 2.9
This section proves the error estimate for the discontinuous Galerkin discretisation (dG(0)) from Theorem 2.9.

Combination of variational inequalities
The sum of (dG(0)) with the new notation v DG = u DG (to highlight that this will be chosen below to approximate the exact solution) and (4.3) provides (4.4) Since (d/dt)J 3 u DG = 0 in I , the fourth term on the right-hand side of (4.4) readŝ Since u DG is constant on I , the binomial formulas shoŵ The approximation properties (2.3) of J 3 and a Cauchy inequality implŷ Since by assumption u ∈ C(0, T ; L 2 ( )) is continuous, and e = u − u DG , the fifth term on the right-hand side of (4.4) readŝ The definition of the average Therefore, (4.7) results in (4.8) The approximation and stability properties (2.3) of J 3 provê (4.9) Since ∇ NC e L 2 (Q ) can be absorbed, the combination of (4.4) with (4.5)-(4.6), (4.8)-(4.9) andQ (4.13)

Choosing u DG
Define u ∈ L 2 ( ) by u(x) = ffl I u dt for a.e. x ∈ . Note that u ∈ H 1 0 ( ) and ∇u = ffl I ∇u dt (from Fubini's theorem). Define u DG | I := I NC u. The integral mean property Jensen's inequality yieldsÎ It remains to estimate (4.10) and (4.11).
A Cauchy inequality implies for the first term The definition of u DG = I NC u and the approximation property of I NC implŷ Since I NC u = ffl I I NC u dt and ∇u = e. x ∈ . Hence, Jensen's inequality proves The combination of the previous inequalities leads tô Since u is continuous, Since u DG = I NC u and u| I ∈ H 1 ( ), (2.1) implies (4.14) Altogether, This concludes the proof of Theorem 2.9.

Generalization to 3D
In this section, let ⊆ R 3 be a three-dimensional bounded polyhedral Lipschitz domain. Define the vector-valued spaces and let H −1 ( ; R 3 ) denote the dual space of H 1 0 ( ; R 3 ). Furthermore, the linearised Green strain tensor ε(u) = (∇u + ∇u )/2 denotes the symmetric part of the gradient of a velocity function u ∈ H 1 0 ( ; R 3 ) and : denotes the scalar product between two matrices A, B ∈ R 3×3 , i.e., A : Given f ∈ L 2 (0, T ; L 2 ( ; R 3 )) and initial data u 0 ∈ H (div, ) with div u 0 = 0, the time-dependent Bingham flow problem in 3D seeks u ∈ L 2 (0, A direct discretisation of the variational inequality with P 1 non-conforming finite elements is not possible, since´ ε NC (•) : ε NC (•) dx is not positive definite on the P 1 non-conforming finite element space: Korn's inequality fails for CR 1 0 (T ; R 3 ) in general. However, for homogeneous Dirichlet boundary conditions, a straightforward calculation reveals, for all v, w ∈ H 1 0 ( ; Since u(t) ∈ Z , the following problem is equivalent to (5.1): The P 1 non-conforming finite element space in this three-dimensional context consists of piecewise affine functions, which are continuous in the midpoints of interior faces and vanish in midpoints of boundary faces. Let div NC and ε NC denote the piecewise versions of div and ε and define the vector-valued spaces , Let U 0 ∈ CR 1 0 (T ; R 3 ) be some approximation of u 0 . With the notation for the partition of the time interval K from Sect. 2, the generalised midpoint rule seeks u h ∈ S 1 (K; Z CR ) such that u h (0) = U 0 and, for all = 1, . . . , N and for all v h ∈ Z CR , it holdŝ Recall the notation of the jump [•] and the left limit (•) + from Sect. 2. The discontinuous Galerkin FEM in 3D seeks u DG ∈ P 0 (K; Z CR ) such that, for all = 1, . . . , N and for all v DG ∈ P 0 (K; The divergence-free condition u h (t) ∈ Z CR and u DG | I ∈ Z CR can be implemented using a piecewise constant Lagrange multiplier as it is common for the Stokes equations. Define H (div, ; R 3×3 ) := {τ = (τ 1 ; τ 2 ; τ 3 ) : → R 3×3 | ∀m = 1, 2, 3, τ m ∈ H (div, )}, Let t ∈ [0, T ] be arbitrary such that (5.2) holds and f (t) −u(t) ∈ H −1 ( ) and u(t) ∈ H 1 0 ( ). Define W 3D (A) := (μ/4)|A| 2 + g|sym A|. The Euler-Lagrange equations of [11,Theorem 5.1] prove the existence of σ (t) ∈ H (div, ; R 3×3 ) and ξ(t) ∈ L 2 0 ( ) with With this σ , the error estimates of Theorems 2.7-2.9 equally hold in this threedimensional situation. The proof follows as in Sects. 3-4, see also [11,Theorem 5.2] for the analysis of the error estimate in space. We want to highlight that it is essential for the proof that the functions in Z CR are pointwise divergence-free.

Numerical experiments
This section is devoted to two numerical experiments for the two discretisations (tB h ) and (dG(0)) in Sects. 6.2-6.3. Sect. 6.1 describes how the solution to the discrete nonlinear problems are approximated, while Sect. 6.4 draws some conclusions of the numerical experiments.

Numerical realisation
The discrete solutions of (tB h ) and (dG(0)) are approximated by the following Uzawa algorithm from [13]. It is recalled here for the sake of completeness.
Set U old := U . end while Set u h (τ ) := U for the approximation of (tB h ) resp. u DG | I := U for the approximation of (dG(0)). end for Output: Approximation u h or u DG of the solution to (tB h ) resp. (dG(0)).
The numerical experiments in Subsects. 6.2-6.3 employ ρ = μ/g 2 ≈ 44.4 and tol = 10 −6 for the discrete solutions (resp. tol = 10 −6 for the computation of the reference solutions) and we neglect the influence of the error in the Usawa algorithm in the discussion.

Numerical experiment on the square
This example approximates the solution to (tB) on the unit square = (0, 1) 2 and the time intervall [0, T ] with T = 0.2 and initial value u(0) = 0, right-hand side f = 5, viscosity μ = 1, and yield limit g = 0.15. The initial triangulation T 0 is depicted in Fig. 1a. The approximations u DG,m, j defined by (dG(0)) and u θ,m, j defined by (tB h ) for θ = θ = 0.5 and θ = θ = 1 for all time steps are computed on a sequence (T m ) m=0,...,5 of red-refined triangulations (see Fig. 2 for a red-refined triangle). The time step sizes are chosen as 3 − j T /2 for j = 0, 1, . . . , 4 and the parameter tol in the Uzawa algorithm as 10 −6 . A discrete rigid region of u θ,m, j or u DG,m, j at some time t is defined as the union of all triangles where the norm of its (spatial) gradient is less than 1‰ of the maximal norm of the gradient at that time. These rigid regions are plotted in Fig. 3 (resp. Fig. 4) for u DG,m, j (resp. u θ,m, j ) for three different times and for different mesh-sizes and time-step sizes to illustrate the behaviour of the discrete   Note that the refinement in time is chosen in such a way that the evaluation points τ with respect to a coarse time step are also evaluation points on a finer time step. The errors e θ,m, j (resp. e DG,m, j ) read for u h = u DG,m, j (resp. u h = u θ,m, j ) and are plotted in Fig. 5 against the numbers of degrees of freedom (ndof) in T m , while the errors for the number of time steps N are plotted in Fig. 6. For a time step size of k = 3 −5 T /2, the errors from (6.1) and (6.2) for all three methods show a convergence rate of ndof −1/2 = O(h). For θ = 0.5 this convergence rate is also achieved for a time step size of k = 3 −4 T /2. In Figs. 7 and 8, the same errors are plotted against the number of time steps. The errors (6.1) show a convergence rate of O(k) on the finest triangulation T 5 for all considered methods for the first time steps. In particular, the observed convergence rate is O(k) instead of the predicted O(k 1/2 ) for the discretisation (dG(0)). For larger time steps, the error in space already seems to dominate the error, so that the convergence is diminished. The errors (6.2) converge with rate between O(k) and O(k 1/2 ) for θ = 1 on T 5 for the first time steps. This is a slower convergence compared to the predicted convergence rate of Theorem 2.7 and it seems that preasymptotic effects diminish the convergence rate. For θ = 0.5, the convergence between O(k) and the predicted O(k 2 ) is visible. It seems that preasymptotic effects or the error in space slow down the convergence.

Numerical experiment on the L-shaped domain
This example approximates the solution of (tB) on the L-shaped domain = (−1, 1) 2 \ ([0, 1] × [0, 1]) and the time interval [0, T ] with T = 0.2 and with initial value u(0) = 0, right-hand side f = 5, viscosity μ = 1, and yield limit g = 0.15. As in Sect. 6.2, the approximations u DG,m, j defined by (dG(0)) and u θ,m, j defined by (tB h ) for θ = θ = 0.5 and θ = θ = 1 for all time steps are computed on a sequence (T m ) m=0,...,5 of red-refined triangulations (see Fig. 2 for a red-refined triangle). The initial triangulation T 0 is depicted in Fig. 1b. The time step sizes are chosen as 3 − j T /2 for j = 0, 1, . . . , 4 and the parameter tol in the Uzawa algorithm as 10 −6 . As in Sect. 6.2, the rigid regions are plotted in Fig. 9 (resp. Fig. 10) for u DG,m, j (resp. u θ,m, j ) for three different times and for different mesh-sizes and time-step sizes to illustrate the behaviour of the discrete solutions.
A reference solution u ref is computed on the triangulation T 6 with time step size 3 −5 T /2 and tol = 10 −7 . The errors from (6.1) are plotted in Fig. 11 against the num-   For coarse time steps, the convergence is limited by the error in time, while for small time steps, the errors for all methods show a convergence rate slightly smaller than O(h). In the stationary case, the experiments of [11] show a convergence rate slightly larger than O(h 2/3 ) due to the re-entrant corner and the resulting singularity of the solution on an L-shaped domain. However, for the time-dependent situation considered here, it seems that (at least preasymptotic) the singularity has merely a marginal influence. Figures 13 and 14 display the same errors plotted against the number of time steps. The observed convergence rates are O(k) for the first refinement in time on the fines triangulation for all three methods and the two errors. The predicted convergence rate of O(k 2 ) for the uniform midpoint rule cannot be observed. The error in space seems to dominate the approximation error as in Sect. 6.2 on coarser triangulations.

Conclusions
Section 2.4 states error estimates for the discretisations (tB h ) and (dG(0)). The numerical experiments of this section illustrate those results. Moreover, they indicate the following results beyond the theoretical statements from Sect. 2.4.
• The convergence history plots show no preasymptotic effect with respect to h. • The empirical convergence for the discontinuous Galerkin discretisation is O(k) rather than O(k 1/2 ) as k → 0 predicted by Theorem 2.9. • The factor h 2 /k 1/2 in the upper bound for the error of the discontinuous Galerkin discretisation cannot be seen in the numerical experiments, i.e., the error do not increase, when k becomes small for large h. • The condition h T L ∞ ( ) < k √ μ/(8 √ T C NC ) that was needed in the error estimate for the uniform midpoint rule in Theorem 2.8 is satisfied only for the pairs (m, j) ∈ {(3, 0), (4, 0), (5, 0), (5, 1)} for both experiments in Subsects. 6.2-6.3. Although Theorem 2.8 only proves a parameter conditional convergence, the numerical experiments indicate that the convergence is independent of that parameter.
Funding Open Access funding enabled and organized by Projekt DEAL.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.