Yongseok Jang
email: yongseok.jang@onera.fr

Simon Shaw
email: simon.shaw@brunel.ac.uk

Discontinuous Galerkin finite element method for dynamic viscoelasticity models of power-law type A PREPRINT

Keywords: Fractional order viscoelasticity, Power-law type stress relaxation, Symmetric interior penalty Galerkin method, A priori analysis

Linear viscoelasticity can be characterized by a stress relaxation function. We consider a power-law type stress relaxation to yield a fractional order viscoelasticity model. The governing equation is a Volterra integral problem of the second kind with a weakly singular kernel. We employ spatially discontinuous Galerkin methods, symmetric interior penalty Galerkin method (SIPG) for spatial discretization, and the implicit finite difference schemes in time, Crank-Nicolson method. Further, in order to manage the weak singularity in the Volterra kernel, we use a linear interpolation technique. We present a priori stability and error analyses without relying on Grönwall's inequality, and so provide high quality bounds that do not increase exponentially in time. This indicates that our numerical scheme is well-suited for long-time simulations. Despite the limited regularity in time, we establish suboptimal fractional order accuracy in time as well as optimal convergence of SIPG. We carry out numerical experiments with varying regularity of exact solutions to validate our error estimates. Finally, we present numerical simulations based on real material data.

Viscoelasticity is a fundamental property exhibited by a wide range of materials, including polymers, gels, biological tissues, and even certain metals, e.g. see [START_REF] Hunter | Mechanics of continuous media[END_REF]. This property indicates the materials' capacity to display combined elastic and viscous behaviour. Unlike purely elastic materials, which deform instantaneously and fully recover their original shape upon removal of the load, viscoelastic materials exhibit time-dependent deformation and dissipate energy during loading and unloading. Various models have been proposed to describe the viscoelastic behaviour of materials, including the Maxwell, Kelvin-Voigt, and Zener models. These rheological models employ different combinations of springs and dashpots to represent the elastic and viscous elements of the material, providing a framework to capture the viscoelastic response. For more details, we refer to [START_REF] Shaw | Some partial differential Volterra equation problems arising in viscoelasticity[END_REF][START_REF] Findley | Creep and relaxation of nonlinear viscoelastic materials[END_REF][START_REF] Drozdov | Viscoelastic structures: mechanics of growth and aging[END_REF][START_REF] Golden | Boundary value problems in linear viscoelasticity[END_REF] and the references therein.

We begin with the momentum balance for a linear homogeneous and isotropic viscoelastic solid material (see e.g. [START_REF] Shaw | Some partial differential Volterra equation problems arising in viscoelasticity[END_REF][START_REF] Golden | Boundary value problems in linear viscoelasticity[END_REF]), given by ρü(t) -∇ • σ(t) = f (t) on Ω × (0, T ], (0.1) where Ω ⊂ R d is an open bounded polytopic domain, T > 0, u is displacement, σ is stress and f is an external body force. Here overdots denote time differentiation so that u is velocity and ü is acceleration, and it is assumed that ρ is the constant mass density of the material. In addition, we suppose a mix of essential and natural boundary conditions so that

u(t) = 0 on Γ D × [0, T ], (0.2) σ(t) • n = g N (t) on Γ N × [0, T ], (0.3)
where Γ D is the Dirichlet boundary (assumed to have positive surface measure), Γ N is the Neumann boundary given by Γ N = ∂Ω\Γ D , n is the outward unit normal vector defined a.e. on Γ N , and g N prescribes surface traction on Γ N . Furthermore, for initial conditions on the displacement and the velocity we take, u(0) = u 0 and u(0) = w 0 (0.4)

for given vector-fields u 0 and w 0 .

The constitutive relation between stress σ and strain ε characterizes the viscoelasticity model. In this paper, we focus on a power-law type constitutive model which is motivated by the intermediate concept of an elastic solid and viscous liquid in continuum mechanics such that the stress is proportional to the strain in solid, e.g. σ ∝ ε, and the stress is proportional to the rate of the strain in Newtonian fluid, e.g. σ ∝ ε. Hence the power-law type constitutive law would follow σ ∝ ∂ α t ε where ∂ α t is a fractional order time differential operator of order α with 0 < α < 1. For example, in [START_REF] Torvik | On the appearance of the fractional derivative in the behavior of real materials[END_REF], the constitutive relation in elastomer 3M-467 exhibits σ ∝ ∂ 0.56 t ε. In this setting we formulate the constitutive equation by

σ(t) = Dε(t) + 0 D α t Dε(t), (0.5)
where D and D are fourth-order tensors, and 0 D α t is a left Riemann-Liouville differential operator of order α in time (see e.g. [START_REF] Shaw | Some partial differential Volterra equation problems arising in viscoelasticity[END_REF][START_REF] Findley | Creep and relaxation of nonlinear viscoelastic materials[END_REF][START_REF] Drozdov | Viscoelastic structures: mechanics of growth and aging[END_REF]). For simplicity, we suppose D and D are piecewise constants, and defined by Dijkl = 2μδ ik δ jl + λδ ij δ kl and Dijkl = 2μδ ik δ jl + λδ ij δ kl for i, j, k, l = 1, . . . , d, where (μ, λ) and (μ, λ) are Lamé parameters, respectively. Using the notation of Cauchy's infinitesimal tensor,

ε ij (v) = 1 2
∂v i ∂x j + ∂v j ∂x i , for i, j = 1, . . . , d,
we define the strain by ε(t) = ε(u(t)) in (0.5) for convenience. For other choice of stress relaxation models, we refer to [START_REF] Shaw | Some partial differential Volterra equation problems arising in viscoelasticity[END_REF], in particular, we refer to [START_REF] Rivière | Discontinuous Galerkin finite element methods for dynamic linear solid viscoelasticity problems[END_REF][START_REF] Jang | Finite element approximation and analysis of a viscoelastic scalar wave equation with internal variable formulations[END_REF][START_REF] Jang | A priori analysis of a symmetric interior penalty discontinuous Galerkin finite element method for a dynamic linear viscoelasticity model[END_REF] for Prony series type constitutive relation with internal variables.

Solving time fractional order integro-differential equations numerically is challenging due to the non-locality and memory effects introduced by the fractional derivatives. The presence of integrals adds computational complexity, while singularities and discontinuities require special treatment. Selecting suitable numerical methods and proving stability and error estimates are additional difficulties. McLean and Thomée [START_REF] Mclean | Numerical solution of an evolution equation with a positive-type memory term[END_REF][START_REF] Mclean | Numerical solution via Laplace transforms of a fractional order evolution equation[END_REF][START_REF] Mclean | Maximum-norm error analysis of a numerical solution via Laplace transformation and quadrature of a fractional-order evolution equation[END_REF] made significant contributions to the field by developing numerical analysis techniques for a fractional order evolution equation corresponding to a scalar analogue to a power-law type fractional order viscoelasticity problem. In their work, they specifically focused on investigating the error analysis associated with the homogeneous Dirichlet boundary condition. Their research provides valuable insights and advancements in understanding the numerical aspects of fractional order evolution problems. However, since their analyses are based on spectral methods, the analyses are limited to the purely homogeneous Dirichlet boundary condition. In the works of [START_REF] Jang | Spatially continuous and discontinuous Galerkin finite element approximations for dynamic viscoelastic problems[END_REF][START_REF] Jang | A priori error analysis for a finite element approximation of dynamic viscoelasticity problems involving a fractional order integro-differential constitutive law[END_REF], the well-posedness and error estimates for the vector-valued fractional order viscoelasticity problem with a mixed boundary condition were established using duality arguments and an L ∞ approach in time, without relying on Grönwall's inequality and a spectral approach. Additionally, for a Mittag-Leffler type fractional order viscoelasticity problem, the works of [START_REF] Larsson | The continuous Galerkin method for an integro-differential equation modeling dynamic fractional order viscoelasticity[END_REF][START_REF] Larsson | Discontinuous Galerkin method for an integro-differential equation modeling dynamic fractional order viscoelasticity[END_REF][START_REF] Saedpanah | Existence and convergence of Galerkin approximation for second order hyperbolic equations with memory term[END_REF] provide relevant contributions to the analysis of such problems.

In this paper, we approximate the dynamic fractional order viscoelasticity model of a power-law type with discontinuous Galerkin finite element method (DGFEM), specifically symmetric interior penalty Garlerkin method (SIPG), for spatial discretization and Crank-Nicolson type finite difference method for temporal discretization. Due to the presence of the singularity in the fractional order Volterra kernel, particularly in the presence of non-smooth initial or boundary conditions, they require special treatment to ensure accuracy and stability. Classical numerical schemes such as standard quadrature rules developed for integer-order integral problems may not be directly applicable or may lose accuracy when applied to fractional order equations. To address this, we incorporate the linear interpolation technique [START_REF] Li | Numerical approaches to fractional calculus and fractional ordinary differential equation[END_REF][START_REF] Linz | Theoretical numerical analysis: an introduction to advanced techniques[END_REF] to handle the weak singularity and to improve accuracy. Stability bounds and spatially optimal error bounds for discrete problems are demonstrated without relying on Grönwall's inequality to avoid exponential growth in time of the so-called generic constants. Furthermore, the regularity of solutions is analyzed to address weak singularities and derive suboptimal and optimal orders of convergence with respect to time.

We would like to highlight that the stability estimates of the fractional order integro-differential equation with a mixed boundary condition can be demonstrated without the use of Grönwall's inequality. Instead, we employ the positivity property in fractional order integration and Markov's inequality to prove stability bounds for semi-discrete and fully discrete problems, respectively. Despite the presence of a weak singularity in the power-law type model and the limited regularity of solutions in time, the fully discrete solutions achieve a higher order of accuracy compared to first-order schemes. This enhanced accuracy is verified through duality arguments and an L ∞ approach in time without Grönwall's inequality and spectral methods. To our knowledge, our study presents, for the first time, stability and a priori error analyses of SIPG for the dynamic viscoelasticity model of power-law type with a mixed boundary condition including the purely elastic response. We can only find certain research works with further assumptions such as imposing Γ D = ∂Ω, vanishing D or problems of Mittag-Leffler type, where the Mittag-Leffler type kernel is asymptotically equivalent to the power-law type as t → 0. The actual computational costs in Mittag-Leffler type simulations are more expensive than the power-law type since the Mittag-Leffler type kernel involves an infinite series. Moreover, the numerical scheme in [START_REF] Larsson | Discontinuous Galerkin method for an integro-differential equation modeling dynamic fractional order viscoelasticity[END_REF] exhibits only first order accuracy in time and [START_REF] Saedpanah | Existence and convergence of Galerkin approximation for second order hyperbolic equations with memory term[END_REF] shows optimal spatial error estimates using Grönwall's inequality without temporal error analysis. Therefore, the novelty of our work is the improved analyses of stability and a priori bounds for more generalized dynamic fractional order viscoelasticity problems where the bounds are non-exponentially increasing in time to give confidence in the long time simulation of viscoelastic response. This article is structured as follows: In Section 1, we introduce the fundamental definitions of fractional calculus, the frameworks of DGFEM, and our notation. Section 2 defines a semi-discrete formulation along with its stability analysis, as well as a fully discrete formulation. The stability analysis and a priori error bounds for the fully discrete problem are stated and proved in Section 3. Numerical experiments using FEniCS (https://fenicsproject.org/) are presented in Section 4. Finally, Section 5 concludes the article.

Preliminary

We use standard notation so that L p (Ω), H s (Ω) and W s p (Ω) (with s and p non-negative) denote the usual Lebesgue, Hilbert and Sobolev spaces. For any normed space X, ∥•∥ X represents the X norm which, for inner product spaces, is always the norm induced by the inner product. For example, ∥•∥ L2(Ω) is the L 2 (Ω) norm, as induced by the L 2 (Ω) inner product denoted-for brevity-by (•, •), but for S ⊂ Ω, we use (•, •) L2(S) for the L 2 (S) inner product. When we denote the Bochner space by L p (0, T ; X), for a time-dependent function f ∈ L p (0, T ; X), the corresponding norm is defined by

∥f ∥ Lp(0,T ;X) = T 0 ∥f (t)∥ p X dt 1/p , for 1 ≤ p < ∞. When p = ∞ this becomes the essential supremum norm: ∥f ∥ L∞(0,T ;X) = ess sup 0≤t≤T ∥f (t)∥ X .
When convenient, we shall often replace the upper limit T in these expressions by some other value t ∈ [0, T ].

For inner products of vector-valued and tensor-valued functions we use the same notation as for the scalar cases. For instance, we have

(v, w) = Ω v • w dΩ, (v, w) = Ω v : w dΩ = d i,j=1 Ω v ij w ij dΩ,
for vector-valued functions v and w, and second order tensors v and w.

We follow the framework of the DGFEM in [START_REF] Rivière | Discontinuous Galerkin methods for solving elliptic and parabolic equations: theory and implementation[END_REF] and refer to it for a detailed explanation. Assume that the closure of Ω is subdivided by E h , where E ∈ E h is a triangle in 2D or a tetrahedron in 3D, and the intersection of any pair of elements is either a vertex, an edge, a face, or empty. We suppose that the subdivision is quasi-uniform, which means that there exists a positive constant C such that h ≤ Ch E for any E where h E is the diameter of E ∈ E h and h is the maximum diameter. Let Γ h be the set of interior edges (in 2D) or faces (in 3D) contained in the subdivision E h . Then for each edge or face element e, we can define a unit normal vector, n e . If e ⊂ ∂Ω, n e is the outward unit normal vector. For an interior edge e such that e ⊂ E i ∩ E j with i < j, the normal vector n e is oriented from

E i to E j .
With the subdivision, we can introduce the broken Sobolev space

H s (E h ) = v ∈ L 2 (Ω) | ∀E ∈ E h , v| E ∈ H s (E)
and endow it with the broken Sobolev norm, |||•||| H s (E h ) , defined by

|||v||| H s (E h ) = E∈E h ∥v∥ 2 H s (E) 1/2 .
These definitions and notations are extended in an obvious way to the the vector field analogue H s (E h ). We can also define the space of polynomials of degree less than or equal to k over E such that

P k (E) = span x i1 1 • • • x i d d | d m=1 i m ≤ k, x ∈ E, i m ∈ N ∪ {0} for each m ,
and then define our DG finite element space as

D k (E h ) = v ∈ H 1 (E h ) v| E ∈ P k (E) for each E ∈ E h .
The analogous vector field is given by D

k (E h ) := [D k (E h )] d .
Next, we want to define an average and a jump for a vector valued function v and a second order tensor v between two elements E e i and E e j sharing the common edge e with i < j by We can now introduce the jump penalty operator,

{v} = (v| E e i )| e + (v| E e j )| e 2 , {v} = (v| E e i )| e + (v| E e j )| e 2 , [v] = (v| E e i )| e -(v| E e j )| e , [v ⊗ n e ] = (v| E e i )|
J γ0,γ1 0 (v, w) = e⊂Γ h ∪Γ D γ 0 |e| γ1 e [v] • [w] de,
where γ 0 and γ 1 are positive constants.

Useful inequalities

We now provide the following inequalities for use later in the a priori analysis.

• Inverse polynomial trace inequalities [START_REF] Warburton | On the constants in hp-finite element trace inverse inequalities[END_REF]:

For any v ∈ P k (E), ∀e ⊂ ∂E, v L2(e) ≤ Ch -1/2 E v L2(E) , and ∇v • n e L2(e) ≤ Ch -1/2 E ∇v L2(E) , (1.1)
where C is a positive constant and is independent of h E but depends on the polynomial degree k. • Poincaré's Inequality [START_REF] Brenner | Poincaré-Friedrichs inequalities for piecewise H 1 functions[END_REF][START_REF] Rivière | Discontinuous Galerkin methods for solving elliptic and parabolic equations: theory and implementation[END_REF]

: If γ 1 (d -1) ≥ 1 and |e| ≤ 1 for every e ⊂ Γ h ∪ Γ D , then, v L2(Ω) ≤ C |||∇v||| 2 H 0 (E h ) + e⊂Γ h ∪Γ D 1 |e| γ1 [v] 2 L2(e) 1/2 , (1.2)
for any v ∈ H 1 (E h ). • Inverse Inequality (or Markov Inequality) [START_REF] Ozisik | On the constants in inverse inequalities in L 2[END_REF][START_REF] Rivière | Discontinuous Galerkin methods for solving elliptic and parabolic equations: theory and implementation[END_REF]: For any E ∈ E h , there is a positive constant C such that

∀v ∈ P k (E), ∇ j v L2(E) ≤ Ch -j E v L2(E) , ∀j ∈ {0, 1, . . . , k}, (1.3)
where

∇ j v = ∇ • ∇ j-1 v for even j, ∇(∇ j-1 v) for odd j, and
∇ 0 v = v.
Note that (1.1)-(1.3) can also be applied to vector-valued functions, componentwisely.

Next, we present the definition of the (left) Riemann-Liouville fractional derivative as well as its properties. Definition (Riemann-Liouville fractional derivative and integral). Let f be a function defined on [a, b] and α ∈ (0, 1). A left Riemann-Liouville derivative of order α and a left fractional integral of order α are defined by for t > a

a D α t f (t) = 1 Γ(1 -α) d dt t a f (t ′ )(t -t ′ ) -α dt ′ and a I α t f (t) = 1 Γ(α) t a f (t ′ )(t -t ′ ) α-1 dt ′
where Γ is the gamma function. We can observe that

a D α t f (t) = d dt a I 1-α t f (t), and a D α t f (t) = f (a)(t -a) -α Γ(1 -α) + a I 1-α t ḟ (t) if f is differentiable.
Furthermore, we have the positive definiteness [START_REF] Mclean | Numerical solution of an evolution equation with a positive-type memory term[END_REF] of the fractional integral of order α ∈ (0, 1) such that

T 0 0 I 1-α t ϕ(t)ϕ(t)dt = 1 Γ(1 -α) T 0 t 0 (t -t ′ ) -α ϕ(t ′ )ϕ(t)dt ′ dt ≥ 0.
(1.4)

Space and time discretization

In this section, we first introduce the SIPG formulation for the momentum equation (0.1) to derive a semi-discrete problem.

Then, using a Crank-Nicolson type finite difference scheme, we can formulate a fully discrete problem. For this, we additionally consider a numerical scheme for the fractional order derivative using the linear interpolation technique in [START_REF] Li | Numerical approaches to fractional calculus and fractional ordinary differential equation[END_REF].

With the use of the power-law type stress relaxation in [START_REF] Shaw | Numerical solution of linear quasistatic hereditary viscoelasticity problems ii: a posteriori estimates[END_REF][START_REF] Jang | A priori error analysis for a finite element approximation of dynamic viscoelasticity problems involving a fractional order integro-differential constitutive law[END_REF], we can derive the power-law type constitutive equation (0.5) as

σ(t) = φ 0 Dε(t) + φ 1 t -α Dε(0) + φ 1 Γ(1 -α) 0 I 1-α t D ε(t), (2.1)
where φ 0 is non-negative, φ 1 is positive and D is a symmetric positive definite piecewise constant fourth order tensor. By substitution of the constitutive law (2.1) into the momentum equation (0.1), the model problem becomes:

ρü(t) -∇ • φ 0 Dε(t) + φ 1 Γ(1 -α) 0 I 1-α t D ε(t) = f (t) + t -α ψ 0 on Ω × (0, T ], (2.2)
where ψ 0 := φ 1 ∇ • Dε(0), with the boundary conditions (0.2) and (0.3), and the initial condition (0.4).

If we set φ 1 = 0 in the constitutive equation, the stress relaxation will disappear. Hence, the corresponding constitutive relation implies linear elasticity. For the purpose of our study in viscoelasticity, we shall assume a positive value for φ 1 . On the other hand, in the absence of elastic response by φ 0 = 0, the model problem exhibits only stress relaxation. Furthermore, the momentum equation can be simplified to a parabolic type integro-differential equation, not a hyperbolic type. For more details on the reduced problem, refer to [START_REF] Jang | A priori error analysis for a finite element approximation of dynamic viscoelasticity problems involving a fractional order integro-differential constitutive law[END_REF].

A semi-discrete problem

To formulate a spatially discrete approximation of (2.2), we first present the SIPG form in the context of our model problem. We define a symmetric DG bilinear form a :

H s (E h ) × H s (E h ) → R for s > 3/2 by a (v, w) = E∈E h E Dε(v) : ε(w) dE - e⊂Γ h ∪Γ D e {Dε(v)} : [w ⊗ n e ] de - e⊂Γ h ∪Γ D e {Dε(w)} : [v ⊗ n e ] de + J γ0,γ1 0 (v, w), (2.3)
for any v, w ∈ H s (E h ). We also define our DG energy norm by

v V = E∈E h E Dε(v) : ε(v) dE + J γ0,γ1 0 (v, v) 1/2 , for v ∈ H s (E h ).
Comparing these we can observe that

a (v, v) = v 2 V -2 e⊂Γ h ∪Γ D e {Dε(v)} : [v ⊗ n e ] de. (2.4)
In the DG bilinear form, the third term is called the "interior penalty" term, while the last term is referred to as the "jump penalty". Depending on the sign of the interior penalty, the bilinear form is either symmetric or nonsymmetric. In this article, we consider only the symmetric DG method and refer to [START_REF] Jang | Spatially continuous and discontinuous Galerkin finite element approximations for dynamic viscoelastic problems[END_REF][START_REF] Rivière | Discontinuous Galerkin finite element methods for dynamic linear solid viscoelasticity problems[END_REF] for an application of the nonsymmetric method for viscoelasticity. The choice to employ SIPG is motivated by the fact that it only requires the standard penalization parameter γ 1 (d -1) ≥ 1 to achieve optimal spatial error estimates. In contrast, the nonsymmetric interior penalty Galerkin method (NIPG) demands a super penalization parameter γ 1 (d -1) ≥ 3 for optimal error estimates. It is worth noting that using super penalization can result in a more ill-conditioned linear system, potentially leading to difficulties when solving the system with iterative solvers. For more comprehensive details, we refer to [START_REF] Jang | Spatially continuous and discontinuous Galerkin finite element approximations for dynamic viscoelastic problems[END_REF]. Remark (Korn's inequality for piecewise H 1 vector fields [START_REF] Brenner | Korn's inequalities for piecewise H 1 vector fields[END_REF][START_REF] Rivière | Discontinuous Galerkin methods for solving elliptic and parabolic equations: theory and implementation[END_REF]). If we have γ 1 (d -1) ≥ 1, then since D is symmetric positive definite and the jump penalty is defined not only on the interior edges but also on the positive measured Dirichlet boundary Γ D , Korn's inequality yields, for any v ∈ H 1 (E h ),

E∈E h ∇v 2 L2(E) ≤ C v 2 V , (2.5)
for some positive C independent of v.

Proposition 1 (DG elliptic projection). The DG elliptic projector, R, is defined for u ∈ H s (E h ) and s > 3/2 by, R :

H s (E h ) → D k (E h ) such that a (u, v) = a (Ru, v) , ∀v ∈ D k (E h
). Note that we have the Galerkin orthogonality such that a (u -Ru, v) = 0 for any v ∈ D k (E h ). Referring to [START_REF] Rivière | Discontinuous Galerkin methods for solving elliptic and parabolic equations: theory and implementation[END_REF][START_REF] Riviére | Discontinuous Galerkin finite element methods for linear elasticity and quasistatic linear viscoelasticity[END_REF][START_REF] Houston | An hp-adaptive mixed discontinuous Galerkin FEM for nearly incompressible linear elasticity[END_REF][START_REF] Wihler | Locking-free adaptive discontinuous Galerkin FEM for linear elasticity problems[END_REF], for example, we recall the following elliptic-error estimates,

u -Ru V ≤ Ch min(k+1,s)-1 |||u||| H s (E h ) and u -Ru L2(Ω) ≤ Ch min(k+1,s) |||u||| H s (E h ) , (2.6)
for u ∈ H s (E h ) with s > 3/2 and for sufficiently large penalty parameters γ 0 > 0 and γ 1 ≥ (d -1) -1 . Here, the positive constant C is independent of u but dependent on the domain, its boundary, and the polynomial degree k. Proposition 2 (Bounds for interior penalty term [START_REF] Jang | A priori analysis of a symmetric interior penalty discontinuous Galerkin finite element method for a dynamic linear viscoelasticity model[END_REF][START_REF] Jang | Spatially continuous and discontinuous Galerkin finite element approximations for dynamic viscoelastic problems[END_REF] ). Suppose γ 0 > 0 and γ 1 (d -1) ≥ 1. For any v, w ∈ D k (E h ), we have

e⊂Γ h ∪Γ D e {Dε(v)} : [w ⊗ n e ] de ≤ C √ γ 0 v 2 V + J γ0,γ1 0 (w, w) , (2.7)
where C is a positive constant independent of v and w but dependent on the inverse polynomial trace inequality's constants and the domain. Proposition 3 (Coercivity and continuity [START_REF] Jang | A priori analysis of a symmetric interior penalty discontinuous Galerkin finite element method for a dynamic linear viscoelasticity model[END_REF][START_REF] Jang | Spatially continuous and discontinuous Galerkin finite element approximations for dynamic viscoelastic problems[END_REF]). Suppose γ 0 > 0 is sufficiently large and γ 1 (d -1) ≥ 1. Then there exist positive constants κ and K such that

κ v 2 V ≤ a (v, v) , and |a (v, w)| ≤ K v V w V , ∀v, w ∈ D k (E h ), (2.8
) where κ and K are independent of v and w.

Hereafter, we assume s ≥ 2, γ 0 > 0 and γ 1 (d -1) ≥ 1 in our paper to fulfill Proposition 1-3. Remark. On account of using the inverse polynomial trace inequality to prove the coercivity and continuity, the DG bilinear form will not be coercive and continuous on the broken Sobolev space. In other words, (2.8) holds only on the finite element space. For the choice of the penalty parameter γ 0 , we refer to [START_REF] Hansbo | Discontinuous Galerkin methods for incompressible and nearly incompressible elasticity by Nitsche's method[END_REF][START_REF] Wihler | Locking-free adaptive discontinuous Galerkin FEM for linear elasticity problems[END_REF]. For instance, we will take γ 0 ∈ [10, 100] in the numerical experiments Section 4.

In the usual way for DGFEM, we follow the standard method of multiplying (2.2) by a test function in H s (E h ) for each E ∈ E h , integrating by parts over the element, summing over all E and then imposing additional penalty terms. This produces a weak formulation of (2.2). Under the assumption that the strong solution to the momentum equation has spatial continuity, the interior penalty and jump penalty terms can be included in our DG formulation. We refer to [START_REF] Jang | Spatially continuous and discontinuous Galerkin finite element approximations for dynamic viscoelastic problems[END_REF] for more details of the DG formulation for the viscoelasticity model problems. Therefore, we can obtain the following semi-discrete problem: find

u h : [0, T ] → D k (E h ) satisfying for any v ∈ D k (E h ), ρ (ü h (t), v) + φ 0 a (u h (t), v) + φ α a 0 I 1-α t uh (t), v = F (t; v) for t > 0, (2.9) a (u h (0), v) = a (u 0 , v) , (2.10) a ( uh (0), v) = a (w 0 , v) , (2.11
) where φ α := φ 1 Γ(1 -α) for convenience and the linear form F (•) is defined by

F (t; v) = (f (t), v) + (g N (t), v) L2(Γ N ) + t -α ψ 0 , v .
Note that the bilinear form is only associated with spatial variables so that the Leibniz integral rule leads to

a 0 I 1-α t uh (t), v = 0 I 1-α t a ( uh (t), v) and t -α a (u 0 , v) = a t -α u 0 , v , for any v ∈ D k (E h ).
The above arguments are used to obtain (2.9). It is easy to show the linear form

F (•) is continuous if f ∈ C(0, T ; L 2 (Ω)), g N ∈ C(0, T ; L 2 (Γ N )), and
ψ 0 ∈ L 2 (Ω). Indeed, we can observe that if u 0 ∈ H 2 (Ω), it implies ψ 0 ∈ L 2 (Ω).
Hereafter, we assume the data terms are bounded and smooth enough to satisfy the continuity condition of the linear form. In addition to the initial data, we suppose

u 0 ∈ H 2 (Ω) ∩ C(Ω) and w 0 ∈ H 2 (Ω) ∩ C(Ω).
For the existence and uniqueness of the solution, we refer to [START_REF] Saedpanah | Existence and convergence of Galerkin approximation for second order hyperbolic equations with memory term[END_REF].

Theorem 1. Let u h be a solution of the semi-discrete problem (2.9)-(2.11). Suppose uh ∈ L ∞ (0, T ; L 2 (Ω)) and u h ∈ L ∞ (0, T ; H s (E h )). If γ 0 is large enough and γ 1 (d -1) ≥ 1, there exists a positive constant C such that ∥ uh ∥ 2 L∞(0,T ;L2(Ω)) + ∥u h ∥ 2 L∞(0,T ;V ) ≤ C ∥w 0 ∥ 2 H 2 (Ω) + ∥u 0 ∥ 2 H 2 (Ω) + T ∥f ∥ 2 L2(0,T ;L2(Ω)) + T h -1 ∥g N ∥ 2 L2(0,T ;L2(Γ N )) + T 2(1-α) ψ 0 2 L2(Ω) .
Here, the constant C is independent of the semi-discrete solution, T and h but dependent on the polynomial degree k, the fractional order α, the domain Ω, its boundary, and the material properties.

Proof. Let v = uh (t) for t ∈ (0, T ] in (2.9) to get,

ρ 2 d dt uh (t) 2 L2(Ω) + φ 0 2 d dt u h (t) 2 V -φ 0 e⊂Γ h ∪Γ D e {Dε(u h (t))} : [ uh (t) ⊗ n e ] de -φ 0 e⊂Γ h ∪Γ D e {Dε( uh (t))} : [u h (t) ⊗ n e ] de + φ α a 0 I 1-α t uh (t), uh (t) = F ( uh (t)).
(2.12)

For 0 < τ ≤ T , time integration yields

ρ 2 uh (τ ) 2 L2(Ω) + φ 0 2 u h (τ ) 2 V + φ α τ 0 a 0 I 1-α t uh (t), uh (t) dt = τ 0 F ( uh (t))dt + ρ 2 uh (0) 2 L2(Ω) + φ 0 2 u h (0) 2 V + φ 0 e⊂Γ h ∪Γ D e ({Dε(u h (τ ))} : [u h (τ ) ⊗ n e ] -{Dε(u h (0))} : [u h (0) ⊗ n e ]) de, (2.13)
by the definition of the bilinear form and integration by parts in time. Since we have

τ 0 a 0 I 1-α t uh (t), uh (t) dt = 1 Γ(1 -α) τ 0 t 0 (t -s) -α a ( uh (t ′ ), uh (t)) dt ′ dt,
and the bilinear form a (•, •) is symmetric positive definite, (1.4) implies τ 0 a 0 I 1-α t uh (t), uh (t) dt ≥ 0.
Hence, we can get

ρ 2 uh (τ ) 2 L2(Ω) + φ 0 2 u h (τ ) 2 V ≤ τ 0 F ( uh (t)) dt + ρ 2 uh (0) 2 L2(Ω) + φ 0 2 u h (0) 2 V + φ 0 e⊂Γ h ∪Γ D e ({Dε(u h (τ ))} : [u h (τ ) ⊗ n e ] -{Dε(u h (0))} : [u h (0) ⊗ n e ]) de.
(2.14)

Next, we shall show that the right-hand side of (2.14) is bounded. Consider the first term of the right-hand side. By the definition of the linear form,

τ 0 F ( uh (t))dt = τ 0 (f (t), uh (t)) dt + τ 0 (g N (t), uh (t)) L2(Γ N ) dt + τ 0 t -α ψ 0 , uh (t) dt.
Using Cauchy-Schwarz inequality, we have

τ 0 F ( uh (t))dt ≤ τ 0 f (t) L2(Ω) uh (t) L2(Ω) dt + τ 0 g N (t) L2(Γ N ) uh (t) L2(Γ N ) dt + τ 0 t -α ψ 0 L2(Ω) uh (t) L2(Ω) dt ≤ τ 0 f (t) L2(Ω) uh (t) L2(Ω) dt + Ch -1/2 τ 0 g N (t) L2(Γ N ) uh (t) L2(Ω) dt + τ 0 t -α ψ 0 L2(Ω) uh (t) L2(Ω) dt,
with the inverse polynomial trace inequality (1.1) on the second term. An L ∞ norm argument over time for uh (t) L2(Ω) implies that

τ 0 F ( uh (t))dt ≤ ∥ uh ∥ L∞(0,T ;L2(Ω)) τ 0 f (t) L2(Ω) + Ch -1/2 g N (t) L2(Γ N ) + t -α ψ 0 L2(Ω) dt .
Thus, by employing Young's inequalities and Cauchy-Schwarz inequalities, we can obtain the bound of the linear form such that

τ 0 F ( uh (t))dt ≤ ϵ a + ϵ b + ϵ c 2 ∥ uh ∥ 2 L∞(0,T ;L2(Ω)) + T 2ϵ a ∥f ∥ 2 L2(0,T ;L2(Ω))
(2.15)

+ CT h -1 2ϵ b ∥g N ∥ 2 L2(0,T ;L2(Γ N )) + T 2(1-α) 2(1 -α)ϵ c ψ 0 2 L2(Ω) , (2.16)
for any positive constants ϵ a , ϵ b and ϵ c . We refer to [START_REF] Riviere | Discontinuous Galerkin finite element methods for dynamic linear solid viscoelasticity problems[END_REF] for the bound

u h (0) 2 V ≤ C ∥u 0 ∥ 2 H 2 (Ω) , (2.17)
and we derive from (2.11) by Korn's inequality (2.5) and Poincaré's inequality (1.2)

uh (0) 2 L2(Ω) ≤ C uh (0) 2 V ≤ C ∥w 0 ∥ 2 H 2 (Ω) . (2.18)
On the other hand, (2.7) leads to

e⊂Γ h ∪Γ D e ({Dε(u h (τ ))} : [u h (τ ) ⊗ n e ] -{Dε(u h (0))} : [u h (0) ⊗ n e ]) ≤ C √ γ 0 u h (τ ) 2 V + u h (0) 2 V ≤ C √ γ 0 ∥u h ∥ 2 L∞(0,T ;V ) + ∥u 0 ∥ 2 H 2 (Ω) , (2.19
2 L2(Ω) + φ 0 2 u h (τ ) 2 V ≤ ϵ a + ϵ b + ϵ c 2 ∥ uh ∥ 2 L∞(0,T ;L2(Ω)) + C √ γ 0 ∥u h ∥ 2 L∞(0,T ;V ) + C ρ 2 ∥w 0 ∥ 2 H 2 (Ω) + C √ γ 0 ∥u 0 ∥ H 2 (Ω) + T 2ϵ a ∥f ∥ 2 L2(0,T ;L2(Ω)) + CT h -1 2ϵ b ∥g N ∥ 2 L2(0,T ;L2(Γ N )) + T 2(1-α) 2(1 -α)ϵ c ψ 0 2 L2(Ω) . (2.20)
Since τ is arbitrary, we can write (2.20) as

ρ 2 ∥ uh ∥ 2 L∞(0,T ;L2(Ω)) + φ 0 2 ∥u h ∥ 2 L∞(0,T ;V ) ≤(ϵ a + ϵ b + ϵ c ) ∥ uh ∥ 2 L∞(0,T ;L2(Ω)) + C √ γ 0 ∥u h ∥ 2 L∞(0,T ;V ) + Cρ ∥w 0 ∥ 2 H 2 (Ω) + C √ γ 0 ∥u 0 ∥ H 2 (Ω) + T ϵ a ∥f ∥ 2 L2(0,T ;L2(Ω)) + CT h -1 ϵ b ∥g N ∥ 2 L2(0,T ;L2(Γ N )) + T 2(1-α) (1 -α)ϵ c ψ 0 2 L2(Ω) ,
and setting

ϵ a = ϵ b = ϵ c = ρ/12 yields ρ 4 ∥ uh ∥ 2 L∞(0,T ;L2(Ω)) + φ 0 2 - C √ γ 0 ∥u h ∥ 2 L∞(0,T ;V ) ≤Cρ ∥w 0 ∥ 2 H 2 (Ω) + C √ γ 0 ∥u 0 ∥ H 2 (Ω) + 12T ρ ∥f ∥ 2 L2(0,T ;L2(Ω)) + CT h -1 ρ ∥g N ∥ 2 L2(0,T ;L2(Γ N )) + 12T 2(1-α) (1 -α)ρ ψ 0 2 L2(Ω) . (2.21)
Requiring a sufficiently large γ 0 to give φ 0 /2 -C/ √ γ 0 > 0, we complete the proof.

If we used Grönwall's inequality to estimate for the time integration, the stability bound would depend on time exponentially, for example, the Grönwall's inequality leads to

τ 0 F ( uh (t))dt ≤C exp (T ) ∥ uh ∥ L2(0,T ;L2(Ω)) τ 0 f (t) L2(Ω) + h -1/2 g N (t) L2(Γ N ) + t -α ψ 0 L2(Ω) dt .
Whereas our stability estimates present non-exponential bounds in time. Remark. Even though the h -1 term appears on the traction in the stability bound, it has no significant effect since h is fixed for the finite element space. While the h -1 term appears from the inverse polynomial trace inequality, the numerical results will only show the weakly imposed Neumann boundary condition, e.g. see [START_REF] Rivière | Discontinuous Galerkin methods for solving elliptic and parabolic equations: theory and implementation[END_REF][START_REF] Riviere | Discontinuous Galerkin finite element methods for dynamic linear solid viscoelasticity problems[END_REF].

A fully discrete formulation

Employing the Crank-Nicolson type time discretization, we can formulate a fully discrete problem for the non-hereditary terms. However, due to the weak singularity in the fractional order integral, it is necessary to also use appropriate numerical techniques for the discrete fractional integral.

A time step size ∆t = T /N > 0 is defined for some N ∈ N and we set t n = n∆t for n = 0, . . . , N . We denote the fully discrete solution for the velocity and the displacement by

W n h ∈ D k (E h ) and U n h ∈ D k (E h )
, respectively, for n = 0, . . . , N (i.e. we consider u(t n ) ≈ W n h and u(t n ) ≈ U n h for each time step). To complete the time discretization, a linear interpolation technique [START_REF] Li | Numerical approaches to fractional calculus and fractional ordinary differential equation[END_REF][START_REF] Jang | A priori error analysis for a finite element approximation of dynamic viscoelasticity problems involving a fractional order integro-differential constitutive law[END_REF] is introduced. A local interpolation operator L n is defined as follows:

L n (w)(t) = - t -t n ∆t w(t n-1 ) + t -t n-1 ∆t w(t n ) for n = 1, . . . , N.
This operator provides a piecewise linear interpolation of w. Proposition 4. If w is of C 2 in time, we can define

E n (t) := w(t) -L n (w)(t) = 1 2 ẅ(ξ t )(t -t n-1 )(t -t n ) for some ξ t ∈ [t n-1 , t n ],
where t ∈ [t n-1 , t n ], by Rolle's theorem. For any t ∈ [t n-1 , t n ], if w(t) ∈ X for a normed space X, it holds that

∥E n (t)∥ X ≤ ∆t 2 2 ∥ ẅ∥ L∞(tn-1,tn;X) .
This inequality also holds on the broken Sobolev norm sense.

Now, we can derive the numerical approximation q n (w) to the fractional order integration satisfying that

0 I 1-α t w(t n ) = 1 Γ(1 -α) n i=1 ti ti-1 (L i (w)(t ′ ) + E i (t ′ )) (t n -t ′ ) -α dt ′ , = ∆t 1-α Γ(3 -α) n i=0 B n,i w(t i ) + 1 Γ(1 -α) n i=1 ti ti-1 E i (t ′ )(t n -t ′ ) -α dt ′ := q n (w) + e n (w), (2.22)
where

B n,i =    n 1-α (2 -α -n) + (n -1) 2-α , i = 0, (n -i -1) 2-α + (n -i + 1) 2-α -2(n -i) 2-α , i = 1, . . . , n -1, 1, i = n.
Note that 0 < B n,i < 2 for any n > 0 and i = 0, . . . , n. By Proposition 4, if w ∈ C 2 (0, T ; X), the numerical approximation error by the linear interpolation is given by, for any n = 1, . . . , N ,

∥e n (w)∥ X = ∥ 0 I 1-α t w(t n ) -q n (w)∥ X ≤ ∆t 2 2Γ(1 -α) ∥ ẅ∥ L∞(0,tn;X) tn 0 (t n -t ′ ) -α dt ′ , ≤ ∆t 2 2Γ(1 -α) ∥ ẅ∥ L∞(0,T ;X) T 1-α . (2.23)
For more details of the linear interpolation technique for fractional integral, we refer to [START_REF] Li | Numerical approaches to fractional calculus and fractional ordinary differential equation[END_REF][START_REF] Jang | Spatially continuous and discontinuous Galerkin finite element approximations for dynamic viscoelastic problems[END_REF][START_REF] Jang | A priori error analysis for a finite element approximation of dynamic viscoelasticity problems involving a fractional order integro-differential constitutive law[END_REF] and the references therein.

In the fully discrete sense, we define

Q n (W h ) by Q n (W h ) = ∆t 1-α Γ(3 -α) n i=0 B n,i W i h .
For simplicity, we assume ψ 0 = 0 which implies that the linear form F (t; •) is well-defined at t = 0. Finally, we can formulate a fully discrete problem by using the Crank-Nicolson type scheme and the approximate fractional integral as follows: find

W n h ∈ D k (E h ) and U n h ∈ D k (E h ) for n = 0, . . . , N such that for any v ∈ D k (E h ), ρ W n+1 h -W n h ∆t , v + φ 0 a U n+1 h + U n h 2 , v + φ α a Q n+1 (W h ) + Q n (W h ) 2 , v + J γ0,γ1 0 W n+1 h -W n h ∆t , v = 1 2 (F (t n+1 ; v) + F (t n ; v)), (2.24)
for n = 0, . . . , N -1, and

a U 0 h , v =a (u 0 , v) , (2.25) a W 0 h , v =a (w 0 , v) , (2.26)
with

W n+1 h + W n h 2 = U n+1 h -U n h ∆t . (2.27)
The jump penalty term for the discrete acceleration (the fourth term in (2.24)) is not necessary for stability analysis. But we will require it for energy error estimation for the velocity later.

Remark. In the previous work [START_REF] Riviere | Discontinuous Galerkin finite element methods for dynamic linear solid viscoelasticity problems[END_REF][START_REF] Jang | Spatially continuous and discontinuous Galerkin finite element approximations for dynamic viscoelastic problems[END_REF][START_REF] Jang | A priori analysis of a symmetric interior penalty discontinuous Galerkin finite element method for a dynamic linear viscoelasticity model[END_REF], the SIPG formulations for linear viscoelasticity problems introduced jump penalty terms for discrete velocity to handle the spatial discontinuity of the discrete velocity. On the other hand, our method introduces the jump penalty term of acceleration in (2.24), providing the capability to regulate the discontinuity of discrete acceleration.

For instance, as detailed in [START_REF] Jang | A priori analysis of a symmetric interior penalty discontinuous Galerkin finite element method for a dynamic linear viscoelasticity model[END_REF], an additional jump penalty term is defined as

J γ0,γ1 0 W n+1 h + W n h 2 , v ,
involving the time-averaged velocity. While this approach allows us to have the boundedness of J γ0,γ1 0 (W n h , W n h ), it may not ensure the management of acceleration's discontinuity. Conversely, our scheme,

J γ0,γ1 0 W n+1 h -W n h ∆t , v ,
can reduce the discontinuity of acceleration by adjusting penalty parameters, immediately.

Stability and convergence analyses

In this section, we present a stability bound as well as an error bound for the fully discrete problem. In the fully discrete problem, the stability bound implies the existence and uniqueness of the solution. To estimate a priori error, we introduce the usual approach using the elliptic projection.

A stability bound

In the semi-discrete problem, we showed the stability analysis without using Grönwall's inequality to yield non-exponentially increasing bounds in time. Using the positive definiteness (1.4), we could deal with the fractional integration term but in the following case of the fully discrete scheme, we need to use a different proof technique based on mathematical induction. For the stability analysis of the fully discrete problem, we also suppose the same smooth data terms. Theorem 2. If γ 1 (d -1) ≥ 1 and γ 0 is sufficiently large, there exists a unique discrete solution to (2.24)-(2.26) that satisfies

max 0≤n≤N W n h 2 L2(Ω) + max 0≤n≤N U n h 2 V + ∆t 2-α N -1 n=0 W n+1 h + W n h 2 V + max 0≤n≤N J γ0,γ1 0 (W n h , W n h ) ≤C ∥w 0 ∥ 2 H 2 (Ω) + ∥u 0 ∥ 2 H 2 (Ω) + ∆t 2-α ∥w 0 ∥ 2 H 2 (Ω) + T ∆t N n=0 f (t n ) 2 L2(Ω) + ∆t N n=0 h -1 g N (t n ) 2 L2(Γ N )
, where C is independent of the solution, ∆t and h.

Proof. Let m be an arbitrary positive integer such that m < N . [START_REF] Rivière | Discontinuous Galerkin finite element methods for dynamic linear solid viscoelasticity problems[END_REF]) and summing all results over n = 0 to n = m -1, we get

Taking v = 2∆t(W n+1 h + W n h ) for 0 ≤ n ≤ m -1 in (2.
2ρ W m h 2 L2(Ω) -W 0 h 2 L2(Ω) + φ 0 ∆t m-1 n=0 a U n+1 h + U n h , W n+1 h + W n h + φ α ∆t m-1 n=0 a Q n+1 (W h ) + Q n (W h ), W n+1 h + W n h + 2 J γ0,γ1 0 (W m h , W m h ) -J γ0,γ1 0 W 0 h , W 0 h =∆t m-1 n=0 f (t n+1 ) + f (t n ), W n+1 h + W n h + ∆t m-1 n=0 g N (t n+1 ) + g N (t n ), W n+1 h + W n h L2(Γ N ) . ( 3
+ W n h ) = 2(U n+1 h -U n h ), we can rewrite (3.1) as 2ρ W m h 2 L2(Ω) + 2φ 0 a (U m h , U m h ) + φ α ∆t m-1 n=0 a Q n+1 (W h ) + Q n (W h ), W n+1 h + W n h + 2J γ0,γ1 0 (W m h , W m h ) =∆t m-1 n=0 f (t n+1 ) + f (t n ), W n+1 h + W n h + ∆t m-1 n=0 g N (t n+1 ) + g N (t n ), W n+1 h + W n h L2(Γ N ) + 2ρ W 0 h 2 L2(Ω) + 2φ 0 a U 0 h , U 0 h + 2J γ0,γ1 0 W 0 h , W 0 h . (3.2)
By the definition of the numerical approximation of the fractional integration, we can expand the third term of (3.2) as

φ α ∆t 2-α Γ(3 -α) m-1 n=0 a n+1 i=0 B n+1,i W i h + n i=0 B n,i W i h , W n+1 h + W n h ,
and we can split it by

φ α ∆t 2-α Γ(3 -α) m-1 n=0 a W n+1 h + W n h , W n+1 h + W n h + φ α ∆t 2-α Γ(3 -α) m-1 n=0 a n i=0 B n+1,i W i h + n-1 i=0 B n,i W i h , W n+1 h + W n h ,
since B n,n = 1 for ∀n. By introducing this expression into (3.2), we can obtain

2ρ W m h 2 L2(Ω) + 2κφ 0 U m h 2 V + κφ α ∆t 2-α Γ(3 -α) m-1 n=0 W n+1 h + W n h 2 V + 2J γ0,γ1 0 (W m h , W m h ) ≤∆t m-1 n=0 f (t n+1 ) + f (t n ), W n+1 h + W n h + ∆t m-1 n=0 g N (t n+1 ) + g N (t n ), W n+1 h + W n h L2(Γ N ) + 2ρ W 0 h 2 L2(Ω) + 2φ 0 a U 0 h , U 0 h + 2J γ0,γ1 0 W 0 h , W 0 h - φ α ∆t 2-α Γ(3 -α) m-1 n=0 a n i=0 B n+1,i W i h + n-1 i=0 B n,i W i h , W n+1 h + W n h , (3.3)
by the coercivity of the bilinear form. Now, we shall present the upper bounds of the first four terms in the right-hand side of (3.3) by following the proof of Theorem 1, e.g. Cauchy-Schwarz inequality, Young's inequality, inverse polynomial trace inequality, etc.

• ∆t m-1 n=0 f (t n+1 ) + f (t n ), W n+1 h + W n h
By the Cauchy-Schwarz inequality and the Young's inequality, we have

∆t m-1 n=0 f (t n+1 ) + f (t n ), W n+1 h + W n h ≤∆t m n=0 2ϵ a f (t n ) 2 L2(Ω) + 2 ϵ a W n h 2 L2(Ω) ≤∆t N n=0 2ϵ a f (t n ) 2 L2(Ω) + 4T ϵ a max 0≤n≤N W n h 2 L2(Ω)
for any positive ϵ a .

• ∆t

m-1 n=0 g N (t n+1 ) + g N (t n ), W n+1 h + W n h L2(Γ N )
Similarly, we can also derive, with the inverse polynomial trace theorem,

∆t m-1 n=0 g N (t n+1 ) + g N (t n ), W n+1 h + W n h L2(Γ N ) ≤ ∆t m n=0 2ϵ b g N (t n ) 2 L2(Γ N ) + ∆t m n=0 2 ϵ b e⊂Γ N W n h 2 L2(e) ≤ ∆t m n=0 2ϵ b g N (t n ) 2 L2(Γ N ) + ∆t m n=0 2Ch -1 ϵ b E∈E h W n h 2 L2(E) ≤ ∆t N n=0 2ϵ b g N (t n ) 2 L2(Γ N ) + Ch -1 4T ϵ b max 0≤n≤N W n h 2 L2(Ω) ,
for any positive ϵ b .

• W 0 h 2 L2(Ω) , a U 0 h , U 0 h and J γ0,γ1 0 W 0 h , W 0 h
As seen in (2.17) and (2.18), (2.25) and (2.26) imply that

a U 0 h , U 0 h ≤ K U 0 h 2 V ≤ C ∥u 0 ∥ 2 H 2 (Ω) , and W 0 h 2 L2(Ω) ≤ C W 0 h 2 V ≤ C ∥w 0 ∥ 2 H 2 (Ω) ,
by the continuity, respectively. Also, the definition of DG energy norm gives

J γ0,γ1 0 W 0 h , W 0 h ≤ W 0 h 2 V ≤ C ∥w 0 ∥ 2 H 2 (Ω) .
Tidying up all results with (3.3), we then have

2ρ W m h 2 L2(Ω) + 2κφ 0 U m h 2 V + κφ α ∆t 2-α Γ(3 -α) m-1 n=0 W n+1 h + W n h 2 V + 2J γ0,γ1 0 (W m h , W m h ) ≤2Cρ ∥w 0 ∥ 2 H 2 (Ω) + Cφ 0 ∥u 0 ∥ 2 H 2 (Ω) + 4T ϵ a + Ch -1 4T ϵ b max 0≤n≤N W n h 2 L2(Ω) + ∆t N n=0 2ϵ a f (t n ) 2 L2(Ω) + ∆t N n=0 2ϵ b g N (t n ) 2 L2(Γ N ) - φ α ∆t 2-α Γ(3 -α) m-1 n=0 a n i=0 B n+1,i W i h + n-1 i=0 B n,i W i h , W n+1 h + W n h :=R - φ α ∆t 2-α Γ(3 -α) m-1 n=0 a n i=0 B n+1,i W i h + n-1 i=0 B n,i W i h , W n+1 h + W n h . (3.4)
Note that in (3.4), R is positive and independent of m. Using mathematical induction, we want to show the following inequality holds such that

2ρ W m h 2 L2(Ω) + 2κφ 0 U m h 2 V + κφ α ∆t 2-α 2Γ(3 -α) m-1 n=0 W n+1 h + W n h 2 V + 2J γ0,γ1 0 (W m h , W m h ) ≤ C(R + ∆t 2-α W 0 h 2 V ). (3.5)
For m = 1 in (3.4), using the continuity of the bilinear form and Young's inequality, we have

2ρ W 1 h 2 L2(Ω) + 2κφ 0 U 1 h 2 V + κφ α ∆t 2-α Γ(3 -α) W 1 h + W 0 h 2 V + 2J γ0,γ1 0 W 1 h , W 1 h ≤R + φ α ∆t 2-α Γ(3 -α) K 2 B 2 1,0 ϵ 2 W 0 h 2 V + 1 2ϵ W 1 h + W 0 h 2 V , for any positive ϵ, since a B 1,0 W 0 h , W 1 h + W 0 h ≤ K 2 B 2 1,0 ϵ 2 W 0 h 2 V + 1 2ϵ W 1 h + W 0 h 2 V .
Taking ϵ = 1/κ, we can observe that (3.5) holds when m = 1. Let us assume that (3.5) holds for j < m. Consider (3.4) for j + 1, particularly the last term. Denoting 0 < G := max 0≤i≤n≤N B n,i < 2, we have

j n=0 a n i=0 B n+1,i W i h + n-1 i=0 B n,i W i h , W n+1 h + W n h ≤ j n=0 n-1 i=1 a B n,i W i+1 h + W i h , W n+1 h + W n h + j n=0 (B n,0 + B n+1,0 + B n+1,1 ) a W 0 h , W n+1 h + W n h + j n=0 B n+1,1 a W 1 h + W 0 h , W n+1 h + W n h ≤ j n=0 n-1 i=1 K 2 G 2 ϵ 2 W i+1 h + W i h 2 V + 1 2ϵ W n+1 h + W n h 2 V + j n=0 K 2 (3G) 2 ε 2 W 0 h 2 V + 1 2ε W n+1 h + W n h 2 V + j n=0 K 2 G 2 ε 2 W 1 h + W 0 h 2 V + 1 2ε W n+1 h + W n h 2 V
, for any positive ϵ, ε and ε. Then the induction assumption (3.5) for j implies the boundedness of

j-1 n=0 W n+1 h + W n h 2 V and j n=0 n-1 i=1 W i+1 h + W i h 2 V .
Therefore, choosing proper Young's constants such as ϵ = ε = ε = 1/(3κ) leads us to prove (3.5) for j + 1. Consequently, (3.5) holds for arbitrary m.

After noting that W 0

h 2 V ≤ C ∥w 0 ∥ 2 H 2 (Ω)
, combining it with the maximum argument such that

a n + b n ≤ C, ∀n, ⇒ max n a n + max n b n ≤ 2C, (3.5) yields 2ρ max 0≤n≤N W n h 2 L2(Ω) + 2κφ 0 max 0≤n≤N U n h 2 V + κφ α ∆t 2-α 2Γ(3 -α) N -1 n=0 W n+1 h + W n h 2 V + 2 max 0≤n≤N J γ0,γ1 0 (W n h , W n h ) ≤4C 2ρ ∥w 0 ∥ 2 H 2 (Ω) + Cφ 0 ∥u 0 ∥ 2 H 2 (Ω) + 4T ϵ a + Ch -1 4T ϵ b max 0≤n≤N W n h 2 L2(Ω) + ∆t N n=0 2ϵ a f (t n ) 2 L2(Ω) + ∆t N n=0 2ϵ b g N (t n ) 2 L2(Γ N ) + C∆t 2-α ∥w 0 ∥ 2 H 2 (Ω) .
In the end, we can complete our proof by taking appropriate Young's constants, for example, ϵ a = 32CT /ρ and ϵ b = 32CT /(ρh).

In the fully discrete case, the stability bound in Theorem 2 implies the existence and uniqueness of the fully discrete problem. The proof uses the maximum argument instead of discrete Grönwall's inequality to estimate the discrete integration in time. As a consequence, the stability bound is not exponentially increasing in time. While we have applied the positive definiteness of a fractional integration to the stability analysis of the semi-discrete problem, by employing mathematical induction, we have proved the stability bound of the fully discrete solution and the (discrete) fractional order integration of the velocity. Indeed, in a similar way to the positive definiteness on the fractional order integral, there exist positivity properties of the quadrature Q n . For instance, it holds that

∆t n i=0 Q n (Φ)Φ n ≥ 0, ∀Φ = (Φ 0 , . . . , Φ n ) T . (3.6)
However, to use the positivity (3.6) in the stability analysis, it is essential to assume the homogeneous Dirichlet boundary condition, since a spectral decomposition is additionally required in the proof. Furthermore, using (3.6) will not provide energy norm bounds of velocity as in Theorem 1. For more details, we refer to [START_REF] Mclean | Numerical solution of an evolution equation with a positive-type memory term[END_REF][START_REF] Mclean | Numerical solution via Laplace transforms of a fractional order evolution equation[END_REF].

Error estimates

Following the usual path for error estimation [START_REF] Jang | A priori analysis of a symmetric interior penalty discontinuous Galerkin finite element method for a dynamic linear viscoelasticity model[END_REF], we want to split spatial and temporal error components by introducing the DG elliptic projection. We define

θ(t) := u(t) -Ru(t), χ n := U n h -Ru(t n ), ϖ n := W n h -R u(t n ), E 1 (t) := ü(t + ∆t) + ü(t) 2 - u(t + ∆t) -u(t) ∆t ,
for t ∈ [0, T ] and n = 0, . . . , N . Then (2.27) gives

χ n+1 -χ n ∆t = ϖ n+1 + ϖ n 2 + E 2 (t n ) + E 3 (t n ), (3.7)
for n = 0, . . . , N -1 where

E 2 (t) := θ(t + ∆t) -θ(t) ∆t - θ(t + ∆t) + θ(t) 2 and E 3 (t) := u(t + ∆t) + u(t) 2 - u(t + ∆t) -u(t) ∆t .
For a three-times time-differentiable function, v(t), when we denote its third time derivative by v (3) , we have

v(t n+1 ) + v(t n ) 2 - v(t n+1 ) -v(t n ) ∆t = 1 2∆t tn+1 tn v (3) (t)(t n+1 -t)(t -t n )dt. (3.8)
Hence, if v (3) ∈ L 2 (t n , t n+1 ; X), the Cauchy-Schwarz inequality gives

v(t n+1 ) + v(t n ) 2 - v(t n+1 ) -v(t n ) ∆t 2 X ≤ ∆t 3 4 ∥v (3) ∥ 2 L2(tn,tn+1;X) . (3.9)
Before evaluating error estimates, we shall consider the regularity of the solution to ensure optimal errors in time. As seen in (3.9), the Crank-Nicolson method requires H 3 smoothness in time. Due to the weak singularity in fractional order integrals, it is not trivial to have second order accuracy in time.

To describe the regularity of solutions, we first introduce a convolution form in fractional integrals and a spatial differential operator. For instance, when we define β 1-α (t) = t -α /Γ(1 -α) and A = ∇ • Dε, and we denote a Laplace convolution by * , in the strong form (2.2), we can rewrite it as

ρü(t) = φ 0 Au(t) + φ α β 1-α * A u(t) + f (t), (3.10)
with the assumption ψ 0 = 0 for simplicity. In fact, to obtain ψ 0 = 0, we need to suppose u 0 ∈ ker(A). When we consider L q norm in time, Young's inequality for convolution yields ∥ρü∥ Lq(0,T ) ≤ φ 0 ∥Au∥ Lq(0,T ) + φ α ∥β 1-α ∥ L1(0,T ) ∥A u∥ Lq(0,T ) + ∥f ∥ Lq(0,T ) , (3.11) for q ≥ 1. In (3.11), ü is L 2 integrable in time if Au, A u and f are L 2 integrable in time. By differentiation of (3.10) in time, we have ρu (3)

(t) = φ 0 A u(t) + φ α β 1-α (t)Aw 0 + φ α β 1-α * Aü(t) + ḟ (t). (3.12)
Since β 1-α is L 1 integrable, so is u (3) with L 1 integrable A u, Aü and ḟ in time. However, β 1-α is only L 2 integrable for α < 1/2. Hence, if α < 1/2 or w 0 ∈ ker(A), u (3) is L 2 integrable in time where A u, Aü, ḟ ∈ L 2 (0, T ). Repeating this process, we can consider the fourth time derivative of u such that

ρu (4) (t) = φ 0 Aü(t) + φ α β1-α (t)Aw 0 + φ α β 1-α (t)Aü(0) + φ α β 1-α * Au (3) (t) + f (t). (3.13)
Note that β1-α / ∈ L q (0, T ) for q ≥ 1, for any 0 < α < 1. Hence, for the L q integrability of u (4) , it is required to assume that w 0 ∈ ker(A) and ü(0) ∈ ker(A).

Remark. For n ≥ 1, since we have u (3) ∈ L 2 (t n , t n+1 ; X) and u (4) ∈ L 2 (t n , t n+1 ; X) with sufficiently smooth initial data and source terms, substitutions of (3.12) and (3.13) into (3.8) lead to (3.9) for u and u, respectively. However, when n = 0, the singularity appears. Thus, in order to take full boundedness in time, we need further attention to manage the bound for n = 0. Lemma 1. Suppose f ∈ H 3 (0, T ; L 2 (Ω)) ∩ H 2 (0, T ; H 2 (Ω)), u 0 ∈ ker(A) and w 0 ∈ H 2 (Ω). Furthermore, we assume that u ∈ H 3 (0, T ; H 2 (Ω)) and Au ∈ H 3 (0, T ; H 2 (Ω)). If w 0 ∈ ker(A), there exists a positive constant C independent of ∆t such that ü(∆t

) + ü(0) 2 - u(∆t) -u(0) ∆t ≤ C∆t 2-α . (3.14)
Moreover, if w 0 ∈ ker(A) and ü(0) ∈ ker(A), we can also obtain .

ü(∆t) + ü(0) 2 - u(∆t) -u(0) ∆t ≤ C∆t 2 . ( 3
On the other hand, since Aü(0) is time independent, we have

1 ∆t ∆t 0 β 1-α (t)Aü(0)p 2 (t)dt = (1 -α)∆t 2-α Γ(4 -α) Aü(0).
Therefore, we can obtain from (3.17)

ü(∆t) + ü(0) 2 - u(∆t) -u(0) ∆t ≤C ∆t 2 + ∆t 2-α ≤ C∆t 2-α , (3.18)
where C depends on u and f but is independent of ∆t.

In addition, if w 0 ∈ ker(A) and ü(0) ∈ ker(A), it immediately implies that the error has second order accuracy.

We refer to [START_REF] Mclean | Numerical solution of an evolution equation with a positive-type memory term[END_REF] for the regularity of solutions of fractional order integro-differential equations, using a spectral decomposition for a representation with Fourier coefficients. Remark. In our model problem, the weak singularity only occurs at t = 0. If the solution has sufficient regularity in time for the second order time discretization schemes, additional regularity properties in space will not be required. Otherwise, the solution needs higher regularity in space such that u (3) (t) ∈ H 4 (Ω) near t = 0 where the regularity conditions rely on initial conditions and f . For more details of regularity results with respect to time and space, and their assumptions, e.g. see [START_REF] Mclean | Numerical solution of an evolution equation with a positive-type memory term[END_REF].

Using the DG elliptic error estimates as in Proposition 1 and Crank-Nicolson type temporal errors (3.9) with Lemma 1, we can derive the following lemma. Lemma 2. Let γ 0 > 0 be large enough and

γ 1 (d -1) ≥ 1. Suppose u ∈ C 2 0, T ; C 2 (Ω) ∩ H s (E h ) ∩ W 3 ∞ (0, T ; H s (E h )),
and the data terms are sufficiently smooth to fulfil Lemma 1. For the fully discrete solutions to

(2.24)-(2.26), (W n h ) N n=0 and (U n h ) N n=0
, there exists a positive constant C such that

max 0≤n≤N ϖ n L2(Ω) + max 0≤n≤N χ n V + ∆t 2-α N -1 n=0 ϖ n+1 + ϖ n 2 V 1/2 + max 0≤n≤N (J γ0,γ1 0 (ϖ n , ϖ n )) 1/2 ≤ CT 2-α h r + ∆t 2-α ,
where r = min(k + 1, s). Here, C is independent of h, ∆t and the numerical solutions. In addition to the condition of smooth data, if ü(0) ∈ ker(A), we have

max 0≤n≤N ϖ n L2(Ω) + max 0≤n≤N χ n V + ∆t 2-α N -1 n=0 ϖ n+1 + ϖ n 2 V 1/2 + max 0≤n≤N (J γ0,γ1 0 (ϖ n , ϖ n )) 1/2 ≤ CT 2-α h r + ∆t 2 .
Proof. Let m be an arbitrary positive integer in [1, N ]. For any 0 ≤ n ≤ m -1, average of (2.9) over t n and t n+1 and subtraction of it from (2.24) give

ρ W n+1 h -W n h ∆t - ü(t n+1 ) + ü(t n ) 2 , v + φ 0 a U n+1 h + U n h 2 - u(t n+1 ) + u(t n ) 2 , v + φ α a Q n+1 (W h ) + Q n (W h ) 2 - 0 I 1-α tn+1 u + 0 I 1-α tn u 2 , v + J γ0,γ1 0 W n+1 h -W n h ∆t - ü(t n+1 ) + ü(t n ) 2 , v =0
ρ ∆t ϖ n+1 -ϖ n , v + φ 0 2 a χ n+1 + χ n , v + φ α 2 a Q n+1 (ϖ) + Q n (ϖ), v + 1 ∆t J γ0,γ1 0 ϖ n+1 -ϖ n , v = ρ ∆t θ(t n+1 ) -θ(t n ), v + φ 0 2 a (θ(t n+1 ) + θ(t n ), v) + φ α 2 a q n+1 ( θ) + q n ( θ), v + φ α 2 a (e n+1 ( u) + e n ( u), v) + ρ (E 1 (t n ), v) + 1 ∆t J γ0,γ1 0 θ(t n+1 ) -θ(t n ), v + J γ0,γ1 0 (E 1 (t n ), v) . (3.20)
The Galerkin orthogonality, the linearity of the bilinear form and the spatial continuity reduce (3.20) to

ρ ∆t ϖ n+1 -ϖ n , v + φ 0 2 a χ n+1 + χ n , v + φ α 2 a Q n+1 (ϖ) + Q n (ϖ), v + 1 ∆t J γ0,γ1 0 ϖ n+1 -ϖ n , v = ρ ∆t θ(t n+1 ) -θ(t n ), v + φ α 2 a (e n+1 ( u) + e n ( u), v) + ρ (E 1 (t n ), v) . (3.21) Since u(t) ∈ C 2 (Ω), e n ( u) is of C 2 in
the spatial domain too. Hence, the continuity implies that a (e n+1 ( u) + e n ( u), v) = (A(e n+1 ( u) + e n ( u)), v) .

To complete the proof, we will follow the same arguments in the stability analysis with the spatial error estimates as well as the time discretization errors. By substitution of v = 2∆t(ϖ n+1 + ϖ n ) into (3.21) with (3.7) and (3.22), summing together from n = 0 to n = m -1 yields

2ρ ϖ m 2 L2(Ω) + 2φ 0 a (χ m , χ m ) + φ α ∆t m-1 n=0 a Q n+1 (ϖ) + Q n (ϖ), ϖ n+1 + ϖ n + 2J γ0,γ1 0 (ϖ m , ϖ m ) =2ρ ϖ 0 2 L2(Ω) + 2φ 0 a χ 0 , χ 0 + 2J γ0,γ1 0 ϖ 0 , ϖ 0 + 2ρ m-1 n=0 θ(t n+1 ) -θ(t n ), ϖ n+1 + ϖ n + φ α ∆t m-1 n=0
A(e n+1 ( u) + e n ( u)),

ϖ n+1 + ϖ n + 2ρ∆t m-1 n=0 E 1 (t n ), ϖ n+1 + ϖ n + 2φ 0 ∆t m-1 n=0 a χ n+1 + χ n , E 2 (t n ) + 2φ 0 ∆t m-1 n=0 a χ n+1 + χ n , E 3 (t n ) . (3.23)
Expanding the discrete fractional integration Q, and using the coercivity in (3.23), we can derive

2ρ ϖ m 2 L2(Ω) + 2κφ 0 χ m 2 V + κφ α ∆t 2-α Γ(3 -α) m-1 n=0 ϖ n+1 + ϖ n 2 V + 2J γ0,γ1 0 (ϖ m , ϖ m ) ≤2ρ ϖ 0 2 L2(Ω) + 2φ 0 a χ 0 , χ 0 + 2J γ0,γ1 0 ϖ 0 , ϖ 0 + 2ρ m-1 n=0 θ(t n+1 ) -θ(t n ), ϖ n+1 + ϖ n + φ α ∆t m-1 n=0
A(e n+1 ( u) + e n ( u)),

ϖ n+1 + ϖ n + 2ρ∆t m-1 n=0 E 1 (t n ), ϖ n+1 + ϖ n + 2φ 0 ∆t m-1 n=0 a χ n+1 + χ n , E 2 (t n ) + 2φ 0 ∆t m-1 n=0 a χ n+1 + χ n , E 3 (t n ) - ∆t 2-α Γ(3 -α) m-1 n=0 a n i=0 B n+1,i ϖ i + n-1 i=0 B n,i ϖ i , ϖ n+1 + ϖ n . (3.24)
Next, we shall show the bounds of the right-hand side of (3.24) except for the last term.

• ϖ 0 2 L2(Ω) , a χ 0 , χ 0 and J γ0,γ1 (2.26) and the Garlerkin orthogonality. Then Poincaré's inequality (1.2) leads us to obtain

0 ϖ 0 , ϖ 0 For any v ∈ D k (E h ), we have a ϖ 0 , v = a W 0 h -w 0 , v + a (w 0 -Rw 0 , v) = 0 by
ϖ 0 2 L2(Ω) ≤ C ϖ 0 2 V = 0
and so J γ0,γ1 0 ϖ 0 , ϖ 0 ≤ ϖ 0 2 V = 0. In this manner, the Galerkin orthogonality and (2.25) gives

a χ 0 , χ 0 = 0. • m-1 n=0 θ(t n+1 ) -θ(t n ), ϖ n+1 + ϖ n
The use of Cauchy-Schwarz inequalities, Young's inequality and elliptic error estimates yields

m-1 n=0 θ(t n+1 ) -θ(t n ), ϖ n+1 + ϖ n = m-1 n=0 tn+1 tn θ(t ′ ), ϖ n+1 + ϖ n dt ′ ≤ ϵ a 2 tm 0 θ(t ′ ) 2 L2(Ω) dt ′ + ∆t 2ϵ a m-1 n=0 ϖ n+1 + ϖ n 2 L2(Ω) ≤ ϵ a 2 θ 2 L2(0,T ;L2(Ω)) + ∆t 2ϵ a 4N max 0≤n≤N ϖ n 2 L2(Ω) ≤C ∥ü∥ 2 L2(0,T ;H s (E h )) ϵ a 2 h 2r + 2T ϵ a max 0≤n≤N ϖ n 2 L2(Ω) ,
for any positive ϵ a and r = min(k + 1, s).

• ∆t m-1 n=0 A(e n+1 ( u) + e n ( u)), ϖ n+1 + ϖ n
In a similar way, we have

∆t m-1 n=0
A(e n+1 ( u) + e n ( u)),

ϖ n+1 + ϖ n ≤ ∆t N -1 n=0 ϵ b 2 A(e n+1 ( u)+e n ( u)) 2 L2(Ω) + 2T ϵ b max 0≤n≤N ϖ n 2 L2(Ω) ,
for any positive ϵ b . After noting that

A(e n+1 ( u) + e n ( u)) 2 L2(Ω) ≤ C ∥e n+1 ( u) + e n ( u)∥ 2 H 2 (Ω) , (2.23) allows us to derive ∆t m-1 n=0 A(e n+1 ( u) + e n ( u)), ϖ n+1 + ϖ n ≤ C∥u∥ 2 W 3 ∞ (0,T ;H 2 (Ω)) T 3-2α ϵ b ∆t 4 + 2T ϵ b max 0≤n≤N ϖ n 2 L2(Ω) . • ∆t m-1 n=0 E 1 (t n ), ϖ n+1 + ϖ n
Recalling Lemma 1 and the time discretization error (3.9) gives

∆t m-1 n=0 E 1 (t n ), ϖ n+1 + ϖ n ≤ ϵ c 8 ∆t 4 ∥u (4) ∥ 2 L2(t1,T ;L2(Ω)) + CT ϵ c ∆t 4-2α + 2T ϵ c max 0≤n≤N ϖ n 2 L2(Ω) ≤ CT ϵ c ∆t 4-2α + 2T ϵ c max 0≤n≤N ϖ n 2 L2(Ω)
by Cauchy-Schwarz inequalities and Young's inequality, for any positive

ϵ c . If ü(0) ∈ ker(A), it holds ∥u∥ H 4 (0,T ;H 2 (Ω)) < ∞ to yield ∆t m-1 n=0 E 1 (t n ), ϖ n+1 + ϖ n ≤ C∥u∥ 2 H 4 (0,T ;H 2 (Ω)) T ϵ c ∆t 4 + 2T ϵ c max 0≤n≤N ϖ n 2 L2(Ω) . • ∆t m-1 n=0 a χ n+1 + χ n , E 2 (t n ) and ∆t m-1 n=0 a χ n+1 + χ n , E 3 (t n )
Following similar arguments but using the continuity of the DG bilinear form rather than Cauchy-Schwarz inequality, we get

∆t m-1 n=0 a χ n+1 + χ n , E 2 (t n ) ≤ ϵ d 8K 2 ∆t 4 ∥θ (3) ∥ 2 L2(0,T ;V ) + 2T ϵ d max 0≤n≤N χ n 2 V ,
for any positive ϵ d . Note that the DG elliptic error estimates such as (2.6) provide spatial error estimates of θ and its time derivatives. For example, (2.6) implies θ (3)

(t) V ≤ C u (3) (t) H 1 (E h ) for any t. Also, we have ∆t m-1 n=0 a χ n+1 + χ n , E 3 (t n ) ≤C ϵ d 8K 2 ∆t 4 ∥u (3) ∥ 2 L2(0,T ;H 1 (E h )) + 2T ϵ d max 0≤n≤N χ n 2 V .
Collecting up all the above results, we can derive a bound for (3.24) as

2ρ ϖ m 2 L2(Ω) + 2κφ 0 χ m 2 V + κφ α ∆t 2-α Γ(3 -α) m-1 n=0 ϖ n+1 + ϖ n 2 V + 2J γ0,γ1 0 (ϖ m , ϖ m ) ≤C ρϵ a h 2r + φ α T 3-2α ϵ b ∆t 4 + ρT ϵ c ∆t 4-2α + φ 0 ϵ d ∆t 4 + 4ρT ϵ a + 2φ α T ϵ b + 2ρT ϵ c max 0≤n≤N ϖ n 2 L2(Ω) + 8φ 0 T ϵ d max 0≤n≤N χ n 2 V - ∆t 2-α Γ(3 -α) m-1 n=0 a n i=0 B n+1,i ϖ i + n-1 i=0 B n,i ϖ i , ϖ n+1 + ϖ n , (3.25)
where C is a positive constant independent of ∆t, h and the discrete solutions. We can observe that the bound in (3.25) except the last term is independent of m. As shown in the proof of Theorem 2, the use of mathematical induction and the maximum argument leads to

2ρ max 0≤n≤N ϖ n 2 L2(Ω) + 2κφ 0 max 0≤n≤N χ n 2 V + κφ α ∆t 2-α 2Γ(3 -α) N -1 n=0 ϖ n+1 + ϖ n 2 V + 2 max 0≤n≤N J γ0,γ1 0 (ϖ n , ϖ n ) ≤4C ρϵ a h 2r + φ α T 3-2α ϵ b ∆t 4 + ρT ϵ c ∆t 4-2α + φ 0 ϵ d ∆t 4 + 2T 2ρ ϵ a + φ α ϵ b + ρ ϵ c max 0≤n≤N ϖ n 2 L2(Ω) + 8φ 0 T ϵ d max 0≤n≤N χ n 2 V . (3.26)
At last, by the setting of coefficients of Young's inequalities,

ϵ a = 48CT, ϵ b = 24φ α CT, ϵ c = 24CT, ϵ d = 32CT /κ, we have ρ max 0≤n≤N ϖ n 2 L2(Ω) + κφ 0 max 0≤n≤N χ n 2 V + κφ α ∆t 2-α 2Γ(3 -α) N -1 n=0 ϖ n+1 + ϖ n 2 V + 2 max 0≤n≤N J γ0,γ1 0 (ϖ n , ϖ n ) ≤CT 4-2α (h 2r + ∆t 4 + ∆t 4-2α ), (3.27
) for some positive C which is independent of h, ∆t, T and the discrete solutions but depends on the strong solution, the data terms, the spatial domain and the material properties. Moreover, if we suppose ü(0) ∈ ker(A), the ∆t 4-2α terms in the bounds will disappear and (3.27) becomes

ρ max 0≤n≤N ϖ n 2 L2(Ω) + κφ 0 max 0≤n≤N χ n 2 V + κφ α ∆t 2-α 2Γ(3 -α) N -1 n=0 ϖ n+1 + ϖ n 2 V + 2 max 0≤n≤N J γ0,γ1 0 (ϖ n , ϖ n )
≤CT 4-2α (h 2r + ∆t 4 ).

(3.28) This completes the proof.

In Lemma 2, under the assumption of H 4 regularity in time, it is clear to see that the constant of the bound C depends on ∥u∥ H 4 (0,T ;H s (E h )) . On the other hand, if the solution has only H 3 regularity in time, Lemma 1 plays an important role to present sub-optimal orders of time discretization errors. Once H 4 regularity of the solution is possessed in time, the smoothness of f is no longer required. Only the condition of the continuous linear form is needed for the stability of discrete formulations.

By Lemma 2, we can prove the following a priori error estimates. N n=0 are the fully discrete solution. Then we can observe optimal orders of L 2 error estimates as well as energy error estimates with respect to h. Although we have sub-optimal 2 -α order accuracy in time, some additional conditions such as ü(0) ∈ ker(A) or H 4 regularity in time lead us to obtain second-order accuracy. Thus, we obtain max

0≤n≤N u(t n ) -W n h L2(Ω) ≤ CT 2-α h r + ∆t 2-α , max 0≤n≤N u(t n ) -U n h V ≤ CT 2-α h r-1 + ∆t 2-α ,
and with higher regularity in time,

max 0≤n≤N u(t n ) -W n h L2(Ω) ≤ CT 2-α h r + ∆t 2 , max 0≤n≤N u(t n ) -U n h V ≤ CT 2-α h r-1 + ∆t 2 ,
where r = min(k + 1, s) and C is a positive constant independent of h and ∆t.

Moreover, we can derive energy norm error estimates for the velocity as well as L 2 norm error estimates for the displacement:

max 0≤n≤N u(t n ) -W n h V ≤ CT 2-α h r-1 + ∆t 2-α , max 0≤n≤N u(t n ) -U n h L2(Ω) ≤ CT 2-α h r + ∆t 2-α ,
and if is H 4 regular in time max 0≤n≤N u(t n ) -W n h V ≤ CT 2-α h r-1 + ∆t 2 , max 0≤n≤N u(t n ) -U n h L2(Ω) ≤ CT 2-α h r + ∆t 2 .
Proof. Let us consider u(t n ) -W n h L2(Ω) for any n = 0, . . . , N . Using the triangular inequality, we have

u(t n ) -W n h L2(Ω) = θ(t n ) -ϖ n L2(Ω) ≤ θ(t n ) L2(Ω) + ϖ n L2(Ω)
. By (2.6) and Lemma 2, this immediately gives

u(t n ) -W n h L2(Ω) ≤ CT 2-α (h r + ∆t 2-α ), and since n is arbitrary, max 0≤n≤N u(t n ) -W n h L2(Ω) ≤ CT 2-α h r + ∆t 2-α . With higher regularity in time, we have max 0≤n≤N u(t n ) -W n h L2(Ω) ≤ CT 2-α h r + ∆t 2 .
In this manner, we take into account the energy norm error for the displacement. By the triangular inequality, the elliptic energy error estimates (2.6) and Lemma (2), we can obtain

u(t n ) -U n h V ≤ θ(t n ) V + χ n V ≤ CT 2-α (h r-1 + ∆t 2-α ), hence max 0≤n≤N u(t n ) -U n h V ≤ CT 2-α (h r-1 + ∆t 2-α ).
In addition, if the strong solution has H 4 regularity in time, it holds

max 0≤n≤N u(t n ) -U n h V ≤ CT 2-α (h r-1 + ∆t 2 ).
On the other hand, to show energy error estimates of the velocity, we need the inverse inequality (1.3). Then, after noting that

v 2 V = E∈E h E Dε(v) : ε(v) dE + J γ0,γ1 0 (v, v) ≤ C|||v||| 2 H 1 (E h ) + J γ0,γ1 0 (v, v) ⇒ v V ≤ 2C|||v||| H 1 (E h ) + 2 J γ0,γ1 0 (v, v),
the inverse inequality implies

u(t n ) -W n h V ≤ θ(t n ) V + ϖ n V ≤ θ(t n ) V + C|||ϖ n ||| H 1 (E h ) + 2 J γ0,γ1 0 (ϖ n , ϖ n ) ≤ θ(t n ) 2 V + Ch -1 ϖ n L2(Ω) + 2 J γ0,γ1 0 (ϖ n , ϖ n ).
Then, by employing (2.6) and Lemma 2, we can obtain the energy norm error bounds for the velocity.

In the case of L 2 error estimation of the displacement, Poincaré inequality (1.2) gives,

u(t n ) -U n h L2(Ω) ≤ θ(t n ) L2(Ω) + χ n L2(Ω) ≤ θ(t n ) L2(Ω) + C χ n V .
Therefore, we can complete the proof using (2.6) and Lemma 2.

Remark. In the context of the stability analysis, the jump penalty term of the discrete velocity in (2.24) may not directly contribute. However, it plays a crucial role in achieving energy error bounds for the velocity. This additional jump penalty term allows us to handle the spatial discontinuity of the error between the numerical velocity and the elliptic projection of the exact solution over the edges. It ensures the energy norm error estimates of the numerical velocity, which are essential for accurate and reliable numerical solutions.

Numerical experiments

Using the open-source finite element method library FEniCS of version 2019.1.0 (https://fenicsproject.org), we conduct numerical simulations to validate our error analysis. We consider two cases to demonstrate the influence of the regularity of solutions in time:

1. Example 1: This case involves a solution that is not of class H 4 in time, showcasing the performance of our method with less regular solutions. 2. Example 2: Here, we consider a smoother case with a solution having higher regularity in time.

Additionally, we provide Example 3 to demonstrate the practical applicability of our method using real-material data. The numerical simulations presented in this manuscript were implemented with code available on author Jang's Git repository (https://github.com/Yongseok7717/visco_frac_dg). We believe in the importance of open and reproducible research, and thus, we encourage readers to access and explore our code for a better understanding of our proposed approach.

Let e n u := u(t n ) -U n h and e n w := u(t n ) -W n h be the numerical errors of displacement and velocity, respectively. On account of the dependency of the DG energy norm on the penalty parameters, γ 0 and γ 1 , we consider the broken H 1 norm of errors instead. Thanks to Korn's inequality, the (broken) H 1 norm error estimates follow the same convergence rates as the DG energy error estimates. Therefore, by Theorem 3, for a solution with H 3 regularity in time and sufficient smoothness of f and initial conditions, the error estimates are as follows: ∀n,

• Displacement errors: |||e n u ||| H 1 (E h ) = O(h r-1 + ∆t 2-α
) and e n u L2(Ω) = O(h r + ∆t 2-α ).

• Velocity errors:

|||e n w ||| H 1 (E h ) = O(h r-1 + ∆t 2-α ) and e n w L2(Ω) = O(h r + ∆t 2-α ).
Here, r is the spatial convergence rate, α is the fractional order of the time derivative, and h and ∆t are the mesh sizes in space and time, respectively. A higher regularity of the solution will lead to the second order accuracy in time as the optimal result of the Crank-Nicolson Scheme. The numerical convergent rate can be computed by the differences between two errors divided by mesh differences in the logarithm. For example, the spatial order of convergence d h is obtained by

d h = log(error of h 1 ) -log(error of h 2 ) log(h 1 ) -log(h 2 ) ,
for different mesh sizes h 1 and h 2 , when the temporal errors are negligible. In this manner, we can derive a numerical order of convergence in time d t as well. This allows us to quantify how the error decreases as we refine the mesh or change the time step size, providing valuable insights into the accuracy and efficiency of our numerical method.

Remark. For the stability of our numerical scheme, we should take sufficiently large penalty parameters, since the coercivity, the continuity, the DG elliptic error estimates and the bounds for interior penalty rely on the penalty parameters. We refer to [START_REF] Jang | A priori analysis of a symmetric interior penalty discontinuous Galerkin finite element method for a dynamic linear viscoelasticity model[END_REF] for the failure of DG simulations when the penalty parameters are not large enough. In our following numerical experiments, we define γ 0 = 20 and γ 1 = 1 in 2D problems. Example 1. Let us consider an exact solution to the primal model problem in the strong form defined by u(t, x, y) = (0.5t 2 + 0.4t 2.5 ) sin(πx) sin(πy)

x(1 -x)y(1 -y) on [0, 1] × Ω,
with Ω = (0, 1) 2 and its boundary splitting in

Γ N := {(x, y) ∈ ∂Ω | x = 0} and Γ D := ∂Ω\Γ N . We set α = 1/2, ρ = 1, φ 0 = 0, φ 1 = 1/Γ(1/2) and Dε = ε so that u solves ü(t) -∇ • 0 I 1/2 t ε(t) = f (t), (4.1)
for some f that can be readily determined analytically (easy to compute fractional integrals of polynomials. e.g. 1/2 order integral of t k is Γ(k + 1)/Γ(k + 3/2)t k+1/2 ). Also, the traction g N (t) can be obtained from the exact solution. Here, we can observe that u ∈ C 3 (0, T ; C ∞ (Ω)), u(0) = u(0) = 0 and ü(0) / ∈ ker(A).

Note that u (4) is not integrable in time hence the numerical errors of Example 1 will follow 1.5 order of accuracy in time, i.e. suboptimal convergence in time. This example is equivalent to [START_REF] Jang | A priori error analysis for a finite element approximation of dynamic viscoelasticity problems involving a fractional order integro-differential constitutive law[END_REF]Example 5.1], where the model problem has reduced to a parabolic type evolution problem of fractional order viscoelasticity. We refer to the reference for the suboptimal numerical results by the continuous finite element method imposed by purely homogeneous Dirichlet boundary conditions. In contrast, our simulation utilizes DGFEM for spatial discretization and imposes a mixed boundary condition of a non-homogeneous Neumann boundary and a homogeneous Dirichlet boundary. Tables 1 and2 indicate numerical errors at the final time in H 1 norm and L 2 norm with fixed fine timesteps for linear and quadratic polynomial bases, respectively. Those errors exhibit optimal orders of convergence with respect to the spatial mesh h, i.e. d h ≈ k in H 1 norm and d h ≈ k + 1 in L 2 norm for both displacement and velocity, where k = 1 or 2. However, due to the weak singularity, the second-order schemes in Example 1 cannot fully exploit their second-order accuracy in time, resulting in suboptimal convergent rates in time as shown in Table 3, where d t ≈ 1.5 regardless of variables and norms. In the next example, we solve the fractional order viscoelasticity problem for smoother solutions than Example 1 so that error estimates will follow optimal convergence rates not only of h but also of ∆t. Example 2. Let us define u(t, x, y) = t 4.5 sin(πx) sin(πy)

x(1 -x)y(1 -y) on [0, 1] × Ω,
with the same parameters and domains setting in Example 1 but φ 0 = 1 so that we suppose u solves

ü(t) -∇ • ε(t) -∇ • 0 I 1/2 t ε(t) = f (t), (4.2)
where data terms are obtained from the exact solution u. Clearly, the strong solution satisfies the regularity in time and space for the optimal error estimation theorem such that u ∈ C 4 (0, T ; C ∞ (Ω)), and u(0) = u(0) = ü(0) = 0.

By following the error estimates theorem for smooth solutions, on account of the regularity of solutions, the loss of accuracy in time discretization will disappear in Example 2. The numerical errors will follow O(h k + ∆t 2 ) in H 1 norm and O(h k+1 + ∆t 2 ) in L 2 norm for both displacement and velocity, respectively. Therefore, we can observe the optimal orders of convergence with respect to h such that d h ≈ 1 or 2 with linear polynomial bases and d h ≈ 2 or 3 with quadratic polynomial bases, depending on the choice of norms, in Tables 4 and5. On the other hand, Table 6 illustrates the second order accuracy in time of our numerical scheme for fixed h with cubic polynomial bases. Example 3. According to the real material data of butyl rubber, butyl 70821, from [START_REF] Bagley | A theoretical basis for the application of fractional calculus to viscoelasticity[END_REF][START_REF] Makris | Three-dimensional constitutive viscoelastic laws with fractional order time derivatives[END_REF], we illustrate fractional order viscoelasticity behaviour of the butyl rubber in 2D. For example, we have the parameters of the material such that ρ = 920 kg m -3 † , α = 0.449, φ 0 = 0.685 MN m -2 , φ 1 = 1.37 MN m -2 , µ = 0.228 and λ = 0.456. We suppose Ω = (0, 2) × (0, 1), T = 0. where A = 1 MN m -2 , small ϵ > 0 and H is the Heaviside step function. To impose non-zero traction only at the beginning of the simulation, we assume ∆t > ϵ.

In this numerical simulation, we want to solve

ρü(t) -φ 0 ∇ • Dε(t) -φ α ∇ • 0 I 1-α t D ε(t) = 0, (4.3)
with the given boundary conditions and initial conditions, where Dε is computed by (Dε) ij = D ijkl ε kl = 2µε ij + λδ ij ε kk .

By the zero body force and initial conditions, the numerical solution will satisfy at least suboptimal order of accuracy in time, i.e. O(∆t 2-α ). In this numerical simulation, the non-zero traction force on Γ left appears only at the first time iteration and then there is no more loading in the system, i.e. zero body force, homogeneous Dirichlet boundary condition and homogeneous Neumann boundary condition.

To compare mechanical responses between elasticity and viscoelasticity, we solve a linear elastic model by setting φ 1 = 0, for simplicity. Hence, the constitutive equation of the linear elasticity model is given by σ elastic (t) = φ 0 Dε(t).

For space and time discretization, we define the piecewise quadratic DG finite element space of 60 × 30 uniform mesh resulting in right-angled triangles and the timestep ∆t = 1/1000.

In Figure 1, the physical properties of wave propagation are observed well with respect to elastic and viscoelastic problems. The left figures of Figure 1 exhibit characteristics of elastic waves, while the solution of the viscoelastic model shows large attenuation. For more details on the attenuation in the unified elastic-viscoelastic model, we refer to [START_REF] Wang | Generalized viscoelastic wave equation[END_REF].

Conclusion

In conclusion, this research presents a rigorous analysis of the discontinuous Galerkin finite element method for addressing complex challenges in fractional order viscoelasticity. We have developed a fully discrete numerical approach that ensures stability and provides reliable numerical solutions, incorporating the Crank-Nicolson time-stepping scheme. Theoretical error estimates have been derived, revealing optimal convergence rates in both space and time for sufficiently smooth solutions. Conversely, solutions lacking high regularity exhibit suboptimal convergence in time. Our extensive numerical experiments affirm the efficiency and effectiveness of the proposed DGFEM method. These numerical findings solidify theoretical error estimates, confirming the numerical reliability of our approach.

The proposed approach has proven to be a robust and efficient numerical tool, capable of accurately predicting viscoelastic behaviour, even for solutions with nonsmooth features and weak singularities. The ability to handle non-homogeneous Neumann boundary conditions adds to its versatility and practicality, making it suitable for a wide range of real-world applications. The method's adaptability to different boundary conditions and solution regularities makes it important in various engineering and scientific domains. Higher order methods for time discretization are of our future work with fast computations of the fractional order integral/differentiation for practical use.

  e ⊗ n e -(v| E e j )| e ⊗ where the normal vector n e is oriented from E e i to E e j and ⊗ is the outer product defined, for vectors a and b, by (a ⊗ b) mn = a m b n for m, n = 1, . . . , d. On the other hand, if e ⊂ ∂Ω and e ⊂ ∂E {v} = v| e , {v} = v| e , [v] = v| e • n e , and [v ⊗ n e ] = v| e ⊗ n e .

Theorem 3 .

 3 Assume that u, f and initial conditions are given satisfying Lemma 2, and (W n h ) N n=0 and (U n h )

  25, Γ D = {(x, y) ∈ ∂Ω | x = 2} and Γ N = ∂Ω\Γ D . We impose zero initial conditions, zero body force, and homogeneous boundary conditions on the Γ D (the right edge) and the Neumann boundary of the top and bottom edges. For the left edge, we define the traction g(x, t) by g(x, t) = A(H(t) -H(t -ϵ) 0 on x ∈ Γ left := {0} × [0, 1],

Figure 1 :

 1 Figure 1: Displacement snapshots in the x-direction: the discrete solution u 1 for the linear elastic problem on the left and the viscoelastic problem on the right.

  .[START_REF] Larsson | The continuous Galerkin method for an integro-differential equation modeling dynamic fractional order viscoelasticity[END_REF] Proof. Consider(3.13) and split it into two parts by L 2 integrable part l(t) and others. Then we have ρu(4) (t) = l(t) + φ α β1-α (t)Aw 0 + φ α β 1-α (t)Aü(0). (3.16) If w 0 ∈ ker(A), (3.16) implies u (4) is L 1 integrable in time. Let p 2 (t) = (∆t -t)t for t ∈ [0, ∆t]. By substitution (3.16) into (3.8), when w 0 ∈ ker(A), we derive the following equation ′ )p 3 (t ′ )dt ′ , where p 3 (t) = ∆t/2t 2 -t 3 /3, and so there exists a positive C such that

	ü(∆t) + ü(0) 2	-	u(∆t) -u(0) ∆t	=	1 2ρ∆t	0	∆t	l(t)p 2 (t)dt +	φ α 2ρ∆t	0	∆t	β 1-α (t)Aü(0)p 2 (t)dt.	(3.17)
	Using integration by parts, we get
			t					t
			0 l(t ∆t l(t)p 2 (t)dt = l(t)p 3 (t) -0
		0	l(t)p 2 (t)dt ≤ C∆t 3 l(∆t) L2(Ω) + ∥ l∥ L1(0,∆t;L2(Ω))

Table 1 :

 1 Example 1; numerical errors and spatial orders of convergence when k = 1 and ∆t = 1/512.

		H 1 norm error				L 2 norm error
	h	e N u	H 1 (E h )	Rate	e N w H 1 (E h )	Rate	e N u L2(Ω)	Rate	e N w L2(Ω)	Rate
	1/2	1.3147e+00		2.8811e+00		1.6518e-01		3.8249e-01
	1/4	7.4445e-01	0.82	1.6378e+00	0.81 5.2814e-02 1.65 1.2549e-01 1.61
	1/8	3.8444e-01	0.95	8.4734e-01	0.95 1.4469e-02 1.87 3.4810e-02 1.85
	1/16	1.9377e-01	0.99	4.2730e-01	0.99 3.7269e-03 1.97 9.0123e-03 1.95
	1/32	9.7095e-02	1.00	2.1416e-01	1.00 9.4061e-04 1.99 2.2769e-03 1.98

Table 2 :

 2 Example 1; numerical errors and spatial orders of convergence when k = 2 and ∆t = 1/512.

		H 1 norm error				L 2 norm error
	h	e N u	H 1 (E h )	Rate	e N w H 1 (E h )	Rate	e N u L2(Ω)	Rate	e N w L2(Ω)	Rate
	1/2	3.9823e-01		8.7083e-01		2.6754e-02		6.0143e-02
	1/4	1.1133e-01	1.84	2.4333e-01	1.84 3.6105e-03 2.89 7.9411e-03 2.92
	1/8	2.8963e-02	1.94	6.3316e-02	1.94 4.6371e-04 2.96 1.0064e-03 2.98
	1/16	7.3409e-03	1.98 1.6054e-022 1.98 5.8663e-05 2.98 1.2665e-04 2.99
	1/32	1.8446e-03	1.99	4.0350e-03	1.99 7.5356e-06 2.96 1.6114e-05 2.97

Table 3 :

 3 Example 1; numerical errors and temporal orders of convergence when k = 3 and h = 1/64.

		H 1 norm error				L 2 norm error
	∆t	e N u	H 1 (E h )	Rate	e N w H 1 (E h )	Rate	e N u L2(Ω)	Rate	e N w L2(Ω)	Rate
	1/8	6.7401e-04		6.3975e-04		1.4615e-04		1.3384e-04
	1/16	2.7884e-04	1.27	2.9409e-04	1.12 6.0717e-05 1.27 6.4152e-05 1.10
	1/32	1.0582e-04	1.40	9.1038e-05	1.69 2.3073e-05 1.40 1.9915e-05 1.69
	1/64	3.8465e-05	1.46	2.7789e-05	1.71 8.3728e-06 1.46 5.9461e-06 1.74
	1/128	1.3947e-05	1.46	1.0314e-05	1.43 2.9817e-06 1.49 1.8027e-06 1.72

Table 4 :

 4 Example 2; numerical errors and spatial orders of convergence when k = 1 and ∆t = 1/512.

		H 1 norm error				L 2 norm error
	h	e N u	H 1 (E h )	Rate	e N w H 1 (E h )	Rate	e N u L2(Ω)	Rate	e N w L2(Ω)	Rate
	1/2	1.4802e+00		6.6002e+00		1.7595e-01		8.1404e-01
	1/4	8.3331e-01	0.83	3.7310e+00	0.82 5.4569e-02 1.69 2.5799e-01 1.66
	1/8	4.2932e-01	0.96	1.9255e+00	0.95 1.4649e-02 1.90 7.0264e-02 1.88
	1/16	2.1627e-01	0.99	9.7035e-01	0.99 3.7427e-03 1.97 1.8056e-02 1.96
	1/32	1.0835e-01	1.00	4.8620e-01	1.00 9.4141e-04 1.99 4.5539e-03 1.99

Table 5 :

 5 Example 2; numerical errors and spatial orders of convergence when k = 2 and ∆t = 1/512. H 1 norm error L 2 norm error

	h	e N u	H 1 (E h )	Rate	e N w H 1 (E h )	Rate	e N u L2(Ω)	Rate	e N w L2(Ω)	Rate
	1/2	4.4695e-01		1.9979e+00		2.9012e-02		1.3276e-01
	1/4	1.2494e-01	1.84	5.5857e-01	1.84 4.0018e-03 2.86 1.8044e-02 2.88
	1/8	3.2497e-02	1.94	1.4530e-01	1.94 5.2061e-04 2.94 2.3277e-03 2.95
	1/16	8.2348e-03	1.98	3.6826e-02	1.98 6.6068e-05 2.98 2.9492e-04 2.98
	1/32	2.0690e-03	1.99	9.2531e-03	1.99 8.4004e-06 2.98 3.7547e-05 2.97

Table 6 :

 6 Example 2; numerical errors and temporal orders of convergence when k = 3 and h = 1/64.

		H 1 norm error				L 2 norm error
	∆t	e N u	H 1 (E h )	Rate	e N w H 1 (E h )	Rate	e N u L2(Ω)	Rate	e N w L2(Ω)	Rate
	1/8	2.3658e-02		9.6445e-02		6.1306e-03		1.9759e-02
	1/16	5.6874e-03	2.01	2.4780e-02	1.96 1.4761e-03 2.05 5.1026e-03 1.95
	1/32	1.3851e-03	2.04	6.2993e-03	1.98 3.5974e-04 2.04 1.2997e-03 1.97
	1/64	3.4020e-04	2.03	1.5921e-03	1.98 8.8378e-05 2.03 3.2885e-04 1.98
	1/128	8.4070e-05	2.02	4.0123e-04	1.99 2.1827e-05 2.02 8.2883e-05 1.99

† https://www.aqua-calc.com/page/density-table/substance/rubber-coma-and-blank-butyl