FEM–BEM coupling for the thermoelastic wave equation with transparent boundary conditions in 3D

We consider the thermoelastic wave equation in three dimensions with transparent boundary conditions on a bounded, not necessarily convex domain. In order to solve this problem numerically, we introduce a coupling of the thermoelastic wave equation in the interior domain with time-dependent boundary integral equations. Here, we want to highlight that this type of problem differs from other wave-type problems that dealt with FEM–BEM coupling so far, e.g., the acoustic as well as the elastic wave equation, since our problem consists of coupled partial differential equations involving a vector-valued displacement field and a scalar-valued temperature field. This constitutes a nontrivial challenge which is solved in this paper. Our main focus is on a coercivity property of a Calderón operator for the thermoelastic wave equation in the Laplace domain, which is valid for all complex frequencies in a half-plane. Combining Laplace transform and energy techniques, this coercivity in the frequency domain is used to prove the stability of a fully discrete numerical method in the time domain. The considered numerical method couples finite elements and the leapfrog time-stepping in the interior with boundary elements and convolution quadrature on the boundary. Finally, we present error estimates for the semi- and full discretization.


Introduction and outline
In this paper, we study the thermoelastic wave equation in R 3 , which describes the interaction of the elastic behavior of a material and its temperature. Initial conditions and inhomogeneities are assumed to have support in a bounded, possibly non-convex domain, such that waves may re-enter the domain. Instead of considering this domain explicitly, we equip the wave equation with transparent boundary conditions, which are nonlocal in space and time. We achieve a stable numerical coupling of the interior and exterior problems by use of a so-called Calderón operator. Thus, it suffices to solve a problem which couples the thermoelastic wave equation to a boundary integral equation which represents the transmission conditions between the interior and exterior domain.
The kind of coupling we intend to use has been applied successfully for other wave-type equations, e.g., in [4] for the acoustic wave equation, in [8] for the elastodynamic wave equation, and in [15] for Maxwell's equations. Further, we want to mention the comparison of the acoustic and elastic FEM-BEM coupling in [10] as well as the implementation of the elastodynamic problem in [9]. However, the functional analytic setting is different as thermoelasticity has as variables a vector-valued 1 displacement field U and a scalar-valued temperature field T . One of the challenges is to include the equation involving the time

The thermoelastic wave equation
It is well known that elastic solids expand if their temperature rises. This interaction also works the other way around, i.e., elastic deformations influence the temperature. Such interactions are summarized in the term thermoelasticity. In a simple, linear, but non-stationary setting, an adequate description in three-dimensional Euclidean space is given by the system (see, e.g., [13], Eq. (3.1.13))

4)
T ( x, 0) = T 0 in R 3 , Here, · 2 is the Euclidean norm in R 3 , the index c.f. means the curl-free part of U , the index d.f. means the divergence-free part of U , and ∂ ∂r denotes a derivative in radial direction. In order to obtain a system in which each PDE is of second order with respect to time, we apply the variable transformation T = ∂ t Q such that Eqs. (2.1)-(2.5) become

12)
Q( x, 0) = Q 0 in R 3 , (2.13) (2.14) This transformation is essential for successfully applying several techniques in the proofs later on in this paper. We also want to highlight that this choice is not trivial and required a deep examination of the underlying problem. In contrast to other equations for which the numerical method we have in mind has been considered, we now have to deal with a vector-valued and a scalar-valued quantity. Moreover, our system involves terms which mix spatial and temporal derivatives. Note that due to the mixed derivative term κ T ref Δ∂ t Q( x, t) in Eq. (2.10), we have a third-order system here.

Transmission conditions between interior and exterior space
To account for bounded support of initial conditions and inhomogeneities, we split our problem: For the interior part of a bounded Lipschitz-domain Ω ⊂ R 3 which contains the supports of initial conditions and right-hand sides, find the solution U − for the system (2.9)-(2.14), whereas for the exterior part Ω + = R 3 \ Ω, where F = 0, G = 0 and initial conditions vanish, determine the solution U + of the homogenous version of (2.9)-(2.14). Both problems are coupled by transmission conditions on the boundary ∂Ω = Γ. Thus, we introduce Dirichlet as well as Neumann traces. We remind the reader that Neumann traces are related to co-normal derivatives, which may differ from normal derivatives. In our case, they even  (2.18) where n is the outer unit normal to Ω, is the elastic stress tensor. Here, we also introduced the elastic strain tensor (2.20) The thermoelastic stress tensor is Based on this, we finally state the transmission conditions We also extend the trace operators to Sobolev spaces in the usual way, such that

Laplace transform of the thermoelastic wave equation
We introduce the Laplace transforms 2 with the parameter s ∈ C. If we assume vanishing initial conditions for the moment, transferring the modified thermoelastic wave equations (2.9)-(2.10) to the Laplace domain yields This leads to the relation

Potentials and Calderón operator
The material of this section is taken from [19], Chapters 6 and 7, and included for the convenience of the reader. We introduce the trace operators For any linear partial differential equation, the single layer potential is given by and the double layer potential by The fundamental solution tensor ⇒ G s ( x, y) can be adapted from [13], Section 3.2(b), or [16], Chapter II, §3. We include it in Appendix A. The layer potentials and related traces fulfill the following mapping properties: The jumps of the traces over Γ are defined as which allows us to state the corresponding jump relations such that with the definitions Moreover, we define the operators where we also introduced the notation {{·}} for averages of inner and outer limits at the boundary. This allows us to compactly summarize the limit and jump relations ⎛ Moreover, we introduce the Calderón operator (3.20)

Coercivity results for the Calderón operator
In order to show stability of the suggested numerical model, we need to show that the Calderón operator satisfies a coercivity estimate, both in the Laplace as well as in the time domain. Our considerations here extend the results for the elastic wave equations in [8] as we consider a vector-valued partial differential equation (second order in time and space) which is coupled with a scalar-valued third-order PDE (second order in time), which stems from the transformation of the original problem as well as the corresponding boundary conditions.

Coercivity of the Laplace transformed Calderón operator B(s)
We need the following estimates which can be found in [11], p. 28, and [19], Theorem 10.2, respectively.

Lemma 1.
Let Ω be a domain with a Lipschitz boundary Γ. Then, it holds the trace inequality and Korn's second inequality With these estimates, we can show coercivity of the Calderón operator. In the following, ·, · Γ denotes the anti-duality in (H − 1 2 (Γ)) 4 × (H 1 2 (Γ)) 4 and (H 1 2 (Γ)) 4 × (H − 1 2 (Γ)) 4 . Furthermore, we define for for Re(s) > 0 and for all Proof. We start with the consideration of If we now substitute the definitions of γ N U and γ N Q from Eqs. (2.28) and (2.30), we can apply a modified version of Green's second identity for poroelasticity from [2], Theorem 4.3., to obtain where d min = min 2μ, λ, ρ, κ T ref , . Next, we take a look at the norms of the boundary densities. It is where we have also used Eqs. (4.1) and (2.9), the triangle inequality, and defined d 1,max := d 1 max(8μ 2 , 2λ 2 , where we used Eqs. (4.1) and (2.10), the triangle inequality, the relation which is obtained by using Korn's inequality (4.2), and By adding the estimates for the norms of the four boundary densities, we find where d max := 3 max{d 1,max , d 2,max , d 3 , d 4 }. By combining the final estimates in Eqs. (4.6) and (4.11), we achieve the desired result with c := dmin dmax .

Coercivity of the time-dependent Calderón operator B(∂ t )
In order to analyze the stability of our numerical scheme and find error estimates, we have to transfer the coercivity result for the Calderón operator from the Laplace domain to the time domain. This is directly possible due to an operator-valued variant of the classical Herglotz theorem and the convolution quadrature as shown in [4]. The following lemmas follow from this Herglotz theorem in a way similar to the corresponding results for the elastic wave equation derived in [8] and from the estimates where M depends only on σ = Re(s). These estimates can be derived in a similar way as in [3,17]. Hence, As in the proof of Lemma 2, we get The Lax-Milgram lemma in the form of [22] (Lemma 2.1.51 with the definition of ellipticity as in (2.43)) then gives With Eq. (4.6), we get (4.26) The Lax-Milgram lemma in the form of [22] (Lemma 2.1.51 with the definition of ellipticity as in (2.43)) then gives for any t end > 0 and for all is a short-hand notation for an integration with respect to time from 0 to t which corresponds to the factor 1 s in the Laplace domain. Lemma 3 can be extended to bound the time behavior of an energy.
where d end = min(t −2 end , t −5 end ) andĖ denotes the derivative of E with respect to t.

Remark 1.
Here, we want to highlight one of the major differences to the acoustic [4], elastic [8] and Maxwell's [15] cases. In both, Lemmas 3 and 4, we have integrals involving 4 . This leads to nontrivial challenges in the energy estimates later on, which are solved in this paper.

Lemma 5.
In the situation of Lemma 3, we have for t end = mΔt and 0 < Δt

33)
for all sequences where d end = min(t −2 end , t −5 end ) and c > 0 depends only on Δt 0 and tends to 1 as Δt 0 tends to zero.

Discretization
As we now have all the necessary theoretical results, let us discuss the discretization of our problem. We do this in two steps, first a semi-discretization in space and then a full discretization including space and time.

Variational formulation
Returning to the time domain, the second relation in Eq. (3.12) translates to the more explicit formulation which yields for the Calderón operator Together with Eqs. (2.9) and (2.10), this leads to the following formulation which is of first order in time: 3 , and χ U ∈ (H 1 2 (Γ)) 3 , we perform integration by parts to obtain a variational formulation. However, we split up some of the terms before performing integration by parts in order to end up with antisymmetries. As it should be clear from context, we only write (·, ·) to denote the L 2 -inner product in each of the spaces H 1 (Ω), (H 1 (Ω)) 3 , and (H 1 (Ω)) 3×3 . By observing that These equations are completed by the relations (5.14) which are valid due to the transmission conditions (2.22).

Choosing as test functions
as well as χ Q = ψ Q , and rearranging terms slightly yields Taking Eq. (5.14) into account, terms given the same color cancel out when we sum up these equations to obtain Hence, we end up with an energy given by dt .

(5.22)
Here, we see the first consequence of the presence of a Laplace operator in Eq. (5.7). This is actually similar to the dependence on U as ⇒ V is the time-integrated symmetrized Jacobi matrix of U and W its time-integrated curl. However, we want to remark that this energy is a mathematical construct and should not be interpreted as a quantity with a physical meaning.

Remark 2.
We want to mention that this special choice of testing (see (5.15)-(5.19)) is required for the construction of the energy (5.22) and is only possible due to our variable substitution T ← ∂ t Q, we introduced in (2.10).
(5.23) 3 be finite dimensional subspaces of the given Sobolev spaces. Q h is a finite element space of continuous piecewise linear function, Z Q,h a boundary element space of continuous piecewise linear functions, and P Q,h a boundary element space of piecewise constant functions. We denote the chosen bases of these spaces by

FEM-BEM spatial semi-discretization. To discretize in space
. Furthermore, Z Q,h and P Q,h contain the traces of Q h , i.e., γ D Q h ⊆ Z Q,h and γ N Q h ⊆ P Q,h . This is not problematic since we assume that the boundary of Ω is approximated in a piecewise polygonal manner. The inclusion also implies γ D U h ⊆ Z U ,h and γ N U h ⊆ P U ,h . The semi-discretized problem then reads: h , and χ U,h ∈ Z U,h , Eqs. (5.9)-(5.14) hold with the respective substitutions. We can turn the result into matrix-vector-form by testing with the according basis functions in the nodes of the discretization) 4 : Here, we have and for the boundary The terms in gray are meant to be perturbations. They do not have any physical meaning, but need to be considered in our stability analysis later on. For Eqs. (5.24) and (5.27), perturbations are considered as parts of the actual right-hand sides F and G, respectively. Further, we divided the boundary term from Eq. (5.12) into two terms and plugged in φ Q,h for one, but ψ Q,h for the other. Thus, terms involving these boundary contributions and the related terms from Eq. (5.28) also cancel out in the semi-discrete and discrete settings.
and, thus, represent the matrix B(s) by . (5.33)

Time discretization via leapfrog and convolution quadrature
In order to apply the leapfrog method, see, e.g., [12], and convolution quadrature, we use the same ideas as in [8], but have to account for the additional terms involving D 1 Q as well as the equation for Q and related boundary conditions. As a consequence, we get a time-stepping scheme which is implicit for Q, in accordance with [12], Section 1.8. Thus, we obtain where the bar above a quantity denotes an average, e.g., and analogously for all others. Further, we introduced where δ > 0 is a stabilization parameter which is discussed in more detail in Section 8 and the dot above a quantity denotes a central difference with respect to time, i.e., and similarly for other quantities. Please note that there is another discretization for the first-order derivative with respect to time in play here, which is denoted by ∂ Δt t . It is used for the discretization of the Calderón operator and given by the backward differentiation formula BDF-2 according to [4], Section 2.3.
Let us compare the above system of equations to the known ones for the acoustic and the elastic wave equations. For both of those, only the system of equations incorporating the Calderón operator is implicit.
Thus, it is possible to express the analog to U decoupling such that a more concrete expression can be given for the (still implicit) equations involving the Calderón operator, for the acoustic wave equation in [4], Section 5.3, and for the elastic wave equation in [8], Section 5.3. Due to the implicit nature of Eqs. (5.36) and (5.37), we cannot do something like this here. Only Eqs. (5.34), (5.35), (5.38), and (5.39) can be treated as explicit time-stepping schemes, whereas all other equations must be considered together as a coupled system of implicit equations. Further, we provide the equations for V and W using a staggered grid, cf. [12], Eq. (1.7), as this is the version used to prove stability results and is obtained by substituting (5.38) and (5.39) into (5.34) and (5.35), respectively, so that We will also need a set of semi-discrete equations which is of third order in time. These can be derived by taking derivatives with respect to time of Eqs. For our stability considerations, we also need a formulation which is closer to a time discretized version of that system which reads where the dot above a quantity is to be interpreted as a central difference in time as in Eq. (5.44) and a bar denotes averaging, see Eq. (5.41). Hence we havë The accompanying equations for the boundary densities are

Stability of the spatial semi-discretization
Let us assume that our bases of U h , V h , Q h , Z h , and P h are orthonormal in L 2 (Ω) 3 , (L 2 (Ω)) 3×3 , L 2 (Ω), (H 1 2 (Γ)) 4 , and (H − 1 2 (Γ)) 4 , respectively. Then the matrices M 0 and M 2 become unit matrices. The matrix M 1 , however, becomes a diagonal matrix whose diagonal elements are all 1 2 . Consequently, M −1 1 is a diagonal matrix with a 2 for each entry on the diagonal.
To analyze the propagation of errors in the spatial discretization, we consider bounds of the Euclidean norms of the solution of these equations, where we take the norm of the perturbations F , G, I, K, ν, ς, η, and ξ into account. Similar to the continuous field energy we introduce a semi-discrete field energy and show its boundedness.

Lemma 7. The semi-discrete field energy
is bounded along the solution by Proof. This lemma is proven in a similar way as the corresponding results for the acoustic wave equation [4] or the elastodynamic wave equation [8] by applying the bounds of the boundary integral operators.
By combining the previous results, we obtain the following theorem.
Theorem 11. Assume that the initial values U (·, 0), ⇒ V (·, 0), W (·, 0), and Q(·, 0) have their support in Ω. Let the initial values for the semi-discretization be chosen as Here P − denotes the L 2 (Ω)-orthogonal projection onto the corresponding finite element spaces. If we assume that the solution of the wave equation is sufficiently smooth, then the error of the FEM-BEM semi-discretization is bounded by

4)
where the constant m(t) grows at most polynomially with t.
Proof. For the first term on the left-hand side, the triangle inequality yields and similar estimates hold for all other quantities. The second term on the right can be estimated according to Lemma 9. For the sum of the first terms, observe that the differences U h (t) − P U h U (t) , and for the boundary densities All that is left to do is apply the stability results from the previous section together with Lemmas 9 and 10 to get the desired result, similar to the proof of Theorem 7.1 in [4] and Theorem 6.1 in [15].

Stability of the full discretization
We want to study the stability of the full discrete scheme under the corresponding CFL condition. Inspecting Eqs. (5.36), (5.37), (5.40), (5.45), and (5.46), we observe that Eq. (5.36) is similar in spirit to the Crank-Nicolson scheme for parabolic PDEs. Thus, there is no CFL-or other stability condition with regard to Q and we only need the CFL condition from elasticity, which reads As a further stability condition, it is sufficient to assume resulting from the considerations of the proof of Lemma 12. This guarantees that the discrete energy E n ≥ 0 for all n = 0, 1, 2, . . ..  by Taylor expansion that then where the first term represents the consistency error and the second one is of order O(Δt 2 ). For the further defects of the full discrete equations, we refer to the equations of the spatial semi-discretization. In addition, we obtain for a temporally smooth solution terms of order O(Δt 2 ). Further, the initial values of the time discretization are chosen as described in [18], so that H − 1 2 is well defined and is of order O(Δt 2 ) as well. Summarizing the discrete stability results and the error bounds of the spatial semi-discretization, we finally obtain the error bound for the full discretization.

Summary and conclusion
In this paper, we presented a numerically stable FEM-BEM coupling for the thermoelastic wave equation with transparent boundary conditions. We started with the introduction of the involved equations and stated the corresponding fundamental solution as well as the required layer potentials. Based on these, we took a look at the boundary integral operators which are needed for the construction of the Calderón operator. The key issue for the stability analysis was the coercivity result of the Calderón operator together with energy techniques. Here, we first considered the semi-discretization followed by the full discretization including time, where we applied a leapfrog scheme in the domain's interior and convolution quadrature on the boundary. Finally, we ended up with an asymptotically optimal convergence result.
In order to manage the coupled system of a vector-valued PDE (second order in time and space) and a third-order scalar-valued PDE (first order in time and second order in space), we introduced a timeintegrated temperature deviation. In particular, this allowed us to handle the coupling term without having to use a factor s in the test functions for the displacement vector, see, e.g., [24], Theorem 3.9. All in all, we understand this contribution as a step-stone to generalize the methods from [4,8] to more general sets of partial differential equations, as, e.g., those which describe wave propagation in poroelastic media [24].

Funding Information Open Access funding enabled and organized by Projekt DEAL.
Open Access. This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.