Staggered explicit-implicit time-discretization for elastodynamics with dissipative internal variables

In computational mechanics one can distinguish two main classes of time-dependent problems, quasistatic and dynamic. Focusing on the latter one, one can further distinguish two other cases: (i) Low-frequency regimes, which are typically related with vibrations of structures and where the energy is not dominantly transmitted through space. In this case implicit time-discretization is relatively efficient, even though large systems of algebraic equations are to be solved at each time step. (ii) High-frequency regimes which arise typically within wave propagation and for which only explicit time-discretizations are reasonably efficient, in particular in three-dimensions. These explicit methods can essentially be used in hyperbolic problems, as mere elastodynamics (treated here) or the elasto-magnetic Maxwell system or some conservation laws. Yet, many applications need to combine the convervative hyperbolic problems with various dissipative processes of a parabolic character, involving typically some internal variables. For parabolic problems, however, explicit methods are known to be problematic due to severe time-step restrictions. Therefore, we propose and analyse an explicit/implicit scheme, using the fractional-step (also called staggered) technique. In fact, the proposed method becomes completely explicit in the absence of gradients of internal variables, as it can be in plasticicity or viscoelasticity as in Sect. 5.1, or an interfacial variant of damage models (so-called delamination) as in [41].

Our starting point is the linear elastodynamic problem: Find the displacement $u:[0,T]\times\varOmega\to{\mathbb{R}}^{d}$ satisfying

		$\displaystyle\varrho\mathchoice{{\buildrel\hskip 1.00006pt\text{\LARGE..}\over{u}}}{{\buildrel\hskip 1.00006pt\text{\Large..}\over{u}}}{{\buildrel\hskip 1.00006pt\text{\large..}\over{u}}}{{\buildrel\hskip 1.00006pt\text{\large..}\over{u}}}-{\rm div}\,\mathbb{C}e(u)=f$	$\displaystyle\text{on }\ \varOmega\ \text{ for }\ t\in(0,T],$		(1.1a)
		$\displaystyle[\mathbb{C}e(u)]\vec{n}+\mathbb{B}u=g$	$\displaystyle\text{on }\ \varGamma\ \text{ for }\ t\in(0,T],$		(1.1b)
		$\displaystyle u|_{t=0}=u_{0},\ \ \ \mathchoice{{\buildrel\hskip 1.00006pt\text{\LARGE.}\over{u}}}{{\buildrel\hskip 1.00006pt\text{\Large.}\over{u}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{u}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{u}}}|_{t=0}=v_{0}$	$\displaystyle\text{on }\ \varOmega,$		(1.1c)

for $T>0$ a fixed time horizon. Here $\varOmega\subset{\mathbb{R}}^{d}$ is a bounded Lipschitz domain, $d=2$ or $3$, $\varGamma$ is its boundary, and $\vec{n}$ the unit outward normal. The dot-notation stands for the time derivative. In (1.1), $\varrho>0$ denotes the mass density, $\mathbb{C}$ is the elasticity tensor which is symmetric and positive definite, $e(u)$ denotes the small-strain tensor defined as $e(u)=\frac{1}{2}(\nabla u)^{\top}{+}\frac{1}{2}\nabla u$, and $\mathbb{B}$ is a symmetric positive semidefinite 2nd-order tensor determining the elastic support on the boundary. The terms appearing in the second member of (1.1) are, the bulk force $f$ and the surface loading $g$. In (1.1c), $u_{0}$ denotes the initial displacement and $v_{0}$ the initial velocity. For a more compact notation, we write the initial-boundary-value problem (1.1) in the following abstract form

$\displaystyle\mathcal{T}^{\prime}\mathchoice{{\buildrel\hskip 1.00006pt\text{\LARGE..}\over{u}}}{{\buildrel\hskip 1.00006pt\text{\Large..}\over{u}}}{{\buildrel\hskip 1.00006pt\text{\large..}\over{u}}}{{\buildrel\hskip 1.00006pt\text{\large..}\over{u}}}+\mathcal{W}^{\prime}u=\mathcal{F}_{u}^{\prime}(t)\ \ \text{ for }\ t\in(0,T],\ \ u|_{t=0}=u_{0},\ \ \ \mathchoice{{\buildrel\hskip 1.00006pt\text{\LARGE.}\over{u}}}{{\buildrel\hskip 1.00006pt\text{\Large.}\over{u}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{u}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{u}}}|_{t=0}=v_{0}.$

(1.2)

Here $\mathcal{T}$ is the kinetic energy, $\mathcal{W}$ is the stored energy, and $\mathcal{F}$ is the external force, while $(\cdot)^{\prime}$ denotes the Gâteaux derivative. In the context of (1.1), we have

$$\mathcal{T}(\color[rgb]{0,0,0}\mathchoice{{\buildrel\hskip 1.00006pt\text{\LARGE.}\over{u}}}{{\buildrel\hskip 1.00006pt\text{\Large.}\over{u}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{u}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{u}}}\color[rgb]{0,0,0})=\int_{\varOmega}\frac{1}{2}\varrho|\color[rgb]{0,0,0}\mathchoice{{\buildrel\hskip 1.00006pt\text{\LARGE.}\over{u}}}{{\buildrel\hskip 1.00006pt\text{\Large.}\over{u}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{u}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{u}}}\color[rgb]{0,0,0}|^{2}\,\mathrm{d}x,\ \quad\mathcal{W}(u)=\int_{\varOmega}\frac{1}{2}\mathbb{C}e(u){:}e(u)\,\mathrm{d}x+\int_{\varGamma}\frac{1}{2}\mathbb{B}u\cdot u\,\mathrm{d}S$$

and

$$\mathcal{F}(t,u)=\int_{\varOmega}f(t){\cdot}u\,\mathrm{d}x+\int_{\varGamma}g(t){\cdot}u\,\mathrm{d}S.$$

Thus $\mathcal{F}_{u}^{\prime}(t)$ is a linear functional, let us denote it shortly by $F(t)$.

For high-frequency wave propagation problems, implicit time discretizations are computationally expensive, especially in the three dimensional setting. Therefore, we focus our attention on explicit methods. The simplest explicit scheme is the following second-order finite difference scheme

$\displaystyle\mathcal{T}^{\prime}\frac{u_{\tau h}^{k+1}-2u_{\tau h}^{k}+u_{\tau h}^{k-1}}{\tau^{2}}+\mathcal{W}_{h}^{\prime}u_{\tau h}^{k}=F_{h}(k\tau).$

(1.3)

Here $\tau>0$ denotes the time step, and $\mathcal{W}_{h}$ (resp. $F_{h}$) denote some discrete approximations of the respective continuous functionals obtained by a suitable finite-element method (FEM) with mesh size $h>0$. For simplicity, we assume the mass density constant (or at least piecewise constant in space) so that the kinetic energy $\mathcal{T}$ does not need any numerical approximation. In particular, a numerical approximation leading to a diagonal the mass matrix $\mathcal{T}^{\prime}$ in (1.3), typically referred to as mass lumping, is an important ingredient so as to obtain efficient explicit methods. Here we will consider that $u$ is discretized in an element-wise constant way so that $\mathcal{T}^{\prime}$ leads to a diagonal matrix even without any approximation. Multiplying (1.3) by $\frac{u_{\tau h}^{k+1}-u_{\tau h}^{k}}{\tau}$ and using the scheme, it is easy to show that the energy preserved is

$$\frac{1}{2}\langle\color[rgb]{0,0,0}{\mathcal{T}}^{\prime}\color[rgb]{0,0,0}\frac{u_{\tau h}^{k+1}-u_{\tau h}^{k}}{\tau},\frac{u_{\tau h}^{k+1}-u_{\tau h}^{k}}{\tau}\rangle+\frac{1}{2}\langle\mathcal{W}_{h}^{\prime}u_{\tau h}^{k+1},u_{\tau h}^{k}\rangle.$$

This is an energy, i.e., a positive quantity under the following Courant-Fridrichs-Lewy (CFL) condition [20]

$\displaystyle\big{\langle}\mathcal{T}^{\prime}u_{h},u_{h}\big{\rangle}\geq\frac{\tau^{2}}{4}\big{\langle}{\cal W}_{h}^{\prime}u_{h},u_{h}\big{\rangle}$

(1.4)

for any $u_{h}$ from the respective finite-dimensional subspace. The CFL typically bounds the time discretization step $\tau=\mathscr{O}(h_{\rm min})$ with $h_{\rm min}$ the smallest element size on a FEM discretization. This method has frequently been used and analysed from various aspects, including comparison with implicit time discretizations, cf. e.g. [34, 33]. However, the form of the discrete stored energy $\frac{1}{2}\langle\mathcal{W}_{h}^{\prime}u_{\tau h}^{k+1},u_{\tau h}^{k}\rangle$ makes this discretization less suitable for the problem that we wish to consider in this paper where the stored energy is enhanced by some internal variables and (possibly) nonlinear processes on them.

Therefore, we use another explicit discretization scheme, the so-called leap-frog scheme. To this end, we first write the velocity/stress formulation of (1.1a). Introducing the velocity, $v=\mathchoice{{\buildrel\hskip 1.00006pt\text{\LARGE.}\over{u}}}{{\buildrel\hskip 1.00006pt\text{\Large.}\over{u}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{u}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{u}}}$ and the stress tensor $\sigma:=\mathbb{C}e(u)$, we get

		$\displaystyle\varrho\mathchoice{{\buildrel\hskip 1.00006pt\text{\LARGE.}\over{v}}}{{\buildrel\hskip 1.00006pt\text{\Large.}\over{v}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{v}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{v}}}-{\rm div}\,\sigma=f\ \ \ \text{ and }\ \ \ \mathchoice{{\buildrel\hskip 1.00006pt\text{\LARGE.}\over{\sigma}}}{{\buildrel\hskip 1.00006pt\text{\Large.}\over{\sigma}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{\sigma}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{\sigma}}}=\mathbb{C}e(v)$	$\displaystyle\text{on }\ \varOmega\ \text{ for }\ t\in(0,T],$		(1.5a)
		$\displaystyle\mathchoice{{\buildrel\hskip 1.00006pt\text{\LARGE.}\over{\sigma}}}{{\buildrel\hskip 1.00006pt\text{\Large.}\over{\sigma}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{\sigma}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{\sigma}}}\vec{n}+\mathbb{B}v=\mathchoice{{\buildrel\hskip 1.00006pt\text{\LARGE.}\over{g}}}{{\buildrel\hskip 1.00006pt\text{\Large.}\over{g}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{g}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{g}}}$	$\displaystyle\text{on }\ \varGamma\ \text{ for }\ t\in(0,T],$		(1.5b)
		$\displaystyle v|_{t=0}=v_{0},\ \ \ \sigma|_{t=0}=\sigma_{0}:=\mathbb{C}e(u_{0})$	$\displaystyle\text{on }\ \varOmega.$		(1.5c)

In the abstract form (1.2), when writing $\mathcal{W}=\mathscr{W}\circ E$ with $E$ denoting the linear operator $u\mapsto(e,w):=(e(u),u|_{\varGamma})$ and with $u|_{\varGamma}$ denoting the trace of $u$ on the boundary $\varGamma$, this reads as

		$\displaystyle\mathcal{T}^{\prime}\mathchoice{{\buildrel\hskip 1.00006pt\text{\LARGE.}\over{v}}}{{\buildrel\hskip 1.00006pt\text{\Large.}\over{v}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{v}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{v}}}+E^{*}\varSigma=F(t)$	$\displaystyle\text{ for }\ t\in(0,T],$	$\displaystyle v|_{t=0}=v_{0},\ \ \ \ \text{ and }\ \$		(1.6a)
		$\displaystyle\mathchoice{{\buildrel\hskip 1.00006pt\text{\LARGE.}\over{\varSigma}}}{{\buildrel\hskip 1.00006pt\text{\Large.}\over{\varSigma}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{\varSigma}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{\varSigma}}}=\mathscr{W}^{\prime}Ev+\mathchoice{{\buildrel\hskip 1.00006pt\text{\LARGE.}\over{G}}}{{\buildrel\hskip 1.00006pt\text{\Large.}\over{G}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{G}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{G}}}(t)$	$\displaystyle\text{ for }\ t\in(0,T],$	$\displaystyle\varSigma|_{t=0}=\varSigma_{0}:=\mathscr{W}^{\prime}Eu_{0},$		(1.6b)

where $E^{*}$ is the adjoint operator to $E$. The stored energy governing (1.5) is

$$\mathscr{W}(e,w)=\int_{\varOmega}\frac{1}{2}\mathbb{C}e{:}e\,\mathrm{d}x+\int_{\varGamma}\frac{1}{2}\mathbb{B}w{\cdot}w\,\mathrm{d}S$$

while the external loading is now split into two parts acting differently, namely

$$\langle F(t),u\rangle=\int_{\varOmega}f(t)\cdot u\,\mathrm{d}x\qquad\text{and}\qquad\langle G(t),w\rangle=\int_{\varGamma}g(t)\cdot w\,\mathrm{d}x.$$

Let us note that (1.6) involves the equation on $\varOmega$, as well as, the equation on $\varGamma$. Thus $\mathcal{T}$ is to be understood as the functional on $\varOmega\times\varGamma$, that is trivial on $\varGamma$ since no inertia is considered on the $(d{-}1$)-dimensional boundary $\varGamma$. In particular, the “generalized” stress $\varSigma=\mathscr{W}^{\prime}Eu=(\mathbb{C}e(u),\mathbb{B}u|_{\varGamma})$ contains, besides the bulk stress tensor, also the traction stress vector. Relying on the linearity of $\mathscr{W}^{\prime}$, we have $\mathchoice{{\buildrel\hskip 1.00006pt\text{\LARGE.}\over{\varSigma}}}{{\buildrel\hskip 1.00006pt\text{\Large.}\over{\varSigma}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{\varSigma}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{\varSigma}}}=\mathscr{W}^{\prime}Ev$ with $v=\mathchoice{{\buildrel\hskip 1.00006pt\text{\LARGE.}\over{u}}}{{\buildrel\hskip 1.00006pt\text{\Large.}\over{u}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{u}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{u}}}$, as used in (1.6b). Let us note that the adjoint operator $E^{*}:(\sigma,\varsigma)\mapsto\mathfrak{F}$ in (1.6a) with the traction force $\varsigma=\mathbb{B}u|_{\varGamma}$ determines a bulk force $\mathfrak{F}$ as a functional on test displacements $u$ by

	$\displaystyle\int_{\varOmega}\mathfrak{F}\cdot u\,\mathrm{d}x=\left\langle\mathfrak{F},u\right\rangle$	$\displaystyle=\left\langle E^{*}(\sigma,\varsigma),u\right\rangle=\left\langle(\sigma,\varsigma),Eu\right\rangle=\left\langle(\sigma,\varsigma),(e(u),u|_{\varGamma})\right\rangle$
		$\displaystyle=\int_{\varOmega}\sigma:e(u)\,\mathrm{d}x+\int_{\varGamma}\varsigma\cdot u\,\mathrm{d}S=\int_{\varGamma}(\sigma\cdot\vec{n}+\varsigma)\cdot u\,\mathrm{d}S-\int_{\varOmega}{\rm div}\,\sigma\cdot u\,\mathrm{d}x\,,$

which clarifies the force $E^{*}\varSigma=\mathfrak{F}$ in (1.6a). The leap-frog time discretization of (1.6) then reads as

$\displaystyle\frac{\varSigma_{\tau h}^{k+1/2}-\varSigma_{\tau h}^{k-1/2}}{\tau}=\mathscr{W}_{h}^{\prime}E_{h}v_{\tau h}^{k}+D_{\tau h}^{k}\ \ \ \ \ \text{ and }\ \ \ \ {\mathcal{T}}^{\prime}\frac{v_{\tau h}^{k+1}{-}v_{\tau h}^{k}}{\tau}+E_{h}^{*}\varSigma_{\tau h}^{k+1/2}=F_{\tau h}^{k+1/2},$

(1.7)

where $\mathscr{W}_{h}$ and $E_{h}$ are suitable FEM discretizations of $\mathscr{W}$ and $E$, cf. (3.1a) and (3.1) below, and

$\displaystyle F_{\tau h}^{k+1/2}:=\frac{1}{\tau}\int_{k\tau}^{(k+1)\tau}\!\!\!\!\!\!\!F_{h}(t)\,\mathrm{d}t\ \ \ \text{ and }\ \ \ D_{\tau h}^{k}:=\frac{1}{\tau}\int_{(k-1/2)\tau}^{(k+1/2)\tau}\!\!\mathchoice{{\buildrel\hskip 1.00006pt\text{\LARGE.}\over{G}}}{{\buildrel\hskip 1.00006pt\text{\Large.}\over{G}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{G}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{G}}}_{h}(t)\,\mathrm{d}t=\frac{G_{h}((k{+}\frac{1}{2})\tau)-G_{h}((k{-}\frac{1}{2})\tau)}{\tau}\,.$

(1.8)

As mentioned before, we assume here that $v$ is discretized in an element-wise constant way so that $\mathcal{T}$ leads to a diagonal matrix. In this case we do not need to employ numerical integration to approximate the mass matrix. For higher-order discretizations, however, mass lumping is necessary so as to obtain explicit discretization schemes. We refer to [10, 31, 50] for details in the case $G\equiv 0$. The proposed FEM leads to a block diagonal matrix for $\mathcal{W}_{h}^{\prime}=E_{h}^{*}\mathscr{W}_{h}^{\prime}E_{h}$, which means that the resulting scheme does not require the solution of a big linear system at each iteration in time. The spatial FEM discretization exploits regularity available in linear elastodynamics, in particular that ${\rm div}\,\sigma$ and $e(v)$ in (1.5a) live in $L^{2}$-spaces. Moreover, the equations in (1.7) are decoupled in the sense that, first, $\varSigma_{\tau h}^{k+1/2}$ is calculated from the former equation and, second, $v_{\tau h}^{k+1}$ is calculated from the latter equation assuming, that $(v_{\tau h}^{k},\varSigma_{\tau h}^{k-1/2})$ is known from the previous time step. For $k=0$, it starts from $v_{\tau h}^{0}=v_{0}$ and from a half time step $\varSigma_{\tau h}^{1/2}=\varSigma_{\tau h}^{0}+\frac{\tau}{2}\mathscr{W}_{h}^{\prime}E_{h}v_{\tau h}^{0}$. For the space discretization, the lower order $Q^{\rm div}_{k+1}-Q_{k}$ finite element is obtained for $k=0$ and in this case the velocity is discretized as piecewise constant on rectangular or cubic elements while the stress is discretized by piecewise bi-linear functions with some continuities. Namely the normal component of the stress is continuous across edges of adjacent elements while the tangential component is allowed to be discontinuous. For more details about the space discretization we refer the interested reader to [10]. Alternative discretizations for the linear elasticity problem have been proposed by D. Arnold and his collaborators who designed mixed finite elements for general rectangular and triangular grids [6, 4, 5]. For tetrahedral leap-frog discretization of the elastodynamics see [21]. In general, the leap-frog scheme has been frequently used in geophysics to calculate seismic wave propagation with the finite differences method, cf. e.g. [13, 25, 51].

When taking the average (i.e. the sum with the weights $\frac{1}{2}$ and $\frac{1}{2}$) of the second equation in (1.7) in the level $k$ and $k{-}1$ tested by $v_{\tau h}^{k}$ and summing it with the first equation in (1.7) tested by $[\mathscr{W}_{h}^{\prime}]^{-1}(\varSigma_{\tau h}^{k+1/2}{+}\varSigma_{\tau h}^{k-1/2})/2$, we obtain

	$\displaystyle\frac{1}{2\tau}\big{\langle}[\mathscr{W}_{h}^{\prime}]^{-1}\varSigma_{\tau h}^{k+1/2},\varSigma_{\tau h}^{k+1/2}\big{\rangle}-\frac{1}{2\tau}\big{\langle}[\mathscr{W}_{h}^{\prime}]^{-1}\varSigma_{\tau h}^{k-1/2},\varSigma_{\tau h}^{k-1/2}\big{\rangle}$
	$\displaystyle\qquad\qquad=\Big{\langle}\frac{\varSigma_{\tau h}^{k+1/2}{+}\varSigma_{\tau h}^{k-1/2}}{2},E_{h}v_{\tau h}^{k}\Big{\rangle}+\Big{\langle}[\mathscr{W}_{h}^{\prime}]^{-1}D_{\tau h}^{k},\frac{\varSigma_{\tau h}^{k+1/2}{+}\varSigma_{\tau h}^{k-1/2}}{2}\Big{\rangle}\ \ \text{ and}$
	$\displaystyle\Big{\langle}{\mathcal{T}}^{\prime}\frac{v_{\tau h}^{k+1}{-}v_{\tau h}^{k-1}}{\tau},v_{\tau h}^{k}\Big{\rangle}+\Big{\langle}\frac{\varSigma_{\tau h}^{k+1/2}{+}\varSigma_{\tau h}^{k-1/2}}{2},E_{h}v_{\tau h}^{k}\Big{\rangle}=\Big{\langle}\frac{F_{\tau h}^{k+1/2}{+}F_{\tau h}^{k-1/2}}{2},v_{\tau h}^{k}\Big{\rangle}\,.$

Summing it up, we get that the following discrete energy is conserved

$\displaystyle\frac{1}{2}\langle\color[rgb]{0,0,0}{\mathcal{T}}{}^{\prime}v_{\tau h}^{k+1},v_{\tau h}^{k}\rangle+\varPhi_{h}(\varSigma_{\tau h}^{k+1/2})\ \ \text{ with }\ \ \varPhi_{h}(\varSigma)=\frac{1}{2}\langle[\mathscr{W}_{h}^{\prime}]^{-1}\varSigma,\varSigma\rangle\,.$

(1.9)

Note that $\varPhi_{h}$ is the discrete stored energy expressed in terms of the generalized stress. In contrast to (1.3), this formulation allows for enhancement of the discrete stored energy by some internal variables. The energy (1.9) is shown to be a positive quantity under the following CFL condition

$\displaystyle\big{\langle}[\mathcal{W}^{\prime}_{h}]^{-1}\varSigma_{h},\varSigma_{h}\big{\rangle}\geq\frac{\tau^{2}}{4}\big{\langle}E^{*}_{h}\varSigma_{h},({\mathcal{T}}^{\prime})^{-1}E^{*}_{h}\varSigma_{h}\big{\rangle}$

(1.10)

for any $\varSigma_{h}$ from the respective finite-dimensional subspace. Moreover, $F=0$ is often considered, which makes the a-priori estimation easier. Let us also note that the adjective “leap-frog” is sometimes used also for the time-discretization (1.3) if written as a two-step scheme, cf. e.g. [19, Sect. 7.1.1.1].

The plan of this article is as follows: In Section 2, we complement the abstract system (1.6) by another equation for some internal variable and cast its weak formulation without relying on any regularity. Then, in Section 3, we extend the two-step leap-frog discrete scheme (1.7) to a suitable three-step scheme, and study the energy properties of the proposed scheme. Then, in Section 4, we prove the numerical stability of the 3-step staggered approximation scheme and its convergence under the CFL condition (4.1). Such an abstract scheme is then illustrated in Section 5 on several examples from continuum mechanics, in particular on models of plasticity, creep, diffusion, and damage. For illustration of computational efficiency, we refer to [43] where this scheme was implemented for another problem, namely a delamination (i.e. interfacial damage).

It should be emphasized that, to the best of our knowledge, a rigorously justified (as far as numerical stability and convergence) combination of the explicit staggered discretization with nonlinear dissipative processes on some internal variables is new, although occasionally some dissipative nonlinear phenomena can be found in literature as in [45] for a unilateral contact, in [13] for a Maxwell viscoelastic rheology, in [47] for electroactive polymers, in [23] for an aeroelastic system, or in [24] for general thermomechanical systems, but without any numerical stability (a-priori estimates) and convergence guaranteed.

The concept of internal variables has a long tradition, cf. [36], and opens wide options for material modelling, cf. e.g. [35, 37] and references therein. Typically, internal variables are governed by a parabolic-type 1st-order evolution. The abstract system (1.2) is thus generalized to

		$\displaystyle\mathcal{T}^{\prime}\mathchoice{{\buildrel\hskip 1.00006pt\text{\LARGE..}\over{u}}}{{\buildrel\hskip 1.00006pt\text{\Large..}\over{u}}}{{\buildrel\hskip 1.00006pt\text{\large..}\over{u}}}{{\buildrel\hskip 1.00006pt\text{\large..}\over{u}}}+\mathcal{W}_{u}^{\prime}(u,z)=F(t)$	$\displaystyle\text{ for }\ t\in(0,T],\ \ u|_{t=0}=u_{0},\ \ \ \mathchoice{{\buildrel\hskip 1.00006pt\text{\LARGE.}\over{u}}}{{\buildrel\hskip 1.00006pt\text{\Large.}\over{u}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{u}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{u}}}|_{t=0}=v_{0},$		(2.1a)
		$\displaystyle\partial\varPsi(\mathchoice{{\buildrel\hskip 1.00006pt\text{\LARGE.}\over{z}}}{{\buildrel\hskip 1.00006pt\text{\Large.}\over{z}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{z}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{z}}})+\mathcal{W}_{z}^{\prime}(u,z)\ni 0$	$\displaystyle\text{ for }\ t\in(0,T],\ \ z|_{t=0}=z_{0}.$		(2.1b)

The inclusion in (2.1b) refers to a possibility that the convex (pseudo)potential $\varPsi$ of dissipative forces may be nonsmooth and then its subdifferential $\partial\varPsi$ can be multivalued.

Combination of the 2nd-order evolution (1.2) with such 1st-order evolution is to be handled carefully. In contrast to the implicit staggered schemes, cf. [42], the constitutive equation is differentiated in time, cf. (1.5a), and it seems necessary to use the staggered scheme so that the internal-variable flow rule can be used without being differentiated in time, even for a quadratic stored energy $\mathcal{W}$.

Moreover, to imitate the leap-frog scheme, it seems suitable (or maybe even necessary) that the stored energy $\mathcal{W}$ may be expressed in terms of the generalized stress as

$\displaystyle\mathcal{W}(u,z)=\varPhi(\varSigma,z)\ \ \text{ with }\ \ \varSigma=\mathfrak{C}Eu,\ \text{ and }\ \ \varPhi(\cdot,z)\ \text{ and }\ \varPhi(\varSigma,\cdot)\text{ quadratic},$

(2.2)

where $\mathfrak{C}$ stands for a “generalized” elasticity tensor and $E$ is an abstract gradient-type operator. Typically $Eu=(e(u),u|_{\varGamma})$ or simply $Eu=e(u)$ are considered here in the context of continuum mechanics at small strains, cf. the examples in Sect. 5. Here, $\varSigma$ may not directly enter the balance of forces and is thus to be called rather as some “proto-stress”, while the actual generalized stress will be denoted by $S$. For a relaxation of the last requirement of (2.2) see Remark 4.4 below.

Then, likewise (1.6), we can write the system (2.1) in the velocity/proto-stress formulation as

		$\displaystyle\mathchoice{{\buildrel\hskip 1.00006pt\text{\LARGE.}\over{\varSigma}}}{{\buildrel\hskip 1.00006pt\text{\Large.}\over{\varSigma}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{\varSigma}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{\varSigma}}}=\mathfrak{C}Ev+\mathchoice{{\buildrel\hskip 1.00006pt\text{\LARGE.}\over{G}}}{{\buildrel\hskip 1.00006pt\text{\Large.}\over{G}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{G}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{G}}}(t)$	$\displaystyle\text{ for }\ t\in(0,T],\ \$		(2.3a)
		$\displaystyle\mathcal{T}^{\prime}\mathchoice{{\buildrel\hskip 1.00006pt\text{\LARGE.}\over{v}}}{{\buildrel\hskip 1.00006pt\text{\Large.}\over{v}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{v}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{v}}}+E^{*}S=F(t)\ \ \text{ with }\ \ S=\mathfrak{C}^{*}\varPhi_{\varSigma}^{\prime}(\varSigma,z)\!\!\!\!\!\!$	$\displaystyle\text{ for }\ t\in(0,T],\ \$		(2.3b)
		$\displaystyle\partial\varPsi(\mathchoice{{\buildrel\hskip 1.00006pt\text{\LARGE.}\over{z}}}{{\buildrel\hskip 1.00006pt\text{\Large.}\over{z}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{z}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{z}}})+\varPhi_{z}^{\prime}(\varSigma,z)\ni 0$	$\displaystyle\text{ for }\ t\in(0,T],$		(2.3c)
		$\displaystyle\varSigma|_{t=0}=\varSigma_{0}:=\mathfrak{C}Eu_{0}+G(0)\,,\ \ \ v|_{t=0}=v_{0}\,,\ \ z|_{t=0}=z_{0}\,.$			(2.3d)

Here $\varPhi_{\varSigma}^{\prime}(\varSigma,z)$ is a “generalized” strain and, when multiplied by $\mathfrak{C}^{*}$, it becomes a generalized stress. Actually, (1.6) is a special case of (2.3) when $\varPhi=\varPhi(\varSigma)=\int_{\varOmega}\frac{1}{2}\mathbb{C}^{-1}\sigma:\sigma\,\mathrm{d}x+\int_{\varGamma}\mathbb{B}^{-1}\varsigma\cdot\varsigma\,\mathrm{d}S$ and $S=\varSigma=(\sigma,\varsigma)$ while $\mathfrak{C}=\mathscr{W}^{\prime}=(\mathbb{C},\mathbb{B})$, so that indeed $S=\mathfrak{C}^{*}\varPhi_{\varSigma}^{\prime}(\varSigma)=(\mathbb{C},\mathbb{B})(\mathbb{C}^{-1}\sigma,\mathbb{B}^{-1}\varsigma)=(\sigma,\varsigma)=\varSigma$.

The energy properties of this system can be revealed by testing the particular equations/inclusions in (2.3) by $\varPhi_{\varSigma}^{\prime}(\varSigma,z)$, $v$, and $\mathchoice{{\buildrel\hskip 1.00006pt\text{\LARGE.}\over{z}}}{{\buildrel\hskip 1.00006pt\text{\Large.}\over{z}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{z}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{z}}}$. Thus, at least formally, we obtain

	$\displaystyle\big{\langle}\varPhi_{\varSigma}^{\prime}(\varSigma,z),\mathchoice{{\buildrel\hskip 1.00006pt\text{\LARGE.}\over{\varSigma}}}{{\buildrel\hskip 1.00006pt\text{\Large.}\over{\varSigma}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{\varSigma}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{\varSigma}}}\big{\rangle}=\big{\langle}\varPhi_{\varSigma}^{\prime}(\varSigma,z),\mathfrak{C}Ev+\mathchoice{{\buildrel\hskip 1.00006pt\text{\LARGE.}\over{G}}}{{\buildrel\hskip 1.00006pt\text{\Large.}\over{G}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{G}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{G}}}\big{\rangle}=\big{\langle}\mathfrak{C}^{*}\varPhi_{\varSigma}^{\prime}(\varSigma,z),Ev\big{\rangle}+\big{\langle}\varPhi_{\varSigma}^{\prime}(\varSigma,z),\mathchoice{{\buildrel\hskip 1.00006pt\text{\LARGE.}\over{G}}}{{\buildrel\hskip 1.00006pt\text{\Large.}\over{G}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{G}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{G}}}\big{\rangle},$		(2.4a)
	$\displaystyle\color[rgb]{0,0,0}\big{\langle}{\mathcal{T}}^{\prime}\mathchoice{{\buildrel\hskip 1.00006pt\text{\LARGE.}\over{v}}}{{\buildrel\hskip 1.00006pt\text{\Large.}\over{v}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{v}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{v}}},v\big{\rangle}\color[rgb]{0,0,0}+\big{\langle}\mathfrak{C}^{*}\varPhi_{\varSigma}^{\prime}(\varSigma,z),Ev\big{\rangle}=\big{\langle}F(t),v\big{\rangle}\,,$		(2.4b)
	$\displaystyle\varXi(\mathchoice{{\buildrel\hskip 1.00006pt\text{\LARGE.}\over{z}}}{{\buildrel\hskip 1.00006pt\text{\Large.}\over{z}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{z}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{z}}})+\big{\langle}\varPhi_{z}^{\prime}(\varSigma,z),\mathchoice{{\buildrel\hskip 1.00006pt\text{\LARGE.}\over{z}}}{{\buildrel\hskip 1.00006pt\text{\Large.}\over{z}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{z}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{z}}}\big{\rangle}\leq 0\ \ \ \ \text{ with }\ \ \ \varXi(\mathchoice{{\buildrel\hskip 1.00006pt\text{\LARGE.}\over{z}}}{{\buildrel\hskip 1.00006pt\text{\Large.}\over{z}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{z}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{z}}}):=\inf\big{\langle}\partial\varPsi(\mathchoice{{\buildrel\hskip 1.00006pt\text{\LARGE.}\over{z}}}{{\buildrel\hskip 1.00006pt\text{\Large.}\over{z}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{z}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{z}}}),\mathchoice{{\buildrel\hskip 1.00006pt\text{\LARGE.}\over{z}}}{{\buildrel\hskip 1.00006pt\text{\Large.}\over{z}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{z}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{z}}}\big{\rangle}\,.$		(2.4c)

The functional $\varXi$ is a dissipative rate and the “inf” in it refers to the fact that the dissipative potential $\varPsi$ can be nonsmooth and thus the subdifferential $\partial\varPsi$ can be multivalued even at $\mathchoice{{\buildrel\hskip 1.00006pt\text{\LARGE.}\over{z}}}{{\buildrel\hskip 1.00006pt\text{\Large.}\over{z}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{z}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{z}}}\neq 0$, otherwise an equality in (2.4c) holds. Summing it up and using the calculus $\frac{\mathrm{d}}{\mathrm{d}t}{\mathcal{T}}(v)=\langle{\mathcal{T}}^{\prime}v,\mathchoice{{\buildrel\hskip 1.00006pt\text{\LARGE.}\over{v}}}{{\buildrel\hskip 1.00006pt\text{\Large.}\over{v}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{v}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{v}}}\rangle=\langle{\mathcal{T}}^{\prime}\mathchoice{{\buildrel\hskip 1.00006pt\text{\LARGE.}\over{v}}}{{\buildrel\hskip 1.00006pt\text{\Large.}\over{v}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{v}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{v}}},v\rangle$ and $\frac{\mathrm{d}}{\mathrm{d}t}\varPhi(\varSigma,z)=\langle\varPhi_{\varSigma}^{\prime}(\varSigma,z),\mathchoice{{\buildrel\hskip 1.00006pt\text{\LARGE.}\over{\varSigma}}}{{\buildrel\hskip 1.00006pt\text{\Large.}\over{\varSigma}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{\varSigma}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{\varSigma}}}\rangle+\langle\varPhi_{z}^{\prime}(\varSigma,z),\mathchoice{{\buildrel\hskip 1.00006pt\text{\LARGE.}\over{z}}}{{\buildrel\hskip 1.00006pt\text{\Large.}\over{z}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{z}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{z}}}\rangle$, we obtain the following inequality for the energy,

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\big{(}\!\!\!\!\begin{array}[t]{c}\begin{array}[t]{c}\underbrace{\mathcal{T}(v)+\varPhi(\varSigma,z)}\vspace*{.5em}\end{array}\vspace*{-1em}\\ _{\mbox{\footnotesize\rm kinetic and}}\vspace*{-.5em}\\ _{\mbox{\footnotesize\rm stored energies}}\end{array}\!\!\!\!\big{)}\ +\!\!\!\!\!\!\begin{array}[t]{c}\begin{array}[t]{c}\underbrace{\varXi(\mathchoice{{\buildrel\hskip 1.00006pt\text{\LARGE.}\over{z}}}{{\buildrel\hskip 1.00006pt\text{\Large.}\over{z}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{z}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{z}}})}\vspace*{.5em}\end{array}\vspace*{-1em}\\ _{\mbox{\footnotesize\rm dissipation}}\vspace*{-.5em}\\ _{\mbox{\footnotesize\rm rate}}\end{array}\!\!\!\!\!\!\color[rgb]{0,0,0}\leq\color[rgb]{0,0,0}\!\!\begin{array}[t]{c}\begin{array}[t]{c}\underbrace{\big{\langle}F(t),v\big{\rangle}+\big{\langle}\varPhi_{\varSigma}^{\prime}(\varSigma,z),\mathchoice{{\buildrel\hskip 1.00006pt\text{\LARGE.}\over{G}}}{{\buildrel\hskip 1.00006pt\text{\Large.}\over{G}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{G}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{G}}}\big{\rangle}}\vspace*{.5em}\end{array}\vspace*{-1em}\\ _{\mbox{\footnotesize\rm power of}}\vspace*{-.5em}\\ _{\mbox{\footnotesize\rm external force}}\end{array}\,.$

(2.17)

Actually, (2.4c) and (2.17) often hold as equalities.

Let us now formulate some abstract functional setting of the system (2.3). For some Banach spaces $\mathcal{S}$, $\mathcal{Z}$, and $\mathcal{Z}_{1}\supset\mathcal{Z}$ and for a Hilbert space $\mathcal{H}$, let $\varPhi:{\mathcal{S}}\times\mathcal{Z}\to{\mathbb{R}}$ be smooth and coercive, $\mathcal{T}:\mathcal{H}\to{\mathbb{R}}$ be quadratic and coercive, and let $\Psi:\mathcal{Z}\to[0,+\infty]$ be convex, lower semicontinuous, and coercive on $\mathcal{Z}_{1}$, cf. (4.2) below. Intentionally, we do not want to rely on any regularity which is usually at disposal in linear problems but might be restrictive in some nonlinear problems. For this reason, we reconstruct the abstract “displacement” and use (2.3a) integrated in time, i.e.

$\displaystyle\varSigma={\mathfrak{C}}Eu+G\ \ \ \text{ with }\ \ u(t):=\int_{0}^{t}\!\!v(t)\,\mathrm{d}t+u_{0}\,.$

(2.18)

Moreover, we still need another Banach space $\mathcal{E}$ and define the Banach space $\mathcal{U}:=\{u\in\mathcal{H};\ Eu\in\mathcal{E}\}$ equipped with the standard graph norm. Then, by definition, we have the continuous embedding $\mathcal{U}\to\mathcal{H}$ and the continuous linear operator $E:\mathcal{U}\to\mathcal{E}$. We assume that $\mathcal{U}$ is embedded into $\mathcal{H}$ densely, so that $\mathcal{H}^{*}\subset\mathcal{U}^{*}$ and that $\mathcal{H}$ is identified with its dual $\mathcal{H}^{*}$, so that we have the so-called Gelfand triple

$$\mathcal{U}\subset\mathcal{H}\cong\mathcal{H}^{*}\subset\mathcal{U}^{*}.$$

We further consider the abstract elasticity tensor $\mathfrak{C}$ as a linear continuous operator $\mathcal{E}\to\mathcal{S}$. Therefore $\mathfrak{C}Eu\in\mathcal{S}$ provided $u\in\mathcal{U}$ so that the equation (2.18) is meant in $\mathcal{S}$ and one needs $G(t)\in\mathcal{S}$. Let us note that $\mathcal{T}^{\prime}:\mathcal{H}\to\mathcal{H}^{*}\cong\mathcal{H}$, $\varPhi_{\varSigma}^{\prime}:{\mathcal{S}}\times\mathcal{Z}\to{\mathcal{S}}^{*}$, $E^{*}:\mathcal{E}^{*}\to\mathcal{U}^{*}$, and $\mathfrak{C}^{*}:\mathcal{S}^{*}\to\mathcal{E}^{*}$, so that $\mathcal{T}^{\prime}v\in\mathcal{H}^{*}$ provided $v\in\mathcal{H}$ and also $S=\mathfrak{C}^{*}\varPhi_{\varSigma}^{\prime}\in\mathcal{E}^{*}$ and $E^{*}S\in\mathcal{H}^{*}$. In particular, the equation (2.3b) can be meant in $\mathcal{H}$ if integrated in time, and one needs $F(t)$ valued in $\mathcal{H}$.

For a Banach space $\mathcal{X}$, we will use the standard notation $L^{p}(0,T;\mathcal{X})$ for Bochner spaces of the Bochner measurable functions $[0,T]\to\mathcal{X}$ whose norm is integrable with the power $p$ or essentially bounded if $p=\infty$, and $W^{1,p}(0,T;\mathcal{X})$ the space of functions from $L^{p}(0,T;\mathcal{X})$ whose distributional time derivative is also in $L^{p}(0,T;\mathcal{X})$. Also, $C^{k}(0,T;\mathcal{X})$ will denote the space of functions $[0,T]\to\mathcal{X}$ whose $k^{\rm th}$-derivative is continuous, and $C_{\rm w}(0,T;\mathcal{X})$ will denote the space of weakly continuous functions $[0,T]\to\mathcal{X}$. Later, we will also use ${\rm Lin}(\mathcal{U},\mathcal{E})$, denoting the space of linear bounded operators $\mathcal{U}\to\mathcal{E}$ normed by the usual sup-norm.

A weak formulation of (2.3b) can be obtained after by-part integration over the time interval $I=[0,T]$ when tested by a smooth function. It is often useful to consider

$\displaystyle\varPhi(\varSigma,z)=\varPhi_{0}(\varSigma,z)+\varPhi_{1}(z)\ \ \ \text{ with }\ [\varPhi_{0}]_{z}^{\prime}:{\mathcal{S}}\times\mathcal{Z}\to\mathcal{Z}_{1}^{*}\ \text{ and }\ \varPhi_{1}^{\prime}:\mathcal{Z}\to\mathcal{Z}^{*}$

(2.19)

and to use integration by-parts for the term $\langle\varPhi_{1}^{\prime}(z),\mathchoice{{\buildrel\hskip 1.00006pt\text{\LARGE.}\over{z}}}{{\buildrel\hskip 1.00006pt\text{\Large.}\over{z}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{z}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{z}}}\rangle$. We thus arrive to the following definition.

Let us note that the remaining initial condition $v(0)=v_{0}$ is contained in (2.20a). Definition 2.1 works successfully for $p>1$, i.e. for rate-dependent evolution of the abstract internal variable $z$, so that $\mathchoice{{\buildrel\hskip 1.00006pt\text{\LARGE.}\over{z}}}{{\buildrel\hskip 1.00006pt\text{\Large.}\over{z}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{z}}}{{\buildrel\hskip 1.00006pt\text{\large.}\over{z}}}\in L^{p}(0,T;\mathcal{Z}_{1})$. For the rate-dependent evolution when $p=1$, we would need to modify it. Here, we restrict ourselves to $p\geq 2$, because of the a-priori estimates in Proposition 4.1.

We derive in this section the leap-frog discretization of (2.3a,b) combined with an implicit discretization for (2.3c), using a fractional-step split (called also a staggered scheme) with a mid-point formula for (2.3c). Instead of a two-step scheme (1.7), we obtain a three-step scheme and therefore, from now on, we abandon the convention of a half-step notation as used in (1.7) and write $k+1$ instead of $k+1/2$.

To this aim, we consider sequences of nested finite-dimensional subspaces $S_{h}\subset\mathcal{S}$, $V_{h}\subset\mathcal{H}$, and $Z_{h}\subset\mathcal{Z}$ where the values of the respective discrete variables $\varSigma_{h}$, $v_{h}$, and $z_{h}$ will be, assuming that their unions are dense in the respective Banach spaces. We will use an interpolation operator $I_{h}:{\rm Lin}(\mathcal{S},S_{h})$ and the embedding operator $J_{h}:Z_{h}\to\mathcal{Z}$; it is important that the collection $\{J_{h}\}_{h>0}$ is uniformly bounded and, since $\bigcup_{h>0}Z_{h}$ is dense in $\mathcal{Z}$, the sequence $\{J_{h}\}_{h>0}$ converges to the identity on $\mathcal{Z}$ strongly. We consider $E_{h}\in{\rm Lin}(V_{h},\mathcal{E})$. Let us note that we allow for a “non-conformal” approximation of $v$, i.e., $V_{h}\subset\mathcal{H}$ is not necessarily a subspace of $\mathcal{U}$. This is in agreement with discretizations of the velocity as in [10, 9, 11, 18, 50].

Considering that we know from previous step $\varSigma_{\tau h}^{k},v_{\tau h}^{k},z_{\tau h}^{k}$, then the proposed discretization scheme is

		1) calculate $\varSigma_{\tau h}^{k+1}$:	$\displaystyle\!\!\!\frac{\varSigma_{\tau h}^{k+1}-\varSigma_{\tau h}^{k}}{\tau}=I_{h}\mathfrak{C}E_{h}v_{\tau h}^{k}+D_{\tau h}^{k},$			(3.1a)
		2) calculate $z_{\tau h}^{k+1}$:	$\displaystyle\!\!\!J_{h}^{*}\partial\varPsi\Big{(}\frac{z_{\tau h}^{k+1}{-}z_{\tau h}^{k}}{\tau}\Big{)}+\varPhi_{z}^{\prime}\Big{(}\varSigma_{\tau h}^{k+1},\frac{z_{\tau h}^{k+1}{+}z_{\tau h}^{k}}{2}\Big{)}\ni 0\,,$			(3.1b)
		3) calculate $v_{\tau h}^{k+1}$:	$\displaystyle\!\!\!{\mathcal{T}}^{\prime}\frac{v_{\tau h}^{k+1}{-}v_{\tau h}^{k}}{\tau}+E_{h}^{\ast}S_{\tau h}^{k+1}=F_{\tau h}^{k+1/2}\ \text{ with }\ S_{\tau h}^{k+1}=\mathfrak{C}^{*}I_{h}^{*}\varPhi_{\varSigma}^{\prime}(\varSigma_{\tau h}^{k+1},z_{\tau h}^{k+1}),$
		and  $u_{\tau h}^{k+1}$:	$\displaystyle\!\!\!u_{\tau h}^{k+1}=u_{\tau h}^{k}+\tau v_{\tau h}^{k+1},$			(3.1c)

where $F_{\tau h}^{k+1/2}$ and $D_{\tau h}^{k}$ are from (1.8). It seems important in the non-linear case to compute the variables in the order given above. We note however, that for the linear viscoelastic problem with Maxwell rheology a scheme with a different ordering has been proposed in [27, Part I, Sect.2]. The potentials $\varPhi$, $\varPsi$, and ${\mathcal{T}}$ are considered restricted on $S_{h}\times Z_{h}$, $Z_{h}$, and $V_{h}$, so that their corresponding (sub)differentials $\varPhi_{\varSigma}^{\prime}$, $\varPhi_{z}^{\prime}$, $\partial\varPsi$, and ${\mathcal{T}}^{\prime}$ are valued in $S_{h}^{*}$, $Z_{h}^{*}$, and $V_{h}^{*}$, respectively. The particular equations/inclusion in (3.1) are thus to be understood in $S_{h}$, $Z_{h}^{*}$, $V_{h}^{*}$, $\mathcal{E}^{*}$, and $V_{h}$, respectively. The only implicit equation is (3.1b). Note however that even this equation becomes explicit if there are no spatial gradients in $\varPhi$ and $\varPsi$. In view of the definition of the convex subdifferential, (3.1b) means the variational inequality

$\displaystyle\varPsi\big{(}\widetilde{z}\big{)}+\bigg{\langle}\varPhi_{z}^{\prime}\Big{(}\varSigma_{\tau h}^{k+1},\frac{z_{\tau h}^{k+1}{+}z_{\tau h}^{k}}{2}\Big{)},\widetilde{z}-\frac{z_{\tau h}^{k+1}{-}z_{\tau h}^{k}}{\tau}\bigg{\rangle}\geq\varPsi\Big{(}\frac{z_{\tau h}^{k+1}{-}z_{\tau h}^{k}}{\tau}\Big{)}$

(3.2)

for any $\widetilde{z}\in V_{h}$. In any case, the equation (3.1b) posseses a potential

$\displaystyle z\mapsto\frac{2}{\tau}\varPhi\Big{(}\varSigma_{\tau h}^{k+1},\frac{z{+}z_{\tau h}^{k}}{2}\Big{)}+\varPsi\Big{(}\frac{z{-}z_{\tau h}^{k}}{\tau}\Big{)}$

(3.3)

which is to be minimized on $Z_{h}$. Therefore the existence of a solution to this inclusion (3.1b), or equivalently of the variational inequality (3.2), can be shown by a direct method, cf. also [42]. The scheme (3.1) is thus to be solved recurrently for $k=0,1,...,T/\tau-1$, starting from the initial conditions (2.3d) assumed, for simplicity, to live in the respective finite-dimensional spaces.

The energy properties of this scheme can be obtained by imitating (2.4)–(2.17). More specifically, we proceed as follows: we test (3.1a) by $\frac{1}{2}\varPhi_{\varSigma}^{\prime}(\varSigma_{\tau h}^{k+1},z_{\tau h}^{k+1})+\frac{1}{2}\varPhi_{\varSigma}^{\prime}(\varSigma_{\tau h}^{k},z_{\tau h}^{k})$, then test (3.1b) by $\frac{z_{\tau h}^{k+1}-z_{\tau h}^{k}}{\tau}$, and eventually test the average of (3.1) at the level $k{+}1$ and ${k}$ by $v_{\tau h}^{k}$. Using that $\varPhi(\cdot,z)$ and $\varPhi(\Sigma,\cdot)$ are quadratic as assumed in (2.2), we have

		$\displaystyle\bigg{\langle}\frac{\varPhi_{\varSigma}^{\prime}(\varSigma_{\tau h}^{k+1},z_{\tau h}^{k+1})+\varPhi_{\varSigma}^{\prime}(\varSigma_{\tau h}^{k},z_{\tau h}^{k})}{2},\frac{\varSigma_{\tau h}^{k+1}\!-\varSigma_{\tau h}^{k}}{\tau}\bigg{\rangle}$
		$\displaystyle\quad=\bigg{\langle}\frac{\varPhi_{\varSigma}^{\prime}(\varSigma_{\tau h}^{k+1},z_{\tau h}^{k})+\varPhi_{\varSigma}^{\prime}(\varSigma_{\tau h}^{k},z_{\tau h}^{k})}{2},\frac{\varSigma_{\tau h}^{k+1}\!-\varSigma_{\tau h}^{k}}{\tau}\bigg{\rangle}$
		$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad+\frac{\tau}{2}\bigg{\langle}\frac{\varPhi_{\varSigma}^{\prime}(\varSigma_{\tau h}^{k+1},z_{\tau h}^{k+1})-\varPhi_{\varSigma}^{\prime}(\varSigma_{\tau h}^{k+1},z_{\tau h}^{k})}{\tau},\frac{\varSigma_{\tau h}^{k+1}\!-\varSigma_{\tau h}^{k}}{\tau}\bigg{\rangle}$
		$\displaystyle\quad=\frac{\varPhi(\varSigma_{\tau h}^{k+1},z_{\tau h}^{k}){-}\varPhi(\varSigma_{\tau h}^{k},z_{\tau h}^{k})}{\tau}+\frac{\tau}{2}\bigg{\langle}\frac{\varPhi_{\varSigma}^{\prime}(\varSigma_{\tau h}^{k+1},z_{\tau h}^{k+1})-\varPhi_{\varSigma}^{\prime}(\varSigma_{\tau h}^{k+1},z_{\tau h}^{k})}{\tau},\frac{\varSigma_{\tau h}^{k+1}\!-\varSigma_{\tau h}^{k}}{\tau}\bigg{\rangle}\,,$		(3.4a)
where we used also (3.1a), and
		$\displaystyle\bigg{\langle}\varPhi_{z}^{\prime}\Big{(}\varSigma_{\tau h}^{k+1},\frac{z_{\tau h}^{k+1}{+}z_{\tau h}^{k}}{2}\Big{)},\frac{z_{\tau h}^{k+1}-z_{\tau h}^{k}}{\tau}\bigg{\rangle}=\frac{\varPhi(\varSigma_{\tau h}^{k+1},z_{\tau h}^{k+1})-\varPhi(\varSigma_{\tau h}^{k+1},z_{\tau h}^{k})}{\tau}\,.$		(3.4b)

Therefore, using again the particular equations/inclusion in (3.1), we get respectively

	$\displaystyle\frac{\varPhi(\varSigma_{\tau h}^{k+1},z_{\tau h}^{k})-\varPhi(\varSigma_{\tau h}^{k},z_{\tau h}^{k})}{\tau}=\bigg{\langle}\frac{\varPhi_{\varSigma}^{\prime}(\varSigma_{\tau h}^{k+1},z_{\tau h}^{k+1})+\varPhi_{\varSigma}^{\prime}(\varSigma_{\tau h}^{k},z_{\tau h}^{k})}{2},I_{h}\mathfrak{C}E_{h}v_{\tau h}^{k}+D_{\tau h}^{k}\bigg{\rangle}$
	$\displaystyle\hskip 140.00021pt-\frac{\tau}{2}\bigg{\langle}\frac{\varPhi_{\varSigma}^{\prime}(\varSigma_{\tau h}^{k+1},z_{\tau h}^{k+1})-\varPhi_{\varSigma}^{\prime}(\varSigma_{\tau h}^{k+1},z_{\tau h}^{k})}{\tau},\frac{\varSigma_{\tau h}^{k+1}\!-\varSigma_{\tau h}^{k}}{\tau}\bigg{\rangle}\,,$		(3.5a)
	$\displaystyle\varXi\Big{(}\frac{z_{\tau h}^{k+1}{-}z_{\tau h}^{k}}{\tau}\Big{)}+\frac{\varPhi(\varSigma_{\tau h}^{k+1},z_{\tau h}^{k+1})-\varPhi(\varSigma_{\tau h}^{k+1},z_{\tau h}^{k})}{\tau}\leq 0\,,\qquad\text{and}$		(3.5b)
	$\displaystyle\big{\langle}{\mathcal{T}}^{\prime}\frac{v_{\tau h}^{k+1}{-}v_{\tau h}^{k-1}}{2\tau},v_{\tau h}^{k}\big{\rangle}+\bigg{\langle}E_{h}^{*}\mathfrak{C}^{*}I_{h}^{*}\frac{\varPhi_{\varSigma}^{\prime}(\varSigma_{\tau h}^{k+1},z_{\tau h}^{k+1}){+}\varPhi_{\varSigma}^{\prime}(\varSigma_{\tau h}^{k},z_{\tau h}^{k})}{2},v_{\tau h}^{k}\bigg{\rangle}=\big{\langle}F_{\tau h}^{k+1/2},v_{\tau h}^{k}\big{\rangle}\,.$		(3.5c)

Let us also note that, if $\varPsi(0)=0$ is assumed, the substitution $\widetilde{z}=0$ into the inequality (3.2) gives $\varPsi(\frac{z_{\tau h}^{k+1}-z_{\tau h}^{k}}{\tau})$ instead of the dissipation rate $\varXi(\frac{z_{\tau h}^{k+1}-z_{\tau h}^{k}}{\tau})$ in (3.5b), which is a suboptimal estimate except if $\varPsi$ is degree-1 positively homogeneous.