Relaxation of regularity for the Westervelt equation by nonlinear damping with application in acoustic-acoustic and elastic-acoustic coupling

In this paper we show local (and partially global) in time existence for the Westervelt equation with several versions of nonlinear damping. This enables us to prove well-posedness with spatially varying $L_\infty$-coefficients, which includes the situation of interface coupling between linear and nonlinear acoustics as well as between linear elasticity and nonlinear acoustics, as relevant, e.g., in high intensity focused ultrasound (HIFU) applications.


Introduction
The Westervelt equation for the acoustic pressure p where k = β a /λ, β a = 1 + B/(2A), λ = ̺c 2 is the bulk modulus, ̺ is the mass density, c is the speed of sound, b is the diffusivity of sound, B/A is the parameter of nonlinearity, is a classical model of nonlinear acoustics, cf. [10,16,20]. In particular, it is widely used for the simulation of high intensity focused ultrasound (HIFU) which has a broad range of technical and medical applications ranging from lithotripsy or thermotherapy to ultrasound cleaning or welding and sonochemistry, see [1,16] and the references therein. Note that the Westervelt equation can be written alternatively in terms of the acoustic velocity potential ψ, where ̺ψ t = p andk = ̺k. An analysis of the Westervelt equation with homogeneous [12] and inhomogeneous [13] Dirichlet and Neumann [15] boundary conditions as well as with boundary instead of interior damping [11] has yielded well-posedness and exponential decay of small and regular (i.e., H 2 (Ω)) solutions. Here Ω ⊆ R d , d ∈ {1, 2, 3} is the spatial domain on which (1.1) is considered.
An important open question to be addressed in this paper is existence of spatially less regular solutions as needed, e.g., for coupling with elastic or acoustic regions exhibiting different material parameters in the simulation of a focusing silicone lens immersed in an acoustic medium.
An important feature of equation (1.1) (similarly also of (1.2)) is the potential degeneracy due to the factor (1 − 2kp) of the second time derivative p tt . In order to avoid degeneracy, it is crucial to obtain an L ∞ -estimate of p in order to bound (1 − 2kp) away from zero. This has been achieved until now by combining a bound of ∆p (obtained by energy estimates) with Sobolev's embedding H 2 (Ω) → L ∞ (Ω). In this paper, we will employ the nonlinear damping for this purpose instead. Note that, indeed, the particular choice of the damping is to some extent left open from the point of view of the physical modeling. We will use this freedom to devise possible damping terms leading to existence of H 1 (Ω)-solutions.
in a smooth domain Ω ⊆ R d , d ∈ {1, 2, 3} with ε ≥ 0, δ ∈ [0, 1] and p, q ≥ 1. The constant parameters b, c will be positive while k will not be assumed to have a particular sign. Equations (1.3) and (1.4) are motivated by (1.1) while (1.5) comes from (1.2). In all equations we changed the notation to the typical mathematical one for solutions of PDEs, i.e., p → u, and ψ → u, respectively.
Note that if ε = 0 or if p = 1 in (1.3), one obtains the classical Westervelt equation; likewise, for δ = 0 or q = 1 in (1.4) and (1.5). Here we will analyze the problem for ε, δ > 0 and p, q > d − 1 , where d is the space dimension such that W 1,p+1 (Ω) and W 1,q+1 (Ω) are continuously embedded in L ∞ (Ω); precise conditions for these inclusions to hold will be discussed later in the paper. In Sections 3 and 5 we will have to to use the stronger assumption q ≥ 3, since there we will need W 1,q+1 (Ω) to be continuosly embedded in W 1,4 (Ω).
In all three cases the damping terms enable us to derive an L ∞ -estimate on u (or u t ) and thus avoid degeneracy of the coefficient 1 − 2ku (or 1 − 2ku t ). Namely, for the damping term of (1.3) we have, using homogeneity of the Dirichlet boundary values of u, dt Ω |∇u(s, y)| p+1 dy ds Finally, for (1.5), we get, replacing u by u t and p by q in (1.6) div |∇u t (s, y)| q−1 ∇u t (s, y) u tt (s, y) dy ds 1 q+1 .
Hence, in all three cases, multiplication of the damping term with u t (or u tt ) and integration over time and space will provide us with the desired L ∞ -estimate on u (or u t ). In this manner, we avoid estimates on ∆u in order to conclude L ∞ -boundedness of u, a strategy that has been used in previous studies of the Westervelt equation [12,13]. A major advantage of this relaxed regularity is the fact that it can be expected to enable solutions of the modified Westervelt equation to be coupled with other equations or with jumping coefficients, a situation that is also of high practical relevance for HIFU devices based on the use of acoustic lenses immersed in a fluid medium [16].
We will first consider an acoustic-acoustic coupling. This can be modeled by the presence of spatially varying, namely piecewise constant coefficients, in a pressure formulation (cf. [2] for the linear case) Here we have emphasized space dependence of the coefficients while suppressing space and time dependence of u in the notation as before. While in the above equations (1.3), (1.4), (1.5) and in the regions of nonlinearity (i.e., k = 0) in (1.9), strong damping b > 0 is needed for ensuring well-posedness, we may set b = 0 in regions where k vanishes. This corresponds to the physically relevant situation of a linearly acoustic (possibly to be considered as approximation to linearly elastic) silicone lens immersed in a nonlinear acoustic fluid.
The physically more relevant model requires a linearly elastic model for the lens. Therefore, we consider a velocity based formulation for elastic-acoustic coupling (see also the displacement based formulation in [4] and the velocity potential formulation in [9], both for the linear case) where u plays the role of the velocity, ψ determines the gradient part in the Helmholtz decomposition of u, u = ∇ψ + ∇ × A , and the first order differential operator B is given by Note that ψ can be determined from u as the solution of −∆ψ = −div u which is unique, e.g., if we imposed homogeneous Dirichlet boundary conditions on ψ, which we will do below. Here we think of Ω being decomposed into an acoustic (fluid) and an elastic (solid) subdomain Ω = Ω f ∪ Ω s , Ω f ∩ Ω s = ∅ . No particular smoothness assumption will have to be imposed on the domains Ω f and Ω s as these subdomains are just characterized as the sets of points where the L ∞ -coefficients ̺, c,k, and b take on certain values that are typical for fluid and solid, respectively. The acoustic region Ω f is characterized as the region of vanishing shear modulus µ = 0 in the tensor (Note that [c] could be set to any symmetric positive definite 6 × 6 matrix with entries in L ∞ (Ω) in the elastic region Ω s , thus allowing for anisotropic elasticity.) The tensor [b](x) is assumed to be symmetric nonnegative definite and to have the same structure as [c](x) in the fluid region, i.e., From the Westervelt equation on the subdomain of nonlinearity we have that the vector potential A may be set to zero in the acoustic region (ψ equals the acoustic velocity potential there) and therewith [c]Bu = λ∆ψe in Ω f , where e = (1, 1, 1, 0, 0, 0) T , so that on Ω f the PDE in (1.10) becomes Multiplying with an arbitrary vector valued test function v = ∇w +∇×W compactly supported in Ω f , integrating by parts on Ω, using the fact that e T Bv = ∆w , div v = ∆w and assuming that ̺ is constant on Ω f such that ̺∇ψ tt = ∇(̺ψ tt ), we arrive at for any smooth compactly supported w. With c 2 = λ ̺ , b =b ̺ we get (1.5) which is (up to the damping term) equivalent to (1.2).
In order to be able to make use of the embeddings 12) in this paper, we will impose zero Dirichlet boundary conditions and assume that Ω ⊂ R d is an open bounded set with Lipschitz boundary, d ∈ {1, 2, 3}.
Actually, (1.11) allows to increase the space dimension even to d = 4. However, this case is not of practical relevance. The proofs below show that the embedding H 1 0 (Ω) → L r (Ω) with r = 4 indeed suffices and we do not need to use the maximal possible exponent r = 6 in R 3 .
Throughout the paper we will use Poincaré's inequality; on convex domains the inequality reads (see for example page 364 in [14]) where u Ω is the average of u over the domain Ω and diam(Ω) denotes the diameter of the domain Ω. The Poincaré constant plays a role in the existence time for the solutions, however, our focus is on establishing local well-posedness of solutions and not on estimating the time of existence.
Moreover, this solution is unique inX and it satisfies the energy estimate and C(ǫ, r) as in (1.14).
Nonlinear p-Laplace damping was in considered in [18], where a nonlinear source term appears too. However, other features of the equation (such as variable α) were not present. For the sake of self-completeness we include below the proof of Proposition 2.1.
Proof. The weak form of (2.1) reads as with initial conditions (u 0 , u 1 ).
Step 1: Smooth approximation of α, f , and g: We consider sequences (α k ) k∈N , (f k ) k∈N and (g k ) k∈N such that (Ω))) * , • f k − 1 2 α k,t L∞(0,T ;L 2 (Ω)) ≤ b, and, for fixed k ∈ N, prove that there exists a solution u (k) of (2.6) with α, f and g replaced by α k , f k and g k , respectively, i.e., with initial conditions (u 0 , u 1 ). Later we will consider limits as k → ∞ to prove well-posedness of (2.6). The existence proof follows the line of the standard approach for linear parabolic or second order hyperbolic PDEs as it can be found, e.g., in [8]. The proof is divided into three subparts: (a) Galerkin approximation, (b) energy estimates and (c) weak limit.
Step 1 (a): Galerkin approximation. We will first show existence and uniqueness of solutions for a finite-dimensional approximation of (2.7). Assume w m = w m (x), m ∈ N are smooth functions such that where Lα k 2 is the weighted L 2 -space based on the inner product Moreover, let V n be the finite dimensional subspace of L α k 2 (Ω)∩W 1,q+1 Based on this, we consider a sequence of discretized versions of (2.7), n (t) ∈ V n and initial conditions (u 0,n , u 1,n ). For each n ∈ N the equality in (2.8) together with initial conditions (u 0,n , u 1,n ) gives an initial value problem for a second order system of ordinary differential equations which has smooth (with respect to time) coefficients and right-hand side. The standard existence theory for ordinary differential equations (cf. [19]) provides us with a unique solution u (k) n ∈ C ∞ (0,T ; V n ) of the finite-dimensional approximation (2.8) of (2.6) for someT ≤ T sufficiently small. By the uniform energy and norm estimates below, we obtaiñ T = T , i.e. there is a unique solution u (k) n ∈ C 1 (0, T ; V n ) of (2.8).
Step 1 (b): Energy estimate. Testing (2.8) with w n = u (k) n,t (t), integrating with respect to time and using the identity as well as Young's inequality in the form (1.13) with (1.14), we obtain which corresponds to the energy estimate (2.2) upon replacement of α, f, g by α k , f k , g k . As, by assumption, g k ∈ (L q+1 (0, T ; W 1,q+1 0 (Ω))) * , we have that u (k) n n∈N is a bounded sequence in the Banach space which is a reflexive Banach space. Furthermore, we obtain that (Ω)). Step n,tt . To show (ii) in case α k,t = 0 and therewithα k = α k , we now prove an analog of (2.3) with α k , f k , g k in place of α, f, g. For this purpose, let v ∈ W 1,q+1 Here holds due to the fact that with a linearly independent set of dual vectors By using orthogonality we obtain the first equality below; by (2.8) we obtain the second equality); finally, by using |. .
Multiplying the resulting inequality with a test function φ ∈ C ∞ (0, T ) and integrating with respect to time yields which shows that for some constant C > 0. Hence with a uniform bound with respect to n.
Step 3: Uniqueness of weak solutions. Uniqueness of weak solutions follows from the fact that the differenceû = u 1 − u 2 between any two weak solutions u 1 , u 2 of (2.1) satisfies Above we used the fact that for Multiplication of (2.19) byû t yields Thereforeû = 0 almost everywhere and the proof of uniqueness is complete.
Step 4: Higher energy estimate. For the proof of (iii), we need a higher energy estimate which can be obtained by testing (2.6) with w = u tt (t) (strictly speaking, we multiply by a smooth approximation of u tt and take weak limits); then integration with respect to time yields where we have used integration by parts with respect to time for the 2 , then adding (2.2) to (2.22) multiplied by λ for any τ, σ > 0 such that Note that by the assumptions α(t, x) ≥ −α > −1, c 2 > 0, (2.2) implies an estimate of the form for the usual lower order energy . Using Proposition 2.1 in a fixed point argument, we now show local well-posedness. and m sufficiently small, and u is unique in W.
Proof. We define the fixed point operator T : which is well-defined by Proposition 2.1. Indeed, by v ∈ W, (2.26), and the penultimate line in (1.7) we have that α ∈ C(0, T ; L ∞ (Ω)), with C P F the constant in the Poincaré-Friedrichs inequality Hence, we can make use of the higher energy estimate (2.4) to conclude that for any m, M > 0 with 2kC Ω we obtain that under the assumption the operator T maps into W. Contractivity is obtained by considering v i ∈ W, u i = T v i ∈ W, i = 1, 2 and subtracting the equations for u 1 and u 2 , which yields (2.28) (2.20). Due to the special form of the nonlinear strong damping term here, we cannot apply Proposition 2.1 directly, but we can proceed analogously to its proof. By multiplication of (2.28) withû t , integration with respect to space and time, and the fact that the bδ term yields a nonnegative contribution on the left hand side, we obtain This by v 1 , v 2 , u 1 , u 2 ∈ W yields where we have used |∇v| 2 L∞(0,T ;L 2 (Ω)) ≤ T |∇v t | 2 L 2 (0,T ;L 2 (Ω)) since ∇v(0) = 0, hence contractivity of T with respect to the norm |||v||| = v t L∞(0,T ;L 2 (Ω)) + ∇v L∞(0,T ;L 2 (Ω)) + ∇v t L 2 (0,T ;L 2 (Ω)) provided m is chosen sufficiently small. Like in Section 2 we assume that Ω ⊆ R d , d ∈ {1, 2, 3}, is an open bounded set with Lipschitz boundary in order to make use of (1.11), (1.12). Again we first consider an equation in which the nonlinearity occurs only through damping, namely Then (3.1) has a weak solution and any solution satisifies the energy estimate for C( bδ 2 , q+1 2 ) according to (1.14).
Proof. Since we deal with an autonomous second order in time PDE, the proof can be done directly via Galerkin discretization, energy estimates and weak limits. We therefore only provide the crucial energy estimates for (i) and (ii), which are obtained by multiplying the discretized version of (3.1) with u n,t and u n,tt , respectively and integrating by parts with respect to space and time.
For (i), from multiplication of the discretized version of (3.1) with u n,t we get 1 2 |u n,t | 2 If q > 1, using (1.13) and (1.14), we further estimate  if q > 1.
Remark 3.2. Note that multiplication with u t (via q ≥ d − 1) according to (3.2) only gives an L 2 (0, T ; L ∞ (Ω)) bound on u t and in order to obtain a C(0, T ; L ∞ (Ω)) bound, multiplication with u tt is required, cf. (1.8). Thus part (i) of Proposition 3.1 is only an intermediate result.
In part (ii) we make use of a positive sign of the term pertaining to the potential energy in the equation. However this leads to an L 4 (Ω) norm term on the right hand side of the energy identity, which we can only dominate by means of the L q+1 (Ω) norm term on the left hand side. Hence, for this part, q + 1 ≥ 4 is needed.
To this end, we define the fixed point operator T : which is well-defined by Proposition 3.1 (ii), since we have Thus by making T and the bound κ on the initial data sufficiently small, by an appropriate choice of m, M we can achieve that u ∈ W. By closedness of W we obtain existence of a solution.
Remark 3.4. To see that contractivity and therefore also uniqueness fails for (1.5), consider two solutions u i = T (v i ), i = 1, 2 as well as their differenceû = u 1 − u 2 , which is a weak solution to (2.20). Upon multiplication withû t and integration with respect to space and time, like in the proof of Proposition 2.1 and Theorem 2.3, the bδ term yields a nonnegative contribution on the left hand side. Therewith, similarly to (3.7), we obtain We see that estimation of the first term on the right hand side (the one containing v 1 tt |∇û| 2 ) would require higher spatial summability of v 1 tt and/or of ∇û. The estimates on these two quantities, following from v 1 ∈ W and from the left hand side of the above estimate onû are not strong enough for this purpose. However, multiplication of (1.5) with higher time or space derivatives of u leads to difficulties involving the strong nonlinear damping term; the same holds for multiplication of (3.13) withû tt .
Since (1.9) can be viewed as a version of (1.4) with spatially varying coefficients, inspection of the proof of Theorem 2.3 immediately yields For any T > 0 there is a κ T > 0 such that for all u 0 , there exists a weak solution u ∈ W (defined as in (2.25) with (2.26) and m sufficiently small) of (1.9) and u is unique in W.
Here the notation for velocity is motivated by the linearized Euler equation ̺v t = −∇p relating the velocity v with the pressure p, which is denoted by u here.

5.
Linear elasticity -nonlinear Westervelt equation coupling via nonlinear strong damping (1.10) In this section we consider the coupled nonlinear acoustic -elastic system We again assume that Ω ⊆ R d is an open bounded set with Lipschitz boundary and impose the following conditions on the coefficients in (5.1) , [c(x)] symmetric positive semidefinite matrices, ∃̺, ̺, c, c, γ > 0 : in Ω nl , where Ω nl := {x ∈ Ω :k(x) = 0} is an open domain and Ω nl is a compact subset of Ω or Ω has C 2 boundary .
In particular, strong linear damping is required in the whole domain, whereas strong nonlinear damping is only needed in the region Ω nl of nonlinear acoustics. Similarly to Theorem 3.3 we obtain the following local existence result Theorem 5.1. Let q ≥ 3 and assumption (5.2) be satisfied. There exist κ > 0, m > 0, a > 0, M > 0, T > 0 such that for all u 0 , u 1 ∈ (H 1 0 (Ω)) d , u 1 | Ω nl ∈ (W 1,q+1 (Ω nl )) d with we have existence of a weak solution u ∈ W ⊆ H 2 (0, T ; (L 2 (Ω)) d ) ∩ C 1 (0, T ; (H 1 0 (Ω)) d ) of (5.1) where Proof. Recall that ψ is related to u via −∆ψ = −div u and that we impose homogeneous Dirichlet boundary conditions on ψ. Like in the proof of Theorem 3.3 we can use a fixed point argument with an operator T mapping v to a weak solution of (5.4) −∆ denotes the Laplace operator on Ω with homogeneous Dirichlet boundary conditions on ∂Ω, and C ∆,Ω nl is the constant in the regularity estimate Existence of a weak solution for fixed v ∈ W can be shown like in Proposition 3.1. In particular by choosing a sufficiently small so that 2 k L∞(Ω nl ) C ∆,Ω nl C H 2 (Ω nl ),L∞(Ω nl ) a < 1 we have that for v ∈ W, provided a is chosen sufficiently small. By multiplying (5.4) with u t + εu tt we obtain We estimate the right-hand side and arrive at , thus by (5.2) using (1.13), (1.14) in Thus by choosing ε, κ sufficiently small, we obtain u ∈ W .
Due to the appearance of positive powers of t in the energy estimate (5.5), also here only a local in time well-posedness result can be expected. Morever, for similar reasons as in Section 3 uniqueness is not likely to hold here, see Remark 3.4.
Proof. The proof can be done analogously to the one of Proposition 2.1 by first showing that (6.5) below holds for smooth approximations α k , f k and g k of α, f and g by means of Galerkin's method, energy estimates and weak limits and then letting α k → α, f k → f and g k → g, respectively. Here, we omit details and only show the key energy estimates. The weak form of (6.1) reads as (6.5) with initial conditions (u 0 , u 1 ). Testing (6.5) with w = u t (t) and integrating with respect to time yields we arrive at (6.2), which directly implies (6.4). Remark 6.2. Note that in the case of p-Laplace damping we do not obtain a higher order energy estimate like in the proof of Proposition 2.1. Testing (6.5) with w = u tt yields but here the first term in the third line cannot be controlled.
Using Proposition 6.1 we now show local existence of solutions.
For this equation we even have global existence and exponential decay of the energy Theorem 6.5. Let c 2 , b, ε > 0, k ∈ R, p > d − 1.
Thus we can split the term b t 0 |∇u t (s)| 2 L 2 (Ω) ds in (6.9) as follows which inserted into (6.9) with λ ∈ (0, 1) sufficiently small implies for some constants c, C > 0 and all t > 0. This relation by a standard argument implies exponential decay of E[u](t).

Conclusions and Remarks
The introduction of nonlinear strong damping in equations of nonlinear acoustics allows us to prove existence of solutions with weaker regularity; in particular, this enables us to show well-posedness of solutions to coupled acoustic-acoustic and acoustic-elastic problems. However, the nonlinear strong damping also introduces additional challenges to the analysis: due to the relatively high order of differentiation they contain, these terms only allow to derive energy estimates for certain low order multipliers. For this reason, for some of the equations under consideration, uniqueness of solutions remains open.
The presence of the linear strong damping term −∆u t would seem to imply that the use of a nonlinear strong damping −div (|∇u t | q−1 ∇u t ) of viscosity type is more natural, it turns out that a p-Laplace damping term −div (|∇u| q−1 ∇u) yields some nicer mathematical properties such as global in time existence and exponential decay; uniqueness, however, remains an open problem for this formulation. Here we have only investigated the p-Laplace damping for the acoustic pressure formulation; we do expect that global existence and exponential decay will also carry over to the velocity potential formulation (1.5) and to the coupled problems (1.9), (1.10) upon replacement of viscosity by p-Laplace damping. This setting together with a choice of different boundary conditions (e.g., practically relevant Neumann as well as absorbing boundary conditions) will be subject of future research. R.B. and B.K. gratefully acknowledge financial support of their research by the FWF (Austrian Science Fund): P24970. The work of P.R. was supported by the NSF Grant DMS 0908435.