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

We propose a new implicit–explicit scheme to address the challenge of modeling 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.

example References 3-8.However, fewer studies concern the use of asymptotic methods for wave propagation in these media-see for instance Reference 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-for example, the so-called semi-analytical finite element methods or the finite difference (FD) schemes and their extensions, see Reference 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-12This approach has proven to be particularly efficient for long distance propagation and accurate rendering of wave packets dispersion, 13-15 and it has been successfully applied to guided wave propagation modeling. 16,17In 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-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 propagation 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  in the following.This leads to an equivalent re-scaled formulation of the problem where the dependency of the differential operators w.r.t  is clearly visible.Applying the standard explicit scheme on this re-scaled formulation shows that its related stability condition is proportionate to .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 for example, the related prior works of References 21-23.For this implicit-explicit scheme we are able to prove that there exists a sufficient stability condition that is independent of .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, that is, 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 article 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 .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 .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 .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.

SEMI-DISCRETE ELASTIC WAVE PROPAGATION PROBLEM
As shown in Figure 1, we consider a plate occupying the open and bounded volume Ω  = Σ × (−, ) ⊂ R 3 , where Σ ⊂ 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 where u  is the time and space dependent displacement associated to the mechanical waves, and (⋅) is the linearized Green-Lagrange tensor.The solid is characterized by its mass density   , and the fourth-order tensor C  representing the linear stress-strain relation, that is, Hooke's law.Denoting by Γ  + the upper boundary of the domain Ω  , relation ( 1) is completed with the following boundary conditions where n is the normal vector field of the boundary.For any point (x 1 , x 2 ) ⊺ ∈ Γ  + , the time-dependent source term g + (x 1 , x 2 , t) typically represents the effect of an ultrasonic actuator.For simplicity we assume that and we also provide null initial conditions in displacement field and velocity, namely As described in Figure 1, we transform each point . Note that Ω is independent of the domain half-thickness, and for this "fixed" geometry we denote by Γ + its upper boundary.Concerning the physical parameters, we assume that there exists a scalar field  and a fourth-order tensor field C defined in Ω such that The scalar field  is a strictly positive bounded function independent of  and we assume that there exists two constants  − and  + such that Concerning the fourth-order tensor field C we assume that it is uniformly coercive and bounded.More precisely, denoting by (R 3 ) the space of tensors-assimilated to linear mappings from R 3 into itself-equipped with its natural inner product and norm and by  s (R 3 ) the subspace of (R 3 ) of symmetric tensors, the mentioned uniform properties read: there exists two constants  > 0 and  > 0 such that Additionally, in what follows we assume that the scalar field  and the tensor field 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. 12

2.1
Re-scaled continuous model problem Lemma 1. Defining V = H 1 (Ω) 3 , the re-scaled weak formulation of problem (1) reads: for any time t ∈ (0, T), find u(t) ∈ V such that for any v ∈ V and Proof.Let us denote by u the re-scaled displacement defined on the fixed geometry.It still depends on  even though we do not make explicit this dependency in order to lighten the presentation.In particular, we have This change of variable transfer the dependency w.r.t  of the geometry onto the differential operators.Namely, the partial derivatives in space of the re-scaled displacement read Let u  be the solution of (1) such that u  (t) ∈ V  = H 1 (Ω  ) 3 , it satisfies the following weak formulation After operating the change of variable from Ω  to Ω, and upon the assumption (4) on the mass density and the stress-strain, we obtain the following equivalent weak formulation which entails the desired result after simplification by .▪ 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 .This will be the fundamental building block for the developments of dedicated time discretizations.For any v = (v 1 , v 2 , v 3 ) ⊺ ∈ V we define the tangential and normal part of the Green-Lagrange tensor as Denoting by and remarking that the re-scaled Green-Lagrange tensor can be expressed as the bilinear form a  (⋅, ⋅) can be rewritten as

Galerkin approximation and semi-discrete weak formulation
Let us denote by V  ⊂ V a Galerkin 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 ∈ (0, T), find u  (t) ∈ V  such that for any 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  (⋅, ⋅) and a   (⋅, ⋅) the resulting approximated forms, the "effective" semi-discrete formulation is and in the next section we will make the following assumption.
Assumption 1.For any  > 0 and for any v  ∈ V  , we have the following inequality 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.

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 .

Explicit, implicit, and implicit-explicit schemes
Let Δt be the time discretization step.We assume that there exists a maximal time step index N T ∈ N * such that T = N T Δt.
For any integer n ∈ ⟦0; N T ⟧, we denote by t n = nΔt and u n  represents the approximation of u  (t n ) the solution of (10).According to the null initial conditions (3) and (2) we have null initial vectors, namely 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  ∈ R a time discretization parameter and

𝔥
another approximation of u  (t n ).The first time scheme considered in our work is the implicit scheme with null initial vectors (12).This scheme is unconditionally stable for any  ≥ 1 4 -which is in itself an interesting robustness property as  vanishes.It can be shown 24 that consistency error is minimized for  = 1 4 .The scheme then reads It is again standard to show that, when (t n , v  ) = 0 one has an energy conservation property, The discrete energy is clearly a positive quadratic functional for any couple (u n+1  , u n  ).This proves the stability of the scheme 12,25 independently of the space discretization parameters, and of .Although this property is highly desirable it comes with the disadvantage that at each iteration the linear system deriving from the bilinear form 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 When mass lumping is achieved 11,12,20,26 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 .To see that, the scheme is rewritten as a second-order perturbation of (14).Using that we obtain that ( 16) is equivalent to As before, when (t n , v  ) = 0 one has an energy conservation property, is defined in (15).The scheme is stable if the energy functional defined above is positive.This can be ensure choosing Δt sufficiently small so that ) 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 The time step Δt ex depends on the space discretization, and is proportional to the small parameter .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  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 References 21-23 and is referred to as an implicit-explicit strategy.It reads 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  + Δt 2  −2 a , )(⋅, ⋅) 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 .
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.

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 , we obtain that ( 19) is equivalent to m with As previously, when (t n , v  ) = 0 one has an energy conservation property, is defined in (15).The stability of the scheme is therefore obtained by ensuring that m (⋅, ⋅) is a positive bilinear form, which leads to the following CFL condition When choosing  = 1 4 one can see that Δt imex depends linearly on √ .This is already an improvement compared to the explicit scheme.However, by choosing  > 1  4 we show below that a sufficient condition on the time step independent of  can be derived-upon Assumption 1.
Theorem 1. Assuming that  > 1  4 , and that Assumption 1 is satisfied, then  Proof.Our starting point is the CFL condition (21).Using Assumption 1, we have for any  > 0 and v  ∈ V  the following inequality Hence, choosing  =  −1 (4 − 1) in the previous expression-which is strictly positive upon the assumption that  > 1 4 -we have We obtain the desired result using (21).▪

Conditioning analysis of the implicit scheme
At each time step of the implicit scheme, a linear system of the following form needs to be solved where and the right hand side depends on the previous steps as 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  ,im (⋅, ⋅) 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 as preconditioner, we obtain a condition number independent of .Additionally, this preconditioning strategy is computationally efficient since the bilinear form a , (⋅, ⋅) is partially lumped, as detailed later in Lemma 4.
Theorem 2. We denote by  ⋆ the condition number of the linear system (23) preconditioned by the operator a  ,imex (⋅, ⋅), namely If Assumption 1 is satisfied and the time step of the implicit scheme satisfies the following condition then we have Proof.We begin by estimating the second term (the supremum) in (24).We have, using Assumption 1, .
The positivity of the bilinear forms m  (⋅, ⋅) and a , (⋅, ⋅) thus entails 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   (⋅, ⋅), one has the following inequality From this estimate we can deduce that for any  ′ > 0 .
Choosing  ′ < 1, from the positivity of the bilinear form  ′ m  (⋅, ⋅) we obtain inf which finally leads to the following estimate for the infimum inf By assumption  < 1, hence one can choose  = √  in the upper bound of the supremum and  ′ = √  in the lower bound of the infimum to obtain the desired result on the condition number.▪ Remark 2. Choosing  = 1 4 , the condition  < 1 in the previous theorem is equivalent to which is the CFL condition for the fully explicit scheme involving only a , .Hence the obtained estimate is valid if the time step satisfies a CFL-like condition independent of .The preconditioning is therefore efficient.Note that the scheme is stable for any choice of Δt, it is however unclear how the conditioning of the underlying system behave when ( 25) is not satisfied.

APPLICATION TO SPACE DISCRETIZATION WITH HIGH ORDER SPECTRAL FINITE ELEMENTS
In this section, we detail a specific construction of the approximation space V  based upon high order spectral finite elements. 10-12On the one hand, we show that using an extruded mesh of the domain Ω one recovers some sort of lumping-referred to as partial lumping-of the matrix associated to the bilinear form a  ,imex (⋅, ⋅).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 Relying on the nature of the geometry Ω, 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 Σ
Let h > 0 be the (planar) space discretization parameter.We consider { h (Σ)} h>0 a family of meshes of the surface Σ.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,  h (Σ) is a set of N h quadrangles generating a tessellation of Σ in the standard sense, without hanging nodes.
For any  ∈ ⟦1; N h ⟧ we denote by   ∈  h (Σ) a quadrangle of the mesh.Each quadrangle   is uniquely defined by its four points {s  j } 3 j=0 ⊂ R 2 , and we denote by f  the unique transformation of the reference quadrangle κ = (0; 1) 2 such that f  (κ) =   .
Introducing some usual definitions-see for example, in Appendix A of Reference 27-we denote by h  = diam(  ) the maximal Euclidean distance between two points in {s  j } 3 j=0 .The space discretization parameter associated to the mesh is then defined as For any j ∈ Z∕4Z, the triangle built from the points s  j−1 , s  j , and s  j+1 is denoted by T  j , and the diameter of the corresponding inscribed circle by   j .Defining the aspect ratio of a quadrangle element is defined as 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 Σ.
Assumption 2. We assume that there exists two constants  > 0 and  > 0 (that are independent of the mesh step) such that for any mesh in { h (Σ)} h>0 we have Note that the left-hand side of (26) entails in particular that h ≤ h  ≤ h, which imposes a mesh composed of elements with similar diameters.Whereas the right-hand side of (26) leads to a bounded   , 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.
With these notations, we give without demonstration the following lemma-see for example, Reference 27 or 28 for more details.Note that in what follows we use the notation ≲ to represent an inequality which is valid up to a constant independent of any element in the family of meshes.
Lemma 2. Denoting j  = det f  , we have For the inverse geometric transformation, similar results hold, namely 4.1.2Quasi-uniform extruded mesh of Ω Let  > 0 be the space discretization parameter in the orthogonal direction of the surface plane.We consider a subdivision of N  edge elements of the interval (−1, 1) such that As for the space discretization of the surface Σ we make the assumption of quasi-uniformity of the subdivision of the interval.
Assumption 3. We assume that there exists a constant  > 0 (independent of ) such that Denoting by  = (h, ) the couple of space discretization parameters of the surface Σ and the interval (−1, 1), we consider a family of space discretizations of the complete domain {  (Ω)}  which are composed of hexahedra built by extrusion of the quadrangles of the surface mesh.Namely, for two fixed mesh parameters h and , and for any  ∈ ⟦1; N h ⟧ and m ∈ ⟦1; N  ⟧ we denote by K ,m ∈   (Ω) an hexahedron of the mesh such that To each hexahedron we associate a unique transformation F ,m of the reference hexahedron K = (0, 1) 3 such that F ,m ( K) = K ,m and defined as from which we can derive the following expression of the Jacobian matrix and its determinant Note that, using this expression of the determinant and the second estimates in ( 27) and ( 28), we have

High-order spectral finite elements
Let us define p  and p  two natural integers representing the order of approximation in the tangent plane and in the orthogonal direction respectively, and n = p  + 1 and n = p  + 1.We define on the reference hexahedron the following polynomial space The dimension of this local polynomial space is denoted by n = n2  n .From this local space, we can build the global conforming approximation space that supports our semi-discretization procedure, namely, In our work, we rely on so-called spectral finite elements. 10-12It 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 { ξ,i } n i=1 and { ξ,i } n i=1 the n and the n Gauss-Lobatto points on the reference edge ê = (0, 1 (see e.g., Reference 12 for an explicit definition and related examples of these points).From these points, we associate the one-dimensional Lagrange polynomials where  k (ê) is the space of polynomials of order lower or equal to k on the reference edge.As Lagrange polynomials they satisfy, by definition, Let us denote by n = ⟦1; n ⟧ 2 × ⟦1; n ⟧, and by i = (i 1 , i 2 , i 3 ) any multi-index such that i ∈ n.Then, on the reference hexahedron K, we construct the Gauss-Lobatto points { ξi } i∈n by tensorization of the one-dimensional points defined on ê, In the same manner, we define the Lagrange polynomials One can verify that they actually are a generating basis of the approximation space, that is,

Numerical integration
We rely on a quadrature formula to numerically compute the integrals involved in the bilinear forms m(⋅, ⋅) and a(⋅, ⋅).
The formula is defined on the reference hexahedron.We denote its associated weights and points respectively by 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 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-12It 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. 19,20,26In the following m  (⋅, ⋅) and a   (⋅, ⋅) introduced in (10) are the approximated bilinear forms where the integral over elements are evaluated using this quadrature formula, namely 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  > 0 and for any v  ∈ V  , we have Proof.Upon the Equation ( 6) of uniform coercivity and continuity on the fourth order tensor C we can insure that for all x ∈ Ω, C(x) defines an equivalent inner product on  s (R 3 ).Hence, using Cauchy-Schwarz and Young inequalities enables us to verify that for any  > 0 we have Using the fact that all the quadrature weights are strictly positive, we obtain the result of the lemma by numerically integrating-that is, using the quadrature formula-over each elements and summing on the elements the previous expression with  =   (v  ) and  =   (v  ).▪

Partial lumping of the discrete stiffness bilinear form a 𝖍,𝝂𝝂
As already mentioned the matrix corresponding to the bilinear form m  (⋅, ⋅) is diagonal.In the remaining part of this section we show that the matrix corresponding to the bilinear form a , (⋅, ⋅) 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  ,imex (⋅, ⋅) 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 ,m ∈   (Ω) and  ∈  0 (K ,m ) be any continuous scalar-valued function, for all ( I ,  J ) ∈ where I = (i, , m), J = (j, , 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 From the tensorized construction of the local Lagrange polynomial, their derivative in the third direction evaluated at the quadrature points reads where i = (i 1 , i 2 , i 3 ) and k = (k 1 , k 2 , k 3 ).Proceeding in the same fashion for  3 φj ( ξk ) with j = (j 1 , j 2 , j 3 ) leads to the desired result.▪ This lemma shows that-locally to each element K ,m -the only test functions for which the integral 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 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  + Δt 2  −2 a , )(⋅, ⋅) 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 presented 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  > 0 and C  > 0 independent of h such that Furthermore, there exists two constants ĉ > 0 and Ĉ > 0 that depends only on the reference element and the local polynomial space Q p  ,p  ( K) such that Proof.In the following, we denote by K any element in   (Ω) of indexes  ∈ ⟦1, N h ⟧ and m ∈ ⟦1, N  ⟧.We decompose the proof by starting with the upper bound and finishing with the lower bound in (33).▪

4.5.2
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 Also, let be the local expression of the mass and stiffness bilinear forms m  (⋅, ⋅) and a , (⋅, ⋅).By construction we have 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  s (R 3 ) and ( 6), we obtain From the definition (8) of the tangential Green-Lagrange tensor, one can verify that , with A  = I − e 3 ⊗ e 3 .Applying the change of variable to the reference element yields ) , and using estimates ( 28) and ( 31), we can insure that there exists a constant C > 0 independent of the element such that where v = (v 1 , v2 , v3 ) ⊺ .Upon the Assumption 2 of quasi-uniformity of the mesh of the surface Σ, this boils down to 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 Using once again the quasi-uniformity Assumption 2, we obtain Combining ( 35) and ( 36) we obtain the upper bound, for any We obtain an upper bound independent of K and v  ∈ V ♭  by taking the supremum over all function v ∈ [Q p  ,p  ( K)] 3 , hence we obtain (33) with which clearly depends exclusively on the reference element and the local polynomial space.

4.5.3
Lower bound in (33) Once again, using the bounds (5) on the mass density and on the coercivity of C in (6), we have In order to provide a lower bound for this supremum we choose a relevant test function.Namely, for any global index I ∈ ⟦1, n  ⟧, we define v ⋆  = (0, 0,  I ) ⊺ .Starting with the stiffness term, for any Performing the change of variable to the reference element leads to where i ∈ n is a local index such that I = (i, , 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 Moving on to the mass term, using the first estimate in (31), there exists a constant c ′ independent of the space discretization such that where j is the local index such that I = (j, , m).Combining the last two inequalities we obtain the lower bound in (33) with Note that, upon the quasi-uniformity Assumptions 2 and 3, the number of elements forming the support of a basis function  I is necessarily bounded uniformly w.r.t the space discretization parameters h and .Mention (as an exercise) of this result can be found at the end of Chapter 3 in the reference book. 28Therefore, we retrieve a constant that depends exclusively on the reference element and the local polynomial space.Lemma 6.There exists two constants c  > 0 and C  > 0 independent of  such that Furthermore, there exists two constants ĉ > 0 and Ĉ > 0 that depends only on the reference element and the local polynomial space Q p  ,p  ( K) such that Proof.The demonstration follows the same principles as the one detailed for Lemma 5. ▪

Upper bound in (41)
Denoting by A  = I − (e 1 ⊗ e 1 + e 2 ⊗ e 2 ) and following the same arguments as for the proof of the upper bound in Lemma 5, we have where K has indexes  ∈ ⟦1, N h ⟧ and m ∈ ⟦1, N  ⟧, and vi = As previously we remark that the matrix product is ) .
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 Combining with the lower bound (36) on the local bilinear form m K  (⋅, ⋅) we obtain the desired result.Note that the constant appearing in the expression of C  is identical to the one obtained in Lemma 5, and is given in (37).

4.5.5
Lower bound in (41) Following the same line of arguments as for the previous lemma, we derive a lower bound by choosing a specific test function defined as v ⋆  = (0, 0,  I ) ⊺ for any global index I ∈ ⟦1, n  ⟧. for any K ,m ∈  h (Ω) such that K ,m ⊂ supp { I }, one can show, in the same way we obtained (38), that where i ∈ n is a local index such that I = (i, , m).Using the estimate just above, estimate (31) for a lower bound of J ,m , and the quasi-uniformity Assumption 2, one can show that there exists a constant c > 0 independent of the space discretization such that Combining this inequality with (39) we obtain the lower bound in (41), where the constant ĉ is identical to (40) up to a factor two.
Lemma 7.There exists two constants c > 0 and C > 0 independent of the discretization parameters and the plate thickness such that Proof.From the first inequality of Lemma 3 with  =  −1 we obtain, for any Hence, using Lemmas 5 and 6 yields the sought upper bound with C = 2 max(C  , C  ).For the lower bound, we carry out the same procedure as for Lemma 5 or 6, and we choose a specific test function v ⋆  = (0, 0,  I ) ⊺ , for any global index I ∈ ⟦1, n  ⟧.Using the uniform coercivity of C, the definition (9) of the Green-Lagrange tensor, and the fact that the product ) .
Using the inequalities (38) and ( 42), there exists a constant ĉ > 0 that only depends on the reference elements and the local polynomial space such that ) .

4.5.6
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 This shows that as the plate thickness decreases the stability condition of the time step becomes more restrictive, that is, Δt = O().Note that, from this estimation we can infer the natural result that choosing  ∼  −1 , that is, 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  < 1  4 , Proof.For  > 1  4 , the time step Δt imex is given by ( 22) and we conclude using Lemma 5.For  = 1 4 , the maximum time step allowed is given by Using the first inequality of Lemma 3 with  = 1 yields from which we can deduce that Using the upper bounds of Lemmas 5 and 6 give the result of the theorem.▪ Remark 4. The estimates of Theorem 3 show that for  = 1 4 a sufficient stability condition is at worst in Δt = O( √ ), 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  > 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.

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.

Comparing convergence rates between explicit and implicit-explicit schemes
In the following illustrations, we consider a domain Ω  = (0,500 mm) 2 × (0, 2), with two different half-thicknesses  ∈ {0.625, 0.078125} mm.The maximal time of interest is set to T = 75 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   = 2.71 g cm −3 ,   = 62 GPa, and   = 25 GPa, where   and   are the Lamé coefficients.We assume that these parameters are constant values over the domain Ω  .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 The volume source term f(x, t) is simply constructed from the product of a Ricker wavelet times a plane Gaussian space function (constant in the thickness direction) where the polarization is p = (1, 1, 0) ⊺ , the origin is (x 0 1 , x 0 2 ) = (250 mm, 250 mm) the standard deviation is  = 30 mm, and the Ricker wavelet is defined as 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 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 ∞ norm in time and L 2 norm in space, while the second type of errors e 1 ij corresponds to a L ∞ in time and semi-H 1 norm in space.More precisely, considering u ⋆ a reference solution for the problem (7), we define In the previous expression u 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  = 2 and N  = 1.The mesh of the mid-surface is composed of regular square-shaped elements.We vary the mesh step in the following manner h ∈ {50, 25, 12.5, 6.25} mm, and, unless specify, we prescribe finite elements of order four that is, p  = 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 Δt sample .This avoids in particular errors coming from a potential time interpolation when computing convergence errors (45).-They follow the same division as the mesh steps of the mid-surface, that is, 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.
Note that the stability condition of each time schemes are computed by estimating the spectral radius (as in Reference 29) of the corresponding global matrices coming from the discretization in space (see also Reference 30 for a local approach at the element level).This is achieved using an iterated power method.
The reference solution u ⋆ appearing in (45) is computed using the explicit scheme ( 16) and a significantly refined mid-surface mesh h ⋆ = 1.5625 mm.We show in Tables 1 and 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.
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).Second, 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. 12In Figure 3, we plot the e 0 and e 1  11 convergence curves for different values of p  .We retrieve the theoretical convergence rates of order p  + 1 for the L ∞ − L 2 norm and p  for the L ∞ − semi H 1 norm.Note for the particular case of order 4, the L ∞ − L 2 convergence is saturated for the lowest mesh step.This is most likely the threshold where the time discretization error dominates.

F I G U R E 3
Convergence curves for the source, explicit scheme, with  = 0.078125 mm, and a variation on the order of the space approximation p  .
In Figure 4, we plot the e 1 i1 and e 1 i2 convergence curves with i ∈ {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 convergence rate equal or greater than two, showing no sign of numerical locking.Note the symmetrical behavior of e 1  11 ∼ e 1 22 and e 1 21 ∼ e 1 12 .In Figure 5, we plot the e 1 i3 convergence curves with i ∈ {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 1  13 and e 1 23 .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 1  13 ∼ e 1 23 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 than 10 −2 for the L ∞ − 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 (43) 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  = 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, with Δt ex = 0.004395 s for the explicit scheme and Δ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.

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 Ω  = (0, 1000 mm) 2 × (0, 3 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   = 1.6 g cm −3 .For the plies with odd or even indexes, the elastic modulii following the Voigt's notations are given by These modulii are representing a unidirectional 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 (44) with a different polarization of p = (1, 0, 0) ⊺ .As for the previous illustration we set h = 25 mm, p  = 4 and for each ply corresponds one element, that is, N  = 20, of order two, p  = 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.

CONCLUSIONS AND PERSPECTIVES
In our work, we have addressed the issue of wave propagation modeling 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 SFEMs.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-that is, 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 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 loose optimal convergence in  3 u 1 and  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 31-33 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 , (⋅, ⋅) 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.

F I G U R E 1
Illustration of the re-scaling of the domain of interest w.r.t the plate half-thickness , and the corresponding mesh elements   ∈  h (Σ) and K ,m ∈   (Ω) of the mid-surface and the re-scaled domain respectively.
To do so, let us introduce the standard Euclidean norm || ⋅ || d in R d , with d = 2 or d = 3, and for any domain D ∈ R d the operator norms

From
this set of points, one can construct the so-called "local-to-global" mapping (i, , m), which to a local point of index i ∈ n within an element of indexes  ∈ ⟦1; N h ⟧ (numbering in the tangent plane) and m ∈ ⟦1; N  ⟧ (numbering in the orthogonal direction) renders the associated global point of index I ∈ ⟦1; n  ⟧.The definition of the global Lagrange basis functions { I } n  I=1 follows U R E 2 L ∞ − 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 Tables1 and 2. (A) Volume source.(B) Surface source.

F
I G U R E 4 L ∞ − 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 1 i1 while the second row to e 1 i2 with i ∈ {1, 2, 3}.

F
I G U R E 5 L ∞ − semi-H 1 convergence curves e 1 i3 with i ∈ {1, 2, 3} for the explicit and implicit-explicit schemes, and both thicknesses of the plate.(First row) Volume source; (second row) surface source.F I G U R E 6Snapshots of the displacement magnitude obtained with the implicit-explicit scheme for the aluminium isotropic plate of half-thickness  = 0.078125 mm.

F I G U R E 8
Snapshots of the displacement magnitude obtained with the implicit-explicit scheme for the stratified anisotropic plat.
Values of the time steps used for each configurations in the convergence curve for the case of  = 0.625 mm.Last column corresponds to the mesh step used for computing the reference solution.Values of the time steps used for each configurations in the convergence curve for the case of  = 0.078125 mm.
Note:TA B L E 2Note: Last column corresponds to the mesh step used for computing the reference solution.
Comparison of performances between explicit and implicit-explicit scheme for the isotropic aluminium plate configuration.