Staggered explicit-implicit time-discretization for elastodynamics with dissipative internal variables
aa r X i v : . [ m a t h . NA ] J un Staggered explicit-implicit time-discretization for elastodynamicswith dissipative internal variables
Tom´aˇs Roub´ıˇcek ∗† and Chrysoula Tsogka ‡§ June 11, 2020
Abstract
An extension of the two-step staggered time discretization of linear elastodynamics in stress-velocity form to systems involving internal variables subjected to a possibly non-linear dissipativeevolution is proposed. The original scheme is thus enhanced by another step for the internalvariables which, in general, is implicit, although even this step might be explicit if no spatialgradients of the internal variables are involved. Using an abstract Banach-space formulation,a-priori estimates and convergence are proved under a CFL condition. The developed three-stepscheme finds applications in various problems of continuum mechanics at small strain. Here,we consider in particular plasticity, viscoelasticity (creep), diffusion in poroelastic media, anddamage.
AMS Subject Classification
Keywords
Elastodynamics, explicit discretization, fractional steps, mixed finite-element method,plasticity, creep, poro-elasticity, damage.
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 theenergy is not dominantly transmitted through space. In this case implicit time-discretization is rel-atively 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 ex-plicit time-discretizations are reasonably efficient, in particular in three-dimensions. These explicitmethods 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 needto combine the convervative hyperbolic problems with various dissipative processes of a paraboliccharacter, involving typically some internal variables. For parabolic problems, however, explicitmethods are known to be problematic due to severe time-step restrictions. Therefore, we propose ∗ Charles University, Mathematical Institute, Sokolovsk´a 83, CZ-186 75 Praha 8, Czech Republic † Institute of Thermomechanics of the Czech Acad. Sci., Dolejˇskova 5, CZ–182 08 Praha 8, Czech Rep. ‡ Department of Applied Mathematics, University of California, Merced, 5200 North Lake road, Merced, CA95343, USA § Institute of Applied and Computational Mathematics, Foundation for Research and Technology - Hellas, Niko-laou Plastira 100, Vassilika Vouton, GR-700 13 Heraklion, Crete, Greece u : [0 , T ] × Ω → R d satisfying ̺ .. u − div C e ( u ) = f on Ω for t ∈ (0 , T ] , (1.1a)[ C e ( u )] ~n + B u = g on Γ for t ∈ (0 , T ] , (1.1b) u | t =0 = u , . u | t =0 = v on Ω, (1.1c)for T > Ω ⊂ R d is a bounded Lipschitz domain, d = 2 or 3, Γ isits boundary, and ~n the unit outward normal. The dot-notation stands for the time derivative. In(1.1), ̺ > C is the elasticity tensor which is symmetric and positivedefinite, e ( u ) denotes the small-strain tensor defined as e ( u ) = ( ∇ u ) ⊤ + ∇ u , and B is a symmetricpositive semidefinite 2nd-order tensor determining the elastic support on the boundary. The termsappearing in the second member of (1.1) are, the bulk force f and the surface loading g . In (1.1c), u denotes the initial displacement and v the initial velocity. For a more compact notation, wewrite the initial-boundary-value problem (1.1) in the following abstract form T ′ .. u + W ′ u = F ′ u ( t ) for t ∈ (0 , T ] , u | t =0 = u , . u | t =0 = v . (1.2)Here T is the kinetic energy, W is the stored energy, and F is the external force, while ( · ) ′ denotesthe Gˆateaux derivative. In the context of (1.1), we have T ( . u ) = Z Ω ̺ | . u | d x, W ( u ) = Z Ω C e ( u ): e ( u ) d x + Z Γ B u · u d S and F ( t, u ) = Z Ω f ( t ) · u d x + Z Γ g ( t ) · u d S. Thus F ′ u ( t ) is a linear functional, let us denote it shortly by F ( t ).For high-frequency wave propagation problems, implicit time discretizations are computation-ally expensive, especially in the three dimensional setting. Therefore, we focus our attention onexplicit methods. The simplest explicit scheme is the following second-order finite differencescheme T ′ u k +1 τh − u kτh + u k − τh τ + W ′ h u kτh = F h ( kτ ) . (1.3)Here τ > W h (resp. F h ) denote some discrete approximations of therespective continuous functionals obtained by a suitable finite-element method (FEM) with meshsize h >
0. For simplicity, we assume the mass density constant (or at least piecewise constant inspace) so that the kinetic energy T does not need any numerical approximation. In particular, anumerical approximation leading to a diagonal the mass matrix T ′ in (1.3), typically referred to asmass lumping, is an important ingredient so as to obtain efficient explicit methods. Here we willconsider that u is discretized in an element-wise constant way so that T ′ leads to a diagonal matrix2ven without any approximation. Multiplying (1.3) by u k +1 τh − u kτh τ and using the scheme, it is easyto show that the energy preserved is12 hT ′ u k +1 τh − u kτh τ , u k +1 τh − u kτh τ i + 12 hW ′ h u k +1 τh , u kτh i . This is an energy, i.e. , a positive quantity under the following
Courant-Fridrichs-Lewy (CFL)condition [20] (cid:10) T ′ u h , u h (cid:11) ≥ τ (cid:10) W ′ h u h , u h (cid:11) (1.4)for any u h from the respective finite-dimensional subspace. The CFL typically bounds the timediscretization step τ = O ( h min ) with h min the smallest element size on a FEM discretization. Thismethod has frequently been used and analysed from various aspects, including comparison withimplicit time discretizations, cf. e.g. [33, 34]. However, the form of the discrete stored energy hW ′ h u k +1 τh , u kτh i makes this discretization less suitable for the problem that we wish to consider inthis paper where the stored energy is enhanced by some internal variables and (possibly) nonlinearprocesses on them.Therefore, we use another explicit discretization scheme, the so-called leap-frog scheme . To thisend, we first write the velocity/stress formulation of (1.1a). Introducing the velocity, v = . u andthe stress tensor σ := C e ( u ), we get ̺ . v − div σ = f and . σ = C e ( v ) on Ω for t ∈ (0 , T ] , (1.5a) . σ~n + B v = . g on Γ for t ∈ (0 , T ] , (1.5b) v | t =0 = v , σ | t =0 = σ := C e ( u ) on Ω. (1.5c)In the abstract form (1.2), when writing W = W ◦ E with E denoting the linear operator u ( e, w ) := ( e ( u ) , u | Γ ) and with u | Γ denoting the trace of u on the boundary Γ , this reads as T ′ . v + E ∗ Σ = F ( t ) for t ∈ (0 , T ] , v | t =0 = v , and (1.6a) . Σ = W ′ Ev + . G ( t ) for t ∈ (0 , T ] , Σ | t =0 = Σ := W ′ Eu , (1.6b)where E ∗ is the adjoint operator to E . The stored energy governing (1.5) is W ( e, w ) = Z Ω C e : e d x + Z Γ B w · w d S while the external loading is now split into two parts acting differently, namely h F ( t ) , u i = Z Ω f ( t ) · u d x and h G ( t ) , w i = Z Γ g ( t ) · w d x. Let us note that (1.6) involves the equation on Ω , as well as, the equation on Γ . Thus T is to beunderstood as the functional on Ω × Γ , that is trivial on Γ since no inertia is considered on the( d − Γ . In particular, the “generalized” stress Σ = W ′ Eu = ( C e ( u ) , B u | Γ )contains, besides the bulk stress tensor, also the traction stress vector. Relying on the linearityof W ′ , we have . Σ = W ′ Ev with v = . u , as used in (1.6b). Let us note that the adjoint operator3 ∗ : ( σ, ς ) F in (1.6a) with the traction force ς = B u | Γ determines a bulk force F as a functionalon test displacements u by Z Ω F · u d x = h F , u i = h E ∗ ( σ, ς ) , u i = h ( σ, ς ) , Eu i = h ( σ, ς ) , ( e ( u ) , u | Γ ) i = Z Ω σ : e ( u ) d x + Z Γ ς · u d S = Z Γ ( σ · ~n + ς ) · u d S − Z Ω div σ · u d x , which clarifies the force E ∗ Σ = F in (1.6a). The leap-frog time discretization of (1.6) then reads as Σ k +1 / τh − Σ k − / τh τ = W ′ h E h v kτh + D kτh and T ′ v k +1 τh − v kτh τ + E ∗ h Σ k +1 / τh = F k +1 / τh , (1.7)where W h and E h are suitable FEM discretizations of W and E , cf. (3.1a) and (3.1c) below, and F k +1 / τh := 1 τ Z ( k +1) τkτ F h ( t ) d t and D kτh := 1 τ Z ( k +1 / τ ( k − / τ . G h ( t ) d t = G h (( k + ) τ ) − G h (( k − ) τ ) τ . (1.8)As mentioned before, we assume here that v is discretized in an element-wise constant way sothat T leads to a diagonal matrix. In this case we do not need to employ numerical integration toapproximate the mass matrix. For higher-order discretizations, however, mass lumping is necessaryso as to obtain explicit discretization schemes. We refer to [10, 31, 50] for details in the case G ≡ W ′ h = E ∗ h W ′ h E h , which means that theresulting scheme does not require the solution of a big linear system at each iteration in time. Thespatial FEM discretization exploits regularity available in linear elastodynamics, in particular thatdiv σ and e ( v ) in (1.5a) live in L -spaces. Moreover, the equations in (1.7) are decoupled in thesense that, first, Σ k +1 / τh is calculated from the former equation and, second, v k +1 τh is calculatedfrom the latter equation assuming, that ( v kτh , Σ k − / τh ) is known from the previous time step. For k = 0, it starts from v τh = v and from a half time step Σ / τh = Σ τh + τ W ′ h E h v τh . For the spacediscretization, the lower order Q div k +1 − Q k finite element is obtained for k = 0 and in this case thevelocity is discretized as piecewise constant on rectangular or cubic elements while the stress isdiscretized by piecewise bi-linear functions with some continuities. Namely the normal componentof the stress is continuous across edges of adjacent elements while the tangential component isallowed to be discontinuous. For more details about the space discretization we refer the interestedreader to [10]. Alternative discretizations for the linear elasticity problem have been proposed byD. Arnold and his collaborators who designed mixed finite elements for general rectangular andtriangular grids [4–6]. For tetrahedral leap-frog discretization of the elastodynamics see [21]. Ingeneral, the leap-frog scheme has been frequently used in geophysics to calculate seismic wavepropagation with the finite differences method, cf. e.g. [13, 25, 51].When taking the average (i.e. the sum with the weights and ) of the second equation in(1.7) in the level k and k − v kτh and summing it with the first equation in (1.7) testedby [ W ′ h ] − ( Σ k +1 / τh + Σ k − / τh ) /
2, we obtain12 τ (cid:10) [ W ′ h ] − Σ k +1 / τh , Σ k +1 / τh (cid:11) − τ (cid:10) [ W ′ h ] − Σ k − / τh , Σ k − / τh (cid:11) = D Σ k +1 / τh + Σ k − / τh , E h v kτh E + D [ W ′ h ] − D kτh , Σ k +1 / τh + Σ k − / τh E and4 T ′ v k +1 τh − v k − τh τ , v kτh E + D Σ k +1 / τh + Σ k − / τh , E h v kτh E = D F k +1 / τh + F k − / τh , v kτh E . Summing it up, we get that the following discrete energy is conserved12 hT ′ v k +1 τh , v kτh i + Φ h ( Σ k +1 / τh ) with Φ h ( Σ ) = 12 h [ W ′ h ] − Σ, Σ i . (1.9)Note that Φ h is the discrete stored energy expressed in terms of the generalized stress. In contrastto (1.3), this formulation allows for enhancement of the discrete stored energy by some internalvariables. The energy (1.9) is shown to be a positive quantity under the following CFL condition (cid:10) [ W ′ h ] − Σ h , Σ h (cid:11) ≥ τ (cid:10) E ∗ h Σ h , ( T ′ ) − E ∗ h Σ h (cid:11) (1.10)for any Σ 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” issometimes 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 onany regularity. Then, in Section 3, we extend the two-step leap-frog discrete scheme (1.7) to asuitable three-step scheme, and study the energy properties of the proposed scheme. Then, inSection 4, we prove the numerical stability of the 3-step staggered approximation scheme andits convergence under the CFL condition (4.1). Such an abstract scheme is then illustrated inSection 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] wherethis 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 asnumerical stability and convergence) combination of the explicit staggered discretization with non-linear dissipative processes on some internal variables is new, although occasionally some dissipa-tive nonlinear phenomena can be found in literature as in [45] for a unilateral contact, in [13] fora Maxwell viscoelastic rheology, in [47] for electroactive polymers, in [23] for an aeroelastic sys-tem, or in [24] for general thermomechanical systems, but without any numerical stability (a-prioriestimates) and convergence guaranteed. The concept of internal variables has a long tradition, cf. [36], and opens wide options for materialmodelling, cf. e.g. [35, 37] and references therein. Typically, internal variables are governed by aparabolic-type 1st-order evolution. The abstract system (1.2) is thus generalized to T ′ .. u + W ′ u ( u, z ) = F ( t ) for t ∈ (0 , T ] , u | t =0 = u , . u | t =0 = v , (2.1a) ∂Ψ ( . z ) + W ′ z ( u, z ) ∋ t ∈ (0 , T ] , z | t =0 = z . (2.1b)The inclusion in (2.1b) refers to a possibility that the convex (pseudo)potential Ψ of dissipativeforces may be nonsmooth and then its subdifferential ∂Ψ can be multivalued.5ombination of the 2nd-order evolution (1.2) with such 1st-order evolution is to be handledcarefully. In contrast to the implicit staggered schemes, cf. [42], the constitutive equation is dif-ferentiated in time, cf. (1.5a), and it seems necessary to use the staggered scheme so that theinternal-variable flow rule can be used without being differentiated in time, even for a quadraticstored energy W .Moreover, to imitate the leap-frog scheme, it seems suitable (or maybe even necessary) that thestored energy W may be expressed in terms of the generalized stress as W ( u, z ) = Φ ( Σ, z ) with Σ = C Eu, and Φ ( · , z ) and Φ ( Σ, · ) quadratic , (2.2)where C stands for a “generalized” elasticity tensor and E is an abstract gradient-type operator.Typically Eu = ( e ( u ) , u | Γ ) or simply Eu = e ( u ) are considered here in the context of continuummechanics at small strains, cf. the examples in Sect. 5. Here, Σ may not directly enter the balanceof forces and is thus to be called rather as some “proto-stress”, while the actual generalized stresswill 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 . Σ = C Ev + . G ( t ) for t ∈ (0 , T ] , (2.3a) T ′ . v + E ∗ S = F ( t ) with S = C ∗ Φ ′ Σ ( Σ, z ) for t ∈ (0 , T ] , (2.3b) ∂Ψ ( . z ) + Φ ′ z ( Σ, z ) ∋ t ∈ (0 , T ] , (2.3c) Σ | t =0 = Σ := C Eu + G (0) , v | t =0 = v , z | t =0 = z . (2.3d)Here Φ ′ Σ ( Σ, z ) is a “generalized” strain and, when multiplied by C ∗ , it becomes a generalized stress.Actually, (1.6) is a special case of (2.3) when Φ = Φ ( Σ ) = R Ω C − σ : σ d x + R Γ B − ς · ς d S and S = Σ = ( σ, ς ) while C = W ′ = ( C , B ), so that indeed S = C ∗ Φ ′ Σ ( Σ ) = ( C , B )( C − σ, B − ς ) =( σ, ς ) = Σ .The energy properties of this system can be revealed by testing the particular equations/inclusionsin (2.3) by Φ ′ Σ ( Σ, z ), v , and . z . Thus, at least formally, we obtain (cid:10) Φ ′ Σ ( Σ, z ) , . Σ (cid:11) = (cid:10) Φ ′ Σ ( Σ, z ) , C Ev + . G (cid:11) = (cid:10) C ∗ Φ ′ Σ ( Σ, z ) , Ev (cid:11) + (cid:10) Φ ′ Σ ( Σ, z ) , . G (cid:11) , (2.4a) (cid:10) T ′ . v, v (cid:11) + (cid:10) C ∗ Φ ′ Σ ( Σ, z ) , Ev (cid:11) = (cid:10) F ( t ) , v (cid:11) , (2.4b) Ξ ( . z ) + (cid:10) Φ ′ z ( Σ, z ) , . z (cid:11) ≤ Ξ ( . z ) := inf (cid:10) ∂Ψ ( . z ) , . z (cid:11) . (2.4c)The functional Ξ is a dissipative rate and the “inf” in it refers to the fact that the dissipativepotential Ψ can be nonsmooth and thus the subdifferential ∂Ψ can be multivalued even at . z = 0,otherwise an equality in (2.4c) holds. Summing it up and using the calculus dd t T ( v ) = hT ′ v, . v i = hT ′ . v, v i and dd t Φ ( Σ, z ) = h Φ ′ Σ ( Σ, z ) , . Σ i + h Φ ′ z ( Σ, z ) , . z i , we obtain the following inequality for theenergy, dd t (cid:0) T ( v ) + Φ ( Σ, z ) | {z } kinetic andstored energies (cid:1) + Ξ ( . z ) | {z } dissipationrate ≤ (cid:10) F ( t ) , v (cid:11) + (cid:10) Φ ′ Σ ( Σ, z ) , . G (cid:11)| {z } power ofexternal force . (2.5)Actually, (2.4c) and (2.5) often hold as equalities.Let us now formulate some abstract functional setting of the system (2.3). For some Banachspaces S , Z , and Z ⊃ Z and for a Hilbert space H , let Φ : S × Z → R be smooth and coercive,6 : H → R be quadratic and coercive, and let Ψ : Z → [0 , + ∞ ] be convex, lower semicontinuous,and coercive on Z , cf. (4.2) below. Intentionally, we do not want to rely on any regularity whichis usually at disposal in linear problems but might be restrictive in some nonlinear problems. Forthis reason, we reconstruct the abstract “displacement” and use (2.3a) integrated in time, i.e. Σ = C Eu + G with u ( t ) := Z t v ( t ) d t + u . (2.6)Moreover, we still need another Banach space E and define the Banach space U := { u ∈ H ; Eu ∈ E} equipped with the standard graph norm. Then, by definition, we have the continuous embedding U → H and the continuous linear operator E : U → E . We assume that U is embedded into H densely, so that H ∗ ⊂ U ∗ and that H is identified with its dual H ∗ , so that we have the so-calledGelfand triple U ⊂ H ∼ = H ∗ ⊂ U ∗ . We further consider the abstract elasticity tensor C as a linear continuous operator E → S . There-fore C Eu ∈ S provided u ∈ U so that the equation (2.6) is meant in S and one needs G ( t ) ∈ S .Let us note that T ′ : H → H ∗ ∼ = H , Φ ′ Σ : S × Z → S ∗ , E ∗ : E ∗ → U ∗ , and C ∗ : S ∗ → E ∗ , so that T ′ v ∈ H ∗ provided v ∈ H and also S = C ∗ Φ ′ Σ ∈ E ∗ and E ∗ S ∈ H ∗ . In particular, the equation(2.3b) can be meant in H if integrated in time, and one needs F ( t ) valued in H .For a Banach space X , we will use the standard notation L p (0 , T ; X ) for Bochner spaces of theBochner measurable functions [0 , T ] → X whose norm is integrable with the power p or essentiallybounded if p = ∞ , and W ,p (0 , T ; X ) the space of functions from L p (0 , T ; X ) whose distributionaltime derivative is also in L p (0 , T ; X ). Also, C k (0 , T ; X ) will denote the space of functions [0 , T ] → X whose k th -derivative is continuous, and C w (0 , T ; X ) will denote the space of weakly continuousfunctions [0 , T ] → X . Later, we will also use Lin( U , E ), denoting the space of linear boundedoperators U → 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 Φ ( Σ, z ) = Φ ( Σ, z ) + Φ ( z ) with [ Φ ] ′ z : S × Z → Z ∗ and Φ ′ : Z → Z ∗ (2.7)and to use integration by-parts for the term h Φ ′ ( z ) , . z i . We thus arrive to the following definition. Definition 2.1 (Weak solution to (2.3) .) The quadruple ( u, Σ, v, z ) ∈ C w (0 , T ; U ) × C w (0 , T ; S ) × C w (0 , T ; H ) × C w (0 , T ; Z ) with Ψ ( . z ) ∈ L ( I ) and . z ∈ L (0 , T ; Z ) will be called a weak solution tothe initial-value problem (2.3) with (2.6) if v = . u in the distributional sense, Σ = C Eu + G holdsa.e. on I , and if Z T (cid:10) Φ ′ Σ ( Σ, z ) , C E e v (cid:11) S ∗ ×S − (cid:10) T ′ v, . e v (cid:11) H ∗ ×H d t = (cid:10) T ′ v , e v (0) (cid:11) H ∗ ×H + Z T (cid:10) F, e v (cid:11) H ∗ ×H d t (2.8a) for any e v ∈ C (0 , T ; H ) ∩ C (0 , T ; U ) with e v ( T ) = 0 , and Z T Ψ ( e z ) + (cid:10) [ Φ ] ′ z ( Σ, z ) , e z − . z (cid:11) Z ∗ ×Z + (cid:10) Φ ′ ( z ) , e z (cid:11) Z ∗ ×Z d t + Φ ( z ) ≥ Φ ( z ( T )) + Z T Ψ ( . z ) d t (2.8b) for any e z ∈ C (0 , T ; Z ) , where indices in the dualities h· , ·i indicate the respective spaces in dualities,and if also u (0) = u , Σ (0) = Σ , and z (0) = z . v (0) = v is contained in (2.8a). Definition 2.1works successfully for p >
1, i.e. for rate-dependent evolution of the abstract internal variable z , sothat . z ∈ L p (0 , T ; Z ). For the rate-dependent evolution when p = 1, we would need to modify it.Here, we restrict ourselves to p ≥
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 dis-cretization for (2.3c), using a fractional-step split (called also a staggered scheme) with a mid-pointformula 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 + 1instead of k + 1 / S h ⊂ S , V h ⊂ H ,and Z h ⊂ Z where the values of the respective discrete variables Σ h , v h , and z h will be, assumingthat their unions are dense in the respective Banach spaces. We will use an interpolation operator I h : Lin( S , S h ) and the embedding operator J h : Z h → Z ; it is important that the collection { J h } h> is uniformly bounded and, since S h> Z h is dense in Z , the sequence { J h } h> convergesto the identity on Z strongly. We consider E h ∈ Lin( V h , E ). Let us note that we allow for a“non-conformal” approximation of v , i.e., V h ⊂ H is not necessarily a subspace of U . This is inagreement with discretizations of the velocity as in [9–11, 18, 50].Considering that we know from previous step Σ kτh , v kτh , z kτh , then the proposed discretizationscheme is1) calculate Σ k +1 τh : Σ k +1 τh − Σ kτh τ = I h C E h v kτh + D kτh , (3.1a)2) calculate z k +1 τh : J ∗ h ∂Ψ (cid:16) z k +1 τh − z kτh τ (cid:17) + Φ ′ z (cid:16) Σ k +1 τh , z k +1 τh + z kτh (cid:17) ∋ , (3.1b)3) calculate v k +1 τh : T ′ v k +1 τh − v kτh τ + E ∗ h S k +1 τh = F k +1 / τh with S k +1 τh = C ∗ I ∗ h Φ ′ Σ ( Σ k +1 τh , z k +1 τh ) , and u k +1 τh : u k +1 τh = u kτh + τ v k +1 τh , (3.1c)where F k +1 / τh and D kτh are from (1.8). It seems important in the non-linear case to compute thevariables in the order given above. We note however, that for the linear viscoelastic problem withMaxwell rheology a scheme with a different ordering has been proposed in [27, Part I, Sect.2]. Thepotentials Φ , Ψ , and T are considered restricted on S h × Z h , Z h , and V h , so that their corresponding(sub)differentials Φ ′ Σ , Φ ′ z , ∂Ψ , and T ′ are valued in S ∗ h , Z ∗ h , and V ∗ h , respectively. The particularequations/inclusion in (3.1) are thus to be understood in S h , Z ∗ h , V ∗ h , E ∗ , and V h , respectively. Theonly implicit equation is (3.1b). Note however that even this equation becomes explicit if thereare no spatial gradients in Φ and Ψ . In view of the definition of the convex subdifferential, (3.1b)means the variational inequality Ψ (cid:0)e z (cid:1) + (cid:28) Φ ′ z (cid:16) Σ k +1 τh , z k +1 τh + z kτh (cid:17) , e z − z k +1 τh − z kτh τ (cid:29) ≥ Ψ (cid:16) z k +1 τh − z kτh τ (cid:17) (3.2)for any e z ∈ V h . In any case, the equation (3.1b) posseses a potential z τ Φ (cid:16) Σ k +1 τh , z + z kτh (cid:17) + Ψ (cid:16) z − z kτh τ (cid:17) (3.3)8hich 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 , , ..., T /τ −
1, starting from the initialconditions (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.5). More specifically,we proceed as follows: we test (3.1a) by Φ ′ Σ ( Σ k +1 τh , z k +1 τh ) + Φ ′ Σ ( Σ kτh , z kτh ), then test (3.1b) by z k +1 τh − z kτh τ , and eventually test the average of (3.1c) at the level k +1 and k by v kτh . Using that Φ ( · , z )and Φ (Σ , · ) are quadratic as assumed in (2.2), we have (cid:28) Φ ′ Σ ( Σ k +1 τh , z k +1 τh ) + Φ ′ Σ ( Σ kτh , z kτh )2 , Σ k +1 τh − Σ kτh τ (cid:29) = (cid:28) Φ ′ Σ ( Σ k +1 τh , z kτh ) + Φ ′ Σ ( Σ kτh , z kτh )2 , Σ k +1 τh − Σ kτh τ (cid:29) + τ (cid:28) Φ ′ Σ ( Σ k +1 τh , z k +1 τh ) − Φ ′ Σ ( Σ k +1 τh , z kτh ) τ , Σ k +1 τh − Σ kτh τ (cid:29) = Φ ( Σ k +1 τh , z kτh ) − Φ ( Σ kτh , z kτh ) τ + τ (cid:28) Φ ′ Σ ( Σ k +1 τh , z k +1 τh ) − Φ ′ Σ ( Σ k +1 τh , z kτh ) τ , Σ k +1 τh − Σ kτh τ (cid:29) , (3.4a)where we used also (3.1a), and (cid:28) Φ ′ z (cid:16) Σ k +1 τh , z k +1 τh + z kτh (cid:17) , z k +1 τh − z kτh τ (cid:29) = Φ ( Σ k +1 τh , z k +1 τh ) − Φ ( Σ k +1 τh , z kτh ) τ . (3.4b)Therefore, using again the particular equations/inclusion in (3.1), we get respectively Φ ( Σ k +1 τh , z kτh ) − Φ ( Σ kτh , z kτh ) τ = (cid:28) Φ ′ Σ ( Σ k +1 τh , z k +1 τh ) + Φ ′ Σ ( Σ kτh , z kτh )2 , I h C E h v kτh + D kτh (cid:29) − τ (cid:28) Φ ′ Σ ( Σ k +1 τh , z k +1 τh ) − Φ ′ Σ ( Σ k +1 τh , z kτh ) τ , Σ k +1 τh − Σ kτh τ (cid:29) , (3.5a) Ξ (cid:16) z k +1 τh − z kτh τ (cid:17) + Φ ( Σ k +1 τh , z k +1 τh ) − Φ ( Σ k +1 τh , z kτh ) τ ≤ , and (3.5b) (cid:10) T ′ v k +1 τh − v k − τh τ , v kτh (cid:11) + (cid:28) E ∗ h C ∗ I ∗ h Φ ′ Σ ( Σ k +1 τh , z k +1 τh )+ Φ ′ Σ ( Σ kτh , z kτh )2 , v kτh (cid:29) = (cid:10) F k +1 / τh , v kτh (cid:11) . (3.5c)Let us also note that, if Ψ (0) = 0 is assumed, the substitution e z = 0 into the inequality (3.2) gives Ψ ( z k +1 τh − z kτh τ ) instead of the dissipation rate Ξ ( z k +1 τh − z kτh τ ) in (3.5b), which is a suboptimal estimateexcept if Ψ is degree-1 positively homogeneous.Summing (3.5) up, we benefit from the cancellation of the terms ± Φ ( Σ k +1 τh , z kτh ), which is theusual goal and attribute of well-designed fractional-split schemes. Thus, using also the simplealgebra (cid:10) T ′ ( v k +1 τh − v k − τh ) , v kτh (cid:11) = (cid:10) T ′ v k +1 τh , v kτh (cid:11) − (cid:10) T ′ v kτh , v k − τh (cid:11) , (3.6)we obtain the analog of (2.5), namely hT ′ v k +1 τh , v kτh i − hT ′ v kτh , v k − τh i τ + Φ ( Σ k +1 τh , z k +1 τh ) − Φ ( Σ kτh , z kτh ) τ + Ξ (cid:16) z k +1 τh − z kτh τ (cid:17) (cid:10) F k +1 / τh , v kτh (cid:11) + (cid:28) Φ ′ Σ ( Σ k +1 τh , z k +1 τh ) + Φ ′ Σ ( Σ kτh , z kτh )2 , D kτh (cid:29) − τ (cid:28) Φ ′ Σ ( Σ k +1 τh , z k +1 τh ) − Φ ′ Σ ( Σ k +1 τh , z kτh ) τ , Σ k +1 τh − Σ kτh τ (cid:29) . (3.7)If Ψ is smooth except possibly at zero, there is even equality in (3.7).Considering some approximate values { z kτh } k =0 ,...,K of the variable z with K = T /τ , we definethe piecewise-constant and the piecewise affine interpolants respectively by z τh ( t ) = z kτh , z τh ( t ) = 12 z kτh + 12 z k − τh , and (3.8a) z τh ( t ) = t − ( k − ττ z kτh + kτ − tτ z k − τh for ( k − τ < t ≤ kτ. (3.8b)Similar meaning is implied for Σ τh , v τh , Σ τh , v τh , F τh , etc. The discrete scheme (3.1) can bewritten in a “compact” form as . Σ τh = I h C E h v τh + . G τh and . u τh = v τh , (3.9a) J ∗ h ∂Ψ (cid:0) . z τh (cid:1) + Φ ′ z (cid:0) Σ τh , z τh (cid:1) ∋ , (3.9b) T ′ . v τh + E ∗ h S τh = F τh with S τh = C ∗ I ∗ h Φ ′ Σ ( Σ τh , z τh ) (3.9c)to be valid a.e. on the time interval [0 , T ]. Because the energy (1.9) involves now also the internal variable, the CFL condition has to bemodified. More specifically, we assume that ∃ η > ∀ Σ h ∈ S h , z h ∈ Z h : Φ ( Σ h , z h ) ≥ τ − η (cid:10) E ∗ h S h , ( T ′ ) − E ∗ h S h (cid:11) H ∗ ×H with S h = C ∗ I ∗ h Φ ′ Σ ( Σ h , z h ) , (4.1)where Σ h and z h are considered from the corresponding finite-dimensional subspaces. Let us stillintroduce the Banach space X := { X ∈ S ∗ ; E ∗ C ∗ X ∈ H ∗ } . We further assume C ∈ Lin( E , S )invertible and that the collection of the interpolation operators { I h : S → S h } h> is bounded inLin( S , S ). Proposition 4.1 (Numerical stability.)
Let F be constant in time, valued in H ∗ , G ∈ W , (0 , T ; S ) , u ∈ U so that Σ = C Eu ∈ S , v ∈ H , z ∈ Z , the functionals T , Φ , and Ψ be coercive and Φ ′ Σ ( Σ, · ) be Lipschitz continuous uniformly for Σ ∈ S in the sense ∃ ǫ > p ≥ ∀ ( Σ, v, z ) ∈ S×H×Z : T ( v ) ≥ ǫ k v k H , Φ ( Σ, z ) ≥ ǫ k Σ k S + ǫ k z k Z , Ψ ( z ) ≥ ǫ k z k p Z , (4.2a) ∃ C ∀ Σ ∈ S , z ∈ Z : (cid:13)(cid:13) Φ ′ Σ ( Σ, z ) (cid:13)(cid:13) S ∗ ≤ C (cid:0) k Σ k S + k z k Z (cid:1) , (4.2b) ∃ ℓ ∈ R ∀ Σ ∈ S , z, e z ∈ Z : (cid:13)(cid:13) Φ ′ Σ ( Σ, z ) − Φ ′ Σ ( Σ, e z ) (cid:13)(cid:13) S ∗ ≤ ℓ k z − e z (cid:13)(cid:13) Z . (4.2c)10 et also the CFL condition (4.1) hold with τ > sufficiently small (in order to make the discreteGronwall inequality effective). Then the following a-priori estimates hold: k u τh k W , ∞ (0 ,T ; H ) ≤ C , (4.3a) k Σ τh k L ∞ (0 ,T ; S ) ≤ C and k . Σ τh k L (0 ,T ; X ∗ ) ≤ C, (4.3b) k v τh k L ∞ (0 ,T ; H ) ≤ C and kT ′ . v τh k L ∞ (0 ,T ; U ∗ ) ≤ C, (4.3c) k z τh k L ∞ (0 ,T ; Z ) ≤ C and k . z τh k L p (0 ,T ; Z ) ≤ C. (4.3d) Proof.
The energy imbalance that we have here is (3.7) which can be re-written as E k +1 h − E kh τ + Ξ (cid:16) z k +1 τh − z kτh τ (cid:17) ≤ (cid:10) F kτh , v kτh (cid:11) H ∗ ×H + (cid:28) Φ ′ Σ ( Σ k +1 τh , z k +1 τh ) + Φ ′ Σ ( Σ kτh , z kτh )2 , D kτh (cid:29) S ∗ ×S − τ (cid:28) Φ ′ Σ ( Σ k +1 τh , z k +1 τh ) − Φ ′ Σ ( Σ k +1 τh , z kτh ) τ , Σ k +1 τh − Σ kτh τ (cid:29) S ∗ ×S (4.4)with an analog of the energy (1.9), namely E k +1 h = 12 hT ′ v k +1 τh , v kτh i H ∗ ×H + Φ ( Σ k +1 τh , z k +1 τh ) . (4.5)We need to show that E k +1 h is indeed a sum of the kinetic and the stored energies at least up tosome positive coefficients. To do so, like e.g. [45, Lemma 4.2] or [50, Sect. 6.1.6], let us write hT ′ v k +1 τh , v kτh i = D T ′ v k +1 τh + v kτh , v k +1 τh + v kτh E − D T ′ v k +1 τh − v kτh , v k +1 τh − v kτh E = D T ′ v k +1 τh + v kτh , v k +1 τh + v kτh E − τ (cid:10) E ∗ h (cid:0) C ∗ I ∗ h Φ ′ Σ ( Σ k +1 τh , z k +1 τh ) − F k +1 / τh (cid:1) , ( T ′ ) − E ∗ h (cid:0) C ∗ I ∗ h Φ ′ Σ ( Σ k +1 τh , z k +1 τh ) − F k +1 / τh (cid:1)(cid:11) , (4.6)where all the duality pairings are between H ∗ and H ; here also (3.1c) has been used. Thus, usingalso T ( v ) = hT ′ v, v i , we can write the energy (4.5) as E k +1 h = T ( v k +1 / τh ) + a k +1 τh Φ ( Σ k +1 τh , z k +1 τh ) + τ (cid:10) ( T ′ ) − E ∗ h C ∗ I ∗ h Φ ′ Σ ( Σ k +1 τh , z k +1 τh ) , F k +1 / τh (cid:11) − τ k F k +1 / τh k H with a k +1 τh := 1 − τ (cid:10) E ∗ h C ∗ I ∗ h Φ ′ Σ ( Σ k +1 τh , z k +1 τh ) , ( T ′ ) − E ∗ h C ∗ I ∗ h Φ ′ Σ ( Σ k +1 τh , z k +1 τh ) (cid:11) Φ ( Σ k +1 τh , z k +1 τh ) ≥ η (4.7)and with v k +1 / τh := v k +1 τh + v kτh . The energy E k +1 h yields a-priori estimates if the coefficient a kτh isnon-negative, which is just ensured by our CFL condition (4.1) used for Σ h = Σ k +1 τh and z h = z k +1 τh .Here η > k = 0 , ..., l ≤ T /τ − ǫ (cid:18)(cid:13)(cid:13) v l +1 / τh (cid:13)(cid:13) H + a l +1 τh (cid:13)(cid:13) Σ l +1 τh (cid:13)(cid:13) S + a l +1 τh (cid:13)(cid:13) z l +1 τh (cid:13)(cid:13) Z + τ l X k =0 (cid:13)(cid:13)(cid:13) z k +1 τh − z kτh τ (cid:13)(cid:13)(cid:13) p Z (cid:19) τ k F l +1 / τh k H − τ (cid:10) ( T ′ ) − E ∗ h C ∗ I ∗ h Φ ′ Σ ( Σ l +1 τh , z l +1 τh ) , F l +1 / τh (cid:11) − τ (cid:10) ( T ′ ) − E ∗ h C ∗ I ∗ h Φ ′ Σ ( Σ τh , z τh ) , F / τh (cid:11) + T ( v − / τh ) + a τh Φ ( Σ τh , z τh ) + τ l X k =0 (cid:18)(cid:10) F k +1 / τh , v kτh (cid:11) + 12 (cid:13)(cid:13) Φ ′ Σ ( Σ k +1 τh , z k +1 τh ) + Φ ′ Σ ( Σ kτh , z kτh ) (cid:13)(cid:13) S ∗ (cid:13)(cid:13) D kτh (cid:13)(cid:13) S + τ ℓ (cid:13)(cid:13)(cid:13) z kτh + 1 − z kτh τ (cid:13)(cid:13)(cid:13) Z (cid:13)(cid:13)(cid:13) Σ k +1 τh − Σ kτh τ (cid:13)(cid:13)(cid:13) S (cid:19) , (4.8)where ǫ , p , ℓ and a l +1 τh come from (4.2) and (4.7). Here we also have used that the collection { I h } h> is bounded. Using (4.2b), we estimate k Φ ′ Σ ( Σ kτh , z kτh ) k S ∗ ≤ C (1 + k Σ kτh k S + k z kτh k Z ) and k Φ ′ Σ ( Σ k +1 τh , z k +1 τh ) k S ∗ ≤ C (1+ k Σ k +1 τh k S + k z k +1 τh k Z ), and then use the summability of k D kτh (cid:13)(cid:13) S neededfor the discrete Gronwall inequality; here the assumption . G ∈ L (0 , T ; S ) is needed. The last termin (4.8) is to be estimated by the Young inequality as τ ℓ (cid:13)(cid:13)(cid:13) z kτh + 1 − z kτh τ (cid:13)(cid:13)(cid:13) Z (cid:13)(cid:13)(cid:13) Σ k +1 τh − Σ kτh τ (cid:13)(cid:13)(cid:13) S ≤ ǫ (cid:13)(cid:13)(cid:13) z kτh + 1 − z kτh τ (cid:13)(cid:13)(cid:13) Z + ℓ ǫ (cid:13)(cid:13) Σ k +1 τh − Σ kτh (cid:13)(cid:13) S ≤ ǫ (cid:13)(cid:13)(cid:13) z kτh + 1 − z kτh τ (cid:13)(cid:13)(cid:13) p Z + C p,ǫ,ℓ (cid:0) k Σ k +1 τh k S + k Σ kτh k S (cid:1) with some C p,ǫ,ℓ depending on p , ǫ , and ℓ . Here we needed p ≥
2; note that this is related withthe specific explicit time discretization due to the last term in (3.7) but not with the problemitself. Then we use the discrete Gronwall inequality to obtain the former estimates in (4.3b,c) andthe estimates (4.3a,d). Using the discrete Gronwall inequality is a bit tricky because of the term k v l +1 / τh k H on the left-hand side of (4.8) while there is v kτh instead of v k +1 / τh on the right-hand sideof (4.8). To cope with it, we have to rely on F being constant (as assumed). We prove the estimatefor l = 1, then we sum up (4.8) for l +1 and l to get (cid:10) F k +1 / τh , v k +1 / τh (cid:11) also on the right-hand side.Note also that, in view of (3.6) for k = 0, we have obtained the term T ( v − / τh ) on the right-handside of (4.8) which, however, can simply be ignored if taking the “fictitious” velocity at level k = − − v τh = − v h .The equation . Σ τh = I h C E h v τh + . G τh gives the latter estimate in (4.3b) by estimating Z T (cid:10) . Σ τh , X (cid:11) X ∗ ×X d t = Z T (cid:10) I h C E h v τh + . G τh , X (cid:11) X ∗ ×X d t = Z T (cid:10) v τh , E ∗ h C ∗ I ∗ h X (cid:11) H×H ∗ d t + Z T (cid:10) . G τh , X (cid:11) X ∗ ×X d t (4.9)for X ∈ L ∞ (0 , T ; X ) and using also the already proved boundedness of v τh in L ∞ (0 , T ; H ) and theassumed boundedness of E h uniform in h >
0; here we used also that . Σ τh ( t ) ∈ S ⊂ X ∗ .Eventually, the already obtained estimates (4.2b) give Φ ′ Σ ( Σ τh , z τh ) bounded in L ∞ (0 , T ; S ∗ ).Therefore S τh = C ∗ I ∗ h Φ ′ Σ ( Σ τh , z τh ) is bounded in L ∞ (0 , T ; E ∗ ), hence E ∗ h S τh is bounded in L ∞ (0 , T ; U ∗ ),so that T ′ . v τh = F τh − E ∗ h S τh gives the latter estimate in (4.3c). (cid:3) We will now need the approximation properties for h → v ∈ U , v h ∈ V h , v h → v in H ⇒ E h v h → Ev in E with E ∈ Lin( U , E ) , (4.10a)Σ ∈ S , Σ h ∈ S h , Σ h → Σ in
S ⇒ I h Σ h → Σ in S . (4.10b)12 roposition 4.2 (Convergence.) Let (2.7) and (4.10a) hold, all the involved Banach spaces beseparable, and the assumptions of Proposition 4.1 hold. Moreover, let ∀ z ∈Z : Φ ′ Σ ( · , z ) continuous linear, and Φ ′ Σ : S×Z → Lin( S , S ∗ ) continuousor Φ ′ Σ : S × Z → S ∗ is continuous linear, (4.11a) ∀ z ∈Z : [ Φ ] ′ z ( · , z ) continuous linear, and [ Φ ] ′ z : S×Z → Lin( Z , Z ) continuousor [ Φ ] ′ z : S × Z → Z ∗ is continuous linear, and (4.11b) Φ ′ : Z → Z ∗ is linear continuous , (4.11c) for some Banach space Z into which Z is embedded compactly, where Φ and Φ are from (2.7) .Then there is a selected subsequence, again denoted { ( u τh , Σ τh , v τh , z τh ) } τ> converging weakly*in the topologies indicated in the estimates (4.3) to some ( u, Σ, v, z ) . Moreover, any ( u, Σ, v, z ) obtained as such a limit is a weak solution according Definition 2.1.Proof. By the Banach selection principle, we can select the weakly* converging subsequence asclaimed; here the separability of the involved Banach spaces is used.Referring to the compact embedding
Z ⊂ Z used in the former option in (4.11a,b) and relyingon a generalization the Aubin-Lions compact-embedding theorem with . z τh being bounded in thespace of the Z -valued measures on I , cf. [39, Corollary 7.9], we have z τh → z strongly in L r (0 , T ; Z )for any 1 ≤ r < + ∞ .Further, we realize that the approximate solution satisfy identities/inequality analogous to whatis used in Definition 2.1. In view of (2.8a), the equations (3.9c) now means Z T (cid:10) Φ ′ Σ ( Σ τh , z τh ) , I h C E h e v (cid:11) S ∗ ×S − (cid:10) T ′ v τh , . e v (cid:11) H ∗ ×H d t = (cid:10) T ′ v , e v (0) (cid:11) H ∗ ×H + Z T (cid:10) F h , e v (cid:11) H ∗ ×H d t (4.12a)for any e v ∈ C (0 , T ; H ) valued in V h and with e v ( T ) = 0. Like in (2.8b), the inclusion (3.9b) means Z T Ψ ( e z ) + (cid:10) [ Φ ] ′ z ( Σ τh , z τh ) , e z − . z τh (cid:11) Z ∗ ×Z + (cid:10) Φ ′ z τh , e z (cid:11) Z ∗ ×Z d t + Φ ( z ) ≥ Φ ( z τh ( T )) + Z T Ψ ( . z τh ) d t (4.12b)for all e z ∈ L (0 , T ; Z ). This is completed by (3.9a).It is further important that the equations in (3.9a) and the first equation in (3.9c) are linear,so that the weak convergence is sufficient for the limit passage there. In particular, we use (4.10a)and the Lebesgue dominated-convergence theorem.As to the weak convergence of (3.9a) integrated in time towards (3.1a) integrated in time, i.e.towards Σ = C Eu + G as used in Definition 2.1, we need to prove that Z T (cid:10) Σ τh − G τh , X (cid:11) S×S ∗ − (cid:10) u τh , E ∗ h C ∗ X (cid:11) H×H ∗ d t → Z T (cid:10) Σ − G, X (cid:11)
S×S ∗ − (cid:10) u, E ∗ C ∗ X (cid:11) H×H ∗ d t (4.13)for any X ∈ L (0 , T ; S ∗ ). By (4.10a), we have also E ∗ h S → E ∗ S in H for any S ∈ E ∗ , in particularfor S = C ∗ X ( t ). Thus certainly E ∗ h C ∗ X → E ∗ C ∗ X in L (0 , T ; H ) strongly. Using the weak*13onvergence u τh → u in L ∞ (0 , T ; H ), we obtain (4.13). Moreover, in the limit Eu = C − ( Σ − G ) ∈ L ∞ (0 , T ; E ) so that u ∈ L ∞ (0 , T ; U ).For the limit passage in (4.12a), we also use Φ ′ Σ ( Σ τh , z τh ) → Φ ′ Σ ( Σ, z ) weakly* in L ∞ (0 , T ; S ∗ )because Φ ′ Σ is continuous in the (weak × strong,weak)-mode, cf. (4.11a), and because of the men-tioned strong convergence of z τh → z .Furthermore, we need to show the convergence [ Φ ] ′ z ( Σ τh , z τh ) → [ Φ ] ′ z ( Σ, z ). For this, we useagain the mentioned generalized Aubin-Lions theorem to have the strong convergence z τh → z in L r (0 , T ; Z ) for any 1 ≤ r < + ∞ and then the continuity of [ Φ ] ′ z in the (weak × strong,weak)-mode,cf. the former option in (4.11b). The limit passage of (4.12b) towards (2.8b) then uses also theweak lower semicontinuity of Φ and the weak convergence z τh ( T ) → z ( T ) in Z ; here for thispointwise convergence in all time instants t and in particular in t = T , we also used that we havesome information about . z τh , cf. (4.3d).So far, we have relied on the former options in (4.11a,b) and the Aubin-Lions compactnessargument as far as the z -variable is concerned. If Φ is quadratic (as e.g. in the examples inSects. 5.1–5.2 below), we can use the latter options in (4.11a,b) and simplify the above arguments,relying merely on the weak convergence z τh → z and z τh → z . (cid:3) Remark 4.3 (
Alternative weak formulation)
Here, we used the weak formulation of (2.3c) con-taining the term h Φ ′ z ( Σ, z ) , . z (cid:11) which often does not have a good meaning since . z may not beregular enough in some applications. This term is thus eliminated by substituting it, after inte-gration over the time interval, by Φ ( Σ ( T ) , z ( T )) − R T h Φ ′ Σ ( Σ, z ) , . Σ (cid:11) d t − Φ ( Σ , z ) or even ratherby Φ ( Σ ( T ) , z ( T )) − R T h Φ ′ Σ ( Σ, z ) , C Ev (cid:11) d t − Φ ( Σ , z ). This however, would bring even more diffi-culties because we would need to prove a strong convergence of Φ ′ Σ ( Σ, z ), or of . Σ , or C Ev in ourexplicit-discretization scheme, which seems not easy. Remark 4.4 (
Nonquadratic Φ ( Σ, · ) ) Some applications use such Φ ( Σ, · ) which is not quadratic.This is still consistent with the explicit leap-frog-type discretization if, instead of Φ ′ z ( Σ, z ), weconsider an abstract difference quotient Φ ◦ z ( Σ, z, e z ) with the properties Φ ◦ z ( Σ, z, z ) = Φ ′ z ( Σ, z ) and (cid:10) Φ ◦ z ( Σ, z, e z ) , z − e z (cid:11) = Φ ( Σ, z ) − Φ ( Σ, e z ) , (4.14)cf. [42]. Then, instead of Φ ′ z ( Σ k +1 τh , z k +1 τh + z kτh ) in (3.1b), we should write Φ ◦ z ( Σ k +1 τh , z k +1 τh , z kτh ). Remark 4.5 (
State-dependent dissipation)
The generalization of Ψ dependent also on z or evenon ( Σ, z ) is easy. Then ∂Ψ is to be replaced by the partial subdifferential ∂ . z Ψ and (3.1b) shoulduse Ψ ( Σ k +1 τh , z kτh , · ) instead of Ψ ( · ). Remark 4.6 (
Spatial numerical approximation)
From the coercivity of the stored energy Φ , wehave Σ kτh ∈ S for any k = 0 , , ... and thus, from (3.1a), E h v kτh ∈ E so that v kτh ∈ U , althoughthe limit v cannot be assumed valued in U in general. Similarly, from (3.1c), one can read that E ∗ h S kτh ∈ H although this cannot be expected in the limit in general. Anyhow, on the time-discretelevel, one can use the FEM discretization similarly as in the linear elastodynamics where regularitycan be employed, cf. [9–11, 50] for a mixed finite-element method and [18] for the more recentlydeveloped staggered discontinuous Galerkin method for elastodynamics.14 emark 4.7 ( Other explicit-implicit schemes)
Combination of explicit and implicit time discretiza-tion might not only be due to parabolic evolution of internal variables but also due to geometricalreasons, e.g. transmission through a thin layer, that lead to a very restrictive CFL condition,cf. [14].
We present three examples from continuum mechanics of deformable bodies at small strains ofdifferent characters to illustrate applicability of the ansatz (2.2) and the above discretization scheme.Various combinations of these examples are possible, too, covering thus a relatively wide variety ofmodels.We use a standard notation concerning function spaces. Beside the Lebesgue L p -spaces, wedenote by H k ( Ω ; R n ) the Sobolev space of functions whose distributional derivatives are from L ( Ω ; R n × d k ). The simplest example with quadratic stored energy and local dissipation potential is the model ofplasticity or creep. The internal variable is then the plastic strain π , valued in the set of symmetrictrace-free matrices R d × d dev = { P ∈ R d × d ; P ⊤ = P, tr P = 0 } . For simplicity, we consider onlyhomogeneous Neumann or Dirichlet boundary conditions, so that simply E = e ( u ) and C = C . Thestored energy in terms of strain e ( u ) is W ( u, π ) = Z Ω C ( e ( u ) − π ):( e ( u ) − π ) d x , (5.1)which is actually a function of the elastic strain e el = e − π . The additive decomposition e ( u ) = e el + π is referred to as Green-Naghdi’s [26] decomposition. This energy leads to Φ ( σ, π ) = Z Ω C − σ : σ − σ : π + 12 C π : π d x with σ = C e ( u ) . (5.2)Let us note that Φ ′ σ ( σ, π ) = C − σ − π = e − π , i.e. the elastic strain e el , and that the proto-stress Σ = σ is indeed different from the actual stress σ − C π .The dissipation potential is standardly chosen as Ψ ( . π ) = Z Ω σ Y | . π | + 12 D . π : . π d x (5.3)with σ Y ≥ D a positive semidefinite viscosity tensor. The dissipationrate is then Ξ ( . π ) = R Ω σ Y | . π | + D . π : . π d x . For D > σ Y = 0, we obtain mere creep model or,in other words, the linear viscoelastic model in the Maxwell rheology . For both D > σ Y > viscoplasticity . For D = 0 and σ Y >
0, we would obtain the rate-independent (perfect)plasticity but our Proposition 4.1 does not cover this case (i.e. p = 1 is not admitted).The functional setting is H = L ( Ω ; R d ), E = S = Z = Z = L ( Ω ; R d × d sym ) where R d × d sym denotessymmetric ( d × d )-matrices. Thus U := { v ∈ L ( Ω ; R d ); e ( v ) ∈ L ( Ω ; R d × d ) } = H ( Ω ; R d ) by Korn’sinequality. 15 modification of the stored energy models an isotropic hardening , enhancing (5.1) as W ( u, π ) = Z Ω C ( e ( u ) − π ):( e ( u ) − π ) + 12 C π : π d x (5.4)so that the energy Φ from (5.2) is modified as Φ ( σ, π ) = Z Ω C − σ : σ − σ : π + 12 ( C + C ) π : π d x . (5.5)In the pure creep variant σ Y = 0, this is actually the standard linear solid (in a so-called Zenerform), considered together with the leap-frog time discretization in [7]. The isochoric constrainttr π = 0 can then be avoided, assuming that C is positive definite.All these models lead to a flow rule which is localized on each element when an element-wiseconstant approximation of π is used, and no large system of algebraic equations need to be solvedso that the combination with the explicit discretization of the other equations leads to a very fastcomputational procedure.Another modification for gradient plasticity by adding terms κ |∇ π | into the stored energy iseasily possible, too. This modification uses Z = H ( Ω ; R d × d sym ) and (2.7) with Φ ( z ) = R Ω κ |∇ π | and makes, however, the flow rule nonlocal but at least one can benefit from that the usual spacediscretization of the proto-stress σ uses the continuous piecewise smooth elements which allows forhandling gradients ∇ π if used consistently also for π . For the quasistatic variant of this model,we refer to the classical monographs [28, 49], while the dynamical model with D = 0 is e.g. in [37,Sect.5.2].Noteworthy, all these models bear time regularity if the loading is smooth and initial conditionsregular enough, which can be advantageously reflected in space FEM approximation, too.The Maxwell visco-elastodynamics was also studied by J.-P. Groby [27, Part I, Sect.2] using aslightly modified time discretization scheme, namely the order of (3.1a) and (3.1c) was exchanged.The CFL condition (4.1) here is actually the same as the standard one (1.4). This is becausethe internal variable actually does not influence the elasticity response and, likewise, the inertia isindependent of the internal variable, so the wave speed is not influenced either. The CFL is thusof the form τ ≤ h p ̺/ | λ max ( C ) | where λ max ( C ) is the maximal eigenvalue of C . Another example with quadratic stored energy but less trivial dissipation potential is a saturated
Darcy or Fick flow of a diffusant in porous media, e.g. water in porous elastic rocks or in concrete,or a solvent in elastic polymers. The most simple model is the classical
Biot model [12], capturingeffects as swelling or seepage. In a one-component flow, the internal variable is then the scalar-valued diffusant content (or concentration) denoted by ζ .As in the previous Section 5.1 , we consider only Neumann or Dirichlet boundary conditions,so that E = e ( u ). Here we use the orthogonal decomposition e = sph e + dev e with the spherical(volumetric) part sph e := (tr e ) I /d and the deviatoric part dev e and confine ourselves to isotropicmaterials where the elastic-moduli tensor C ijkl = Kδ ij δ kl + G ( δ ik δ jl + δ il δ jk − δ ij δ kl /d ) with K the bulk modulus and G the shear modulus (= the second Lam´e constant), which is the standardnotation hopefully without any confusion with the notation used in (1.6). Then the proto-stress Σ = σ = C e = K sph e + 2 G dev e . In particular, sph σ = K sph e so that tr e = K − tr σ .16dopting the gradient theory for this internal variable ζ , the stored energy in terms of strainis considered W ( u, ζ ) = Z Ω C e ( u ): e ( u ) + 12 M ( β tr e ( u ) − ζ ) + 12 L ( ζ − ζ eq ) + κ |∇ ζ | d x = Z Ω (cid:16) K + β d M (cid:17) | sph e ( u ) | + G | dev e ( u ) | − βM ζ tr e ( u ) + 12 M ζ + 12 L ( ζ − ζ eq ) + κ |∇ ζ | d x which, in terms of the (here partial) stress σ = C e , reads as R Ω ( K + β dK M ) | sph σ | + G | dev σ | − βζ MK tr σ + M ζ + L ( ζ − ζ eq ) d x . Here M > β > κ > ζ eq is a given equilibrium content.From (5.6), we arrive at the overall stored energy as: Φ ( σ, ζ ) = Z Ω (cid:16) K + β dK M (cid:17) | sph σ | + 1 G | dev σ | − βζ MK tr σ d x + Z Ω M ζ + 12 L ( ζ − ζ eq ) + κ |∇ ζ | d x , | {z } =: Φ ( ζ ) . (5.6)Let us note that Φ ′ σ ( σ, ζ ) = C − σ + βMdK ( β sph σ − ζK I ), i.e. the elastic strain, and that the proto-stress Σ = σ indeed differs from an actual stress by the spherical pressure part βMdK ( β sph σ − ζK I ).The driving force for the diffusion is the chemical potential µ = Φ ′ ζ ( σ, ζ ), i.e. here µ = ( M + L ) ζ − β MK tr σ − Lζ eq − κ ∆ ζ . (5.7a)The diffusion equation is . ζ − div( M ∇ µ ) = 0 (5.7b)with M denoting the diffusivity tensor. The system (5.7) is called the Cahn-Hilliard equation ,here combined with elasticity so that the flow of the diffusant is driven both by the gradientof concentration (Fick’s law) and the gradient of the mechanical pressure (Darcy’s law). Thedissipation potential in terms of ∇ µ , let us denote it by R behind this system, is R ( µ ) = Z Ω M ∇ µ ·∇ µ d x . (5.8)For the analysis cf. e.g. [35, Sect. 7.6].One would expect the dissipation potential as a function of the rate of internal variables, as in(2.3c). In fact, the system (5.7) turns into the form (2.3c) if one takes the dissipation potential Ψ = Ψ ( . ζ ) as Ψ ( . ζ ) = R ∗ ( . ζ ) (5.9)with R ∗ denoting the convex conjugate of R . Now, Ψ is nonlocal. The functional setting is as inthe previous example but now Z = H ( Ω ) and Z = H ( Ω ) ∗ . For a discretization of the type(3.1b), see [40]. 17ften, the diffusivity is considered dependent on ζ . Or even one can think about M = M ( σ, ζ ). Then the modification in Remark 4.5 is to be applied. In particular, R ( σ, ζ, µ ) = R Ω M ( σ, ζ ) ∇ µ ·∇ µ d x and Ψ ( σ, ζ, . ζ ) = [ R ( σ, ζ, · )] ∗ ( . ζ ).For this Biot model in the dynamical variant, the reader is also referred to the books [1,16,17,48]or also [35,37]. In any case, the diffusion involves gradients and in the implicit discretization it leadsto large systems of algebraic equations, which inevitably slows down the fast explicit discretizationof the mechanical part itself.For this case also the CFL condition (4.1) is the same as the standard one and leads to arestriction of the form τ ≤ Ch/V max where V max denotes the maximal speed with which wavespropagate in the medium. We note that the pressure velocity which is the maximal speed ofpropagation in isotropic solids is enhanced in the Biot model. The stability analysis of the discretescheme is quite technical and does not always lead to a practical CFL condition. A. Ezziani in hisPh.D thesis [22] studied a discretization of Biot’s model similar to (5.6) but the stability analysisof the discrete scheme is very nontrivial, cf. [22, Formula (7.4.11)] and, as he points out, cannot betranslated into a practical condition. Therefore, in practice he proposes to use τ ≤ a r h/V pf where V pf is the speed of the fast wave and a r is a constant depending on the order of the particularspace discretization used. The attenuation caused by diffusion causes also some dispersion of wavevelocities which stay however bounded from above by a high-frequency limit, cf. also [22, Fig. 5.2.1],so the CFL condition expectedly holds uniformly like for the pure elastodynamics. The simplest examples of nonconvex stored energy are models of damage. The most typical mod-els use as an internal variable the scalar-valued bulk damage α having the interpretation as aphenomenological volume fraction of microcracks or microvoids manifested macroscopically as acertain weakening of the elastic response. This concept was invented by L.M. Kachanov [32] andYu.N. Rabotnov [38].Considering gradient theories, the stored energy in terms of the strain and damage is hereconsidered as W ( e, α ) = Z Ω γ ( α ) C e : e + φ ( α ) + κ |∇ α | + ε ∇ ( C e ): ∇ e d x , (5.10)where φ ( · ) is an energy of damage which gives rise to an activation threshold for damage evolutionand may also lead to healing (if allowed). The last term is mainly to facilitate the mathematicstowards convergence and existence of a weak solution in such purely elastic materials withoutinvolving any viscosity, cf. [35, Sect. 7.5.3]. This regularization can also control dispersion of elasticwaves. More specifically, the 4th-order term resulted in the momentum-equilibrium equation fromthe ε -term in (5.10) causes an anomalous dispersion, i.e. waves with shorter wavelength propagatefaster than longer wavelength ones, cf. e.g. [35, Remark 6.3.6]. The ∇ α -term also facilitates theanalysis and controls the internal length-scale of damage profiles.Let us consider the “generalized” elasticity tensor C = C independent of x . As in the previousexamples, Eu = e ( u ) and G = 0. According (2.3a), the proto-stress Σ = C Eu + G , denoted by σ ,now looks as C e =: σ ; in damage mechanics, the proto-stress is also called an effective stress with18 specific mechanical interpretation, cf. [38]. In terms of σ , the stored energy is then Φ ( σ, α ) = Z Ω γ ( α ) C − σ : σ + ε ∇ C − σ : ∇ σ d x + Z Ω φ ( α ) + κ |∇ α | d x . | {z } =: Φ ( α ) (5.11)Then Φ ′ σ = γ ( α ) C − σ − div( ε ∇ ( C − σ )) and the true stress S = C ∗ Φ ′ σ is then γ ( α ) σ − div( ε ∇ σ ) pro-vided C is constant and symmetric. The damage driving force (energy) is Φ ′ α ( σ, α ) = γ ′ ( α ) C − σ : σ + φ ′ ( α ) − div( κ ∇ α ). When γ ′ (0) = 0 and φ ′ (0) ≤
0, then always α ≥ α ≥ Z = L ( Ω )with p ≥ Ψ ( . α ) = (R Ω ε . α d x + ∞ or (R Ω ε . α d x if . α ≤ Ω, R Ω . α /ε d x otherwise (5.12)with some (presumably small) coefficient ε >
0. The former option corresponds to a unidirectional(i.e. irreversible) damage not allowing any healing (as used in engineering) while the latter optionallows for (presumably slow) healing as used in geophysical models on large time scales.Since σ appears nonlinearly in Φ ′ α ( σ, α ), the strong convergence σ τh → σ in L ( Q ; R d × d ) isneeded. For this, the strain-gradient term with ε > L /ǫ (0 , T ; L d/ ( d − − ǫ ( Ω ; R d × d ))for arbitrarily small ǫ > . σ τh is bounded in some norm, which can be shown byusing . σ τh = C e ( v τh ) and the Green formula (cid:13)(cid:13) . σ τh k L ∞ (0 ,T ; H − ( Ω ; R d × d )) = sup k e e k L ,T ; H Ω ; R d × d )) ≤ Z T Z Ω . σ τh : e e d x d t = sup k e e k L ,T ; H Ω ; R d × d )) ≤ Z T Z Ω C e ( v τh ): e e d x d t = sup k e e k L ,T ; H Ω ; R d × d )) ≤ − Z T Z Ω v τh · div( C e e ) d x d t ≤ C k v τh k L ∞ (0 ,T ; L ( Ω ; R d )) with C depending on | C | . Cf. also the abstract estimation (4.9).When γ or φ are not quadratic but continuously differentiable, one can use the abstract differencequotient (4.14) defined, in the classical form, as Φ ◦ z ( Σ, α, e α ) = γ ( α ) − γ ( e α ) α − e α C − σ : σ + φ ( α ) − φ ( e α ) α − e α − κ ∆ α + e α α = e α . γ ′ ( α ) C − σ : σ + φ ′ ( α ) − κ ∆ α where α = e α . (5.13)Of course, rigorously, the ∆-operator in (5.13) is to be understood in the weak form when using itin (3.1b).Due to the gradient κ -term in (5.11), the implicit incremental problem (3.1b) leads to an al-gebraic problem with a full matrix, which may substantially slow down the otherwise fast explicitscheme. Like in the previous model the capillarity, now this gradient theory controls the length-scale of the damage profile and also serves as a regularization to facilitate mathematical analysis.19ometimes, a nonlocal “fractional” gradient can facilitate the analysis, too. Then, some waveletequivalent norm can be considered to accelerate the calculations, cf. also [3]. As far as the stress-gradient term, it is important that the discretization of the proto-stress in the usual implementationof the leap-frog method is continuous piecewise smooth, so that ∇ σ has a good sense in the dis-cretization without need to use higher-order elements. Here we use that the latter relation in(3.1c) is to be understood in the weak form, namely R Ω S k +1 τh : ˜ E h d x = h Φ ′ Σ ( Σ k +1 τh , z k +1 τh ) , C ˜ E h i for˜ E h = ˜ e h = e (˜ u h ), which means Z Ω S k +1 τh : ˜ E h d x = Z Ω γ ( α k +1 τh ) C − σ k +1 : C ˜ e h + ǫ ∇ C − σ k +1 : ∇ C ˜ e h d x for any ˜ e h from the corresponding finite-dimensional subspace of H ( Ω ; R d × d sym ). Thus we indeeddo not need higher-order elements, and also we do not need to specify explicitly homogeneousboundary conditions in this boundary-value problem.The functional setting is H = L ( Ω ; R d ), E = S = H ( Ω ; R d × d sym ), Z = H ( Ω ), and Z = Z = L ( Ω ). Then U = H ( Ω ; R d ), and E = e ( · ) is understood as an operator H ( Ω ; R d ) → H ( Ω ; R d × d sym ), and C ∗ ∼ = C ⊤ = C is understood as an a operator from H ( Ω ; R d × d sym ) to itself.A particular case of this model is a so-called phase-field fracture , taking as a basic choice γ ( α ) := ε /ε + α , φ ( α ) := g c (1 − α ) /ε, and κ := ε g c (5.14)with g c denoting the energy of fracture and with ε controlling a “characteristic” width of the phase-field fracture . The physical dimension of ε as well as of ε is m (meters) while the physical dimensionof g c is J/m . This is known as the so-called Ambrosio-Tortorelli functional [2]. In the dynamicalcontext, only various implicit discretization schemes seems to be used so far, cf. [15, 29, 44, 46].There are a lot of improvements of this basic model, approximating a mode-sensitive fracture,or ε -insensitive models (with ε referring to (5.14)), or ductile fracture, cf. [40]. This last variantcombines this model with the plasticity as in Sect. 5.1.As mentioned above in this case we have anomalous dispersion, i.e. the high frequencies prop-agate faster, cf. e.g. [35, Remark 6.3.6]. The resulting CFL condition is a combination of the usualCFL (1.4) for the 2nd-order elastodynamic model with the CFL condition for 4th-order plate asin [8]. More specifically, the speed of elastic waves in such combined model is like v ∼ v p ε/λ with v the speed in the elastodynamic case (i.e. ε = 0) and with λ the wavelength, cf. [35, Re-mark 6.3.6] for a one-dimensional analysis. For particular space discretisations, implementablewavelengths λ are bounded from below just by h . This yields to a CFL condition of the type τ ≤ C h p ε/h . (5.15)Asymptotically, for h → τ is to be small as O ( ε − / h ). For fixed ε >
0, this isactually very restrictive like in the explicit discretization of the heat equation where it practicallyprevents from efficient usage of explicit discretizations. However, here the role of ε is primarilyto facilitate rigorous existence of weak solutions of this model and can be assumed to be small.Then the influence of this 4th-order term and this restrictive asymptotics is presumably small, andthe usual CFL condition resulting from (5.15) with ε = 0 will dominate except on very fine spacediscretizations.Let us eventually remark that better asymptotics of the type τ ∼ O ( ε − / h δ ) for h → δ ) for some20 > Acknowledgments
This research has been partly supported from the grants 17-04301S (especially as far as the focuson the dissipative evolution of internal variables concerns) and 19-04956S (especially as far as thefocus on the dynamic and nonlinear behaviour concerns) of the Czech Sci. Foundation, and by theinstitutional support RVO: 61388998 ( ˇCR). T.R. also acknowledges the hospitality of FORTH inHeraklion, Crete.The authors are deeply thankful to Christos G. Panagiotopoulos for fruitful discussions. Alsothe careful and critical reading of the original version and valuable suggestions by two anonymousreferees are acknowledged.
References [1] Y.N. Abousleiman, A.H.-D. Cheng, and F.-J. Ulm, editors.
Poromechanics III: Biot Centennial(1905–2005) , London, 2005. Taylor & Francis.[2] L. Ambrosio and V.M. Tortorelli. Approximation of functional depending on jumps via byelliptic functionals via Γ-convergence.
Comm. Pure Appl. Math. , 43:999–1036, 1990.[3] M. Arndt, M. Griebel, and T.Roub´ıˇcek. Modelling and numerical simulation of martensitictransformation in shape memory alloys.
Continuum Mech. Thermodyn. , 15:463–485, 2003.[4] D.N. Arnold and G. Awanou. Rectangular mixed finite elements for elasticity.
Math. ModelsMethods Appl. Sci. , 15:1417–1429, 2005.[5] D.N. Arnold, G. Awanou, and R. Winther. Finite elements for symmetric tensors in threedimensions.
Math. Comput. , 77:1229–1251, 2008.[6] D.N. Arnold and R. Winther. Mixed finite elements for elasticity.
Numerische Mathematik ,92:401–419, 2002.[7] E. B´ecache, A. Ezziani, and P. Joly. A mixed finite element approach for viscoelastic wavepropagation.
Computational Geosciences , 8:255–299, 2004.[8] E. B´ecache, G.Derveaux, and P. Joly. An efficient numerical method for the resolution of theKirchhoff-Love dynamic plate equation.
Numer. Meth. P. D. E. , 21:323–348, 2005.[9] E. B´ecache, P. Joly, and C. Tsogka. Fictitious domains, mixed finite elements and perfectlymatched layers for 2D elastic wave propagation.
J. of Comput. Acoustics , 9(3):1175–1202,2001.[10] E. B´ecache, P. Joly, and C. Tsogka. A new family of mixed finite elements for the linearelastodynamic problem.
SIAM J. Numer. Anal. , 39:2109–2132, 2002.[11] E. B´ecache, J. Rodr´ıguez, and C. Tsogka. Convergence results of the fictitious domain methodfor a mixed formulation of the wave equation with a Neumann boundary condition.
ESAIM:Math. Model. Numer. Anal. , 43:377–398, 2009.[12] M.A. Biot. General theory of three-dimensional consolidation.
J. Appl. Phys. , 12:155–164,1941. 2113] T. Bohlen. Parallel 3-D viscoelastic finite difference seismic modelling.
Computers & Geo-sciences , 28:887–899, 2002.[14] M. Bonnet, A. Burel, M. Durufl´e, and P. Joly. Effective transmission conditions for thin-layertransmission problems in elastodynamics. The case of a planar layer model.
ESAIM: Math.Model. Numer. Anal. , 50:43–75, 2016.[15] M.J. Borden, C.V. Verhoosel, M.A. Scott, T.J.R. Hughes, and C.M. Landis. A phase-fielddescription of dynamic brittle fracture.
Comput. Meth. Appl. Mech. Engr. , 217-220:77–95,2012.[16] J.M. Carcione.
Wave Fields in Real Media, Wave Propagation in Anisotropic, Anelastic,Porous and Electromagnetic Media . Elsevier, Amsterdam, 2015.[17] A.H.-D. Cheng.
Poroelasticity . Springer, Switzerland, 2016.[18] E.T. Chung, C.Y. Lam, and J. Qian. A staggered discontinuous Galerkin method for thesimulation of seismic waves with surface topography.
Geophysics , 80:T119–T135, 2015.[19] G. Cohen and S. Pernet.
Finite Element and Discontinuous Galerkin Methods for TransientWave Equations . Springer, Dordrecht, 2017.[20] R. Courant, K. Friedrichs, and H. Lewy. ¨Uber die partiellen Differenzengleichungen der math-ematischen Physik.
Math. Annalen , 100:32–74, 1928.[21] S. Delcourte and N. Glinsky. Analysis of a high-order space and time discontinuous Galerkinmethod for elastodynamic equations. Application to 3D wave propagation.
ESAIM: Math.Model. Numer. Anal. , 49:1085–1126, 2015.[22] A. Ezziani.
Mod´elisation math´ematique et num´erique de la propagation d’ondes dans les milieuxvisco´elastiques et poro´elastiques . PhD thesis, Univ. Paris IX Dauphine, 2005.[23] C. Farhat, M. Lesoinne, and N. Maman. Mixed explicit/implicit time integration of coupledaeroelastic problems: Threefield formulation, geometric conservation and distributed solution.
Intl. J. Numer. Meth. Fluids , 21:807–835, 1995.[24] C.A. Felippa, K.C. Park, and C. Farhat. Partitioned analysis of coupled mechanical systems.
Comput. Methods Appl. Mech. Engrg. , 190:3247–3270, 2001.[25] R.W. Graves. Simulating seismic wave propagation in 3D elastic media using staggered-gridfinite differences.
Bull. Seismological Soc. Amer. , 86:1091–1106, 1996.[26] A.E. Green and P.M. Naghdi. A general theory of an elastic-plastic continuum.
Arch. RationalMech. Anal. , 18:251–281, 1965.[27] J.-P. Groby.
Mod´elisation de la propagation des ondes ´elastiques g´en´er´ees par un s´eisme procheou ´eloign´e `a l’int´erieur d’une ville . PhD thesis, Universit de la M´editerran´ee - Aix-MarseilleII, 2005.[28] W. Han and B.D. Reddy.
Plasticity . Springer, New York, 1999.[29] M. Hofacker and C. Miehe. Continuum phase field modeling of dynamic fracture: variationalprinciples and staggered FE implementation.
Intl. J. Fract. , 178:113–129, 2012.[30] M. Jir´asek. Nonlocal theories in continuum mechanics.
Acta Polytechnica , 44:16–34, 2004.[31] P. Joly and C. Tsogka.
Finite Element Methods with Discontinuous Displacement , chapter 11.Chapman & Hall/CRC, Boca Raton, FL, 2008.2232] L.M. Kachanov. Time of rupture process under creep conditions.
Izv. Akad. Nauk SSSR , 8:26,1958.[33] R. Kolman, J. Pleˇsek, J. ˇCerv, and D. Gabriel. Grid dispersion analysis of plane squarebiquadratic serendipity finite elements in transient elastodynamics.
Int. J. Numer. Meth.Engng. , 96:1–28, 2013.[34] R. Kolman, J. Pleˇsek, J. ˇCerv, M. Okrouhl´ık, and P. Paˇr´ık. Temporal-spatial dispersion andstability analysis of finite element method in explicit elastodynamics.
Intl. J. Numer. Meth.Engr. , 106:113–128, 2016.[35] M. Kruˇz´ık and T. Roub´ıˇcek.
Mathematical Methods in Continuum Mechanics of Solids .Springer, Switzeland, 2019.[36] G.A. Maugin. The saga of internal variables of state in continuum thermo-mechanics (1893-2013).
Mechanics Research Communications , 69:79–86, 2015.[37] A. Mielke and T. Roub´ıˇcek.
Rate-Independent Systems – Theory and Application . Springer,New York, 2015.[38] Yu.N. Rabotnov.
Creep Problems in Structural Members . North-Holland, Amsterdam, 1969.[39] T. Roub´ıˇcek.
Nonlinear Partial Differential Equations with Applications . Birkh¨auser, Basel,2nd edition, 2013.[40] T. Roub´ıˇcek. An energy-conserving time-discretisation scheme for poroelastic media withphase-field fracture emitting waves and heat.
Disc. Cont. Dynam. Syst. S , 10:867–893, 2017.[41] T. Roub´ıˇcek, M. Kruˇz´ık, V. Mantiˇc, C.G. Panagiotopoulos, R. Vodiˇcka, and J. Zeman. De-lamination and adhesive contacts, their mathematical modeling and numerical treatment. InV.Mantiˇc, editor,
Math. Methods and Models in Composites , chapter 11. Imperial CollegePress, 2nd edition.[42] T. Roub´ıˇcek and C.G. Panagiotopoulos. Energy-conserving time-discretisation of abstractdynamical problems with applications in continuum mechanics of solids.
Numer. Funct. Anal.Optim. , 38:1143–1172, 2017.[43] T. Roub´ıˇcek, C.G. Panagiotopoulos, and C. Tsogka. Explicit time-discretisation of elastody-namics with some inelastic processes at small strains. Preprint arXiv no.1903.11654, 2019.[44] T. Roub´ıˇcek and R. Vodiˇcka. A monolithic model for phase-field fracture and waves in solid-fluid media towards earthquakes.
Intl. J. Fracture , 219, 2019.[45] G. Scarella.
Etude th´eorique et num´erique de la propagation d’ondes en pr´esence de contactunilat´eral dans un milieu fissur´e . PhD thesis, Univ. Paris Dauphine, 2004.[46] A. Schl¨uter, A. Willenb¨ucher, C. Kuhn, and R. M¨uller. Phase field approximation of dynamicbrittle fracture.
Comput Mech , 54:1141–1161, 2014.[47] S. Seifi, K.C. Park, and H.S. Park. A staggered explicit-implicit finite element formulationforelectroactive polymers.
Comput. Methods Appl. Mech. Engrg. , 337:150–164, 2018.[48] B. Straughan.
Stability and Wave Motion in Porous Media . Springer, New York, 2008.[49] R. Temam.
Mathematical Problems in Plasticity . Gauthier-Villars, Paris, 1985. (Frenchoriginal in 1983).[50] C. Tsogka.
Modelisation math´ematique et num´erique de la propagation des ondes ´elastiquestridimensionnelles dans des milieux fissur´es . PhD thesis, Univ. Paris IX Dauphine, 1999.[51] J. Virieux. SH-wave propagation in heterogeneous media: Velocity-stress finite-differencemethod.