# Locally implicit and stabilized explicit time schemes for transient visco-elastic wave propagation problems

Vinduja Vasanthan, Alexandre Imperiale, Sébastien Imperiale

<figure><img src="image_1.png" /></figure>

## To cite this version:

Vinduja Vasanthan, Alexandre Imperiale, Sébastien Imperiale. Locally implicit and stabilized explicit time schemes for transient visco-elastic wave propagation problems. 2025. hal-05003766

HAL Id: hal-05003766

https://inria.hal.science/hal-05003766v1

Preprint submitted on 24 Mar 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Locally implicit and stabilized explicit time schemes for transient visco-elastic
wave propagation problems

Vinduja Vasanthan<sup>a</sup>, Alexandre Imperiale<sup>a,*</sup>, Sébastien Imperiale<sup>b</sup>

<sup>a</sup>Université Paris-Saclay, CEA, List, F-91120, Palaiseau, France

<sup>b</sup>Inria, Ecole Polytechnique, CNRS, Palaiseau, France

**Abstract**

In the context of numerical methods for time-domain wave propagation problems, combining high-order lumped finite elements with an explicit time scheme is a popular approach for either inviscid or visco-elastic models. This strategy has proven to be efficient in numerous cases. However, when dealing with non-uniform meshes or high-contrast materials, the stability condition on the time step becomes drastically stringent. One can encounter such configurations when meshing “unfortunate” CAD input file, *e.g.* when dealing with heterogeneous materials where neighboring heterogeneities produce very small elements in-between them, or when considering materials with high and localized wave velocities. To address efficiently these configurations, we propose to adapt the locally implicit and stabilized leapfrog methods to the Kelvin-Voigt, Maxwell and Zener visco-elastic models. We prove using energy arguments that the global stability condition of these schemes can be much more favorable compared to a fully-explicit scheme, decreasing the number of iterations for a fixed time window. We illustrate our approaches with 2D and 3D numerical test cases related to ultrasonic non-destructive testing experiments.

*Keywords:* Time-domain wave problems, Visco-elasticity, Implicit-Explicit Schemes, Stabilized leapfrog schemes

# 1. Introduction

In the context of Non-Destructive Evaluation (NDE), Ultrasonic Testing (UT) is popular across numerous industrial fields and consists in examining the internal state of parts without damaging their integrity by propagating elastic and acoustic waves within them. To assist this process, building efficient numerical schemes for high-frequency transient wave propagation to model the UT experiments is relevant [1] and can be challenging depending on the considered configuration. In particular, visco-elastic wave propagation modeling is used in NDE to represent either intrinsic attenuation or multiple reflections occurring at sub-scales. Such phenomena are typically present in micro-structured materials, *e.g.* composite materials or welded metals. Different approaches exist to model such events, such as asymptotic methods (ray tracing) or numerical methods. In this work, we focus on a popular numerical approach: the time domain Spectral Element Method (SEM) [2, 3, 4], which we consider for visco-elastic wave propagation problems, namely the Maxwell (MX), Zener (ZN) and Kelvin-Voigt (KV) models – see *e.g.* the text book [5] and references therein. Since these are large scale propagation problems, *i.e.* the number of wavelengths within the propagation domain is (very) high, traditional time domain numerical methods resort to explicit time discretization. Combined with the so-called mass lumping strategy, the SEM is a very efficient approach and has successfully been applied in the context of visco-elastic problems [6].

In our work we consider cases where the propagation domain $\Omega$ is separated into two non-overlapping domains $\Omega_0$ and $\Omega_*$, such that $\bar{\Omega} = \bar{\Omega}_0 \cup \bar{\Omega}_*$. The restriction of the differential operators onto $\Omega_*$ is a

<footnote>*Corresponding author
Email address: alexandre.imperiale@cea.fr (Alexandre Imperiale)</footnote>
significantly stiffer part of the operators appearing in the visco-elastic propagation models. For example, this configuration can occur in the presence of an inclusion $\Omega_*$ with higher contrast compared to a background material occupying the volume $\Omega_0$. Another typical example consists in configurations where the mesh supporting the space discretization procedure is significantly more refined or non-uniform in the domain $\Omega_*$ – due, for example, to a specific shape of the domain (thin layer, or complex geometric features). It is well-known that explicit discretizations for hyperbolic problems have a stability condition that depends on the spectral radius of the operator [7, 8]. Hence, in the aforementioned configurations, the performance of the explicit time marching algorithm decreases drastically due to very restrictive time steps.

In order to overcome this difficulty, we propose in our work two strategies that we develop for the MX, ZN and KV models. A first one consists in using an implicit time scheme exclusively on the stiffest part of the operators. This leads to a Locally Implicit (LI) [9, 10, 11, 12, 13] time scheme with a stability condition independent of the restrictions of the operators on $\Omega_*$, hence a potentially great improvement on the time step. However, this gain is to be balanced with the cost of solving a linear system at each time step. The dimension of this linear system is equal to the number of degrees of freedom within $\Omega_*$, which can become quite large, especially for three-dimensional problems.

As an alternative strategy, we extend the stabilized explicit approach proposed in [14, 15] to visco-elastic problems. In essence, it aims at considering the mortar element domain decomposition approach [16, 17, 18] between $\Omega_0$ and $\Omega_*$, where in $\Omega_*$ the stiffness operator is replaced with a Chebychev polynomial. This method can lead to a gain in the time step proportionate to the order of this stabilization polynomial. This gain also needs to be balanced with the number of application of the stiffness operator in the final time marching algorithm, as well as the inversion of a linear system at the interface between $\Omega_0$ and $\Omega_*$.

This paper is outlined as follows. In Section 2, the groundwork to obtain a space-time discretization in a spectral finite element framework is laid for each considered visco-elastic model. In Section 3, we propose locally implicit schemes for each visco-elastic models, and we analyze their corresponding stability conditions using energy arguments. We show in particular that the stability condition on the time step is independent of the stiffness operator in $\Omega_*$. In Section 4, for each visco-elastic models, we extend an alternative approach based on the combination of stabilized leapfrog schemes and the mortar domain decomposition method. The stability conditions of the obtained schemes are analyzed, again using energy arguments. In Section 5 and 6, we study the precision of the schemes through a numerical convergence illustration, we present 2D and 3D numerical applications of each method linked to UT experiments and analyze the performance by quantifying the gain in CPU time. Finally, we conclude and give some perspectives in Section 7.

## 2. Visco-elastic models, associated explicit schemes and positioning of the problem

### 2.1. Maxwell, Zener and Kelvin-Voigt models

We consider visco-elastic wave propagation within a material occupying the open volume $\Omega \subset \mathbb{R}^d$ (with $d=2$ in 2D and $d=3$ in 3D). The boundary of the domain is denoted by $\partial\Omega$ and its outer normal vector field by $n$. Denoting by $u$ the displacement vector field generated from a low-amplitude high-frequency source $f$, the wave propagation phenomena are classically represented through the following equilibrium relation

$$ \rho \partial_{tt}^2 u - \nabla \cdot \sigma = f, \quad \text{in } \Omega, \quad \sigma \cdot n = 0, \quad \text{in } \partial\Omega, \qquad (1) $$

where $\rho$ is the mass density of the material, and $\sigma$ is the Cauchy stress tensor field. This relation is completed with given initial conditions in displacement $u|_{t=0} = u_0$ and velocity $\partial_t u|_{t=0} = v_0$, and a constitutive stress-strain relation.

In this regard, we assume that the effective visco-elasticity behavior of the material has been characterized – e.g. through experimental studies – at a specific frequency of interest $\omega_*$, yielding complex valued elastic modulii $C_* + iD_*$, namely

$$ \hat{\sigma}(\omega_*) = (C_* + iD_*)\hat{\varepsilon}(\omega_*) \qquad (2) $$

The real part $C_*$ is a symmetric strictly positive fourth-order tensor. It represents the purely elastic component of the constitutive law. The imaginary part $D_*$ is also a symmetric strictly positive fourth-order
tensor. It embodies the effective viscosity at wavelength scale. Usually, materials of sufficient quality factor $Q_*$ = $D_*^{-1}C_*$ bear a dominant elastic part, i.e. $C_* \gg D_*$, allowing for actual propagation of waves.

In our work we consider the three most popular MX, ZN and KV visco-elastic models – see e.g. textbook [5] and references therein for more details. We briefly recall hereafter the associated constitutive laws of these models. Also, following a rather straightforward identification process in the frequency domain, we review how each model can be tuned to exactly recover the complex constitutive law $C_*+iD_*$ at $\omega_*$ – see [6] for more details. By doing so, each model respects the same quality factor at $\omega_*$. This parametrization process is of capital importance when trying to adequately compare the three models in the same configuration setting.

**The Maxwell model.** The MX model is defined through the following constitutive relation

$$S_{\text{mx}} \partial_t \boldsymbol{\sigma} + \mathbf{H}_{\text{mx}} \boldsymbol{\sigma} = \partial_t \varepsilon(\mathbf{u}),$$

where $\varepsilon(\mathbf{u}) = \frac{1}{2}(\nabla\mathbf{u} + \nabla\mathbf{u}^\mathrm{T})$ is the linearized Green-Lagrange tensor, and $S_{\text{mx}}$ and $H_{\text{mx}}$ are homogeneous to a compliance tensor and the inverse of a viscosity tensor respectively. In other words, we define $C_{\text{mx}}$ and $D_{\text{mx}}$ two symmetric strictly positive fourth-order tensors such that

$$S_{\text{mx}} = (C_{\text{mx}})^{-1}, \quad \text{and} \quad H_{\text{mx}} = (D_{\text{mx}})^{-1}.$$

In the frequency domain, the constitutive relation reads

$$(i\omega S_{\text{mx}} + H_{\text{mx}})\hat{\boldsymbol{\sigma}} = i\omega\hat{\boldsymbol{\varepsilon}},$$

where $\hat{\boldsymbol{\sigma}}$ and $\hat{\boldsymbol{\varepsilon}}$ are the time-domain Fourier transform of $\boldsymbol{\sigma}$ and $\boldsymbol{\varepsilon}$ respectively. Identifying this relation with (2) leads to the relations

$$C_{\text{mx}} = C_* + D_* C_*^{-1} D_*, \quad \text{and} \quad D_{\text{mx}} = \frac{1}{\omega_*} (D_* + C_* D_*^{-1} C_*). \qquad (3)$$

As depicted in [6], Appendix A, the MX model follows a linear dependency of the quality factor w.r.t. the frequency, namely $Q_{\text{mx}}(\omega) = \omega D_{\text{mx}}^{-1}C_{\text{mx}}$, and such that $Q_{\text{mx}}(\omega_*) = Q_*$. The Zener model is defined through the following constitutive law

**The Zener model.** The ZN model is defined through the following constitutive law

$$\boldsymbol{\sigma} + \tau_{\text{zn}} \partial_t \boldsymbol{\sigma} = C_{\text{zn}} \boldsymbol{\varepsilon}(\mathbf{u}) + \tau_{\text{zn}} D_{\text{zn}} \partial_t \boldsymbol{\varepsilon}(\mathbf{u}),$$

where $\tau_{\text{zn}} > 0$ is a relaxation time, and $C_{\text{zn}}$ and $D_{\text{zn}}$ are symmetric strictly positive fourth-order tensors. In the frequency domain, the constitutive law is

$$(1 + i\omega\tau_{\text{zn}})\hat{\boldsymbol{\sigma}} = (C_{\text{zn}} + i\omega\tau_{\text{zn}}D_{\text{zn}})\hat{\boldsymbol{\varepsilon}},$$

and matching with (2) leads to

$$\tau_{\text{zn}} = \frac{1}{\omega_*}, \quad C_{\text{zn}} = C_* - D_*, \quad D_{\text{zn}} = C_* + D_* \qquad (4)$$

The associated quality factor is $Q_{\text{zn}}(\omega) = (D_{\text{zn}} - C_{\text{zn}})^{-1}((\omega\tau_{\text{zn}})^{-1}C_{\text{zn}} + \omega\tau_{\text{zn}}D_{\text{zn}})$ – see [5, 6] – and one recovers $Q_{\text{zn}}(\omega_*) = Q_*$. An interesting property of the Zener model, particularly sought in modeling viscosity for earth sciences [19], is that, in the high frequency regime, the quality factor is nearly constant at the frequency of interest

$$d_{\omega} Q_{\text{zn}}(\omega) = \tau_{\text{zn}} (D_{\text{zn}} - C_{\text{zn}})^{-1} (D_{\text{zn}} - (\omega \tau_{\text{zn}})^{-2} C_{\text{zn}}), \quad d_{\omega} Q_{\text{zn}}(\omega_*) = \frac{1}{\omega_*} I \underset{\omega_* \to +\infty}{\longrightarrow} 0,$$
where **I** is the fourth-order identity tensor. As detailed in [20], an alternate expression of the ZN model – typically more suited for energy analysis – is obtained by considering the purely viscous part of the stress tensor **s** = **σ** – **C**<sub>zn</sub>**ε**(**u**). Written in the variables (**u**, **s**) the complete ZN model reads

$$
\varrho \partial_{tt}^2 \mathbf{u} - \nabla \cdot \mathbf{s} - \nabla \cdot (\mathbf{C}_{\text{zn}} \varepsilon(\mathbf{u})) = \mathbf{f}, \quad \mathbf{S}_{\text{zn}} \partial_t \mathbf{s} + \mathbf{H}_{\text{zn}} \mathbf{s} = \partial_t \varepsilon(\mathbf{u}), \quad \text{in } \Omega,
$$

with adequate boundary and initial conditions, where $\mathbf{S}_{\text{zn}}$ and $\mathbf{H}_{\text{zn}}$ are defined through

$$
\mathbf{S}_{\text{zn}} = (\mathbf{D}_{\text{zn}} - \mathbf{C}_{\text{zn}})^{-1}, \quad \mathbf{H}_{\text{zn}} = \frac{1}{\tau_{\text{zn}}} \mathbf{S}_{\text{zn}}. \tag{5}
$$

The Kelvin-Voigt model. The KV model is the simplest visco-elastic model of the three considered in our work, and it actually benefits from a relative popularity when modelling visco-elastic waves propagation. The constitutive law is defined as

$$
\sigma = C_{kv} \epsilon(u) + D_{kv} \partial_t \epsilon(u),
$$

or, in frequency domain, $\hat{\sigma} = (C_{kv} + i\omega D_{kv})\hat{\varepsilon}$, where $C_{kv}$ and $D_{kv}$ are symmetric strictly positive fourth-order tensors. Matching with (2) leads to

$$
C_{kv} = C_{\star}, \quad D_{kv} = \frac{1}{\omega_{\star}} D_{\star}.
$$

The corresponding quality factor $Q_{kv}(\omega) = \omega^{-1} D_{kv}^{-1} C_{kv}$ is proportionate to the inverse of the frequency and satisfies $Q_{kv}(\omega_*) = Q_\star$. In practice, the stress unknown is discarded, yielding only the displacement unknown satisfying

$$
\varrho \partial_{tt}^2 \mathbf{u} - \nabla \cdot (\mathbf{C}_{\text{kv}} \varepsilon(\mathbf{u})) - \nabla \cdot (\mathbf{D}_{\text{kv}} \partial_t \varepsilon(\mathbf{u})) = \mathbf{f}, \quad \text{in } \Omega,
$$

with initial and boundary conditions.

2.2. Conform approximations and space-discrete operators

In our work, we consider exclusively primal-dual formulations [7], meaning that the maximal regularity is supported by the primal variable, namely the displacement vector field. More precisely, for any time $t \ge 0$, the solution of problem (1) is looked in the following functional spaces

$$
u(t) \in V = [H^1(\Omega)]^d \subset H, \quad \sigma(t) \in W = [L^2(\Omega)]_{\text{sym}}^{d \times d},
$$

where $H = [L^2(\Omega)]^d$, and $W$ is the space of square-integrable symmetric tensor fields. We define the scalar product in $H$ and $W$ as

$$
\left|
\begin{array}{ll}
(u, \tilde{u})_H = \displaystyle \int_{\Omega} \varrho \mathbf{u} \cdot \tilde{\mathbf{u}} \, \mathrm{d}\Omega, & \forall u, \tilde{u} \in H, \\
(\sigma, \tilde{\sigma})_W = \displaystyle \int_{\Omega} \boldsymbol{\sigma} : \tilde{\boldsymbol{\sigma}} \, \mathrm{d}\Omega, & \forall \sigma, \tilde{\sigma} \in W.
\end{array}
\right.
$$

Note that we differentiate the vector and tensor fields from their corresponding functions using bold font.
In our work we consider conform finite dimensional approximation spaces of the above defined functional
spaces, namely

$$
V_h \subset V, \quad W_h \subset W.
$$

We denote by $\mathcal{L}(W_h)$ the space of linear and bounded operators on $W_h$. For all symmetric strictly positive fourth-order tensor $\mathbf{S}$, we define the weighted mass operator $M_h(\mathbf{S}) \in \mathcal{L}(W_h)$ as

$$
(M_h(\mathbf{S})\sigma_h, \tilde{\sigma}_h)_W = \int_{\Omega} \mathbf{S}\sigma_h : \tilde{\sigma}_h \, d\Omega, \quad \forall \sigma_h, \tilde{\sigma}_h \in W_h.
$$

Note that, in practice, **S** typically represents a compliance tensor. $M_h(\mathbf{S})$ is a symmetric, strictly positive
finite dimensional operator, thus invertible, and we can define its square root as

$$
(M_h(\mathbf{S})^{\frac{1}{2}}\sigma_h, M_h(\mathbf{S})^{\frac{1}{2}}\tilde{\sigma}_h)_W = (M_h(\mathbf{S})\sigma_h, \tilde{\sigma}_h)_W, \quad \forall \sigma_h, \tilde{\sigma}_h \in W_h.
$$
We also define the divergence operator $B_h \in \mathcal{L}(V_h, W_h)$ as

$$ (B_h u_h, \sigma_h)_W = \int_{\Omega} \varepsilon(\mathbf{u}_h) : \sigma_h \, d\Omega, \quad \forall u_h \in V_h, \forall \sigma_h \in W_h, \qquad (6) $$

and we introduce $B_h^t \in \mathcal{L}(W_h, V_h)$ the (discrete) adjoint operator satisfying

$$ (B_h u_h, \sigma_h)_W = (u_h, B_h^t \sigma_h)_H. $$

For any symmetric strictly positive fourth-order tensor $C$, typically representing elastic modulii, we define the stiffness operator $A_h(C) \in \mathcal{L}(V_h)$ as

$$ A_h(C) = B_h^t M_h (C^{-1})^{-1} B_h. \qquad (7) $$

This definition is in fact a general form of a discrete stiffness operator, as proved in the following lemma.

**Lemma 1.** Let $C$ be a symmetric strictly positive fourth-order tensor, assuming that

$$ \forall u_h \in V_h, \quad \exists \gamma_h \in W_h, \quad \gamma_h = C\varepsilon(\mathbf{u}_h), \quad \text{a.e. } \Omega, \qquad (8) $$

then

$$ (A_h(C)u_h, \tilde{u}_h)_H = \int_{\Omega} C\varepsilon(\mathbf{u}_h) : \varepsilon(\tilde{u}_h) \, d\Omega, \quad \forall u_h, \tilde{u}_h \in V_h. $$

*Proof.* Given $u_h \in V_h$, we define $\sigma_h \in W$ such that $\sigma_h = C\varepsilon(\mathbf{u}_h)$ almost everywhere in $\Omega$. Upon assumption (8), $\sigma_h$ is an element of $W_h$ satisfying

$$ \int_{\Omega} C^{-1} \sigma_h : \tilde{\sigma}_h \, d\Omega = \int_{\Omega} \varepsilon(\mathbf{u}_h) : \tilde{\sigma}_h \, d\Omega, \quad \forall \tilde{\sigma}_h \in W_h, $$

or, equivalently, $M_h(C^{-1})\sigma_h = B_h u_h$. Thus, for any element $\tilde{u}_h \in V_h$, we have

$$ (B_h^t \sigma_h, \tilde{u}_h)_H = (B_h^t M_h(C^{-1})^{-1} B_h u_h, \tilde{u}_h)_H = (A_h(C)u_h, \tilde{u}_h)_H. $$

We conclude using the definition of the discrete adjoint divergence operator, namely

$$ (B_h^t \sigma_h, \tilde{u}_h)_H = (\sigma_h, B_h \tilde{u}_h)_W = \int_{\Omega} \sigma_h : \varepsilon(\tilde{u}_h) \, d\Omega = \int_{\Omega} C\varepsilon(\mathbf{u}_h) : \varepsilon(\tilde{u}_h) \, d\Omega. $$

□

**Remark 1.** We consider in our work continuous finite element approximation spaces based on parallelotopes. Namely, denoting by $\mathcal{T}_h$ a conform mesh of $\Omega$ without hanging nodes, and by $\mathbf{F}_K$ the transformation of the reference element $\tilde{K}$ (e.g. unit square in 2D or unit cube in 3D) into an element $K$, we define the following finite element space

$$ V_h = \{ v_h \in C^0(\bar{\Omega}, \mathbb{R}^d) \mid \forall K \in \mathcal{T}_h, \exists \hat{v} \in [\mathbb{Q}_k(\tilde{K})]^d, v_h|_K = \hat{v} \circ \mathbf{F}_K^{-1} \} \subset V, $$

where $k \ge 1$, and $\mathbb{Q}_k(\tilde{K})$ is the space of multivariate polynomials with maximal order smaller than $k$. For $W_h$ we simply consider the corresponding discontinuous finite element space, namely

$$ W_h = \{ \sigma_h \in W \mid \forall K \in \mathcal{T}_h, \exists \hat{\sigma} \in [\mathbb{Q}_k(\tilde{K})]_{\text{sym}}^{d \times d}, \sigma_h|_K = \hat{\sigma} \circ \mathbf{F}_K^{-1} \} \subset W. $$

With these definitions of $V_h$ and $W_h$, assumption (8) is satisfied e.g. when each $\mathbf{F}_K$ is an affine transformation of the reference element, and $C$ is piece-wise constant.
The definition (7) of the stiffness operator does not imply, in general, linearity w.r.t the tensor $C$. It is however expected from such operator to bear linear dependency, i.e.

$$A_h(C_0 + C_1) = A_h(C_0) + A_h(C_1), \qquad (9)$$

which essentially boils down to the following assumption on the weighted mass operator

$$M_h((C_0 + C_1)^{-1})^{-1} = M_h(C_0^{-1})^{-1} + M_h(C_1^{-1})^{-1}, \qquad (10)$$

for any symmetric strictly positive fourth-order tensor fields $C_0$ and $C_1$. Note that, if (8) holds for every $C$ in a given vectorial space $C$ of symmetric strictly positive fourth-order tensor fields (e.g. piece-wise constant tensor fields), the linearity is obviously true in that space – in virtue of Lemma 1. We give in the next remark another practical case where this property is satisfied.

**Remark 2.** In addition to these conform space discretizations, we assume that a suitable quadrature formula is used to approximate the integrals appearing in the expression of the inner product $(\cdot, \cdot)_H$ and the weighted mass operators $M_h(\cdot)$, rendering diagonal mass matrices. This approach is the so-called mass lumping technique, available for either parallelotopes [2, 8, 4, 21] or simplicial elements [22, 3, 23]. It is of particular interest for high-frequency wave propagation problems in time domain since it can lead to fully explicit scheme after time discretization. Since the mass-lumping procedure entails diagonal mass operators, the assumption (10) is necessarily satisfied, entailing linearity of the discrete stiffness operators w.r.t. the elastic modulii.

### 2.3. Semi-discrete formulations of visco-elastic models

With these definitions of the discrete operators, the space discretizations of the visco-elastic wave propagation problems read

$$f = d_{tt}^2 u_h + \begin{cases} B_h^t \sigma_h, & \text{for MX,} \\ A_h(C_{zn})u_h + B_h^t s_h, & \text{for ZN,} \\ A_h(C_{kv})u_h + A_h(D_{kv})d_t u_h, & \text{for KV.} \end{cases} \qquad (11)$$

For the MX and ZN models, this relation is completed with the space discretization of their respective constitutive laws, namely

$$\left| \begin{array}{ll} M_h(\mathbf{S}_{\text{mx}})d_t \sigma_h + M_h(\mathbf{H}_{\text{mx}})\sigma_h = B_h d_t u_h, & \text{for MX,} \\ M_h(\mathbf{S}_{\text{zn}})d_t s_h + M_h(\mathbf{H}_{\text{zn}})s_h = B_h d_t u_h & \text{for ZN.} \end{array} \right. \qquad (12)$$

**Remark 3.** In the following, in order to simplify the presentation, we will denote by

$$S_{\text{mx},h} = M_h(\mathbf{S}_{\text{mx}}), \quad H_{\text{mx},h} = M_h(\mathbf{H}_{\text{mx}}), \quad \text{and} \quad S_{\text{zn},h} = M_h(\mathbf{S}_{\text{zn}}), \quad H_{\text{zn},h} = M_h(\mathbf{H}_{\text{zn}}),$$

the weighted mass operators.

### 2.4. Fully-discrete time explicit discretizations

Explicit time discretization with mass lumping is a well-known combination since it leads to particularly efficient time marching algorithm. Indeed, thanks to the lumped mass, no linear system needs to be solved at each time step. Furthermore, the stability condition – or CFL condition – of standard explicit scheme for hyperbolic problem is usually acceptable, taking the form of $\Delta t = \mathcal{O}(h/c_{\max})$, with $c_{\max}$ the maximal wave velocity in the material.

Since they are the starting point of our work, we review in this section the most simple explicit schemes for the space-discrete MX, ZN and KV models, and give their associated stability conditions. Readers may
refer to [6] for more details. Let us denote by $u_h^n$ the time approximation of $u_h(n\Delta t)$. Using the same notation for the stress unknowns and the source term, the explicit time discretization of (11) is

$$f^n = \frac{u_h^{n+1} - 2u_h^n + u_h^{n-1}}{\Delta t^2} + \begin{cases} B_h^t \sigma_h^n, & \text{for MX,} \\ A_h(C_{\text{zn}})u_h^n + B_h^t s_h^n, & \text{for ZN,} \\ A_h(C_{\text{kv}})u_h^n + A_h(D_{\text{kv}}) \frac{u_h^n - u_h^{n-1}}{\Delta t}, & \text{for KV.} \end{cases} \quad (13)$$

Particular emphasis should be given to the uncentered discretization of the first order time derivative in the KV case, in order to keep a fully explicit scheme, since a stiffness operator cannot be lumped. This operation however effectively reduces the approximation order of the time discretization procedure for the KV case. For the MX and ZN case, the constitutive relations (12) are discretized using an implicit scheme centered at time $t^{n+1/2} = (n + 1/2)\Delta t$, i.e.

$$\left| \begin{aligned} & S_{\text{mx},h} \frac{\sigma_h^{n+1} - \sigma_h^n}{\Delta t} + H_{\text{mx},h} \frac{\sigma_h^{n+1} + \sigma_h^n}{2} = B_h \frac{u_h^{n+1} - u_h^n}{\Delta t}, && \text{for MX,} \\ & S_{\text{zn},h} \frac{s_h^{n+1} - s_h^n}{\Delta t} + H_{\text{zn},h} \frac{s_h^{n+1} + s_h^n}{2} = B_h \frac{u_h^{n+1} - u_h^n}{\Delta t}, && \text{for ZN,} \end{aligned} \right.$$

without breaking the explicit nature of the final algorithm since the displacement unknown can be computed in a prior step, and the weighted mass operators $S_{\text{mx},h}$, $H_{\text{mx},h}$, $S_{\text{zn},h}$ and $H_{\text{zn},h}$ are lumped – see Remark 2. For any symmetric strictly positive fourth-order tensor $C$, we denote the spectral radius of the associated discrete stiffness operator by

$$\mathcal{R}_C = \sup_{\{v_h \in V_h, v_h \neq 0\}} \frac{(A_h(C)v_h, v_h)_H}{\|v_h\|_H^2},$$

so that the stability conditions for the explicit time schemes are, see [6],

$$\Delta t \le 2\mathcal{R}_{C_{\text{mx}}}^{-\frac{1}{2}}, \quad \Delta t \le 2\mathcal{R}_{D_{\text{zn}}}^{-\frac{1}{2}}, \quad \Delta t \le 2\mathcal{R}_{C_{\text{kv}}}^{-\frac{1}{2}}\left(\left(\frac{\mathcal{R}_{D_{\text{kv}}}^2}{4\mathcal{R}_{C_{\text{kv}}}} + 1\right)^{\frac{1}{2}} - \frac{\mathcal{R}_{D_{\text{kv}}}}{2\sqrt{\mathcal{R}_{C_{\text{kv}}}}}\right), \quad (14)$$

for MX, ZN and KV respectively.

## 2.5. Positioning of the problem

In our work, we consider cases where the propagation domain can be decomposed into two parts: $\Omega_0$, which may be referred to as the background domain, and $\Omega_*$, a perturbation domain. Naturally, the integral in the definition of the divergence operator (6) can be decomposed as a sum of integrals over both domains. Hence, we define the operators

$$(B_{h,0}u_h, \sigma_h)_W = \int_{\Omega_0} \varepsilon(\mathbf{u}_h) : \sigma_h \, d\Omega, \quad (B_{h,\star}u_h, \sigma_h)_W = \int_{\Omega_*} \varepsilon(\mathbf{u}_h) : \sigma_h \, d\Omega, \quad (15)$$

for any $u_h \in V_h$, and $\sigma_h \in W_h$, so that $B_h = B_{h,0} + B_{h,\star}$. Furthermore, we assume that the conform approximation $W_h$ of the dual space $W$ can be written as the direct sum of two sub-spaces

$$W_h = W_{h,0} \oplus W_{h,\star}, \quad \left| \begin{array}{ll} \forall \sigma_{h,0} \in W_{h,0}, & \text{supp } \sigma_{h,0} \subset \bar{\Omega}_0, \\ \forall \sigma_{h,\star} \in W_{h,\star}, & \text{supp } \sigma_{h,\star} \subset \bar{\Omega}_*. \end{array} \right. \quad (16)$$

**Property 1.** Upon the assumption (16), and the definition (15) of the operators $B_{h,0}$ and $B_{h,\star}$, we have
i) For any $\sigma_h = \sigma_{h,0} + \sigma_{h,\star} \in W_h$, with $\sigma_{h,0} \in W_{h,0}$ and $\sigma_{h,\star} \in W_{h,\star}$,

$$B_h^t \sigma_h = B_{h,0}^t \sigma_{h,0} + B_{h,\star}^t \sigma_{h,\star}. \qquad (17)$$

ii) For any symmetric strictly positive fourth-order tensor field $C$, defining

$$A_{h,0}(C) = B_{h,0}^t M_h(C^{-1})^{-1} B_{h,0}, \quad A_{h,\star}(C) = B_{h,\star}^t M_h(C^{-1})^{-1} B_{h,\star}, \qquad (18)$$

then the following decomposition holds

$$A_h(C) = A_{h,0}(C) + A_{h,\star}(C).$$

**Proof.** i) From the definition (15), and using the fact that $\text{supp } \sigma_{h,\star} \subset \bar{\Omega}_\star$, we have, for any $u_h \in V_h$,

$$(B_{h,0}u_h, \sigma_h)_W = \int_{\Omega_0} \varepsilon(\mathbf{u}_h) : (\boldsymbol{\sigma}_{h,0} + \boldsymbol{\sigma}_{h,\star}) \, d\Omega = \int_{\Omega_0} \varepsilon(\mathbf{u}_h) : \boldsymbol{\sigma}_{h,0} \, d\Omega = (B_{h,0}^t \sigma_{h,0}, u_h)_H.$$

Proceeding equivalently for $B_{h,\star}$ leads to

$$(B_h^t \sigma_h, u_h)_H = (\sigma_h, B_{h,0}u_h)_W + (\sigma_h, B_{h,\star}u_h)_W = (B_{h,0}^t \sigma_{h,0} + B_{h,\star}^t \sigma_{h,\star}, u_h)_H,$$

ii) Remarking that $\text{Im } B_{h,0} \subset (W_{h,\star})^{\perp} = W_{h,0}$, and reciprocally that $\text{Im } B_{h,\star} \subset W_{h,\star}$, leads to

$$\text{Im } B_{h,0} \subset (W_{h,\star})^{\perp} \subset (\text{Im } B_{h,\star})^{\perp} = \text{Ker } B_{h,\star}^t.$$

Let $u_h \in V_h$, and consider $\sigma_{h,0} = B_{h,0}u_h \in \text{Im } B_{h,0}$, and $\mu_h = M_h(C^{-1})^{-1}\sigma_{h,0} \in W_h$. We have

$$\forall \sigma'_{h,\star} \in W_{h,\star}, \quad \left\| \begin{array}{l} (\sigma_{h,0}, \sigma'_{h,\star})_W = (M_h(C^{-1})\mu_h, \sigma'_{h,\star})_W = \int_{\Omega} C^{-1} \mu_h : \sigma'_{h,\star} \, d\Omega, \\ (\sigma_{h,0}, \sigma'_{h,\star})_W = (B_{h,0}u_h, \sigma'_{h,\star})_W = \int_{\Omega_0} \varepsilon(\mathbf{u}_h) : \sigma'_{h,\star} \, d\Omega = 0, \end{array} \right.$$

Combined with the bijectivity of $C$, it implies that $\mu_h \in (W_{h,\star})^{\perp}$. In other words,

$$\text{Im } (M_h(C^{-1})^{-1} B_{h,0}) \subset (W_{h,\star})^{\perp} \subset \text{Ker } B_{h,\star}^t,$$

i.e. $B_{h,\star}^t M_h(C^{-1})^{-1} B_{h,0} = 0$. We conclude using (7), and $B_h = B_{h,0} + B_{h,\star}$.

□

**Remark 4.** Since, in practice, we consider discontinuous approximation functions to build $W_h$ – see Remark 1 – the assumption (16) is satisfied if the mesh $T_h$ can be decomposed into two conform sub-meshes of $\Omega_0$ and $\Omega_\star$.

In the context of our work, we assume that the stiffness operator within $\Omega_\star$ – namely $A_{h,\star}$ – has a significantly greater spectral radius than the stiffness operator within $\Omega_0$ – namely $A_{h,0}$. More precisely, denoting the suprema of the Rayleigh quotients by

$$\mathcal{R}_{C,0} = \sup_{\{v_h \in V_h, v_h \neq 0\}} \frac{(A_{h,0}(C)v_h, v_h)_H}{\|v_h\|_H^2}, \quad \mathcal{R}_{C,\star} = \sup_{\{v_h \in V_h, v_h \neq 0\}} \frac{(A_{h,\star}(C)v_h, v_h)_H}{\|v_h\|_H^2}, \qquad (19)$$

we assume that

$$\mathcal{R}_{C,0} \ll \mathcal{R}_{C,\star}.$$

This is typically the case when the elastic modulii (thus the associated wave velocities) are significantly greater within the perturbation zone $\Omega_\star$ compared to the background medium or when the mesh is much finer or composed of distorted elements. In any case, the stability conditions (14) and the relation $\mathcal{R}_{C,\star} \le \mathcal{R}_C$, imply that the perturbation zone has a potentially dramatic effect on the time step of the explicit time schemes. For a fixed global time window, this can lead to a substantial increase of the overall computation cost. To circumvent this difficulty we propose two strategies:
* We apply an implicit scheme exclusively on the operator penalizing the stability conditions. This is in essence a locally implicit approach, and can be understood as an extension for visco-elastic models of the prior works of [11, 13, 12]. By doing so, we recover stability conditions independent of the perturbation zone, at a cost of solving (localized) linear systems at each time step.

* We extend to visco-elastic problems the stabilization strategy proposed in [14]. This approach consists in applying the mortar element method [16, 17, 18] to separate the background and perturbed domains, and applying a stabilized leapfrog scheme, in the vein of [15, 24], within $\Omega_{\star}$. This procedure is explicit, and renders a CFL condition which can be increased by a factor proportional to the order of the polynomial used in the stabilized scheme. Compared to other local time-stepping approaches [25, 24] the main difference resides in the use of a domain decomposition procedure, which provides nonetheless an efficient framework for insuring stability using energy arguments at the cost of additional multiplications by stiffness operators in $\Omega_{\star}$ and solving a linear system at the interfaces between the subdomains.

### 3. Locally implicit time schemes for visco-elastic wave problems

In this section, we derive locally implicit schemes for each visco-elastic models, and we analyze their stability. To ensure stability, we need to show that there exists a uniform estimate in a space-time norm (typically $L^{\infty}$ in time and $L^2$ or $H^1$ in space) of the discrete solution w.r.t the initial conditions and/or source terms. In [3, 7, 26], such estimates for classes of wave problems are provided using so-called energy techniques. In these works, authors show that the conservation of a positive discrete energy functional is a sufficient condition to obtain these stability estimates. Therefore, following this line of work, we assume in the rest of this paper that conservation or decay of a positive energy functional, potentially upon a CFL condition, is sufficient to ensure stability of the discrete scheme.

#### 3.1. Preliminary identities for time-discrete constitutive relations

Let us start by establishing some simple but relatively technical identities that are used in the next sections to obtain energy relations. To do so we consider $S$ and $H$ two symmetric strictly positive fourth-order tensor fields.

Time-discrete constitutive relation centered at $t^{n+\frac{1}{2}}$. Let $\{u_h^n\}_{n \ge 0} \subset V_h$ and $\{\sigma_h^n\}_{n \ge 0} \subset W_h$ satisfy the following time-discrete constitutive relation

$$M_h(S) \frac{\sigma_h^{n+1} - \sigma_h^n}{\Delta t} + M_h(H) \frac{\sigma_h^{n+1} + \sigma_h^n}{2} = B_h \frac{u_h^{n+1} - u_h^n}{\Delta t}. \quad (20)$$

**Property 2.** Associated to the stress and displacement unknown satisfying (20), we define the following functionals

$$\left\{
\begin{aligned}
\mathcal{J}^{n+\frac{1}{2}} &= \frac{1}{4} \left\| M_h(S)^{\frac{1}{2}} \sigma_h^{n+1} \right\|_W^2 + \frac{1}{4} \left\| M_h(S)^{\frac{1}{2}} \sigma_h^n \right\|_W^2, \\
\mathcal{B}^{n+\frac{1}{2}} &= \left( \frac{\sigma_h^{n+1} - \sigma_h^n}{\Delta t}, B_h \frac{u_h^{n+1} - u_h^n}{\Delta t} \right)_W, \\
\mathcal{D}^n &= \frac{1}{2} \left\| M_h(H)^{\frac{1}{2}} \frac{\sigma_h^{n+1} + \sigma_h^n}{2} \right\|_W^2 + \frac{1}{2} \left\| M_h(H)^{\frac{1}{2}} \frac{\sigma_h^n + \sigma_h^{n-1}}{2} \right\|_W^2.
\end{aligned}
\right.
\quad (21)$$

These functionals satisfy the following relation

$$\frac{\mathcal{J}^{n+\frac{1}{2}} - \mathcal{J}^{n-\frac{1}{2}}}{\Delta t} - \frac{\Delta t^2}{4} \frac{\mathcal{B}^{n+\frac{1}{2}} - \mathcal{B}^{n-\frac{1}{2}}}{\Delta t} + \mathcal{D}^n = \left( \sigma_h^n, B_h \frac{u_h^{n+1} - u_h^{n-1}}{2\Delta t} \right)_W.$$
*Proof.* Taking the scalar product in W of the equation (20) with $\frac{1}{2}(\sigma_h^{n+1} + \sigma_h^n)$ and then computing the average of the result at times $t^{n+\frac{1}{2}}$ and $t^{n-\frac{1}{2}}$ yields the following

$$ \frac{1}{4\Delta t} \left( \left\| M_h(S)^{\frac{1}{2}} \sigma_h^{n+1} \right\|_W^2 - \left\| M_h(S)^{\frac{1}{2}} \sigma_h^{n-1} \right\|_W^2 \right) + \frac{1}{2} \left\| M_h(H)^{\frac{1}{2}} \frac{\sigma_h^{n+1} + \sigma_h^n}{2} \right\|_W^2 + \frac{1}{2} \left\| M_h(H)^{\frac{1}{2}} \frac{\sigma_h^n + \sigma_h^{n-1}}{2} \right\|_W^2 \\ - \frac{1}{2} \left( \frac{\sigma_h^{n+1} + \sigma_h^n}{2}, B_h \frac{u_h^{n+1} - u_h^n}{\Delta t} \right)_W - \frac{1}{2} \left( \frac{\sigma_h^n + \sigma_h^{n-1}}{2}, B_h \frac{u_h^n - u_h^{n-1}}{\Delta t} \right)_W = 0. \quad (22) $$

Remarking that

$$ -\frac{1}{2} \left( \frac{\sigma_h^{n+1} + \sigma_h^n}{2}, B_h \frac{u_h^{n+1} - u_h^n}{\Delta t} \right)_W - \frac{1}{2} \left( \frac{\sigma_h^n + \sigma_h^{n-1}}{2}, B_h \frac{u_h^n - u_h^{n-1}}{\Delta t} \right)_W \\ = - \left( \frac{\sigma_h^{n+1} - \sigma_h^n}{2}, B_h \frac{u_h^{n+1} - u_h^n}{2\Delta t} \right)_W + \left( \frac{\sigma_h^n - \sigma_h^{n-1}}{2}, B_h \frac{u_h^n - u_h^{n-1}}{2\Delta t} \right)_W - \left( \sigma_h^n, B_h \frac{u_h^{n+1} - u_h^{n-1}}{2\Delta t} \right)_W. $$

leads to the desired result.

<figure><img src="image_3.png" /></figure>

**Property 3.** Denoting by $C = S^{-1}$, the first two functionals defined in (21) satisfy

$$ \mathcal{P}^{n+\frac{1}{2}} - \frac{\Delta t^2}{4} \mathcal{B}^{n+\frac{1}{2}} = -\frac{\Delta t^2}{8} \left( A_h(C) \frac{u_h^{n+1} - u_h^n}{\Delta t}, \frac{u_h^{n+1} - u_h^n}{\Delta t} \right)_H + \frac{1}{2} \left\| M_h(S)^{\frac{1}{2}} \frac{\sigma_h^{n+1} + \sigma_h^n}{2} \right\|_W^2 \\ + \frac{\Delta t^2}{8} \left\| M_h(S)^{\frac{1}{2}} \frac{\sigma_h^{n+1} - \sigma_h^n}{\Delta t} - M_h(S)^{-\frac{1}{2}} B_h \frac{u_h^{n+1} - u_h^n}{\Delta t} \right\|_W^2 $$

*Proof.* We can express $\mathcal{P}^{n+\frac{1}{2}}$ as

$$ \mathcal{P}^{n+\frac{1}{2}} = \frac{1}{2} \left\| M_h(S)^{\frac{1}{2}} \frac{\sigma_h^{n+1} + \sigma_h^n}{2} \right\|_W^2 + \frac{\Delta t^2}{8} \left\| M_h(S)^{\frac{1}{2}} \frac{\sigma_h^{n+1} - \sigma_h^n}{\Delta t} \right\|_W^2, $$

and see that, using the decomposition (7) of the stiffness operator, we have

$$ \frac{\Delta t^2}{8} \left\| M_h(S)^{\frac{1}{2}} \frac{\sigma_h^{n+1} - \sigma_h^n}{\Delta t} \right\|_W^2 - \frac{\Delta t^2}{4} \mathcal{B}^{n+\frac{1}{2}} = -\frac{\Delta t^2}{8} \left( A_h(C) \frac{u_h^{n+1} - u_h^n}{\Delta t}, \frac{u_h^{n+1} - u_h^n}{\Delta t} \right)_H \\ + \frac{\Delta t^2}{8} \left( \left\| M_h(S)^{\frac{1}{2}} \frac{\sigma_h^{n+1} - \sigma_h^n}{\Delta t} \right\|_W^2 - 2 \left( \frac{\sigma_h^{n+1} - \sigma_h^n}{\Delta t}, B_h \frac{u_h^{n+1} - u_h^n}{\Delta t} \right)_W + \left( A_h(C) \frac{u_h^{n+1} - u_h^n}{\Delta t}, \frac{u_h^{n+1} - u_h^n}{\Delta t} \right)_H \right), $$

from which we recognize in the last three terms the square of a binomial, yielding the expected result.

<figure><img src="image_4.png" /></figure>

*Time-discrete constitutive relation centered at t^n.* Let us now consider the following time-discrete constitutive relation

$$ M_h(S) \frac{\sigma_h^{n+\frac{1}{2}} - \sigma_h^{n-\frac{1}{2}}}{\Delta t} + M_h(H) \frac{\sigma_h^{n+\frac{1}{2}} + \sigma_h^{n-\frac{1}{2}}}{2} = B_h \frac{u_h^{n+1} - u_h^{n-1}}{2\Delta t}, \quad (23) $$

where $\{u_h^n\}_{n \ge 0} \subset V_h$ and $\{(\sigma_h^{n+1})\}_{n \ge 0} \subset W_h$.

**Property 4.** Associated to the stress and displacement unknowns satisfying (23), we define the following functionals

$$ \mathcal{P}^{n+\frac{1}{2}} = \frac{1}{2} \left\| M_h(S)^{\frac{1}{2}} \sigma_h^{n+\frac{1}{2}} \right\|_W^2, \quad \mathcal{D}^n = \left\| M_h(H)^{\frac{1}{2}} \frac{\sigma_h^{n+\frac{1}{2}} + \sigma_h^{n-\frac{1}{2}}}{2} \right\|_W^2. $$

These functionals satisfy the following relation

$$ \frac{\mathcal{P}^{n+\frac{1}{2}} - \mathcal{P}^{n-\frac{1}{2}}}{\Delta t} + \mathcal{D}^n = \left( \frac{\sigma_h^{n+\frac{1}{2}} + \sigma_h^{n-\frac{1}{2}}}{2}, B_h \frac{u_h^{n+1} - u_h^{n-1}}{2\Delta t} \right)_W. $$
Proof. Taking the scalar product of (23) with $(\sigma_h^{n+\frac{1}{2}} + \sigma_h^{n-\frac{1}{2}})/2$ leads to

$$ \frac{1}{2\Delta t} \left( \left\| M_h(S)^{\frac{1}{2}} \sigma_h^{n+\frac{1}{2}} \right\|_W^2 - \left\| M_h(S)^{\frac{1}{2}} \sigma_h^{n-\frac{1}{2}} \right\|_W^2 \right) + \mathcal{D}^n = \left( \frac{\sigma_h^{n+\frac{1}{2}} + \sigma_h^{n-\frac{1}{2}}}{2}, B_h \frac{u_h^{n+1} - u_h^{n-1}}{2\Delta t} \right)_W. $$

□

### 3.2. Maxwell model

Upon on the assumption (16) and the property (17) of the divergence operator, we rely on a splitting of the dual variable. For any time $t > 0$, we consider

$$ \sigma_{h,0}(t) \in W_{h,0} \quad \text{and} \quad \sigma_{h,\star}(t) \in W_{h,\star} $$

and $u_h(t) \in V_h$ the solution of

$$ \begin{cases} d_{tt}^2 u_h + B_{h,0}^t \sigma_{h,0} + B_{h,\star}^t \sigma_{h,\star} = f, \\ S_{\mathrm{mx},h} d_t \sigma_{h,0} + H_{\mathrm{mx},h} \sigma_{h,0} = B_{h,0} d_t u_h, \\ S_{\mathrm{mx},h} d_t \sigma_{h,\star} + H_{\mathrm{mx},h} \sigma_{h,\star} = B_{h,\star} d_t u_h. \end{cases} $$

One can check, by linearity and in virtue of (17), that $\sigma_h = \sigma_{h,0} + \sigma_{h,\star}$. Upon this formulation, the locally implicit time scheme for the Maxwell model reads

$$ \left\{ \begin{aligned} & \frac{u_h^{n+1} - 2u_h^n + u_h^{n-1}}{\Delta t^2} + B_{h,0}^t \sigma_{h,0} + B_{h,\star}^t \frac{\sigma_{h,\star}^{n+\frac{1}{2}} + \sigma_{h,\star}^{n-\frac{1}{2}}}{2} = f^n, \\ & S_{\mathrm{mx},h} \frac{\sigma_{h,0}^{n+1} - \sigma_{h,0}^n}{\Delta t} + H_{\mathrm{mx},h} \frac{\sigma_{h,0}^{n+1} + \sigma_{h,0}^n}{2} = B_{h,0} \frac{u_h^{n+1} - u_h^n}{\Delta t}, \\ & S_{\mathrm{mx},h} \frac{\sigma_{h,\star}^{n+\frac{1}{2}} - \sigma_{h,\star}^{n-\frac{1}{2}}}{\Delta t} + H_{\mathrm{mx},h} \frac{\sigma_{h,\star}^{n+\frac{1}{2}} + \sigma_{h,\star}^{n-\frac{1}{2}}}{2} = B_{h,\star} \frac{u_h^{n+1} - u_h^{n-1}}{2\Delta t}. \end{aligned} \right. \quad (24) $$

Remark that in the first equation, the stress unknown in the perturbed domain $\Omega_*$ is treated using an implicit scheme. Additionally, in the sound domain $\Omega_0$ we use a time-discrete constitutive relation of the form of (20), while in the perturbed domain $\Omega_*$ we use (23).

**Lemma 2.** Associated to the fully-discrete solution of (24), let us define the following

* Kinetic energy functional:

$$ \mathcal{K}^{n+\frac{1}{2}} = \frac{1}{2} \left\| \frac{u_h^{n+1} - u_h^n}{\Delta t} \right\|_H^2, $$

* Potential energy functionals:

$$ \mathcal{P}_0^{n+\frac{1}{2}} = \frac{1}{4} \left\| S_{\mathrm{mx},h}^{\frac{1}{2}} \sigma_{h,0}^{n+1} \right\|_W^2 + \frac{1}{4} \left\| S_{\mathrm{mx},h}^{\frac{1}{2}} \sigma_{h,0}^n \right\|_W^2, \quad \mathcal{P}_*^{n+\frac{1}{2}} = \frac{1}{2} \left\| S_{\mathrm{mx},h}^{\frac{1}{2}} \sigma_{h,\star}^{n+\frac{1}{2}} \right\|_W^2, $$

* Cross-term functional:

$$ \mathcal{B}^{n+\frac{1}{2}} = \left( \frac{\sigma_{h,0}^{n+1} - \sigma_{h,0}^n}{\Delta t}, B_{h,0} \frac{u_h^{n+1} - u_h^n}{\Delta t} \right)_W, $$

* Energy dissipation functional:

$$ \mathcal{D}^n = \frac{1}{2} \left\| H_{\mathrm{mx},h}^{\frac{1}{2}} \frac{\sigma_{h,0}^{n+1} + \sigma_{h,0}^n}{2} \right\|_W^2 + \frac{1}{2} \left\| H_{\mathrm{mx},h}^{\frac{1}{2}} \frac{\sigma_{h,0}^n + \sigma_{h,0}^{n-1}}{2} \right\|_W^2 + \left\| H_{\mathrm{mx},h}^{\frac{1}{2}} \frac{\sigma_{h,\star}^{n+\frac{1}{2}} + \sigma_{h,\star}^{n-\frac{1}{2}}}{2} \right\|_W^2, $$
• Total energy functional:

$$ \varepsilon^{n+\frac{1}{2}} = \mathcal{K}^{n+\frac{1}{2}} + \mathcal{P}_{\star}^{n+\frac{1}{2}} + \mathcal{P}_{0}^{n+\frac{1}{2}} - \frac{\Delta t^2}{4} \mathcal{B}^{n+\frac{1}{2}}, $$

then, upon the assumption of null source, i.e. $f^n = 0$, we have the following decay property,

$$ \frac{\varepsilon^{n+\frac{1}{2}} - \varepsilon^{n-\frac{1}{2}}}{\Delta t} = -\mathcal{D}^n \leq 0. $$

*Proof.* By taking the scalar product of the first equation in (24) with $\frac{1}{2\Delta t}(u_h^{n+1} - u_h^{n-1})$, we obtain a first relation:

$$ \frac{\mathcal{K}^{n+\frac{1}{2}} - \mathcal{K}^{n-\frac{1}{2}}}{\Delta t} + \left( B_{h,0} \frac{u_h^{n+1} - u_h^{n-1}}{2\Delta t}, \sigma_{h,0}^n \right)_W + \left( B_{h,\star} \frac{u_h^{n+1} - u_h^{n-1}}{2\Delta t}, \frac{\sigma_{h,\star}^{n+\frac{1}{2}} + \sigma_{h,\star}^{n-\frac{1}{2}}}{2} \right)_W = 0. $$

Using Property 2 and Property 4 for the second and third equation in (24) respectively, and combining the two obtained expressions with the one above leads to

$$ \frac{\mathcal{K}^{n+\frac{1}{2}} - \mathcal{K}^{n-\frac{1}{2}}}{\Delta t} + \frac{\mathcal{P}_{\star}^{n+\frac{1}{2}} - \mathcal{P}_{\star}^{n-\frac{1}{2}}}{\Delta t} + \frac{\mathcal{P}_{0}^{n+\frac{1}{2}} - \mathcal{P}_{0}^{n-\frac{1}{2}}}{\Delta t} - \frac{\Delta t^2}{4} \frac{\mathcal{B}^{n+\frac{1}{2}} - \mathcal{B}^{n-\frac{1}{2}}}{\Delta t} = -\mathcal{D}^n \leq 0. $$

□

**Lemma 3.** The total energy functional associated to the time discrete scheme (24) is positive on the following condition

$$ \Delta t \leq 2\mathcal{R}_{C_{\mathrm{mx},0}}^{-\frac{1}{2}}. \qquad (25) $$

where $C_{\mathrm{mx}}$ is defined in (3), and $\mathcal{R}_{C_{\mathrm{mx},0}}$ in (19).

*Proof.* Injecting into the energy functional the expression obtained when using Property 3 to $\mathcal{P}_{0}^{n+\frac{1}{2}} - \frac{\Delta t^2}{4} \mathcal{B}^{n+\frac{1}{2}}$ gives the following

$$ \varepsilon^{n+\frac{1}{2}} = \tilde{\varepsilon}^{n+\frac{1}{2}} + \frac{1}{2} \left( \left( I_h - \frac{\Delta t^2}{4} A_{h,0}(C_{\mathrm{mx}}) \right) \frac{u_h^{n+1} - u_h^n}{\Delta t}, \frac{u_h^{n+1} - u_h^n}{\Delta t} \right)_H, $$

where $A_{h,0}$ is defined in (18), and

$$ \tilde{\varepsilon}^{n+\frac{1}{2}} = \frac{1}{2} \left\| S_{\mathrm{mx},h}^{\frac{1}{2}} \frac{\sigma_{h,0}^{n+1} + \sigma_{h,0}^{n}}{2} \right\|_{W}^{2} + \frac{1}{2} \left\| S_{\mathrm{mx},h}^{\frac{1}{2}} \sigma_{h,\star}^{n+\frac{1}{2}} \right\|_{W}^{2} + \frac{\Delta t^2}{8} \left\| S_{\mathrm{mx},h}^{\frac{1}{2}} \frac{\sigma_{h,0}^{n+1} - \sigma_{h,0}^{n}}{\Delta t} - S_{\mathrm{mx},h}^{-\frac{1}{2}} B_h \frac{u_h^{n+1} - u_h^n}{\Delta t} \right\|_{W}^{2}. $$

is a positive functional. Hence, the total energy functional is positive if $I_h - \frac{\Delta t^2}{4} A_{h,0}$ is positive definite, equivalently if $\Delta t \leq 2\mathcal{R}_{C_{\mathrm{mx},0}}^{-\frac{1}{2}}$. □

**Remark 5.** Owing to the fact that

$$ \forall v_h \in V_h, \quad \|v_h\|_H^2 \geq (\rho v_h, v_h)_{L^2(\Omega_0)^d}, $$

a sufficient condition for the positivity of the total energy functional, hence the stability of the scheme, is

$$ \Delta t \leq 2 \left( \sup_{\{v_h \in V_h, v_h \neq 0\}} \frac{(A_{h,0}(C_{\mathrm{mx}})v_h, v_h)_H}{(\rho v_h, v_h)_{L^2(\Omega_0)^d}} \right)^{-\frac{1}{2}}, $$

which is independent of the material parameters $C_{\mathrm{mx}}$ and $\rho$ in $\Omega_\star$, or the mesh elements occupying $\Omega_\star$.
### 3.3. Zener model

To express the locally implicit schemes for the Zener model, we introduce the so-called $\theta$-schemes, where $\theta \in \mathbb{R}$. Namely, for any $\{u_h^n\}_{n \ge 0} \subset V_h$, we define

$$ \{u_h^n\}_\theta = \theta u_h^{n+1} + (1-2\theta)u_h^n + \theta u_h^{n-1} \quad (26) $$

To simplify the presentation, we denote by $A_h^e = A_h(C_{zn})$ the purely elastic stiffness operator. With these notations, and following the same decomposition strategy as for the Maxwell case, the fully-discrete locally implicit scheme for the ZN model defined in (11) reads

$$ \left\{ \begin{array}{ll} \displaystyle \frac{u_h^{n+1} - 2u_h^n + u_h^{n-1}}{\Delta t^2} + A_{h,0}^e u_h^n + A_{h,\star}^e \{u_h^n\}_\theta + B_{h,0}^t s_{h,0}^n + B_{h,\star}^t \frac{s_{h,\star}^{n+1/2} + s_{h,\star}^{n-1/2}}{2} & = f^n, \\ \\ S_{zn,h} \frac{s_{h,0}^{n+1} - s_{h,0}^n}{\Delta t} + H_{zn,h} \frac{s_{h,0}^{n+1} + s_{h,0}^n}{2} & = B_{h,0} \frac{u_h^{n+1} - u_h^n}{\Delta t}, \\ \\ S_{zn,h} \frac{s_{h,\star}^{n+1/2} - s_{h,\star}^{n-1/2}}{\Delta t} + H_{zn,h} \frac{s_{h,\star}^{n+1/2} + s_{h,\star}^{n-1/2}}{2} & = B_{h,\star} \frac{u_h^{n+1} - u_h^{n-1}}{2\Delta t}. \end{array} \right. \quad (27) $$

**Lemma 4.** Associated to the fully-discrete solution of (27), let us define the following

* Kinetic energy functional:

$$ \mathcal{K}^{n+\frac{1}{2}} = \frac{1}{2} \left( \left( I_h - \frac{\Delta t^2}{4} A_{h,0}^e + \Delta t^2 \left(\theta - \frac{1}{4}\right) A_{h,\star}^e \right) \frac{u_h^{n+1} - u_h^n}{\Delta t}, \frac{u_h^{n+1} - u_h^n}{\Delta t} \right)_H, $$

* Potential energy functionals:

$$ \left\{ \begin{aligned} \mathcal{P}_e^{n+\frac{1}{2}} &= \frac{1}{2} \left( A_h^e \frac{u_h^{n+1} + u_h^n}{2}, \frac{u_h^{n+1} + u_h^n}{2} \right)_H, \\ \mathcal{P}_0^{n+\frac{1}{2}} &= \frac{1}{4} \left\| S_{zn,h}^{\frac{1}{2}} s_{h,0}^{n+1} \right\|_W^2 + \frac{1}{4} \left\| S_{zn,h}^{\frac{1}{2}} s_{h,0}^n \right\|_W^2, \\ \mathcal{P}_\star^{n+\frac{1}{2}} &= \frac{1}{2} \left\| S_{zn,h}^{\frac{1}{2}} s_{h,\star}^{n+\frac{1}{2}} \right\|_W^2, \end{aligned} \right. $$

* Cross-term functional:

$$ \mathcal{B}^{n+\frac{1}{2}} = \left( \frac{s_{h,0}^{n+1} - s_{h,0}^n}{\Delta t}, B_{h,0} \frac{u_h^{n+1} - u_h^n}{\Delta t} \right)_W, $$

* Energy dissipation functional:

$$ \mathcal{D}^n = \frac{1}{2} \left\| H_{zn,h}^{\frac{1}{2}} \frac{s_{h,0}^{n+1} + s_{h,0}^n}{2} \right\|_W^2 + \frac{1}{2} \left\| H_{zn,h}^{\frac{1}{2}} \frac{s_{h,0}^n + s_{h,0}^{n-1}}{2} \right\|_W^2 + \left\| H_{zn,h}^{\frac{1}{2}} \frac{s_{h,\star}^{n+1/2} + s_{h,\star}^{n-1/2}}{2} \right\|_W^2, $$

* Total energy functional:

$$ \mathcal{E}^{n+\frac{1}{2}} = \mathcal{K}^{n+\frac{1}{2}} + \mathcal{P}_e^{n+\frac{1}{2}} + \mathcal{P}_\star^{n+\frac{1}{2}} + \mathcal{P}_0^{n+\frac{1}{2}} - \frac{\Delta t^2}{4} \mathcal{B}^{n+\frac{1}{2}}, $$

then, upon the assumption of null source, i.e. $f^n = 0$, we have the following decay property:

$$ \frac{\mathcal{E}^{n+\frac{1}{2}} - \mathcal{E}^{n-\frac{1}{2}}}{\Delta t} = -\mathcal{D}^n \le 0. $$
*Proof.* By remarking that $\{u_h^n\}_{\frac{1}{4}}$ can be used as a pivot since

$$ \left\{
\begin{aligned}
\{u_h^n\}_\theta &= \{u_h^n\}_{\frac{1}{4}} + (\theta - \frac{1}{4})(u_h^{n+1} - 2u_h^n + u_h^{n-1}), \\
u_h^n &= \{u_h^n\}_{\frac{1}{4}} - \frac{1}{4}(u_h^{n+1} - 2u_h^n + u_h^{n-1}),
\end{aligned}
\right. \qquad (28) $$

the first equation of the scheme can be rewritten as

$$ \left( I_h - \frac{\Delta t^2}{4} A_{h,0}^e + \Delta t^2 \left(\theta - \frac{1}{4}\right) A_{h,\star}^e \right) \frac{u_h^{n+1} - 2u_h^n + u_h^{n-1}}{\Delta t^2} + A_h^e \{u_h^n\}_{\frac{1}{4}} + B_{h,0}^t s_{h,0}^n + B_{h,\star}^t \frac{s_{h,\star}^{n+\frac{1}{2}} + s_{h,\star}^{n-\frac{1}{2}}}{2} = 0. $$

By taking the scalar product of this equation with $\frac{1}{2\Delta t}(u_h^{n+1} - u_h^{n-1})$, we obtain a first relation

$$ \frac{\mathcal{K}^{n+\frac{1}{2}} - \mathcal{K}^{n-\frac{1}{2}}}{\Delta t} + \frac{\mathcal{P}_e^{n+\frac{1}{2}} - \mathcal{P}_e^{n-\frac{1}{2}}}{\Delta t} + \left(B_{h,0} \frac{u_h^{n+1} - u_h^{n-1}}{2\Delta t}, s_{h,0}^n\right)_W + \left(B_{h,\star} \frac{u_h^{n+1} - u_h^{n-1}}{2\Delta t}, \frac{s_{h,\star}^{n+\frac{1}{2}} + s_{h,\star}^{n-\frac{1}{2}}}{2}\right)_W = 0. $$

Using Property 2 and Property 4 for the second and third equation in (27) respectively, and combining the two obtained expressions with the one above leads to

$$ \frac{\mathcal{K}^{n+\frac{1}{2}} - \mathcal{K}^{n-\frac{1}{2}}}{\Delta t} + \frac{\mathcal{P}_e^{n+\frac{1}{2}} - \mathcal{P}_e^{n-\frac{1}{2}}}{\Delta t} + \frac{\mathcal{P}_\star^{n+\frac{1}{2}} - \mathcal{P}_\star^{n-\frac{1}{2}}}{\Delta t} + \frac{\mathcal{P}_0^{n+\frac{1}{2}} - \mathcal{P}_0^{n-\frac{1}{2}}}{\Delta t} - \frac{\Delta t^2}{4} \frac{\mathcal{B}^{n+\frac{1}{2}} - \mathcal{B}^{n-\frac{1}{2}}}{\Delta t} = -\mathcal{D}^n \le 0. $$

<figure><img src="image_5.png" /></figure>

**Lemma 5.** Assuming that $\theta \ge \frac{1}{4}$, the energy functional associated to the time discrete scheme (27) is positive on the following condition

$$ \Delta t \le 2\mathcal{R}_{D_{zn,0}}^{-\frac{1}{2}}. \qquad (29) $$

where $D_{zn}$ is defined in (4), and $\mathcal{R}_{D_{zn,0}}$ in (19).

*Proof.* Using Property 3 we have

$$ \mathcal{K}^{n+\frac{1}{2}} + \mathcal{P}_0^{n+\frac{1}{2}} - \frac{\Delta t^2}{4} \mathcal{B}^{n+\frac{1}{2}} = \frac{1}{2} \left( \left( I_h - \frac{\Delta t^2}{4} (A_{h,0}^e + A_{h,0}(\mathbf{S}_{zn}^{-1})) + \Delta t^2 \left(\theta - \frac{1}{4}\right) A_{h,\star}^e \right) \frac{u_h^{n+1} - u_h^n}{\Delta t}, \frac{u_h^{n+1} - u_h^n}{\Delta t} \right)_H + \frac{1}{2} \left\| S_{zn,h}^{\frac{1}{2}} \frac{s_{h,0}^{n+1} + s_{h,0}^n}{2} \right\|_W^2 + \frac{\Delta t^2}{8} \left\| S_{zn,h}^{\frac{1}{2}} \frac{s_{h,0}^{n+1} - s_{h,0}^n}{\Delta t} - S_{zn,h}^{-\frac{1}{2}} B_h \frac{u_h^{n+1} - u_h^n}{\Delta t} \right\|_W^2 $$

Hence, every other terms being positive and upon the assumption that $\theta \ge \frac{1}{4}$, the total energy functional is positive upon the condition that

$$ I_h - \frac{\Delta t^2}{4} (A_{h,0}^e + A_{h,0}(\mathbf{S}_{zn}^{-1})) \ge 0. $$

The expression (5) gives us $\mathbf{S}_{zn}^{-1} = \mathbf{D}_{zn} - \mathbf{C}_{zn}$, and, from linearity of the stiffness operator w.r.t. the elastic modulii – see (9) – we obtain

$$ A_{h,0}^e + A_{h,0}(\mathbf{S}_{zn}^{-1}) = A_{h,0}(\mathbf{C}_{zn}) + A_{h,0}(\mathbf{S}_{zn}^{-1}) = A_{h,0}(\mathbf{D}_{zn}), $$

which concludes the proof.

<figure><img src="image_6.png" /></figure>

**Remark 6.** In the same vein of Remark 5, the CFL condition (29) is independent of material parameters or the mesh elements occupying $\Omega_\star$.
### 3.4. Kelvin-Voigt model

To simplify the presentation, we denote by $A_h^e = A_h(C_{kv})$ and $A_h^v = A_h(D_{kv})$ the purely elastic and purely viscous stiffness operators respectively. The fully-discrete locally implicit scheme for the KV model defined in (11) reads

$$ \frac{u_h^{n+1} - 2u_h^n + u_h^{n-1}}{\Delta t^2} + A_{h,0}^e u_h^n + A_{h,\star}^e \{u_h^n\}_\theta + A_{h,0}^v \frac{u_h^n - u_h^{n-1}}{\Delta t} + A_{h,\star}^v \frac{u_h^{n+1} - u_h^{n-1}}{2\Delta t} = f^n. \quad (30) $$

This expression of the locally implicit scheme differs from the two previous cases in the treatment of the viscous part. Here, we merely treat the dissipative term in $\Omega_\star$ using an (implicit) discrete time derivative centered at $t^n$.

**Lemma 6.** Associated to the fully-discrete solution of (30), let us define the following

* Kinetic energy functional:

$$ \mathcal{K}^{n+\frac{1}{2}} = \frac{1}{2} \left( \left( I_h - \frac{\Delta t^2}{4} A_{h,0}^e + \Delta t^2 \left( \theta - \frac{1}{4} \right) A_{h,\star}^e \right) \frac{u_h^{n+1} - u_h^n}{\Delta t}, \frac{u_h^{n+1} - u_h^n}{\Delta t} \right)_H, $$

* Potential energy functionals:

$$ \mathcal{P}_e^{n+\frac{1}{2}} = \frac{1}{2} \left( A_h^e \frac{u_h^{n+1} + u_h^n}{2}, \frac{u_h^{n+1} + u_h^n}{2} \right)_H, \quad \mathcal{P}_v^{n+\frac{1}{2}} = \left( A_{h,0}^v \frac{u_h^{n+1} - u_h^n}{\Delta t}, \frac{u_h^{n+1} - u_h^n}{\Delta t} \right)_H, $$

* Energy dissipation functional:

$$ \mathcal{D}^n = \left( A_h^v \frac{u_h^{n+1} - u_h^{n-1}}{2\Delta t}, \frac{u_h^{n+1} - u_h^{n-1}}{2\Delta t} \right)_H, $$

* Total energy functional:

$$ \mathcal{E}^{n+\frac{1}{2}} = \mathcal{K}^{n+\frac{1}{2}} + \mathcal{P}_e^{n+\frac{1}{2}} - \frac{\Delta t}{4} \mathcal{P}_v^{n+\frac{1}{2}}, $$

then, upon the assumption of null source, i.e. $f^n = 0$, we have the following decay property:

$$ \frac{\mathcal{E}^{n+\frac{1}{2}} - \mathcal{E}^{n-\frac{1}{2}}}{\Delta t} = -\mathcal{D}^n \le 0. $$

*Proof.* Using (28), the scheme reads

$$ \left( I_h - \frac{\Delta t^2}{4} A_{h,0}^e + \Delta t^2 \left( \theta - \frac{1}{4} \right) A_{h,\star}^e \right) \frac{u_h^{n+1} - 2u_h^n + u_h^{n-1}}{\Delta t^2} + A_h^e \{u_h^n\}_\frac{1}{4} + A_{h,0}^v \frac{u_h^n - u_h^{n-1}}{\Delta t} + A_{h,\star}^v \frac{u_h^{n+1} - u_h^{n-1}}{2\Delta t} = 0. $$

Taking the scalar product of this equation with $\frac{1}{2\Delta t}(u_h^{n+1} - u_h^{n-1})$ leads to

$$ \frac{\mathcal{K}^{n+\frac{1}{2}} - \mathcal{K}^{n-\frac{1}{2}}}{\Delta t} + \frac{\mathcal{P}_e^{n+\frac{1}{2}} - \mathcal{P}_e^{n-\frac{1}{2}}}{\Delta t} + \left( A_{h,\star}^v \frac{u_h^{n+1} - u_h^{n-1}}{2\Delta t}, \frac{u_h^{n+1} - u_h^{n-1}}{2\Delta t} \right)_H + \left( A_{h,0}^v \frac{u_h^{n+1} - u_h^{n-1}}{2\Delta t}, \frac{u_h^n - u_h^{n-1}}{\Delta t} \right)_H = 0. $$

Using the following relation

$$ \left( A_{h,0}^v \frac{u_h^{n+1} - u_h^{n-1}}{2\Delta t}, \frac{u_h^n - u_h^{n-1}}{\Delta t} \right)_H = \left( A_{h,0}^v \frac{u_h^{n+1} - u_h^{n-1}}{2\Delta t}, \frac{u_h^{n+1} - u_h^{n-1}}{2\Delta t} \right)_H - \frac{\Delta t}{4} \frac{\mathcal{P}_v^{n+\frac{1}{2}} - \mathcal{P}_v^{n-\frac{1}{2}}}{\Delta t}, $$

along with $A_{h,\star}^v + A_{h,0}^v = A_h^v$ by linearity of the stiffness operator w.r.t. the elastic modulii, see (9) - leads to the desired result. $\square$
**Lemma 7.** Assuming that $\theta \ge \frac{1}{4}$, the energy functional associated to the time discrete scheme (30) is positive on the following condition

$$ \Delta t \le 2\mathcal{R}_{0,\mathcal{C}_{kv}}^{-\frac{1}{2}} \left( \left( \frac{\mathcal{R}_{0,\mathcal{D}_{kv}}^2}{4\mathcal{R}_{0,\mathcal{C}_{kv}}} + 1 \right)^{\frac{1}{2}} - \frac{\mathcal{R}_{0,\mathcal{D}_{kv}}}{2\sqrt{\mathcal{R}_{0,\mathcal{C}_{kv}}}} \right). $$

*Proof.* The energy functional can be rewritten as follows

$$ \varepsilon^{n+\frac{1}{2}} = \frac{1}{2} \left( \left( I_h - \frac{\Delta t^2}{4} A_{h,0}^e - \frac{\Delta t}{2} A_{h,0}^v + \Delta t^2 \left( \theta - \frac{1}{4} \right) A_{h,\star}^e \right) \frac{u_h^{n+1} - u_h^n}{\Delta t}, \frac{u_h^{n+1} - u_h^n}{\Delta t} \right)_H + \frac{1}{2} \left( A_h^e \frac{u_h^{n-1} + u_h^n}{2}, \frac{u_h^{n-1} + u_h^n}{2} \right)_H, $$

which, assuming that $\theta \ge \frac{1}{4}$, is positive if $I_h - \frac{\Delta t^2}{4} A_{h,0}^e - \frac{\Delta t}{2} A_{h,0}^v$ is positive, equivalently if

$$ \frac{\Delta t^2}{4} \mathcal{R}_{0,\mathcal{C}_{kv}} + \frac{\Delta t}{2} \mathcal{R}_{0,\mathcal{D}_{kv}} - 1 \le 0. $$

Computing the positive root of this second-order polynomial inequality leads to the desired result. $\square$

### 3.5. Time marching algorithm

In this section, we present the algorithms performed to obtain, given the previous steps, the solution at the next time step for the three locally implicit visco-elastic schemes presented previously. To this effect, we denote by $\{\varphi_I\}_{I=1}^N$ and $\{\psi_J\}_{J=1}^{N'}$ the basis functions of the approximation spaces, with $N$ and $N'$ the dimensions of $V_h$ and $W_h$ respectively. Using these notations, we introduce the following mass in $V_h$ and weighted mass in $W_h$

$$ (\mathbb{M}_{\varrho})_{IJ} = \left( \int_{\Omega} \varrho \varphi_I \cdot \varphi_J \, d\Omega \right)_{I,J=1}^{N}, \quad (\mathbb{M}_{\mathcal{S}})_{IJ} = \left( \int_{\Omega} \mathcal{S} \psi_I : \psi_J \, d\Omega \right)_{I,J=1}^{N'} $$

for any symmetric strictly positive fourth-order tensor $\mathcal{S}$. As explained in Remark 2, these matrices become diagonal using the mass-lumping technique. We introduce the finite element approximation of the divergence operator in $\Omega_0$ as

$$ (\mathbb{B}_0)_{IJ} = \int_{\Omega_0} \psi_I : \varepsilon(\varphi_J) \, d\Omega, \quad \forall I = 1, \dots, N', \quad \forall J = 1, \dots, N, $$

and define $\mathbb{B}_*$ in the same manner. The finite element approximation of the divergence operator on the complete domain $\Omega$ is simply defined as $\mathbb{B} = \mathbb{B}_0 + \mathbb{B}_*$, so that the stiffness matrices read

$$ \mathbb{K}_C = \mathbb{B}^T \mathbb{M}_S^{-1} \mathbb{B}, \quad \mathbb{K}_{C,0} = \mathbb{B}_0^T \mathbb{M}_S^{-1} \mathbb{B}_0, \quad \mathbb{K}_{C,\star} = \mathbb{B}_*^T \mathbb{M}_S^{-1} \mathbb{B}_*, $$

with $C = S^{-1}$.

We denote by $\vec{u} \in \mathbb{R}^N$, $\vec{\sigma} \in \mathbb{R}^{N'}$ and $\vec{s} \in \mathbb{R}^{N'}$ the vectors of degrees of freedom of the discrete displacement field, stress and purely viscous stress respectively. As in (??), we denote the decomposition of the vector of degrees of freedom of a stress-type variable by $\vec{\sigma} = \vec{\sigma}_0 + \vec{\sigma}_*$. We show in the following that for each visco-elastic propagator, the locally implicit time scheme yields a linear system of the form

$$ (\mathbb{M}_{\varrho} + \beta \mathbb{K}_{A,\star}) \vec{u}^{n+1} = \vec{b}, \quad (31) $$

where $\beta \in \mathbb{R}$ is function of $\Delta t$ and $\theta$. Systems of the form (31) can be reduced to a small linear system of size of the number of degrees of freedom in $\Omega_*$. In fact, modulo a renumbering process, due to mass lumping, the matrix $\mathbb{M}_{\varrho} + \beta \mathbb{K}_{A,\star}$ is diagonal everywhere except on the degrees of freedom corresponding to basis functions with support non-vanishing in $\Omega_*$.
**Maxwell scheme.** In order to obtain the solution at the next time step for the MX case, we use the Schur complement method. Let us define

$$T^{\text{mx}} = (M_{S_{\text{mx}}} + \frac{\Delta t}{2} M_{H_{\text{mx}}})^{-1}, \quad Q^{\text{mx}} = T^{\text{mx}}(M_{S_{\text{mx}}} - \frac{\Delta t}{2} M_{H_{\text{mx}}}). \qquad (32)$$

First, introducing an auxiliary vector $\vec{\mu}_{\star} \in \mathbb{R}^{N'}$, we compute, using the second and third equation of (24)

$$\vec{\sigma}_{0}^{n} = Q^{\text{mx}}\vec{\sigma}_{0}^{n-1} + T^{\text{mx}}\mathbb{B}_{0}(\vec{u}^{n} - \vec{u}^{n-1}), \quad \vec{\mu}_{\star} = Q^{\text{mx}}\vec{\sigma}_{\star}^{n-\frac{1}{2}} - \frac{1}{2}T^{\text{mx}}\mathbb{B}_{\star}\vec{u}^{n-1},$$

such that $\vec{\sigma}_{\star}^{n+\frac{1}{2}} = \vec{\mu}_{\star} + \frac{1}{2}T^{\text{mx}}\mathbb{B}_{\star}\vec{u}^{n+1}$, which we inject into the first equation of (24) to obtain

$$(M_{\varrho} + \frac{\Delta t^2}{4} \mathbb{B}_{\star}^{\text{T}} T^{\text{mx}} \mathbb{B}_{\star}) \vec{u}^{n+1} = \Delta t^2 \vec{f}^n + M_{\varrho}(2\vec{u}^n - \vec{u}^{n-1}) - \Delta t^2 \mathbb{B}_{0}^{\text{T}} \vec{\sigma}_{0}^{n} - \frac{\Delta t^2}{2} \mathbb{B}_{\star}^{\text{T}} (\vec{\mu}_{\star} + \vec{\sigma}_{\star}^{n-\frac{1}{2}}),$$

and recover a system of the form of (31) since $\mathbb{B}_{\star}^{\text{T}} T^{\text{mx}} \mathbb{B}_{\star}$ can be interpreted as a stiffness matrix.

**Zener scheme.** We define $T^{\text{zn}}$ and $Q^{\text{zn}}$ similarly to (32). Introducing the auxiliary vector $\vec{\nu}_{\star} \in \mathbb{R}^{N'}$, we compute, using the second and third equation of (27)

$$\vec{s}_{0}^{n} = Q^{\text{zn}}\vec{s}_{0}^{n-1} + T^{\text{zn}}\mathbb{B}_{0}(\vec{u}^{n} - \vec{u}^{n-1}), \quad \vec{\nu}_{\star} = Q^{\text{zn}}\vec{s}_{\star}^{n-\frac{1}{2}} - \frac{1}{2}T^{\text{zn}}\mathbb{B}_{\star}\vec{u}^{n-1},$$

such that $\vec{s}_{\star}^{n+\frac{1}{2}} = \vec{\nu}_{\star} + \frac{1}{2}T^{\text{zn}}\mathbb{B}_{\star}\vec{u}^{n+1}$, which we inject into the first equation of (27) to obtain

$$(M_{\varrho} + \theta \Delta t^2 \mathbb{K}_{C_{\text{zn},\star}} + \frac{\Delta t^2}{4} \mathbb{B}_{\star}^{\text{T}} T^{\text{zn}} \mathbb{B}_{\star}) \vec{u}^{n+1} = \Delta t^2 \vec{f}^n + M_{\varrho}(2\vec{u}^n - \vec{u}^{n-1}) \\ - \Delta t^2 \mathbb{K}_{C_{\text{zn},0}} \vec{u}^n - \Delta t^2 \mathbb{K}_{C_{\text{zn},\star}}((1-2\theta)\vec{u}^n + \theta \vec{u}^{n-1}) \\ - \Delta t^2 \mathbb{B}_{0}^{\text{T}} \vec{s}_{0}^{n} - \frac{\Delta t^2}{2} \mathbb{B}_{\star}^{\text{T}} (\vec{\nu}_{\star} + \vec{s}_{\star}^{n-\frac{1}{2}}),$$

which is again a linear system of the form of (31).

**Kelvin-Voigt model.** In order to obtain the solution at the next time step, we rewrite the scheme (30) as follows

$$(M_{\varrho} + \theta \Delta t^2 \mathbb{K}_{C_{\text{kv},\star}} + \frac{\Delta t}{2} \mathbb{K}_{D_{\text{kv},\star}}) \vec{u}^{n+1} = \Delta t^2 \vec{f}^n + M_{\varrho}(2\vec{u}^n - \vec{u}^{n-1}) \\ - \Delta t^2 \mathbb{K}_{C_{\text{kv},0}} \vec{u}^n - \Delta t^2 \mathbb{K}_{C_{\text{kv},\star}}((1-2\theta)\vec{u}^n + \theta \vec{u}^{n-1}) - \Delta t \mathbb{K}_{D_{\text{kv},0}}(\vec{u}^n - \vec{u}^{n-1}) + \frac{\Delta t}{2} \mathbb{K}_{D_{\text{kv},\star}} \vec{u}^{n-1},$$

which is of the form of (31).

For all three models, we notice the implicit treatment of the penalizing operator. In fact, we need to solve a non-trivial linear system at each time step and the corresponding matrix, which needs to be inverted, is not diagonal even after mass-lumping. However, the size of the linear system can be reduced to the number of degrees of freedom corresponding to basis functions with support non-vanishing in $\Omega_{\star}$, which is expected to be considerably smaller than in $\Omega_{0}$. Thus, one can hope that the gain on the time step is not overly impacted by the cost of the inversion.

**Remark 7.** In the MX and ZN locally implicit schemes, the time derivative of the stress variable in $\Omega_{\star}$ is given using half time steps $n + 1/2$ and $n - 1/2$. A different choice could have been made using full time steps, e.g. $(\vec{\sigma}_{\star}^{n+1} - \vec{\sigma}_{\star}^{n-1})/2\Delta t$, also leading to a stable scheme. However in practice, using full time steps leads to the memory storage of an additional variable during the implementation. In fact, we need $\vec{\sigma}_{\star}^{n+1}$, $\vec{\sigma}_{\star}^{n}$ and $\vec{\sigma}_{\star}^{n-1}$, whereas the choice of half time steps ensures that only two variables need to be stored in memory, i.e. $\vec{\sigma}_{\star}^{n+1/2}$ and $\vec{\sigma}_{\star}^{n-1/2}$.
4. Stabilized explicit time schemes with the mortar method for visco-elastic wave problems

In essence, stabilized schemes consist in introducing a modified stiffness operator, which, by construction, renders a consistent and stable scheme with an optimized, i.e. potentially increased, CFL condition. As for the choice of the operator, Chebyshev polynomials are a good candidate to construct such modified stiffness operator. This method has been applied to second-order [27, 28] and higher order [29, 26, 15, 24] time schemes for hyperbolic problems, and also combined with domain decomposition approach. As explained previously, this last approach consists in applying the mortar element method to separate Ω₀ from Ω*, in which a stabilized leapfrog scheme – as in [14] – is used. This leads to the conservation of the discrete energy, and ultimately to a gain in the time step. In our work, we extend this approach to the visco-elastic case and present stabilized discretizations for the KV, MX and ZN models.

To start with, we consider the simpler inviscid case and explain the construction of stabilized schemes. Let us consider the homogeneous semi-discretized wave equation with non-zero initial conditions

$$
\frac{d^2 u_h}{dt^2} + A_h(C)u_h = 0, \quad u_h(0) = u_{0,h}, \quad \frac{du_h}{dt} = u_{1,h},
$$

for any symmetric strictly positive tensor C. The stabilized time discretization of this system is written as
follows

$$
\mathcal{P}^{-1}(\Delta t^2 A_h(\mathbf{C})) \frac{u_h^{n+1} - 2u_h^n + u_h^{n-1}}{\Delta t^2} + A_h(\mathbf{C}) u_h^n = 0, \quad (33)
$$

with $u_h^0$ and $u_h^1$ given, and where the polynomial operator $\mathcal{P}$ is instrumental in rendering a consistent and
stable time scheme. To this effect, let us suppose the following as in [14]

$$
\forall x \in [0; \Delta t^2 \mathcal{R}_c], \quad \mathcal{P}(x) > 0, \tag{34}
$$

$$
\mathcal{P}(0) = 1, \tag{35}
$$

$$
\exists \alpha > 0 \text{ s.t. } \alpha = \sup_{\beta \in \mathbb{R}} \{ \beta \ge 0 \mid \forall x \in [0, 4\beta^2], 0 \le x \mathcal{P}(x) \le 4 \}. \quad (36)
$$

The assumption (34) implies in particular that $\mathcal{P}(\Delta t^2 A_h(C))$ is invertible. Also, if (35) holds, we retrieve the leapfrog discretization when $\Delta t^2 \to 0$. Thus, at least formally, the scheme (33) is consistent of order two, and can be seen as a consistent perturbation of the standard leapfrog scheme.

Let us now derive the discrete energy and find the condition upon which it is positive – a necessary
condition for the stability of the scheme. We take the scalar product of (33) with $\frac{1}{2\Delta t}(u_h^{n+1} - u_h^{n-1})$, which
yields

$$
\epsilon^{n+\frac{1}{2}} = \frac{1}{2} \left( \mathcal{P}^{-1}(\Delta t^2 A_h(\mathbf{C})) \frac{u_h^{n+1} - u_h^n}{\Delta t}, \frac{u_h^{n+1} - u_h^n}{\Delta t} \right)_H + \frac{1}{2} (A_h(\mathbf{C}) u_h^n, u_h^{n+1})_H .
$$

As in [6], remarking that

$$
(A_h(C)u_h^n, u_h^{n+1})_H = \left( A_h(C) \frac{u_h^{n+1} + u_h^n}{2}, \frac{u_h^{n+1} + u_h^n}{2} \right)_H - \frac{\Delta t^2}{4} \left( A_h(C) \frac{u_h^{n+1} - u_h^n}{\Delta t}, \frac{u_h^{n+1} - u_h^n}{\Delta t} \right)_H,
$$

we can rewrite the discrete energy as follows

$$
\epsilon^{n+\frac{1}{2}} = \frac{1}{2} \left( \left( \mathcal{P}^{-1}(\Delta t^2 A_h(C)) - \frac{\Delta t^2}{4} A_h(C) \right) \frac{u_h^{n+1} - u_h^n}{\Delta t}, \frac{u_h^{n+1} - u_h^n}{\Delta t} \right)_H \right. \\
\left. + \frac{1}{2} \left( A_h(C) \frac{u_h^{n+1} + u_h^n}{2}, \frac{u_h^{n+1} + u_h^n}{2} \right)_H \right).
$$

The positivity of the discrete energy is guaranteed if $\mathcal{P}^{-1}(\Delta t^2 A_h(C)) - \frac{\Delta t^2}{4} A_h(C)$ is positive definite, which is equivalent to the following condition

$0 \le x \mathcal{P}(x) \le 4, \quad \forall x \in [0; \Delta t^2 \mathcal{R}_c].$
Upon assumption (36) on the polynomial $\mathcal{P}$, the above condition is true if

$$0 \le \Delta t^2 \mathcal{R}_C \le 4\alpha^2.$$

We then have the following CFL-like condition on the time step, which guarantees the positivity of the discrete energy,

$$\Delta t \le 2\alpha \mathcal{R}_C^{-\frac{1}{2}}, \qquad (37)$$

where $\alpha$ depends on $\mathcal{P}$ and, according to (36), $\alpha > 0$. We can expect that for specific choices of $\mathcal{P}$ – thus specific choices of $\alpha$ –, the CFL condition will be greater than the one of the classical leapfrog scheme $2\mathcal{R}_C^{-1/2}$.

### 4.1. Choice of modified stiffness operator $\mathcal{P}$

In [27, 14], the following polynomial $\mathcal{P}_l$ is proposed as the modified stiffness operator

$$\mathcal{P}_l(x) = \frac{2}{x} \left( 1 - \mathcal{T}_{l+1} \left( 1 - \frac{x}{2(l+1)^2} \right) \right)$$

where $\mathcal{T}_{l+1}$ represents the Chebyshev polynomial of the first kind of order $l+1$. In order for this scheme to be consistent and stable, $\mathcal{P}_l$ needs to satisfy the conditions (35) and (36). By using the explicit expression of $\mathcal{T}_{l+1}(x)$, the consistency condition (35) is verified. In [27], it is proven that $\mathcal{P}_l$ satisfies the stability condition (36). More precisely, it is shown that

$$0 \le x \mathcal{P}_l(x) \le 4 \quad \forall x \in [0, 4(l+1)^2].$$

Thus, we have the following CFL-like condition

$$\Delta t \le 2(l+1)\mathcal{R}_C^{-\frac{1}{2}}.$$

However, $\mathcal{P}_l$ may vanish thus fails to satisfy (34). Another choice for the modified stiffness operator is proposed in [29, 15, 24]

$$\mathcal{P}_l^{\nu}(x) = \frac{2}{x} \left[ 1 - \frac{\mathcal{T}_{l+1}(\delta_l^{\nu} - \frac{x}{\omega_l^{\nu}})}{\mathcal{T}_{l+1}(\delta_l^{\nu})} \right], \qquad (38)$$

with $\nu$ a small positive stabilization parameter. The shifting and stretching parameters $\delta_l^{\nu}$ and $\omega_l^{\nu}$ are positive and chosen following [29]

$$\delta_l^{\nu} = 1 + \frac{\nu}{(l+1)^2} \quad \text{and} \quad \omega_l^{\nu} = 2 \frac{\mathcal{T}_{l+1}'(\delta_l^{\nu})}{\mathcal{T}_{l+1}(\delta_l^{\nu})}.$$

This choice ensures that the conditions (34), (35) and (36) are fulfilled. In particular, this polynomial is strictly positive.

### 4.2. Maxwell model

Similarly as in the previous subsection, we construct stabilized schemes for the visco-elastic models. Thus, for the Maxwell model we propose the following stabilized scheme

$$\begin{cases} \mathcal{P}^{-1}(\Delta t^2 A_h) \frac{u_h^{n+1} - 2u_h^n + u_h^{n-1}}{\Delta t^2} + B_h^t \sigma_h^n = f^n, \\ S_{\text{mx},h} \frac{\sigma_h^{n+1} - \sigma_h^n}{\Delta t} + H_{\text{mx},h} \frac{\sigma_h^{n+1} + \sigma_h^n}{2} = B_h \frac{u_h^{n+1} - u_h^n}{\Delta t}, \end{cases} \qquad (39)$$

where, for simplicity, we denote $A_h = A_h(C_{\text{mx}})$.
**Lemma 8.** Let us define the following functionals

* Kinetic energy functional:

$$ \mathcal{K}^{n+\frac{1}{2}} = \frac{1}{2} \left( \mathcal{P}^{-1}(\Delta t^2 A_h) \frac{u_h^{n+1} - u_h^n}{\Delta t}, \frac{u_h^{n+1} - u_h^n}{\Delta t} \right)_H, $$

* Potential energy functional:

$$ \mathcal{P}^{n+\frac{1}{2}} = \frac{1}{4} \left\| S_{\mathrm{mx},h}^{\frac{1}{2}} \sigma_h^{n+1} \right\|_W^2 + \frac{1}{4} \left\| S_{\mathrm{mx},h}^{\frac{1}{2}} \sigma_h^n \right\|_W^2, $$

* Cross-term functional:

$$ \mathcal{B}^{n+\frac{1}{2}} = \left( B_h \frac{u_h^{n+1} - u_h^n}{\Delta t}, \frac{\sigma_h^{n+1} - \sigma_h^n}{\Delta t} \right)_W, $$

* Energy dissipation functional:

$$ \mathcal{D}^n = \frac{1}{2} \left\| H_{\mathrm{mx},h}^{\frac{1}{2}} \frac{\sigma_h^{n+1} + \sigma_h^n}{2} \right\|_W^2 + \frac{1}{2} \left\| H_{\mathrm{mx},h}^{\frac{1}{2}} \frac{\sigma_h^n + \sigma_h^{n-1}}{2} \right\|_W^2, $$

* Total energy functional:

$$ \mathcal{E}^{n+\frac{1}{2}} = \mathcal{K}^{n+\frac{1}{2}} + \mathcal{P}^{n+\frac{1}{2}} - \frac{\Delta t^2}{4} \mathcal{B}^{n+\frac{1}{2}}, $$

then, upon the assumption of null source, i.e. $f^n = 0$, we have the following decay property:

$$ \frac{\mathcal{E}^{n+\frac{1}{2}} - \mathcal{E}^{n-\frac{1}{2}}}{\Delta t} = -\mathcal{D}^n \le 0. $$

*Proof.* Taking the scalar product of the first equation of the scheme with $\frac{1}{2\Delta t}(u_h^{n+1} - u_h^{n-1})$, we obtain a first relation

$$ \frac{\mathcal{K}^{n+\frac{1}{2}} - \mathcal{K}^{n-\frac{1}{2}}}{\Delta t} + \left( B_h \frac{u_h^{n+1} - u_h^{n-1}}{2\Delta t}, \sigma_h^n \right)_W = 0. $$

Using Property 2 to the second equation of the scheme, and combining it to the above expression leads to

$$ \frac{\mathcal{K}^{n+\frac{1}{2}} - \mathcal{K}^{n-\frac{1}{2}}}{\Delta t} + \frac{\mathcal{P}^{n+\frac{1}{2}} - \mathcal{P}^{n-\frac{1}{2}}}{\Delta t} - \frac{\Delta t^2}{4} \frac{\mathcal{B}^{n+\frac{1}{2}} - \mathcal{B}^{n-\frac{1}{2}}}{\Delta t} = -\mathcal{D}^n \le 0. $$

□

**Lemma 9.** Upon (36) on $\mathcal{P}$, the energy functional associated to the time discrete scheme (39) is positive on the following condition

$$ \Delta t \le 2\alpha \mathcal{R}_{C_{\mathrm{mx}}}^{-\frac{1}{2}}. \qquad (40) $$

*Proof.* Injecting into the energy functional the expression obtained when using Property 3 to $\mathcal{P}^{n+\frac{1}{2}} - \frac{\Delta t^2}{4}\mathcal{B}^{n+\frac{1}{2}}$ gives the following

$$ \mathcal{E}^{n+\frac{1}{2}} = \frac{1}{2} \left( \left( \mathcal{P}^{-1}(\Delta t^2 A_h) - \frac{\Delta t^2}{4} A_h \right) \frac{u_h^{n+1} - u_h^n}{\Delta t}, \frac{u_h^{n+1} - u_h^n}{\Delta t} \right)_H + \frac{\Delta t^2}{8} \left\| S_{\mathrm{mx},h}^{\frac{1}{2}} \frac{\sigma_h^{n+1} - \sigma_h^n}{\Delta t} - S_{\mathrm{mx},h}^{-\frac{1}{2}} B_h \frac{u_h^{n+1} - u_h^n}{\Delta t} \right\|_W^2 + \frac{1}{2} \left\| S_{\mathrm{mx},h}^{\frac{1}{2}} \frac{\sigma_h^{n+1} + \sigma_h^n}{2} \right\|_W^2. $$

So, the energy functional is positive if $I_h - \frac{1}{4}\mathcal{P}(\Delta t^2 A_h)\Delta t^2 A_h$ is positive definite. Upon (36) on $\mathcal{P}$, this is satisfied if (40) is satisfied (see the beginning of Section 4 for more details). □

We obtain a stability condition similar to the one in the inviscid case (37). Thus, we can draw the same conclusion that the CFL condition can be increased for specific choices of $\mathcal{P}$ and $\alpha$. We show in the next section that a similar result can be obtained for the ZN case.
### 4.3. Zener Model

In the same manner as previously, we propose the following stabilized scheme for the Zener model

$$ egin{cases} \mathcal{P}^{-1}(\Delta t^2 A_h^v) \frac{u_h^{n+1} - 2u_h^n + u_h^{n-1}}{\Delta t^2} + A_h^e u_h^n + B_h^t s_h^n = f^n, \\ S_{zn,h} \frac{s_h^{n+1} - s_h^n}{\Delta t} + H_{zn,h} \frac{s_h^{n+1} + s_h^n}{2} - B_h \frac{u_h^{n+1} - u_h^n}{\Delta t} = 0, \end{cases} \quad (41) $$

where, for simplicity, we denote $A_h^e = A_h(C_{zn})$ and $A_h^v = A_h(D_{zn})$.

**Lemma 10.** Let us define the following functionals

* Kinetic energy functional:

$$ \mathcal{K}^{n+\frac{1}{2}} = \frac{1}{2} \left( \mathcal{P}^{-1}(\Delta t^2 A_h^v) \frac{u_h^{n+1} - u_h^n}{\Delta t}, \frac{u_h^{n+1} - u_h^n}{\Delta t} \right)_H, $$

* Potential energy functionals:

$$ \mathcal{P}_e^{n+\frac{1}{2}} = \frac{1}{2} (A_h^e u_h^n, u_h^{n+1})_H, \quad \mathcal{P}^{n+\frac{1}{2}} = \frac{1}{4} \left\| S_{zn,h}^{\frac{1}{2}} s_h^{n+1} \right\|_W^2 + \frac{1}{4} \left\| S_{zn,h}^{\frac{1}{2}} s_h^n \right\|_W^2, $$

* Cross-term functional:

$$ \mathcal{B}^{n+\frac{1}{2}} = \left( B_h \frac{u_h^{n+1} - u_h^n}{\Delta t}, \frac{s_h^{n+1} - s_h^n}{\Delta t} \right)_W, $$

* Energy dissipation functional:

$$ \mathcal{D}^n = \frac{1}{2} \left\| H_{zn,h}^{\frac{1}{2}} \frac{s_h^{n+1} + s_h^n}{2} \right\|_W^2 + \frac{1}{2} \left\| H_{zn,h}^{\frac{1}{2}} \frac{s_h^n + s_h^{n-1}}{2} \right\|_W^2, $$

* Total energy functional:

$$ \mathcal{E}^{n+\frac{1}{2}} = \mathcal{K}^{n+\frac{1}{2}} + \mathcal{P}_e^{n+\frac{1}{2}} + \mathcal{P}^{n+\frac{1}{2}} - \frac{\Delta t^2}{4} \mathcal{B}^{n+\frac{1}{2}}, $$

then, upon the assumption of null source, i.e. $f^n = 0$, we have the following decay property:

$$ \frac{\mathcal{E}^{n+\frac{1}{2}} - \mathcal{E}^{n-\frac{1}{2}}}{\Delta t} = -\mathcal{D}^n \le 0. $$

*Proof.* By remarking that $u_h^n = \{u_h^n\}_{\frac{1}{4}} - \frac{1}{4}(u_h^{n+1} - 2u_h^n + u_h^{n-1})$, the first equation of the scheme can be rewritten as

$$ \left( \mathcal{P}^{-1}(\Delta t^2 A_h^v) - \frac{\Delta t^2}{4} A_h^e \right) \frac{u_h^{n+1} - 2u_h^n + u_h^{n-1}}{\Delta t^2} + A_h^e \{u_h^n\}_{\frac{1}{4}} + B_h^t s_h^n = 0. $$

Taking the scalar product of this equation with $\frac{1}{2\Delta t}(u_h^{n+1} - u_h^{n-1})$, we obtain a first relation:

$$ \frac{\mathcal{K}^{n+\frac{1}{2}} - \mathcal{K}^{n-\frac{1}{2}}}{\Delta t} + \frac{\mathcal{P}_e^{n+\frac{1}{2}} - \mathcal{P}_e^{n-\frac{1}{2}}}{\Delta t} + \left( B_h \frac{u_h^{n+1} - u_h^{n-1}}{2\Delta t}, s_h^n \right)_W = 0. $$

Using Property 2 to the second equation in (41), and combining it to the above expression leads to

$$ \frac{\mathcal{K}^{n+\frac{1}{2}} - \mathcal{K}^{n-\frac{1}{2}}}{\Delta t} + \frac{\mathcal{P}_e^{n+\frac{1}{2}} - \mathcal{P}_e^{n-\frac{1}{2}}}{\Delta t} + \frac{\mathcal{P}^{n+\frac{1}{2}} - \mathcal{P}^{n-\frac{1}{2}}}{\Delta t} - \frac{\Delta t^2}{4} \frac{\mathcal{B}^{n+\frac{1}{2}} - \mathcal{B}^{n-\frac{1}{2}}}{\Delta t} = -\mathcal{D}^n \le 0. $$

□
**Lemma 11.** Upon (36) on $\mathcal{P}$, the energy functional associated to the time discrete scheme (41) is positive on the following condition

$$ \Delta t \leq 2\alpha \mathcal{R}_{D_{\text{zn}}}^{-\frac{1}{2}}. \qquad (42) $$

*Proof.* Injecting into the energy functional the expression obtained when using Property 3 to $\mathcal{P}^{n+\frac{1}{2}} - \frac{\Delta t^2}{4}\mathcal{B}^{n+\frac{1}{2}}$ entails

$$ \begin{aligned} \mathcal{E}^{n+\frac{1}{2}} = \frac{1}{2} \| S_{\text{zn},h}^{\frac{1}{2}} \frac{s_h^{n+1} + s_h^n}{2} \|_{W}^2 &+ \frac{1}{2} \left( \left( \mathcal{P}^{-1}(\Delta t^2 A_h^v) - \frac{\Delta t^2}{4} A_h^v \right) \frac{u_h^{n+1} - u_h^n}{\Delta t}, \frac{u_h^{n+1} - u_h^n}{\Delta t} \right)_H \\ &+ \frac{\Delta t^2}{8} \left\| S_{\text{zn},h}^{\frac{1}{2}} \frac{s_h^{n+1} - s_h^n}{\Delta t} - S_{\text{zn},h}^{-\frac{1}{2}} \frac{u_h^{n+1} - u_h^n}{\Delta t} \right\|_{W}^2 + \frac{1}{2} \left( A_h^e \frac{u_h^{n-1} + u_h^n}{2}, \frac{u_h^{n-1} + u_h^n}{2} \right)_H. \end{aligned} $$

So, the energy functional is positive if $I_h - \frac{1}{4}\mathcal{P}(\Delta t^2 A_h^v)\Delta t^2 A_h^v$ is positive definite. Upon (36) on $\mathcal{P}$, this is satisfied if (42) is satisfied. $\square$

### 4.4. Kelvin-Voigt model

In the same manner as previously, we consider the following stabilized scheme for the Kelvin-Voigt model

$$ \mathcal{P}^{-1}(\Delta t^2 A_h^e + 2\Delta t A_h^v) \frac{u_h^{n+1} - 2u_h^n + u_h^{n-1}}{\Delta t^2} + A_h^e u_h^n + A_h^v \frac{u_h^n - u_h^{n-1}}{\Delta t} = f^n, \qquad (43) $$

where, for simplicity, we denote $A_h^e = A_h(C_{kv})$ and $A_h^v = A_h(D_{kv})$. As in the previous cases, this polynomial satisfies (35) and renders a consistent scheme. However, the order of consistency here is lower than for the MX and ZN stabilized schemes. In fact, we see that $\mathcal{P}(\Delta t^2 A_h^e + 2\Delta t A_h^v)$ goes to zero as $\Delta t \to 0$, hence, at least formally, this scheme is consistent with order one.

**Lemma 12.** Let us define the following functionals

* Kinetic energy functional:

$$ \mathcal{K}^{n+\frac{1}{2}} = \frac{1}{2} \left( \mathcal{P}^{-1}(\Delta t^2 A_h^e + 2\Delta t A_h^v) \frac{u_h^{n+1} - u_h^n}{\Delta t}, \frac{u_h^{n+1} - u_h^n}{\Delta t} \right)_H, $$

* Potential energy functionals:

$$ \mathcal{P}_e^{n+\frac{1}{2}} = \frac{1}{2} (A_h^e u_h^{n+1}, u_h^n)_H, \qquad \mathcal{P}_v^{n+\frac{1}{2}} = \left( A_h^v \frac{u_h^{n+1} - u_h^n}{\Delta t}, \frac{u_h^{n+1} - u_h^n}{\Delta t} \right)_H, $$

* Energy dissipation functional:

$$ \mathcal{D}^n = \left( A_h^v \frac{u_h^{n+1} - u_h^{n-1}}{2\Delta t}, \frac{u_h^{n+1} - u_h^{n-1}}{2\Delta t} \right)_H, $$

* Total energy functional:

$$ \mathcal{E}^{n+\frac{1}{2}} = \mathcal{K}^{n+\frac{1}{2}} + \mathcal{P}_e^{n+\frac{1}{2}} - \frac{\Delta t}{4} \mathcal{P}_v^{n+\frac{1}{2}}, $$

then, upon the assumption of null source, i.e. $f^n = 0$, we have the following decay property:

$$ \frac{\mathcal{E}^{n+\frac{1}{2}} - \mathcal{E}^{n-\frac{1}{2}}}{\Delta t} = -\mathcal{D}^n \leq 0. $$
*Proof.* By remarking that $u_h^n = \{u_h^n\}_{\frac{1}{4}} - \frac{1}{4}(u_h^{n+1} - 2u_h^n + u_h^{n-1})$ as in the ZN case, the scheme reads

$$ \left( \mathcal{P}^{-1}(\Delta t^2 A_h^e + 2\Delta t A_h^v) - \frac{\Delta t^2}{4} A_h^e \right) \frac{u_h^{n+1} - 2u_h^n + u_h^{n-1}}{\Delta t^2} + A_h^e \{u_h^n\}_{\frac{1}{4}} + A_h^v \frac{u_h^n - u_h^{n-1}}{\Delta t} = 0. $$

Taking the scalar product of this equation with $\frac{1}{2\Delta t}(u_h^{n+1} - u_h^{n-1})$ leads to

$$ \frac{\mathcal{K}^{n+\frac{1}{2}} - \mathcal{K}^{n-\frac{1}{2}}}{\Delta t} + \frac{\mathcal{P}_e^{n+\frac{1}{2}} - \mathcal{P}_e^{n-\frac{1}{2}}}{\Delta t} + \left( A_h^v \frac{u_h^{n+1} - u_h^{n-1}}{2\Delta t}, \frac{u_h^n - u_h^{n-1}}{\Delta t} \right)_H = 0. $$

Using the following relation

$$ \left( A_h^v \frac{u_h^{n+1} - u_h^{n-1}}{2\Delta t}, \frac{u_h^n - u_h^{n-1}}{\Delta t} \right)_H = \left( A_h^v \frac{u_h^{n+1} - u_h^{n-1}}{2\Delta t}, \frac{u_h^{n+1} - u_h^{n-1}}{2\Delta t} \right)_H - \frac{\Delta t}{4} \frac{\mathcal{P}_v^{n+\frac{1}{2}} - \mathcal{P}_v^{n-\frac{1}{2}}}{\Delta t} $$

leads to the desired result.

<figure><img src="image_7.png" /></figure>

**Lemma 13.** Upon (36) on $\mathcal{P}$, the energy functional associated to the time discrete scheme (30) is positive on the following condition:

$$ \Delta t \le 2 \mathcal{R}_{C_{kv}}^{-\frac{1}{2}} \left( \left( \frac{\mathcal{R}_{D_{kv}}^2}{4 \mathcal{R}_{C_{kv}}} + \alpha^2 \right)^{\frac{1}{2}} - \frac{\mathcal{R}_{D_{kv}}}{2 \sqrt{\mathcal{R}_{C_{kv}}}} \right). \quad (44) $$

*Proof.* The energy functional can be rewritten as follows

$$ \mathcal{E}^{n+\frac{1}{2}} = \frac{1}{2} \left( \left( \mathcal{P}^{-1}(\Delta t^2 A_h^e + 2\Delta t A_h^v) - \frac{\Delta t^2}{4} A_h^e - \frac{\Delta t}{2} A_h^v \right) \frac{u_h^{n+1} - u_h^n}{\Delta t}, \frac{u_h^{n+1} - u_h^n}{\Delta t} \right)_H + \frac{1}{2} \left( A_h^e \frac{u_h^{n-1} + u_h^n}{2}, \frac{u_h^{n-1} + u_h^n}{2} \right)_H, $$

which is positive if $I_h - \frac{1}{4}\mathcal{P}(\Delta t^2 A_h^e + 2\Delta t A_h^v) (\Delta t^2 A_h^e + 2\Delta t A_h^v)$ is positive. Upon (36) on $\mathcal{P}$, this is satisfied when

$$ \Delta t^2 \mathcal{R}_{C_{kv}} + 2\Delta t \mathcal{R}_{D_{kv}} - 4\alpha^2 \le 0 $$

is satisfied. Computing the positive root of this second-order polynomial inequality w.r.t. $\Delta t$ leads to the desired result.

<figure><img src="image_8.png" /></figure>

Contrary to the MX and ZN case, the stability condition we obtain for the KV scheme is less favorable and differs from (37). In Figure 1, we plot the CFL condition of each stabilized propagator w.r.t. $\alpha$ in a toy configuration. We observe that, for the same choice of operator $\mathcal{P}$, the stability condition of the MX and ZN propagator are linear w.r.t. $\alpha$, whereas in the KV case, it is smaller up to a certain polynomial order ($\approx 20$), after which it exceeds the MX and ZN stability conditions. Though the stability condition can be more favorable as the order increases, in practice we rarely go above order 8 and in this range, the KV CFL remains less favorable.

### 4.5. Coupling through mortar element method

In practice, we use the polynomial $P_l$ (see Section 4.1) as the modified stiffness operator $\mathcal{P}$. The inverse of $P_l$ is, in fact, never needed since we directly apply $P_l$. With this choice, we observe that the CFL conditions (40), (42) and (44) increase along with the order of $P_l$. So, one would expect to obtain a considerably reduced overall CPU time as we increase $l$ as well. However, as $l$ increases, the number of applications of the stiffness operator increases as well ($l+1$ applications) and this operation is the most expensive one in this framework, since the mass matrix is lumped and the inversion of a diagonal matrix is trivial. Thus in practice, the gain in the stability condition is evened out by the cost of the multiple applications of
<figure><img src="image_9.png" /><figcaption>Figure 1: Comparison of the CFL condition of the stabilized scheme for the MX, ZN and KV propagators with the material properties $\rho = 1.0 \text{ g} \cdot \text{cm}^{-3}$, $\lambda_e = \lambda_v = 0.4 \text{ GPa}$, $\mu_e = \mu_v = 0.4 \text{ GPa}$ and the same calibration as in [6] on a regular ten by ten grid with $h = 0.1 \text{ mm}$ and fourth-order spectral elements.</figcaption></figure>

the stiffness matrix. The stabilized scheme is then equivalent to the standard leapfrog scheme in terms of precision and performance. In the case of the KV model, the (much smaller) gain in the CFL condition can even be overridden by the cost of the stiffness applications (see Figure 1). In the end, the gain in CPU time is close to zero in this context and the stabilized scheme (as such, at least) does not seem quite advantageous to our initial aim of relaxing the stability condition in order to reduce the simulation time.

In order to benefit from the increased CFL condition of stabilized schemes and reduce the computational time, the authors in [14] propose to combine stabilized schemes with mortar domain decomposition. The idea consists in localizing the stabilization in a subdomain by decomposing the domain as in the locally implicit case $\Omega = \Omega_0 \cup \Omega_*$, where $\Omega_*$ represents the most constraining elements for the stability condition. Then, a classical leapfrog time discretization is used in $\Omega_0$, whereas a stabilized scheme is used in $\Omega_*$. In fact, we expect $\Omega_*$ to be quite small compared to $\Omega_0$ and the overall cost of applying multiple sparse stiffness matrices w.r.t. the gain on the time step would be reduced. The mortar element method allows to couple different kind of formulations and weakly enforces the continuity at the interface by the introduction of a Lagrange multiplier. It renders an interesting stability condition that is also increased and less constrained by $\Omega_*$ as $l$ increases. However, the gain is limited by the inversion of a linear system at the interfaces between the subdomains, as detailed in this section.

*Coupled visco-elastic formulations.* To illustrate the method, let us consider the non-overlapping two-domain problem, i.e. $\Omega = \Omega_0 \cup \Omega_*$ with $\Omega_0 \cap \Omega_* = \emptyset$, where both subdomains support a visco-elastic Maxwell formulation (11). In $\Omega_0$, the problem is discretized in time with a classical leapfrog scheme and in $\Omega_*$, we use a stabilized scheme. The mortar element method introduces a Lagrange multiplier $\lambda$, which represents the normal trace. The continuity of the traces is weakly enforced at the interface $\Gamma = \partial\Omega_0 \cap \Omega_* - \text{assumed to be a closed boundary. Defining the following functional spaces}$

$$V_0 = [H^1(\Omega_0)]^d, \quad W_0 = [L^2(\Omega_0)]_{\text{sym}}^{d \times d}, \quad V_* = [H^1(\Omega_*)]^d, \quad W_* = [L^2(\Omega_*)]_{\text{sym}}^{d \times d}, \quad Z = [H^{-\frac{1}{2}}(\Gamma)]^d,$$
the weak formulation of this problem reads: for any time $t \ge 0$, find $(u_0(t), u_*(t), \sigma_0(t), \sigma_*(t), \lambda(t)) \in V_0 \times V_* \times W_0 \times W_* \times Z$ solution of

$$ egin{cases} 
\displaystyle \int_{\Omega_0} \rho \partial_{tt}^2 u_0 \cdot v d\Omega_0 + \int_{\Omega_0} \sigma_0 : \varepsilon(v) d\Omega_0 - \int_{\Gamma} \lambda \cdot v d\Gamma = \int_{\Omega_0} f \cdot v d\Omega_0 \\[1.5em] 
\displaystyle \int_{\Omega_*} \rho \partial_{tt}^2 u_* \cdot w d\Omega_* + \int_{\Omega_*} \sigma_* : \varepsilon(w) d\Omega_* + \int_{\Gamma} \lambda \cdot w d\Gamma = 0 \\[1.5em] 
\displaystyle \int_{\Gamma} (\partial_t u_0 - \partial_t u_*) \cdot \mu d\Gamma = 0 
\end{cases} 
\quad (45) $$

for all test functions $(v, w, \mu) \in V_0 \times V_* \times Z$ and completed with the MX constitutive relations. Using the same notation as in Section 2, the conform finite element approximation space of $Z$ is denoted $Z_h$ and is constructed with the traces of the basis functions of $V_h$. We define the coupling operator $C_{h,0}$ as

$$ (C_{h,0}\lambda_h, v_h)_{L^2(\Gamma)^d} = \int_{\Gamma} \lambda_h \cdot v_h d\Gamma, \quad \forall \lambda_h \in Z_h, \quad \forall v_h \in V_{h,0}, $$

and $C_{h,\star}$ is defined similarly. The space discretization of (45) reads

$$ egin{cases} 
\displaystyle d_{tt}^2 u_{h,0} + B_{h,0}^t \sigma_{h,0} - C_{h,0} \lambda_h = 0, & \text{in } \Omega_0, \\[1.5em] 
S_{h,0} d_t \sigma_{h,0} + H_{h,0} \sigma_{h,0} - B_{h,0} d_t u_{h,0} = 0, & \text{in } \Omega_0, \\[1.5em] 
\displaystyle d_{tt}^2 u_{h,\star} + B_{h,\star}^t \sigma_{h,\star} + C_{h,\star} \lambda_h = 0, & \text{in } \Omega_*, \\[1.5em] 
S_{h,\star} d_t \sigma_{h,\star} + H_{h,\star} \sigma_{h,\star} - B_{h,\star} d_t u_{h,\star} = 0, & \text{in } \Omega_*, \\[1.5em] 
C_{h,0}^t d_t u_{h,0} - C_{h,\star}^t d_t u_{h,\star} = 0, & \text{in } \Gamma. 
\end{cases} 
\quad (46) $$

We apply to (46) a leapfrog time discretization in $\Omega_0$ and a stabilized time discretization in $\Omega_*$, which reads

$$ egin{cases} 
\displaystyle \frac{u_{h,0}^{n+1} - 2u_{h,0}^n + u_{h,0}^{n-1}}{\Delta t^2} + B_{h,0}^t \sigma_{h,0}^n - C_{h,0} \frac{\lambda_h^{n+1} + \lambda_h^{n-1}}{2} = 0, & \text{in } \Omega_0, \\[1.5em] 
S_{h,0} \frac{\sigma_{h,0}^{n+1} - \sigma_{h,0}^n}{\Delta t} + H_{h,0} \frac{\sigma_{h,0}^{n+1} + \sigma_{h,0}^n}{2} - B_{h,0} \frac{u_{h,0}^{n+1} - u_{h,0}^n}{\Delta t} = 0, & \text{in } \Omega_0, \\[1.5em] 
\mathcal{P}^{-1}(\Delta t^2 A_{h,\star}) \frac{u_{h,\star}^{n+1} - 2u_{h,\star}^n + u_{h,\star}^{n-1}}{\Delta t^2} + B_{h,\star}^t \sigma_{h,\star}^n + C_{h,\star} \frac{\lambda_h^{n+1} + \lambda_h^{n-1}}{2} = 0, & \text{in } \Omega_*, \\[1.5em] 
S_{h,\star} \frac{\sigma_{h,\star}^{n+1} - \sigma_{h,\star}^n}{\Delta t} + H_{h,\star} \frac{\sigma_{h,\star}^{n+1} + \sigma_{h,\star}^n}{2} - B_{h,\star} \frac{u_{h,\star}^{n+1} - u_{h,\star}^n}{\Delta t} = 0, & \text{in } \Omega_*, \\[1.5em] 
C_{h,0}^t \frac{u_{h,0}^{n+1} - u_{h,0}^{n-1}}{2\Delta t} - C_{h,\star}^t \frac{u_{h,\star}^{n+1} - u_{h,\star}^{n-1}}{2\Delta t} = 0, & \text{in } \Gamma. 
\end{cases} 
\quad (47) $$

where $A_{h,\star} = A_{h,\star}(C_{\text{mx}})$.

*Stability analysis.* Let us now give an energy conservation result using energy methods as previously.

**Lemma 14.** The energy functional associated to the considered time discrete scheme is positive on the following condition

$$ \Delta t \le \min(2\mathcal{R}_{C_{\text{mx},0}}^{-\frac{1}{2}}, 2\alpha\mathcal{R}_{C_{\text{mx},*}}^{-\frac{1}{2}}) \quad (48) $$

*Proof.* Proceeding in the same manner as previously, we take the scalar product of the first and third equation of the scheme (47) with $\frac{1}{2\Delta t}(u_{h,0}^{n+1} - u_{h,0}^{n-1})$ and $\frac{1}{2\Delta t}(u_{h,\star}^{n+1} - u_{h,\star}^{n-1})$ respectively. Then, we use
Property 2 to the second and fourth equation in (47), and combine it to the previous expressions. Finally, summing the results leads to

$$ \frac{\varepsilon_0^{n+\frac{1}{2}} - \varepsilon_0^{n-\frac{1}{2}}}{\Delta t} - \frac{1}{2\Delta t} \left( C_{h,0} \frac{\lambda_h^{n+1} + \lambda_h^{n-1}}{2}, u_{h,0}^{n+1} - u_{h,0}^{n-1} \right)_{L^2(\Gamma)^d} \\ + \frac{\varepsilon_*^{n+\frac{1}{2}} - \varepsilon_*^{n-\frac{1}{2}}}{\Delta t} - \frac{1}{2\Delta t} \left( C_{h,\star} \frac{\lambda_h^{n+1} + \lambda_h^{n-1}}{2}, u_{h,\star}^{n+1} - u_{h,\star}^{n-1} \right)_{L^2(\Gamma)^d} = -D_0^n - D_*^n. $$

Using the continuity relation, which is the fifth equation of (47), leads to the following decay property

$$ \frac{\varepsilon_0^{n+\frac{1}{2}} - \varepsilon_0^{n-\frac{1}{2}}}{\Delta t} + \frac{\varepsilon_*^{n+\frac{1}{2}} - \varepsilon_*^{n-\frac{1}{2}}}{\Delta t} = -D_0^n - D_*^n. $$

Thus, the stability is guaranteed if the condition (48) is satisfied. □

In the expression (48), we can see that the stabilization is localized in $\Omega_\star$. Indeed in our framework, we expect $2\mathcal{R}_{C_{\text{mx},\star}}^{-1/2}$ to be the most constraining condition, since $\Omega_\star$ represents the most limiting cells for the stability condition. Thus, (48) will increase along with $\alpha$. Since $\Omega_\star$ is expected to be quite small in comparison to $\Omega_0$, the cost of the application of the stiffness operator $A_{h,\star}(C_{\text{mx}})$ would be smaller compared to the application of $A_{h,0}(C_{\text{mx}})$.

### 4.6. Time marching algorithm

In this section, we present the algorithm performed to obtain, given the previous steps, the solution at the next time step for the stabilized Maxwell scheme with mortars (47) presented above. To this effect, we denote by $\{\mu_I\}_{I=1}^{N_\Gamma}$ the basis functions of $Z_h$ with $N_\Gamma$ the dimension of the approximation space and we use the same matrix-vector notations as in Section 3.5. This leads to the following system

$$ \left( \begin{array}{ccc} \frac{1}{\Delta t^2} \mathbb{M}_{\varrho,0} & 0 & -\frac{1}{2} \mathbb{C}_0 \\ 0 & \frac{1}{\Delta t^2} \mathbb{M}_{\varrho,\star} & \frac{1}{2} \mathcal{P}(\Delta t^2 \mathbb{K}_{C_{\text{mx},\star}} \mathbb{M}_{\varrho,\star}^{-1}) \mathbb{C}_\star \\ \frac{1}{2\Delta t} \mathbb{C}_0^\Gamma & -\frac{1}{2\Delta t} \mathbb{C}_\star^\Gamma & 0 \end{array} \right) \left( \begin{array}{c} \vec{u}_0^{n+1} \\ \vec{u}_\star^{n+1} \\ \vec{\lambda}^{n+1} \end{array} \right) = \left( \begin{array}{c} \vec{F}_0 \\ \vec{F}_\star \\ \vec{G} \end{array} \right), \quad (49) $$

where we define the transmission matrices

$$ (C_0)_{IJ} = \int_{\Gamma} \varphi_{0,I} \cdot \boldsymbol{\mu}_J d\Gamma, \quad (C_\star)_{IJ} = \int_{\Gamma} \varphi_{\star,I} \cdot \boldsymbol{\mu}_J d\Gamma, \quad \forall I = 1, \dots, N, \quad \forall J = 1, \dots, N_\Gamma. $$

The various right-hand sides are defined as

$$ \left\{ \begin{array}{l} \vec{F}_0 = \frac{1}{\Delta t^2} \mathbb{M}_{\varrho,0} (2 \vec{u}_0^n - \vec{u}_0^{n-1}) - \mathbb{B}_0^\Gamma \vec{\sigma}_0^n + \frac{1}{2} \mathbb{C}_0 \vec{\lambda}^{n-1}, \\ \vec{F}_\star = \frac{1}{\Delta t^2} \mathbb{M}_{\varrho,\star} (2 \vec{u}_\star^n - \vec{u}_\star^{n-1}) - \mathcal{P}(\Delta t^2 \mathbb{K}_{C_{\text{mx},\star}} \mathbb{M}_{\varrho,\star}^{-1}) (\mathbb{B}_\star^\Gamma \vec{\sigma}_\star^n + \frac{1}{2} \mathbb{C}_\star \vec{\lambda}^{n-1} - \vec{f}^n), \\ \vec{G} = \frac{1}{2\Delta t} (\mathbb{C}_0^\Gamma \vec{u}_0^{n-1} - \mathbb{C}_\star^\Gamma \vec{u}_\star^{n-1}). \end{array} \right. $$

Then, we define the following matrices

$$ D = \begin{pmatrix} \frac{1}{\Delta t^2} \mathbb{M}_{\varrho,0} & 0 \\ 0 & \frac{1}{\Delta t^2} \mathbb{M}_{\varrho,\star} \end{pmatrix}, \quad C_P = \begin{pmatrix} -\frac{1}{2} \mathbb{C}_0 \\ \frac{1}{2} \mathcal{P}(\Delta t^2 \mathbb{K}_{C_{\text{mx},\star}} \mathbb{M}_{\varrho,\star}^{-1}) \mathbb{C}_\star \end{pmatrix}, \quad C^\Gamma = \left( \frac{1}{2\Delta t} \mathbb{C}_0^\Gamma - \frac{1}{2\Delta t} \mathbb{C}_\star^\Gamma \right) $$

representing the block matrices of (49). For solving such a linear system, the Schur complement consists in, first, computing the matrix $S = -C^\Gamma D^{-1} C_P$ using a LU factorization. Then, the following steps are repeated over time:
$$ \textit{Step 1. Preprocess:} \quad \begin{pmatrix} \vec{u}_0 \\ \vec{u}_* \end{pmatrix} = \mathbb{D}^{-1} \begin{pmatrix} \vec{F}_0 \\ \vec{F}_* \end{pmatrix}, $$

$$ \textit{Step 2. Solve } \mathbb{S} \vec{\lambda}^{n+1} = G - \mathbb{C} \begin{pmatrix} \vec{u}_0 \\ \vec{u}_* \end{pmatrix}, $$

$$ \textit{Step 3. Postprocess:} \quad \begin{pmatrix} \vec{u}_0^{n+1} \\ \vec{u}_*^{n+1} \end{pmatrix} = \begin{pmatrix} \vec{u}_0 \\ \vec{u}_* \end{pmatrix} - \mathbb{D}^{-1} \mathbb{C}_P \vec{\lambda}^{n+1}. $$

Using the mass-lumping technique, we are able to obtain diagonal mass matrices, thus the application of $\mathbb{D}^{-1}$ is trivial. In a classical leapfrog scheme without stabilization, one would have

$$ \mathbb{C}_P|_{P=Id} = \begin{pmatrix} -\frac{1}{2}\mathbb{C}_0 \\ \frac{1}{2}\mathbb{C}_* \end{pmatrix}. $$

As for the mass matrices, the volume-to-interface matrices $\mathbb{C}_0$ and $\mathbb{C}_*$ can be lumped in the conform case so that the application of $\mathbb{C}_P|_{P=Id}$, the inversion of $\mathbb{S}$ and its application are also quite effortless. In the stabilized framework, we notice that $\mathbb{C}_*$ is multiplied by $\mathcal{P}(\Delta t^2 \mathbb{K}_{C_{\text{mx},*}} \mathbb{M}_{\theta, *}^{-1})$, which is an operator made of stiffness matrices and cannot be lumped. Thus, the inversion of $\mathbb{S}$ becomes more complex and time-consuming than in the classical leapfrog framework. Hence, the overall gain in CPU-time induced by (48) is potentially reduced because of this inversion. However, this operation is only applied at the interfaces between subdomains and so, we can expect the loss to be limited.

**Remark 8.** We propose an alternate scheme that we call the volume-stabilized mortar scheme – as opposed to the previously introduced fully-stabilized mortar scheme – which consists in replacing the third equation of (47) by the following

$$ \frac{u_{h, \star}^{n+1} - 2u_{h, \star}^{n} + u_{h, \star}^{n-1}}{\Delta t^2} + \mathcal{P}(\Delta t^2 A_{h, \star}) B_{h, \star}^{t} \sigma_{h, \star}^{n} + C_{h, \star} \frac{\lambda_{h}^{n+1} + \lambda_{h}^{n-1}}{2} = 0, \quad \text{in } \Omega_{\star}, \quad (50) $$

In this system, the transmission operator $C_{h, \star}$ is not multiplied by the polynomial $\mathcal{P}$, thus there is no linear system to invert at the interface between $\Omega_0$ and $\Omega_*$ similarly to the classical leapfrog formulation without stabilization. This change does not affect the stability of the scheme and the efficiency and accuracy of each scheme will be investigated and compared in Section 6.

# 5. Numerical results and illustrations of the locally implicit scheme

## 5.1. Illustration of numerical convergence of the locally implicit schemes

In order to give a numerical illustration of convergence, we compute the relative $L^2$-error of the solution

$$ e = \max_n \frac{\|u_{h, \star}^n - \mathcal{I}_* \vec{u}^n\|_{L^2}}{\|u_{h, \star}^n\|_{L^2}} \quad (51) $$

where $\vec{u}_*^n$ stands for the reference solution, obtained by solving the MX, ZN or KV equations with an explicit scheme on a refined mesh $\mathcal{T}_h^*$ and $\mathcal{I}_* \vec{u}^n$ stands for the interpolation operator of the locally implicit solution on $\mathcal{T}_h^*$. This error is computed for diminishing values of the space step $h$, and so, the time step $\Delta t$ through the CFL condition. We consider the locally implicit schemes for the MX (24), ZN (27) and KV (30) models with $\theta = 1/2$, supplemented by Neumann boundary conditions on a domain $\Omega = (0, 1) \times (0, 1)$ mm with $\Omega_* = (0.4, 0.6) \times (0.4, 0.6)$ mm. We mesh the domain with regular quadrilaterals and with seventh-order spectral elements. This configuration is depicted in Figure 2. The mechanical properties of the material
in $\Omega$ are represented by the following mass density $\rho$, Lamé parameters of elasticity $\lambda_C, \mu_C$, and viscosity $\lambda_D, \mu_D$

$$ \rho = 1 \text{ g} \cdot \text{cm}^{-3}, \quad \lambda_C = 0.4 \text{ GPa}, \quad \mu_C = 0.4 \text{ GPa}, \quad \lambda_D = 800 \text{ GPa}, \quad \mu_D = 800 \text{ GPa}. $$

In this configuration, the quality factor $Q_*$ bears a dominant viscous part. This choice is made in order to emphasize the difference in convergence rate between the three visco-elastic propagators. Namely, we would like to highlight the fact that the KV model has only a convergence order of one, contrary to the MX and ZN case, which are of order two. We use the calibrated expressions of the constitutive laws for each model given in Section 2. We consider a Gaussian source point centered at (0.5, 0.5) mm with an amplitude of 1 and standard deviation of 0.024. The simulation runs during 600 iterations. First, we study the behavior of

<figure><img src="image_10.png" /><figcaption>(a) Definition of $\Omega_0$ and $\Omega_*$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(b) Mesh with seventh-order spectral elements.</figcaption></figure>

<caption>Figure 2: Configuration of the domain for the convergence study of the locally implicit scheme.</caption>

the CFL condition of each of the visco-elastic models with respect to their space step $h$ in log-scale in Figure 3a. Similarly to the 1D results in [6], the CFL condition in the KV case is quadratic with respect to $h$ and much lower than the CFL conditions for the MX and ZN propagators. Then, we represent the error (51) for each visco-elastic locally implicit scheme in Table 1 and Figure 3b. We observe the convergence rates of each scheme and retrieve an expected result: the KV propagator has a convergence rate of one, whereas the MX and ZN propagators have a convergence rate of two. Moreover, we remark that to obtain an equivalent precision as the two other models, the mesh for the KV case needs to be considerably refined and thus, has a higher computational time.

<table><caption>Table 1: Error of displacement field between locally implicit and explicit scheme for varying mesh steps.</caption><thead><tr><th>$h_*$</th><th>3.125e-3</th><th>$h$</th><th>6.25e-3</th><th>1.25e-2</th><th>2.5e-2</th><th>5e-2</th></tr></thead><tbody><tr><td>$\Delta t_{\text{mx},*}$</td><td>5.51765e-8</td><td>$\Delta t_{\text{mx}}$</td><td>1.10353e-7</td><td>2.20705e-7</td><td>4.41411e-7</td><td>8.82821e-7</td></tr><tr><td></td><td></td><td>$e$</td><td>0.00136701</td><td>0.00444655</td><td>0.0444575</td><td>0.317385</td></tr><tr><td>$\Delta t_{\text{zn},*}$</td><td>2.466945e-6</td><td>$\Delta t_{\text{zn}}$</td><td>4.93389e-6</td><td>9.86778e-6</td><td>1.97356e-5</td><td>3.94711e-5</td></tr><tr><td></td><td></td><td>$e$</td><td>0.00141822</td><td>0.00457739</td><td>0.0456138</td><td>0.315874</td></tr><tr><td>$\Delta t_{\text{kv},*}$</td><td>4.356575e-11</td><td>$\Delta t_{\text{kv}}$</td><td>1.74263e-10</td><td>6.53488e-10</td><td>2.78821e-9</td><td>1.11093e-8</td></tr><tr><td></td><td></td><td>$e$</td><td>0.00494642</td><td>0.0304165</td><td>0.15493</td><td>0.673148</td></tr></tbody></table>

## 5.2. Overview of the cells' selection process through Irons theorem

In order to construct a locally implicit scheme, we first need to determine $\Omega_*$ and $\Omega_0$, i.e. find which cells are the most constraining for the CFL condition. To do so, we perform a local spectral analysis of $\Omega$
<figure><img src="image_12.png" /><figcaption>Figure 3: Comparing the evolution of the stability condition and the convergence of the KV, MX and ZN propagators.</figcaption></figure>

and identify the elements with the highest spectral radius, which will then correspond to $\Omega_*$ and simply retrieve $\Omega_0 = \Omega \setminus \Omega_*$. We recall the following theorem, see [30, 31, 32, 33]

**Theorem 15** (Irons and Treharne). All eigenvalues $\{\lambda_i\}_{i=1,N}$ of $M_{\varrho}^{-1}K_C$ are bounded by

$$ \min_{K \in \mathcal{T}_h} \lambda_K^{\min} \le \lambda_i \le \max_{K \in \mathcal{T}_h} \lambda_K^{\max} \quad \forall i = 1, \dots, N, $$

where $M_{\varrho,K}$ and $K_{C,K}$ represent the local mass and stiffness matrices for each element $K \in \mathcal{T}_h$, $\lambda_K^{\min}$ and $\lambda_K^{\max}$ represent the minimal and maximal eigenvalues of $M_{\varrho,K}^{-1}K_{C,K}$ respectively.

Thus, a sufficient stability condition on the time step is

$$ \Delta t \le 2 \left( \max_{K \in \mathcal{T}_h} \lambda_K^{\max} \right)^{-\frac{1}{2}} \le 2 R_C^{-\frac{1}{2}}, $$

and we expect that the CFL condition would increase if $\max_{K \in \mathcal{T}_h} \lambda_K^{\max}$ decreases. To this purpose, we remove the elements with the highest spectral radius from $\varrho$ (which will then constitute $\Omega_*$) and treat them with an implicit scheme. Thus, we compute local spectral radii for each mesh element and set a threshold $\epsilon$ above which the cells are in $\Omega_*$ (and symmetrically, the ones below are in $\Omega_0$). We propose to set $\epsilon$ using simple statistical tools such as the mean value $\mu_{\mathcal{T}_h}$ and the standard deviation $\sigma_{\mathcal{T}_h}$ of the set of spectral radii defined as follows

$$ \mu_{\mathcal{T}_h} = \frac{1}{\mathcal{T}_h} \sum_{K \in \mathcal{T}_h} \mathcal{R}_{C,K}, \quad \sigma_{\mathcal{T}_h} = \sqrt{\frac{1}{\mathcal{T}_h} \sum_{K \in \mathcal{T}_h} (\mathcal{R}_{C,K} - \mu_{\mathcal{T}_h})^2}, \quad \epsilon = \mu_{\mathcal{T}_h} + \alpha \sigma_{\mathcal{T}_h}, $$

where $\mathcal{R}_{C,K}$ represents the local spectral radius on an element $K$, $\mathcal{T}_h$ is the number of elements in $\mathcal{T}_h$, and $\alpha \in \mathbb{R}$ is a constant to be defined for each configuration of interest. In practice, we run several simulations with a varying $\alpha$ in order to find the optimal ratio of implicit to explicit elements leading to the lowest CPU time.

Let us illustrate this selection process with an example. We consider a two-dimensional welded component. The considered chamfer is enclosed in a 60 by 30 mm rectangle and is composed of an isotropic aluminium material, represented by the following mass density $\rho$, Lamé parameters of elasticity $\lambda_C, \mu_C$, and viscosity $\lambda_D, \mu_D$

$$ \rho = 7.8 \text{ g} \cdot \text{cm}^{-3}, \quad \lambda_C = 57 \text{ GPa}, \quad \mu_C = 26 \text{ GPa}, \quad \lambda_D = 0.057 \text{ GPa}, \quad \mu_D = 0.026 \text{ GPa}. $$
Thus, the quality factor equals $Q_*=10^3I$. The domain is meshed with quadrilaterals as is depicted in Figure 4a, which often produces distorted elements when confronted with complex CADs. In particular, one quadrilateral is very distorted because of the chamfer geometry, due to which the CFL condition drastically decreases. In Figure 4b, we plot the spectral radius of each element and represent the threshold $\epsilon$ that varies with $\alpha$. We notice that one element stands out and has a much higher spectral radius than the rest. This corresponds, in fact, to the distorted quadrangle, which we achieve to select through our threshold and constitutes $\Omega_*$.

<figure><img src="image_13.png" /><figcaption>Figure 4: Spectral analysis of 2D chamfer configuration.</figcaption></figure>

### 5.3. Application to a 2D welded component

We propose to numerically illustrate the method by modeling the UT of the same two-dimensional welded component as in Section 5.2. We consider a Gaussian source point with an amplitude of 1 and a standard deviation of 0.014 located at 10 mm from the bottom-left corner. The simulation runs during $T = 100~\mu s$ with $\alpha = 10$. This configuration is depicted in Figure 5a. We observe in Figure 6 snapshots of the displacement field through time. In Figure 5b, we plot the Euclidean norm of the displacement field obtained with the explicit or the locally implicit scheme and observe a good agreement between the two solutions. We then compare the performances of a strictly explicit scheme to the locally implicit scheme. In order to evaluate the performance of the locally implicit scheme compared to the explicit one, we compare the ratio of time steps $\Delta t_{\text{IMEX}}/\Delta t_{\text{EX}}$, which represents the ideal optimal gain, to the ratio of CPU times $\text{CPU}_{\text{EX}}/\text{CPU}_{\text{IMEX}}$, which is the effective gain we observe. The simulations are run on a laptop PC, 12th Gen Intel(R) Core(TM) i7-12800H CPU @ 2.40GHz and the results are presented in Table 2. We observe a considerable gain on the CPU time. In fact, this configuration is particularly inclined to be solved with a locally implicit scheme as there is one very distorted mesh element, thus a high contrast in the mesh elements.

### 5.4. Application to a 3D plate with a cracked through-hole

As a more realistic numerical illustration, we propose to model a three-dimensional (3D) thin plate with a cracked through-hole as can be seen in Figure 7a. The considered domain is of size 500 mm × 500 mm × 5 mm and is composed of a visco-elastic Maxwell material, with the same mechanical properties as in Section 5.2.
<figure><img src="image_14.png" /><figcaption>Figure 5: Configuration of the 2D welded component and comparison of the displacement fields between the explicit and locally implicit solution.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 6: Snapshots of the displacement field at different times obtained with a locally implicit scheme.</figcaption></figure>

<table><caption>Table 2: Comparing the performance of the simulation between the explicit (Ex.) and locally implicit case for the 2D chamfer configuration.</caption><thead><tr><th>Scheme</th><th>Δt (μs)</th><th>CPU (s)</th><th>Δt<sub>LI</sub> / Δt<sub>Ex.</sub></th><th>CPU<sub>Ex.</sub> / CPU<sub>LI</sub></th></tr></thead><tbody><tr><td>Ex.</td><td>2.79e<sup>-6</sup></td><td>7282.6</td><td></td><td></td></tr><tr><td>LI</td><td>2.76e<sup>-3</sup></td><td>8.2009</td><td>989.4</td><td>888.0</td></tr></tbody></table>

and consider the calibrated expressions given in Section 2 for $C_{mx}$ and $D_{mx}$. The plate is meshed with hexahedra and thus, this leads to the same issues as with the quadrilaterals in 2D, which motivates the wish to relax the CFL condition. The mesh and the selected elements to be treated implicitly are depicted in Figure 8a. We observe in Figure 8b snapshots of the displacement field through time. We compare the Euclidean norm of the displacement fields obtained with the locally implicit or the explicit scheme through time. We can see in Figure 7b that both solutions are in good agreement. As previously, we present the performance in Table 3 and we observe a gain in the CPU time for a simulation ran during 10μs. We remark that as $\alpha$ diminishes (i.e. as we select more and more elements to be treated with an implicit scheme) the gap between the ideal and the effective gain widens, which is to be expected. However, this gap is considerably wider than we observed in other 2D applications. This is most probably due to the fact that the linear system to invert is much larger due to the additional dimension.
<figure><img src="image_16.png" /><figcaption>(a) Representation of the domain of interest: plate with through-hole and crack.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(b) Comparing the magnitude of the displacement field at the observation point (200, 200, 4.999) mm for the explicit and LI MX model.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>(a) Mesh of the plate with fourth-order spectral elements and $\Omega_*$ depicted in magenta.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>(b) Snapshots of the displacement field at different times.</figcaption></figure>

<caption>Figure 8: Mesh configuration of the considered plate with through-hole and crack and snapshots of the displacement field obtained with a locally implicit scheme.</caption>

# 6. Numerical results and illustrations of the stabilized explicit scheme

## 6.1. Illustration of numerical convergence of the stabilized explicit schemes

In order to give a numerical illustration of convergence, we compute the relative $L^2$-error between the solution of a stabilized and explicit scheme. We consider the stabilized schemes for the MX (39), ZN (41) and KV (43) models, supplemented by Neumann boundary conditions on a domain $\Omega = (0, 1) \times (0, 1)$ mm divided into two subdomains: $\Omega_0 = (0, 0.5) \times (0, 1)$ mm and $\Omega_* = (0.5, 1) \times (0, 1)$ mm. We mesh the domain
<table><caption>Table 3: Comparing the performance of the simulation between the explicit (Ex.) and implicit case with different values for α.</caption><thead><tr><th>Scheme (α)</th><th>Δt (μs)</th><th>CPU (s)</th><th>Δt<sub>LI</sub></th><th>CPU<sub>Ex.</sub></th></tr></thead><tbody><tr><td>Ex.</td><td>1.60 e<sup>-4</sup></td><td>873.6</td><td></td><td></td></tr><tr><td>LI (4)</td><td>4.42 e<sup>-4</sup></td><td>587.5</td><td>2.7</td><td>1.4</td></tr><tr><td>LI (1)</td><td>7.65 e<sup>-4</sup></td><td>500.3</td><td>4.7</td><td>1.7</td></tr><tr><td>LI (0)</td><td>2.43 e<sup>-3</sup></td><td>250.7</td><td>15</td><td>3.4</td></tr><tr><td>LI (-0.1)</td><td>6.62 e<sup>-3</sup></td><td>177.7</td><td>41</td><td>4.9</td></tr><tr><td>LI (-0.11)</td><td>1.19 e<sup>-2</sup></td><td>178.4</td><td>74</td><td>4.8</td></tr></tbody></table>

with regular quadrilaterals and use seventh-order spectral elements. This configuration is depicted in Figure 9a. We consider in $\Omega$ a material with the same mechanical properties as in Section 5.1 and a Gaussian source point centered at (0.25, 0.25) mm with an amplitude of 1 and standard deviation of 0.014907. The simulation runs during 600 iterations. We represent the error (51) for each visco-elastic stabilized scheme in Table 4 and Figure 9b. We obtain satisfying convergence rates for each scheme, i.e. a second order convergence for the MX, ZN models and a first order convergence for the KV model.

<table><caption>Table 4: Error of displacement field between stabilized and explicit scheme for varying mesh steps for $l = 2$.</caption><thead><tr><th>h<sub>*</sub></th><th>3.125e<sup>-3</sup></th><th>h</th><th>6.25e<sup>-3</sup></th><th>1.25e<sup>-2</sup></th><th>2.5e<sup>-2</sup></th><th>5e<sup>-2</sup></th></tr></thead><tbody><tr><td>Δt<sub>mx,*</sub></td><td>5.51765e<sup>-8</sup></td><td>Δt<sub>mx</sub></td><td>1.10353e<sup>-7</sup></td><td>2.20705e<sup>-7</sup></td><td>4.41411e<sup>-7</sup></td><td>8.82821e<sup>-7</sup></td></tr><tr><td></td><td></td><td>e</td><td>0.00123612</td><td>0.00370529</td><td>0.00863152</td><td>0.0187048</td></tr><tr><td>Δt<sub>zn,*</sub></td><td>2.466945e<sup>-6</sup></td><td>Δt<sub>zn</sub></td><td>4.93389e<sup>-6</sup></td><td>9.86778e<sup>-6</sup></td><td>1.97356e<sup>-5</sup></td><td>3.94711e<sup>-5</sup></td></tr><tr><td></td><td></td><td>e</td><td>0.00126264</td><td>0.00378311</td><td>0.00880516</td><td>0.0190392</td></tr><tr><td>Δt<sub>kV,*</sub></td><td>4.356575e<sup>-11</sup></td><td>Δt<sub>kV</sub></td><td>1.74263e<sup>-10</sup></td><td>6.53488e<sup>-10</sup></td><td>2.78821e<sup>-9</sup></td><td>1.11093e<sup>-8</sup></td></tr><tr><td></td><td></td><td>e</td><td>0.00491722</td><td>0.0229545</td><td>0.107838</td><td>0.494968</td></tr></tbody></table>

<figure><img src="image_20.png" /><figcaption>(a) Configuration of the domain for the convergence study of (b) Convergence rates of the stabilized scheme for visco-elastic models with polynomial order $l = 1$ and $l = 2$.</figcaption></figure>

<caption>Figure 9: Illustration of numerical convergence of the stabilized scheme with mortar elements.</caption>

## 6.2. Overview of the domain decomposition process

In the following numerical illustrations, we consider a domain $\Omega$ divided into $N_{sub}$ subdomains $\{\Omega_i\}_{i=1,N_{sub}}$. We need to determine which subdomains constitute $\Omega_*$, i.e., are the most penalizing for the CFL condition
and are then treated with a stabilized scheme. These subdomains are denoted by {$\Omega_{*,i}$}$_{i=1,N_{*,sub}}$ with $N_{*,sub}$ the number of subdomains in $\Omega_{*}$. These notations naturally translate to $\Omega_{0}$ where an explicit scheme is used. In the present work, we select the penalizing subdomains qualitatively by examining the mesh of the considered domain. The use of simple statistical tools as in the locally implicit case (see Section 5.2) will be considered in a future work but in the following, we assume known the decomposition on $\Omega_{0}$ and {$\Omega_{*,i}$}$_{i=1,N_{*,sub}}$. Once that $\Omega_{0}$ and $\Omega_{*}$ are determined, we solve an explicit scheme in $\Omega_{0}$ and a stabilized one in $\Omega_{*}$, as in (39), (41) and (43), with $\mathcal{P}_{l}$ as the modified stiffness operator. As for the order $l$ of this polynomial, it can either be set as the same for all subdomains {$\Omega_{*,i}$}$_{i=1,N_{*,sub}}$ or can be allowed to vary from one to another, which is referred to as the "automatic mode" in the following. For this latter choice, we propose the following polynomial order {$l_{\Omega_{*,i}}$}$_{i=1,N_{*,sub}}$

$$l_{\Omega_{*,i}} = \frac{(\mathcal{R}_{0,C_{\text{mx}}})^{-\frac{1}{2}}}{(\mathcal{R}_{\Omega_{*,C_{\text{mx}}}}^{i})^{-\frac{1}{2}}}, \quad \forall i = \{1, \dots, N_{*,\text{sub}}\},$$

where $\mathcal{R}_{\Omega_{*,i}}^{i}$ represents the spectral radius of a stabilized subdomain $\Omega_{*,i}$ and

$$\mathcal{R}_{0,C_{\text{mx}}} = \min_{j=1,N_{0,\text{sub}}} \mathcal{R}_{0,C_{\text{mx}}}^{j}.$$

### 6.3. Application to a 3D plate with stiffener and through-hole

We propose to numerically illustrate the stabilization method by modeling the UT of a 3D thin plate with a stiffener and a through-hole. This configuration is depicted in Figure 11a. The plate is of size 300 mm × 300 mm × 1 mm, while the stiffener has the dimensions given in Figure 11b and is placed on top of the plate. The domain is composed of a visco-elastic Maxwell material with the same mechanical properties as in Section 5.2. We consider the calibrated expressions given in Section 2 for $C_{\text{mx}}$ and $D_{\text{mx}}$. The mesh of the plate and stiffener is composed of hexahedra and is subdivided into $N_{\text{sub}} = 168$ subdomains. The close proximity of the through-hole to the stiffener produces very thin elements, which are likely to have a high spectral radius and impact the CFL condition severely. Moreover, the area around the source is usually difficult to mesh and produces distorted mesh elements. Thus, we select the subdomains containing the thin elements along the stiffener and the distorted elements around the source to be in $\Omega_{*}$, i.e. be treated with a stabilized scheme, as depicted in Figure 12a. We observe snapshots of the displacement field through time in Figure 13. We confirm that using $\mathcal{P}_{l}$, which is not invertible, instead of $\mathcal{P}_{l}^{\nu}$ does not affect the results. (see Section 4.1). In fact, we compare the Euclidean norm of the displacement fields obtained with the stabilized and the explicit scheme through time. We can see in Figure 12b that both solutions are in good agreement. Using the same markers as previously, we present the performance of the stabilized scheme compared to the explicit one in Table 5. We test out the stabilized scheme with an unique order $l = 1$ for all subdomains in $\Omega_{*}$ and also with the automatic mode. In the latter case, the computed optimal order for the thin hexahedra along the stiffener is $l = 3$ and for the area around the source, we obtain $l = 1$. We observe a gain, however it is not very close to the ideal gain. This can be explained by the cost of the inversion of a system at the interfaces, which escalates quickly in 3D. This remark led us to test out the alternate method referred to as the volume-stabilized mortar scheme (see Remark 8). We run the same simulation as above with this new system and compare both solutions along with their performance. In Figure 10, we plot a reference solution obtained with an explicit mortar scheme without stabilization and the solutions resulting from both mortar fully-stabilized and volume-stabilized schemes. We compare the cross-section of each solution along a predefined axis at final time. We observe that, while the solution of the fully-stabilized discretization is in good agreement with the reference one, the volume-stabilized solution, though sufficiently close to the reference, loses precision over time. In fact, this corresponds to a lack of consistency which is to be expected with the change we operated. However in terms of performance, the volume-stabilized scheme is much more efficient as can be seen in Table 6. The gain we obtain is, in fact, much closer to the ideal one – greater in some cases – which is explained by the lack of linear system to invert. Thus, the volume-stabilized scheme can be a good alternative in specific low-frequency configurations, where the precision does not suffer as much and the time gain is considerable.
<figure><img src="image_21.png" /><figcaption>Figure 10: Comparing the magnitude of the displacement field at the observation point (71, 65, 0.999) mm for the explicit, fully-stabilized and volume-stabilized MX model with mortars.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Figure 11: Configuration of the domain of interest.</figcaption></figure>

<table><caption>Table 5: Comparing the performance of the simulation between the explicit (Ex.) and fully-stabilized (Ful.-stab.) schemes.</caption><thead><tr><th>Scheme (order)</th><th>Δt (μs)</th><th>CPU (s)</th><th>Δt<sub>Stab.</sub> / Δt<sub>Ex.</sub></th><th>CPU<sub>Ex.</sub> / CPU<sub>Stab.</sub></th></tr></thead><tbody><tr><td>Ex.</td><td>5.741e<sup>-3</sup></td><td>889.0</td><td></td><td></td></tr><tr><td>Ful.-stab. (2)</td><td>1.148e<sup>-2</sup></td><td>686.8</td><td>1.99</td><td>1.29</td></tr><tr><td>Ful.-stab. (auto)</td><td>1.829e<sup>-2</sup></td><td>543.6</td><td>3.18</td><td>1.63</td></tr></tbody></table>
<figure><img src="image_23.png" /><figcaption>Figure 12: Construction of $\Omega_*$ and comparison of the explicit and stabilized solution of the Maxwell model for the plate with stiffener and through-hole.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Figure 13: Snapshots of the displacement field at different times.</figcaption></figure>

## 7. Conclusion and perspectives

In this work, we introduced numerical schemes that allow the CFL condition to be less restrictive when confronted with stiffener parts in visco-elastic problems. These methods are meant to work best on problems with a high contrast between stiff and non-stiff operators, albeit not a general tool for performance. In the first part, we have proposed a locally implicit approach, where the most constraining elements for the CFL condition are treated with an implicit scheme, and thus, do not impact the stability condition anymore.
<table><caption>Table 6: Comparing the performance of the simulation between the explicit (Ex.) and volume-stabilized (Vol.-stab.) schemes.</caption><thead><tr><th>Scheme (order)</th><th>Δt (μs)</th><th>CPU (s)</th><th>Δt<sub>Stab.</sub></th><th>CPU<sub>Ex.</sub></th></tr></thead><tbody><tr><td>Ex.</td><td>5.741e<sup>-3</sup></td><td>690.4</td><td></td><td></td></tr><tr><td>Vol.-stab. (2)</td><td>1.148e<sup>-2</sup></td><td>344.0</td><td>1.99</td><td>2.00</td></tr><tr><td>Vol.-stab. (auto)</td><td>1.829e<sup>-2</sup></td><td>210.1</td><td>3.18</td><td>3.28</td></tr></tbody></table>

We have shown that the locally implicit schemes for the MX, ZN and KV models are stable upon a better CFL condition than in the fully explicit scheme. In 2D, we have given a numerical convergence illustration and we have proposed to model the UT of a welded component along with a performance study. In 3D, we have considered a thin plate with a cracked through-hole. We observe an improvement in the computational cost, whether in 2D or in 3D. However, there is a gap between the effective gain and ideal (theoretical) gain in CPU time, which is mostly due to the inversion of a linear system. This effect is more noticeable in 3D, as the number of elements drastically increase. So, as the number of elements treated with an implicit scheme increase, we can expect this gap to widen and the computational cost to increase. This constitutes a limitation of this method. In the second part, we have proposed a stabilized approach with mortar domain decomposition. In this case, the most constraining mesh elements are separated from the rest of the domain through a mortar domain decomposition and are treated with a stabilized scheme. In the latter, the idea is to apply a polynomial stiffness operator, which characteristics make it so that the CFL condition is less restrictive as the order of the polynomial increases. In fact, the gain in the CFL condition is less in the KV case. We have shown using energy arguments that the stabilized schemes for the MX, ZN and KV models are stable, and that the introduction of the mortar method (imperative to exploit the gain in CFL of the stabilized schemes) does not affect this stability. In 2D, we have given a numerical convergence illustration and in 3D, we have proposed to model the UT of a thin plate with a stiffener and a through-hole, along with a performance study. We observe an improvement in the computational cost compared to the fully explicit numerical scheme, close to the one ideally expected. Moreover, the inversion of a linear system at the mortar interfaces does not seem to affect the result too much for this configuration. However, this method can only be used in specific configurations with a possible domain decomposition.

Several perspectives following this work can be considered, whether in a more practical point of view or in the sense of broadening the proposed methods to further notions. Currently, there is a parameter $\alpha$ that needs to be set in the locally implicit framework. This parameter needs to be tuned in order to find the most optimal distribution of elements in $\Omega_*$ and $\Omega_0$. Ideally, $\alpha$ would be computed using the given information on the mesh and the propagator. But the process of this computation is still unclear and an interesting perspective would be to go into this in depth. In the case of the stabilized method, we select the subdomains that seem more constraining to the CFL condition and tag them to be treated with a stabilized scheme. The same process as in the locally implicit case could be established, i.e., compute the spectral radius of each subdomain and select the ones above a threshold parameterized by a constant $\alpha$. Then, will be left the same issue as in the previous case, which is the choice of $\alpha$. But, this constitutes a better automation than the actual process. As explained above, each of these methods work best in specific configurations. So, we could imagine identifying which scheme (implicit, locally implicit, stabilized or explicit) is best compatible with which kind of configuration (UT of welded components, SHM, ...). Another axis to explore would be the locally stabilized schemes [25, 34, 15, 24] (or local time stepping) applied to visco-elastic propagators. This approach does not involve domain decomposition, which eliminates the cost of inverting a linear system at the interfaces between $\Omega_0$ and $\Omega_*$. However, this method does not have a priori stability result [25, 29]. This could be an interesting approach for configurations such as the ones used in the locally implicit case, though with an added cost as $\Omega_*$ may need to be enriched to guarantee stability, it could constitute an additional tool to enhance the CFL condition. For further prospects, we could consider high-order methods, which are often used to increase the efficiency of the algorithm or the precision of the solution. In the case of the locally implicit scheme, the authors in [35] have proposed a high-order locally implicit scheme for linear ODEs and
in particular semi-discretized wave propagation models. In the case of the stabilized approach, high-order schemes are proposed in [29] for wave problems using the modified equation approach. Another interesting axis of research would be considering non-conform meshes. In the case of the stabilized schemes, the mortar domain decomposition method would have to be adapted to take into account the non-conformity.

References

[1] http://www.extende.com/civa-in-a-few-words.

[2] D. Komatitsch, J.-P. Vilotte, The spectral element method: an efficient tool to simulate the seismic response of 2d and 3d geological structures, Bull. Seismol. Soc. Am. 88 (2) (1998) 368–392.

[3] G. Cohen, P. Joly, J. E. Roberts, N. Tordjman, Higher order triangular finite elements with mass lumping for the wave equation, SIAM J. Numer. Anal. 38 (6) (2001) 2047–2078.

[4] P. Joly, Numerical Methods for Elastic Wave Propagation, in: Waves in Nonlinear Pre-Stressed Materials, no. 495 in CISM Courses and Lectures, Springer Vienna, 2007, pp. 181–281.

[5] J. M. Carcione, Wave fields in real media: Wave propagation in anisotropic, anelastic, porous and electromagnetic media, Vol. 38, Elsevier, 2007.

[6] A. Imperiale, N. Leymarie, E. Demaldent, Numerical modeling of wave propagation in anisotropic viscoelastic laminated materials in transient regime: Application to modeling ultrasonic testing of composite structures, Int. J. Numer. Methods Eng. 121 (15) (2020) 3300–3338.

[7] P. Joly, Variational Methods for Time-Dependent Wave Propagation Problems, in: Topics in Computational Wave Propagation, no. 31 in Notes Comput. Sci. Eng., Springer Berlin Heidelberg, 2003, pp. 201–264.

[8] G. Cohen, Higher-Order Numerical Methods for Transient Wave Equations, Scientific computation, Springer, 2002.

[9] S. J. Ruuth, Implicit-explicit methods for time-dependent pde's, Ph.D. thesis, University of British Columbia (1993).

[10] J. Frank, W. Hundsdorfer, J. G. Verwer, On the stability of implicit-explicit linear multistep methods, Appl. Numer. Math. 25 (2-3) (1997) 193–205.

[11] T. Rylander, A. Bondeson, Stability of explicit-implicit hybrid time-stepping schemes for maxwell's equations, J. Comput. Phys. 179 (2) (2002) 426–438.

[12] M. Hochbruck, A. Sturm, Error analysis of a second-order locally implicit method for linear maxwell's equations, SIAM J. Numer. Anal. 54 (5) (2016) 3167–3191.

[13] S. Descombes, S. Lanteri, L. Moya, Locally implicit time integration strategies in a discontinuous galerkin method for maxwell's equations, J. Sci. Comput. 56 (2013) 190–218.

[14] J. Chabassier, S. Imperiale, Construction and convergence analysis of conservative second order local time discretisation for linear wave equations, ESAIM: Math. Model. Numer. Anal. 55 (4) (2021) 1507–1543.

[15] C. Carle, M. Hochbruck, A. Sturm, On leapfrog-chebyshev schemes, SIAM J. Numer. Anal. 58 (4) (2020) 2404–2433.

[16] F. Ben Belgacem, Y. Maday, The mortar element method for three dimensional finite elements, ESAIM: Math. Model. Numer. Anal. 31 (2) (1997) 289–302.
[17] F. Casadei, E. Gabellini, G. Fotia, F. Maggio, A. Quarteroni, A mortar spectral/finite element method for complex 2d and 3d elastodynamic problems, Comput. Methods in Appl. Mech. Eng. 191 (45) (2002) 5119–5148.

[18] P. Hauret, P. Le Tallec, A discontinuous stabilized mortar method for general 3d elastic problems, Comput. Methods in Appl. Mech. Eng. 196 (49–52) (2007) 4881–4900.

[19] E. Blanc, D. Komatitsch, E. Chaljub, B. Lombard, Z. Xie, Highly accurate stability-preserving optimization of the Zener viscoelastic model, with application to wave propagation in the presence of strong attenuation, Geophys. J. Int. 205 (1) (2016) 427–439.

[20] E. Bécache, A. Ezziani, P. Joly, A mixed finite element approach for viscoelastic wave propagation, Computat. Geosci. 8 (3) (2005) 255–299.

[21] M. Durufle, P. Grob, P. Joly, Influence of Gauss and Gauss-Lobatto quadrature rules on the accuracy of a quadrilateral finite element method in the time domain, Numer. Methods Partial Differential Equ. 25 (3) (2009) 526–551.

[22] M. J. S. Chin-Joe-Kong, W. A. Mulder, M. Van Veldhuizen, Higher-order triangular and tetrahedral finite elements with mass lumping for solving the wave equation, J. Eng. Math. 35 (4) (1999) 405–426.

[23] W. A. Mulder, R. Shamasundar, Performance of continuous mass-lumped tetrahedral elements for elastic wave propagation with and without global assembly, Geophys. J. Int. 207 (1) (2016) 414–421.

[24] M. Grote, S. Michel, S. Sauter, Stabilized leapfrog based local time-stepping method for the wave equation, Math. Comput. 90 (332) (2021) 2603–2643.

[25] J. Diaz, M. J. Grote, Energy conserving explicit local time stepping for second-order wave equations, SIAM J. Sci. Comput. 31 (3) (2009) 1985–2014.

[26] J. Chabassier, S. Imperiale, Space/time convergence analysis of a class of conservative schemes for linear wave equations, C. R. Math. 355 (3) (2017) 282–289.

[27] J. C. Gilbert, P. Joly, Higher order time stepping for second order hyperbolic problems and optimal cfl conditions, in: Partial Differential Equations: Modeling and Numerical Simulation, Springer, 2008, pp. 67–93.

[28] P. Joly, J. Rodríguez, Optimized higher order time discretization of second order hyperbolic problems: Construction and numerical study, J. Comput. Appl. Math. 234 (6) (2010) 1953–1961.

[29] S. Imperiale, Etude mathématique et numérique de capteurs piézoélectriques, Ph.D. thesis, PhD Thesis, Université Paris Dauphine (2012).

[30] I. Fried, Bounds on the extremal eigenvalues of the finite element stiffness and mass matrices and their spectral condition number, J. Sound Vib. 22 (4) (1972) 407–418.

[31] A. Wathen, An analysis of some element-by-element techniques, Comput. Methods Appl. Mech. Eng. 74 (3) (1989) 271–287.

[32] J. I. Lin, An element eigenvalue theorem and its application for stable time steps, Comput. Methods Appl. Mech. Eng. 73 (3) (1989) 283–294.

[33] R. Cottereau, R. Sevilla, Stability of an explicit high-order spectral element method for acoustics in heterogeneous media based on local element stability criteria, Int. J. Numer. Methods Eng. 116 (4) (2018) 223–245.
[34] J. Chabassier, J. Diaz, S. Imperiale, Construction and analysis of fourth order, energy consistent, family of explicit time discretizations for dissipative linear wave equations, ESAIM: Math. Model. Numer. Anal. 54 (3) (2020) 845-878.

[35] H. Barucq, M. Duruflé, M. N'diaye, High-order locally a-stable implicit schemes for linear odes, J. Sci. Comput. 85 (2) (2020) 31.