Some comparisons and analyses of time or space discontinuous Galerkin methods applied to elastic wave propagation in anisotropic and heterogeneous media

This research work presents some comparisons and analyses of the time discontinuous space–time Galerkin method and the space discontinuous Galerkin method applied to elastic wave propagation in anisotropic and heterogeneous media. Mechanism of both methods to ensure their stability using time or space discontinuities of unknown fields is analyzed and compared. The most general case of anisotropic and heterogeneous media with physical interfaces of discontinuous material properties is considered, especially for the space discontinuous Galerkin method. A new stability result is proved. Numerical applications to different elastic media, more particularly polycrystalline materials containing a large number of physical interfaces, are also presented to confirm theoretical analyses.

Nowadays, as a classical problem, the numerical modeling of elastic wave propagation can be done reliably and efficiently in a large number of cases. However, for heterogeneous media and when the domain of propagation is large compared to the involved wavelengths and to the characteristic length of heterogeneities, accurate and efficient simulations of wave propagation phenomena still remain a challenging task. The case where the geometries of interior physical interfaces of heterogeneities are complicated, for instance the polycrystalline materials, and the 3D case are even more problematic. Therefore, the research field for the development of effective solvers is still active and it is also important to compare and assess the existing methods that have already proved their effectiveness in simpler cases.

The purpose of the paper is to present some comparisons and analyses between a space discontinuous Galerkin (dG) method and a time discontinuous space–time Galerkin method, hereafter named the time dG method. Both methods use discontinuities of unknown fields to ensure the stability but in different ways leading to different consequences.

The time dG method considered herein is a displacement and velocity two-fields method, which is a specific case of a general time discontinuous space–time finite element method for the second order hyperbolic elastodynamic equations developed by Hughes and Hulbert [1,2,3]. In the weak formulation of the general method, additional stabilizing terms of least-squares type are introduced so it can be viewed as a Petrov–Galerkin method. However, these least-squares terms are not necessary for the stability of the method and can be omitted. Without these terms, we get a time discontinuous Galerkin method. In this case, the only appropriate definition of upwind fluxes (due to the causality), together with the use of the elastic energy inner product to enforce the displacement–velocity compatibility, is sufficient to make the time dG method stable (see “Displacement–velocity two fields time dG method” section). By choosing specific space–time elements, the time dG method can be finally written as a time-stepping scheme and results in an implicit solver. As one important advantage, the method provides an appropriate framework to develop adaptive computing as it remains unconditionally stable even if the space discretization changes over time [4,5,6,7,8]. From one time step to another, the energies that are dissipated in the jumps in time of both displacement and velocity fields guarantee the unconditional stability.

The space dG method^{Footnote 1} is based on the use of spatially element-wise discontinuous finite element basis functions. It has been widely used for computational fluid dynamics, and its application to solid elastodynamics is more recent. We refer to [9,10,11,12,13,14,15,16,17,18] for a non-exhaustive review of related research works. To apply the space dG method, the second-order elastic wave equations are transformed to a first-order hyperbolic system, both velocity–stress [10, 11] and velocity–strain [15] formulations can be used. One important advantage of the space dG method is that the use of discontinuous finite element basis functions allows element-wise solving and makes straightforward the parallel data structures and computing.

For the space dG method, developing appropriate numerical fluxes on element interfaces is the key point for its success. For this purpose, exact solving of the Riemann problem defined on element interfaces is usually recommended. Herein, the upwind numerical fluxes that are exact solutions of the Riemann problem are called exact upwind numerical fluxes and those are only approximate solutions are called approximate upwind numerical fluxes. When the space dG method is applied to elastic media, due to the involvement of the fourth order elastic tensor, the exact upwind numerical fluxes can be easily defined only in the case of continuous material properties [11, 15]. In [15], the exact upwind numerical fluxes were also developed for 2D isotropic elastic media with discontinuous material properties for the velocity–strain formulation.

Recently, in the most general case of multidimensional anisotropic elastic media with discontinuous material properties, a unified and wave oriented variational framework has been proposed by the author [17], which allows a systematic development of the exact upwind numerical fluxes for both velocity-stress and velocity–strain formulations [18]. Owing to the compact and intrinsic forms expressed in terms of tensors, the derivation of upwind numerical fluxes can be done in a well structured way, which allows better understanding of the physical meaning of the developed upwind numerical fluxes. The consistency of the exact numerical fluxes can be easily proved [15, 18]. However, the stability analysis of the space dG method is much more difficult. Wilcox et al. has proved the stability of the velocity–strain dG method in 2D isotropic elastic media including physical interfaces. Their result shows that the energies dissipated in the jumps in space of the velocity and strains fields guarantee the stability for the velocity–strain dG formulation that is not yet discretized in time [15].

In the present paper, the velocity–stress space dG formulation is considered in the general case of anisotropic elastic media with physical interfaces and its stability is analyzed. The result allows explaining the following phenomena previously observed in [17]: when the degree of discontinuity at physical interfaces is high, the use of approximate numerical fluxes leads to instability problems, which cannot be resolved by simply reducing the time step used by the time discretization scheme.

Several numerical examples are presented to illustrate the analytical results. More particularly, polycrystalline materials that present an interesting case of heterogeneous media with a large number of physical interfaces are considered. For this, ultrasonic wave propagation in a single-phase and untextured polycristal composed of elliptic grains is simulated.

We consider the wave propagation in an elastic medium \(\Omega \subset \mathbb {R}^d \) of space dimension d (\(d=1,2,3\)) and in a time interval [0, T]. \(\Omega \) is submitted to dynamic body and surface forces \(\varvec{f}\) and \(\varvec{g}\). \(\varvec{g}\) is applied on a part of \(\Omega \)’s boundary \(\partial \Omega _N \subset \partial \Omega \).

Within the framework of the displacement–velocity two fields time dG method, the elastic wave propagation problem is governed by the following first-order in time equations with as primary unknowns the displacement field \(\varvec{u}(\varvec{x}, t)\) and the velocity field \(\varvec{v}(\varvec{x},t)\): \(\forall (\varvec{x}, t) \in \Omega \times ]0,T[\)

where \(\rho \) denotes the density and \(\varvec{\sigma }\) the second-order Cauchy stress tensor. \(\varvec{\sigma }\) is linked to the second-order infinitesimal strain tensor \(\varvec{\varepsilon }\) and the displacement \(\varvec{u}\) by the Hooke’s law and the definition of \(\varvec{\varepsilon }\) under the hypothesis of small deformations:

In (2), \(\varvec{C}\) is the fourth-order elasticity Hooke tensor, \((\varvec{\cdot })^T\) denotes the transpose operator and “ : ” the usual double dot product between two tensors defined as \((\varvec{C}:\varvec{\varepsilon })_{ij} = \sum _{k,l} C_{ijkl} \varepsilon _{kl} \equiv C_{ijkl} \varepsilon _{kl}\). Herein, the Einstein summation convention is systematically used and all the vectors and tensors are denoted using bold letters. The partial derivative operator with respect to time is denoted \(\partial _t\).

Otherwise, it is useful to recall the usual space gradient and divergence operators defined using an orthonormal basis \((\varvec{e}_i)_{i=1,\ldots ,d}\):

with “\(\otimes \)” the usual tensor product between two vectors: \((\varvec{a}\otimes \varvec{b})_{ij} = a_i b_j\) and “\(\cdot \)” the usual dot product between a tensor and a vector: \((\varvec{A} \cdot \varvec{a})_i = A_{ij} a_j\). In the following, the symmetrized tensor product “\(\otimes _s\)” will also be used and it is defined as: \((\varvec{a}\otimes _s\varvec{b}) = \frac{1}{2}(\varvec{a}\otimes \varvec{b} + \varvec{b}\otimes \varvec{a})\).

Among the two equations of the governing system (1), the first one (1a) expresses the elastodynamic equilibrium and the second one (1b) the compatibility condition between the displacement and velocity fields, which are treated as independent unknowns. However, (1b) is less classical due to the use of the space derivative operator \({{\,\mathrm{\mathbf {Div}}\,}}_{\varvec{x}} (\varvec{\sigma }(\varvec{\cdot }))\), which allows, as will be indicated later, proving the stability of the time dG method in terms of kinetic and elastic energies. It is worth noticing that, due to the special form of the equation (1b), the extension of the time dG method to nonlinear problems is not straightforward, since the meaning of the term \({{\,\mathrm{\mathbf {Div}}\,}}_{\varvec{x}} (\varvec{\sigma }(\partial _t \varvec{u}- \varvec{v}))\) in (1b) becomes unclear. For example, in the case of an elastoplastic problem, if only the elastic part of the stress operator is taken into account to have (1b) well-defined, then the stability analysis of the obtained solver would be difficult.

For the time-dependent problem (1), the following initial conditions are necessary: \(\forall \varvec{x} \in \Omega \)

with \(\mathbf {u}_0\) and \(\varvec{v}_0\) respectively the initial displacement and velocity fields. Finally, to complete the definition of the elastic wave propagation problem, the following boundary conditions are prescribed:

The main idea of the displacement–velocity two fields time dG method is to establish a space–time variational formulation of the strong formulation (1), which is furthermore already discretized in time.

To do this, the whole space–time domain \(S = \Omega \, \times \, ]0,T[\) is subdivided into N space–time slabs: \(\{S_n = \Omega \times \Delta T_n\}_{n = 1,\cdots ,N}\), with \(\Delta T_n =\, ]t_{n-1},t_n[\). Between two successive space–time slabs, the unknown fields \(\mathbf {u}\) and \(\varvec{v}\) can be discontinuous. Then, by considering the balance of the virtual works integrated in the space-time domain, the variational formulation in each space–time slab \(S_n\) expressing the dynamic equilibrium and the displacement–velocity compatibility reads as: \(\forall (\varvec{w}_u, \varvec{w}_v) \in V(S_n) \times V(S_n)\)

where \(V(S_n)\) is the space of kinematically admissible fields \((\varvec{w}_u, \varvec{w}_v)\), \([\varvec{\cdot }(t_{n-1})] = (\varvec{\cdot })(t_{n-1}^+) - (\varvec{\cdot })(t_{n-1}^-)\) denotes the jump quantity in time at \(t_{n-1}\) and \((\varvec{\cdot },\varvec{\cdot })_D\) denotes the inner product between two vector fields or two tensor fields integrated over a domain D, e.g.

It is worth recalling that to obtain (6), the integration by parts is also made in time, then a numerical flux between two successive space–time slabs should be chosen, as the fluxes are not continuous, e.g. at \(t_n\), we have \((\varvec{v}(t_n^-), \varvec{\sigma }(\mathbf {u}(t_n^-))) \ne (\varvec{v}(t_n^+), \varvec{\sigma }(\mathbf {u}(t_n^+)))\). Due to causality, a natural choice is to take the numerical fluxes equal to the upwind fluxes \((\varvec{v}(t_n^-), \varvec{\sigma }(\mathbf {u}(t_n^-)))\).

Now, let us consider the total energy functional \(\mathcal {E}_\Omega (\mathbf {u}(t),\varvec{v}(t))\) over the whole studied domain \(\Omega \) that is defined as follows:

The following result concerning the stability of the time dG method (6) can be proved.

Theorem 1

By assuming that there is no source term inside \(\Omega \) (i.e. \(\varvec{f} = \mathbf {0}\)) and the homogeneous Neumann and Dirichlet boundary conditions (i.e. \(\varvec{g} = \mathbf {0}\) and \(\varvec{u}_D = \mathbf {0}\)), the time dG method (6) is stable in the sense that the following equations hold for each space–time slab \(S_n\) and for the whole space–time domain S:

In the sense that (9) holds also for finite element solutions whatever the space–time mesh used to discretize each space–time slab \(S_n\), the time dG method (6) is unconditionally stable.

Proof

Taking \((\varvec{w}_v,\varvec{w}_u) = (\varvec{v}, \mathbf {u})\) in (6), (9a) is straightforward for each space–time slab \(S_n\). Then summing up all the space–time slabs and noticing that \((\mathbf {u}(\varvec{x}, 0^-),\varvec{v}(\varvec{x}, 0^-)) = (\mathbf {u}_0(\varvec{x}),\varvec{v}_0(\varvec{x}))\), (9b) is immediate. \(\square \)

Hence, the kinetic and elastic energies dissipated in the displacement and velocity jumps in time between two successive space–time slabs guarantee the unconditional stability of the time dG method.

Otherwise, we recall that space–time meshes generally used for the time dG method are composed of a classical finite element mesh in space and one linear element in time in each space–time slab and so an implicit solver is actually obtained [4,5,6,7,8, 17].

The governing equations of the same model problem of elastic wave propagation defined in the preceding are firstly given within the framework of the velocity–stress space dG method. The unified and wave oriented variational framework proposed in [17] is recalled. Then the exact upwind numerical fluxes developed in [17, 18] are given before the presentation of the main result of the present work: the stability of the space dG method in the general case of anisotropic heterogeneous media with physical interfaces of discontinuous material properties.

Strong and variational formulations

We consider the same model problem of elastic wave propagation defined in the preceding section. However, within the framework of the velocity–stress space dG method proposed in [10, 11, 19] (see Eqs. (22.14–22.15) in [19]), its governing equations are put into the following strong first-order hyperbolic system with as primary unknowns the velocity and stress fields \((\varvec{v},\varvec{\sigma })\): \(\forall (\varvec{x}, t) \in \Omega \times ]0,T[\)

where the space derivation operator \(\varvec{A}^{\partial _{\varvec{x}}}\) and its transpose (useful hereafter) are defined as follows: \(\forall \varvec{W} = (\varvec{w} \ \varvec{\tau })^T\)

We note that in this section no body force is considered in the equilibrium equation without loss of generality of the purpose of the present work.

The tensorial compact form in (10) was proposed in [17]: the generalized unknown \(\varvec{U}(\varvec{x},t) = (\varvec{v}(\varvec{x},t) \ \varvec{\sigma }(\varvec{x},t))^T\) composed of \(\varvec{v}\) and \(\varvec{\sigma }\) is a field in \(\mathbb {R}^d \times \mathbb {R}^{d\times d_{sym}}\) and defined over the space–time domain \(S = \Omega \times ]0, T[\), with \(\mathbb {R}^{d\times d_{sym}}\) indicating that \(\varvec{\sigma }\) is a \(d \times d\) symmetric second-order tensor. The inner product in the vectorial space \(\mathbb {R}^d \times \mathbb {R}^{d\times d_{sym}}\) is defined as: \(\forall \varvec{W}_i = (\varvec{w}_i \ \varvec{\tau }_i)^T, (i=1,2)\),

According to the definition of the space derivative operators (3), it is easy to show that, on the boundary \(\partial D\) of any subdomain \(D \subseteq \Omega \), the flux operator \(\varvec{F}_{\varvec{n}}\) (with \(\varvec{n} = n_i\varvec{e}_i\) the outward unit normal vector defined on \(\partial D\)) associated to the first-order system (10) is in fact equal to \(\varvec{A}_{\varvec{n}}\) the Jacobian operator in the \(\varvec{n}\) direction, that is: \(\forall \varvec{W} = (\varvec{w} \ \varvec{\tau })^T\),

In (13), the subscript index “\(\varvec{n}\)” indicates the dependency of \(\varvec{F}_{\varvec{n}}\) and \(\varvec{A}_{\varvec{n}}\) on \(\varvec{n}\). In the following, the local orthonormal basis defined on \(\partial D\) will be denoted by \((\varvec{n}, \{\varvec{t}_\alpha \}_{\alpha = 1,\ldots , d-1})\).

Furthermore, it is useful to recall one important result given in [17]: \(\varvec{A}_n\) can be decomposed using its two eigenbases as follows:

where \((\lambda ^-_{\varvec{n},k}, \varvec{R}^-_{\varvec{n},k}, \varvec{L}^-_{\varvec{n},k})_{k = qL, \{qT_\alpha \}_{\alpha = 1, d-1}}\) and \((\lambda ^+_{\varvec{n},k}, \varvec{R}^+_{\varvec{n},k}, \varvec{L}^+_{\varvec{n},k})_{k = qL, \{qT_\alpha \}_{\alpha = 1, d-1}}\) are respectively the strictly negative and positive eigenvalues and the corresponding right and left eigenvectors of \(\varvec{A}_n\). The left eigenvectors of \(\varvec{A}_n\) are the right eigenvectors of \(\varvec{A}_n^T\).

We recall that, among the \(m = d + d(d+1)\slash 2\) eigenvalues of \(\varvec{A}_{\varvec{n}}\), there are d strictly negative eigenvalues \((\lambda ^-_{\varvec{n},k} = -c_{\varvec{n},k})_{k = qL, \{qT_\alpha \}_{\alpha = 1, \ldots , d-1}}\) and d strictly positive eigenvalues \((\lambda ^+_{\varvec{n},k} = c_{\varvec{n},k})_{k = qL, \{qT_\alpha \}_{\alpha = 1, \ldots , d-1}}\), \(c_{\varvec{n},qL}\) and \(c_{\varvec{n},qT_\alpha }\) being respectively the velocity of quasi longitudinal and quasi transverse wave modes propagating in the \(\varvec{n}\) direction. The subscript indices “qL” and “qT” respectively refer to terms “quasi longitudinal” and “quasi transverse”.

According to [17], the right and left eigenmodes corresponding to the nonzero eigenvalues of \(\varvec{A}_{\varvec{n}}\) read as: \(\forall k = qL, \{qT_\alpha \}_{\alpha = 1, \ldots , d-1}\)

with \(z^\pm _{\varvec{n},k} = \rho \lambda ^\pm _{\varvec{n},k}\) the acoustic impedance, \(\varvec{w}_{\varvec{n}, k} = \frac{1}{\sqrt{2}}\varvec{\gamma }_{\varvec{n},k}\) , \((\varvec{\gamma }_{\varvec{n},k})_{k = qL, \{qT_\alpha \}_{\alpha = 1, \ldots , d-1}}\) the unit eigenvectors of the usual eigensystem of the Christoffel tensor \(\varvec{\Gamma }_{\varvec{n}}\):

The definition of the Christoffel tensor \(\varvec{\Gamma }_{\varvec{n}}\) is recalled in the following:

We recall that the word “quasi” means that in the general anisotropic case we have neither pure longitudinal wave mode verifying \(\varvec{\gamma }_{\varvec{n},qL} \parallel \varvec{n}\) nor pure transverse waves modes verifying \(\varvec{\gamma }_{\varvec{n},qT} \perp \varvec{n}\), in contrast to the isotropic case.

To develop the variational formulation for the space dG method applied to the strong formulation (10), the main idea is to seek an approximated solution \(\varvec{U}_h = (\varvec{v}_h \ \varvec{\sigma }_h)^T\) that are discontinuous in space, so we need to use a finite element mesh and the obtained space dG formulation is already discretized in space.

Let \(\mathcal {M}_h = \{\Omega _k\}_k\) denote a finite element mesh of \(\Omega \). In the following, any element \(\Omega _k\) of the mesh \(\mathcal {M}_h\) will be denoted by E and any of its neighboring elements by \(E'\). The discontinuous solutions in E and \(E'\) are respectively denoted by \(\varvec{U}_h\) and \(\varvec{U}_h'\). Then the space dG variational formulation of the elastic wave model problem (10) for any elemnt E can be put into the following two equivalent forms: \(\forall \varvec{W}_h(\varvec{x}) = (\varvec{w}_h(\varvec{x}) \ \varvec{\tau }_h(\varvec{x}))^T\)

In (18), \(\hat{\varvec{F}}_{\varvec{n}}(\varvec{U}_h, \varvec{U}'_h)\) denotes a numerical flux used to replace the discontinuous flux \(\varvec{F}_{\varvec{n}}(\varvec{U}_h)\) on the element boundary \(\partial E\) when the variational formulation (18a) is established. Hence, an appropriate choice of the numerical flux is essential for the success of the space dG method.

In the following, is firstly recalled the definition of numerical fluxes on the interior element boundary \(\partial E_{int} = \partial E \backslash (\partial E \cap \partial \Omega )\), which should also depend on the solution \(\varvec{U}'_h\) in the neighboring elements \(E'\) of E. As for the exterior element boundary \(\partial E_{ext} = \partial E \cap \partial \Omega \), the definition of the numerical flux \(\hat{\varvec{F}}_{\varvec{n}}(\varvec{U}_h)\) should take into account the boundary conditions \(\varvec{U}'_h\). It was given in [17] and will be recalled in Lemma 1.

Exact upwind numerical fluxes

Let us consider the interface of two adjacent elements E and \(E'\) having respectively \((\rho , \varvec{C}, \varvec{U}_h)\) and \(( \rho ', \varvec{C}', \varvec{U}'_h)\) as densities, elastic moduli and initial states (Fig. 1). The two outward unit normal vectors of E and \(E'\) on their interface are respectively \(\varvec{n}\) and \(\varvec{n}'\) and verify \(\varvec{n} + \varvec{n}' = \mathbf {0}\).

Sketch illustrating the Rankine–Hugoniot jump conditions in the Riemann problem in the 3D case

Full size image

To define exact upwind numerical fluxes, we should consider the Riemann problem governing the states that are results of the propagation of the discontinuity \(\varvec{U}_h - \varvec{U'}_h\) [15, 17,18,19,20]. In the following all the equations are written in the 3D case without loss of generality.

As discussed in [17], solving exactly the Riemann problem at element interfaces is not a sufficient condition to get physically sound numerical solutions. Indeed, different equivalent strong forms of the elastic wave problem give rise to different forms of the Riemann problem and consequently to different interface conditions. For example, the Riemann problem directly derived from the strong form (10) leads to the following interface conditions:

Obviously, on a physical interface with \(\rho \ne \rho '\) or/and \(\varvec{\Gamma }_{\varvec{n}} \ne \varvec{\Gamma }'_{\varvec{n}'}\), (19) are generally not coherent with the classical interface conditions, i.e. the following velocity and stress vector continuities:

In the present work, only the exact upwind numerical fluxes developed in [18] that are coherent with (20) are considered.

Recalling the results given in [18], the Rankine–Hugoniot jump conditions leading to the interface conditions (20) read as:

It is worth noticing that the Riemann problem (21) corresponds to another equivalent form of the first-order velocity–stress system (10):

with

Furthermore, by recalling (15) the definition of \((\varvec{R}^-_{\varvec{n},k},\varvec{L}^-_{\varvec{n},k})\) and (14) the decomposition of \(\varvec{A}_{\varvec{n}}\), it can be shown that:

Solving of the Riemann problem (21) leads to the determination of the six unknown states \(\{\varvec{U}^a, \varvec{U}^b, \varvec{U}^c, \varvec{U}^{a'}, \varvec{U}^{b'}, \varvec{U}^{c'}\}\), i.e. the six characteristic coefficients \(\{\tilde{\alpha }_k, \tilde{\alpha }'_{k}\}_{k=qL, qT_1, qT_2}\). Then, the exact numerical fluxes are defined in the following way:

We remark that (25b) gives another flux definition that is equivalent to (25a). It will be used in the proof of Lemma 1.

In [18], it has been proved that \(\{ \tilde{\alpha }_k, \tilde{\alpha }'_k\}_{k=qL, qT_1, qT_2}\) are the solution of the following linear system of equations:

In (26), \([I_d]\) is the \(d \times d\) identity matrix and \([\tilde{B}]\) and \([\tilde{B}']\) are \(d \times d\) matrices with zero diagonal terms and the following extra-diagonal terms:

with \(\delta z^-_{\varvec{n},kl} = z^-_{\varvec{n},k} - z^-_{\varvec{n},l}\), \(\delta z^{-'}_{\varvec{n}',kl} = z^{-'}_{\varvec{n}',k} - z^{-'}_{\varvec{n}',l}\) and,

\(\overline{z^-_{\varvec{n},k}}^{R}\) and \(\overline{z^-_{\varvec{n},k}}^{V}\) respectively denote the harmonic and arithmetic means of the acoustic impedances of the k-th eigenvector. \(\{ \tilde{\,}\varvec{L}^-_{\varvec{n},k}, \tilde{\,}\varvec{L}^{-'}_{\varvec{n}',k} \}\) are the perturbed left eigenmodes of \(\{ \varvec{A}_{\varvec{n}}, \varvec{A}'_{\varvec{n}'} \}\) calculated by using the material properties of the adjacent element \(E'\) in the following way:

Two simpler but interesting cases to consider are:

• Continuous case

  On an element interface with continuous material properties, it is straightforward that:

  $$\begin{aligned} \tilde{\alpha }_k = \varvec{L}^-_{\varvec{n},k} \cdot (\varvec{U}_h - \varvec{U}'_h) \;,\;\tilde{\alpha }'_k = \varvec{L}^{-'}_{\varvec{n}',k} \cdot (\varvec{U}'_h - \varvec{U}_h) \end{aligned}$$

  (30)

• Discontinuous isotropic case

  On an element interface separating two different isotropic materials, by simply recalling that in the isotropic case \(\varvec{\gamma }_{\varvec{n},qL} = \varvec{n}\), \(\varvec{\gamma }_{\varvec{n},qT_1} = \varvec{t}_1\) and \(\varvec{\gamma }_{\varvec{n},qT_2} = \varvec{t}_2\), it is straightforward that:

  $$\begin{aligned} \tilde{\alpha }_k = \tilde{\,}\varvec{L}^-_{\varvec{n},k} \cdot (\varvec{U}_h - \varvec{U}'_h) \;,\;\tilde{\alpha }'_k = \tilde{\,}\varvec{L}^{-'}_{\varvec{n}',k} \cdot (\varvec{U}'_h - \varvec{U}_h) \end{aligned}$$

  (31)

We note that, only in the general case of an element interface separating two different anisotropic materials, we have \([\tilde{B}] \ne [0]\) and \([\tilde{B}'] \ne [0]\). Otherwise, it is interesting to note that the perturbed left eigenmodes \(\{ \tilde{\,}\varvec{L}^-_{\varvec{n},k}, \tilde{\,}\varvec{L}^{-'}_{\varvec{n}',k} \}\) take into account the coupling between the wave modes of the same type, e.g. both qL-modes, of the adjacent elements E and \(E'\), while the matrices [B] and \([B']\) take into account the coupling between two wave modes of different type, e.g. a qL wave mode from E with a qT wave mode from \(E'\).

Stability analysis

Now the main result concerning the stability of the space dG method (18) using the exact upwind numerical fluxes (25)–(26) is given. Firstly, the stability of the space dG method is proved in both continuous case and discontinuous isotropic case whatever the space dimension. Then, in the general case of anisotropic media with physical interfaces separating different materials, the proof, which is much more complicated, is given only for the space dimension equal to 2.

Hereafter, \(\Gamma _{EE'} = \partial E \cap \partial E'\) denotes the interface between any pair of adjacent elements E and \(E'\) and the following notations are used for space jumps across it:

Let us consider the total energy functional \(\mathcal {E}_\Omega (t)\) over the whole studied domain \(\Omega \) that is defined within the framework of the space dG method as follows:

Firstly, the following lemma can be proved.

Lemma 1

By assuming that there is no source term inside \(\Omega \) (i.e. \(\varvec{f} = \mathbf {0}\)) and the homogeneous Neumann and Dirichlet boundary conditions (i.e. \(\varvec{g} = \mathbf {0}\) and \(\varvec{u}_D = \mathbf {0}\)), we have

Proof

Let us consider any element E. When \(\partial E_{ext} \ne \varnothing \), let us assume, as in [17] (“Exact upwind numerical fluxes” section), that the Dirichlet boundary conditions are strongly imposed by eliminating the prescribed degrees of freedom of velocity in the final system to solve. On the other hand, as also discussed in [17], it is natural to choose the following numerical flux on \(\partial \Omega _N\) where the Neumann boundary conditions are prescribed:

Then, when the strategy of weak verification of the Neumann boundary conditions is chosen, the two equivalent space dG variational formulations (18) become:

and, when the strategy of strong verification of the Neumann boundary conditions is chosen, i.e. \(\varvec{\sigma }_h\cdot \varvec{n} = \varvec{g}\), they become:

Now let us sum up the two equivalent space dG variational formulations (36) or (37), divide the obtained equation by 2 and take

then, by remarking that the following equations are true:

we finally get:

Now, by summing up all the elements, remarking that:

and taking into the definition (25b) of the numerical fluxes \(\hat{\tilde{\varvec{F}}}_{\varvec{n}}(\varvec{U}_h, \varvec{U}'_h)\) and (24a), the Eq. (34) can be finally proved. \(\square \)

Now, using Lemma 1, the following results concerning the stability of the space dG method can be proved.

Theorem 2

In the continuous case, the space dG method is stable in the sense that:

with \(|\tilde{\varvec{A}}_{\varvec{n}}| = |z^\pm _{\varvec{n},k}| \varvec{L}^\pm _{\varvec{n},k} \otimes \varvec{L}^\pm _{\varvec{n},k}\).

Proof

In the continuous case, \(\{\tilde{\alpha }_k, \tilde{\alpha }'_{k}\}_{k=qL, qT_1, qT_2}\) are given by (30) and the following equations hold:

Using the decomposition (24b) of \(\tilde{\varvec{A}}_{\varvec{n}}\), the equation in (42) is straightforward. As \(|\tilde{\varvec{A}}_{\varvec{n}}|\) is symmetric and positive-definite, the inequation in (42) is obvious. \(\square \)

Theorem 3

In the discontinuous isotropic case, the space dG method is stable in the sense that:

Proof

In this case, \(\{\tilde{\alpha }_k, \tilde{\alpha }'_{k}\}_{k=qL, qT_1, qT_2}\) are given in (31), then we get:

with \(\varvec{A}^* = |z^-_{\varvec{n},k}|\varvec{L}^-_{\varvec{n},k}\otimes \tilde{\,}\varvec{L}^-_{\varvec{n},k} + |z^{-'}_{\varvec{n}',k}|\varvec{L}^{-'}_{\varvec{n}',k}\otimes \tilde{\,}\varvec{L}^{-'}_{\varvec{n}',k}\).

According to (15), (29) and (28), the definitions of \((\varvec{L}^-_{\varvec{n},k}, \varvec{L}^{-'}_{\varvec{n}',k})\), \((\tilde{\,}\varvec{L}^-_{\varvec{n},k}, \tilde{\,}\varvec{L}^{-'}_{\varvec{n}',k})\) and \((C^-_{z,k}, C^{-'}_{z,k})\), we get:

Herein, the orthonormal basis on \(\Gamma _{EE'}\) chosen for E (resp. \(E'\)) is \((\varvec{n}, \{\varvec{t}_l\}_{l=1,d-1})\) (resp. \((\varvec{n}', \{\varvec{t}'_l\}_{l=1,d-1})\)) with \(\varvec{t}_l + \varvec{t}'_l = \mathbf {0}\). Then by taking into account the fact that in the isotropic case the eigenvector \(\{\varvec{\gamma }_{\varvec{n},k}\}_{k=1,d}\) (resp. \(\{\varvec{\gamma }'_{\varvec{n}',k}\}_{k=1,d}\)) can be taken equal to the elements of the corresponding basis, we have \(\varvec{\gamma }'_{\varvec{n}',k} = - \varvec{\gamma }_{\varvec{n},k}\) (\(\forall k\)). Putting the latter and (46) in (45), it is easy to get (44). \(\square \)

We note that a similar result for the velocity–strain space dG formulation was proved by Wilcox et al. in [15], similar roles played by the harmonic and arithmetic means of the acoustic impedances \(\overline{|z^-_{\varvec{n},k}|}^{R}\) and \(\overline{|z^-_{\varvec{n},k}|}^{V}\) were highlighted.

It can be noticed that the proof of the stability in the discontinuous isotropic case is more tricky than the one in the continuous case. However, thanks to the fact that the relation \(\varvec{\gamma }'_{\varvec{n}',k} = - \varvec{\gamma }_{\varvec{n},k}\) (\(\forall k\)) remains true, the proof remains relatively simple, contrary to the stability analysis in the discontinuous anisotropic case. In the following, only the 2D case is considered and the stability is proved only under a sufficient condition.

In the 2D case, let us assume that four orthonormal bases are used: \((\varvec{n}, \varvec{t})\) and \((\varvec{\gamma }_{\varvec{n},1}, \varvec{\gamma }_{\varvec{n},2})\) for E, and \((\varvec{n}', \varvec{t}')\) and \((\varvec{\gamma }'_{\varvec{n}',1}, \varvec{\gamma }'_{\varvec{n}',2})\) for \(E'\), with \(\varvec{e}_3=\varvec{n}\wedge \varvec{t} = \varvec{n}'\wedge \varvec{t}' = \varvec{\gamma }_{\varvec{n},1}\wedge \varvec{\gamma }_{\varvec{n},2} = \varvec{\gamma }'_{\varvec{n}',1}\wedge \varvec{\gamma }'_{\varvec{n}',2}\). Otherwise, \(\zeta \) defined as follows is also used:

We remark that \(\zeta = 0\) in both continuous and discontinuous isotropic cases.

Theorem 4

In the 2D discontinuous anisotropic case, the space dG method is stable in the sense that:

under the following sufficient condition:

In (48), \(G(C_1, C_2, \Delta ; \varvec{w})\) is defined as follows: \(\forall C_1 \ge 0\), \(C_2 \ge 0\) and \(\varvec{w}\)

with

Proof

See Appendix A. \(\square \)

We remark that in the discontinuous isotropic case, as \(\varvec{\gamma }'_{\varvec{n}',1} = - \varvec{\gamma }_{\varvec{n},1}\) and \(\varvec{\gamma }'_{\varvec{n}',2} = - \varvec{\gamma }_{\varvec{n},2}\), so (48) is identical to (44).

In the discontinuous anisotropic case, Theorem 4 shows that the stability of the space dG method is much more complicated to establish, even in the 2D case. The difficulty is essentially due to the interactions between wave modes of different types, i.e. \(\varvec{\gamma }_{\varvec{n},1}\) with \(\varvec{\gamma }'_{\varvec{n}',2}\) and \(\varvec{\gamma }_{\varvec{n},2}\) with \(\varvec{\gamma }'_{\varvec{n}',1}\), coming from two adjacent elements separated by a physical interface. The kinetic and elastic energies dissipated in the space jumps of velocity \([\varvec{v}_h]\) and of stress vector \([\varvec{\sigma }_h \cdot \varvec{n}]\) can still guarantee the stability, but under a sufficient condition. For the moment, it has not been possible either to remove this condition or prove that it is also a necessary condition. However, we remark that all equations developed in Theorem 4 and in its proof are perfectly symmetric with respect to all quantities coming from E and \(E'\), hence, we believe that it would be difficult to further improve the condition proposed in Theorem 4.

Concerning the sufficient condition (49), it depends obviously on the degree of discontinuity across a physical interface, as \(\Delta ^R_{12}\), \(\Delta ^V_{12}\) and \(\zeta \) are proportional to the latter. Future theoretical and numerical investigations are necessary to better understand this condition. However, a first analysis suggests that it should be satisfied for relatively high degrees of discontinuity.

In conclusion, for both time dG and space dG methods, the mechanism to ensure the stability is the energy dissipation of velocity, displacement or stress vector jumps in time or in space. However, contrary to the stability of the time dG method that is unconditional whatever the space and time discretizations used, the stability of the space dG method is proved for its variational formulation discretized in space but not yet discretized in time. When a time discretization scheme is applied, its own stability condition should also be taken into account. For example, in the present work, the standard four-stage fourth-order Runge–Kutta (RK) method is used for the time discretization of the space dG method, it is an explicit scheme submitted to Courant–Friedrichs–Lewy (CFL) type stability condition [17]. Inversely, even when the stability condition is satisfied by the time discretization method, if the stability is not proved at the continuous level for the space dG variational formulation, typically in the case of the use of approximated numerical fluxes, the stability of discretized solutions cannot be guaranteed.

This section presents several numerical examples of comparison between both time and space dG methods. As previously indicated, concerning the time dG method, a space–time finite element mesh composed of a classical finite element mesh in space and one unique linear element in time is used to discretize each space–time slab. For the comparison, a same space finite element mesh is used by both methods, but the space dG method uses the standard four-stage fourth-order RK algorithm that is explicit and conditionally stable.

The strategy adopted in our works to choose space and time discretization parameters was detailed in [17]. In summary, the element size \(h^E\) is determined by the shortest wavelength of interest, which is defined by the highest frequency of interest. Concerning the time step \(\Delta t\), it is determined by the element size and the following CFL (Courant–Friedrichs–Levy) type conditions: for the space dG method it is a stability condition, while for the time dG method it is a condition to prevent high frequency modes of interest from numerical damping:

In (52), \(N_p\) is the order of FE basis function and \(CFL_{sDG}\) the CFL number in classical sense.

For the numerical examples presented hereafter, only linear (1D) or bilinear (2D, quadrilateral element with 4 nodes) elements are used, and we take \(CFL_{sDG} \approx C_{tDG-damping}\) with \(N_p = 1\). So the time step used for the time dG method is three times larger than the one used for the space dG method.

As far as the external loadings are concerned, their dependency upon the time t is chosen to be a ricker signal, defined as follows:

with \(T_r\) the period of the ricker signal. We recall that the ricker signal (53) provides a perfectly controlled frequency content, its center frequency is \(f_{max} = 2T^{-1}_r\) and its cutoff frequency can be reasonably considered as \(f_c = 3f_{max}\) (see Figure 2 in [17]).

1D case with a physical interface

A 1D elastic rod \(\Omega = \Omega _1 \cup \Omega _2\) composed of two parts of different materials is studied. It was already considered in [17]. The density \(\rho = 2500\; \text {kg}\, \text {m}^{-3}\) is uniform in the entire rod, while the Young’s modulus in the two subdomains is respectively (from left to right) \(E_1 = 22.5\; \text {GPa}\) and \(E_2 = a^2 E_1\). Hence, the phase velocity in the two subdomains is respectively equal to \(c_{L1} = 3000\; \text {m}\,\text {s}^{-1}\) and \(c_{L2} = a\,c_{L1}\). The total length of the rod is \(L = L_1 + L_2 = (1+a)L_1\), with \(L_1 = 900\;\text {m}\). The impedance contrast between the two subdomains is \(\sqrt{E_2\rho _2} \slash \sqrt{E_1\rho _1} = a\). It is obvious that the smaller the value of a, the higher the degree of discontinuity. Three degrees of discontinuity are considered, corresponding to three values of a: 1 (homogeneous case), 0.5 (discontinuous case) and 0.01 (highly discontinuous case).

A pressure loading g(t) is applied on the left end of the rod and the free boundary condition is prescribed on the right end. g(t) is a ricker signal with period \(T_r = 0.04\,\text {s}\) and cutoff frequency \(f_c = 125\,\text {Hz}\). The corresponding shortest wavelength to consider is \(\ell _L(f_c) = 24\,\text {m}\). The time that waves need to propagate from the left end to the physical interface is equal to \(0.2\,\text {s}\), so is the time from the physical interface to the right end. The total time of analysis is taken equal to \(1.6\,\text {s}\), which corresponds to the time required for 2 round trips.

Firstly, the following element sizes \(h^E_1 = 0.6\,\text {m}\) and \(h^E_2 = a\,h^E_1\) are used in the two subdomains, which results in \(N_e = \frac{\lambda _{min}}{h^E}= 40\) elements in the shortest longitudinal wavelength of interest. The time step \(\Delta t = 62.5\,\upmu \text {s}\) is used by the space dG method and \(\Delta t = 20\,\upmu \text {s}\) by the time dG method, leading to \(CFL_{sDG} \approx C_{tDG-damping} \approx 1\).

Decreasing over time of the total energy in the entire rod is studied using two quantities: the first one is the total loss accumulated up to a time step: \(\frac{\mathcal {E}(t)-max(\mathcal {E})}{max(\mathcal {E})}\) (Fig. 2b), and the second one is the incremental loss between two successive time steps: \(\frac{\mathcal {E}(t_n)-\mathcal {E}(t_{n-1})}{max(\mathcal {E})}\) (Fig. 2b).

Decreasing of the total energy over time by both space and time dG methods with different degrees of discontinuity. a Total loss accumulated w.r.t. time; b incremental loss w.r.t. time

Full size image

Figure 2a leads to a first conclusion: with only one physical interface, there is no significant influence of the degree of discontinuity, i.e. of a, on the numerical damping rate over time. Moreover, with our choice of discretization parameters \((h^E, \Delta t)\), i.e. \(CFL_{sDG} \approx C_{tDG-damping} \approx 1\), same levels of damping are observed for the two methods.

However, detailed analyses of the energy decreasing between two successive time steps (Figs. 2b,  3a) show that, at the physical interface, perturbations are observed only for the time dG method, which increase with decreasing value of a, while the space dG method shows a perfect behavior at the physical interface. Otherwise, due to the weak verification of the free boundary conditions, the energy damping of both methods is perturbed and such perturbations do not depend on a and are much more important for the space dG method than for the time dG method (Figs. 2b,  3b). It will be shown in the next section that such perturbations may become problematic in the 2D case.

a At the physical interface, influence of a on the energy damping by the time dG method; b influence of boundary conditions on the energy damping

Full size image

In conclusion, compared to the time dG method, the behavior of the space dG method across the physical interface is very interesting especially in cases involving a large number of physical interfaces. But future studies are necessary to understand and improve its behavior at external boundaries.

Secondly, influences of space–time discretization parameters \((h^E, \Delta t)\) are investigated. For both methods, the energy damping increases rapidly with mesh coarsening in space and in time (Fig. 4a). As expected, the energy damping of the space dG method depends principally on the element size in space (Fig. 4b), while the one of the time dG method depends principally on the time step (Fig. 4c).

Influence of \((h^E =\frac{\lambda _{min}}{N_e}, \Delta t)\) on the numerical damping for a both space and time dG methods with \(CFL_{sDG} \approx C_{tDG-damping} \approx 1\); b for the space dG method; c for the time dG method

Full size image

2D case with physical interfaces

For the 2D case, the studied 2D rectangular elastic domain \(\Omega \) is presented in Fig. 5a and the width of the whole domain is \(L_x = 5\;\text {mm}\). A similar but four times larger (in each space dimension) domain was already considered in [17]. The following four cases are simulated:

• (1) Homogeneous isotropic case

  In this case, the density, the Young’s moduli and the Poisson’s ratio of the elastic material are \(\rho = 4428\;\text {kg}\,\text {m}^{-3}\), \(E_1 = 73.95\;\text {GPa}\) and \(\nu = 0.395\). The phase velocities of the longitudinal and the transverse waves are respectively \(c_{L1} = 5879\;\text {m}\,\text {s}^{-1}\) and \(c_{T1} = 2446\;\text {m}\,\text {s}^{-1}\).

• (2) Heterogeneous isotropic case

  In this case, the domain \(\Omega = \Omega _1 \cup \Omega _2\) contains two vertical interfaces as shown in Fig. 5a. Only the Young’s modulus differs in \(\Omega _1\) and \(\Omega _2\) with \(E_1\) already given in the case (1) and \(E_2 = (1.5)^2 E_1\). The phase velocities of the longitudinal and the transverse waves in \(\Omega _2\) are respectively \(c_{L2} = 1.5\,c_{L1} = 8819\;\text {m}\,\text {s}^{-1}\) and \(c_{T2} = 1.5\,c_{T1} = 3670\;\text {m}\,\text {s}^{-1}\)

• (3) Homogeneous anisotropic case

  This case corresponds to a monocrystalline material, a single crystal of cubic symmetry from titanium alloy. The Euler angles of its crystallographic orientation are equal to \((0^{\circ }, 0^{\circ }, 0^{\circ })\). The elastic moduli of the Hooke tensor \(\varvec{C}\) in the local material anisotropic basis \((\varvec{a}_1, \varvec{a}_2, \varvec{a}_3)\) of the single crystal are: \(C_{iiii} = 134.0\; \text {GPa}\), \(C_{iijj, i\ne j} = 110.0\;\text {GPa}\) and \(C_{ijij, i \ne j} = 36.0 \;\text {GPa}\). The density is the same as the isotropic material given in the cases (1) and (2).

• (4) Heterogeneous polycrystalline case

  This case corresponds to a single-phase polycrystalline material with 85 elliptic grains of size \(480\;\upmu \text {m} \times 240\;\upmu \text {m}\) (Fig. 5b). Each grain has the same Hooke tensor \(\varvec{C}\) given in the case (3) but a randomly chosen crystallographic orientation. The layer surrounding the 85 grains is the single crystal of the case (3). The maximum phase velocity of quasi longitudinal waves is \(c_{qL, max} = 6123\;\text {m}\,\text {s}^{-1}\).

Studied 2D rectangular elastic domain. a Case with two vertical physical interfaces; b case of polycrystalline material with elliptic grains

Full size image

Elastic waves are generated by applying a pressure loading \(g(x,t) \varvec{e}_y\) in the vertical direction \(\varvec{e}_y\) on an extremely short segment \(L_e\) of length \(50\;\upmu \text {m}\) at the center of the external top boundary of \(\Omega \), while free surface boundary conditions are prescribed on the remaining parts (Fig. 5a). The pressure \(g(x,t) \varvec{e}_y\) is uniformly distributed in space and its period of ricker signal in time is \(T_r = 0.64\; \upmu \text {s}\), resulting in \(f_{max} = 3.1\;\text {MHz}\) and \(f_c = 7.8\;\text {MHz}\).

Bilinear square finite elements of size \(h^E = 25 \;\upmu \text {m}\) are used to discretize the whole domain \(\Omega \). For the two first cases, this choice corresponds to respectively 12 and 18 elements in the shortest transverse wavelengths of interest \(\ell _{T1}(f_c) = 314\, \upmu \text {m}\) in \(\Omega _1\) and \(\ell _{T2}(f_c) = 470\; \upmu \text {m}\) in \(\Omega _2\). For the polycrystalline material, it corresponds to 16 elements in the shortest quasi transverse wavelengths of interest \(\ell _{qT,min}(f_c) = 422\; \upmu \text {m}\). With respect to each elliptic grain in the polycrystalline case, it corresponds to about 13 elements in the major axis and about 10 element in the minor axis. This choice of element size is based on our previous studies on the mesh convergence of the time dG solver in terms of the attenuation and scattering (noise level) coefficients in the case of polycrystalline materials [21, 22]. In particular, it has been shown that the mesh convergence depends not only on the wavelength to element size ratio but also on the grain size to element size ratio.

The time step used for the space dG method is \(\Delta t_{sDG} = 0.533\,\text {ns}\) and the one for the time dG method is \(\Delta t_{tDG} = 3 \Delta t_{sDG} = 1.6\,\text {ns}\), leading to the value of \(CFL_{sDG} = C_{tDG-damping}\) equal to 0.37 for the case (1), to 0.56 for the case (2) and to 0.42 for the last two cases.

Evolution over time of the total, the kinetic and the elastic energies calculated by both space and time dG methods are compared in Figs. 6,  7. Both methods give very close curves for the kinetic energy, which is coherent with de comparison presented in [17] done only with the data of the velocity unknown field \(\varvec{v}_h\).

Energy evolution over time. a Comparison between both space and time dG methods in the heterogeneous isotropic case; b comparison between homogenous and heterogeneous isotropic cases

Full size image

Energy evolution over time. a Comparison between both space and time dG methods in the polycrystalline case; b comparison between monocrystalline and polycrystalline cases

Full size image

However, as the space dG method implemented in the present work uses the strategy of weak verification of the Neumann boundary conditions, the elastic energy calculated by the space dG method is highly perturbed by the non verification of those boundary conditions. Indeed, according to (40), for example, when the free Neumann boundary conditions are not exactly verified, the quantity \((\varvec{v}_h, \varvec{g})_{\partial \Omega _N}\), with \(\varvec{g} = \varvec{\sigma }_h\cdot \varvec{n}\), should be added in the energy evolution. In Fig. 6b, it is shown that this perturbation increases in the domain \(\Omega _2\) with larger elastic moduli. In Fig. 7b, it is shown that this perturbation highly increases in the polycrystalline case, probably because the grains constantly reflect waves in the boundary layer and thus reinforces the perturbation. Future studies are therefore necessary to improve the stress field solution of the space dG method with respect to the Neumann boundary condition, by adopting for example the strategy of strong verification of them. Afterward, quantitative comparison of the damping behavior between both space and dG methods would be performed.

Comparisons and analyses of the time discontinuous space–time Galerkin method and the space discontinuous Galerkin method applied to the elastic wave propagation in anisotropic and heterogeneous media were presented. More particularly, the stability of both methods was investigated. It has been shown that the mechanism to ensure the stability of both methods is similar, in the sense that it is based on the energies dissipated in time or space discontinuities of unknown fields, but essentially different. Indeed, contrary to the stability of the time dG method that is unconditional whatever the space and time discretizations used, the stability of the space dG method is proved for its variational formulation discretized in space but not yet discretized in time. When a time discretization scheme is applied, its own stability condition should also be taken into account.

As the main contribution, a new result of stability of the velocity–stress space dG method was given. The stability of the space dG method was proved at the continuous level of time in the case of elastic media with continuous properties and the case of isotropic elastic media with discontinuous properties, whatever the space dimension. However, the proof of the stability was given only in the 2D case for anisotropic elastic media with discontinuous properties under a sufficient condition. Future theoretical and numerical investigations are necessary to better understand this condition. However, a first analysis suggests that it should be satisfied for relatively high degrees of discontinuity. Otherwise, the existence of such a sufficient condition in the 2D case suggests the existence of a similar condition in the 3D case. It should be important to develop it even if the proof in the 3D case would be very complicated.

Otherwise, the stability result of the space dG method was proved with the use of the exact numerical fluxes. It allows explaining previously observed phenomena [17]: when the degree of discontinuity at physical interfaces is high, the use of approximate numerical fluxes leads to an instability problem, which can be considered as intrinsic to the space dG method, in the sense that it cannot be resolved by simply reducing the time step used by the time discretization scheme.

Numerical applications to different elastic media with physical interfaces were presented and the numerical damping behavior of both space and time dG methods was investigated. In the 1D case, the quantitative comparison between the two methods was complete and has shown, at physical interfaces, the perfect behavior of the space dG method using the exact numerical fluxes. In the 2D case, the necessity to well satisfy the Neumann boundary conditions by the stress field calculated by the space dG method was highlighted. Future studies should be done to complete the quantitative comparison between the two dG methods and to assess the convergence rate of the space dG solver in terms of the energy and L2 norms and also of the physical quantities of interest, such as the attenuation and noise levels in the case of polycrystalline materials.

1. It is worth noticing that the space dG method considered is usually called discontinuous Galerkin method. The word “space” is added in the present work in order to distinguish it from the considered time dG method.

Authors’ contributions

The author read and approved the final manuscript.

Acknowledgements

This study was the continuation of a work granted by the ANR (Agence Nationale de la Recherche)-MAPIE project, ANR-13-MONU-002.

Competing interests

The author declares no competing interests.

Availability of data and materials

All numerical setting data are presented in the article.

Authors and Affiliations

1. Laboratory MSSMat (UMR8579), CNRS, CentraleSupélec, Université Paris-Saclay, 8-10 rue Joliot-Curie, 91190, Gif-sur-Yvette, France

   Bing Tie

Corresponding author

Correspondence to Bing Tie.

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Appendix A. Proof of Theorem 4

Firstly, the following Lemma can be proved.

Lemma 2

Let us consider \(\Phi _\pm (C_1, C_2, \Delta , \zeta ; x, y)\) two second order polynomial functions of (x, y) defined as follows:

with \(C_1 \ge 0\), \(C_2 \ge 0\). Under the following condition:

they can be written in the following form:

where \(\beta _\pm (C_1, C_2, \Delta ) \) are two real numbers defined by (51c) in Theorem  4.

Proof The proof is given for \(\Phi _-\), as it can be done in the same way for \(\Phi _+\). The idea is to find an appropriate quadruplet \((a_1, b_1, a_2, b_2)\) so that the following decomposition of \(\Phi _-\) holds:

It is obvious that the quadruplet should verify the following equations:

with \(r_1\) and \(r_2\) defined as

Eliminating \((b_1, b_2)\), the first two equations of (58) become

Taking into account the symmetry between x and y, an intuitive and reasonable choice is \(r_1 = r_2 = \dfrac{1}{4}\) and \(a^2_1 = a^2_2\), then finally both equations of (60) become

It is obvious that, for both solutions of (61) to be real, the condition (55) should be verified. Then, together with (59), the following results are straightforward:

Putting (62) into (57) gives rise to the final result to prove. \(\square \)

Then the proof of Theorem 4 is organized in eight parts.

Proof of Theorem 4

• (1) Inverse of the system (26)

  To get \(\{\tilde{\alpha }_k, \tilde{\alpha }'_{k}\}_{k=1, 2}\), with \(k=1\) and 2 corresponding respectively to qL and qT, we need to calculate the inverse of the \(2d \times 2d\) matrix, denoted as \([\tilde{R}]\), of the system of linear equations (26). The inverse matrix is:

  $$\begin{aligned}{}[\tilde{R}]^{-1} = \begin{bmatrix} [\tilde{D}]&[\tilde{H}] \\ [\tilde{H}']&[\tilde{D}'] \end{bmatrix} = \begin{bmatrix} \tilde{D}_{11}&0&0&-\tilde{B}_{12}\tilde{D}_{11} \\ 0&\tilde{D}_{22}&-\tilde{B}_{21}D_{22}&0 \\ 0&-\tilde{B}'_{12}\tilde{D}'_{11}&\tilde{D}'_{11}&0 \\ -\tilde{B}'_{21}\tilde{D}'_{22}&0&0&\tilde{D}'_{22} \\ \end{bmatrix} \end{aligned}$$

  (63)

  with

  $$\begin{aligned} \tilde{D}_{11} = \tilde{D}'_{22} = \dfrac{1}{1-\tilde{B}_{12}\tilde{B}'_{21}} \; , \; \tilde{D}_{22} = \tilde{D}'_{11} = \dfrac{1}{1-\tilde{B}_{21}\tilde{B}'_{12}} \end{aligned}$$

  (64)

  We recall that four orthonormal bases are used: \((\varvec{n}, \varvec{t})\), \((\varvec{n}', \varvec{t}')\), \((\varvec{\gamma }_{\varvec{n},1}, \varvec{\gamma }_{\varvec{n},2})\) and \((\varvec{\gamma }'_{\varvec{n}',1}, \varvec{\gamma }'_{\varvec{n}',2})\) with \(\varvec{e}_3=\varvec{n}\wedge \varvec{t} = \varvec{n}'\wedge \varvec{t}' = \varvec{\gamma }_{\varvec{n},1}\wedge \varvec{\gamma }_{\varvec{n},2} = \varvec{\gamma }'_{\varvec{n}',1}\wedge \varvec{\gamma }'_{\varvec{n}',2}\), then it can be shown that the following equations hold:

  $$\begin{aligned} \zeta&= \varvec{\gamma }_{\varvec{n},1}\cdot \varvec{\gamma }'_{\varvec{n}',2} = - \varvec{\gamma }_{\varvec{n},2}\cdot \varvec{\gamma }'_{\varvec{n}',1} = \gamma _{\varvec{n},1\varvec{n}} \gamma '_{\varvec{n}',1\varvec{t}'} - \gamma _{\varvec{n},1\varvec{t}} \gamma '_{\varvec{n}',1\varvec{n}'} \end{aligned}$$

  (65a)
  $$\begin{aligned} \tilde{B}_{12}&= -\dfrac{C^-_{z,1}}{2} \dfrac{\delta z^{-'}_{\varvec{n}',12}}{z^{-'}_{\varvec{n}',1}} \zeta \;,\; \tilde{B}'_{12} = \dfrac{C^{-'}_{z,1}}{2} \dfrac{\delta z^-_{\varvec{n},12}}{z^-_{\varvec{n},1}} \zeta = -\dfrac{\delta z^-_{\varvec{n},12}}{\delta z^{-'}_{\varvec{n}',12}} \tilde{B}_{12} \end{aligned}$$

  (65b)
  $$\begin{aligned} \tilde{B}_{21}&= -\dfrac{C^-_{z,2}}{2} \dfrac{\delta z^{-'}_{\varvec{n}',12}}{z^{-'}_{\varvec{n}',2}} \zeta \;,\; \tilde{B}'_{21} = \dfrac{C^{-'}_{z,2}}{2} \dfrac{\delta z^-_{\varvec{n},12}}{z^-_{\varvec{n},2}} \zeta = -\dfrac{\delta z^-_{\varvec{n},12}}{\delta z^{-'}_{\varvec{n}',12}} \tilde{B}_{21} \end{aligned}$$

  (65c)

  with \(\gamma _{\varvec{n},1\varvec{n}} = \varvec{\gamma }_{\varvec{n},1}\cdot \varvec{n}\), \(\gamma _{\varvec{n},1\varvec{t}} = \varvec{\gamma }_{\varvec{n},1}\cdot \varvec{t}\), \(\gamma '_{\varvec{n}',1\varvec{n}'} = \varvec{\gamma }'_{\varvec{n}',1}\cdot \varvec{n}'\) and \(\gamma _{\varvec{n}',1\varvec{t}'} = \varvec{\gamma }_{\varvec{n}',1}\cdot \varvec{t}'\). We note that \(\gamma _{\varvec{n},2\varvec{n}} = - \gamma _{\varvec{n},1\varvec{t}}\) and \(\gamma _{\varvec{n},2\varvec{t}} = \gamma _{\varvec{n},1\varvec{n}}\).

  Otherwise, we recall that \(\lambda ^-_{\varvec{n},1} \le \lambda ^-_{\varvec{n},2}\) and \(\lambda ^{-'}_{\varvec{n}',1} \le \lambda ^{-'}_{\varvec{n}',2}\), hence we get:

  $$\begin{aligned} \tilde{D}_{11} = \tilde{D}_{22} = \tilde{D}'_{11} = \tilde{D}'_{22} =\dfrac{1}{1+\tilde{B}_{12}\tilde{B}_{21}\dfrac{\delta z^-_{\varvec{n},12}}{\delta z^{-'}_{\varvec{n}',12}}} > 0 \end{aligned}$$

  (66)

• (2) Decomposition of \(\partial _t\mathcal {E}_\Omega \) into different parts to analyze

  According to (26), we have:

  $$\begin{aligned} \tilde{\alpha }_k&= \tilde{D}_{kl} \,\tilde{\,}\varvec{L}^-_{\varvec{n},l} \cdot [\varvec{U}_h] - \tilde{H}_{kl} \, \tilde{\,}\varvec{L}^{-'}_{\varvec{n}',l} \cdot [\varvec{U}_h] \end{aligned}$$

  (67a)
  $$\begin{aligned} \tilde{\alpha }'_k&= -\tilde{D}'_{kl} \, \tilde{\,}\varvec{L}^{-'}_{\varvec{n}',l} \cdot [\varvec{U}_h] + \tilde{H}'_{kl} \, \tilde{\,}\varvec{L}^-_{\varvec{n},l} \cdot [\varvec{U}_h] \end{aligned}$$

  (67b)

  Using Lemma 1, we get:

  $$\begin{aligned} \partial _t\mathcal {E}_\Omega = \sum _{\Gamma _{EE'}} \dfrac{\tilde{D}_{11}}{2} \big ([\varvec{U}_h], (\varvec{P}_{EE'} + \varvec{Q}_{EE'} + \varvec{Q}'_{EE'})\cdot [\varvec{U}_h] \big )_{\Gamma _{EE'}} \end{aligned}$$

  (68)

  with

  $$\begin{aligned} \varvec{P}_{EE'}&= z^-_{\varvec{n},1} \varvec{L}^-_{\varvec{n},1} \otimes \tilde{\,}\varvec{L}^-_{\varvec{n},1} + z^-_{\varvec{n},2} \varvec{L}^-_{\varvec{n},2} \otimes \tilde{\,}\varvec{L}^-_{\varvec{n},2} \nonumber \\&\qquad + z^{-'}_{\varvec{n}',1} \varvec{L}^{-'}_{\varvec{n}',1} \otimes \tilde{\,}\varvec{L}^{-'}_{\varvec{n}',1} + z^{-'}_{\varvec{n}',2} \varvec{L}^{-'}_{\varvec{n}',2} \otimes \tilde{\,}\varvec{L}^{-'}_{\varvec{n}',2} \end{aligned}$$

  (69a)
  $$\begin{aligned} \varvec{Q}_{EE'}&= \tilde{B}_{12} z^-_{\varvec{n},1} \varvec{L}^-_{\varvec{n},1} \otimes \tilde{\,}\varvec{L}^{-'}_{\varvec{n}',2} + \tilde{B}_{21} z^-_{\varvec{n},2} \varvec{L}^-_{\varvec{n},2} \otimes \tilde{\,}\varvec{L}^{-'}_{\varvec{n}',1} \end{aligned}$$

  (69b)
  $$\begin{aligned} \varvec{Q}'_{EE'}&= \tilde{B}'_{12} z^{-'}_{\varvec{n},1} \varvec{L}^{-'}_{\varvec{n}',1} \otimes \tilde{\,}\varvec{L}^-_{\varvec{n},2} + \tilde{B}'_{21} z^{-'}_{\varvec{n}',2} \varvec{L}^{-'}_{\varvec{n}',2} \otimes \tilde{\,} \varvec{L}^-_{\varvec{n},1} \end{aligned}$$

  (69c)

• (3) Calculation of the term concerning \(\varvec{P}_{EE'}\)

  According to (15) and (29) the definitions of \((\varvec{L}^-_{\varvec{n},k}, \varvec{L}^{-'}_{\varvec{n}',k})\) and \((\tilde{\,}\varvec{L}^-_{\varvec{n},k}, \tilde{\,}\varvec{L}^{-'}_{\varvec{n}',k})\) and recalling that \(\varvec{w}_{\varvec{n}, k} = \frac{1}{\sqrt{2}}\varvec{\gamma }_{\varvec{n},k}\), we get:

  $$\begin{aligned}&[\varvec{U}_h]\cdot \varvec{P}_{EE'}\cdot [\varvec{U}_h] \nonumber \\ =&-\sum _k \dfrac{\overline{|z^-_{\varvec{n},k}|}^{R}}{2}\bigg ([\varvec{v}_h]\cdot \varvec{\gamma }_{\varvec{n},k} + \dfrac{[\varvec{\sigma }_h\cdot \varvec{n}]\cdot \varvec{\gamma }_{\varvec{n},k}}{|z^-_{\varvec{n},k}|} \bigg )\bigg ([\varvec{v}_h]\cdot \varvec{\gamma }_{\varvec{n},k} + \dfrac{[\varvec{\sigma }_h\cdot \varvec{n}]\cdot \varvec{\gamma }_{\varvec{n},k}}{|z^{-'}_{\varvec{n}',k}|} \bigg ) \nonumber \\&-\sum _k\dfrac{\overline{|z^-_{\varvec{n},k}|}^{R}}{2}\bigg ([\varvec{v}_h]\cdot \varvec{\gamma }'_{\varvec{n}',k} - \dfrac{[\varvec{\sigma }_h\cdot \varvec{n}]\cdot \varvec{\gamma }'_{\varvec{n}',k}}{|z^{-'}_{\varvec{n}',k}|} \bigg )\bigg ([\varvec{v}_h]\cdot \varvec{\gamma }'_{\varvec{n}',k} - \dfrac{[\varvec{\sigma }_h\cdot \varvec{n}]\cdot \varvec{\gamma }'_{\varvec{n}',k}}{|z^-_{\varvec{n},k}|} \bigg ) \end{aligned}$$

  (70a)
  $$\begin{aligned} =&-\dfrac{1}{2}\sum _k \overline{|z^-_{\varvec{n},k}|}^{R} \bigg ( ([\varvec{v}_h]\cdot \varvec{\gamma }_{\varvec{n},k})^2 + ([\varvec{v}_h]\cdot \varvec{\gamma }'_{\varvec{n}',k})^2 \bigg ) \nonumber \\&- \dfrac{1}{2} \sum _k \dfrac{ ([\varvec{\sigma }_h \cdot \varvec{n}] \cdot \varvec{\gamma }_{\varvec{n},k})^2 + ([\varvec{\sigma }_h \cdot \varvec{n}] \cdot \varvec{\gamma }'_{\varvec{n}',k})^2}{\overline{|z^-_{\varvec{n},k}|}^{V}} \nonumber \\&+ ([\varvec{v}_h]\otimes [\varvec{\sigma }_h\cdot \varvec{n}]):\bigg (\sum _k \varvec{\gamma }_{\varvec{n},k} \otimes \varvec{\gamma }_{\varvec{n},k} - \sum _k \varvec{\gamma }'_{\varvec{n}',k} \otimes \varvec{\gamma }'_{\varvec{n}',k}\bigg ) \end{aligned}$$

  (70b)

  We remark that the last term of (70b) equals to zero as \(\sum _k \varvec{\gamma }_{\varvec{n},k} \otimes \varvec{\gamma }_{\varvec{n},k} = \sum _k \varvec{\gamma }'_{\varvec{n}',k} \otimes \varvec{\gamma }'_{\varvec{n}',k} = \varvec{I}\). Finally (70) is written in the following form by decomposing it into two independent parts, related respectively to \([\varvec{v}_h]\) and \([\varvec{\sigma }_h\cdot \varvec{n}]\):

  $$\begin{aligned}{}[\varvec{U}_h]\cdot \varvec{P}_{EE'}\cdot [\varvec{U}_h] = P_{[\varvec{v}_h]} + P_{[\varvec{\sigma }_h\cdot \varvec{n}]} \end{aligned}$$

  (71)

• (4) Calculation of the term concerning \(\varvec{Q}_{EE'}\)

  According to (15), (29) and (28), the definitions of \((\varvec{L}^-_{\varvec{n},k}, \varvec{L}^{-'}_{\varvec{n}',k})\), \((\tilde{\,}\varvec{L}^-_{\varvec{n},k}, \tilde{\,}\varvec{L}^{-'}_{\varvec{n}',k})\) and \((C^-_{z,k}, C^{-'}_{z,k})\) and to (65), we get:

  $$\begin{aligned} \begin{array}{rcl} &{} &{}[\varvec{U}_h]\cdot \varvec{Q}_{EE'}\cdot [\varvec{U}_h] \\ &{}= &{}- \dfrac{\overline{|z^-_{\varvec{n},1}|}^{R}\overline{|z^-_{\varvec{n},2}|}^{R}}{2} \big |\dfrac{\delta z^{-'}_{\varvec{n}',12}}{z^{-'}_{\varvec{n}',1}z^{-'}_{\varvec{n}',2}}\big | \zeta \\ &{}&{} \bigg ( ([\varvec{v}_h] \otimes [\varvec{v}_h] - \dfrac{[\varvec{\sigma }_h\cdot \varvec{n}] \otimes [\varvec{\sigma }_h\cdot \varvec{n}]}{|z^{-}_{\varvec{n},1}||z^{-}_{\varvec{n},2}|}\bigg ):(\varvec{\gamma }_{\varvec{n},1}\otimes \varvec{\gamma }'_{\varvec{n}',2} + \varvec{\gamma }_{\varvec{n},2}\otimes \varvec{\gamma }'_{\varvec{n}',1}) + C_Q \big ) \end{array} \end{aligned}$$

  (72)

  with

  $$\begin{aligned} \begin{array}{rcl} C_Q &{}= &{}\dfrac{[\varvec{v}_h] \otimes [\varvec{\sigma }_h\cdot \varvec{n}]}{|z^-_{\varvec{n},1}|} (\varvec{\gamma }'_{\varvec{n}',2}\otimes \varvec{\gamma }_{\varvec{n},1} - \varvec{\gamma }_{\varvec{n},2}\otimes \varvec{\gamma }'_{\varvec{n}',1}) \\ &{} \qquad +&{} \dfrac{[\varvec{v}_h] \otimes [\varvec{\sigma }_h\cdot \varvec{n}]}{|z^-_{\varvec{n},2}|} (- \varvec{\gamma }_{\varvec{n},1} \otimes \varvec{\gamma }'_{\varvec{n}',2} + \varvec{\gamma }'_{\varvec{n}',1} \otimes \varvec{\gamma }_{\varvec{n},2}) \big ) \end{array} \end{aligned}$$

  (73)

  Then by taking into account that:

  $$\begin{aligned} (\varvec{\gamma }'_{\varvec{n}',2}\otimes \varvec{\gamma }_{\varvec{n},1} - \varvec{\gamma }_{\varvec{n},2}\otimes \varvec{\gamma }'_{\varvec{n}',1}) = (\gamma _{\varvec{n},1\varvec{n}} \gamma '_{\varvec{n}',1\varvec{t}'} - \gamma _{\varvec{n},1\varvec{t}} \gamma '_{\varvec{n}',1\varvec{n}'} )(\varvec{n}\otimes \varvec{n}+ \varvec{t}\otimes \varvec{t}) = \zeta \varvec{I} \end{aligned}$$

  (74)

  it can be proved that:

  $$\begin{aligned} \begin{array}{l} [\varvec{U}_h]\cdot \varvec{Q}_{EE'}\cdot [\varvec{U}_h] \\ \quad =- \dfrac{\overline{|z^-_{\varvec{n},1}|}^{R}\overline{|z^-_{\varvec{n},2}|}^{R}}{2} \big |\dfrac{\delta z^{-'}_{\varvec{n}',12}}{z^{-'}_{\varvec{n}',1}z^{-'}_{\varvec{n}',2}}\big | \zeta \left( [\varvec{v}_h] \otimes [\varvec{v}_h] - \dfrac{[\varvec{\sigma }_h\cdot \varvec{n}] \otimes [\varvec{\sigma }_h\cdot \varvec{n}]}{|z^{-}_{\varvec{n},1}||z^{-}_{\varvec{n},2}|}\right) :\varvec{A}_Q\\ \qquad +\; \dfrac{|\delta z^-_{\varvec{n},12}||\delta z^{-'}_{\varvec{n}',12}| \zeta ^2}{2\overline{|z^-_{\varvec{n},1}|}^{V}\overline{|z^-_{\varvec{n},2}|}^{V}} [\varvec{\sigma }_h\cdot \varvec{n}] \cdot [\varvec{v}_h] \end{array} \end{aligned}$$

  (75)

  In (75), the following notation is adopted:

  $$\begin{aligned} \varvec{A}_Q = \varvec{\gamma }_{\varvec{n},1}\otimes \varvec{\gamma }'_{\varvec{n}',2} + \varvec{\gamma }_{\varvec{n},2}\otimes \varvec{\gamma }'_{\varvec{n}',1} \end{aligned}$$

  (76)

  Otherwise, it can be proved that \(\varvec{A}_Q\) is symmetric, as we have

  $$\begin{aligned} \begin{array}{rcl} \varvec{A}_Q &{}=&{} (\gamma _{\varvec{n},1\varvec{n}} \gamma '_{\varvec{n}',1\varvec{t}'} + \gamma _{\varvec{n},1\varvec{t}} \gamma '_{\varvec{n}',1\varvec{n}'} )(\varvec{n}\otimes \varvec{n}- \varvec{t}\otimes \varvec{t}) \\ &{}&{}-\; (\gamma _{\varvec{n},1\varvec{n}} \gamma '_{\varvec{n}',1\varvec{n}'} - \gamma _{\varvec{n},1\varvec{t}'} \gamma '_{\varvec{n}',1\varvec{n}'} )(\varvec{n}\otimes \varvec{t} + \varvec{t}\otimes \varvec{n}) \end{array} \end{aligned}$$

  (77)

• (5) Calculation of the term concerning \(\varvec{Q}'_{EE'}\)

  In the same way as for \(\varvec{Q}_{EE'}\), we get:

  $$\begin{aligned}&[\varvec{U}_h]\cdot \varvec{Q}'_{EE'}\cdot [\varvec{U}_h] \nonumber \\&\quad = \dfrac{\overline{|z^-_{\varvec{n},1}|}^{R}\overline{|z^-_{\varvec{n},2}|}^{R}}{2} \left| \dfrac{\delta z^-_{\varvec{n},12}}{z^-_{\varvec{n},1}z^-_{\varvec{n},2}} \right| \zeta \left( [\varvec{v}_h] \otimes [\varvec{v}_h] - \dfrac{[\varvec{\sigma }_h\cdot \varvec{n}] \otimes [\varvec{\sigma }_h\cdot \varvec{n}]}{|z^{-'}_{\varvec{n}',1}||z^{-'}_{\varvec{n}',2}|}\right) :\varvec{A}_{Q'} \nonumber \\&\qquad - \dfrac{|\delta z^-_{\varvec{n},12}||\delta z^{-'}_{\varvec{n}',12}| \zeta ^2}{2\overline{|z^-_{\varvec{n},1}|}^{V}\overline{|z^-_{\varvec{n},2}|}^{V}} [\varvec{\sigma }_h\cdot \varvec{n}] \cdot [\varvec{v}_h] \end{aligned}$$

  (78)

  with

  $$\begin{aligned} \varvec{A}'_Q = \varvec{\gamma }'_{\varvec{n}',1}\otimes \varvec{\gamma }_{\varvec{n},2} + \varvec{\gamma }'_{\varvec{n}',2}\otimes \varvec{\gamma }_{\varvec{n},1} \end{aligned}$$

  (79)

  It can be proved that \(\varvec{A}'_Q = (\varvec{A}_Q)^T = \varvec{A}_Q\). Otherwise, we remark that the sum of the two last terms in (75) and (78) is equal to zero.

• (6) Calculation of the two terms concerning \(\varvec{Q}_{EE'}\) and \(\varvec{Q}'_{EE'}\)

  The two terms that need further analysis are:

  $$\begin{aligned} \begin{array}{lll} &{}&{}[\varvec{U}_h]\cdot (\varvec{Q}_{EE'}+\varvec{Q}'_{EE'})\cdot [\varvec{U}_h] \\ &{}&{}\quad = - \;\dfrac{\overline{|z^-_{\varvec{n},1}|}^{R}\overline{|z^-_{\varvec{n},2}|}^{R}}{2} \left| \dfrac{\delta z^{-'}_{\varvec{n}',12}}{z^{-'}_{\varvec{n}',1}z^{-'}_{\varvec{n}',2}}\right| \zeta \left( [\varvec{v}_h] \otimes [\varvec{v}_h] - \dfrac{[\varvec{\sigma }_h\cdot \varvec{n}] \otimes [\varvec{\sigma }_h\cdot \varvec{n}]}{|z^{-}_{\varvec{n},1}||z^{-}_{\varvec{n},2}|}\right) :\varvec{A}_Q \\ &{}&{}\qquad + \; \dfrac{\overline{|z^-_{\varvec{n},1}|}^{R}\overline{|z^-_{\varvec{n},2}|}^{R}}{2} \left| \dfrac{\delta z^-_{\varvec{n},12}}{z^-_{\varvec{n},1}z^-_{\varvec{n},2}}\right| \zeta \left( [\varvec{v}_h] \otimes [\varvec{v}_h] - \dfrac{[\varvec{\sigma }_h\cdot \varvec{n}] \otimes [\varvec{\sigma }_h\cdot \varvec{n}]}{|z^{-'}_{\varvec{n}',1}||z^{-'}_{\varvec{n}',2}|}\right) :\varvec{A}_{Q'} \\ &{}&{} \quad \equiv Q_{[\varvec{v}_h]} + Q_{[\varvec{\sigma }_h\cdot \varvec{n}]} \end{array} \end{aligned}$$

  (80)

• (7) Calculation of the terms \(P_{[\varvec{v}_h]}\) and \(Q_{[\varvec{v}_h]}\)

  Using (51a), (54), (70), (71), (80) the definition of the polynomial functions \(\Phi _\pm \) given in Lemma 2 and (50) the definition of \(G(C_1, C_2, \Delta ; \varvec{w})\) in Theorem 4, it can be shown that:

  $$\begin{aligned} P_{[\varvec{v}_h]} + Q_{[\varvec{v}_h]}= & {} -\dfrac{1}{2} \big (\Phi _-(\overline{|z^-_{\varvec{n},1}|}^{R}, \overline{|z^-_{\varvec{n},2}|}^{R}, \Delta ^R_{12}, \zeta ; [\varvec{v}_h] \cdot \varvec{\gamma }_{\varvec{n},1}, [\varvec{v}_h] \cdot \varvec{\gamma }'_{\varvec{n}',2}) \nonumber \\&+ \; \Phi _+(\overline{|z^-_{\varvec{n},2}|}^{R}, \overline{|z^-_{\varvec{n},1}|}^{R}, \Delta ^R_{12}, \zeta ; [\varvec{v}_h] \cdot \varvec{\gamma }_{\varvec{n},2}, [\varvec{v}_h] \cdot \varvec{\gamma }'_{\varvec{n}',1}) \big ) \nonumber \\= & {} - \dfrac{1}{4} G(\overline{|z^-_{\varvec{n},1}|}^{R}, \overline{|z^-_{\varvec{n},1}|}^{R}, \Delta ^R_{12}; [\varvec{v}_h]) \end{aligned}$$

  (81)

  with

  $$\begin{aligned} \Delta ^R_{12} \equiv \bigg |\dfrac{\delta z^-_{\varvec{n},12}}{z^-_{\varvec{n},1}z^-_{\varvec{n},2}}\bigg | - \bigg |\dfrac{\delta z^{-'}_{\varvec{n}',12}}{z^{-'}_{\varvec{n}',1}z^{-'}_{\varvec{n}',2}}\bigg | = \bigg | \dfrac{1}{z^-_{\varvec{n},1}} - \dfrac{1}{z^-_{\varvec{n},2}} \bigg | - \bigg | \dfrac{1}{z^{-'}_{\varvec{n}',1}} - \dfrac{1}{z^{-'}_{\varvec{n}',2}} \bigg | \end{aligned}$$

  (82)

  Furthermore, according to (55), it is easy to obtain (49a) the first part of the sufficient condition (49).

• (8) Calculation of the terms \(P_{[\varvec{\sigma }_h\cdot \varvec{n}]}\) and \(Q_{[\varvec{\sigma }_h\cdot \varvec{n}]}\)

  Using (51b), (54), (70), (71), (80) the definition of the polynomial functions \(\Phi _\pm \) given in Lemma 2 and (50) the definition of \(G(C_1, C_2, \Delta ; \varvec{w})\) in Theorem 4, it can be shown that:

  $$\begin{aligned} \begin{array}{rl} P_{[{\varvec{\sigma }_h\cdot \varvec{n}}]} + Q_{[{\varvec{\sigma }_h\cdot \varvec{n}}]} = -\dfrac{1}{2} &{} \Bigg (\Phi _+\left( \dfrac{1}{\overline{|z^-_{\varvec{n},1}|}^{V}}, \dfrac{1}{\overline{|z^-_{\varvec{n},2}|}^{V}}, \Delta ^V_{12}, \zeta ; [\varvec{\sigma }_h \cdot \varvec{n}] \cdot \varvec{\gamma }_{\varvec{n},1}, [\varvec{\sigma }_h \cdot \varvec{n}] \cdot \varvec{\gamma }'_{\varvec{n}',2}\right) \\ + &{} \Phi _-\left( \dfrac{1}{\overline{|z^-_{\varvec{n},2}|}^{V}}, \dfrac{1}{\overline{|z^-_{\varvec{n},1}|}^{V}}, \Delta ^V_{12}, \zeta ; [\varvec{\sigma }_h \cdot \varvec{n}] \cdot \varvec{\gamma }_{\varvec{n},2}, [\varvec{\sigma }_h \cdot \varvec{n}] \cdot \varvec{\gamma }'_{\varvec{n}',1}\right) \Bigg ) \\ = - \dfrac{1}{4} &{} G\left( \dfrac{1}{\overline{|z^-_{\varvec{n},1}|}^{{V}}}, \dfrac{1}{\overline{|z^-_{\varvec{n},1}|}^{V}}, \Delta ^V_{12}; [\varvec{\sigma }_h \cdot \varvec{n}] \right) \end{array} \end{aligned}$$

  (83)

  with

  $$\begin{aligned} \Delta ^V_{12} = |\delta z^-_{\varvec{n},12}| - |\delta z^{-'}_{\varvec{n}',12}| \end{aligned}$$

  (84)

  Furthermore, according to (55), it is easy to obtain (49b) the second part of the sufficient condition (49).

Finally, taking into account (68), (81) and (83) leads to (48) of Theorem 4. \(\square \)

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Reprints and permissions