# An implicit-explicit time discretization for elastic wave propagation problems in plates

Hajer Methenni, Alexandre Imperiale, Sébastien Imperiale

► To cite this version:

Hajer Methenni, Alexandre Imperiale, Sébastien Imperiale. An implicit-explicit time discretization for elastic wave propagation problems in plates. International Journal for Numerical Methods in Engineering, 2024, 125, pp.e7393. 10.1002/nme.7393. hal-04291680

HAL Id: hal-04291680

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

Submitted on 17 Nov 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# An implicit-explicit time discretization for elastic wave propagation problems in plates

Hajer Methenni<sup>‡,a,b</sup>, Alexandre Imperiale<sup>b,*</sup>, Sébastien Imperiale<sup>c</sup>

<sup>a</sup>POEMS (CNRS – Inria – ENSTA Paristech), Institut Polytechnique de Paris, 91128 Palaiseau, France

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

<sup>c</sup>Inria — LMS, Ecole Polytechnique, CNRS — Université Paris-Saclay

## Abstract

We propose a new implicit-explicit scheme to address the challenge of modelling wave propagation within thin structures using the time-domain finite element method. Compared to standard explicit schemes, our approach renders a time marching algorithm with a time step independent of the plate thickness and its associated discretization parameters (mesh step and order of approximation). Relying on the standard three dimensional elastodynamics equations, our strategy can be applied to any type of material, either isotropic or anisotropic, with or without discontinuities in the thickness direction. Upon the assumption of an extruded mesh of the plate-like geometry, we show that the linear system to be solved at each time step is partially lumped thus efficiently treated. We provide numerical evidence of an adequate convergence behavior, similar to a reference solution obtained using the well-known leapfrog scheme. Further numerical investigations show significant speed up factors compared to the same reference scheme, proving the efficiency of our approach for the configurations of interest.

**Keywords:** Elastodynamics, plates, spectral finite elements, implicit-explicit time discretisation

## 1. Introduction

Mechanical wave propagation is a physical phenomenon upon which rely various industrial applications. Examples are numerous, across different scientific and technical fields: soil imaging in geophysics, medical imaging in healthcare applications, or ultrasonic testing experiments. The main advantage of mechanical waves is their capacity to propagate over large distances, and during their course to be sensitive to the fluctuations of the material characteristics (variations in sound velocity or in geometry, presence of defects...). After propagation, signals bearing capital information on the material at hand are recorded and analyzed. Recovering quantitative data from this analyzing phase is a complex task. A typical family of configurations generating complex signals is guided wave propagation in plate-like geometries, where the propagation domain occupies a volume with a small thickness compared to the transversal wavelengths. Uses and applications associated to these waves can be found in the non-destructive testing community [1], where metallic or composite panels are inspected for flaws or mechanical weaknesses. In these cases, the signals are complex mainly because of the dispersive behavior of the so-called Lamb modes, and their analysis often requires modelling and simulation [2]. Modelling wave propagation in thin structures can be carried out using different approaches. A first relevant approach is to derive – prior to any discretization procedure – an asymptotic model, potentially taking the form of a specific set of partial differential equations. The study of plate deformation is a very standard subject on which many references can be found – see for example [3, 4, 5, 6, 7, 8]. However, fewer studies concern the use of asymptotic methods for wave propagation in

<footnote>*Corresponding author: alexandre.imperiale@cea.fr, ‡ now at Touch Sensity, Technopole Bordeaux Technowest, 33700 Mérignac.</footnote>
these media – see for instance [9]. A second possibility is to directly address the standard three dimensional elastodynamics equation and to propose an efficient mean for approximating its solution. Even though numerous methods fall into this category – e.g. the so-called Semi-Analytical Finite Element (SAFE) methods or the Finite Difference (FD) schemes and their extensions, see [2] for a review on other potential strategies – we focus in the scope of our work on the time-domain Spectral Finite Element Method (SFEM) [10, 11, 12]. This approach has proven to be particularly efficient for long distance propagation and accurate rendering of wave packets dispersion [13, 14, 15], and it has been successfully applied to guided wave propagation modelling [16, 17]. In essence, it combines a finite element space discretization based upon high-order basis functions defined on Gauss-Lobatto points with an explicit time difference scheme – the most popular scheme being the leapfrog scheme. Using the so-called mass-lumping technique [11, 12, 18, 19, 20], this approach leads to a very efficient algorithm, benefiting from the precision of spectral methods and the efficiency of explicit FD schemes. However, as any explicit time discretization scheme, in the context of propagation in thin structures it suffers from a stability condition on the time step – also referred to as the CFL condition – that becomes more stringent as the thickness of the plate decreases. This aspect becomes critical for very thin plates, or for plates with stratified materials, which is typically the case for composite structures.

The goal of the presented work is to derive a time-domain finite element method for solving wave prop-
agation problems in thin plates that is both robust w.r.t to the thickness and efficient in terms of precision
and computational performances. Broadly, our approach is structured around the following elements. We
identify the part of the stiffness bilinear form that strongly penalizes the stability condition on the time
step of the explicit scheme. This is done by re-scaling the problem into a geometry independent of the
plate half-thickness, denoted by $\delta$ in the following. This leads to an equivalent re-scaled formulation of the
problem where the dependency of the differential operators w.r.t. $\delta$ is clearly visible. Applying the standard
explicit scheme on this re-scaled formulation shows that its related stability condition is proportionate to $\delta$.
After revisiting this well-known constraint on the explicit scheme, and comparing it with the unconditionally
stable implicit scheme, we apply an implicit-explicit scheme – see e.g. the related prior works of [21, 22, 23].
For this implicit-explicit scheme we are able to prove that there exists a sufficient stability condition that is
independent of $\delta$. In fact, after careful analysis of this stability condition, we show that it only depends on the
discretization parameters (mesh steps, mesh elements ratio, order of approximation, ...) of the mid-surface.
We also verify that upon some assumptions on the mesh of the plate (mainly extrusion of the mid-surface
mesh) we conserve an efficient scheme, i.e. the computational cost overhead of solving the implicit part of
the scheme is substantially limited. This aspect is due to the partial lumping of the part of the stiffness
bilinear form appearing in the linear system to be solved at each time step of the implicit-explicit scheme.
Combining robustness of the scheme as the plate thickness tends to zero with efficiency enables us to obtain
significant improvements in terms of overall performances compared to the standard leapfrog scheme, used
as a reference algorithm in what follows.

The structure of the paper is as follows. In Section 2 we present the weak form of the linear elastic wave propagation problem and its re-scaled formulation to a geometry independent of the half-thickness $\delta$. We also introduce the associated semi-discrete formulation based upon a Galerkin approximation of the appropriate solution space. The precise description of this approximation space using finite elements is postponed to a subsequent section. Upon this general semi-discrete re-scaled formulation, we review in Section 3 the standard implicit and explicit schemes, and we introduce a new implicit-explicit scheme. In this section we derive a sufficient stability condition of this scheme, defined using only the partial derivatives in the tangent plane – thus independent of $\delta$. Additionally, as a collateral result, we are able to envision a relevant preconditioning strategy for the implicit scheme based upon the matrix to be inverted as each time step of the implicit-explicit scheme. In Section 4 we specify the nature of the Galerkin approximation space. We start by introducing the precise notations for describing the mid-surface mesh and its extrusion, leading to a three dimensional mesh of the re-scaled geometry. We also recall the standard construction of high-order spectral finite elements and their associated notions of numerical integration. With these notions at hand, we are able to show the partial lumping of the part of the stiffness bilinear form appearing in the linear system of the implicit-explicit scheme. This is a key element to insure good performances of this scheme, in addition to its robustness w.r.t. $\delta$. We also provide precise estimates on the Rayleigh quotients appearing
in the stability conditions of both the explicit and implicit-explicit schemes. This shows in particular that the latter benefits from a time step that depends only on the discretization parameters of the mid-surface mesh. In the last Section 5 we present some numerical results. As first numerical investigations, we compare the convergence of the implicit-explicit scheme to the leapfrog scheme both in $L^2$ and semi-$H^1$ norms, and observe similar convergence behaviors. On a second set of experiments, we compare the performances of the two schemes on a simple isotropic plate and on a stratified anisotropic plate. In these illustrations we clearly observer a significant gain in the computational performances of the implicit-explicit scheme, thus proving the interest of our approach in such configurations. Finally, we give some conclusions and perspectives of our work in Section 6.

## 2. Semi-discrete elastic wave propagation problem

<figure><img src="image_2.png" /><figcaption>Figure 1: Illustration of the re-scaling of the domain of interest w.r.t. the plate half-thickness $\delta$, and the corresponding mesh elements $\kappa_{\ell} \in T_h(\bar{\Sigma})$ and $K_{\ell,m} \in T_h(\bar{\Omega})$ of the mid-surface and the re-scaled domain respectively.</figcaption></figure>

As shown in Figure 1, we consider a plate occupying the open and bounded volume $\Omega^\delta = \Sigma \times (-\delta, \delta) \subset \mathbb{R}^3$, where $\Sigma \subset \mathbb{R}^2$ is the mid-surface. We also consider a maximal time $T > 0$, so that $(0, T)$ represents the time interval of interest. Within this space-time domain, we assume that the propagation of high-frequency waves are governed by the standard linear elastodynamics equations

$$ \rho^\delta \partial_{tt}^2 \mathbf{u}^\delta - \mathrm{div}(\mathbf{C}^\delta \boldsymbol{\varepsilon}(\mathbf{u}^\delta)) = 0, \quad \text{in } \Omega^\delta \times (0, T), \qquad (1) $$

where $\mathbf{u}^\delta$ is the time and space dependent displacement associated to the mechanical waves, and $\boldsymbol{\varepsilon}(\cdot)$ is the linearized Green-Lagrange tensor. The solid is characterized by its mass density $\rho^\delta$, and the fourth-order tensor $\mathbf{C}^\delta$ representing the linear stress - strain relation, i.e. Hooke's law. Denoting by $\Gamma_+^\delta$ the upper boundary of the domain $\Omega^\delta$, relation (1) is completed with the following boundary conditions

$$ (\mathbf{C}^\delta \boldsymbol{\varepsilon}(\mathbf{u}^\delta)) \cdot \mathbf{n} = \delta \mathbf{g}_+, \quad \text{on } \Gamma_+^\delta, \quad (\mathbf{C}^\delta \boldsymbol{\varepsilon}(\mathbf{u}^\delta)) \cdot \mathbf{n} = 0, \quad \text{on } \partial\Omega \setminus \Gamma_+^\delta, $$
where **n** is the normal vector field of the boundary. For any point $(x_1, x_2)^T \in \Gamma_+^\delta$, the time-dependent source term $g_+(x_1, x_2, t)$ typically represents the effect of an ultrasonic actuator. For simplicity we assume that

$$g_+(x_1, x_2, 0) = 0, \quad \forall (x_1, x_2)^T \in \Gamma_+^\delta, \qquad (2)$$

and we also provide null initial conditions in displacement field and velocity, namely

$$\mathbf{u}^{\delta}(\mathbf{x}, 0) = \partial_t \mathbf{u}^{\delta}(\mathbf{x}, 0) = 0, \quad \forall \mathbf{x} \in \Omega^{\delta}. \qquad (3)$$

As described in Figure 1, we transform each point $(x_1, x_2, \delta x_3)^T \in \Omega^\delta$ in a point $(x_1, x_2, x_3)^T \in \Omega = \Sigma \times (-1, 1)$. Note that $\Omega$ is independent of the domain half-thickness, and for this “fixed” geometry we denote by $\Gamma_+$ its upper boundary. Concerning the physical parameters, we assume that there exists a scalar field $\rho$ and a fourth-order tensor field $\mathbf{C}$ defined in $\Omega$ such that

$$\rho^{\delta}(x_1, x_2, x_3) = \rho(x_1, x_2, \frac{x_3}{\delta}), \quad \mathbf{C}^{\delta}(x_1, x_2, x_3) = \mathbf{C}(x_1, x_2, \frac{x_3}{\delta}), \quad \forall (x_1, x_2, x_3)^T \in \Omega^{\delta}. \qquad (4)$$

The scalar field $\rho$ is a strictly positive bounded function independent of $\delta$ and we assume that there exists two constants $\rho_-$ and $\rho_+$ such that

$$0 < \rho_- < \rho_-(x) < \rho_+ < +\infty, \quad \forall x \in \bar{\Omega}. \qquad (5)$$

Concerning the fourth-order tensor field $\mathbf{C}$ we assume that it is uniformly coercive and bounded. More precisely, denoting by $\mathcal{L}(\mathbb{R}^3)$ the space of tensors – assimilated to linear mappings from $\mathbb{R}^3$ into itself – equipped with its natural inner product and norm

$$\forall \boldsymbol{\sigma}, \boldsymbol{\varepsilon} \in \mathcal{L}(\mathbb{R}^3), \quad \boldsymbol{\sigma} : \boldsymbol{\varepsilon} = \sum_{i,j=1}^{3} \sigma_{ij} \varepsilon_{ij}, \quad |\boldsymbol{\varepsilon}| = \sqrt{\boldsymbol{\varepsilon} : \boldsymbol{\varepsilon}},$$

and by $\mathcal{L}_s(\mathbb{R}^3)$ the subspace of $\mathcal{L}(\mathbb{R}^3)$ of symmetric tensors, the mentioned uniform properties read: there exists two constants $\alpha > 0$ and $\beta > 0$ such that

$$\forall \boldsymbol{\varepsilon} \in \mathcal{L}_s(\mathbb{R}^3), \quad \alpha |\boldsymbol{\varepsilon}|^2 \le \mathbf{C}(\boldsymbol{x}) \boldsymbol{\varepsilon} : \boldsymbol{\varepsilon} \le \beta |\boldsymbol{\varepsilon}|^2, \quad \forall \boldsymbol{x} \in \bar{\Omega}. \qquad (6)$$

Additionally, in what follows we assume that the scalar field $\rho$ and the tensor field $\mathbf{C}$ are either globally continuous functions over the whole domain, or at least piece-wise continuous over mesh elements. With these assumptions – along with standard regularity hypothesis on the source term and the domain geometry – one can verify that the problem is well-posed and a unique smooth solution exists (see [12]).

## 2.1. Re-scaled continuous model problem

**Lemma 1.** Defining $\mathbf{V} = H^1(\Omega)^3$, the re-scaled weak formulation of problem (1) reads: for any time $t \in (0, T)$, find $\mathbf{u}(t) \in \mathbf{V}$ such that for any $\mathbf{v} \in \mathbf{V}$

$$\frac{d^2}{dt^2} m(\mathbf{u}, \mathbf{v}) + a^{\delta}(\mathbf{u}, \mathbf{v}) = \ell(t; \mathbf{v}), \qquad (7)$$

with

$$\left\{
\begin{aligned}
m(\mathbf{w}, \mathbf{v}) &= \int_{\Omega} \rho \mathbf{w} \cdot \mathbf{v} \, d\Omega, \\
a^{\delta}(\mathbf{w}, \mathbf{v}) &= \int_{\Omega} \mathbf{C} \boldsymbol{\varepsilon}^{\delta}(\mathbf{w}) : \boldsymbol{\varepsilon}^{\delta}(\mathbf{v}) \, d\Omega, \quad \forall \mathbf{w}, \mathbf{v} \in \mathbf{V} \\
\ell(t; \mathbf{v}) &= \int_{\Gamma_+} \mathbf{g}_+(t) \cdot \mathbf{v} \, d\Gamma_+,
\end{aligned}
\right.$$

and

$$\boldsymbol{\varepsilon}^{\delta}(\mathbf{v}) = \frac{1}{2} \begin{pmatrix} 2\partial_1 v_1 & \cdot & \cdot \\ \partial_2 v_1 + \partial_1 v_2 & 2\partial_2 v_2 & \cdot \\ \delta^{-1}\partial_3 v_1 + \partial_1 v_3 & \delta^{-1}\partial_3 v_2 + \partial_2 v_3 & 2\delta^{-1}\partial_3 v_3 \end{pmatrix}, \quad \forall \mathbf{v} = (v_1, v_2, v_3)^T \in \mathbf{V}.$$
*Proof.* Let us denote by **u** the re-scaled displacement defined on the fixed geometry. It still depends on $\delta$ even though we do not make explicit this dependency in order to lighten the presentation. In particular, we have

$$ \mathbf{u}(x_1, x_2, x_3, t) = \mathbf{u}^{\delta}(x_1, x_2, \delta x_3, t), \quad \forall (x_1, x_2, x_3)^{\intercal} \in \Omega, \ t \in (0, T). $$

This change of variable transfer the dependency w.r.t $\delta$ of the geometry onto the differential operators. Namely, the partial derivatives in space of the re-scaled displacement read

$$ \partial_t \mathbf{u}^{\delta}(x_1, x_2, x_3, t) = \partial_t \mathbf{u}(x_1, x_2, x_3, t), \quad i = 1, 2, \quad \partial_3 \mathbf{u}^{\delta}(x_1, x_2, x_3, t) = \delta^{-1} \partial_3 \mathbf{u}(x_1, x_2, x_3, t). $$

Let $\mathbf{u}^{\delta}$ be the solution of (1) such that $\mathbf{u}^{\delta}(t) \in \mathbf{V}^{\delta} = H^1(\Omega^{\delta})^3$, it satisfies the following weak formulation

$$ \forall \mathbf{v}^{\delta} \in \mathbf{V}^{\delta} \quad \int_{\Omega^{\delta}} \rho^{\delta} \frac{\partial^2 \mathbf{u}^{\delta}}{\partial t^2} \cdot \mathbf{v}^{\delta} \, d\Omega^{\delta} + \int_{\Omega^{\delta}} \mathbf{C}^{\delta} \varepsilon(\mathbf{u}^{\delta}) \cdot \varepsilon(\mathbf{v}^{\delta}) \, d\Omega^{\delta} = \delta \int_{\Gamma_{+}^{\delta}} \mathbf{g}_{+} \cdot \mathbf{v}^{\delta} \, d\Gamma_{+}^{\delta}. $$

After operating the change of variable from $\Omega^{\delta}$ to $\Omega$, and upon the assumption (4) on the mass density and the stress – strain, we obtain the following equivalent weak formulation

$$ \forall \mathbf{v} \in \mathbf{V} \quad \delta \int_{\Omega} \rho \frac{\partial^2 \mathbf{u}}{\partial t^2} \cdot \mathbf{v} \, d\Omega + \delta \int_{\Omega} \mathbf{C} \varepsilon^{\delta}(\mathbf{u}) : \varepsilon^{\delta}(\mathbf{v}) \, d\Omega = \delta \int_{\Gamma_{+}} \mathbf{g}_{+} \cdot \mathbf{v} \, d\Gamma_{+}, $$

which entails the desired result after simplification by $\delta$. $\square$

From the expression (7) it is possible to provide the explicit form of the dependency of the stiffness bilinear form in terms of power of $\delta$. This will be the fundamental building block for the developments of dedicated time discretizations. For any $\mathbf{v} = (v_1, v_2, v_3)^{\intercal} \in \mathbf{V}$ we define the tangential and normal part of the Green-Lagrange tensor as

$$ \varepsilon_{\tau}(\mathbf{v}) = \frac{1}{2} \begin{pmatrix} 2\partial_1 v_1 & \cdot & \cdot \\ \partial_2 v_1 + \partial_1 v_2 & 2\partial_2 v_2 & \cdot \\ \partial_1 v_3 & \partial_2 v_3 & 0 \end{pmatrix}, \quad \varepsilon_{\nu}(\mathbf{v}) = \frac{1}{2} \begin{pmatrix} 0 & \cdot & \cdot \\ 0 & 0 & \cdot \\ \partial_3 v_1 & \partial_3 v_2 & 2\partial_3 v_3 \end{pmatrix}. \quad (8) $$

Denoting by

$$ a_{\alpha\beta}(\mathbf{w}, \mathbf{v}) = \int_{\Omega} \mathbf{C} \varepsilon_{\alpha}(\mathbf{w}) : \varepsilon_{\beta}(\mathbf{v}) \, d\Omega, \quad \forall \mathbf{w}, \mathbf{v} \in \mathbf{V}, \quad \forall \alpha \in \{\tau, \nu\}, \quad \forall \beta \in \{\tau, \nu\}, $$

and remarking that the re-scaled Green-Lagrange tensor can be expressed as

$$ \varepsilon^{\delta}(\mathbf{v}) = \varepsilon_{\tau}(\mathbf{v}) + \delta^{-1} \varepsilon_{\nu}(\mathbf{v}), \quad (9) $$

the bilinear form $a^{\delta}(\cdot, \cdot)$ can be rewritten as

$$ a^{\delta}(\mathbf{w}, \mathbf{v}) = a_{\tau\tau}(\mathbf{w}, \mathbf{v}) + \delta^{-1}(a_{\tau\nu}(\mathbf{w}, \mathbf{v}) + a_{\nu\tau}(\mathbf{w}, \mathbf{v})) + \delta^{-2}a_{\nu\nu}(\mathbf{w}, \mathbf{v}). $$

## 2.2. Galerkin approximation and semi-discrete weak formulation

Let us denote by $\mathbf{V}_h \subset \mathbf{V}$ a Galerking approximation of the solution space. A construction of this finite dimensional space using high order spectral finite elements is detailed in Section 4. The semi-discrete counter part of the weak formulation (7) reads: for any time $t \in (0, T)$, find $\mathbf{u}_h(t) \in \mathbf{V}_h$ such that for any $\mathbf{v}_h \in \mathbf{V}_h$

$$ \frac{d^2}{dt^2} m(\mathbf{u}_h, \mathbf{v}_h) + a^{\delta}(\mathbf{u}_h, \mathbf{v}_h) = \ell(t; \mathbf{v}_h). $$

In practice, the integrals appearing in the bilinear forms are approximated using an appropriate quadrature formula – leading in particular to the so called mass-lumping technique. Denoting by $m_h(\cdot, \cdot)$ and $a_h^{\delta}(\cdot, \cdot)$ the resulting approximated forms, the “effective” semi-discrete formulation is

$$ \frac{d^2}{dt^2} m_h(\mathbf{u}_h, \mathbf{v}_h) + a_h^{\delta}(\mathbf{u}_h, \mathbf{v}_h) = \ell(t; \mathbf{v}_h), \quad (10) $$

and in the next section we will make the following assumption.
**Assumption 1.** For any $\epsilon > 0$ and for any $\mathbf{v}_h \in \mathbf{V}_h$, we have the following inequality

$$2a_{h,\tau\nu}(\mathbf{v}_h, \mathbf{v}_h) \leq \frac{1}{\epsilon}a_{h,\tau\tau}(\mathbf{v}_h, \mathbf{v}_h) + \epsilon a_{h,\nu\nu}(\mathbf{v}_h, \mathbf{v}_h). \quad (11)$$

Note that one can verify – see proof of Lemma 3 – that the inequality (11) is satisfied by the exact bilinear forms. Hence, satisfying Assumption 1 exclusively depends on the choice of the quadrature formula used to approximate the integrals.

### 3. Robust time discretizations w.r.t. the plate thickness

The goals of this section are twofold: (1) we aim at emphasizing the performance bottlenecks that arise from a standard – and commonly used – explicit or implicit time discretizations for wave propagation in thin plates. In particular we are interested to configurations where the plate thickness tends to zero; (2) relying on the decomposition of the stiffness bilinear form appearing the semi-discrete formulation (10), we propose a time scheme that provide robustness w.r.t. the plate thickness and potentially alleviate the aforementioned bottlenecks. A stability analysis is provided to illustrate in which sense the scheme is robust w.r.t. $\delta$.

#### 3.1. Explicit, implicit, and implicit-explicit schemes

Let $\Delta t$ be the time discretization step. We assume that there exists a maximal time step index $N_T \in \mathbb{N}^*$ such that $T = N_T \Delta t$. For any integer $n \in [0; N_T]$, we denote by $t^n = n\Delta t$ and $\mathbf{u}_h^n$ represents the approximation of $\mathbf{u}_h(t^n)$ the solution of (10). According to the null initial conditions (3) and (2) we have null initial vectors, namely

$$\mathbf{u}_h^1 = \mathbf{u}_h^0 = 0. \quad (12)$$

Below we present three time discretizations of (10). An implicit and explicit time discretizations, that correspond to standard time discretizations. We present their advantages and weaknesses in order to motivate the third scheme – an implicit-explicit time discretization – that aims at avoiding the weaknesses of the fully implicit or explicit schemes while preserving their strengths.

**Implicit scheme.** Let us define $\theta \in \mathbb{R}$ a time discretization parameter and $\{\mathbf{u}_h^n\}_\theta = \theta \mathbf{u}_h^{n+1} + (1-2\theta)\mathbf{u}_h^n + \theta \mathbf{u}_h^{n-1}$ another approximation of $\mathbf{u}_h(t^n)$. The first time scheme considered in our work is the implicit scheme

$$m_h(\frac{\mathbf{u}_h^{n+1} - 2\mathbf{u}_h^n + \mathbf{u}_h^{n-1}}{\Delta t^2}, \mathbf{v}_h) + a_h^\delta(\{\mathbf{u}_h^n\}_\theta, \mathbf{v}_h) = \ell(t^n, \mathbf{v}_h), \quad (13)$$

with null initial vectors (12). This scheme is unconditionally stable for any $\theta \ge \frac{1}{4}$ – which is in itself an interesting robustness property as $\delta$ vanishes. It can be shown (see [24]) that consistency error is minimized for $\theta = \frac{1}{4}$. The scheme then reads

$$m_h(\frac{\mathbf{u}_h^{n+1} - 2\mathbf{u}_h^n + \mathbf{u}_h^{n-1}}{\Delta t^2}, \mathbf{v}_h) + a_h^\delta(\{\mathbf{u}_h^n\}_{\frac{1}{4}}, \mathbf{v}_h) = \ell(t^n, \mathbf{v}_h). \quad (14)$$

It is again standard to show that, when $\ell(t^n, \mathbf{v}_h) = 0$ one has an energy conservation property,

$$\mathcal{E}_{\text{im},h}^{n+\frac{1}{2}} = \mathcal{E}_{\text{im},h}^{n-\frac{1}{2}}$$

with

$$\mathcal{E}_{\text{im},h}^{n+\frac{1}{2}} = \frac{1}{2}m_h\left(\frac{\mathbf{u}_h^{n+1}-\mathbf{u}_h^n}{\Delta t}, \frac{\mathbf{u}_h^{n+1}-\mathbf{u}_h^n}{\Delta t}\right) + \mathcal{P}_h^{n+\frac{1}{2}} \quad \text{with} \quad \mathcal{P}_h^{n+\frac{1}{2}} = \frac{1}{2}a_h^\delta\left(\frac{\mathbf{u}_h^{n+1}+\mathbf{u}_h^n}{2}, \frac{\mathbf{u}_h^{n+1}+\mathbf{u}_h^n}{2}\right). \quad (15)$$

The discrete energy is clearly a positive quadratic functional for any couple $(\mathbf{u}_h^{n+1}, \mathbf{u}_h^n)$. This proves the stability of the scheme (see [12, 25]) independently of the space discretization parameters, and of $\delta$. Although this property is highly desirable it comes with the disadvantage that at each iteration the linear system deriving from the bilinear form $(m_h + \frac{\Delta t^2}{4}a_h^\delta)(\cdot, \cdot)$ must be solved. To avoid such problem a common strategy is to use a fully explicit scheme
*Explicit scheme.* The second time scheme that we consider is the second order, centered explicit scheme, a.k.a. the leap-frog scheme, namely

$$m_h(\frac{\mathbf{u}_h^{n+1} - 2\mathbf{u}_h^n + \mathbf{u}_h^{n-1}}{\Delta t^2}, \mathbf{v}_h) + a_h^\delta(\mathbf{u}_h^n, \mathbf{v}_h) = \ell(t^n; \mathbf{v}_h). \quad (16)$$

When mass lumping is achieved (see [11, 26, 12, 20]) this scheme is fully explicit, in the sense that no linear system should be solved at each time step. However its performance deteriorates as the thickness tends to zero since the stability condition associated to this scheme is proportional to $\delta$. To see that, the scheme is rewritten as a second-order perturbation of (14). Using that

$$\mathbf{u}_h^n = \{\mathbf{u}_h^n\}_{\frac{1}{4}} - \frac{1}{4}(\mathbf{u}_h^{n+1} - 2\mathbf{u}_h^n + \mathbf{u}_h^{n+1})$$

we obtain that (16) is equivalent to

$$(m_h - \frac{\Delta t^2}{4} a_h^\delta) (\frac{\mathbf{u}_h^{n+1} - 2\mathbf{u}_h^n + \mathbf{u}_h^{n-1}}{\Delta t^2}, \mathbf{v}_h) + a_h^\delta(\{\mathbf{u}_h^n\}_{\frac{1}{4}}, \mathbf{v}_h) = \ell(t^n; \mathbf{v}_h). \quad (17)$$

As before, when $\ell(t^n, \mathbf{v}_h) = 0$ one has an energy conservation property,

$$\mathcal{E}_{\text{ex},h}^{n+\frac{1}{2}} = \mathcal{E}_{\text{ex},h}^{n-\frac{1}{2}} \quad \text{with} \quad \mathcal{E}_{\text{ex},h}^{n+\frac{1}{2}} = \frac{1}{2}(m_h - \frac{\Delta t^2}{4} a_h^\delta) \left( \frac{\mathbf{u}_h^{n+1} - \mathbf{u}_h^n}{\Delta t}, \frac{\mathbf{u}_h^{n+1} - \mathbf{u}_h^n}{\Delta t} \right) + \mathcal{P}_h^{n+\frac{1}{2}},$$

where $\mathcal{P}_h^{n+\frac{1}{2}}$ is defined in (15). The scheme is stable if the energy functional defined above is positive. This can be ensure choosing $\Delta t$ sufficiently small so that $(m_h - \frac{\Delta t^2}{4} a_h^\delta)(\cdot, \cdot)$ is a positive bilinear form. This gives a condition on the time step that is the so-called CFL condition. In details this condition reads

$$\Delta t \le \bar{\Delta} t_{\text{ex}} \quad \text{with} \quad \bar{\Delta} t_{\text{ex}} = 2 \left( \sup_{\mathbf{v}_h \in \mathbf{V}_h} \frac{a_h^\delta(\mathbf{v}_h, \mathbf{v}_h)}{m_h(\mathbf{v}_h, \mathbf{v}_h)} \right)^{-\frac{1}{2}}. \quad (18)$$

The time step $\bar{\Delta} t_{\text{ex}}$ depends on the space discretization, and is proportional to the small parameter $\delta$. This is rigorously shown in the next section where a space discretization with high order spectral finite elements is considered and analyzed. This last point is a major drawback in situations where $\delta$ is small, and is overcome by using the scheme presented below.

*Implicit-explicit scheme.* The last time scheme we consider consists in treating implicitly the part of the stiffness bilinear form that depends solely on the derivative in the thickness direction while the remaining part is treated explicitly. This approach is linked to the one presented in [21, 22, 23] and is referred to as an implicit-explicit strategy. It reads

$$m_h(\frac{\mathbf{u}_h^{n+1} - 2\mathbf{u}_h^n + \mathbf{u}_h^{n-1}}{\Delta t^2}, \mathbf{v}_h) + (a_h^\delta - \delta^{-2}a_{\nu\nu,h})(\mathbf{u}_h^n, \mathbf{v}_h) + \delta^{-2}a_{\nu\nu,h}(\{\mathbf{u}_h^n\}_\theta, \mathbf{v}_h) = \ell(t^n; \mathbf{v}_h), \quad (19)$$

completed with null initialization (12). As we will show in the next section, this scheme is stable upon a condition on the time step that does not depend on the thickness of the plate. We remark that, at each time step, a linear system emanating from $(m_h + \theta \Delta t^2 \delta^{-2} a_{\nu\nu,h})(\cdot, \cdot)$ has to be solved. It is however, as we will show in Lemma 4, a partially lumped bilinear form. Hence, one can expect this scheme to bear the lowest cost overhead per iteration compared to the explicit scheme while having the desired robustness w.r.t. $\delta$.

**Remark 1.** We also propose in Section 3.3 an iterative resolution strategy for the implicit scheme with an efficient preconditioning that relies on the same linear system appearing in each iterations of the implicit-explicit scheme.
### 3.2. Stability analysis of the implicit-explicit scheme

As for the explicit scheme, we start by rewriting (19) as a perturbation of the scheme (14). Using $\{\mathbf{u}_h^n\}_\theta = \mathbf{u}_h^n + \theta(\mathbf{u}_h^{n+1} - 2\mathbf{u}_h^n + \mathbf{u}_h^{n+1})$, we obtain that (19) is equivalent to

$$ \tilde{m}_h \left( \frac{\mathbf{u}_h^{n+1} - 2\mathbf{u}_h^n + \mathbf{u}_h^{n-1}}{\Delta t^2}, \mathbf{v}_h \right) + a_h^\delta (\{\mathbf{u}_h^n\}_\frac{1}{4}, \mathbf{v}_h) = \ell(t^n; \mathbf{v}_h). \quad (20) $$

with

$$ \tilde{m}_h(\mathbf{u}_h, \mathbf{v}_h) = \left( m_h - \frac{\Delta t^2}{4} \left( a_{\tau\tau,h} + 2\delta^{-1}a_{\tau\nu,h} + \delta^{-2}(1-4\theta)a_{\nu\nu,h} \right) \right) (\mathbf{u}_h, \mathbf{v}_h). $$

As previously, when $\ell(t^n, \mathbf{v}_h) = 0$ one has an energy conservation property,

$$ \mathcal{E}_{\text{imex},h}^{n+\frac{1}{2}} = \mathcal{E}_{\text{imex},h}^{n-\frac{1}{2}} \quad \text{with} \quad \mathcal{E}_{\text{imex},h}^{n+\frac{1}{2}} = \frac{1}{2} \tilde{m}_h \left( \frac{\mathbf{u}_h^{n+1} - \mathbf{u}_h^n}{\Delta t}, \frac{\mathbf{u}_h^{n+1} - \mathbf{u}_h^n}{\Delta t} \right) + \mathcal{P}_h^{n+\frac{1}{2}}, $$

where $\mathcal{P}_h^{n+\frac{1}{2}}$ is defined in (15). The stability of the scheme is therefore obtained by ensuring that $\tilde{m}_h(\cdot, \cdot)$ is a positive bilinear form, which leads to the following CFL condition

$$ \Delta t \le \overline{\Delta t}_{\text{imex}} \quad \text{with} \quad \overline{\Delta t}_{\text{imex}} = 2 \left( \sup_{\mathbf{v}_h \in \mathbf{V}_h} \frac{(a_{\tau\tau,h} + 2\delta^{-1}a_{\tau\nu,h} + (1-4\theta)\delta^{-2}a_{\nu\nu,h})(\mathbf{v}_h, \mathbf{v}_h)}{m_h(\mathbf{v}_h, \mathbf{v}_h)} \right)^{-\frac{1}{2}}. \quad (21) $$

When choosing $\theta = \frac{1}{4}$ one can see that $\overline{\Delta t}_{\text{imex}}$ depends linearly on $\sqrt{\delta}$. This is already an improvement compared to the explicit scheme. However, by choosing $\theta > \frac{1}{4}$ we show below that a sufficient condition on the time step independent of $\delta$ can be derived – upon Assumption 1.

**Theorem 1.** Assuming that $\theta > \frac{1}{4}$, and that Assumption 1 is satisfied, then $\mathcal{E}_{\text{imex}}^{n+\frac{1}{2}}$ is a positive functional if $\Delta t \le \overline{\Delta t}_{\text{imex}}$ with

$$ \overline{\Delta t}_{\text{imex}} = 2 \left( 1 + \frac{1}{4\theta - 1} \right)^{-\frac{1}{2}} \left( \sup_{\mathbf{v}_h \in \mathbf{V}_h} \frac{a_{h,\tau\tau}(\mathbf{v}_h, \mathbf{v}_h)}{m_h(\mathbf{v}_h, \mathbf{v}_h)} \right)^{-\frac{1}{2}}. \quad (22) $$

*Proof.* Our starting point is the CFL condition (21). Using Assumption 1, we have for any $\epsilon > 0$ and $\mathbf{v}_h \in \mathbf{V}_h$ the following inequality

$$ (2\delta^{-1}a_{h,\tau\nu} + (1-4\theta)\delta^{-2}a_{h,\nu\nu})(\mathbf{v}_h, \mathbf{v}_h) \le (\epsilon^{-1}\delta^{-1}a_{h,\tau\tau} + (\epsilon\delta^{-1} + (1-4\theta))\delta^{-2}a_{h,\nu\nu})(\mathbf{v}_h, \mathbf{v}_h). $$

Hence, choosing $\epsilon = \delta^{-1}(4\theta - 1)$ in the previous expression – which is strictly positive upon the assumption that $\theta > \frac{1}{4}$ – we have

$$ (a_{h,\tau\tau} + 2\delta^{-1}a_{h,\tau\nu} + (1-4\theta)\delta^{-2}a_{h,\nu\nu})(\mathbf{v}_h, \mathbf{v}_h) \le (1 + \frac{1}{4\theta - 1})a_{h,\tau\tau}(\mathbf{v}_h, \mathbf{v}_h). $$

We obtain the desired result using (21).

<figure><img src="image_3.png" /></figure>
### 3.3. Conditionning analysis of the implicit scheme

At each time step of the implicit scheme, a linear system of the following form needs to be solved

$$ \forall \mathbf{v}_h \in \mathbf{V}_h, \quad a_{h,\mathrm{im}}^{\delta}(\mathbf{u}_h^{n+1}, \mathbf{v}_h) = \tilde{\ell}_{h}^{n,n-1}(\mathbf{v}_h), \qquad (23) $$

where

$$ a_{h,\mathrm{im}}^{\delta}(\cdot, \cdot) = (m_h + \theta \Delta t^2 a_h^{\delta})(\cdot, \cdot), $$

and the right hand side depends on the previous steps as

$$ \tilde{\ell}_{h}^{n,n-1}(\mathbf{v}_h) = \Delta t^2 \ell(t^n; \mathbf{v}_h) + (2m_h - \Delta t^2(1-2\theta)a_h^{\delta})(\mathbf{u}_h^n, \mathbf{v}_h) - a_{h,\mathrm{im}}^{\delta}(\mathbf{u}_h^{n-1}, \mathbf{v}_h). $$

In the context of three dimensional high-frequency wave propagation, a strategy based upon a direct solver for computing the iterates is usually not recommended due to large memory requirements. Alternate approaches rely on iterative solvers, more specifically on Conjugate Gradient (CG) iterations since the bilinear form $a_{h,\mathrm{im}}^{\delta}(\cdot, \cdot)$ is symmetric and coercive. However, the amount of iterations required to obtain a precise solution is known to be sensitive to the conditioning of the underlying matrix. Unfortunately, the condition number is deteriorating as the plate thickness decreases, even when resorting to simple – and commonly accepted – preconditioning strategy based upon the mass operator. In the next result, we show that using the bilinear form

$$ a_{h,\mathrm{imex}}^{\delta}(\cdot, \cdot) = (m_h + \theta \Delta t^2 \delta^{-2} a_{h,\nu\nu})(\cdot, \cdot) $$

as pre-conditioner, we obtain a condition number independent of $\delta$. Additionally, this pre-conditioning strategy is computationally efficient since the bilinear form $a_{h,\nu\nu}(\cdot, \cdot)$ is partially lumped, as detailed later in Lemma 4.

**Theorem 2.** We denote by $\chi^*$ the condition number of the linear system (23) pre-conditioned by the operator $a_{h,\mathrm{imex}}^{\delta}(\cdot, \cdot)$, namely

$$ \chi^* = \left( \inf_{\mathbf{v}_h \in \mathbf{V}_h} \frac{a_{h,\mathrm{im}}^{\delta}(\mathbf{v}_h, \mathbf{v}_h)}{a_{h,\mathrm{imex}}^{\delta}(\mathbf{v}_h, \mathbf{v}_h)} \right)^{-1} \left( \sup_{\mathbf{v}_h \in \mathbf{V}_h} \frac{a_{h,\mathrm{im}}^{\delta}(\mathbf{v}_h, \mathbf{v}_h)}{a_{h,\mathrm{imex}}^{\delta}(\mathbf{v}_h, \mathbf{v}_h)} \right). \qquad (24) $$

If Assumption 1 is satisfied and the time step of the implicit scheme satisfies the following condition

$$ \alpha = \Delta t^2 \theta \sup_{\mathbf{v}_h \in \mathbf{V}_h} \frac{a_{h,\tau\tau}(\mathbf{v}_h, \mathbf{v}_h)}{m_h(\mathbf{v}_h, \mathbf{v}_h)} < 1, $$

then we have

$$ \chi^* < \frac{(\sqrt{\alpha} + 1)^2}{(\sqrt{\alpha} - 1)^2}. $$

*Proof.* We begin by estimating the second term (the supremum) in (24). We have, using Assumption 1,

$$ \sup_{\mathbf{v}_h \in \mathbf{V}_h} \frac{a_{h,\mathrm{im}}^{\delta}(\mathbf{v}_h, \mathbf{v}_h)}{a_{h,\mathrm{imex}}^{\delta}(\mathbf{v}_h, \mathbf{v}_h)} \le \sup_{\mathbf{v}_h \in \mathbf{V}_h} \frac{(m_h + \Delta t^2 \theta(1 + \epsilon^{-1})a_{h,\tau\tau} + \Delta t^2 \theta(1 + \epsilon)\delta^{-2}a_{h,\nu\nu})(\mathbf{v}_h, \mathbf{v}_h)}{a_{h,\mathrm{imex}}^{\delta}(\mathbf{v}_h, \mathbf{v}_h)}. $$

The positivity of the bilinear forms $m_h(\cdot, \cdot)$ and $a_{h,\nu\nu}(\cdot, \cdot)$ thus entails

$$ \sup_{\mathbf{v}_h \in \mathbf{V}_h} \frac{a_{h,\mathrm{im}}^{\delta}(\mathbf{v}_h, \mathbf{v}_h)}{a_{h,\mathrm{imex}}^{\delta}(\mathbf{v}_h, \mathbf{v}_h)} \le (1 + \epsilon) + (1 + \epsilon^{-1})\alpha. $$

We now come back to the estimation of the first term in (24). Once again from Assumption 1, and using the definition of the bilinear form $a_h^{\delta}(\cdot, \cdot)$, one has the following inequality

$$ a_h^{\delta}(\mathbf{v}_h, \mathbf{v}_h) \ge (1 - \frac{1}{\epsilon})a_{h,\tau\tau}(\mathbf{v}_h, \mathbf{v}_h) + (1 - \epsilon)\delta^{-2}a_{h,\nu\nu}(\mathbf{v}_h, \mathbf{v}_h). $$
From this estimate we can deduce that for any $\epsilon' > 0$

$$a_{h,im}^{\delta}(\mathbf{v}_h, \mathbf{v}_h) \geq ((1 - \epsilon')a_{h,imex}^{\delta} + \epsilon'm_h + \Delta t^2 \theta(1 - \frac{1}{\epsilon'})a_{h,\tau\tau})(\mathbf{v}_h, \mathbf{v}_h).$$

Therefore

$$\inf_{\mathbf{v}_h \in \mathbf{V}_h} \frac{a_{h,im}^{\delta}(\mathbf{v}_h, \mathbf{v}_h)}{a_{h,imex}^{\delta}(\mathbf{v}_h, \mathbf{v}_h)} \geq (1 - \epsilon') + \inf_{\mathbf{v}_h \in \mathbf{V}_h} \frac{(\epsilon' m_h + \Delta t^2 \theta(1 - \epsilon'^{-1}) a_{h,\tau\tau})(\mathbf{v}_h, \mathbf{v}_h)}{a_{h,imex}^{\delta}(\mathbf{v}_h, \mathbf{v}_h)}.$$

Choosing $\epsilon' < 1$, from the positivity of the bilinear form $\epsilon'm_h(\cdot, \cdot)$ we obtain

$$\inf_{\mathbf{v}_h \in \mathbf{V}_h} \frac{(\epsilon' m_h + \Delta t^2 \theta(1 - \epsilon'^{-1}) a_{h,\tau\tau})(\mathbf{v}_h, \mathbf{v}_h)}{a_{h,imex}^{\delta}(\mathbf{v}_h, \mathbf{v}_h)} \geq (1 - \frac{1}{\epsilon'}) \Delta t^2 \theta \inf_{\mathbf{v}_h \in \mathbf{V}_h} \frac{a_{h,\tau\tau}(\mathbf{v}_h, \mathbf{v}_h)}{a_{h,imex}^{\delta}(\mathbf{v}_h, \mathbf{v}_h)} \geq (1 - \frac{1}{\epsilon'}) \alpha,$$

which finally leads to the following estimate for the infimum

$$\inf_{\mathbf{v}_h \in \mathbf{V}_h} \frac{a_{h,im}^{\delta}(\mathbf{v}_h, \mathbf{v}_h)}{a_{h,imex}^{\delta}(\mathbf{v}_h, \mathbf{v}_h)} \geq (1 - \epsilon') + (1 - \frac{1}{\epsilon'}) \alpha.$$

By assumption $\alpha < 1$, hence one can choose $\epsilon = \sqrt{\alpha}$ in the upper bound of the supremum and $\epsilon' = \sqrt{\alpha}$ in the lower bound of the infimum to obtain the desired result on the condition number. $\square$

**Remark 2.** Choosing $\theta = \frac{1}{4}$, the condition $\alpha < 1$ in the previous Theorem is equivalent to

$$\Delta t < 2 \left( \sup_{\mathbf{v}_h \in \mathbf{V}_h} \frac{a_{h,\tau\tau}(\mathbf{v}_h, \mathbf{v}_h)}{m_h(\mathbf{v}_h, \mathbf{v}_h)} \right)^{-\frac{1}{2}}, \qquad (25)$$

which is the CFL condition for the fully explicit scheme involving only $a_{h,\tau\tau}$. Hence the obtained estimate is valid if the time step satisfies a CFL-like condition independent of $\delta$. The preconditioning is therefore efficient. Note that the scheme is stable for any choice of $\Delta t$, it is however unclear how the conditionning of the underlying system behave when (25) is not satisfied.

## 4. Application to space discretization with high order spectral finite elements

In this section we detail a specific construction of the approximation space $\mathbf{V}_h$ based upon high order spectral finite elements [10, 11, 12]. On the one hand, we show that using an extruded mesh of the domain $\Omega$ one recovers some sort of lumping – referred to as partial lumping – of the matrix associated to the bilinear form $a_{h,imex}^{\delta}(\cdot, \cdot)$. This is a key point to ensure the efficiency of the implicit-explicit scheme. On the other hand, we estimate the CFL conditions of the schemes presented in the previous section in order to make explicit their dependencies w.r.t discretization parameters and the thickness of the plate.

### 4.1. Quasi-uniform and extruded mesh of the domain $\Omega$

Relying on the nature of the geometry $\Omega$, we assume in the following that the mesh is extruded and satisfies additional quasi-uniformity assumptions – see Figure 1 for an illustration.
Quasi-uniform surface mesh of $\Sigma$. Let $h > 0$ be the (planar) space discretization parameter. We consider $\{\mathcal{T}_h(\bar{\Sigma})\}_{h>0}$ a family of meshes of the surface $\Sigma$. The mesh step $h$ tends to zero as the surface mesh is refined, and we provide a precise definition of this parameter later on. For a fixed $h$, $\mathcal{T}_h(\bar{\Sigma})$ is a set of $N_h$ quadrangles generating a tessellation of $\bar{\Sigma}$ in the standard sense, without hanging nodes. For any $\ell \in \llbracket 1; N_h \rrbracket$ we denote by $\kappa_\ell \in \mathcal{T}_h(\bar{\Sigma})$ a quadrangle of the mesh. Each quadrangle $\kappa_\ell$ is uniquely defined by its four points $\{\mathbf{s}_j^\ell\}_{j=0}^3 \subset \mathbb{R}^2$, and we denote by $\mathbf{f}_\ell$ the unique transformation of the reference quadrangle $\hat{\kappa} = (0; 1)^2$ such that $\mathbf{f}_\ell(\hat{\kappa}) = \kappa_\ell$.

Introducing some usual definitions – see e.g. in Appendix A of [27] – we denote by $h_\ell = \text{diam}(\kappa_\ell)$ the maximal euclidean distance between two points in $\{\mathbf{s}_j^\ell\}_{j=0}^3$. The space discretization parameter associated to the mesh is then defined as

$$h = \max_{\ell \in \llbracket 1; N_h \rrbracket} \{h_\ell\}.$$

For any $j \in \mathbb{Z}/4\mathbb{Z}$, the triangle built from the points $\mathbf{s}_{j-1}^\ell$, $\mathbf{s}_j^\ell$ and $\mathbf{s}_{j+1}^\ell$ is denoted by $T_j^\ell$, and the diameter of the corresponding inscribed circle by $\rho_j^\ell$. Defining

$$\rho_\ell = 2 \min_{j \in \llbracket 0; 3 \rrbracket} \{\rho_j^\ell\},$$

the aspect ratio of a quadrangle element is defined as

$$\sigma_\ell = \frac{h_\ell}{\rho_\ell},$$

which is equal to one for square-shaped quadrangle and grows larger as the element flatten. With these definitions, we now give a precise statement of the quasi-uniformity assumption satisfied by the family of meshes of the surface $\Sigma$.

**Assumption 2.** We assume that there exists two constants $\sigma > 0$ and $\tau > 0$ (that are independent of the mesh step) such that for any mesh in $\{\mathcal{T}_h(\bar{\Sigma})\}_{h>0}$ we have

$$\tau h \le h_\ell \le \sigma \rho_\ell, \quad \forall \ell \in \llbracket 1; N_h \rrbracket. \qquad (26)$$

Note that the left-hand side of (26) entails in particular that $\tau h \le h_\ell \le h$, which imposes a mesh composed of elements with similar diameters. Whereas the right-hand side of (26) leads to a bounded $\sigma_\ell$, in other words the aspect ratio of each elements in any mesh of the surface does not degenerates.

Associated to each transformation of the elements, one can provide estimates for Jacobian matrix and its determinant. To do so, let us introduce the standard euclidean norm $\| \cdot \|_d$ in $\mathbb{R}^d$, with $d = 2$ or $d = 3$, and for any domain $D \in \mathbb{R}^d$ the operator norms

$$\left\| M \right\|_{\infty, D} = \sup_{\mathbf{x} \in D} \left\| M(\mathbf{x}) \right\|, \quad \left\| M(\mathbf{x}) \right\| = \sup_{\mathbf{v} \in \mathbb{R}^d \setminus \{0\}} \frac{\| M(\mathbf{x})\mathbf{v} \|_d}{\| \mathbf{v} \|_d}, \quad \forall M \in \mathcal{C}^0(D, \mathbb{R}^{d \times d}), \\ \left| \alpha \right|_{\infty, D} = \sup_{\mathbf{x} \in D} \left| \alpha(\mathbf{x}) \right|, \quad \forall \alpha \in \mathcal{C}^0(D, \mathbb{R}).$$

With these notations, we give without demonstration the following Lemma – see e.g. [27] or [28] for more details. Note that in what follows we use the notation $\lesssim$ to represent an inequality which is valid up to a constant independent of any element in the family of meshes.

**Lemma 2.** Denoting $j_\ell = \det \nabla \mathbf{f}_\ell$, we have

$$\| \nabla \mathbf{f}_\ell \|_{\infty, \hat{\kappa}} \lesssim h_\ell, \quad |j_\ell|_{\infty, \hat{\kappa}} \lesssim h_\ell^2, \quad \forall \ell \in \llbracket 1; N_h \rrbracket. \qquad (27)$$

For the inverse geometric transformation, similar results hold, namely

$$\| \nabla \mathbf{f}_\ell^{-1} \|_{\infty, \hat{\kappa}} \lesssim \frac{h_\ell}{\rho_\ell^2}, \quad |j_\ell|_{\infty, \hat{\kappa}} \lesssim \frac{1}{\rho_\ell^2}, \quad \forall \ell \in \llbracket 1; N_h \rrbracket. \qquad (28)$$
Quasi-uniform extruded mesh of $\Omega$. Let $\eta > 0$ be the space discretization parameter in the orthogonal direction of the surface plane. We consider a subdivision of $N_\eta$ edge elements of the interval $(-1, 1)$ such that

$$1 = z_{N_{\eta}+1} > z_{N_{\eta}} > \cdots > z_1 = -1, \quad \eta = \max_{m \in \llbracket 1; N_{\eta} \rrbracket} \{z_{m+1} - z_m\}.$$

As for the space discretization of the surface $\Sigma$ we make the assumption of quasi-uniformity of the subdivision of the interval.

**Assumption 3.** We assume that there exists a constant $\mu > 0$ (independent of $\eta$) such that

$$\mu\eta \le z_{m+1} - z_m \le \eta, \quad \forall m \in \llbracket 1; N_{\eta} \rrbracket. \qquad (29)$$

Denoting by $\mathfrak{h} = (h, \eta)$ the couple of space discretization parameters of the surface $\Sigma$ and the interval $(-1, 1)$, we consider a family of space discretizations of the complete domain $\{\mathcal{T}_{\mathfrak{h}}(\bar{\Omega})\}_{\mathfrak{h}}$ which are composed of hexahedra built by extrusion of the quadrangles of the surface mesh. Namely, for two fixed mesh parameters $h$ and $\eta$, and for any $\ell \in \llbracket 1; N_h \rrbracket$ and $m \in \llbracket 1; N_\eta \rrbracket$ we denote by $K_{\ell,m} \in \mathcal{T}_{\mathfrak{h}}(\bar{\Omega})$ an hexahedron of the mesh such that

$$K_{\ell,m} = \kappa_{\ell} \times (z_m, z_{m+1}).$$

To each hexahedron we associate a unique transformation $\mathbf{F}_{\ell,m}$ of the reference hexahedron $\hat{\mathbf{K}} = (0, 1)^3$ such that $\mathbf{F}_{\ell,m}(\hat{\mathbf{K}}) = K_{\ell,m}$ and defined as

$$\mathbf{F}_{\ell,m}(\hat{x}_1, \hat{x}_2, \hat{x}_3) = \begin{pmatrix} \mathbf{f}_{\ell}(\hat{x}_1, \hat{x}_2) \\ (1 - \hat{x}_3)z_m + \hat{x}_3 z_{m+1} \end{pmatrix}, \quad \forall (\hat{x}_1, \hat{x}_2, \hat{x}_3)^{\top} \in \hat{\mathbf{K}},$$

from which we can derive the following expression of the Jacobian matrix and its determinant

$$\nabla \mathbf{F}_{\ell,m} = \begin{pmatrix} \nabla \mathbf{f}_{\ell} & 0 \\ 0 & z_{m+1} - z_m \end{pmatrix}, \quad J_{\ell,m} = \det \nabla \mathbf{F}_{\ell,m} = j_{\ell}(z_{m+1} - z_m). \qquad (30)$$

Note that, using this expression of the determinant and the second estimates in (27) and (28), we have

$$|J_{\ell,m}|_{\infty, \hat{\mathbf{K}}} \lesssim h_{\ell}^2 \eta, \quad |J_{\ell,m}^{-1}|_{\infty, \hat{\mathbf{K}}} \lesssim \frac{1}{\rho_{\ell}^2 \mu \eta}. \qquad (31)$$

## 4.2. High-order spectral finite elements

Let us define $p_\tau$ and $p_\nu$ two natural integers representing the order of approximation in the tangent plane and in the orthogonal direction respectively, and $\hat{n}_\tau = p_\tau + 1$ and $\hat{n}_\nu = p_\nu + 1$. We define on the reference hexahedron the following polynomial space

$$Q_{p_{\tau}, p_{\nu}}(\hat{\mathbf{K}}) = \operatorname{span} \{ \hat{x}_1^i \hat{x}_2^j \hat{x}_3^k, \quad \forall (i, j, k) \in \mathbb{N} : \max\{i, j\} \le p_{\tau}, k \le p_{\nu} \}.$$

The dimension of this local polynomial space is denoted by $\hat{n} = \hat{n}_{\tau}^2 \hat{n}_{\nu}$. From this local space, we can build the global conforming approximation space that supports our semi-discretization procedure, namely,

$$V_{\mathfrak{h}} = \{ v_{\mathfrak{h}} \in \mathcal{C}^0(\bar{\Omega}) \mid \forall K_{\ell,m} \in \mathcal{T}_{\mathfrak{h}}(\bar{\Omega}), \exists \hat{v} \in Q_{p_{\tau},p_{\nu}}(\hat{\mathbf{K}}), v_{\mathfrak{h}}|_{K_{\ell,m}} = \hat{v} \circ \mathbf{F}_{\ell,m}^{-1} \}, \quad \mathbf{V}_{\mathfrak{h}} = V_{\mathfrak{h}}^3 \subset \mathbf{V}.$$

In our work, we rely on so-called spectral finite elements [10, 11, 12]. It is a high-order conform finite element methods relying on Gauss-Lobatto points in order to avoid Runge's phenomenon and to render spectral-like convergence. In the remaining part of this subsection, we recall some standard notations on the construction of the finite element basis functions of the approximation space. Let us define $\{\hat{\xi}_{\tau,i}\}_{i=1}^{\hat{n}_{\tau}}$ and $\{\hat{\xi}_{\nu,i}\}_{i=1}^{\hat{n}_{\nu}}$ the $\hat{n}_{\tau}$ and the $\hat{n}_{\nu}$ Gauss-Lobatto points on the reference edge $\hat{e} = (0, 1)$ (see e.g. [12] for an explicit definition and related examples of these points). From these points, we associate the one-dimensional Lagrange polynomials

$$\{\hat{\varphi}_{\tau,i}\}_{i=1}^{\hat{n}_{\tau}} \subset \mathcal{P}_{p_{\tau}}(\hat{e}) \quad \text{and} \quad \{\hat{\varphi}_{\nu,i}\}_{i=1}^{\hat{n}_{\nu}} \subset \mathcal{P}_{p_{\nu}}(\hat{e}),$$
where $\mathcal{P}_k(\hat{e})$ is the space of polynomials of order lower or equal to $k$ on the reference edge. As Lagrange polynomials they satisfy, by definition,

$$ \hat{\varphi}_{\alpha,i}(\hat{\xi}_{\alpha,j}) = \delta_{ij}, \quad \forall i, j = 1, \dots, \hat{n}_{\alpha}, \quad \forall \alpha \in \{\tau, \nu\}. $$

Let us denote by $\hat{n} = [1; \hat{n}_{\tau}]^2 \times [1; \hat{n}_{\nu}]$, and by $i = (i_1, i_2, i_3)$ any multi-index such that $i \in \hat{n}$. Then, on the reference hexahedron $\hat{K}$, we construct the Gauss-Lobatto points $\{\hat{\xi}_i\}_{i \in \hat{n}}$ by tensorisation of the one-dimensional points defined on $\hat{e}$, namely

$$ \forall i \in \hat{n}, \quad \hat{\xi}_i = (\hat{\xi}_{\tau,i_1}, \hat{\xi}_{\tau,i_2}, \hat{\xi}_{\nu,i_3})^T \in \hat{K}. $$

In the same manner, we define the Lagrange polynomials $\{\hat{\varphi}_i\}_{i \in \hat{n}} \subset Q_{p_{\tau}, p_{\nu}}(\hat{K})$ on the reference element by

$$ \forall i \in \hat{n}, \quad \hat{\varphi}_i(\hat{x}) = \hat{\varphi}_{\tau,i_1}(\hat{x}_1)\hat{\varphi}_{\tau,i_2}(\hat{x}_2)\hat{\varphi}_{\nu,i_3}(\hat{x}_3), $$

which satisfy $\hat{\varphi}_i(\hat{\xi}_j) = \delta_{ij}$. The set of $n_h$ interpolation points is the largest set $\{\hat{\xi}_I\}_{I=1}^{n_h}$ of distinct points such that

$$ \{\hat{\xi}_I\}_{I=1}^{n_h} \subset \bigcup_{K_{\ell,m} \in \mathcal{T}_h(\bar{\Omega})} \bigcup_{i=1}^{\hat{n}} \mathbf{F}_{\ell,m}(\hat{\xi}_i) \subset \bar{\Omega}. $$

From this set of points, one can construct the so-called “local-to-global” mapping $\mathcal{LG}(i, \ell, m)$, which to a local point of index $i \in \hat{n}$ within an element of indexes $\ell \in [1; N_h]$ (numbering in the tangent plane) and $m \in [1; N_{\eta}]$ (numbering in the orthogonal direction) renders the associated global point of index $I \in [1; n_h]$. The definition of the global Lagrange basis functions $\{\varphi_I\}_{I=1}^{n_h}$ follows

$$ \forall I \in [1; n_h], \quad \forall K_{\ell,m} \in \mathcal{T}_h(\bar{\Omega}), \quad \varphi_I|_{K_{\ell,m}} = \begin{cases} \hat{\varphi}_i \circ \mathbf{F}_{\ell,m}^{-1}, & \text{if } I = \mathcal{LG}(i, \ell, m), \\ 0, & \text{otherwise.} \end{cases} $$

One can verify that they actually are a generating basis of the approximation space, i.e. $V_h = \text{span}\{\varphi_I\}_{I=1}^{n_h}$.

### 4.3. Numerical integration

We rely on a quadrature formula to numerically compute the integrals involved in the bilinear forms $m(\cdot, \cdot)$ and $a(\cdot, \cdot)$. The formula is defined on the reference hexahedron. We denote its associated weights and points respectively by

$$ (\{\hat{\omega}_i\}_{i \in \hat{n}}, \{\hat{\xi}_i\}_{i \in \hat{n}}). $$

For any function $f$ at least continuous over each elements in the mesh, we introduce the following notations to differentiate the exact integration form the numerical one

$$ \left| \begin{array}{l} \displaystyle \int_{K_{\ell,m}} f \, dK \approx \int_{K_{\ell,m}}^{Q} f \, dK = \int_{\hat{K}}^{Q} (f \circ \mathbf{F}_{\ell,m}) |J_{\ell,m}| \, d\hat{K} = \sum_{k=1}^{\hat{n}} \hat{\omega}_k f(\mathbf{F}_{\ell,m}(\hat{\xi}_k)) |J_{\ell,m}(\hat{\xi}_k)|, \\ \\ \displaystyle \int_{\Omega} f \, d\Omega \approx \int_{\Omega}^{Q} f \, d\Omega = \sum_{K_{\ell,m} \in \mathcal{T}_h(\bar{\Omega})} \int_{K_{\ell,m}}^{Q} f \, dK. \end{array} \right. $$

The quadrature points are the Gauss-Lobatto points and the quadrature weights are the strictly positive Gauss-Lobatto weights. Matching the quadrature points with the points associated to the local Lagrange basis functions is referred to as “mass lumping” [10, 11, 12]. It leads to a diagonal approximated mass matrix. It has been the subject of particular interest in the context of finite element methods applied to time domain wave propagation problems because it enables fully explicit time discretization schemes. In order to conserve stability and accuracy, specific conditions must be enforced on the quadrature points and weights. For quadrangles and hexahedra, these conditions are met at any order of approximations by the
Gauss-Lobatto rules [12, 18], while for triangles and tetrahedra dedicated construction of the quadrature formula is needed [26, 19, 20]. In the following $m_h(\cdot, \cdot)$ and $a_h^\delta(\cdot, \cdot)$ introduced in (10) are the approximated bilinear forms where the integral over elements are evaluated using this quadrature formula, namely

$$ \left\| \begin{array}{l} m(\mathbf{w}_h, \mathbf{v}_h) \approx m_h(\mathbf{w}_h, \mathbf{v}_h) = \int_{\Omega}^{Q} \rho \mathbf{w}_h \cdot \mathbf{v}_h \, d\Omega, \\ a^\delta(\mathbf{w}_h, \mathbf{v}_h) \approx a_h^\delta(\mathbf{w}_h, \mathbf{v}_h) = \int_{\Omega}^{Q} \mathbf{C} \varepsilon^\delta(\mathbf{w}_h) : \varepsilon^\delta(\mathbf{v}_h) \, d\Omega, \end{array} \right. \quad \forall \mathbf{v}_h, \mathbf{w}_h \in \mathbf{V}_h. $$

From this specific choice of quadrature formula, we verify that Assumption 1 holds, which is a direct consequence of the positivity of the quadrature weights.

**Lemma 3.** For any $\epsilon > 0$ and for any $\mathbf{v}_h \in \mathbf{V}_h$, we have

$$ 2a_{h,\tau\nu}(\mathbf{v}_h, \mathbf{v}_h) \leq \frac{1}{\epsilon}a_{h,\tau\tau}(\mathbf{v}_h, \mathbf{v}_h) + \epsilon a_{h,\nu\nu}(\mathbf{v}_h, \mathbf{v}_h). $$

*Proof.* Upon the equation (6) of uniform coercivity and continuity on the fourth order tensor $\mathbf{C}$ we can insure that for all $\mathbf{x} \in \bar{\Omega}$, $\mathbf{C}(\mathbf{x})$ defines an equivalent inner product on $L_s(\mathbb{R}^3)$. Hence, using Cauchy-Schwarz and Young inequalities enables us to verify that for any $\epsilon > 0$ we have

$$ \forall \varepsilon, \sigma \in \mathcal{L}_s(\mathbb{R}^3), \quad 2|\mathbf{C}(\mathbf{x})\varepsilon : \sigma| \leq \frac{1}{\epsilon}(\mathbf{C}(\mathbf{x})\varepsilon : \varepsilon) + \epsilon(\mathbf{C}(\mathbf{x})\sigma : \sigma). $$

Using the fact that all the quadrature weights are strictly positive, we obtain the result of the lemma by numerically integrating – i.e. using the quadrature formula – over each elements and summing on the elements the previous expression with $\varepsilon = \varepsilon_\tau(\mathbf{v}_h)$ and $\sigma = \varepsilon_\nu(\mathbf{v}_h)$. $\square$

### 4.4. Partial lumping of the discrete stiffness bilinear form $a_{h,\nu\nu}$

As already mentioned the matrix corresponding to the bilinear form $m_h(\cdot, \cdot)$ is diagonal. In the remaining part of this section we show that the matrix corresponding to the bilinear form $a_{h,\nu\nu}(\cdot, \cdot)$ is block diagonal where only degrees of freedom corresponding to interpolation points with the same tangent coordinates are coupled. To show this we only study the simpler scalar case as the extension to the vectorial case is straightforward. This partial lumping has an important consequence since the matrix corresponding to the bilinear form $a_{h,\text{imex}}^\delta(\cdot, \cdot)$ inherits this block diagonal property. It is therefore easy to invert by direct solver and in parallel for each set of degrees of freedom with the same coordinates in the tangent plane.

**Lemma 4.** Let $K_{\ell,m} \in \mathcal{T}_h(\bar{\Omega})$ and $\gamma \in C^0(K_{\ell,m})$ be any continuous scalar-valued function, for all $(\varphi_I, \varphi_J) \in \mathbf{V}_h \times \mathbf{V}_h$ non-vanishing over $K_{\ell,m}$

$$ \int_{K_{\ell,m}}^{Q} \gamma \partial_3 \varphi_I \partial_3 \varphi_J \, d\Omega = \frac{\delta_{i_1 j_1} \delta_{i_2 j_2}}{(z_{m+1} - z_m)^2} \left( \sum_{k_3=1}^{\hat{n}_\nu} \hat{\omega}_k \gamma(\mathbf{F}_{\ell,m}(\hat{\boldsymbol{\xi}}_k)) |J_{\ell,m}(\hat{\boldsymbol{\xi}}_k)| \hat{\varphi}'_{\nu, i_3}(\hat{\xi}_{\nu, k_3}) \hat{\varphi}'_{\nu, j_3}(\hat{\xi}_{\nu, k_3}) \right), \quad (32) $$

where $I = \mathcal{LG}(i, \ell, m)$, $J = \mathcal{LG}(j, \ell, m)$ and $k = (i_1, i_2, k_3)$.

*Proof.* By applying the change variable to the reference element and the quadrature formula the approximated integral is given by

$$ (z_{m+1} - z_m)^{-2} \sum_{k \in \hat{n}} \hat{\omega}_k \gamma(\mathbf{F}_{\ell,m}(\hat{\boldsymbol{\xi}}_k)) |J_{\ell,m}(\hat{\boldsymbol{\xi}}_k)| \partial_3 \hat{\varphi}_i(\hat{\boldsymbol{\xi}}_k) \partial_3 \hat{\varphi}_j(\hat{\boldsymbol{\xi}}_k). $$
From the tensorized construction of the local Lagrange polynomial, their derivative in the third direction evaluated at the quadrature points reads

$$ \partial_3 \hat{\varphi}_i(\hat{\xi}_k) = \hat{\varphi}_{\tau,i_1}(\hat{\xi}_{\tau,k_1}) \hat{\varphi}_{\tau,i_2}(\hat{\xi}_{\tau,k_2}) \hat{\varphi}'_{\nu,i_3}(\hat{\xi}_{\nu,k_3}) = \hat{\varphi}'_{\nu,i_3}(\hat{\xi}_{\nu,k_3}) \delta_{i_1 k_1} \delta_{i_2 k_2}, $$

where $i = (i_1, i_2, i_3)$ and $k = (k_1, k_2, k_3)$. Proceeding in the same fashion for $\partial_3 \hat{\varphi}_j(\hat{\xi}_k)$ with $j = (j_1, j_2, j_3)$ leads to the desired result. $\square$

This lemma shows that – locally to each element $K_{\ell,m}$ – the only test functions for which the integral

$$ \int_{K_{\ell,m}}^{Q} \gamma \partial_3 \varphi_I \partial_3 \varphi_J d\Omega $$

does not vanish correspond to degrees of freedom sharing the same indexes in the tangent plane (in the local numbering). As a direct consequence this property translates to

$$ \int_{\Omega}^{Q} \gamma \partial_3 \varphi_I \partial_3 \varphi_J d\Omega. $$

**Remark 3.** Lemma 4 is of major importance for the practical application of the presented implicit-explicit scheme. In essence, it states that the linear system associated to the bilinear form $(m_h + \theta \Delta t^2 \delta^{-2} a_{\nu\nu,h})(\cdot, \cdot)$ and to be solved at each time step is block diagonal. Each blocks are symmetric positive definite matrices, their dimension is equal to the number of degrees of freedom in the thickness and they are independent from one degrees of freedom in the tangent plane to another. In practice – and in the numerical illustrations pretested in Section 5 – we compute the Cholesky decomposition of these blocks in a initializing phase and we apply the forward and backward substitutions within a parallel loop over the degrees of freedom of the tangent plane. By doing so, we limit drastically the cost of solving this linear system, thus benefiting the most from the gain on the time step brought by the implicit-explicit scheme.

### 4.5. Estimates on the stability conditions

Estimates on Rayleigh quotients w.r.t the discretization parameters and the plate thickness. In this subsection, we gather estimates on different Rayleigh quotients. These estimates will provide tools to analyzing in the next sections the dependency of the stability conditions of the explicit and implicit-explicit schemes w.r.t the discretization parameters and the plate thickness.

**Lemma 5.** There exists two constants $c_\tau > 0$ and $C_\tau > 0$ independent of $h$ such that

$$ c_{\tau}h^{-2} \leq \sup_{\mathbf{v}_h \in \mathbf{V}_h} \frac{a_{h,\tau\tau}(\mathbf{v}_h, \mathbf{v}_h)}{m_h(\mathbf{v}_h, \mathbf{v}_h)} \leq C_{\tau}h^{-2}. \quad (33) $$

Furthermore, there exists two constants $\hat{c} > 0$ and $\hat{C} > 0$ that depends only on the reference element and the local polynomial space $Q_{p_\tau, p_\nu}(\hat{K})$ such that

$$ c_{\tau} = \hat{c} \frac{\alpha}{\varrho_{+}} \frac{\mu}{\sigma^2}, \quad C_{\tau} = \hat{C} \frac{\beta}{\varrho_{-}} \frac{\sigma^6}{\mu \tau^2}. $$

*Proof.* In the following, we denote by $K$ any element in $\mathcal{T}_h(\bar{\Omega})$ of indexes $\ell \in [1, N_h]$ and $m \in [1, N_\eta]$. We decompose the proof by starting with the upper bound and finishing with the lower bound in (33).
Upper bound in (33). We introduce the (scalar and vectorial) spaces of square integrable functions that are polynomials on each elements, without any continuity constraint, namely

$$V_h^b = \{v_h \in L^2(\Omega) \mid \forall K \in \mathcal{T}_h(\bar{\Omega}), \exists \hat{v} \in Q_{p_\tau, p_\nu}(\hat{K}), v_h|_K = \hat{v} \circ \mathbf{F}_{\ell,m}^{-1}\}, \quad \mathbf{V}_h^b = V_h^b \times V_h^b \times V_h^b.$$

Also, let

$$\forall \mathbf{v}_h, \mathbf{w}_h \in \mathbf{V}_h^b, \quad m_h^K(\mathbf{v}_h, \mathbf{w}_h) = \int_K^Q \varrho \mathbf{v}_h \cdot \mathbf{w}_h \, dK, \quad a_h^K_{\mathbf{v}_h, \tau\tau}(\mathbf{w}_h, \mathbf{v}_h) = \int_K^Q \mathbf{C} \varepsilon_\tau(\mathbf{w}_h) : \varepsilon_\tau(\mathbf{v}_h) \, dK,$$

be the local expression of the mass and stiffness bilinear forms $m_h(\cdot, \cdot)$ and $a_h, \tau\tau(\cdot, \cdot)$. By construction we have $\mathbf{V}_h \subset \mathbf{V}_h^b$, implying in particular that

$$\sup_{\mathbf{v}_h \in \mathbf{V}_h} \frac{a_h, \tau\tau(\mathbf{v}_h, \mathbf{v}_h)}{m_h(\mathbf{v}_h, \mathbf{v}_h)} \le \sup_{\mathbf{v}_h \in \mathbf{V}_h^b} \frac{a_h, \tau\tau(\mathbf{v}_h, \mathbf{v}_h)}{m_h(\mathbf{v}_h, \mathbf{v}_h)} \le \max_{K \in \mathcal{T}_h(\bar{\Omega})} \left( \sup_{\mathbf{v}_h \in \mathbf{V}_h^b} \frac{a_{h, \tau\tau}^K(\mathbf{v}_h, \mathbf{v}_h)}{m_h^K(\mathbf{v}_h, \mathbf{v}_h)} \right), \quad (34)$$

where, for the second inequality, we used the decomposition of the global linear forms over each elements. We begin by finding an upper bound on the local stiffness term. Using successively the triangle inequality, the Cauchy-Schwarz inequality for the inner product in $L_s(\mathbb{R}^3)$ and (6), we obtain

$$\forall \mathbf{v}_h \in \mathbf{V}_h^b, \quad a_{h, \tau\tau}^K(\mathbf{v}_h, \mathbf{v}_h) \le \beta \int_K^Q |\varepsilon_\tau(\mathbf{v}_h)|^2 \, dK.$$

From the definition (8) of the tangential Green-Lagrange tensor, one can verify that

$$\forall \mathbf{v}_h = (v_{h,1}, v_{h,2}, v_{h,3})^\intercal \in \mathbf{V}_h^b, \quad |\varepsilon_\tau(\mathbf{v}_h)|^2 \le \sum_{i=1}^3 \|\mathbf{A}_\tau \nabla v_{h,i}\|_3^2,$$

with $\mathbf{A}_\tau = \mathbf{I} - \mathbf{e}_3 \otimes \mathbf{e}_3$. Applying the change of variable to the reference element yields

$$\forall \mathbf{v}_h \in \mathbf{V}_h^b, \quad a_{h, \tau\tau}^K(\mathbf{v}_h, \mathbf{v}_h) \le \beta \sum_{i=1}^3 \int_{\hat{K}}^Q \|\mathbf{A}_\tau \nabla F_{\ell,m}^{-\intercal} \nabla \hat{v}_i\|_3^2 |J_{\ell,m}| \, d\hat{K},$$

where $\hat{v}_i = (v_{h,i} \circ \mathbf{F}_{\ell,m}) \in Q_{p_\tau, p_\nu}(\hat{K})$. Remarking that

$$\mathbf{A}_\tau \nabla \mathbf{F}_{\ell,m}^{-\intercal} = \begin{pmatrix} \nabla \mathbf{f}_\ell^{-\intercal} & 0 \\ 0 & 0 \end{pmatrix},$$

and using estimates (28) and (31), we can insure that there exists a constant $C > 0$ independent of the element such that

$$\forall \mathbf{v}_h \in \mathbf{V}_h^b, \quad a_{h, \tau\tau}^K(\mathbf{v}_h, \mathbf{v}_h) \le C\beta (h_\ell^2 \eta) \left(\frac{h_\ell}{\rho_\ell^2}\right)^2 \int_{\hat{K}}^Q |\nabla \hat{v}|^2 \, d\hat{K},$$

where $\hat{v} = (\hat{v}_1, \hat{v}_2, \hat{v}_3)^\intercal$. Upon the Assumption 2 of quasi-uniformity of the mesh of the surface $\Sigma$, this boils down to

$$\forall \mathbf{v}_h \in \mathbf{V}_h^b, \quad a_{h, \tau\tau}^K(\mathbf{v}_h, \mathbf{v}_h) \le C\beta\eta\sigma^4 \int_{\hat{K}}^Q |\nabla \hat{v}|^2 \, d\hat{K}. \quad (35)$$

Now moving on to the estimate of the local mass term in (34). Using the bounds (5) on the mass density and the last estimate in (31), there exists a constant $C' > 0$ independent of the element such that for any $\mathbf{v}_h \in \mathbf{V}_h^b$

$$m_h^K(\mathbf{v}_h, \mathbf{w}_h) \ge \varrho_- \int_K^Q \|\mathbf{v}_h\|_3^2 dK \ge \varrho_- \inf_{\hat{\mathbf{x}} \in \hat{K}} |J_{\ell,m}(\hat{\mathbf{x}})| \int_{\hat{K}}^Q \|\hat{\mathbf{v}}\|_3^2 d\hat{K} \ge C' \varrho_- (\rho_\ell^2 \mu\eta) \int_{\hat{K}}^Q \|\hat{\mathbf{v}}\|_3^2 d\hat{K}.$$
Using once again the quasi-uniformity Assumption 2, we obtain

$$ \forall \mathbf{v}_h \in \mathbf{V}_h^b, \quad m_h^K(\mathbf{v}_h, \mathbf{v}_h) \ge C' \rho_-(\frac{\tau^2 h^2}{\sigma^2} \mu \eta) \int_{\hat{K}}^{Q} \|\hat{\mathbf{v}}\|_3^2 d\hat{K}. \quad (36) $$

Combining (35) and (36) we obtain the upper bound, for any $\mathbf{v}_h \in \mathbf{V}_h^b$ which cleary depends exclusively on the reference element and the local polynomial space.

Lower bound in (33). Once again, using the bounds (5) on the mass density and on the coercivity of C in (6), we have

$$ \sup_{\mathbf{v}_h \in \mathbf{V}_h} \frac{a_{h, \tau}(\mathbf{v}_h, \mathbf{v}_h)}{m_h(\mathbf{v}_h, \mathbf{v}_h)} \ge \left(\frac{\alpha}{\rho_+}\right) \sup_{\mathbf{v}_h \in \mathbf{V}_h} \left( \int_{\Omega}^{Q} |\varepsilon_{\tau}(\mathbf{v}_h)|^2 d\Omega \right) \left( \int_{\Omega}^{Q} \|\mathbf{v}_h\|_3^2 d\Omega \right)^{-1}. $$

In order to provide a lower bound for this supremum we choose a relevant test function. Namely, for any global index $I \in [1, n_h]$, we define $\mathbf{v}_h^* = (0, 0, \varphi_I)^T$. Starting with the stiffness term, for any $K_{\ell,m} \in T_h(\Omega)$ such that $K_{\ell,m} \subset \text{supp}\{\varphi_I\}$ we have

$$ \int_{\Omega}^{Q} |\varepsilon_{\tau}(\mathbf{v}_h^*)|^2 d\Omega = \frac{1}{2} \int_{\Omega}^{Q} \|\mathbf{A}_{\tau} \nabla \varphi_I\|_3^2 d\Omega \ge \frac{1}{2} \int_{K_{\ell,m}}^{Q} \|\mathbf{A}_{\tau} \nabla \varphi_I\|_3^2 dK. $$

Performing the change of variable to the reference element leads to

$$ \int_{\Omega}^{Q} |\varepsilon_{\tau}(\mathbf{v}_h^*)|^2 d\Omega \ge \frac{1}{2} \inf_{\hat{\mathbf{x}} \in \hat{K}} |J_{\ell,m}(\hat{\mathbf{x}})| \|\nabla \mathbf{f}_{\ell}\|_{\infty, \hat{K}}^{-2} \int_{\hat{K}}^{Q} \|\nabla \hat{\varphi}_i\|_3^2 d\hat{K}, $$

where $i \in \hat{n}$ is a local index such that $I = \mathcal{LG}(i, \ell, m)$. Using the last estimate in (31), the first estimate in (27), and upon the quasi-uniformity Assumption 2, one can verify that there exists a constant $c > 0$ independent of the element such that

$$ \int_{\Omega}^{Q} |\varepsilon_{\tau}(\mathbf{v}_h^*)|^2 d\Omega \ge c (\rho_{\ell}^2 \mu \eta) h_{\ell}^{-2} \int_{\hat{K}}^{Q} \|\nabla \hat{\varphi}_i\|_3^2 d\hat{K} \ge c \left(\frac{\mu \eta}{\sigma^2}\right) \int_{\hat{K}}^{Q} \|\nabla \hat{\varphi}_i\|_3^2 d\hat{K}. \quad (37) $$

Moving on to the mass term, using the first estimate in (31), there exists a constant $c'$ independent of the space discretization such that

$$ \int_{\Omega}^{Q} \|\mathbf{v}_h^*\|_3^2 d\Omega = \int_{\Omega}^{Q} \varphi_I^2 d\Omega \le c' h^2 \eta \sum_{K_{\ell,m} \subset \text{supp}\{\varphi_I\}} \int_{\hat{K}}^{Q} \hat{\varphi}_j d\hat{K}, \quad (38) $$

where $j$ is the local index such that $I = \mathcal{LG}(j, \ell, m)$. Combining the last two inequalities we obtain the lower bound in (33) with

$$ \hat{c} = \frac{c}{c'} \left( \int_{\hat{K}}^{Q} \|\nabla \hat{\varphi}_i\|_3^2 d\hat{K} \right) \left( \sum_{K_{\ell,m} \subset \text{supp}\{\varphi_I\}} \int_{\hat{K}}^{Q} \hat{\varphi}_j d\hat{K} \right)^{-1}. \quad (39) $$

Note that, upon the quasi-uniformity Assumptions 2 and 3, the number of elements forming the support of a basis function $\varphi_I$ is necessarily bounded uniformly w.r.t. the space discretization parameters $h$ and $\eta$. Mention (as an exercise) of this result can be found at the end of Chapter 3 in the reference book [28]. Therefore, we retrieve a constant that depends exclusively on the reference element and the local polynomial space. $\square$
**Lemma 6.** There exists two constants $c_\nu > 0$ and $C_\nu > 0$ independent of $\eta$ such that

$$c_\nu \eta^{-2} \le \sup_{\mathbf{v}_h \in \mathbf{V}_h} \frac{a_{h,\nu\nu}(\mathbf{v}_h, \mathbf{v}_h)}{m_h(\mathbf{v}_h, \mathbf{v}_h)} \le C_\nu \eta^{-2}. \quad (40)$$

Furthermore, there exists two constants $\hat{c} > 0$ and $\hat{C} > 0$ that depends only on the reference element and the local polynomial space $Q_{p_\tau, p_\nu}(\hat{K})$ such that

$$c_\nu = \hat{c} \frac{\alpha}{\varrho_+} \frac{\mu \tau^2}{\sigma^2}, \quad C_\nu = \hat{C} \frac{\beta}{\varrho_-} \frac{\sigma^2}{\mu^3 \tau^2}.$$

*Proof.* The demonstration follows the same principles as the one detailed for Lemma 5.

*Upper bound in (40).* Denoting by $\mathbf{A}_\nu = \mathbf{I} - (\mathbf{e}_1 \otimes \mathbf{e}_1 + \mathbf{e}_2 \otimes \mathbf{e}_2)$ and following the same arguments as for the proof of the upper bound in Lemma 5, we have

$$\forall \mathbf{v}_h = (v_{h,1}, v_{h,2}, v_{h,3})^\top \in \mathbf{V}_h^\flat, \quad \forall K \in \mathcal{T}_h(\overline{\Omega}), \quad a_{h,\nu\nu}^K(\mathbf{v}_h, \mathbf{v}_h) \le \beta \sum_{i=1}^3 \int_{\hat{K}}^Q \|\mathbf{A}_\nu \nabla F_{\ell,m}^{-\top} \nabla \hat{v}_i\|_3^2 |J_{\ell,m}| \mathrm{d}\hat{K},$$

where $K$ has indexes $\ell \in [1, N_h]$ and $m \in [1, N_\eta]$, and $\hat{v}_i = (v_{h,i} \circ \mathbf{F}_{\ell,m}) \in Q_{p_\tau, p_\nu}(\hat{K})$. As previously we remark that the matrix product is

$$\mathbf{A}_\nu \nabla \mathbf{F}_{\ell,m}^{-\top} = \begin{pmatrix} 0 & 0 \\ 0 & (z_{m+1} - z_m)^{-1} \end{pmatrix}.$$

Therefore using the first estimate in (31) and the Assumption 3 of quasi-uniformity of the discretization in the thickness of the plate, there exists a constant $C > 0$ independent of the element such that

$$a_{h,\nu\nu}^K(\mathbf{v}_h, \mathbf{v}_h) \le C \beta (h_\ell^2 \eta) \left(\frac{1}{\mu\eta}\right)^2 \int_{\hat{K}}^Q |\nabla \hat{v}|^2 \mathrm{d}\hat{K} \le C\beta \left(\frac{h^2}{\mu^2\eta}\right) \int_{\hat{K}}^Q |\nabla \hat{v}|^2 \mathrm{d}\hat{K},$$

where $\hat{v} = (\hat{v}_1, \hat{v}_2, \hat{v}_3)^\top$. Combining with the lower bound (36) on the local bilinear form $m_h^K(\cdot, \cdot)$ we obtain the desired result. Note that the constant appearing in the expression of $C_\nu$ is identical to the one obtained in Lemma 5, and is given in (22).

*Lower bound in (40).* Following the same line of arguments as for the previous Lemma, we derive a lower bound by choosing a specific test function defined as $\mathbf{v}_h^* = (0, 0, \varphi_I)^\top$ for any global index $I \in [1, n_h]$. for any $K_{\ell,m} \in \mathcal{T}_h(\overline{\Omega})$ such that $K_{\ell,m} \subset \mathrm{supp}\{\varphi_I\}$,

$$\int_{\Omega}^{Q} |\varepsilon_{\nu}(\mathbf{v}_{h}^{*})|^2 d\Omega \geq \int_{K_{\ell,m}}^{Q} \|\mathbf{A}_{\nu} \nabla \varphi_{I}\|_{3}^{2} dK \geq \inf_{\hat{\mathbf{x}} \in \hat{K}} |J_{\ell,m}(\hat{\mathbf{x}})| \eta^{-2} \int_{\hat{K}}^{Q} \|\nabla \hat{\varphi}_{i}\|_{3}^{2} d\hat{K},$$

where $i \in \hat{\mathbf{n}}$ is a local index such that $I = \mathcal{L}\mathcal{G}(i, \ell, m)$.

$$\int_{\Omega}^{Q} |\varepsilon_{\nu}(\mathbf{v}_{h}^{*})|^2 d\Omega \geq c (\rho_{\ell}^2 \mu \eta^{-1}) \int_{\hat{K}}^{Q} \|\nabla \hat{\varphi}_{i}\|_{3}^{2} d\hat{K} \geq c \left(\frac{h^2 \tau^2 \mu}{\sigma^2}\right) \eta^{-1} \int_{\hat{K}}^{Q} \|\nabla \hat{\varphi}_{i}\|_{3}^{2} d\hat{K}. \quad (41)$$

Combining this inequality with (38) we obtain the lower bound in (40), where the constant $\hat{c}$ is identical to (39) up to a factor two. $\square$

**Lemma 7.** There exists two constants $c > 0$ and $C > 0$ independent of the discretization parameters and the plate thickness such that

$$c(h^{-2} + \delta^{-2}\eta^{-2}) \le \sup_{\mathbf{v}_h \in \mathbf{V}_h} \frac{a_h^\delta(\mathbf{v}_h, \mathbf{v}_h)}{m_h(\mathbf{v}_h, \mathbf{v}_h)} \le C(h^{-2} + \delta^{-2}\eta^{-2}).$$
Proof. From the first inequality of Lemma 3 with $\epsilon = \delta^{-1}$ we obtain, for any $\mathbf{v}_h \in \mathbf{V}_h$

$$a_h^\delta(\mathbf{v}_h, \mathbf{v}_h) = (a_{h,\tau\tau} + 2\delta^{-1}a_{h,\tau\nu} + \delta^{-2}a_{h,\nu\nu})(\mathbf{v}_h, \mathbf{v}_h) \le 2(a_{h,\tau\tau} + \delta^{-2}a_{h,\nu\nu})(\mathbf{v}_h, \mathbf{v}_h).$$

Hence, using Lemma 5 and Lemma 6 yields the sought upper bound with $C = 2\max(C_\tau, C_\nu)$. For the lower bound, we carry out the same procedure as for Lemma 5 or Lemma 6, and we choose a specific test function $\mathbf{v}_h^* = (0, 0, \varphi_I)^T$, for any global index $I \in \llbracket 1, n_h \rrbracket$. Using the uniform coercivity of $\mathbf{C}$, the definition (9) of the Green-Lagrange tensor, and the fact that the product $\varepsilon_\tau(\mathbf{v}_h^*) : \varepsilon_\nu(\mathbf{v}_h^*)$ is zero, we obtain

$$a_h^\delta(\mathbf{v}_h^*, \mathbf{v}_h^*) \ge \alpha \left( \int_{\Omega} |\varepsilon_\tau(\mathbf{v}_h^*)|^2 d\Omega + \delta^{-2} \int_{\Omega} |\varepsilon_\nu(\mathbf{v}_h^*)|^2 d\Omega \right).$$

Using the inequalities (37) and (41), there exists a constant $\hat{c} > 0$ that only depends on the reference elements and the local polynomial space such that

$$a_h^\delta(\mathbf{v}_h^*, \mathbf{v}_h^*) \ge \hat{c} \alpha \left( \frac{\mu\eta}{\sigma^2} + \delta^{-2} \frac{h^2\tau^2\mu}{\sigma^2} \eta^{-1} \right).$$

We conclude using (38).

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

The stability condition of the implicit-explicit scheme. Using the estimates provided in the previous subsection we are now able to analyze in details the stability conditions of the different schemes. To start with, let us consider the explicit scheme. From the expression of its stability condition provided by (18) and Lemma 7, one can verify that there exists two constants $c > 0$ and $C > 0$ independent of the plate thickness and the mesh steps such that

$$c(h^{-2} + \delta^{-2}\eta^{-2})^{-\frac{1}{2}} \le \overline{\Delta t}_{\text{ex}} \le C(h^{-2} + \delta^{-2}\eta^{-2})^{-\frac{1}{2}}.$$

This shows that as the plate thickness decreases the stability condition of the time step becomes more restrictive, i.e. $\Delta t = O(\delta)$. Note that, from this estimation we can infer the natural result that choosing $\eta \sim \delta^{-1}$, i.e. a coarse space discretization in the thickness, would render a reasonable stability condition. However, this is not always feasible for at least two types of configurations: when building reference solution where convergence of the discrete solution is insured, or in the presence of stratified material where the mesh step needs to follow the discontinuities of the material. In the next estimates, we show that the implicit-explicit scheme (19) behaves more adequately when the plate thickness decreases.

**Theorem 3.** There exists two constants $c > 0$ and $C > 0$ independent of the plate thickness and the mesh steps such that for $\theta < \frac{1}{4}$,

$$c \left(1 + \frac{1}{4\theta - 1}\right)^{-\frac{1}{2}} h \le \widetilde{\Delta t}_{\text{imex}} \le C \left(1 + \frac{1}{4\theta - 1}\right)^{-\frac{1}{2}} h,$$

whereas for $\theta = \frac{1}{4}$

$$\overline{\Delta t}_{\text{imex}} \le C(h^{-2} + \delta^{-1}(h^{-2} + \eta^{-2}))^{-\frac{1}{2}}.$$

Proof. For $\theta > \frac{1}{4}$, the time step $\widetilde{\Delta t}_{\text{imex}}$ is given by (22) and we conclude using Lemma 5. For $\theta = \frac{1}{4}$, the maximum time step allowed is given by

$$\overline{\Delta t}_{\text{imex}} = 2 \left( \sup_{\mathbf{v}_h \in \mathbf{V}_h} \frac{(a_{\tau\tau,h} + 2\delta^{-1}a_{\tau\nu,h})(\mathbf{v}_h, \mathbf{v}_h)}{m_h(\mathbf{v}_h, \mathbf{v}_h)} \right)^{-\frac{1}{2}}.$$
Using the first inequality of Lemma 3 with $\epsilon = 1$ yields

$$2a_{h,\tau\nu}(\mathbf{v}_h, \mathbf{v}_h) \leq (a_{h,\tau\tau} + a_{h,\nu\nu})(\mathbf{v}_h, \mathbf{v}_h),$$

from which we can deduce that

$$(a_{h,\tau\tau} + 2\delta^{-1}a_{h,\tau\nu})(\mathbf{v}_h, \mathbf{v}_h) \leq ((1 + \delta^{-1})a_{h,\tau\tau} + \delta^{-1}a_{h,\nu\nu})(\mathbf{v}_h, \mathbf{v}_h).$$

Using the upper bounds of Lemma 5 and Lemma 6 give the result of the theorem. □

**Remark 4.** The estimates of Theorem 3 show that for $\theta = \frac{1}{4}$ a sufficient stability condition is at worst in $\Delta t = O(\sqrt{\delta})$, which is – compared to the explicit case – already an improvement. We have also proven that one can obtain a better behavior of the stability condition w.r.t. the plate thickness if $\theta > \frac{1}{4}$. This result entails in particular that there exists a sufficient stability condition of the implicit-explicit scheme that is independent of the plate thickness – and only depends on the discretization parameters in the tangent plane.

## 5. Numerical results

In this section we provide numerical examples in configurations related to ultrasonic testing based on guided wave propagation in plates. A first goal of these results is to verify the order of convergence of the implicit-explicit scheme, and compare it with the second order space-time convergence of the standard leapfrog scheme. Trying to undercover hints of potential numerical locking, particular interests are given to the behavior of the convergence rate as the thickness of the plate decreases. A second goal of these illustrations is to compare the performances of the implicit-explicit scheme to the explicit scheme, a fair comparison being carried out by insuring equivalent precision of the two schemes for the configuration of interests.

### 5.1. Comparing convergence rates between explicit and implicit-explicit schemes

In the following illustrations we consider a domain $\Omega^\delta = (0, 500 \text{ mm})^2 \times (0, 2\delta)$, with two different half-thicknesses $\delta \in \{0.625, 0.078125\} \text{ mm}$. The maximal time of interest is set to $T = 75 \text{ µs}$ – which provides enough time for the fastest wave packet to rebound on the plate vertical boundaries. The material occupying the domain is an aluminium material represented by an isotropic constitutive law. The associated parameters are $\rho^\delta = 2.71 \text{ g} \cdot \text{cm}^{-3}$, $\lambda^\delta = 62 \text{ GPa}$, and $\mu^\delta = 25 \text{ GPa}$, where $\lambda^\delta$ and $\mu^\delta$ are the Lamé coefficients. We assume that these parameters are constant values over the domain $\Omega^\delta$.

Concerning the source term, we consider two cases: a volume source and a surface source. In the first case, the right-hand side of (7) is modified, and defined as

$$\forall \mathbf{v} \in \mathbf{V}, \quad \ell(t, \mathbf{v}) = \int_{\Omega} \mathbf{f}(t) \cdot \mathbf{v} \, d\Omega.$$

The volume source term $\mathbf{f}(\mathbf{x}, t)$ is simply constructed from the product of a Ricker wavelet times a plane Gaussian space function (constant in the thickness direction)

$$\mathbf{f}(\mathbf{x}, t) = \mathbf{f}(x_1, x_2, t) = h(t) \exp \left( - \frac{|x_1 - x_1^0|^2 + |x_2 - x_2^0|^2}{\sigma^2} \right) \mathbf{p}, \quad \forall (\mathbf{x}, t) \in \Omega \times (0, T), \quad (42)$$

where the polarization is $\mathbf{p} = (1, 1, 0)^T$, the origin is $(x_1^0, x_2^0) = (250 \text{ mm}, 250 \text{ mm})$ the standard deviation is $\sigma = 30 \text{ mm}$, and the Ricker wavelet is defined as

$$h(t) = (2\pi^2 (ft - 1)^2 - 1) \exp(-\pi^2 (ft - 1)^2), \quad \forall t \in (0, T),$$
with a central frequency of $f = 100$ kHz. In the surface source case, the right-hand side of (7) is untouched and the space-time dependency of the surface load reads

$$g_+(x_1, x_2, t) = f(x_1, x_2, t), \quad \forall (x_1, x_2, t) \in \Gamma_+ \times (0, T). \quad (43)$$

Note that the condition (2) is not exactly satisfied, but the initial value of the source term is three orders of magnitude smaller than its peak value – which is assumed to be small enough.

We consider two different types of space-time errors. The first one, denoted by $e^0$, corresponds to a $L^\infty$ norm in time and $L^2$ norm in space, while the second type of errors $e_{ij}^1$ corresponds to a $L^\infty$ in time and semi-$H^1$ norm in space. More precisely, considering $\mathbf{u}^*$ a reference solution for the problem (7), we define

$$e^0 = \frac{\sup_{n \in \llbracket 0; N_T \rrbracket} \|\mathbf{u}^*(t^n) - \mathbf{u}_h^n\|_{L^2(\Omega)^3}}{\sup_{n \in \llbracket 0; N_T \rrbracket} \|\mathbf{u}^*(t^n)\|_{L^2(\Omega)^3}}, \quad e_{ij}^1 = \frac{\sup_{n \in \llbracket 0; N_T \rrbracket} \|\partial_j u_i^*(t^n) - \partial_j u_{i,h}^n\|_{L^2(\Omega)}}{\sup_{n \in \llbracket 0; N_T \rrbracket} \|\partial_j u_i^*(t^n)\|_{L^2(\Omega)}}, \quad \forall i, j = 1, \dots, 3. \quad (44)$$

In the previous expression $\mathbf{u}_h^n$ is either a solution of the explicit time scheme (16) or a solution of the implicit-explicit time scheme (19). In terms of space discretization, we prescribe a fixed discretization in the thickness with one element of order two, namely $p_\nu = 2$ and $N_\eta = 1$. The mesh of the mid-surface is composed of regular square-shaped elements. We vary the mesh step in the following manner $h \in \{50, 25, 12.5, 6.25\}$ mm, we prescribe finite elements of order four i.e. $p_\tau = 4$. For each mesh step $h$ the time steps of each schemes are chosen to satisfy the following requirements:

- They are an integer division of a sampling time step $\Delta t_{\text{sample}}$. This avoids in particular errors coming from a potential time interpolation when computing convergence errors (44).

- They follow the same division as the mesh steps of the mid-surface, i.e. from one configuration to another they are divided by two. This allows to theoretically modify at the same rate the constant in the error coming from time and space discretization.

- They are the maximal time steps satisfying the above requirements along with the stability condition of each time schemes. In particular, the time steps of the implicit-explicit scheme are always greater than the time steps of the explicit scheme since the former is independent of the plate thickness.

The reference solution $\mathbf{u}^*$ appearing in (44) is computed using the explicit scheme (16) and a significantly refined mid-surface mesh $h^* = 1.5625$ mm. We show in Table 1 and Table 2 the time steps for each values of the mesh step $h$, for each schemes, and for the two values of the half-thickness. A first observation here is that the time steps are different (and greater for the implicit-explicit scheme). As expected, the time step for the explicit scheme deteriorates when the thickness of the plate decreases, while the time step of the implicit explicit-scheme remains the same.

<table><caption>Table 1: Values of the time steps used for each configurations in the convergence curve for the case of $\delta = 0.625$ mm. Last column corresponds to the mesh step used for computing the reference solution.</caption><tr><td>$h$ (mm)</td><td>50</td><td>25</td><td>12.5</td><td>6.25</td><td>1.5625</td></tr><tr><td>$\Delta t_{\text{ex}}$ (µs)</td><td>0.067397</td><td>0.033698</td><td>0.016849</td><td>0.008425</td><td>0.002106</td></tr><tr><td>$\Delta t_{\text{imex}}$ (µs)</td><td>0.2021936</td><td>0.101097</td><td>0.050548</td><td>0.025274</td><td>—</td></tr></table>

In Figure 2 we plot the $e^0$ convergence curves. We first remark that these results are similar when changing from one type of source to another (from volume to surface). Secondly, we see that, due to its very low time step, the explicit scheme shows an almost optimal convergence rate driven mainly by the order of the space discretization. For the implicit-explicit scheme, after a pre-asymptotic regime, it recovers the convergence rate of order 2 – which is the minimal convergence rate according to a priori error estimation [12].
<table><caption>Table 2: Values of the time steps used for each configurations in the convergence curve for the case of $\delta = 0.078125$ mm. Last column corresponds to the mesh step used for computing the reference solution.</caption><tr><td>h (mm)</td><td>50</td><td>25</td><td>12.5</td><td>6.25</td><td>1.5625</td></tr><tr><td>Δt<sub>ex</sub> (μs)</td><td>0.008791</td><td>0.004395</td><td>0.002198</td><td>0.001099</td><td>0.000275</td></tr><tr><td>Δt<sub>imex</sub> (μs)</td><td>0.2021936</td><td>0.101097</td><td>0.050548</td><td>0.025274</td><td>—</td></tr></table>

<figure><img src="image_5.png" /><figcaption>Figure 2: $L^\infty-L^2$ convergence curves for volume and surface sources, explicit and implicit-explicit schemes, and both thicknesses of the plate. For each point in the curve, values of the time step are provided in Table 1 and 2.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 3:</figcaption></figure>

In Figure 4 we plot the $e_{i1}^1$ and $e_{i2}^1$ convergence curves with $i \in \{1, 2, 3\}$ for a volume source exclusively – the surface source case being very similar. Both the explicit and the implicit-explicit schemes exhibit a
<figure><img src="image_7.png" /><figcaption>Figure 4: $L^\infty$ - semi-$H^1$ convergence curves for a volume source, explicit and implicit-explicit schemes, and both thicknesses of the plate. The first row corresponds to $e_{i1}^1$ while the second row to $e_{i2}^1$ with $i \in \{1, 2, 3\}$.</figcaption></figure>

convergence rate equal or greater than two, showing no sign of numerical locking. Note the symmetrical behavior of $e_{11}^1 \sim e_{22}^1$ and $e_{21}^1 \sim e_{12}^1$. In Figure 5 we plot the $e_{i3}^1$ convergence curves with $i \in \{1, 2, 3\}$ for both volume (first row) and surface (second row) sources. The former case, with convergence rates equal or greater than two, shows no sign of numerical locking. The latter on the other hand shows a strong deterioration of the convergence rate in $e_{13}^1$ and $e_{23}^1$. This phenomenon seems to worsen as the thickness of the plate. It is observed for both the explicit and the implicit-explicit schemes. Hence we cannot conclude that it is due to the specific nature of the implicit part in (19), but rather from a potential lack of regularity of the exact solution of the problem (7) in the case of surface source terms. Note once again the symmetrical behavior $e_{13}^1 \sim e_{23}^1$ for both volume and surface sources.

## 5.2. Comparing performances between explicit and implicit-explicit scheme

*Illustration on a homogeneous isotropic plate.* As a first illustration of performances we start with an aluminium isotropic plate. The configuration parameters are the same as the one presented in the previous section. The previous convergence study shows less that $10^{-2}$ for the $L^\infty - L^2$ relative error for both the explicit and implicit-explicit scheme when $h = 25$ mm. Hence we choose this configuration with the volume source (42) as a comparison point for the performances of each scheme “at equivalent precision”. To render more eloquent results, we increase the maximal time to $T = 300$ µs.

We show in Figure 6 a series of six snapshots showing the magnitude of the displacement field, obtained with the implicit-explicit scheme for $\delta = 0.078125$ mm. Additionally, we show in Figure 7 extractions of the first component of the displacement at three observation points. These plots show an excellent agreement between the explicit and implicit-explicit schemes even though the time steps are very different,
<figure><img src="image_8.png" /><figcaption>Figure 5: $L^\infty$ - semi-$H^1$ convergence curves $e_{i3}^1$ with $i \in \{1, 2, 3\}$ for the explicit and implicit-explicit schemes, and both thicknesses of the plate. (First row) volume source; (Second row) surface source.</figcaption></figure>

with $\Delta t_{ex} = 0.004395$ µs for the explicit scheme and $\Delta t_{imex} = 0.101097$ µs for the implicit-explicit scheme. In Table 3 we show the performances of the different simulation runs. Overall we see a clear improvement brought by the independence of the time step of the implicit-explicit scheme w.r.t the plate thickness. As expected the speed up grows as the plate of the thickness decreases. It should be noted that this factor is not exactly equivalent to the ratio of the time steps since, at each iteration of the implicit-explicit scheme a set of linear systems are solved. Computing at an initializing step the Cholesky decomposition – since the linear systems are symmetric positive definite – we apply the forward and backward substitutions at each time step. Note that since each linear systems are independent from one degree of freedom in the tangent plane to another, the forward and backward substitutions are performed in parallel. Overall, the gain on the time step is larger than the cost of solving these linear systems and we recover a significant speed up factor.

<table><caption>Table 3: Comparison of performances between explicit and implicit-explicit scheme for the isotropic aluminium plate configuration. Computations performed on a laptop PC, Intel(R) Core(TM) i7-8850H CPU @ 2.60 GHz</caption><thead><tr><th></th><th>Δt (μs)</th><th>CPU* (s)</th><th>Δt (μs)</th><th>CPU* (s)</th></tr></thead><tbody><tr><td>EX.</td><td>0.033698</td><td>7</td><td>0.004395</td><td>50</td></tr><tr><td>IM. - EX.</td><td>0.101097</td><td>3</td><td>0.101097</td><td>3</td></tr><tr><td>ratio</td><td>(3)<sup>-1</sup></td><td>2.33</td><td>(23)<sup>-1</sup></td><td>16.66</td></tr><tr><td></td><td colspan="2">δ = 0.625 mm</td><td colspan="2">δ = 0.078125 mm</td></tr></tbody></table>
<figure><img src="image_9.png" /><figcaption>Figure 6: Snapshots of the displacement magnitude obtained with the implicit-explicit scheme for the aluminium isotropic plate of half-thickness $\delta = 0.078125$ mm.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 7: First component of the displacement field over time at three observation points for the explicit and implicit-explicit schemes for the aluminium isotropic plate of half-thickness $\delta = 0.078125$ mm.</figcaption></figure>

*Illustration on a stratified anisotropic plate.* We finally provide a numerical illustration more related to the scope of application of our work. We consider a plate occupying the domain $\Omega^\delta = (0, 1000 \text{ mm})^2 \times (0, 3 \text{ mm})$ and composed of a stacking or stratification of twenty plies – a configuration commonly encountered with composite materials in aeronautics for instance. Each ply has a thickness of 0.15 mm, and is made of an anisotropic material. The mass density of the plies is set to $\rho^\delta = 1.6 \text{ g} \cdot \text{cm}^{-3}$. For the plies with odd or even
indexes, the elastic modulii following the Voigt's notations are given by

$$ C_{\text{odd}} = \begin{pmatrix} 143.2 & 7.5 & 7.5 \\ 7.5 & 15.8 & 8.2 \\ 7.5 & 8.2 & 15.8 \end{pmatrix}_{3.8 \ 7 \ 7} \quad \text{GPa, } C_{\text{even}} = \begin{pmatrix} 15.8 & 7.5 & 8.2 \\ 7.5 & 143.2 & 7.5 \\ 8.2 & 7.5 & 15.8 \end{pmatrix}_{7 \ 3.8 \ 7} \quad \text{GPa.} $$

These modulii are representing a uni-directional fibered material with a fiber orientation in the x-direction for odd indexes or in the y-direction for even indexes. We choose a surface source term as in (43) with a different polarization of $p = (1,0,0)^{\top}$. As for the previous illustration we set $h = 25$ mm, $p_{\tau} = 4$ and for each ply corresponds one element, i.e. $N_{\eta} = 20$, of order two, $p_{\nu} = 2$. The maximal time is also set to $T = 300$ µs. The time steps of the explicit and implicit-explicit schemes are directly obtained from the CFL condition (18) and (22) respectively.

We present in Figure 8 a series of six snapshots showing the magnitude of the displacement field obtained using the implicit-explicit scheme. We clearly observe on this figure the effect of the anisotropic layers of the material. Additionally, we plot in Figure 9 the first component of the displacement field at three observation points, obtained using the explicit and implicit-explicit schemes. Once again we observe an excellent agreement between the two schemes, even though the time step are very different. In Table 4 we show the ratio between the time step of both schemes and the overall speed up gain obtained through the implicit-explicit scheme. Once again we observe that the gain on the time step greatly overcome the cost overhead of the linear systems to solve at each time step of the implicit-explicit scheme.

<table><caption>Table 4: Comparison of performances between explicit and implicit-explicit scheme for the stratified anisotropic plate configuration. Computations performed on a laptop PC, Intel(R) Core(TM) i7-8850H CPU @ 2.60 GHz</caption><thead><tr><th></th><th>Δt (μs)</th><th>CPU* (s)</th></tr></thead><tbody><tr><td>EX.</td><td>0.0175419</td><td>964</td></tr><tr><td>IM. - EX.</td><td>0.171809</td><td>156</td></tr><tr><td>ratio</td><td>(9.7)<sup>-1</sup></td><td>6.2</td></tr></tbody></table>

## 6. Conclusions and perspectives

In our work, we have addressed the issue of wave propagation modelling in thin plate-like geometries, without any particular assumptions on the material, thus including simple homogeneous isotropic materials or stratified anisotropic materials. The presented approach lies within the context of time-domain spectral finite element methods. By identifying the part of the stiffness operator that is the most penalizing for the time step of the explicit leapfrog scheme, we have been able to derive a new implicit-explicit scheme. After analyzing its related stability through standard energy arguments, we have proven that there exists a sufficient stability condition that is independent of the thickness, its potential sub-division (in the case of stratified materials) and its discretization parameters (number of elements and order of approximation). This implicit-explicit scheme requires however to solve at each time step a linear system. Upon the assumption of an extruded mesh, we have shown that this system can be partially lumped to obtain a block-diagonal matrix. Each block are symmetric positive definite sub-matrices, involving only the degrees of freedom in the thickness and independent from one degree of freedom in the tangent plan to another. Therefore, the cost overhead of solving this system is limited. In practice, we store the Cholesky decomposition in an initializing phase, and perform the forward and backward substitutions within a parallel loop over each block – i.e. over each degrees of freedom in the tangent plane. In terms of numerical illustrations, we first have presented evidence of an adequate convergence behavior in $L^2$ and semi-$H^1$ norms. The convergence of the implicit-explicit scheme appears to be similar to the one observed for the explicit scheme. They both show robustness as the plate thickness decreases, with no sign of numerical locking in the case of a volume
<figure><img src="image_11.png" /><figcaption>Figure 8: Snapshots of the displacement magnitude obtained with the implicit-explicit scheme for the stratified anisotropic plate.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 9: First component of the displacement field over time at three observation points for the explicit and implicit-explicit schemes for the stratified anisotropic plate.</figcaption></figure>

source. However, a specific point needs to be emphasized in the case of a surface load. It appears that both the explicit and implicit-explicit schemes lose optimal convergence in $\partial_3 u_1$ and $\partial_3 u_2$ as the plate thickness decreases. Since both schemes suffer from the same loss, one cannot target the implicit-explicit scheme as the sole culprit. Further investigations should be carried out – typically on the explicit scheme alone – to better understand this behavior. In a second set of numerical illustrations we have compared the performances of the implicit-explicit scheme to the performances of the explicit scheme in a simple homogeneous isotropic case and in a more complex case of a stratified anisotropic material. In both cases we have observed a significant speed up factor, showing that the gain on the time step is far greater than the cost of solving the (small) linear systems at each time step. This speed up factor becomes more important as the plate
thickness decreases, thus proving the interest of our approach for the type of configurations considered in the
scope of this article. By combining robustness w.r.t. the plate thickness, efficiency and consistency with the
three dimensional elastodynamics equation, the application of this new scheme ranges from direct numerical
solving for industrial studies, to generating reference solutions for the validation of other approximated or
asymptotic plate models.

A natural perspective of the presented work is to extend the implicit-explicit scheme to the case of a
curved geometries. Compared to the plane case, in this type of configurations the Jacobian matrix (30) has
non zero extra-diagonal terms. This implies in particular that, after the change of variables to the reference
element, bilinear forms of the type of (32) implicate all three derivatives of the basis functions – and not
only the derivative in the thickness direction as for the plane case. As a consequence, one looses the partial
lumping of the linear system of the implicit-explicit scheme defined as is. A potential workaround is to
isolate the contributions of the tangential derivatives – emanating from the curvature of the plate – and to
relocate them within the explicit part of the scheme. By doing so we would recover a partially lumped linear
system, the price to pay is a stability condition more restrictive than for the plane case – depending on the
amplitude of the local curvature. Another interesting prospect could be to use stabilized leapfrog-Chebyshev
schemes [29, 30, 31] to build a fully explicit scheme with a stability condition that is less sensitive to the plate
thickness. One could expect from this approach to apply multiple times the operator associated to $a_{h,\nu\nu}(\cdot, \cdot)$
within a single time step. Since this operator is cheaper to apply than the complete stiffness operator, one
could obtain some performance benefits depending on the gain on the time step. Compared to the presented
implicit-explicit scheme, this approach would have simpler interactions with standard explicit schemes, for
instance in configurations where the thin structure is embedded within large body of other materials.

## Acknowledgements

Authors would like to thank Prof. Sonia Fliss (POEMS, CNRS, Inria, ENSTA Paris, Institut Polytech-
nique de Paris) for the fruitful discussions during the construction of the presented approach, and for the
comments leading to improvements on the structure of this communication.

## Data availability statement

Data sharing is not applicable to this article as no datasets were generated or analysed during the current
study.

## References

[1] D. Chimenti, Guided waves in plates and their use in materials characterization.

[2] C. Willberg, S. Duczek, J. M. Vivar-Perez, Z. B. Ahmad, Simulation methods for guided wave-based structural health monitoring: a review, Applied Mechanics Reviews 67 (1).

[3] M. Dauge, E. Faou, Z. Yosibash, Plates and Shells: Asymptotic Expansions and Hierarchic Models, American Cancer Society, 2004, Ch. 8.

[4] M. Dauge, I. Djurdjevic, A. Rössle, Full asymptotic expansions for thin elastic free plates, Comptes Rendus de l'Académie des Sciences - Series I - Mathematics 326 (10) (1998) 1243-1248.

[5] M. Dauge, I. Gruais, Asymptotics of arbitrary order for a thin elastic clamped plate, i. optimal error estimates, Asymptotic Analysis 13 (2) (1996) 167-197.

[6] P. G. Ciarlet, Plates and junctions in elastic multi-structures : an asymptotic analysis, Recherches en mathématiques appliquées 14., Masson, Springer-Verlag, 1990.

[7] P. G. Ciarlet, Mathematical Elasticity, Volume 2: Theory of Plates, Studies in Mathematics and its Applications, Elsevier, Academic Press, 1997.

[8] D. Chapelle, B. K.-J., The finite element analysis of shells - fundamentals, 2nd Edition, Computational Fluid and Solid Mechanics, Springer-Verlag Berlin Heidelberg, 2011.

[9] G. Cohen, P. Grob, Mixed higher order spectral finite elements for reissner-mindlin equations, SIAM Journal on Scientific Computing 29 (3) (2007) 986-1005.

[10] D. Komatitsch, J.-P. Vilotte, R. Vai, J. M. Castillo-Covarrubias, F. J. Sanchez-Sesma, The spectral element method for elastic wave equations-application to 2-D and 3-D seismic problems, International Journal for numerical methods in engineering 45 (9) (1999) 1139-1164.
[11] G. Cohen, Higher-Order Numerical Methods for Transient Wave Equations, Scientific computation, Springer, 2002.

[12] 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.

[13] G. Seriani, E. Priolo, Spectral element method for acoustic wave simulation in heterogeneous media, Finite elements in analysis and design 16 (3-4) (1994) 337–348.

[14] J. D. De Basabe, M. K. Sen, Grid dispersion and stability criteria of some common finite-element methods for acoustic and elastic wave equations, Geophysics 72 (6) (2007) T81–T95.

[15] G. Seriani, S. P. Oliveira, Dft modal analysis of spectral element methods for acoustic wave propagation, Journal of Computational Acoustics 16 (04) (2008) 531–561.

[16] W. Ostachowicz, P. Kudela, M. Krawczuk, A. Zak, Guided waves in structures for SHM: the time-domain spectral element method, John Wiley & Sons, 2011.

[17] O. Mesnil, A. Recoquillay, T. Druet, V. Serey, H. T. Hoang, A. Imperiale, E. Demaldent, Experimental validation of transient spectral finite element simulation tools dedicated to guided wave-based structural health monitoring, Journal of Nondestructive Evaluation, Diagnostics and Prognostics of Engineering Systems 4 (4) (2021) 041003.

[18] 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 Eq. 25 (3) (2009) 526–551.

[19] 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.

[20] 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.

[21] T. Rylander, A. Bondeson, Stability of explicit-implicit hybrid time-stepping schemes for maxwell's equations, Journal of Computational Physics 179 (2) (2002) 426–438.

[22] S. Descombes, S. Lanteri, L. Moya, Locally implicit time integration strategies in a discontinuous galerkin method for maxwell's equations, Journal of Scientific Computing 56 (1) (2013) 190–218.

[23] A. B. Hamad, G. Beck, S. Imperiale, P. Joly, An efficient numerical method for time domain electromagnetic wave propagation in co-axial cables, Computational Methods in Applied Mathematics.

[24] J. Chabassier, S. Imperiale, Fourth-order energy-preserving locally implicit time discretization for linear wave equations, International Journal for numerical methods in Engineering 106 (8) (2016) 593–622.

[25] J. Chabassier, S. Imperiale, Space/time convergence analysis of a class of conservative schemes for linear wave equations, Comptes Rendus Mathematique 355 (3) (2017) 282–289.

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

[27] V. Girault, P.-A. Raviart, Finite element methods for Navier-Stokes equations: theory and algorithms, Vol. 5, Springer Science & Business Media, 2012.

[28] P. G. Ciarlet, The finite element method for elliptic problems, SIAM, 2002.

[29] C. Carle, M. Hochbruck, A. Sturm, On leapfrog-chebyshev schemes, SIAM Journal on Numerical Analysis 58 (4) (2020) 2404–2433.

[30] M. J. Grote, S. Michel, S. A. Sauter, Stabilized leapfrog based local time-stepping method for the wave equation, Mathematics of Computation.

[31] J. Chabassier, S. Imperiale, Construction and convergence analysis of conservative second order local time discretisation for wave equations based on domain decomposition, working paper or preprint (Aug. 2019).