THE SHAPE DERIVATIVE FOR AN OPTIMIZATION PROBLEM IN LITHOTRIPSY

. In this paper we consider a shape optimization problem motivated by the use of high intensity focused ultrasound in lithotripsy. This leads to the problem of designing a Neumann boudary part in the context of the Westervelt equation, which is a common model in nonlinear acoustics. Based on regularity results for solutions of this equation and its linearization, we rigorously compute the shape derivative for this problem, relying on the variational framework from [9].


1.
Introduction. High intensity focused ultrasound (HIFU) has numerous applications ranging from industrial (ultrasound cleaning, welding, sonschemistry) to medical (thermotherapy, lithotripsy) ones. We will here concentrate on the latter, where high amplitude focused pressure waves are used to destroy kidney stones. The aim of achieving high pressure at the location of the stone, while keeping the pressure low in the surrounding tissue to avoid lesions, naturally leads to mathematical optimization problems, with some nonlinear acoustic wave equation for modeling HIFU as a constraint.
One of the most commonly used models in nonlinear acoustics is the Westervelt equation ((1 − ku)u) − c 2 ∆u − b∆u = 0 (1) where u is the acoustic pressure, c > 0 the speed of sound, b > 0 the diffusivity of sound, > 0 the mass density, and k > 0 a parameter quantifying the nonlinearity. Here and below a prime denotes the derivative with respect to time.
The two most widely used excitation and focusing principles are based on either a) excitation by a magnetomechanical principle and focusing by an acoustic lens, see Figure 1, left. b) excitation by an array of piezoelectric transducers arranged in a spherical calotte and aimed at the focus point (self-focusing), see Figure 1, right.

Figure 1. Excitation and focusing principles
Here we shall concentrate on the latter principle (b), whereas the former one (a), requiring modeling by a coupled system of PDEs for the lens and the surrounding fluid, is investigated in [10].
Excitation by the piezoelectric transducers is modeled by the Neumann boundary conditions ∂u ∂ν = g on (0, T ) × Γ n , where Γ n represents the surface composed by the piezoelectric transducers. In order to be able to restrict the open domain problem of (nonlinear) acoustics to a smooth and bounded domain Ω ⊆ R 3 (which has obvious computational advantages), we employ first order absorbing boundary conditions on the rest of the boudary Γ a = ∂Ω \ Γ n 1 c u + ∂u ∂ν = 0 on (0, T ) × Γ a For higher order absorbing boundary conditions for the Westervelt equation we refer to [11]. The boundary optimal control problem dealing with an optimal choice of the excitation pattern g in (2) has been considered in [5] and [4]. Here we will, for a given prescribed g, consider the problem of finding an optimal shape of the piezoarray, i.e., of Γ n , which is of practical relevance due to the fact that the range of achievable excitation patterns g is limited in some devices. As in [5] and [10] the goal of the optimization is to achieve a prescribed pressure distribution y d , i.e., considering the domain Ω as a design variable, we seek to minimize the cost functional with the PDE model (1), (2), (3) as a constraint. The remainder of the paper is organized as follows. Section 2 contains some results on the nonlinear PDE as well as its linearization, that will be made use of later on. In Section 3 we will provide the framework for domain variation and some useful tools for this purpose. Our main result is contained in Section 4, where we carry out a shape sensitivity analysis for the optimization problem described above SHAPE DERIVATIVE IN LITHOTRIPSY 401 and provide an expression for the Eulerian derivative of the cost functional (4) with respect to variations of the domain.
2. Some properties of the PDE model. In this section we consider the Westervelt equation on a fixed bounded C 1,1 -domain Ω ⊂ R 3 with boundary conditions (2), (3) under the following assumptions For a given function y we associate with (5) the following linear initial boundary value problem which can be regarded as a linearization of (5) with an additional inhomogeneity f . Well-posedness of (5) and (6) has already been stated in [5]. However, the argument leading from L 2 regularity of ∆u to H 3/2+s regularity of u (which together with the embedding H 3/2+s (Ω) → L ∞ (Ω) is crucial for avoiding degeneracy of the coefficient (1 − 2ku) of the second time derivative) had been missing in [5]. Theorem 1. Let T ∈ (0, ∞) and let the following conditions on the data , L ∞ (Ω)) hold and assume that y is small in the sense that i) y L ∞ (0,T ;L ∞ (Ω)) ≤ 1 4k 402 BARBARA KALTENBACHER AND GUNTHER PEICHL ii) y L ∞ (0,T ;L 2 (Ω)) ≤ min(b,c) Then the initial boundary value problem (6) has a unique weak solution and ξ satisfies the energy estimates with constants C 1 , C 2 indepentent of T . If additionally for some s ∈ (0, 1 2 ) and the compatibility conditions as well as (H4) are satisfied, then we have Proof. The proof is divided into five steps. We start with a Galerkin approximation of the weak form of the PDE by restriction to appropriate finite dimensional subspaces of H 1 (Ω), i.e., we do a semidiscretization with respect to space and establish well-posedness of the resulting system of ODEs. For these solutions, we derive energy estimates of first (Step 2) and second (Step 3) order with respect to time. These estimates enable to extract a weakly convergent subsequence in appropriate spaces, whose limit we prove to solve the weak form of the PDE in Step 4. Finally, we prove higher order regularity in space, which is crucial for avoiding degeneracy in the nonlinear equation later on. For convenience we shall use the notation a = 1 − 2ky. for which by assumption i) holds for all σ ∈ [0, T ]. Then the weak form of (6) is given by We mention that the term Γa ξ v dγ should be interpreted as Γa (tr Γa ξ ) v dγ where tr Γa is the trace operator on Γ a .
We note that all these quantities and in particular A kj (σ) and B kj (σ) are well defined. This is a consequence of a ∈ C 1 ([0, T ], L 2 (Ω)) and the embedding of H 1 (Ω) into L 4 (Ω) (with embedding constant C Ω,4 > 0) which, by the generalized Hölder inequality, implies and similarly for |(a (σ)w k , w i )| for every σ ∈ [0, T ]. For simplicity we shall write α, β, γ, F instead of α n , β n , γ n , and F n .
We note that A is positive definite for every σ ∈ [0, T ]: For any σ ∈ [0, T ], y ∈ R n \ {0} and u = n i=1 y i w i we find by (12) and the linear independence of {w 1 , . . . w n }. Thus α satisfies the following system of ordinary differential equations Since N is nonnegative, the matrix A(σ) + b c N is invertible for every σ ∈ [0, T ] and the above system is equivalent to This equation has a unique solution α ∈ H 2 ((0, T ), R n ). Since the coefficients of the equation and the right hand side at least belong to C((0, T ), R n ) we even find α ∈ C 2 ((0, T ), R n ). Hence we obtain that the equation has a unique solution ξ n ∈ C 2 (H 1 ) Step 2. lower order energy estimate. We insert w = α i w i into (17), take the sum over i and integrate from zero to s to obtain s 0 ((aξ n ) , ξ n ) + c 2 (∇ξ n , ∇ξ n ) + b|∇ξ n | 2 In view of the identity (aξ n ) ξ n = aξ n ξ n + a (ξ n ) 2 = 1 2 (a(ξ n ) 2 ) + 1 2 a (ξ n ) 2 one finds which is equivalent to which holds for all s ∈ [0, T ]. From (15) we infer . Since ξ n 0 , ξ n 1 converge to ξ 0 respectively ξ 1 in H 1 (Ω) we may assume for n sufficiently large . In view of (12) and the estimate s 0 Γn as well as the fact that defines an equivalent norm on H 1 (Ω) ([2, p 298 and Example 7.3.15], hence, for some constant C eq,Γa > 0, we obtain The choice ε = min(b,c) 4C 2 eq,Γa and assumption ii) implies

BARBARA KALTENBACHER AND GUNTHER PEICHL
), hence κ 2ε ≤ κ 2 , one eventually obtains the estimate Step 3. higher order energy estimate. In order to derive estimates for the second time derivative we test (17) with v = ξ n , which after integration by parts with respect to time of some of the terms gives (12) and the estimates which holds for all s ∈ [0, T ]. In particular, since we can bound the full H 1 norm of ξ n , so altogether (ξ n ) is bounded in . Therefore there exist a subsequence, again denoted by (ξ n ), and where the second line follows from the first one due to boundedness and linearity, hence weak continuity, of the trace operator H 1 (Ω) → H 1/2 (Γ a ), and the third line from compactness of the embedding of H 1 (H 1 ) into C(L 2 ). Additionally, it can be seen that the limit in the fourth line coincides with the derivative of the trace of ξ . Namley, denoting by (·, ·) Γa the inner product in L 2 (Γ a ), we get, for any The last equality is justified by the fact that (tr Γa ξ n ) converges weakly in L 2 (L 2 (Γ a )) to tr Γa ξ . On the other hand, by the last line in (22), Step 4. weak limits and verification of PDE. Define for arbitrary but fixed Here we have used the fact that for a ∈ L ∞ (L 2 ) ∩ H 1 (L 3 ) and ξ ∈ H 1 (L 2 ) ∩ L 2 (L 6 ) the distributional product rule is applicable and yields aξ + a ξ = (aξ ) , since with an approximating sequence ( By (22) and (23) we also have Now for any φ ∈ C([0, T ]), i ∈ N, we multiply (14) by φ and integrate over [0, T ]. This yields for any n ≥ i The preceding discussion shows that one can pass to the limit to show that ξ satisfies the equation which holds for all i ∈ N and all φ ∈ C([0, T ]). Since the elements w i span H 1 (Ω), we conclude that for any v ∈ H 1 (Ω) Next we discuss the initial condition for ξ. In view of ξ n → ξ strongly in C(L 2 ) we deduce Since ξ ∈ H 2 (L 2 ) ⊆ C 1 (L 2 ) we can conclude that ξ (0) exists in L 2 (Ω). For any φ ∈ C 1 ([0, T ]) such that φ(T ) = 0, integration by parts in the first term of (26) gives

BARBARA KALTENBACHER AND GUNTHER PEICHL
Inserting this into (26) and passing to the limit results in Integrating by parts in (27) leads to The pointwise bounds in (20) together with (21) can be interpreted in the sense that ξ n (s) and ξ n (s) are bounded in L ∞ (H 1 ) = (L 1 ((H 1 ) * )) * . Thus there exists a subsequence and elements z,z ∈ L ∞ (H 1 ) such that ξ n * z , ξ n * z .
In view of ξ n ξ in L 2 (H 1 ) (cf. (22)) we conclude z = ξ. To show that alsoz = ξ holds, for arbitrary v ∈ H 1 (Ω), φ ∈ C ∞ 0 (0, T ) we infer φ ⊗ v ∈ L 1 (H 1 ) and therefore On the other hand ξφ dσ which impliesz = ξ . Thus we also obtain ξ ∈ L ∞ (H 1 ), hence ξ satisfies the regularity (7) and the a priori bound (9). Since (28) holds for all φ ∈ C([0, T ]) and the integrand defines a function in L 2 (0, T ) for v ∈ H 1 (Ω), we deduce that ξ solves (13), i.e., (6) has a weak solution with the additional regularity (7). The a priori estimate (9) ensures the uniqueness of the solution due to linearity and the fact that the right hand side vanishes for vanishing data. Therefore not only the subsequence but also the original sequence of Galerkin approximations (ξ n ) n∈N converges to ξ in the sense of (22).

BARBARA KALTENBACHER AND GUNTHER PEICHL
This yields Elliptic regularity results imply that (37) admits a solution z(σ) ∈ H s+3/2 (Ω) which is unique up to a constant (which may depend on σ). This constant is determined by the fact that ξ is a solution of (36), hence by (31) we have z(σ) = ξ(σ). Moreover, we have the estimate where F is determined by the right hand side in (37) and where the constant C ell > 0 does not depend on σ. Indeed, we can employ the Exact Interpolation Theorem [1,Theorem 7.23] together with boundedness of the solution operator S : to conclude that S is also bounded between the interpolated spaces X s = H s (∂Ω) and Y s = H 3/2+s (Ω). The estimate (39) then follows by decomposition of ξ as ξ(σ) = Sϕ(σ)+ ξ F (σ), where ξ F (σ) solves (37) with homogeneous Neumann boundary conditions (and a σ dependent constant such that (36) is satisfied).
Since H s+3/2 (Ω) embeds continuously into C(Ω) with constant C Ω,s+ 3 2 we obtain the estimate ξ L ∞ (C(Ω)) ≤ C Ω,s+ 3 2 C ell ( F L ∞ (L 2 (Ω)) + ϕ L ∞ (H s (∂Ω)) ) Using (12), one can verify the estimate Combining the above estimates with (9) and (21) one eventually finds is sufficiently small. The solution u satisfies the a priori estimates with X as in (7) Proof. We have shown that for every y satisfying the assumptions of Theorem 1 the linearized system (6) has a unique solution ξ = Ty. In view of it is clear that a fixed point of T is a solution of (5). We now verify the assumptions of the contraction mapping principle. Let y andỹ satisfy the assumptions of Theorem 1 and denote the corresponding solutions by ξ andξ. Then the difference z =ξ − ξ satisfies whereã = 1 − 2kỹ. Since the right hand side of this equation just belongs to L 2 ((H 1 (Ω)) * ), we can only use the lower order energy estimate (8) to obtain For almost all σ ∈ [0, T ] we obtain Similarly one obtains which leads to ). Combining this estimate with (44) and one can in particular derive the estimate (45) This suggests to work with the Banach spacê Then, provided y satisfies the conditions of the first part of Theorem 1, i.e., y ∈ , estimate (45) and the a priori estimate (9) with (21) give rise to H 1 (H −1/2 (Γn)) ) which shows that T is a contraction if we impose some fixed bound r > 0 on ỹ L ∞ (L 3 ) and if is sufficiently small. This argument suggests to define the set as the domain of the fixed point operator T. For y ∈ W the a priori bounds (9), (21) and (40) imply (1 + T )(1 + r 2 )r 2 0 , provided additionally u 0 H 2 is sufficiently small. Hence by boundedness of the embedding H 1 (Ω) → L 3 (Ω), the set W is invariant under T if r 0 is sufficiently small.

SHAPE DERIVATIVE IN LITHOTRIPSY 417
For h ∈ D and t ∈ R sufficiently small define Then F t is a C 1,1 -diffeomorphism which preserves the regularity of Ω. One defines the Eulerian derivative of J by where u and u t satisfy the PDE on the original and on the perturbed domain, respectively E(u, Ω) = 0, E(u t , Ω t ) = 0. The constraint E(u, Ω) = 0 stands for the weak form of the PDE with initial and boundary conditions, which is in our case The expression dJ(Ω) is called shape derivative of J if dJ(Ω)h exists for all h ∈ D and dJ(Ω) ∈ D * . From the Delfour-Hadamard-Zolesio structure theorem we expect a representation of the form to hold, which we will indeed establish in Theorem 8 below. The difficulty ensuing from the fact that the difference quotient in (47) involves functions defined on different domains can be overcome with the method of mapping. Let For convenience we recall some rules for differentiation and integration of the mapped functions, see, e.g., [6].
4. If f ∈ C((−τ, τ ), W 2,1 (U )) and f t (0) exists in W 1,1 (U ), then Replacing Ω by Ω t in (48) defines a family of perturbed problems which by Theorem 2 have a unique solution u t ∈ L ∞ (0, T ; L ∞ (Ω t )) ∩ H 2 (0, T ; L 2 (Ω t )) ∩ W 1,∞ (0, T ; H 1 (Ω t )). Using the above transformation rules and observing ω t ≡ 1 on Γ a one can verify that the transformed solution u t satisfies Theorem 5. Equation (49) has a unique solution for t sufficiently small. This family of solutions satisfies the a priori bounds (41) and (42) uniformly in t. Moreover we have for u 0 , u 1 and g sufficiently small where u is the solution of (5) and Proof. The first part of the Theorem is a straightforward consequence of the fact that u t is the diffeomorphic image of the corresponding solution of the Westervelt equation (5) defined on Ω t . Since the constants appearing in the a priori estimates for u t can be bounded uniformly with respect to t one can argue that the a priori bounds also apply to the family of transformed solutions u t for t sufficiently small. Next we turn to the analysis of the dependence of u t on t. Let z t stand for the difference Subtracting (48) from (49) one can verify that z t satisfies the equation i.e., the weak form (13) of (6) with By Theorem 2, the assumptions of Theorem 1 are satisfied, hence estimates (8) and (21) yield Here we can estimate for any v ∈ H 1/2 (Γ n ) and ) are sufficiently small and therewith, by (41), we get (50) from (52).
Finally, we provide an identity that will be used in the analysis in the next section.
Lemma 6. For u, v ∈ M where M is the linear space induced by the norm for some s ∈ (0, 1 2 ), the following identity holds

Proof.
A proof of this result for u, v ∈ H 2 (Ω) can be found in [3]. The assertion for u, v ∈ M follows from the following result on density of C ∞ (Ω) in M. Let ε > 0 and w ∈ M be fixed. then we can choose f ∈ C ∞ (Ω), g ∈ C 0,1 (∂Ω) = W 1,∞ (∂Ω) ⊆ H 1/2 (∂Ω) (note that the smoothness of g is limited by the smoothness of ∂Ω) such that with c 1 , c 2 > 0 to be chosen independently of ε below. We define z as the solution to the elliptic Neumann problem Elliptic regularity and the fact that f ∈ L 2 (Ω), g ∈ H 1/2 (∂Ω), yields z ∈ H 2 (Ω), so that we can choose ψ ∈ C ∞ (Ω) such that Moreover, since the difference w − z satisfies again by elliptic regularity, as well as (53), we get where C ∆ Ω denotes the constant in the elliptic regularity estimate and C L 2 →H −1/2+s Ω the embedding constant of L 2 (Ω) → H −1/2+s (Ω). For the difference of the Laplace operator values we get for the difference of the normal derivatives we have These estimates together with (54) and a sufficiently small choice of c 1 , c 2 , c 3 > 0 yield w − ψ M < ε. Moreover, by elliptic regularity, for any w ∈ M we have for some positive constant C independent of ε. Thus due to Φ 1 (ũ,ṽ) − Φ 2 (ũ,ṽ) = 0 ))| < 2Cε, hence, since ε > 0 was arbitrary, Φ 1 (u, v) = Φ 2 (u, v).

4.
The shape derivative of J. In this section we turn to the sensitivity analysis of the cost functional with the PDE model (5) as a constraint. In accordance with reality and to simplify the presentation we assume u(0) = u (0) = 0.
Furthermore we also assume that the excitation of the piezoelectric transducers on Γ n is determined by the trace of a globally defined function for some s ∈ (0, 1 2 ) and y d ∈ L 2 (0, T ; L 2 (U )) ∩ H 1 (0, T ; (H 1 (U )) * ) Following the strategy in [9] we introduce the adjoint system which is the weak form (13) of the linearized system (6) considered in Section 2, with a right hand side depending onp. Hence, using Theorem (1) and a fixed point argument, one can show well-posedness of the adjoint equation.
Corollary 7. Under the assumptions of Theorem 2, the adjoint equation (56) has a unique solution p ∈ X with X as in (7).
For further reference we note that the state u t and p also satisfy the integrated equations (c 2 g t + b(g t ) )φω t dγ dσ, φ ∈ L 2 (0, T ; H 1 (Ω)), Note that u satisfies (58) with t = 0. and prove q a ∈ H s (D) with q a H s (D) ≤ q a H s (Da) , and likewise for q n , q n , which by the triangle inequality will yield the assertion. Fix > 0. By [8,Corollary 1.4.4.5.], which states that H s 0 (D a ) = H s (D a ) for s ∈ (0, 1 2 ) and D a an open Lipschitz domain, there exists a sequence (q a,k ) k∈N ⊆ C ∞ 0 (D a ) such that q a,k − q a H s (Da) → 0 as k → ∞.
For (q a k ) k∈N defined by q a k = q a,k on D a 0 on D n we have q a k H s (D) = q a,k H s (Da) ≤ q a H s (Da) + for k sufficiently large, hence there exists a subsequence, again denoted by (q a k ) k∈N , and an elementq a ∈ H s (D) such that q a k H s q a and q a k L 2 q a .
Since {χ ∈ C ∞ 0 (D) : supp(χ) ⊆ D a } + {χ ∈ C ∞ 0 (D) : supp(χ) ⊆ D n } is dense in L 2 (D), this implies q a =q a . As > 0 was arbitrary, this proves the assertion. The proof is the same for q n , with the only difference that due to the assumption q n ∈ H s 0 (D n ), we need not (and cannot, since D n in general cannot be open as well) employ [8,Corollary 1.4.4.5.].