A priori error estimates for the finite element approximation of Westervelt's quasilinear acoustic wave equation

TL;DR

This paper provides a priori error estimates for finite element approximations of Westervelt's quasilinear acoustic wave equation, demonstrating optimal convergence rates through theoretical analysis and numerical experiments.

Contribution

It introduces a novel a priori error analysis for finite element discretization of Westervelt's equation using Banach fixed-point theorem and inverse estimates.

Findings

Optimal convergence rates in $L^2$ norms for small data and mesh size.

Numerical experiments confirm theoretical error estimates.

Method avoids degeneracy issues in semi-discrete equations.

Abstract

We study the spatial discretization of Westervelt's quasilinear strongly damped wave equation by piecewise linear finite elements. Our approach employs the Banach fixed-point theorem combined with a priori analysis of a linear wave model with variable coefficients. Degeneracy of the semi-discrete Westervelt equation is avoided by relying on the inverse estimates for finite element functions and the stability and approximation properties of the interpolation operator. In this way, we obtain optimal convergence rates in $L^2$-based spatial norms for sufficiently small data and mesh size and an appropriate choice of initial approximations. Numerical experiments in a setting of a 1D channel as well as for a focused-ultrasound problem illustrate our theoretical findings.

Figures

Equations

$$(1-2ku)u_{tt}-c^2\Delta u-b\Delta u_t=2k{u}_t^2,$$

$$u_{tt}-\sum _{i=1}^n\frac{\partial }{\partial x_i}{A}_i(x,\nabla u)=f(x,t,u,\nabla u),$$

$$u_{tt}-\Delta b(u)+\frac{\partial }{\partial t}(a(u))=f(x,t),$$

$$\begin{multlined}\|u-u_{h}\|_{L^{\infty}(0,T;L^{2}(\Omega))}+\left\|u_{t}-u_{h,t}\right\|_{L^{\infty}(0,T;L^{2}(\Omega))}+\left\|u_{tt}-u_{h,tt}\right\|_{L^{2}(0,T;L^{2}(\Omega))}\\[5.69054pt] +h\|\nabla(u-u_{h})\|_{L^{\infty}(0,T;L^{2}(\Omega))}+h\left\|\nabla u_{t}-u_{h,t}\right\|_{L^{\infty}(0,T;L^{2}(\Omega))}\leq\,Ch^{s}.\end{multlined}\|u-u_{h}\|_{L^{\infty}(0,T;L^{2}(\Omega))}+\left\|u_{t}-u_{h,t}\right\|_{L^{\infty}(0,T;L^{2}(\Omega))}+\left\|u_{tt}-u_{h,tt}\right\|_{L^{2}(0,T;L^{2}(\Omega))}\\[5.69054pt] +h\|\nabla(u-u_{h})\|_{L^{\infty}(0,T;L^{2}(\Omega))}+h\left\|\nabla u_{t}-u_{h,t}\right\|_{L^{\infty}(0,T;L^{2}(\Omega))}\leq\,Ch^{s}.$$

$$|fgh{|}_{L^1}\le$$

$$xy\le \epsilon x^2+\frac{1}{4\epsilon }{y}^2,\text{where}x,y>0,\epsilon >0;$$

$$u(t)+v(t)\le C_1+C_2{\int }_0^{t}u(s)\text{d}s\text{ for all }t\in [0,T].$$

$$u(t)+v(t)\le C_1{e}^{C_{2}T}\text{ for all }t\in [0,T];$$

$$f\in H_0^1(\Omega )$$

$$S_h=\{u_h\in H_0^1(\Omega ):u_h{|}_K\in P_1(K),\forall K\in {\mathcal{T}}_h\}.$$

$$c_{1}h\le h_K\le c_2{\varrho }_K,K\in {\mathcal{T}}_h,$$

$$|u-\mathcal{U}{|}_{L^2}\le$$

$$|\nabla (u-\mathcal{U}){|}_{L^2}\le$$

$$|\chi {|}_{L^{\infty }}\le C_{\text{inv}}{h}^{-d/p}|\chi {|}_{L^p},1\le p<\infty ,$$

$$|v-I_{h}v{|}_{L^2}\le$$

$$|I_{h}v{|}_{L^{\infty }}\le$$

$$\{\begin{array}{l}{u}_{tt}-c^2\Delta u-b\Delta u_t=2k\left(u{u}_{tt}+u_t^2\right)\text{in }\Omega \times (0,T),\text{vskip}3.0ptplus1.0ptminus1.0pt\\ u=0\text{ on }\partial \Omega \times (0,T),\text{vskip}3.0ptplus1.0ptminus1.0pt\\ \left(u,u_t\right)=(u_0,u_1)\text{ on }\Omega \times \{t=0\}.\end{array}$$

$$\{\begin{array}{l}\left((1-2ku)u_{tt},\varphi \right)+c^2(\nabla u,\nabla \varphi )+b\left(\nabla u_t,\nabla \varphi \right)-2k\left(u_t^2,\varphi \right)=0,\text{vskip}3.0ptplus1.0ptminus1.0pt\\ \text{for all }\varphi \in H_0^1(\Omega )\text{ a.e. in time},\text{vskip}3.0ptplus1.0ptminus1.0pt\\ \left(u(0),u_t(0)\right)=(u_0,u_1).\end{array}$$

$$|u_0{|}_{H^2}^2+|u_1{|}_{H^2}^2\le {\varrho }_T,$$

$$u\in \mathcal{B}=\{u\in L^{\infty }(\Omega \times (0,T)):$$

$$\Vert u_t{\Vert }_{C{H}^1}\le M,(u(0),u_t(0))=(u_0,u_1)\},$$

$$\{\begin{array}{l}\alpha (x,t)u_{tt}-c^2\Delta u-b\Delta u_t+\beta (x,t)u_t=f(x,t)\text{in }\Omega \times (0,T),\text{vskip}3.0ptplus1.0ptminus1.0pt\\ u=0\text{ on }\partial \Omega \times (0,T),\text{vskip}3.0ptplus1.0ptminus1.0pt\\ \left(u,u_t\right)=(u_0,u_1)\text{ on }\Omega \times \{t=0\},\end{array}$$

$$u_h(x,t)=\sum _{i=1}^{N_h}{\xi }_i(t)w_i(x),$$

$$\left({\alpha }_h{u}_{h,tt},\varphi \right)+c^2(\nabla u_h,\nabla \varphi )+b\left(\nabla u_{h,t},\nabla \varphi \right)+\left({\beta }_h{u}_{h,t},\varphi \right)=(f_h,\varphi ),$$

$$\left(u_h(0),u_{h,t}(0)\right)=(u_{h,0},u_{h,1}),$$

$${\Vert u_{h,tt}\Vert }_{L^2{L}^2}^2+\Vert \nabla u_h{\Vert }_{L^{\infty }{L}^2}^2+{\Vert \nabla u_{h,t}\Vert }_{L^{\infty }{L}^2}^2+{\Vert \nabla u_{h,t}\Vert }_{L^2{L}^2}^2$$

$$\le$$

$$C({\alpha }_h,{\beta }_h,T)=C_1\text{exp}(C_2(\Vert {\alpha }_h{\Vert }_{L^{\infty }{L}^3}^2+\Vert {\beta }_h{\Vert }_{L^{\infty }{L}^3}^2+1)T).$$

$$\{\begin{array}{l}{M}_h(t){\xi }_{h,tt}+K_h{\xi }_h+C_h(t){\xi }_{h,t}=F_h,\\ {\xi }_h(0)={\xi }_{h,0},\\ {\xi }_{h,t}(0)={\xi }_{h,1},\end{array}$$

$$M_h(t)=$$

$$K_h=$$

$$C_h(t)=$$

$$|({\alpha }_h{w}_i,w_j)|\le |{\alpha }_h{|}_{L^2}|w_i{|}_{L^4}|w_j{|}_{L^4},$$

Full text

{AMS}

35L05, 65M15, 65M60

†††Technical University of Munich, Department of Mathematics, Chair of Numerical Mathematics, Boltzmannstraße 3, 85748 Garching, Germany (, ).

1 Introduction

The goal of the present work is to analyze a spatial discretization by piecewise linear finite elements in nonlinear acoustics. To this end, we study a discretization of Westervelt’s wave equation for the acoustic pressure $u$

[TABLE]

which represents a classical model for nonlinear ultrasound propagation through thermoviscous fluids [53]. Our research is motivated by a rising number of nonlinear ultrasound applications in medicine and industry [4, 15, 37, 39, 41]. In (1), the constant $c$ denotes the speed of sound, $b$ is the sound diffusivity, and $k=\beta_{a}/(\varrho c^{2})$, where $\varrho$ is the mass density and $\beta_{a}$ the coefficient of nonlinearity of the medium.

Westervelt’s equation is a strongly damped quasilinear wave equation with potential degeneracy due to the factor $1-2ku$ next to the second time derivative. For its derivation and the theoretical foundations of nonlinear acoustics, we refer to [9, 13, 20, 53], while results on the existence of smooth solutions of (1) can be found in [25, 26, 34]. Efficient simulation of the Westervelt equation and, in general, nonlinear sound propagation by the finite element method has been an active area of research. We refer to, e.g., [16, 22, 24, 29, 35, 36, 50, 52], which all focus on algorithmic aspects of finite element discretizations without any a priori analysis.

Error analysis for the standard finite element discretization of linear wave equations is an extensively studied topic; see, e.g., [1, 2, 3, 11, 18, 32, 49] and the references given therein. In particular, we single out the work on a priori analysis in [1] which provides $L^{\infty}(0,T;L^{2})$ error estimates for the undamped linear wave equation and the results on error bounds for strongly damped linear wave equations [32, 49]. Results on a class of nonlinear wave equations of a divergent type are also well-established. In [10], error analysis is provided for a semi-discretization of nonlinear wave equations of the form

[TABLE]

with a monotonicity condition on the corresponding bilinear form; cf. [10, Theorem 3.2]. In [47], semi-discretization for the following damped model is considered

[TABLE]

with $b^{\prime}(u)\geq M_{0}>0$. In [12], convergence of a full discretization for a class of nonlinear second-order in time evolution equations is provided where the operator acting on the first time derivative is assumed to be hemicontinuous, monotone, coercive, and to fulfill a certain growth condition. Moreover, the operator acting on the solution is assumed to be linear, bounded, symmetric, and strongly positive. We also mention the results in [33] on a class of problems of nonlinear elastodynamic and in [40] on the discontinuous Galerkin methods for a class of divergent-type nonlinear hyperbolic equations.

This work contributes to the finite element analysis of Westervelt’s equation in two ways. We first prove that, coupled with non-zero initial conditions and homogeneous Dirichlet data, its semi-discretization by piecewise linear finite elements has a unique solution which remains bounded in an appropriately chosen norm. Secondly, we derive an optimal a priori error estimate that has the form

[TABLE]

where $\max\{1,d/2\}<s\leq 2$. Our results are intended to enhance the numerical analysis of strongly damped quasilinear wave equations where the nonlinearities in the equation involve the time derivatives of the solution. We note that a particular feature of the present quasilinear equation is that the non-degeneracy is not a priori given. In our proofs we have to ensure that the factor $1-2ku_{h}$ next to the second time derivative remains positive.

In the continuous analysis of the Westervelt equation, non-degeneracy is typically achieved by a higher-regularity result for the solution and the use of an embedding, e.g., $H^{2}(\Omega)\hookrightarrow L^{\infty}(\Omega)$; see [25, Theorem 3.1]. Such a strategy is not possible here since we use piecewise linear basis functions. Instead we employ inverse estimates for finite element functions and the stability and approximation properties of the Scott–Zhang interpolation operator [44].

Our analysis relies on the Banach fixed-point theorem combined with error estimates for a linear wave equation with variable coefficients. Therefore, in this work, we also obtain error estimates for strongly damped variable coefficient wave equations that take coefficient error into account as relevant, e.g., in optimal control problems in nonlinear acoustics [7, 28, 36].

The rest of the paper is organized as follows. Section 2 introduces the notation and lays out the most important theoretical results in Sobolev and finite element spaces that we often use in the analysis. In Section 3, we discuss the continuous problem and its well-posedness. In Section 4, we then study a linearized Westervelt equation with variable coefficients and prove that its semi-discretization has a unique solution. Section 5 focuses on the a priori analysis of this linear model. In Section 6, we show well-posedness and derive convergence rates for the semi-discrete Westervelt equation. Finally, Section 7 contains numerical examples that illustrate our theory.

2 Theoretical preliminaries

We begin by setting the notation and summarizing some auxiliary properties of Sobolev and finite element spaces that we will frequently use in the analysis.

2.1 Notation

We denote the standard $L^{2}$ inner product by $(\cdot,\cdot)$. The norms in Sobolev spaces $L^{p}(\Omega)$ and $W^{q,p}(\Omega)$ are denoted by $|\cdot|_{L^{p}}$ and $|\cdot|_{W^{q,p}}$, respectively, where $1\leq p\leq\infty$, $1\leq q<\infty$. The norms in Bochner spaces $W^{q,p}(0,T;W^{r,s}(\Omega))$ are denoted by $\|\cdot\|_{W^{q,p}W^{r,s}}$, where $0\leq q,r<\infty$, $1\leq p,s\leq\infty$. We also introduce the spaces $\dot{H}^{s}(\Omega)=H_{0}^{1}(\Omega)\cap H^{s}(\Omega)$, for $1\leq s\leq 2$.

The constants $0<C_{i}<\infty$, $i\in\mathbb{N}$, appearing in the estimates denote generic constants that might depend on the coefficients in the equation and the domain $\Omega$, but not on the mesh size. Throughout the paper, we assume $T>0$ to be a fixed time horizon.

2.2 Auxiliary inequalities

Let $\Omega\subset\mathbb{R}^{d}$, where $d\in\{1,2,3\}$, be a bounded domain with Lipschitz regular boundary. The nonlinear terms appearing in the Westervelt equation are of a quadratic type, so after variational testing, we often have to employ Hölder’s inequality for a product of three functions. In particular, we frequently make use of the following three special cases of Hölder’s inequality:

[TABLE]

with $(p,q,r)\in\{(2,4,4),(3,6,2),(\infty,2,2)\}$. We also often employ a special case of Young’s $\varepsilon$-inequality in the form

[TABLE]

see [14, Appendix B]. Let $u$ and $v$ be non-negative continuous functions and $C_{1},C_{2}<\infty$ non-negative constants such that

[TABLE]

Then the following modification of Gronwall’s inequality holds

[TABLE]

see [17, Lemma 3.1]. Finally, we recall the Sobolev embeddings

[TABLE]

for $1\leq p\leq 6$, where $C_{H_{0}^{1},L^{p}}<\infty$, noting that $d\leq 3$.

2.3 Finite element spaces

We consider the discretization in space by continuous piecewise linear finite elements that vanish on the boundary. Let $\Omega\subset\mathbb{R}^{d}$, $d\in\{2,3\}$, be a convex polygonal domain. For $h\in(0,\overline{h}]$, let $\mathcal{T}_{h}$ be a triangulation of $\Omega$ made of triangles (in $\mathbb{R}^{2}$) or of tetrahedrons (in $\mathbb{R}^{3}$) so that $\Omega=\cup_{K\in\mathcal{T}_{h}}K$. We denote by $P_{1}(K)$ the space of polynomials on $K$ of degree no greater than $1$. We introduce the finite element space as

[TABLE]

We assume that $\{\mathcal{T}_{h}\}_{0<h\leq\overline{h}}$ is a quasiuniform family: there are constants $0<c_{1},c_{2}<\infty$ such that

[TABLE]

where $h_{K}$ denotes the diameter of the triangle (tetrahedron) K, $\varrho_{K}$ stands for the diameter of the greatest ball (sphere) included in K, and $h=\max_{K\in\mathcal{T}_{h}}h_{K}.$

It is known that there exists $\mathcal{U}\in S_{h}$ and $0<C<\infty$ such that

[TABLE]

for $u\in\dot{H}^{s}(\Omega)$, $1\leq s\leq 2$; see [19].

Inverse estimates. Under the assumptions made above on the family $\{S_{h}\}_{0<h\leq\overline{h}}$, there is a $0<C_{\textup{inv}}<\infty$ such that

[TABLE]

for every $\chi\in S_{h}$; see [5, Theorem 4.5.11]. We will need the special cases $p=2$ and $p=4$ in the proofs.

Bounds for the interpolation error. In our analysis, we will employ an interpolant $I_{h}:W^{l,p}(\Omega)\rightarrow S_{h}$ of Scott–Zhang type, where $0\leq l\leq 1$, $1\leq p\leq\infty$; cf. [5, 44]. The following approximation and stability properties hold:

[TABLE]

where $0<C_{\textup{app}},C_{\textup{sta}}<\infty$; see [5, Theorem 4.8.12 and Corollary 4.8.15].

3 The continuous problem

We start from the following initial-boundary value problem for the Westervelt equation

[TABLE]

The weak form of the problem is then given by

[TABLE]

This problem is known to be well-posed for small data.

Theorem 3.1**.**

[25, Theorem 3.1]* Let $T>0$, $b,k,c^{2}>0$, $0<m<\frac{1}{4k}$, and $M>0$ be arbitrary. Assume that*

[TABLE]

with $\varrho_{T}$ sufficiently small. Then there exists a unique solution $u$ of (10) such that

[TABLE]

*and such that $\Delta u,\ u_{tt},\ \nabla u_{t}\in L^{\infty}(0,T;L^{2}(\Omega)),\ \nabla u_{tt}\in L^{2}(0,T;L^{2}(\Omega)).$ *

We also refer to [34] where the results of Theorem 3.1 are generalized by employing the maximal $L^{p}$ regularity approach. Results on the existence of very smooth solutions for a reformulation of the problem in terms of the acoustic velocity potential $\psi$, where $u=\varrho\psi_{t}$, can be found in [27, 30].

Note that the well-posedness holds for sufficiently small data which by continuity implies smallness of $u$ in the appropriate norms. The condition $\|u\|_{L^{\infty}(\Omega\times(0,T))}\leq m<1/(4k)$ in (11) ensures that the equation does not degenerate. For the well-posedness of the semi-discrete problem, we will also need smallness of data and a bound on the approximate solution that guarantees non-degeneracy. It is also worth noting that the strong damping (i.e., $b>0$) is needed for the continuous problem to be well-posed and the same will hold for the semi-discrete equation.

Going forward, we assume that (9) has a unique solution. We will impose additional conditions on the regularity of $u$ when needed for the convergence results.

4 Finite element approximation of the linearized Westervelt equation with variable coefficients

We first provide numerical analysis of an initial-boundary value problem for a linear wave equation with variable coefficients which can be interpreted as a linearization of the Westervelt equation. We study the following initial boundary value problem for a non-degenerate equation:

[TABLE]

where $0<\alpha_{0}\leq\alpha(x,t)\leq\alpha_{1}$ a.e. in $\Omega\times(0,T)$. Analysis of the linearization (12) allows to later define an iterative map on which we will apply the Banach fixed-point theorem. However, finite element approximation of the partial differential equation in (12) is also of independent interest. For example, this model with $b=f=0$ appears in [6] and is motivated by the study of the transonic gas dynamics. The adjoint problems for the Westervelt equation which arise in the optimal control and shape optimization works [7, 28, 36] have (after time reversal) the form of this PDE as well.

We refer to [27, Proposition 7.2] for the sufficient conditions under which problem (12) has a unique solution $(u,u_{t})$ in $C([0,T];H^{j}(\Omega)\times H^{j-1}(\Omega))$, where $j\in\{2,4\}$. We therefore proceed with the assumption that problem (12) has a unique solution. The conditions on the regularity of $u$ are specified when needed for the a priori estimates. It is implicitly assumed that the coefficients $\alpha$ and $\beta$, the initial data $(u_{0},u_{1})$, and the source term $f$ are sufficiently smooth for such a regularity to hold.

Results on the error estimates for special cases of (12) with constant coefficients are available in the literature. Analysis of the Galerkin approximation of (12) for the case $\alpha=1$, $b=0$, and $\beta=0$ is performed in [1]. The case of a strongly damped wave equation (i.e. with a fixed positive constant $b$) and with $\alpha=1$, $\beta=f=0$ is analyzed in [32, 46, 49].

Let $\{S_{h}\}_{0<h\leq\overline{h}}$ be a family of subspaces of $H_{0}^{1}(\Omega)$ defined in (5) with basis $\{w_{i}\}_{i=1}^{N_{h}}$. We consider Galerkin approximations in space

[TABLE]

where $\xi_{i}:(0,T)\to\mathbb{R}$ are coefficient functions for $i\in[1,N_{h}]$. Let $\alpha_{h}$, $\beta_{h}$, and $f_{h}$ be approximations of functions $\alpha$, $\beta$, and $f$, respectively, in $S_{h}$. {assumption} We assume that the approximate coefficients and the source term satisfy the following conditions

• •

$\alpha_{h}\in L^{\infty}(0,T;L^{\infty}(\Omega))$,  $\exists\ \alpha_{0}:\,\alpha_{h}\geq\alpha_{0}>0$ a.e. in $\Omega\times(0,T)$,

• •

$\beta_{h}\in L^{\infty}(0,T;L^{3}(\Omega))$,

• •

$f_{h}\in L^{2}(0,T;L^{2}(\Omega))$.

For a given $h\in(0,\overline{h}]$, we next study a semi-discretization of (12) in $S_{h}$ and prove that it has a unique solution.

Theorem 4.1**.**

Let $c^{2}$, $b>0$ and let Assumption 4 hold. For each $h\in(0,\overline{h}]$, there exists a unique function $u_{h}\in H^{2}(0,T;S_{h})$ which satisfies

[TABLE]

for all $\phi\in S_{h}$, a.e. in time, and

[TABLE]

where $u_{h,0}$ and $u_{h,1}$ are approximations of $u_{0}$ and $u_{1}$ in $S_{h}$. Moreover, the following a priori bound holds

[TABLE]

The constant above is given by

[TABLE]

Proof 4.2**.**

*The proof follows a general framework of the well-posedness proofs for the linearizations of the classical nonlinear acoustic equations that are based on the Galerkin approximations in space. In particular, we refer to [28, Theorem 1] and [16, Proposition 1]. However, for the continuous problem, the basis functions have to be in $H^{2}(\Omega)$ to later guarantee the non-degeneracy of the nonlinear model via the embedding $H^{2}(\Omega)\hookrightarrow L^{\infty}(\Omega)$. Our basis functions are only $H^{1}$ regular which changes the a priori estimates that we will derive.

Step 1: Existence of a solution. We denote by $\xi_{h,0}=[\xi_{1,0}\ldots\xi_{N_{h},0}]^{T}$ and $\xi_{h,1}=[\xi_{1,1}\ldots\xi_{N_{h},1}]^{T}$ the components of the given initial approximations $u_{h,0}$ and $u_{h,1}$, respectively. Then our semi-discrete problem is to find $\xi_{h}=[\xi_{1}\ldots\xi_{N_{h}}]^{T}$ such that*

[TABLE]

where the matrices are given by

[TABLE]

and the source term is given by $F_{h}=[F_{1}\ldots F_{N_{h}}]^{T}$, $F_{j}=(f_{h},w_{j})$, with $1\leq i,j\leq N_{h}$. Note that the matrices and the right-hand side vector are all well-defined since

[TABLE]

a.e. in time. Furthermore, the matrix $M_{h}$(t) is invertible for a.e. $t\in[0,T]$; cf. [28, Theorem 1]. The statement follows from the fact that $M_{h}(t)$ is positive definite. Indeed, for any $z\in\mathbb{R}^{N_{h}}\setminus\{0\}$, we have

[TABLE]

for a.e. $t\in[0,T]$. Thanks to the fact that $M_{h}$ is invertible, the matrix equation in (17) can be rewritten as

[TABLE]

*Now the existence of a solution $u_{h}\in H^{2}(0,T_{h};S_{h})$ follows from the standard ODE theory; see, for example, [43, Chapter 1]. To extend the existence interval to $[0,T]$, we next show that $u_{h}$ remains bounded on $[0,T]$ in the appropriate norms.

Step 2: A priori estimate. We want to derive an priori bound for $u_{h}$. To this end, we test our problem with two different test functions. We first test (13) with $\phi=\lambda u_{h,t}\in S_{h}$, where $\lambda>0$, and integrate with respect to time from [math] to $t$, $t\leq T_{h}$. After some standard manipulations, this action results in*

[TABLE]

where $\varepsilon>0$ and $\lambda>0$ will be conveniently chosen. To be able to bound the term $\displaystyle\left\|u_{h,tt}\right\|^{2}_{L^{2}L^{2}}$ that appears on the right-hand side above, we next test (13) with $\phi=u_{h,tt}\in S_{h}$. After integrating over $(0,t)$, this action yields the second inequality

[TABLE]

Above, we have estimated the $c^{2}$ term by first integrating by parts with respect to time and then employing Hölder’s inequality and Young’s $\varepsilon$-inequality with $\varepsilon\in\{b/4,1/2\}$:

[TABLE]

To absorb $\displaystyle\frac{c^{4}}{b}|\nabla u_{h}(t)|^{2}_{L^{2}}$ by the corresponding term on the left side in (18), we need to choose $\lambda>0$ sufficiently large so that $\lambda c^{2}/2>c^{4}/b$. By adding the derived inequalities (18) and (19), we then obtain

[TABLE]

We choose $\lambda=4c^{2}/b$ and $\varepsilon=\alpha_{0}/6$, apply Gronwall’s inequality to (21), and take the essential supremum over $t\in(0,T_{h})$. In this way, we obtain

[TABLE]

*The right-hand side of (22) does not depend on $T_{h}$, so we can show by an argument of contradiction that we are allowed to extend the existence interval of $u_{h}$ to $[0,T]$; i.e., $T_{h}=T$ and estimate (15) holds. *

5 A priori estimates for the linearized Westervelt equation with variable coefficients

We now focus on proving a priori estimates for the linearized Westervelt equation that also take into account approximation error of the coefficients and the source term. We wish to estimate $u-u_{h}$. We follow the usual approach in the finite element analysis and split this difference into

[TABLE]

where $R_{h}$ denotes the elliptic projection; cf. [5, 48]. The idea is to rely on the existing results on elliptic projectors to bound $\varrho=u-R_{h}u$, whereas $\theta=R_{h}u-u_{h}$ will be seen as a solution of a wave PDE with a source term. By deriving estimates for this PDE, we will find a bound for $\theta$.

5.1 Auxiliary results for the elliptic projection

We first recall two useful results for an auxiliary elliptic problem for $\varrho$. We employ the Ritz projection $R_{h}$, i.e., the orthogonal projection with respect to the product $(\nabla u,\nabla\phi)$.

Lemma 5.1**.**

[1, Lemma 2.1]* Let $u$ be the solution of (12). Then there exists a unique mapping $R_{h}u\in L^{2}(0,T;S_{h})$ which satisfies*

[TABLE]

Let $1\leq p\leq\infty$. If for some integer $k\geq 0$, $\frac{\partial^{k}u}{\partial t^{k}}\in L^{p}(0,T;H^{s}(\Omega)),$ then $\frac{\partial^{k}R_{h}u}{\partial t^{k}}\in L^{p}(0,T;S_{h})$, and

[TABLE]

*for some constant $C>0$ independent of $h$ and $u$, and $1\leq s\leq 2$. *

We also need a bound on the gradient of $\varrho$ to be able to later derive $H^{1}$ bounds for $u-u_{h}$.

Lemma 5.2**.**

Let $u$ be the solution of (12) such that $u\in L^{p}(0,T;H^{s}(\Omega))$, where $1\leq s\leq 2$, $1\leq p\leq\infty$. Then it holds

[TABLE]

Proof 5.3**.**

*The estimate follows directly from the Galerkin orthogonality of $u-R_{h}u$ a.e. in time, Céa’s lemma and the approximation property (6). *

We note that analogous bounds can be obtained for $\varrho_{t}$ and $\varrho_{tt}$ by differentiating (23) once and twice with respect to time.

5.2 Bounds for $\theta=R_{h}u-u_{h}$

Since we are able to estimate $\varrho$, we now focus on deriving two a priori bounds for $\theta$.

At this point, we choose the approximate initial data as Ritz projections of $u_{0}$, $u_{1}$ in order to have $\theta(0)=\theta_{t}(0)=0$.

Proposition 5.4**.**

Let $c^{2}$, $b>0$. Let $u$ be the solution of (12) which satisfies

[TABLE]

where $1\leq s\leq 2$. Let Assumption 4 hold and $\alpha_{h,t}\in L^{\infty}(0,T;L^{3}(\Omega))$. Furthermore, let $(u_{h,0},u_{h,1})=(R_{h}u_{0},R_{h}u_{1})$. Then there exists a positive constant $C=C(\alpha_{h},\beta_{h},T)$ such that

[TABLE]

Proof 5.5**.**

The main idea of the proof is to see $\theta$ as a solution of a wave equation with variable coefficients and a source term and then test that equation with a suitable test function. Compared to the similar results for solutions of linear wave equations with constant coefficients **[1, 32*]**, we also take into account the error of the varying coefficients. *

By subtracting the weak forms for $u$ and $u_{h}$ and recalling the definition of $R_{h}u$, we find that $\theta$ solves

[TABLE]

for all $\phi\in S_{h}$ a.e. in time. We next want to test (26) with $\displaystyle\phi=\theta_{t}$, noting that $\theta_{t}\in S_{h}$ a.e. in time. To get optimal error estimates for $u_{h}$, it is important to only employ the $L^{2}$ spatial norm of $\varrho_{t}$ and $\varrho_{tt}$ in our estimates. We also have to pay special attention to estimating the last two terms on the right-hand side. Having in mind the nonlinear problem where we will know that

[TABLE]

we should only have the terms $\alpha-\alpha_{h}$ and $\beta-\beta_{h}$ in the $L^{2}$ spatial norm in our estimates to ensure optimal error rates for the nonlinear Westervelt equation. Testing (26) with $\theta_{t}$, integrating over $(0,t)$, and employing Hölder’s inequality then results in

[TABLE]

for a.e. $t\in[0,T]$. Above, we have used the identity

[TABLE]

and the fact that $\theta(0)=\theta_{t}(0)=0$. By further employing the embedding results (4) and Young’s $\varepsilon$-inequality (2) with $\varepsilon\in\{b/8,1\}$ to handle the product terms, we get

[TABLE]

Applying Gronwall’s inequality (3) to the above estimate and taking the essential supremum over $(0,T)$ then leads to

[TABLE]

with $\overline{C}(\alpha_{h},\beta_{h},T)=C_{3}\,\textup{exp}\left(C_{4}\left(\left\|\alpha_{h,t}\right\|^{2}_{L^{\infty}L^{3}}+\|\beta_{h}\|^{2}_{L^{\infty}L^{3}}+1\right)T\right)\,$. Thanks to the results on the elliptic projector stated in Lemma 5.1 which provide the upper bounds on $\|\varrho_{t}\|_{L^{2}L^{2}}$ and $\|\varrho_{tt}\|_{L^{2}L^{2}}$, we then obtain the final bound (25), where the constant is given by

[TABLE]

Remark 5.6**.**

We note that if $\alpha_{h}\in C([0,T];L^{\infty}(\Omega))$, the same error order can be obtained if we choose any $u_{1,h}$ such that $|u_{1}-u_{1,h}|_{L^{2}}\leq Ch^{s}$. This is due to the fact that we would then have an additional term on the right-hand side of (27) that is of order $h^{s}$:

[TABLE]

To be able to later employ a fixed-point approach and derive an a priori estimate for the nonlinear model, we also need to bound $\|\theta_{tt}\|_{L^{2}L^{2}}$.

Proposition 5.7**.**

Let $c^{2}$, $b>0$ and let $u$ to be the solution of (12) that satisfies

[TABLE]

where $s\in[1,2]$. Let Assumption 4 hold and let $\beta_{h}\in L^{2}(0,T;L^{\infty}(\Omega))$ and $\alpha_{h,t}\in L^{\infty}(0,T;L^{3}(\Omega))$. Assume that $(u_{h,0},u_{h,1})=(R_{h}u_{0},R_{h}u_{1})$. Then there exists a positive constant $C=C(\alpha_{h},\beta_{h},T)$ such that

[TABLE]

Proof 5.8**.**

To obtain the higher order estimate, we additionally test (26) with $\displaystyle\phi=\theta_{tt}$. After integrating over $(0,t)$ and recalling that $\theta(0)=\theta_{t}(0)=0$, we find

[TABLE]

for $t\in[0,T]$. Above we have estimated $\displaystyle c^{2}\int_{0}^{t}\int_{\Omega}\nabla\theta\cdot\nabla\theta_{tt}\,\textup{d}x$ in the same manner as (20). Since $\theta(0)=0$, we can further infer that

[TABLE]

Note that, compared to Proposition 5.4, we had to introduce additional assumptions (28) on the $L^{\infty}$ regularity of $u_{t}$ and $u_{tt}$. This is again due to the fact that we do not want to have higher than $L^{2}$ spatial norms of $\alpha-\alpha_{h}$ and $\beta-\beta_{h}$ on the right-hand side of (30).

We choose $\varepsilon<\alpha_{0}/6$ and add (30) to estimate (27) multiplied by $\lambda>0$ such that $\lambda c^{2}/2>c^{4}/b$. We then apply Gronwall’s inequality to the resulting estimate to obtain

[TABLE]

where we have also used that $|v|_{L^{3}}\leq C_{L^{\infty},L^{3}}|v|_{L^{\infty}}$, for $v\in L^{\infty}(\Omega)$. Employing the bounds on $\|\varrho_{t}\|_{L^{\infty}L^{2}}$ and $\|\varrho_{tt}\|_{L^{2}L^{2}}$ then leads to the estimate (29), where the constant is given by

[TABLE]

5.3 A priori estimate for the linear equation

We can now state the a priori estimate for the linearized Westervelt equation with variable coefficients.

Theorem 5.9**.**

Let the assumptions of Proposition 5.7 hold. Then the following a priori estimate is satisfied

[TABLE]

where $u_{h}$ solves (13), (14). The constant appearing above is given by

[TABLE]

Proof 5.10**.**

The estimate follows directly by splitting the difference $u-u_{h}$ into the $\theta$ and $\varrho$ terms, and then employing Proposition 5.7, Lemma 5.1, and Lemma 5.2, and the fact that

[TABLE]

*for some constant $C>0$ independent of $h$ and $u$. *

We note that Theorem 5.9 also includes, as a special case where $\alpha=1$, $\beta=f=0$, the a priori estimate for strongly damped linear wave equations with constant coefficients; this result corresponds to [32, Theorem 3.4]. If we do not have to take the coefficient error into account, regularity conditions (28) for $u$ can be relaxed to

[TABLE]

since we lose the last two terms in the estimate (31).

6 Finite element approximation of Westervelt’s equation

We are now ready to study discretization in space of the initial-boundary value problem (9) for the Westervelt equation. We want to prove that it has a unique solution in a neighbourhood of $u$. To this end, we rely on the Banach fixed-point theorem.

Theorem 6.1**.**

[A priori error estimate]* Let $c^{2}$, $b$, $k>0$, and $T>0$. Assume that the initial-boundary value problem (10) for the Westervelt equation has a unique solution which satisfies*

[TABLE]

where $\max\{1,d/2\}<s\leq 2$. Then for sufficiently small

[TABLE]

and $h$, there exists a unique $u_{h}\in H^{2}(0,T;S_{h})$ in a neighbourhood of $u$ which satisfies equation

[TABLE]

for all $\phi\in S_{h}$ a.e. in time, and $\left(u_{h}(0),u_{h,t}(0)\right)=\left(R_{h}u_{0},R_{h}u_{1}\right).$ Furthermore, there exists a positive constant $C$ that depends on $m$, $M$, and $T$, but not on $h$, such that

[TABLE]

Proof 6.2**.**

The main idea of the proof is to define an iterative map on which we will apply the Banach fixed-point theorem while relying on the results for the linearized problem. We first introduce the set

[TABLE]

the constant $L>0$ is independent of $h$ and $h\leq\overline{h}$. Note that the set $\mathcal{B}_{h}$ is non-empty since $R_{h}u\in\mathcal{B}_{h}$. For $v_{h}\in\mathcal{B}_{h}$, we then consider the following linearization of our semi-discrete problem

[TABLE]

*We introduce an iterative map $\mathcal{F}:v_{h}\mapsto u_{h}$, noting that $\mathcal{F}$ will be well-defined thanks to the uniqueness result from Theorem 4.1. Clearly a fixed point of $\mathcal{F}$ would solve (33), so we proceed with verifying the conditions of the Banach fixed-point theorem.

$(\mathcal{B}_{h},d)$** is a complete metric space.** $\mathcal{B}_{h}$ is closed in $\mathcal{X}=W^{1,\infty}(0,T;H_{0}^{1}(\Omega))\cap H^{2}(0,T;L^{2}(\Omega))$ with respect to the metric $d$ induced by*

[TABLE]

*from which completeness follows.

$\mathcal{F}$** is a self-mapping.** Take any $v_{h}\in\mathcal{B}_{h}$. We want to show that $u_{h}=\mathcal{F}v_{h}\in\mathcal{B}_{h}$. We set*

[TABLE]

and check that the conditions of Theorem 4.1 and Theorem 5.9 are satisfied. We first show that the non-degeneracy condition on $\alpha_{h}$ is fulfilled. Employing the identity $v_{h}=v_{h}-I_{h}u+I_{h}u$ and relying on the inverse estimate (7) yields

[TABLE]

Then by additionally using the identity $v_{h}-I_{h}u=v_{h}-u+u-I_{h}u$ and the stability and approximation properties (8) of the interpolant, we conclude that

[TABLE]

By choosing $\overline{h}$, $M$, and $m$ sufficiently small so that

[TABLE]

we can guarantee that $\|v_{h}\|_{L^{\infty}L^{\infty}}\leq m_{0}<1/(2k)$. In this way, the non-degeneracy condition in Assumption 4 is fulfilled since

[TABLE]

We next want to bound $\|\beta_{h}\|_{L^{2}L^{\infty}}$, $\|\beta_{h}\|_{L^{\infty}L^{3}}$, $\|\alpha_{h}\|_{L^{\infty}L^{3}}$, and $\|\alpha_{h,t}\|_{L^{\infty}L^{3}}$ uniformly with respect to $h$. Similarly to (36), we derive the following estimate

[TABLE]

from which we have $\left\|\beta_{h}\right\|_{L^{2}L^{\infty}}\leq 2k(\sqrt{T}C_{\textup{inv}}L\overline{h}^{s-d/2}+C_{\textup{inv}}C_{\textup{app}}\overline{h}^{s-d/2}M+M)$. Furthermore, it holds that

[TABLE]

and

[TABLE]

According to Theorem 4.1, there exists a unique solution $u_{h}\in X_{h}$ of (35). Thanks to the a priori bound for the linearized Westervelt equation stated in Theorem 5.9, we have

[TABLE]

where the constant appearing above is computed according to (32) and the derived uniform bounds on $\alpha_{h}$ and $\beta_{h}$:

[TABLE]

Therefore, from (38) and the fact that $v_{h}\in\mathcal{B}_{h}$ we infer that

[TABLE]

*for sufficiently small $m$, $M$, and $\overline{h}$. We can conclude that $u_{h}\in\mathcal{B}_{h}$ and, therefore, $\mathcal{F}(\mathcal{B}_{h})\subset\mathcal{B}_{h}$.

$\mathcal{F}$** is a contraction.** Let $v_{h}^{(1)},v_{h}^{(2)}\in\mathcal{B}_{h}$ and $u^{(1)}_{h}=\mathcal{F}v^{1}_{h}$, $u^{(2)}_{h}=\mathcal{F}v^{(2)}_{h}$. We want to show that*

[TABLE]

We note that the difference $\psi_{h}=u^{(1)}_{h}-u^{(2)}_{h}$ satisfies

[TABLE]

with zero initial conditions, for all $\phi\in S_{h}$. The equation can be seen as a special case of the PDE we considered in Theorem 4.1 if we choose

[TABLE]

and zero initial conditions. Owing to Theorem 4.1, we then have the a priori bound

[TABLE]

where the constant above is computed according to (16),

[TABLE]

We next show that $\|u^{(1)}_{h,t}\|_{L^{\infty}L^{4}}$ and $\|u^{(2)}_{h,tt}\|_{L^{2}L^{4}}$ in (39) can be made small so that $\mathcal{F}$ is a contraction. On account of the inverse properties (7) of $\{S_{h}\}_{0<h\leq\overline{h}}$ , we first find that

[TABLE]

We can then bound $\varrho_{tt}=u_{tt}-R_{h}u_{tt}$ by employing Lemma 5.1 to get

[TABLE]

Above, we have also made use of the fact that $\|R_{h}u_{tt}\|_{L^{2}L^{4}}\leq C\|u_{tt}\|_{L^{2}H^{s}}$ for some $C>0$ independent of $h$ and $u$. Furthermore, since $u_{h}^{(1)}\in\mathcal{B}_{h}$, we can bound the term $\|u^{(1)}_{h,t}\|_{L^{\infty}L^{4}}$ as follows

[TABLE]

*where $C>0$ is independent of $h$ and $u$. Altogether from (39), (40), and (41), for sufficiently small $M$ and $\overline{h}$, we can conclude that $\mathcal{F}$ is contractive with respect to the topology induced by $|||\cdot|||$. The statement now follows by applying the Banach fixed-point theorem. *

We note that due to the presence of the strong damping in the model, the assumed regularity of the solution $u$ is to be realistically expected. The higher the sound diffusivity $b$ is, the less pronounced nonlinear effects are in the propagation, such as the steepening of the wave. We refer the interested reader to, e.g., [16, Section 7.1.1] for a detailed discussion on how the $b$-damping influences the behavior of the nonlinear acoustic models.

Remark 6.3**.**

If $\Omega$ is a bounded interval in $\mathbb{R}$, we can rely on the embedding $H_{0}^{1}(\Omega)\hookrightarrow L^{\infty}(\Omega)$ to avoid degeneracy of the semi-discrete Westervelt equation and then we do not need to employ the inverse properties of $\{S_{h}\}_{0<h\leq\overline{h}}$. Indeed, the $L^{\infty}$ bounds (36) and (37) in the proof of Theorem 6.1 can be replaced by

[TABLE]

respectively, for $v_{h}\in\mathcal{B}_{h}$. Similarly, when proving contractivity, we can replace estimate (39) by

[TABLE]

*Since $u^{(1)}_{h}$, $u^{(2)}_{h}\in\mathcal{B}_{h}$, we can easily show that $\mathcal{F}$ is a contraction for sufficiently small $M$ and $\overline{h}$. Therefore, the a priori bound (34) holds in $1$D as well with $1<s\leq 2$. *

Although beyond the scope of the present work, we expect that the error analysis of a fully discrete scheme would rely on similar theoretical tools and an analogous fixed-point approach. We refer to [42, Chapter 8] which may serve as a first step in the direction of the error analysis for the Newmark time-stepping scheme that is often used in practice.

7 Numerical results

To illustrate the theory, we conduct two numerical experiments using a MATLAB implementation based on the GeoPDEs package [51].

7.1 Propagation in a channel

We first perform an experiment in a 1D channel setting. For the medium, we choose water with the parameter values

[TABLE]

cf. [29]. Following [36, Algorithm 1], to resolve the nonlinearities, we employ a fixed-point iteration with respect to the second time derivative. The tolerance is set to $\textup{TOL}=10^{-8}$. Time stepping is performed by employing the Newmark method [38] with the parameters $(\beta,\gamma)=(0.45,0.75)$. The values are chosen in this way since they were shown to provide good results in simulations of the nonlinear acoustic equations; see [16, 36]. For all spatial refinements, we have $2001$ grid points in time with the final time set to $T=37\,\mu\textup{s}$. We conduct this experiment with the initial data

[TABLE]

where $A_{1}=1.2\cdot 10^{8}\,\textup{Pa}$, $A_{2}=-10^{11}\,\textup{Pa}$, $\sigma_{1}=0.015$, $\sigma_{2}=0.02$, and $\mu=0.1$. We note that in this setting the upper bound for the acoustic pressure that guarantees non-degeneracy amounts to $1/(2k)\approx 214\,$MPa.

We use piecewise linear elements in space to compute solutions on different discretization levels. The numerical solution is computed on level $N$, $N\in[1,6]$, by employing $100\cdot 2^{N-1}$ elements for the channel length of $0.2\,\textup{m}$. The reference solution is taken to be the solution on a very fine mesh, where $N=8$. After obtaining the numerical solution on a coarse mesh, we interpolate it linearly to the mesh on level $N=8$ and compare to the reference. Figure 1 displays snapshots of the reference pressure wave as it propagates.

Let $e_{N}$ denote the error in a certain norm on level $N$. We determine the order of convergence on this level via

[TABLE]

The numerical error orders obtained in the experiments are stated in Table 1 and Table 2. We see that they agree with our theoretically predicted bounds in Theorem 6.1.

[TABLE]

[TABLE]

7.2 Focused ultrasound

In our second example, we consider a more application-oriented setting of a focused-ultrasound problem. In ultrasound applications, the sound is often excited by transducers arranged on a spherical surface [8, 29]. The wave then self-focuses as it propagates; see Figure 2. This approach results in localized high-pressure values, and is therefore often used in non-invasive treatments of kidney stones and certain types of cancer [21, 23, 31].

In the present experiment, our computational domain is a rectangle $[0,0.04]\,\textup{m}\times[0,0.05]\,\textup{m}$ with a curved bottom side that belongs to the circle centered at $(0.02\,\textup{m},0.04\,\textup{m})$ with radius $R^{2}=0.002\,\textup{m}^{2}$. On this bottom side we impose Neumann boundary conditions as a modulated sinusoidal wave:

[TABLE]

The angular frequency is taken to be $w=2\pi f$ with $f=60\,\textup{kHz}$, and we set $g_{0}=10^{7}\,\textup{Pa/m}$. On the rest of the domain sides, we impose linear absorbing boundary conditions, i.e.

[TABLE]

We mention that nonlinear absorbing conditions for the Westervelt equation have also been derived and investigated in [35, 45]. Discretization in time is performed with $3500$ time steps for the final time $T=40\,\mu\textup{s}$ and we again employ the Newmark scheme for time stepping with the same parameters as before. For the discretization is space, we employ a quadrilateral mesh which on the discretization level $N$ has $2^{N-1}\cdot 35$ elements in the propagation direction and $2^{N-1}\cdot 20$ elements in the other direction, where $N\in\{1,\ldots,5\}$.

On each discretization level $N$ we compute

[TABLE]

noting that $e_{N}=|q(u)-q(u^{N}_{h})|\leq\,\|u-u^{N}_{h}\|_{L^{\infty}L^{2}}.$ Thus we have

[TABLE]

Figure 3 displays how the function $q$ changes with respect to $h$. The data has been then fitted to a curve $\alpha+\beta h^{\gamma}$ by employing the nonlinear least-squares solver lsqcurvefit in MATLAB with the starting point $(1,1,2)$. We obtain $\gamma\approx 1.82$ for the order of convergence. In Table 3, we have the orders of discretization errors for $|q(u)-q(u_{h})|$ if we take the value on the highest level $N=5$ as the reference, i.e., $u=u_{h}^{5}$, which are again around $2$.

0.2$$0.6$$1$$1.4$$1.8$$2.2$$2.6$$3$$\cdot 10^{-3}$$1{,}176.5$$1{,}177.5$$1{,}178.5$$1{,}179.5$$1{,}180.5$$h$$q(u_{h})\ \,[\textup{Pa}]DataFitted curve

level $N$ $\,|q(u)-q(u^{N}_{h})|\,$

2 1.8386

3 2.0179

4 2.3046

Table 3: Order of error $e_{N}$.

Acknowledgments

We would like to thank Markus Muhr for helpful discussions about the numerical simulation of the problem. The funds provided by the Deutsche Forschungsgemeinschaft under the grant number WO 671/11-1 are gratefully acknowledged.