On the First-Order Shape Derivative of the Kohn-Vogelius Cost Functional of the Bernoulli Problem

Abstract

The exterior Bernoulli free boundary problem is being considered. The solution to the problem is studied via shape optimization techniques. The goal is to determine a domain having a specific regularity that gives a minimum value for the Kohn-Vogelius-type cost functional while simultaneously solving two PDE constraints: a pure Dirichlet boundary value problem and a Neumann boundary value problem. This paper focuses on the rigorous computation of the first-order shape derivative of the cost functional using the Hölder continuity of the state variables and not the usual approach which uses the shape derivatives of states.

1. Introduction

The Bernoulli problem is the prototype of a stationary free boundary problem. It arises in various applications such as electrochemical machining, potential flow in fluid mechanics, tumor growth, optimal insulation, molecular diffusion, and steel and glass production [1–6]. A characteristic feature of free boundary problems is that not only the state variable is unknown but also the domain on which the state equation is posed. This represents a significant theoretical as well as numerical challenge. One can characterize the Bernoulli problem, at least along general lines, by finding a connected domain as well as a function which is harmonic on this domain. One component on the boundary is known. The other one is determined by a set of overdetermined boundary conditions (a Dirichlet condition and a Neumann condition) for the state. If the free boundary component is strictly exterior to the fixed part of the boundary, the problem is called exterior Bernoulli problem and interior Bernoulli problem otherwise. For more discussions related to interior and exterior Bernoulli problems, we refer the reader to [1, 4, 7–10].

Recent strategies to compute a numerical solution are based on reformulating the Bernoulli problem as a shape optimization problem. This can be achieved in several ways. For a given domain, one can choose one of the boundary conditions on the free boundary to obtain a well-posed state equation. The domain is determined by the requirement that the other condition on the free boundary is satisfied in a least squares sense (cf. [11–13]). Alternatively, one can compute on a given domain two auxiliary states: 𝑢_𝐷 which satisfies the Dirichlet condition and 𝑢_𝑁 which satisfies the Neumann condition on the free boundary. The underlying domain is selected such that the difference 𝐽(Ω)=|𝑢_𝐷−𝑢_𝑁|²_𝐻1(Ω) is as small as possible. In fact, if 𝐽(Ω)=0 for a domain Ω then 𝑢_𝐷=𝑢_𝑁 and (𝑢_𝐷,Ω) is a solution of the Bernoulli problem. Sometimes 𝐽 is called Kohn-Vogelius functional since Kohn and Vogelius were among the first who used such a functional in the context of inverse problems [14]. Standard algorithms to minimize 𝐽 require some gradient information. So in this paper, the first-order sensitivity analysis is carried out for the functional 𝐽 for the exterior Bernoulli problem. The main contribution in this paper is the application of a shape optimization technique that leads to the explicit expression for the shape derivative of the cost functional. This is done through variational means similar to the techniques developed in [9, 10, 13], wherein we use the Hölder continuity of the state variables satisfying the Dirichlet and Neumann problems but we do not introduce any adjoint variables. In our approach, we also bypass the use of the material derivatives of the states (which was done in [1]) and the use of states’ shape derivatives.

The rest of the paper is structured as follows. Section 2 presents the Bernoulli free boundary problem and its shape optimization formulations. Section 3 provides a list of shape optimization tools that are needed in the analysis for the shape derivatives of the Kohn-Vogelius cost functional 𝐽. Section 4 presents an exhaustive discussion on the first-order shape derivative of 𝐽. Finally, Section 5 draws conclusion and observation.

2. The Bernoulli Problem

The exterior Bernoulli free boundary problem is formulated as follows. Given a bounded and connected domain 𝐴⊂ℝ² with a fixed boundary Γ∶=𝜕𝐴 and a constant 𝜆<0, one needs to find a bounded connected domain 𝐵⊂ℝ² with a free boundary Σ, containing the closure of 𝐴, and an associated state function 𝑢∶Ω→ℝ, where Ω=𝐵⧵𝐴, such that the overdetermined conditions are satisfied: −Δ𝑢=0inΩ,𝑢=1onΓ,𝑢=0onΣ,𝜕𝑢𝜕𝐧=𝜆onΣ. (1)

On the other hand, the interior Bernoulli free boundary problem has the following formulation. Given a bounded and connected domain 𝐵⊂ℝ² with a fixed boundary Γ∶=𝜕𝐵 and a constant 𝜆>0, one determines a bounded connected domain 𝐴⊂𝐵 with a free boundary Σ and an associated state function 𝑢∶Ω→ℝ, where Ω=𝐵⧵𝐴, subject to the following constraints: −Δ𝑢=0inΩ,𝑢=0onΓ,𝑢=1onΣ,𝜕𝑢𝜕𝐧=𝜆onΣ. (2) In both problems 𝐧 is the outward unit normal vector to Σ. The difference in the domains of these two types of Bernoulli problems is depicted in Figure 1.

(a)

(b)

Figure 1 

The domain Ω for the interior Bernoulli problem (a) and exterior problem (b).

Methods of shape optimization can be employed in solving the exterior Bernoulli free boundary problem (1). As we observe, this boundary problem is ill-posed due to the fact that we have overdetermined conditions on the free boundary Σ. So to overcome the difficulty of solving it, one can reformulate it as one of the following shape optimization problems which involves now a well-posed state equation.(1)Tracking Neumann data [11, 12] as min_Σ𝐽₁(Σ)≡min_Σ12_Σ𝜕𝑢_𝐷𝜕𝐧−𝜆²𝑑𝑠, (3) where the state function 𝑢_𝐷 is the solution to the Dirichlet problem −Δ𝑢_𝐷𝑢=0inΩ,_𝐷𝑢=1onΓ,_𝐷=0onΣ. (4)(2)Tracking Dirichlet data [11, 13] as min_Σ𝐽₂(Σ)≡min_Σ12_Σ𝑢²_𝑁𝑑𝑠, (5) where the state function 𝑢_𝑁 is the solution to the Neumann problem −Δ𝑢_𝑁𝑢=0inΩ,_𝑁=1onΓ,𝜕𝑢_𝑁𝜕𝐧=𝜆onΣ. (6)(3)Minimizing the Kohn-Vogelius type cost functional [12, 15] as min_Ω𝐽(Ω)≡min_Ω12_Ω∇𝑢_𝐷−𝑢_𝑁²𝑑𝑥, (7) where state functions 𝑢_𝐷 and 𝑢_𝑁 satisfy (4) and (6), respectively.

In this paper, we are just interested in the study of minimizing the Kohn-Vogelius functional 𝐽.

3. Tools in Shape Optimization

3.1. Feasible Domain Ω

In this work, we are interested in 𝐶^𝑘,1-domains, where 𝑘≥0. Aside from being 𝐶^𝑘,1 we also assume that these are bounded and connected subsets of a bigger set 𝑈 which is also a bounded connected 𝐶^𝑘,1 domain. This 𝑈 is called the universal or the hold-all domain. The smoothness of these domains can be defined in the following sense (cf. [16]).

Consider the standard unit orthonormal basis {𝑒₁,𝑒₂,…,𝑒_𝑛} in ℝ^𝑛. For a point 𝑥=(𝑥₁,𝑥₂,…,𝑥_𝑛)∈ℝ^𝑛, let 𝑥^′=(𝑥₁,𝑥₂,…,𝑥_𝑛−1)∈ℝ^𝑛−1 so as to write 𝑥=(𝑥^′,𝑥_𝑛). Consider the unit ball 𝐵(0,1) and introduce the subsets 𝐵⁺(0,1)∶=𝑥∈𝐵(0,1)∶𝑥_𝑛,𝐵>0₋(0,1)∶=𝑥∈𝐵(0,1)∶𝑥_𝑛,𝐵<0₀(0,1)∶=𝑥∈𝐵(0,1)∶𝑥_𝑛=0. (8)

Definition 1. A domain Ω⊂ℝ^𝑛 with a nonempty boundary 𝜕Ω is called a 𝐶^𝑘,𝑙-domain, where 0≤𝑘,  0<𝑙≤1, if for every 𝑦∈𝜕Ω there exists a neighborhood 𝑁_𝑦 of 𝑦 and a 𝐶^𝑘,𝑙 diffeomorphism 𝑓_𝑦∶𝑁_𝑦→𝐵(0,1) such that (a)𝑓_𝑦(𝑁_𝑦∩Ω)=𝐵⁺(0,1), (b)𝑓_𝑦(𝑁_𝑦∩𝜕Ω)=𝐵⁰(0,1), and (c)𝑓_𝑦(𝑁_𝑦∩Ω^𝑐)=𝐵⁻(0,1).

To illustrate this for 𝑛=2 and 𝑘=𝑙=1, see Figure 2.

Figure 2 

A 𝐶𝑘,1-domain Ω, where 𝑓𝑦 is a diffeomorphism from the neighborhood 𝑁𝑦 to the ball 𝐵(0,1).

Note that if Ω is a bounded, open, connected set with a 𝐶^0,1 boundary, then intΩ=Ω. This was given in [17] and we prove it as follows.

Theorem 2 (see [17]). If Ω is a bounded open connected subset of ℝ^𝑛 with Lipschitz continuous boundary, then intΩ=Ω.

Proof. The interior of Ω is the largest open set contained in the set Ω. Moreover, Ω⊆Ω. It follows that Ω⊆intΩ. Next, we show that intΩ⊆Ω. Clearly, intΩ⊆Ω. We now show that if 𝑥∈intΩ, then 𝑥∉𝜕Ω.
Suppose 𝑥∈𝜕Ω and 𝑥∈intΩ. We need to show that any open set 𝑁 containing 𝑥 contains an element not in Ω. We first note that by definition of 𝐶^0,1 domain, there exists a neighborhood 𝑁_𝑥 of 𝑥∈𝜕Ω and a diffeomorphism 𝑓_𝑥∶𝑁_𝑥→𝐵(0,1). Let 𝑁 be an open set containing 𝑥 with 𝑁⊂𝑁_𝑥. It follows that 𝑓_𝑥(𝑁) is an open set containing 0 and this set is contained in 𝐵(0,1). Hence, there exists 𝑧∈𝑓_𝑥(𝑁) such that 𝑧∈𝐵⁻(0,1). This implies that 𝑓_𝑥⁻¹(𝑧)∈𝑁∩Ω^𝑐. Thus, 𝑁 contains an element not in Ω, which is a contradiction. Therefore, 𝑥∉intΩ. We have proven that if 𝑥∈𝜕Ω, then 𝑥∉intΩ. Taking the contrapositive of this statement we get that if 𝑥∈intΩ, then 𝑥∉𝜕Ω. Since 𝑥∈Ω but 𝑥∉𝜕Ω, we conclude 𝑥∈Ω. Thus, intΩ⊆Ω. We have shown that Ω⊆intΩ and intΩ⊆Ω. Therefore, intΩ=Ω.

3.2. The Perturbation of Identity Technique

Given bounded connected domains Ω and 𝑈 of ℝ², where Ω⊆𝑈, and a linear space Θ of vector fields 𝐕, one can deform Ω via the perturbation of identity operator 𝑇_𝑡∶𝑈⟶ℝ²,𝑇_𝑡(𝑥)=𝑥+𝑡𝐕(𝑥),𝑥∈𝑈, (9) where 𝐕∈Θ. For a given 𝑡 we denote the deformed domain to be Ω_𝑡, which is the image of Ω under 𝑇_𝑡.

Throughout the paper, we use the usual infinity norms in the spaces 𝐶(𝑋;ℝ), 𝐶(𝑋;ℝ²), and 𝐶(𝑋;ℝ^2×2), where 𝑋 is a compact subset of ℝ². In addition to this, we also denote the Frobenius norm of 𝑀(𝑥) to be |𝑀(𝑥)|_𝐹=²_𝑖,𝑗=1|𝑀_𝑖𝑗(𝑥)|^21/2. (10) This norm and the infinity norm of the matrix 𝑀 can be related as |𝑀(𝑥)|_𝐹≤2|𝑀|_∞,𝑥∈𝑋. (11) This can be shown easily. One can also show that if 𝑀∈𝐶(𝑋;ℝ^2×2) and 𝑦∈𝐿²(𝑋;ℝ²), then the vector 𝑀𝑦 is bounded in 𝐿²(𝑋;ℝ²). In fact, 𝑀𝑦_𝐿2(𝑋)≤2|𝑀|_∞𝑦_𝐿2(𝑋) (12) and the proof is trivial. Finally, the symbols |⋅| or |⋅|₂ will refer to the usual Euclidean norm.

The Perturbed Domain Ω_𝑡. The domains Ω_𝑡 that are considered in this work are of annulus type with boundary 𝜕Ω_𝑡, which is the union of two disjoint sets Γ_𝑡 and Σ_𝑡, referred to as the fixed and free boundaries, respectively. These domains are obtained through the operator defined in (9), where 𝐕 belongs to Θ, which is defined as Θ=𝐕∈𝐶^1,1𝑈,ℝ²∶𝐕|_{Γ∪𝜕𝑈}=0. (13) For 𝑡=0, we obtain the reference domain Ω∶=Ω₀, with a fixed boundary Γ∶=Γ₀ and a free boundary Σ∶=Σ₀.

The main objective in this subsection is to show that 𝑇_𝑡 is a diffeomorphism from Ω to Ω_𝑡 for sufficiently small 𝑡. To verify this, we need the following results, which are given and proven in [17].

Theorem 3. If Ω is a bounded, open, connected set in ℝ^𝑛 such that intΩ=Ω and 𝜑 is a continuous injective mapping from Ω to ℝ^𝑛, then 𝜑Ω=𝜑(Ω),𝜑(Ω)=int𝜑Ω,𝜑(𝜕Ω)=𝜕𝜑(Ω)=𝜕𝜑Ω. (14)

Theorem 4. Suppose(1)Ω is a bounded, open, connected set in ℝ^𝑛 such that intΩ=Ω,(2)𝜑₀∈𝐶(Ω,ℝ²) where 𝜑₀ is injective,(3)𝜑∈𝐶(Ω,ℝ²)∩𝐶¹(Ω,ℝ²) such that (𝑥)𝜑(𝑥)det𝐷𝜑>0∀𝑥∈Ω,=𝜑₀(𝑥)∀𝑥∈𝜕Ω. (15)Then (i)𝜑∶Ω→𝜑(Ω) is a homeomorphism (i.e., 𝜑 is a bijection, 𝜑 is continuous, and 𝜑⁻¹ is continuous),(ii)𝜑∶Ω→𝜑(Ω) is a 𝐶¹-diffeomorphism (i.e., 𝜑 is a bijection, 𝜑∈𝐶¹(Ω), 𝜑⁻¹∈𝐶¹(𝜑(Ω))),(iii)𝜑(Ω)=𝜑₀(Ω), 𝜑(Ω)=𝜑₀(Ω).

We also consider the following property of a domain, which is also found in [17, page 52].

Lemma 5. If Ω is a bounded, open, connected subset of ℝ^𝑛 having a Lipschitz continuous boundary, then there is a number 𝑐(Ω) such that, for any given points 𝑥,𝑦∈Ω, one can find a finite sequence of points 𝑦_𝑘, 𝑘=1,…,𝑙+1, satisfying the following properties: (a)𝑦₁=𝑥, 𝑦_𝑘∈Ω for 2≤𝑘≤𝑙, 𝑦_𝑙+1=𝑦,(b)(𝑦_𝑘,𝑦_𝑘+1)⊂Ω for 1≤𝑘≤𝑙,(c)∑^𝑙_𝑘=1|𝑦_𝑘+1−𝑦_𝑘|≤𝑐(Ω)|𝑥−𝑦|.

We also recall the useful property of the determinant of the Jacobian of 𝑇_𝑡 which is given in the next lemma. Here we use the notation 𝐼_𝑡(𝑥)=det𝐷𝑇_𝑡(𝑥),𝑥∈𝑈. (16)

Lemma 6 (see [9, 13]). Consider the operator 𝑇_𝑡 defined by (9), where 𝐕∈Θ, which is described by (13). Then (i)𝐼_𝑡=1+𝑡div𝐕+𝑡²det𝐷𝐕,(ii)there exist 𝑡_𝑉,𝛼₁,𝛼₂>0 such that 0<𝛼₁≤𝐼_𝑡(𝑥)≤𝛼₂, for |𝑡|≤𝑡_𝑉, 𝑥∈𝑈.

Proof. In general, for 𝑛-dimensional case, the Jacobian of 𝑇_𝑡 is given by 𝐷𝑇_𝑡=(𝑎_𝑖𝑗), where 𝑎_𝑖𝑗=𝑡(𝜕𝑉_𝑖/𝜕𝑥_𝑗) if, 𝑖≠𝑗, and 𝑎_𝑖𝑖=1+𝑡(𝜕𝑉_𝑖/𝜕𝑥_𝑖). By definition of the determinant, we can write det𝐷𝑇_𝑡 as det𝐷𝑇_𝑡=_{𝜎∈𝑆𝑛}sgn(𝜎)𝑎_1,𝜎(1)𝑎_2,𝜎(2)⋯𝑎_{𝑛,𝜎(𝑛)}=^𝑛_𝑖=11+𝑡𝜕𝑉_𝑖𝜕𝑥_𝑖+_{𝜎∈𝐹𝑛⧵{𝑖}}sgn(𝜎)𝑎_1,𝜎(1)𝑎_2,𝜎(2)⋯𝑎_{𝑛,𝜎(𝑛)}+_{𝜎∈𝑆𝑛⧵𝐹𝑛}sgn(𝜎)𝑎_1,𝜎(1)𝑎_2,𝜎(2)⋯𝑎_{𝑛,𝜎(𝑛)}=∶𝐼+𝐼𝐼+𝐼𝐼𝐼, (17) where 𝑆_𝑛 refers to the set of all permutations of {1,2,…,𝑛}, 𝑖 is the identity permutation, 𝐹_𝑛={𝜎∈𝑆_𝑛∶𝜎(𝑘)=𝑘forsomepositiveinteger𝑘≤𝑛}, and sgn(𝜎) is either 1 (if the number of inversions is even) or −1 (if the number of inversions is odd). We observe that the expression 𝐼 can be written as 𝐼=1+𝑡div𝑉+𝑡²𝑅₁(𝑡,𝐕), where 𝑅₁∈𝐶(ℝ,𝐶^0,1(𝑈)). We also observe that, for 𝑛≥2, each term of the expression 𝐼𝐼 has at least 2 factors that are of the form 𝑡(𝜕𝑉_𝑖/𝜕𝑥_𝑗), 𝑖≠𝑗. Hence we can write 𝐼𝐼=𝑡²𝑅₂(𝑡,𝐕), where 𝑅₂ is in 𝐶(ℝ,𝐶^0,1(𝑈)). All terms of 𝐼𝐼𝐼 have factors of the form 𝑡(𝜕𝑉_𝑖/𝜕𝑥_𝑗),  𝑖≠𝑗, and thus we have 𝐼𝐼𝐼=𝑡^𝑛𝑅₃(𝐕), which can be written as 𝑡²𝑅₃(𝑡,𝐕), where 𝑅₃∈𝐶(ℝ,𝐶^0,1(𝑈)). Combining 𝐼, 𝐼𝐼, and 𝐼𝐼𝐼, we get det𝐷𝑇_𝑡=1+𝑡div𝐕+𝑡²𝑅(𝑡,𝐕) (18) with 𝑅∈𝐶(ℝ,𝐶^0,1(𝑈)). In particular, for 𝑛=2, the determinant is computed as follows: det𝐷𝑇_𝑡=1+𝑡div𝐕+𝑡²det𝐷𝐕. (19)
This verifies (i). To show (ii) we first get the lower bound for 𝐼_𝑡(𝑥). Take 𝑡_𝑉1=min2|div𝐕|_𝐶(𝑈)+|det𝐷𝐕|_𝐶(𝑈),1. (20) For |𝑡|≤𝑡_𝑉, we obtain 𝐼_𝑡(𝑥)≥1−|𝑡||div𝐕|_𝐶(𝑈)−|𝑡|²|det𝐷𝐕|_𝐶(𝑈)≥1−|𝑡||div𝐕|_𝐶(𝑈)−|𝑡||det𝐷𝐕|_𝐶(𝑈)=1−|𝑡||div𝐕|_𝐶(𝑈)+|det𝐷𝐕|_𝐶(𝑈)≥12. (21) On the other hand, by triangle inequality we have 𝐼_𝑡(𝑥)≤1+|𝑡||div𝐕|_𝐶(𝑈)+|𝑡|²|det𝐷𝐕|_𝐶(𝑈)≤1+𝑡_𝑉|div𝐕|_𝐶(𝑈)+𝑡²_𝑉|det𝐷𝐕|_𝐶(𝑈)∀𝑥∈𝑈. (22) Hence, we have shown that there are positive constants 𝛼₁=1/2 and 𝛼₂=1+𝑡_𝑉|div𝐕|_𝐶(𝑈)+𝑡²_𝑉|det𝐷𝐕|_𝐶(𝑈) such that 𝛼₁≤𝐼_𝑡(𝑥)≤𝛼₂ for 𝑥∈𝑈.

Considering the theorems and lemmas presented beforehand, we are now ready to prove the following theorem.

Theorem 7. Let Ω and 𝑈 be nonempty bounded open connected subsets of ℝ² with Lipschitz continuous boundaries, such that Ω⊆𝑈, and 𝜕Ω is the union of two disjoint boundaries Γ and Σ. Let 𝑇_𝑡 be defined as in (9) where 𝐕 belongs to Θ, defined as (13).
Then for sufficiently small 𝑡, (1)𝑇_𝑡∶𝑈→𝑈 is a homeomorphism,(2)𝑇_𝑡∶𝑈→𝑈 is a 𝐶^1,1 diffeomorphism, and in particular, 𝑇_𝑡∶Ω→Ω_𝑡 is a 𝐶^1,1 diffeomorphism,(3)Γ_𝑡=𝑇_𝑡(Γ)=Γ,(4)𝜕Ω_𝑡=Γ∪𝑇_𝑡(Σ).

Proof. First, because 𝑈 is a 𝐶^0,1 domain, it follows that int𝑈=𝑈 by Theorem 2. Second, 𝑇₀=𝐼∈𝐶(𝑈,ℝ²), and it is injective. Third, it is evident that 𝑇_𝑡 is 𝐶^1,1 because 𝐕 is 𝐶^1,1. For 𝑥∈𝜕𝑈, 𝑇_𝑡(𝑥)=𝑥 because 𝐕 vanishes on 𝜕𝑈. For 𝑥∈𝑈, the determinant of the Jacobian of the perturbation of identity operator 𝑇_𝑡 is given by (19). By Lemma 6, there exists a 𝑡_𝑉>0, given by (20), such that det𝐷𝑇_𝑡(𝑥)>0 for all 𝑥∈𝑈 and for |𝑡|≤𝑡_𝑉. Hence, by applying Theorem 4, we conclude that 𝑇_𝑡(𝑈)=𝑈 and 𝑇_𝑡(𝑈)=𝑈 for all |𝑡|<𝑡_𝑉, and 𝑇_𝑡∶𝑈→𝑈 is a homeomorphism. Furthermore, by Theorem 4, we find that 𝑇_𝑡∶𝑈→𝑈 is a 𝐶¹ diffeomorphism. To show that 𝑇_𝑡 is a 𝐶^1,1 diffeomorphism, we are left to show that 𝐷𝑇_𝑡⁻¹ is Lipschitz continuous. To verify this we use Lemma 5.
Given any two points 𝑢,𝑣∈𝑈 we choose {𝑦_𝑘} such that properties (a)–(c) of Lemma 5 are satisfied. For fixed |𝑡|<𝑡_𝑉, differentiating the identities 𝑇_𝑡∘𝑇_𝑡⁻¹=𝐼 and 𝑇_𝑡⁻¹∘𝑇_𝑡=𝐼 will lead to 𝐷𝑇_𝑡𝑇_𝑡⁻¹(𝑧)𝐷𝑇_𝑡⁻¹(𝑧)=𝐼,𝐷𝑇_𝑡⁻¹(𝑧)𝐷𝑇_𝑡𝑇_𝑡⁻¹(𝑧)=𝐼, (23) for all 𝑧∈𝑈. Thus, 𝐷𝑇_𝑡𝑇_𝑡⁻¹𝑦_𝑘+1𝐷𝑇_𝑡⁻¹𝑦_𝑘+1=𝐼,𝐷𝑇_𝑡⁻¹𝑦_𝑘𝐷𝑇_𝑡𝑇_𝑡⁻¹𝑦_𝑘=𝐼. (24)
This implies 𝐷𝑇_𝑡⁻¹(𝑢)−𝐷𝑇_𝑡⁻¹(𝑣≤_𝑙𝑘=1)𝐷𝑇_𝑡⁻¹𝑦_𝑘−𝐷𝑇_𝑡⁻¹𝑦_𝑘+1=^𝑙_𝑘=1𝐷𝑇_𝑡⁻¹𝑦_𝑘×𝐷𝑇_𝑡𝑇_𝑡⁻¹𝑦_𝑘+1−𝐷𝑇_𝑡𝑇_𝑡⁻¹𝑦_𝑘×𝐷𝑇_𝑡⁻¹𝑦_𝑘+1. (25)
Applying the infinity norm in the space 𝐶(𝑋;ℝ²) we have 𝐷𝑇_𝑡⁻¹(𝑢)−𝐷𝑇_𝑡⁻¹(𝑣≤_𝑙𝑘=1)𝐷𝑇_𝑡⁻¹𝑦_𝑘𝐷𝑇_𝑡⁻¹𝑦_𝑘+1×𝐷𝑇_𝑡𝑇_𝑡⁻¹𝑦_𝑘+1−𝐷𝑇_𝑡𝑇_𝑡⁻¹𝑦_𝑘≤max_𝑦∈𝑈𝐷𝑇_𝑡⁻¹𝑦²×^𝑙_𝑘=1𝐷𝑇_𝑡𝑇_𝑡⁻¹𝑦_𝑘+1−𝐷𝑇_𝑡𝑇_𝑡⁻¹𝑦_𝑘. (26)
Since 𝐷𝑇_𝑡 is Lipschitz continuous, we have 𝐷𝑇_𝑡⁻¹(𝑢)−𝐷𝑇_𝑡⁻¹(𝑣)=max_𝑦∈𝑈𝐷𝑇_𝑡⁻¹𝑦²𝐿₁^𝑙_𝑘=1𝑇_𝑡⁻¹𝑦_𝑘+1−𝑇_𝑡⁻¹𝑦_𝑘, (27) where 𝐿₁ is the maximum of all Lipschitz constants of 𝐷𝑇_𝑡 for all |𝑡|<𝑡_𝑉. Then finally, using the mean value theorem and property (c) in Lemma 5, we obtain 𝐷𝑇_𝑡⁻¹(𝑢)−𝐷𝑇_𝑡⁻¹(𝑣)≤max_𝑦∈𝑈𝐷𝑇_𝑡⁻¹𝑦³𝐿₁^𝑙_𝑘=1𝑦_𝑘+1−𝑦_𝑘≤max_𝑦∈𝑈𝐷𝑇_𝑡⁻¹𝑦³𝐿₁𝑐(𝑈)|𝑢−𝑣|. (28) Hence 𝐷𝑇_𝑡⁻¹ is Lipschitz continuous which shows that 𝑇_𝑡∶𝑈→𝑈 is a 𝐶^1,1 diffeomorphism for sufficiently small |𝑡|. Restricting to Ω, this proves that 𝑇_𝑡∶Ω→Ω_𝑡 is a 𝐶^1,1 diffeomorphism. (2) is clear because the fixed boundary is invariant under 𝑇_𝑡; that is, Γ_𝑡∶=𝑇_𝑡(Γ)=Γ since 𝐕 vanishes on Γ. Lastly, using Theorem 3, definition of 𝑇_𝑡, (1), and (2), we obtain (3).

Corollary 8. Let Ω and 𝑈 be two domains of ℝ² with 𝐶^1,1 boundary. Then for |𝑡|<𝑡_𝑉, where 𝑡_𝑉 is given by (20), the perturbed domain Ω_𝑡∶=𝑇_𝑡(Ω) is also of class 𝐶^1,1.

Proof. Given 𝑦∈𝜕Ω_𝑡, we let 𝑥=𝑇_𝑡⁻¹(𝑦)∈𝜕Ω. Then there exists a neighborhood 𝑁_𝑥 of 𝑥 and a diffeomorphism 𝑔_𝑥∈𝐶^1,1(𝑁_𝑥,𝐵(0,1)) such that 𝑔_𝑥(𝑁_𝑥∩Ω)=𝐵⁺(0,1),  𝑔_𝑥(𝑁_𝑥∩𝜕Ω)=𝐵⁰(0,1), and 𝑔_𝑥(𝑁_𝑥∩Ω^𝑐)=𝐵⁻(0,1). We have also shown that 𝑇_𝑡 defined in Theorem 7 is a 𝐶^1,1 diffeomorphism. Since 𝑇_𝑡⁻¹ is continuous, 𝑇_𝑡(𝑁_𝑥)=(𝑇_𝑡⁻¹)⁻¹(𝑁_𝑥)=𝑁_𝑦 is a neighborhood of 𝑦 in 𝑈. Define 𝑔_𝑦∶=𝑔_𝑥∘𝑇_𝑡⁻¹. This is bijective because 𝑔_𝑥 and 𝑇_𝑡⁻¹ are bijective. 𝑔_𝑦∈𝐶^1,1(𝑁_𝑦,𝐵(0,1)) because 𝑇_𝑡⁻¹∈𝐶^1,1(𝑈,𝑈) (hence 𝑇_𝑡⁻¹∈𝐶^1,1(𝑁_𝑦,𝑁_𝑥)) and 𝑔_𝑥∈𝐶^1,1(𝑁_𝑥,𝐵(0,1)). Also, 𝑔_𝑦⁻¹=𝑇_𝑡∘𝑔_𝑥⁻¹∈𝐶^1,1(𝐵(0,1),𝑁_𝑦).
Next, we note that 𝜕Ω_𝑡∩𝑁_𝑦=𝑇_𝑡(𝜕Ω)∩𝑁_𝑦=𝑇_𝑡(𝜕Ω)∩𝑇_𝑡(𝑁_𝑥). Since 𝑇_𝑡 is injective, we have 𝜕Ω_𝑡∩𝑁_𝑦=𝑇_𝑡(𝜕Ω∩𝑁_𝑥). Thus by definition of 𝑔_𝑦 we get 𝑔_𝑦𝜕Ω_𝑡∩𝑁_𝑦=𝑔_𝑦𝑇_𝑡𝜕Ω∩𝑁_𝑥=𝑔_𝑥𝜕Ω∩𝑁_𝑥=𝐵⁰(0,1). (29) We also observe the following: 𝑔_𝑦Ω_𝑡∩𝑁_𝑦=𝑔_𝑦𝑇_𝑡Ω∩𝑁_𝑥=𝑔_𝑥Ω∩𝑁_𝑥=𝐵⁺(),𝑔0,1_𝑦Ω^𝑐_𝑡∩𝑁_𝑦=𝑔_𝑦𝑇_𝑡Ω^𝑐∩𝑁_𝑥=𝑔_𝑥Ω^𝑐∩𝑁_𝑥=𝐵⁻(0,1). (30)
This shows that Ω_𝑡 is indeed of class 𝐶^1,1.

Remark 9. Theorem 7 and Corollary 8 tell us that the reference Ω and the perturbed domain Ω_𝑡 have the same topological structure and regularity under the perturbation of identity operator 𝑇_𝑡 for sufficiently small 𝑡. See Figure 3 for illustration.

Figure 3 

The action of 𝑇𝑡 on a 𝐶1,1-domain.

Properties of 𝑇_𝑡. In addition to (16) we also use the following notations throughout the work: 𝑀_𝑡(𝑥)=𝐷𝑇_𝑡(𝑥)^−𝑇,𝑥∈𝐴𝑈,_𝑡(𝑥)=𝐼_𝑡(𝑥)𝑀^𝑇_𝑡(𝑥)𝑀_𝑡(𝑥),𝑥∈𝑤𝑈,_𝑡(𝑥)=𝐼_𝑡(𝑥)𝐷𝑇_𝑡(𝑥)^−𝑇𝐧(𝑥),𝑥∈Σ. (31)

Remark 10. We note the following observations for fixed, sufficiently small 𝑡. (1)𝐼_𝑡∈𝐶^0,1(𝑈).(2)𝑀_𝑡,𝑀^𝑇_𝑡∈𝐶(𝑈;ℝ^2×2).(3)𝐴_𝑡∈𝐶(𝑈;ℝ^2×2).(4)𝑤_𝑡∈𝐶(Σ;ℝ).(5)𝐕∈𝐶^1,1(𝑈;ℝ²) implies that |𝐕|_∞ and |𝐷𝐕|_∞ are both finite.
We now provide several properties of 𝑇_𝑡.

Lemma 11 (see [9, 13, 16, 18]). Consider the transformation 𝑇_𝑡, where the fixed vector field 𝐕 belongs to Θ, defined in (13). Then there exists 𝑡_𝑉>0 such that 𝑇_𝑡 and the functions in (16) and (31) restricted to the interval 𝐼_𝑉=(−𝑡_𝑉,𝑡_𝑉) have the following regularity and properties.(1)𝑡↦𝑇_𝑡∈𝐶¹(𝐼_𝑉,𝐶^1,1(𝑈,ℝ²)).(2)𝑡↦𝐼_𝑡∈𝐶¹(𝐼_𝑉,𝐶^0,1(𝑈)).(3)𝑡↦𝑇_𝑡⁻¹∈𝐶(𝐼_𝑉,𝐶¹(𝑈,ℝ²)).(4)𝑡↦𝑤_𝑡∈𝐶¹(𝐼_𝑉,𝐶(Σ)).(5)𝑡↦𝐴_𝑡∈𝐶(𝐼_𝑉,𝐶(𝑈,ℝ^2×2)).(6)There is 𝛽>0 such that 𝐴_𝑡(𝑥)≥𝛽𝐼 for 𝑥∈𝑈.(7)(𝑑/𝑑𝑡)𝑇_𝑡|_𝑡=0=𝐕.(8)(𝑑/𝑑𝑡)𝑇_𝑡⁻¹|_𝑡=0=−𝐕.(9)(𝑑/𝑑𝑡)𝐷𝑇_𝑡|_𝑡=0=𝐷𝐕.(10)(𝑑/𝑑𝑡)(𝐷𝑇_𝑡)⁻¹|_𝑡=0=−𝐷𝐕.(11)(𝑑/𝑑𝑡)𝐼_𝑡|_𝑡=0=div𝐕.(12)(𝑑/𝑑𝑡)𝐴_𝑡|_𝑡=0=(div𝐕)𝐼−(𝐷𝐕+(𝐷𝐕)^𝑇)≡𝐴.(13)lim_𝑡→0𝑤_𝑡=1.(14)(𝑑/𝑑𝑡)𝑤_𝑡|_𝑡=0=div_Σ𝐕,where the surface divergence div_Σ is defined by div_Σ𝐕=div𝐕|_Σ−(𝐷𝐕𝐧)⋅𝐧. (32)

We provide proofs for properties (3) and (8). The rest can be seen in [19].

Proof. (3) Suppose 𝑥∈𝑈, 𝑡,ℎ∈𝐼_𝑉, and 𝑦=𝑇_ℎ⁻¹(𝑥). Then 𝑇_𝑡⁻¹(𝑥)−𝑇_ℎ⁻¹(𝑥)=𝑇_𝑡⁻¹𝑇_ℎ𝑦−𝑇_ℎ⁻¹𝑇_ℎ𝑦=𝑇_𝑡⁻¹𝑇_ℎ𝑦−𝑇_𝑡⁻¹𝑇_𝑡𝑦. (33) Using Lemma 5, we connect 𝑇_ℎ(𝑦) and 𝑇_𝑡(𝑦) by a chain 𝑦_𝑘, 𝑘=1,…,𝑙+1, satisfying (i)𝑦₁=𝑇_𝑡(𝑦), 𝑦_𝑘∈𝑈 for 2≤𝑘≤𝑙, 𝑦_𝑙+1=𝑇_ℎ(𝑦),(ii)(𝑦_𝑘,𝑦_𝑘+1)⊂𝑈 for 1≤𝑘≤𝑙,(iii)∑^𝑙_𝑘=1|𝑦_𝑘+1−𝑦_𝑘|≤𝑐(𝑈)|𝑇_ℎ(𝑦)−𝑇_𝑡(𝑦)|,and then we get 𝑇_𝑡⁻¹(𝑥)−𝑇_ℎ⁻¹(𝑥≤_𝑙𝑘=1)𝑇_𝑡⁻¹𝑦_𝑘−𝑇_𝑡⁻¹𝑦_𝑘+1≤max_𝑢∈𝑈𝐷𝑇_𝑡⁻¹(𝑢)^𝑙_𝑘=1𝑦_𝑘−𝑦_𝑘+1≤𝑐(𝑈)max_𝑢∈𝑈𝐷𝑇_𝑡⁻¹(𝑢)𝑇_ℎ𝑦−𝑇_𝑡𝑦,≤𝑐(𝑈)max_𝑢∈𝑈𝐷𝑇_𝑡⁻¹(𝑢)𝐕𝑦|𝑡−ℎ|. (34) Thus, max_𝑥∈𝑈𝑇_𝑡⁻¹(𝑥)−𝑇_ℎ⁻¹(𝑥)≤𝑐(𝑈)max_𝑢∈𝑈𝐷𝑇_𝑡⁻¹(𝑢)|𝐕|_∞|𝑡−ℎ|. (35) By reducing 𝑡_𝑉 if necessary, we can assume without loss of generality that |𝑡𝐷𝐕|_∞<1 for 𝑡∈(−𝑡_𝑉,𝑡_𝑉). This allows us to represent (𝐷𝑇_𝑡)⁻¹ as a Neumann series: 𝐷𝑇_𝑡⁻¹(𝑥)=𝐷𝑇_𝑡⁻¹∘𝑇_𝑡⁻¹(𝑥)=^∞_𝑘=0(−𝑡)^𝑘(𝐷𝐕)^𝑘𝑇_𝑡⁻¹(𝑥), (36) and its norm is estimated as follows: 𝐷𝑇_𝑡⁻¹(𝑥)_∞≤^∞_𝑘=0(−𝑡)^𝑘(𝐷𝐕)^𝑘𝑇_𝑡⁻¹(𝑥)_∞=11−|𝑡||𝐷𝐕|_∞<11−𝑡_𝑉|𝐷𝐕|_∞=∶𝑀. (37) This shows uniform convergence in 𝑡∈𝐼_𝑉 and 𝑥∈𝑈. Hence, for every 𝜖>0 one can choose a 𝛿∶=𝜖/𝑐(𝑈)𝑀|𝐕|_∞>0 which implies that, for every 𝑡,ℎ∈𝐼_𝑉, |𝑇_𝑡⁻¹−𝑇_ℎ⁻¹|_∞<𝜖 whenever |𝑡−ℎ|<𝛿. In other words, 𝑡↦𝑇_𝑡⁻¹∈𝐶(𝐼_𝑉,𝐶(𝑈,ℝ²)).
To show that 𝑡↦𝑇_𝑡⁻¹ is continuous from 𝐼_𝑉 to 𝐶¹(𝑈,ℝ²), we only need to show that for every 𝑡∈𝐼_𝑉, |𝐷𝑇_𝑡⁻¹−𝐷𝑇_ℎ⁻¹|_∞<𝜖 whenever |𝑡−ℎ|<𝛿 and ℎ∈𝐼_𝑉. Let 𝑧∈𝑈. Using (37), estimate |𝐷𝑇_𝑡⁻¹(𝑧)−𝐷𝑇_ℎ⁻¹(𝑧)| as follows: 𝐷𝑇_𝑡⁻¹(𝑧)−𝐷𝑇_ℎ⁻¹(𝑧=)𝐷𝑇_𝑡⁻¹𝑇_𝑡⁻¹(𝑧)−𝐷𝑇_ℎ⁻¹𝑇_ℎ⁻¹(𝑧)=𝐷𝑇_𝑡⁻¹𝑇_𝑡⁻¹(𝑧)×𝐷𝑇_ℎ𝑇_ℎ⁻¹(𝑧)−𝐷𝑇_𝑡𝑇_𝑡⁻¹(𝑧)×𝐷𝑇_ℎ⁻¹𝑇_ℎ⁻¹(𝑧)≤𝑀²𝐷𝑇_ℎ𝑇_ℎ⁻¹(𝑧)−𝐷𝑇_𝑡𝑇_𝑡⁻¹(𝑧). (38) Using the definition of Jacobian of a transformation and the regularity of 𝐕, we further simplify (38) as follows: 𝐷𝑇_𝑡⁻¹(𝑧)−𝐷𝑇_ℎ⁻¹(𝑧)≤𝑀²ℎ𝐷𝐕𝑇_ℎ⁻¹(𝑧)−𝑡𝐷𝐕𝑇_𝑡⁻¹(𝑧)≤𝑀²|ℎ−𝑡|𝐷𝐕𝑇_𝑠⁻¹(𝑧)+𝑀²|𝑡|𝐷𝐕𝑇_ℎ⁻¹(𝑧)−𝐷𝐕𝑇_𝑡⁻¹(𝑧)≤𝑀²𝐾|ℎ−𝑡|+𝑀²𝑡_𝑉𝐿𝑇_ℎ⁻¹(𝑧)−𝑇_𝑡⁻¹(𝑧), (39) where 𝐿 is the Lipshitz constant for 𝐷𝐕 and 𝐾 is upper bound for |𝐷𝐕|_∞. Taking the maximum of both sides of the inequality for all 𝑧∈𝑈 and using (34) we get 𝐷𝑇_𝑡⁻¹−𝐷𝑇_ℎ⁻¹_∞≤𝛼|ℎ−𝑡|, (40) where 𝛼=𝑀²𝐾+𝑀²𝑡_𝑉𝐿𝑐(𝑈)𝑀|𝐕|_∞. Thus, for any 𝜖>0, we choose 𝛿=𝜖/𝛼, so that if 0<|ℎ−𝑡|<𝛿, then |𝐷𝑇_𝑡⁻¹−𝐷𝑇_ℎ⁻¹|_∞<𝜖. Therefore, 𝑡↦𝑇_𝑡⁻¹∈𝐶(𝐼_𝑉,𝐶¹(𝑈,ℝ²)).
Proof of property (8) in Lemma 11 is as follows. Given 𝑥∈𝑈, we have 𝑇_𝑡(𝑇_𝑡⁻¹(𝑥))=𝑥. This implies that 𝑇_𝑡+ℎ𝑇⁻¹_𝑡+ℎ(𝑥)−𝑇_𝑡𝑇_𝑡⁻¹(𝑥)=0. (41) Manipulating the left hand side of (41), we get 𝑇_𝑡+ℎ𝑇⁻¹_𝑡+ℎ(𝑥)−𝑇_𝑡𝑇⁻¹_𝑡+ℎ(𝑥+𝑇_𝑡)𝑇⁻¹_𝑡+ℎ(𝑥)−𝑇_𝑡𝑇_𝑡⁻¹(𝑥)=∶𝐴+𝐵. (42) We first work on 𝐴. Applying the definition of 𝑇_𝑡, we get 𝐴=𝑇⁻¹_𝑡+ℎ(𝑥)+(𝑡+ℎ)𝐕𝑇⁻¹_𝑡+ℎ(𝑥)−𝑇⁻¹_𝑡+ℎ(𝑥)−𝑡𝐕𝑇⁻¹_𝑡+ℎ(𝑥)=ℎ𝐕𝑇⁻¹_𝑡+ℎ(𝑥). (43) Thus, lim_ℎ→01ℎ𝐴=𝐕𝑇_𝑡⁻¹(𝑥). (44) Similarly, we can write 𝐵 as follows: 𝐵=𝑇⁻¹_𝑡+ℎ(𝑥)−𝑇_𝑡⁻¹(𝑥)+𝑡𝐕𝑇⁻¹_𝑡+ℎ(𝑥)−𝑡𝐕𝑇_𝑡⁻¹(𝑥). (45) Hence, we have 1ℎ1𝐵=ℎ𝑇⁻¹_𝑡+ℎ(𝑥)−𝑇_𝑡⁻¹(𝑥1)+𝑡ℎ𝐕𝑇⁻¹_𝑡+ℎ(𝑥)−𝐕𝑇_𝑡⁻¹(𝑥). (46) Suppose 𝑣 is a coordinate function of 𝐕. By the mean value theorem, we observe that 1ℎ𝑣𝑇⁻¹_𝑡+ℎ(𝑥)−𝑣𝑇_𝑡⁻¹(𝑥(𝑐)1)=𝐷𝑣ℎ𝑇⁻¹_𝑡+ℎ(𝑥)−𝑇_𝑡⁻¹(𝑥), (47) where 𝑐 is a point on the segment joining 𝑇_𝑡⁻¹(𝑥) and 𝑇⁻¹_𝑡+ℎ(𝑥), and as ℎ tends to infinity, (47) tends to (𝐷𝑣)(𝑑/𝑑𝑡)𝑇_𝑡⁻¹(𝑥). Thus, lim_ℎ→01ℎ𝑑𝐵=𝑇𝑑𝑡^𝑡−1(𝑥)+𝑡𝐷𝐕𝑇_𝑡⁻¹(𝑥)𝑑𝑇𝑑𝑡^𝑡−1(𝑥). (48) Combining (44) and (48), we get 𝐕𝑇_𝑡⁻¹(𝑥)+𝐼+𝑡𝐷𝐕𝑇_𝑡⁻¹(𝑥)𝑑𝑇𝑑𝑡^𝑡−1(𝑥)=0, (49) which implies that 𝑑𝑇𝑑𝑡^𝑡−1(𝑥)=−𝐼+𝑡𝐷𝐕𝑇_𝑡⁻¹(𝑥)⁻¹𝐕𝑇_𝑡⁻¹(𝑥). (50) Evaluating (50) at 𝑡=0, we get (𝑑/𝑑𝑡)𝑇_𝑡⁻¹|_𝑡=0=−𝐕.

3.3. The Method of Mapping

If 𝑢 is defined in Ω and 𝑢_𝑡 is defined in Ω_𝑡, then the direct comparison of 𝑢_𝑡 with 𝑢 is generally not possible since the functions are defined on different domains. To overcome this difficulty, one maps 𝑢_𝑡 back to Ω by composing it with 𝑇_𝑡; that is, one defines 𝑢_𝑡∘𝑇_𝑡∶Ω→ℝ. With this new mapping one can define the material and the shape derivatives of states, the domain and boundary integral transformations, and derivatives of integrals, as well as the Eulerian derivative of the shape functional. This technique is called the method of mapping.

Material and Shape Derivatives. The material and shape derivatives of state variables are defined as follows [20, 21].

Definition 12. Let 𝑢 be defined in [0,𝑡_𝑉]×𝑈. An element 𝑢̇∈𝐻^𝑘(Ω), called the material derivative of 𝑢, is defined as 𝑢̇(𝑥)∶=𝑢̇(0,𝑥)∶=lim_𝑡→0+𝑢𝑡,𝑇_𝑡(𝑥)−𝑢(0,𝑥)𝑡=𝑑𝑑𝑡𝑢(𝑡,𝑥+𝑡𝐕(𝑥))_𝑡=0 (51)

if the limit exists in (𝐻^𝑘(Ω)).

Remark 13. The material derivative can be written as 𝑢̇(𝑥)=lim_𝑡→0+𝑢_𝑡∘𝑇_𝑡(𝑥)−𝑢(𝑥)𝑡=𝑑𝑑𝑡(𝑢_𝑡∘𝑇_𝑡(𝑥))_𝑡=0. (52) It characterizes the behavior of the function 𝑢 at 𝑥∈Ω⊂𝑈 in the direction 𝐕(𝑥).

Definition 14. Let 𝑢 be defined in [0,𝑡_𝑉]×𝑈. An element 𝑢^′∈𝐻^𝑘(Ω) is called the shape derivative of 𝑢 at Ω in the direction 𝐕, if the following limit exists in 𝐻^𝑘(Ω): 𝑢^′(𝑥)∶=𝑢^′(0,𝑥)∶=lim_𝑡→0+𝑢(𝑡,𝑥)−𝑢(0,𝑥)𝑡. (53)

Remark 15. The shape derivative of 𝑢 is also defined as follows: 𝑢^′𝜕(𝑥)∶=𝑢𝜕𝑡^𝑡(𝑥)_𝑡=0. (54) We note that if 𝑢̇ and ∇𝑢⋅𝑉 exist in 𝐻^𝑘(Ω), then the shape derivative can be written as 𝑢^′(𝑥)=𝑢̇(𝑥)−(∇𝑢⋅𝐕)(𝑥). (55) In general, if 𝑢̇(𝑥) and ∇𝑢⋅𝐕(𝑥) both exist in 𝑊^𝑚,𝑝(Ω), then 𝑢^′(𝑥) also exists in that space.
Domain and Boundary Transformations

Lemma 16 (see [18]). (1) Let 𝜑_𝑡∈𝐿¹(Ω_𝑡). Then 𝜑_𝑡∘𝑇_𝑡∈𝐿¹(Ω) and _Ω𝑡𝜑_𝑡𝑑𝑥_𝑡=_Ω𝜑_𝑡∘𝑇_𝑡𝐼_𝑡𝑑𝑥. (56)
(2) Let 𝜑_𝑡∈𝐿¹(𝜕Ω_𝑡). Then 𝜑_𝑡∘𝑇_𝑡∈𝐿¹(𝜕Ω) and _𝜕Ω𝑡𝜑_𝑡𝑑𝑠_𝑡=_𝜕Ω𝜑_𝑡∘𝑇_𝑡𝑤_𝑡𝑑𝑠, (57) where 𝐼_𝑡 and 𝑤_𝑡 are defined in (31).

Proofs can be found in [13, 18].

Domain and Boundary Differentiation. We recall some results concerning the derivative of integrals with respect to the domain of integration. For the first theorem, it is sufficient to have 𝐶^0,1 domains while the second theorem requires 𝐶^1,1 domains. For proofs, see [18].

Theorem 17 (domain differentiation formula). Let 𝑢∈𝐶(𝐼_𝑉,𝑊^1,1(𝑈)) and suppose 𝑢̇(0,⋅)∶=(𝑑/𝑑𝑡)𝑢(𝑡,𝑇_𝑡(⋅))|_𝑡=0 exists in 𝐿¹(𝑈). Then 𝑑𝑑𝑡_Ω𝑡𝑢(𝑡,𝑥)𝑑𝑥_𝑡=0=_Ω𝑢^′(0,𝑥)𝑑𝑥+_Σ𝑢(0,𝑠)𝐕⋅𝐧𝑑𝑠. (58)

Theorem 18 (boundary differentiation formula). Let 𝑢 be defined in a neighborhood of Γ. If 𝑢∈𝐶(𝐼_𝑉,𝑊^2,1(𝑈)) and 𝑢̇(0,⋅)∈𝑊^1,1(𝑈), then 𝑑𝑑𝑡_Σ𝑡𝑢(𝑡,𝑠)𝑑𝑠_𝑡=0=_Σ𝑢^′(0,𝑠)𝑑𝑠+_Σ𝜕𝑢()𝜅𝜕𝐧+𝑢0,𝑠𝐕⋅𝐧𝑑𝑠, (59) where 𝜅 is the mean curvature of the free boundary Σ.
The First-Order Eulerian Derivative

Definition 19. The Eulerian derivative of the shape functional 𝐽∶Ω→ℝ defined in (7) at the domain Ω in the direction of the deformation field 𝐕∈Θ is given by 𝑑𝐽(Ω;𝐕)∶=lim_𝑡→0+𝐽Ω_𝑡−𝐽(Ω)𝑡, (60) if the limit exists.

Remark 20. 𝐽 is said to be shape differentiable at Ω if 𝑑𝐽(Ω;𝐕) exists for all 𝐕∈Θ and is linear and continuous with respect to 𝐕.

4. Main Result

In this section we derive in a rigorous manner the first-order shape derivative of the Kohn-Vogelius functional 𝐽, defined by (7), subject to the Dirichlet and Neumann boundary value problems (BVPs) (4) and (6), respectively. Our strategy bypasses the material or shape derivatives of states. In the derivation, we have employed techniques used in [9, 10, 13] but there is no need to use adjoint variables.

This section discusses the variational forms of the PDEs, the state variables in the perturbed domains, the Hölder continuity of the state variables, and the higher regularity of the solutions to the BVPs. The rest of the proof is presented in the last part of this section.

4.1. Variational Forms of the Dirichlet and Neumann Problems

We recall that we are considering the shape optimization problem (7) where 𝑢_𝐷 solves the pure Dirichlet problem (4) and 𝑢_𝑁 solves the Neumann problem (6). As in [13], we consider the Hilbert space 𝐻¹_Γ,0=𝜑∈𝐻¹(Ω)∶𝜑|_Γ=0, (61) which is endowed with the norm 𝜑²₁=_Ω∇𝜑²d𝑥, (62) and a linear manifold defined by 𝐻¹_Γ,𝑣=𝜑∈𝐻¹(Ω)∶𝜑|_Γ=𝑣 (63) for 𝑣∈𝐻^1/2(Γ).

First, we determine the variational equations for the Dirichlet and the Neumann problems. The variational form of the Dirichlet problem (4) is given by the following.

Find 𝑢_𝐷∈𝐻¹(Ω) such that ∇𝑢_𝐷,∇𝜓_Ω=0∀𝜓∈𝐻¹₀(Ω)𝑢_𝐷𝑢=1onΓ,_𝐷=0onΣ. (64) Equation (64) can be shown to have a unique solution using Theorem 2.4.2.5 of [22]. Similarly, the variational form of the Neumann problem (6) is formulated as follows.

Find 𝑢_𝑁∈𝐻¹(Ω) such that ∇𝑢_𝑁,∇𝜑_Ω−𝜆,𝜑_Σ=0∀𝜑∈𝐻¹_Γ,0(Ω)𝑢_𝑁=1onΓ. (65) It is also well known that (65) has a unique solution.

4.2. Analysis of State Variables in Deformed Domains

We now consider the class of perturbed problems: min_Ω𝑡𝐽Ω_𝑡=12_Ω𝑡∇𝑢_𝐷,𝑡−𝑢_𝑁,𝑡²𝑑𝑥, (66) where 𝑢_𝐷,𝑡 solves the pure Dirichlet problem −Δ𝑢_𝐷,𝑡=0inΩ_𝑡,𝑢_𝐷,𝑡𝑢=1onΓ,_𝐷,𝑡=0onΣ_𝑡, (67) and 𝑢_𝑁,𝑡 solves the Neumann problem −Δ𝑢_𝑁,𝑡=0inΩ_𝑡,𝑢_𝑁,𝑡=1onΓ,𝜕𝑢_𝑁,𝑡𝜕𝐧_𝑡=𝜆onΣ_𝑡. (68) Here, 𝐧_𝑡 is the outward unit normal to the deformed free boundary Σ_𝑡. The variational form of (67) is formulated as follows.

Find 𝑢_𝐷,𝑡∈𝐻¹(Ω_𝑡) such that ∇𝑢_𝐷,𝑡,∇𝜓_𝑡Ω𝑡=0∀𝜓_𝑡∈𝐻¹₀Ω_𝑡𝑢_𝐷,𝑡𝑢=1onΓ,_𝐷,𝑡=0onΣ_𝑡. (69) It is known that (69) has a unique solution.

Remark 21. The function 𝑢_𝑡∶Ω_𝑡→ℝ can be referred to as the reference domain by composing 𝑢_𝑡 with 𝑇_𝑡; that is, 𝑢^𝑡=𝑢_𝑡∘𝑇_𝑡∶Ω⟶ℝ (70) and by chain rule of differentiation, we get ∇𝑢_𝑡∘𝑇_𝑡=𝐷𝑇_𝑡^−𝑇∇𝑢^𝑡=∶𝑀_𝑡∇𝑢^𝑡. (71)

Let 𝑢_𝐷,𝑡 be the solution of (69). Applying Lemma 16 for all 𝜓_𝑡∈𝐻¹₀(Ω_𝑡) we have ∇𝑢_𝐷,𝑡,∇𝜓_𝑡Ω𝑡=_Ω𝑡∇𝑢_𝐷,𝑡⋅∇𝜓_𝑡=_Ω𝐼_𝑡∇𝑢_𝐷,𝑡∘𝑇_𝑡⋅∇𝜓_𝑡∘𝑇_𝑡. (72) Applying (71) and noting that 𝜓^𝑡∈𝐻¹₀(Ω) because 𝜓_𝑡∈𝐻¹₀(Ω_𝑡), we obtain ∇𝑢_𝐷,𝑡,∇𝜓_𝑡Ω𝑡=_Ω𝐼_𝑡𝑀_𝑡∇𝑢^𝑡_𝐷⋅𝑀_𝑡∇𝜓^𝑡=_Ω𝐼_𝑡𝑀^𝑇_𝑡𝑀_𝑡∇𝑢^𝑡_𝐷⋅∇𝜓^𝑡=𝐴_𝑡∇𝑢^𝑡_𝐷,∇𝑣_Ω, (73) where 𝑣=𝜓^𝑡∈𝐻¹₀(Ω). Hence, if 𝑢_𝐷,𝑡 solves the variational equation (69), then 𝑢^𝑡_𝐷=𝑢_𝐷,𝑡∘𝑇_𝑡 satisfies the variational equation, 𝐴_𝑡∇𝑢^𝑡_𝐷,∇𝑣_Ω=0 (74) for all 𝑣∈𝐻¹₀(Ω),  𝑢^𝑡_𝐷=1 on Γ, and 𝑢^𝑡_𝐷=0 on Σ.

Now we show that 𝑢^𝑡_𝐷 is the unique solution of (74) in 𝐻¹(Ω). First, we show that 𝑦^𝑡=𝑢^𝑡_𝐷−𝑢_𝐷∈𝐻¹₀(Ω) is the unique solution to 𝐴_𝑡∇𝑦^𝑡,∇𝑣_Ω=−𝐴_𝑡∇𝑢_𝐷,∇𝑣_Ω (75) for all 𝑣∈𝐻¹₀(Ω). The bilinear form 𝑏_𝑡(⋅,⋅)∶𝐻¹₀(Ω)×𝐻¹₀(Ω)→ℝ defined by 𝑏_𝑡𝑦^𝑡,𝑣=_Ω𝐴_𝑡∇𝑦^𝑡⋅∇𝑣,∀𝑦^𝑡,𝑣∈𝐻¹₀(Ω) (76) is continuous, because 𝑏_𝑡𝑦^𝑡,𝑣=_Ω𝐴_𝑡∇𝑦^𝑡≤𝐴⋅∇𝑣_𝑡∇𝑦^𝑡_𝐿2(Ω)|∇𝑣|_𝐿²_(Ω)≤2𝐴_𝑡∞∇𝑦^𝑡_𝐿2(Ω)|∇𝑣|_𝐿²_(Ω)=𝛽𝑦^𝑡_𝐻1(Ω)|𝑣|_𝐻¹_0(Ω). (77) The bilinear form is also coercive. To show this we recall that lim_𝑡→0𝐴_𝑡=𝐼 uniformly on Ω. This is equivalent to the statement ∀𝜖>0∃𝛿>0suchthat|𝑡|<𝛿⟹𝐴_𝑡−𝐼_∞<𝜖. (78) Let 𝜖=1/2. So for sufficiently small 𝑡, |𝐴_𝑡−𝐼|_∞<1/2, and 𝑏_𝑡𝑦^𝑡,𝑦^𝑡=_Ω𝐴_𝑡∇𝑦^𝑡⋅∇𝑦^𝑡=_Ω𝐴_𝑡−𝐼∇𝑦^𝑡⋅∇𝑦^𝑡+∇𝑦^𝑡2≥∇𝑦^𝑡2_𝐿2(Ω)−𝐴_𝑡−𝐼_∞∇𝑦^𝑡2_𝐿2(Ω)≥12𝑦^𝑡2_𝐻10(Ω). (79) So 𝑏_𝑡(𝑦^𝑡,𝑣) is coercive.

Next, we show that the functional 𝑓∶𝐻¹₀(Ω)→ℝ is bounded: 𝑓,𝑣_{𝐻−1×𝐻10}=−_Ω𝐴_𝑡∇𝑢_𝐷𝐴⋅∇𝑣≤2_𝑡∞∇𝑢_{𝐷𝐿2(Ω)}|∇𝑣|_𝐿²_(Ω)=2𝐴_𝑡∞𝑢_{𝐷𝐻10(Ω)}|𝑣|_𝐻¹_0(Ω)≤𝐶|𝑣|_𝐻¹_0(Ω). (80) Therefore, by the Lax-Milgram lemma, 𝑦^𝑡=𝑢^𝑡_𝐷−𝑢_𝐷 is the unique solution to the variational equation 𝐴_𝑡∇𝑦^𝑡,∇𝑣_Ω=−𝐴_𝑡∇𝑢_𝐷,∇𝑣_Ω,∀𝑣∈𝐻¹₀(Ω). (81) This implies the existence of a unique solution 𝑢^𝑡_𝐷 of (74) as verified below.

Let 𝑢^𝑡_𝐷=𝑦^𝑡+𝑢_𝐷. Using (81) we obtain 𝐴_𝑡∇𝑢^𝑡_𝐷,∇𝑣_Ω=𝐴_𝑡∇𝑦^𝑡,∇𝑣_Ω+𝐴_𝑡∇𝑢_𝐷,∇𝑣_Ω∀𝑣∈𝐻¹₀(Ω)=−𝐴_𝑡∇𝑢_𝐷,∇𝑣_Ω+𝐴_𝑡∇𝑢_𝐷,∇𝑣_Ω=0. (82) Thus (74) is satisfied. The boundary conditions are also satisfied because on Γ, 𝑦^𝑡=0 and 𝑢_𝐷=1 and on Σ, both 𝑦^𝑡 and 𝑢_𝐷 are zero. To show uniqueness, we let 𝑢^𝑡_𝐷 and 𝑢^𝑡_𝐷 be solutions of (74). This implies that there exist 𝑦^𝑡 and 𝑦^𝑡 such that 𝑢^𝑡_𝐷=𝑦^𝑡+𝑢_𝐷 and 𝑢^𝑡_𝐷=𝑦^𝑡+𝑢_𝐷, where 𝑦^𝑡 and 𝑦^𝑡 are solutions to (81). Taking the difference of 𝑢^𝑡_𝐷 and 𝑢^𝑡_𝐷 and considering that solution to (81) is unique, we get 𝑢^𝑡_𝐷=𝑢^𝑡_𝐷.

Next, we consider (68) whose variational form is formulated as follows.

Find 𝑢_𝑁,𝑡∈𝐻¹(Ω_𝑡) such that ∇𝑢_𝑁,𝑡,∇𝜑_𝑡Ω𝑡−𝜆,𝜑_𝑡Σ𝑡=0∀𝜑_𝑡∈𝐻¹_Γ,0Ω_𝑡,𝑢_𝑁,𝑡=1onΓ. (83) Similarly, if 𝑢_𝑁,𝑡 solves the variational problem (83), then 𝑢^𝑡_𝑁 solves the variational equation 𝐴_𝑡∇𝑢^𝑡_𝑁,∇𝜙_Ω−𝑤_𝑡𝜆,𝜙_Σ=0,∀𝜙∈𝐻¹_Γ,0(Ω), (84) where 𝑢^𝑡_𝑁=1 on Γ.

As shown before, the bilinear form 𝑏_𝑡(⋅,⋅)∶𝐻¹_Γ,0(Ω)×𝐻¹_Γ,0(Ω)→ℝ defined by 𝑏_𝑡𝑧^𝑡,𝑣=_Ω𝐴_𝑡∇𝑧^𝑡⋅∇𝑣,∀𝑧^𝑡,𝑣∈𝐻¹_Γ,0(Ω), (85) is coercive and continuous. The linear functional 𝑓∶𝐻¹_Γ,0(Ω)→ℝ defined by ⟨𝑓,𝜙⟩=∫_Σ𝜆𝑤_𝑡𝜙 is continuous on 𝐻¹_Γ,0(Ω) because 𝑓,𝜙_{𝐻−1×𝐻10}=_Σ𝜆𝑤_𝑡𝜙≤|𝜆|𝑤_𝑡∞Σ𝜙=|𝜆|𝑤_𝑡∞|Σ|^1/2𝛾𝜙_𝐿2(Σ)≤|𝜆|𝑤_𝑡∞|Σ|^1/2𝛾𝜙_𝐻1(Ω). (86) By the Lax-Milgram lemma, 𝑧^𝑡=𝑢^𝑡_𝑁−1∈𝐻¹_Γ,0(Ω) is the unique solution in 𝐻¹_Γ,0(Ω) of _Ω𝐴_𝑡∇𝑧^𝑡⋅∇𝜙=𝜆_Σ𝑤_𝑡𝜙,∀𝜙∈𝐻¹_Γ,0(Ω). (87) Let 𝑢^𝑡_𝑁=𝑧^𝑡+1. Then, by (87), we get 𝐴_𝑡∇𝑢^𝑡_𝑁,∇𝜙_Ω−𝑤_𝑡𝜆,𝜙_Σ=𝐴_𝑡∇𝑧^𝑡,∇𝜙_Ω−𝑤_𝑡𝜆,𝜙_Σ=0𝜙∈𝐻¹_Γ,0(Ω). (88) Since 𝑧^𝑡∈𝐻¹_Γ,0(Ω), 𝑢^𝑡_𝑁=𝑧^𝑡+1=1 on Γ. Uniqueness of 𝑢^𝑡_𝑁 follows from the uniqueness of 𝑧^𝑡. Therefore, 𝑢^𝑡_𝑁 is the unique solution of the variational problem (84) in 𝐻¹(Ω).

4.3. Hölder Continuity of the States

We show that 𝑢^𝑡_𝐷 and 𝑢^𝑡_𝑁 are Hölder continuous on 𝑡.

Theorem 22 (see [13]). The solutions 𝑢^𝑡_𝐷 of (74) are uniformly bounded in 𝐻¹(Ω) for 𝑡∈(−𝑡_𝑉,𝑡_𝑉) and lim_𝑡→01𝑡𝑢^𝑡_𝐷−𝑢_𝐷²_𝐻1(Ω)=0, (89) where 𝑢_𝐷 is the weak solution of (4).

Proof. We first prove the uniform boundedness of 𝑢^𝑡_𝐷 in 𝐻¹(Ω) for 𝑡∈(−𝑡_𝑉,𝑡_𝑉). Since 𝑢^𝑡_𝐷−𝑢_𝐷∈𝐻¹₀(Ω), by using coercivity of 𝐴_𝑡 we get 𝐴_𝑡∇𝑢^𝑡_𝐷−𝑢_𝐷,∇𝑢^𝑡_𝐷−𝑢_𝐷Ω≥12∇𝑢^𝑡_𝐷−𝑢_𝐷²_𝐻10(Ω). (90) Also, by applying (81), we have 𝐴_𝑡∇𝑢^𝑡_𝐷−𝑢_𝐷,∇𝑢^𝑡_𝐷−𝑢_𝐷Ω=−𝐴_𝑡∇𝑢_𝐷,∇𝑢^𝑡_𝐷−𝑢_𝐷Ω≤𝐴_𝑡∞𝑢_{𝐷𝐻1(Ω)}𝑢^𝑡_𝐷−𝑢_{𝐷𝐻10(Ω)}. (91) Therefore, 𝑢^𝑡_{𝐷𝐻1(Ω)}≤𝑢^𝑡_𝐷−𝑢_{𝐷𝐻10(Ω)}+𝑢_{𝐷𝐻1(Ω)}≤𝑡_𝑉+𝑢_{𝐷𝐻1(Ω)}. (92)
Now we take the difference between the weak form of (4) and the variational equation (74), to get 𝐴_𝑡∇𝑢^𝑡_𝐷,∇𝑣_Ω−∇𝑢_𝐷,∇𝑣_Ω=0,∀𝑣∈𝐻¹₀(Ω)𝐴_𝑡∇1𝑡(𝑢^𝑡_𝐷−𝑢_𝐷),∇𝑣_Ω=−(𝐴_𝑡−𝐼)𝑡∇𝑢_𝐷,∇𝑣_Ω. (93) Note that 𝑢_𝐷,𝑡∈𝐻¹_Γ,1(Ω_𝑡). This implies that 𝑢^𝑡_𝐷=𝑢_𝐷,𝑡∘𝑇_𝑡 is in 𝐻¹_Γ,1(Ω). Note also that 𝑢_𝐷∈𝐻¹_Γ,1(Ω). Let 𝑦^𝑡=(1/𝑡)(𝑢^𝑡_𝐷−𝑢_𝐷). So for sufficiently small 𝑡, 𝑦^𝑡∈𝐻¹₀(Ω). Now choosing 𝑣=𝑦^𝑡 as a test function and by the uniform coercivity of 𝐴_𝑡 one obtains 𝛼∇𝑦^𝑡2_𝐿2(Ω)≤𝐴_𝑡∇𝑦^𝑡,∇𝑦^𝑡_Ω1=−𝑡𝐴_𝑡−𝐼∇𝑢_𝐷,∇𝑦^𝑡_Ω≤1𝑡𝐴_𝑡−𝐼_∞∇𝑢_{𝐷𝐿2(Ω)}∇𝑦^𝑡_𝐿2(Ω). (94)
If ∇𝑦^𝑡=0, then the inequality above holds. For ∇𝑦^𝑡≠0 we have 𝛼∇𝑦^𝑡2_𝐿2(Ω)≤1𝑡𝐴_𝑡−𝐼∇𝑢_{𝐷𝐿2(Ω)}∇𝑦^𝑡_𝐿2(Ω)≤1𝑡𝐴_𝑡−𝐼_∞∇𝑢_{𝐷𝐿2(Ω)}∇𝑦^𝑡_𝐿2(Ω). (95) Hence, 1𝑡∇(𝑢^𝑡_𝐷−𝑢_𝐷)_𝐿2(Ω)≤1𝛼1𝑡𝐴_𝑡−𝐼_∞∇𝑢_{𝐷𝐿2(Ω)}. (96) Squaring and multiplying 𝑡 on both sides of the inequality give us 1𝑡∇𝑢^𝑡_𝐷−𝑢_𝐷²_𝐿2(Ω)1≤𝑡𝛼²1𝑡𝐴_𝑡−𝐼²_∞∇𝑢_𝐷²_𝐿2(Ω). (97) Since |(1/𝑡)(𝐴_𝑡−𝐼)|²_∞ is uniformly bounded in 𝑡 by Lemma 11, lim_𝑡→01𝑡𝑢^𝑡_𝐷−𝑢_𝐷²_𝐻1(Ω)=lim_𝑡→01𝑡∇𝑢^𝑡_𝐷−𝑢_𝐷²_𝐿2(Ω)=0 (98) follows.

Theorem 23. The solutions 𝑢^𝑡_𝑁 of (84) are uniformly bounded in 𝐻¹(Ω) for 𝑡∈(−𝑡_𝑉,𝑡_𝑉) and lim_𝑡→01𝑡𝑢^𝑡_𝑁−𝑢_𝑁²_𝐻1(Ω)=0, (99) where 𝑢_𝑁 is the solution of (65).

Proof. Subtracting (65) from (84) for all 𝜑∈𝐻¹_Γ,0(Ω) we get 0=𝐴_𝑡∇𝑢^𝑡_𝑁,∇𝜑_Ω−𝑤_𝑡𝜆,𝜑_Σ−∇𝑢_𝑁,∇𝜑_Ω+𝜆,𝜑_Σ=𝐴_𝑡∇𝑢^𝑡_𝑁−∇𝑢^𝑡_𝑁+∇𝑢^𝑡_𝑁−∇𝑢_𝑁,∇𝜑_Ω−𝑤_𝑡𝜆−𝜆,𝜑_Σ. (100) Hence ∇𝑢^𝑡_𝑁−𝑢_𝑁,∇𝜑_Ω=−𝐴_𝑡−𝐼∇𝑢^𝑡_𝑁,∇𝜑_Ω+𝑤_𝑡𝜆−𝜆,𝜑_Σ. (101) Note that 𝑢^𝑡_𝑁−𝑢_𝑁 belongs to 𝐻¹_Γ,0(Ω). Hence 𝑢^𝑡_𝑁−𝑢_𝑁²_𝐻1(Ω)=−𝐴_𝑡−𝐼∇𝑢^𝑡_𝑁,∇𝑢^𝑡_𝑁−𝑢_𝑁Ω+𝜆𝑤_𝑡−1,𝑢^𝑡_𝑁−𝑢_𝑁Σ. (102) By Cauchy-Schwarz inequality, we obtain 𝑢^𝑡_𝑁−𝑢_𝑁²_𝐻1(Ω)≤𝐴_𝑡−𝐼∇𝑢^𝑡_{𝑁𝐿2(Ω)}∇𝑢^𝑡_𝑁−𝑢_{𝑁𝐿2(Ω)}+|𝜆|𝑤_𝑡−1_∞𝑢^𝑡_𝑁−𝑢_{𝑁𝐿2(Σ)}. (103) Furthermore, by trace theorem we have |𝑢^𝑡_𝑁−𝑢_𝑁|_𝐿2(Σ)≤𝑐|𝑢^𝑡_𝑁−𝑢_𝑁|_𝐻1(Ω). Therefore, 𝑢^𝑡_𝑁−𝑢_𝑁²_𝐻1(Ω)≤𝐶𝐴_𝑡−𝐼∇𝑢^𝑡_{𝑁𝐿2(Ω)}+|𝜆|𝑤_𝑡−1_∞𝑢^𝑡_𝑁−𝑢_{𝑁𝐻1(Ω)} (104) holds, where 𝐶=max{1,𝑐}. This implies 𝑢^𝑡_𝑁−𝑢_{𝑁𝐻1(Ω)}≤𝐶𝐴_𝑡−𝐼∇𝑢^𝑡_{𝑁𝐿2(Ω)}+𝑤_𝑡𝜆−𝜆_𝐿2(Σ), (105) which entails 1√𝑡𝑢^𝑡_𝑁−𝑢_{𝑁𝐻1(Ω)}1≤𝐶√𝑡𝐴_𝑡−𝐼_𝐿2(Ω)∇𝑢^𝑡_{𝑁𝐿2(Ω)}+|𝜆||Σ|^1/21√𝑡𝑤_𝑡−1_∞=√1𝑡𝐶𝑡𝐴_𝑡−𝐼_𝐿2(Ω)∇𝑢^𝑡_{𝑁𝐿2(Ω)}+|𝜆||Σ|^1/21𝑡𝑤_𝑡−1_∞. (106)
We now show that ∇𝑢^𝑡_𝑁 is bounded uniformly in 𝑡 in 𝐿²(Ω). Since 𝑢^𝑡_𝑁−1∈𝐻¹_Γ,0(Ω), we have 𝛼∇𝑢^𝑡_𝑁²_𝐿2(Ω)=𝛼∇𝑢^𝑡_𝑁−1²_𝐿2(Ω)≤_Ω𝐴_𝑡∇𝑢^𝑡_𝑁−1⋅∇𝑢^𝑡_𝑁=−1_Ω𝐴_𝑡∇𝑢^𝑡_𝑁⋅∇𝑢^𝑡_𝑁=|𝜆|−1_Σ𝑤_𝑡𝑢^𝑡_𝑁|𝜆𝑤−1≤𝑐|_𝑡𝐿∞|Σ|^1/2𝛾∇𝑢^𝑡_{𝑁𝐿2(Ω)}. (107) Consequently, ∇𝑢^𝑡_{𝑁𝐿2(Ω)}≤|𝜆|𝛼𝑐𝑤_𝑡𝐿∞|Σ|^1/2𝛾, (108) and this shows that ∇𝑢^𝑡_𝑁 is uniformly bounded in 𝐿²(Ω) because |𝑤_𝑡|_∞ is bounded. In addition, 𝐴_𝑡 and 𝑤_𝑡 are differentiable at 𝑡=0 by Lemma 11. Therefore, lim_𝑡→01𝑡𝑢^𝑡_𝑁−𝑢_𝑁²_𝐻1(Ω)=0. (109)

4.4. Higher Regularity of the Solutions

In this section we will show that the solutions to the PDEs (4) and (6) have higher regularity. We begin by considering the state variable 𝑢_𝐷,𝑢_𝑁∈𝐻¹(Ω). For 𝐶^1,1 domains, we show that these solutions also exist in 𝐻²(Ω) and more generally in 𝐻^𝑘+2(Ω) if domains are of class 𝐶^𝑘+1,1, 𝑘≥0.

To prove higher regularity of 𝑢_𝐷, we require the following two theorems, which are proven in [22].

Theorem 24 (see [22, page 124]). Let Ω be a bounded open subset of ℝ^𝑛 with a 𝐶^1,1 boundary. Consider the Dirichlet boundary value problem: 𝐴𝑢=𝑓𝑖𝑛Ω,𝑢=𝑔𝑜𝑛𝜕Ω, (110) where 𝐴=^𝑛_𝑖,𝑗=1𝐷_𝑖𝑎_𝑖𝑗𝐷_𝑗𝑢+^𝑛_𝑖=1𝑎_𝑖𝐷_𝑖𝑢+𝑎₀𝑢. (111) Let 𝑎_𝑖,𝑗 be uniformly Lipschitz functions and let 𝑎_𝑖 be bounded measurable functions such that 𝑎_𝑗,𝑖=𝑎_𝑖,𝑗, 1≤𝑖, 𝑗≤𝑛 and that there exists 𝛼>0 with ^𝑛_𝑖,𝑗=1𝑎_𝑖,𝑗(𝑥)𝜁_𝑖𝜁_𝑗≤−𝛼𝜁² (112) for all 𝜁∈ℝ^𝑛 and for almost every 𝑥∈Ω. Assume in addition that either (i)𝑎_𝑖=0,   1≤𝑖≤𝑛 and 𝑎₀≥0 a.e. or(ii)𝑎₀≥𝛽>0 a.e.Then for every 𝑓∈𝐿^𝑝(Ω) and every 𝑔∈𝑊^{2−1/𝑝,𝑝}(𝜕Ω), there exists a unique 𝑢∈𝑊^2,𝑝(Ω) that solves (110).

Theorem 25 (see [22, page 128]). Let Ω be a bounded open subset of ℝ^𝑛 with a 𝐶^𝑘+1,1 boundary. Consider the operator 𝐴 defined by (110) with 𝑎_𝑖𝑗=𝑎_𝑗𝑖∈𝐶^𝑘,1(Ω), and assume that there exists 𝛼>0 such that (112) holds for all 𝜁∈ℝ^𝑛 and for every 𝑥∈Ω. Also, consider a real boundary operator 𝐵 which is either the identity operator (𝑜𝑟𝑑𝑒𝑟𝑑=0) or 𝐵𝑢=∑^𝑛_𝑖=1𝑏_𝑖𝐷_𝑖𝑢 with 𝑏_𝑖∈𝐶^𝑘,1(Ω), 1≤𝑖≤𝑛(𝑜𝑟𝑑𝑒𝑟𝑑=1) and ∑^𝑛_𝑖=1𝑏_𝑖𝐧^𝑖≠0(𝐧^𝑖𝑖𝑠𝑡ℎ𝑒𝑖𝑡ℎ𝑐𝑜𝑚𝑝𝑜𝑛𝑒𝑛𝑡𝑜𝑓𝐧) everywhere on Γ∶=𝜕Ω. Furthermore, assume that 𝑢∈𝑊^2,𝑝(Ω) satisfies 𝐴𝑢=𝑓∈𝑊^𝑘,𝑝(Ω) and 𝛾𝐵𝑢=𝑔∈𝑊^{2+𝑘−𝑑−1/𝑝,𝑝}(Γ). Then 𝑢∈𝑊^𝑘+2,𝑝(Ω).

We will also justify the higher regularity of 𝑢_𝑁. We use the following results whose proofs are given in the corresponding texts.

Theorem 26 (see [23, page 316]). Let Ω be a bounded open subset of ℝ^𝑛. Suppose 𝑢∈𝐻¹(Ω) is a weak solution of the PDE 𝐿𝑢=𝑓𝑖𝑛Ω, (113) where 𝐿𝑢=−^𝑛_𝑖,𝑗=1𝑎^𝑖𝑗(𝑥)𝑢_{𝑥𝑖𝑥𝑗}+^𝑛_𝑖=1𝑏^𝑖(𝑥)𝑢_𝑥𝑖+𝑐(𝑥)𝑢. (114) Assume furthermore that 𝑎^𝑖𝑗,𝑏^𝑖,𝑐∈𝐶^∞(Ω) for 𝑖,𝑗=1,2,…,𝑛, and 𝑓∈𝐶^∞(Ω). Then 𝑢∈𝐶^∞(Ω).

Theorem 27 (see [24, page 12]). Let Ω⊂ℝ^𝑛 be a bounded domain with 𝐶^𝑘+1,1 boundary 𝜕Ω for some nonnegative integer 𝑘. Suppose the data 𝑓 and 𝑔 of the problem −Δ𝑢=𝑓𝑖𝑛Ω,𝑢=𝑔𝑜𝑛𝜕Ω (115) are in 𝑊^𝑘,𝑝(Ω) and 𝑊^{𝑘+2−1/𝑝,𝑝}(𝜕Ω), respectively, for some real number 𝑝 with 1<𝑝<∞. Then 𝑢∈𝑊^𝑘+2,𝑝(Ω).

For proof, see [22].

Theorem 28 (see [22, page 84]). Let Ω⊂ℝ^𝑛 be a bounded domain with 𝐶^1,1 boundary Γ and 1<𝑝<∞. Consider the Neumann problem −Δ𝑢+𝑎₀𝛾𝑢=𝑓𝑖𝑛Ω,𝜕𝑢𝜕𝐧=𝑔𝑜𝑛Γ. (116) If 0<𝑎₀∈𝐿^∞(Ω), 𝑓∈𝐿^𝑝(Ω), and 𝑔∈𝑊^{1−1/𝑝,𝑝}(𝜕Ω), then the weak solution 𝑢 to (116) exists in 𝑊^2,𝑝(Ω).

For proof, see [22].

Using the theorems presented above, we will now prove our claim that the solutions to the PDEs (4) and (6) have indeed higher regularity. This result is given in the following theorem.

Theorem 29. Let Ω be a bounded domain with boundary of class 𝐶^1,1. Let 𝑢_𝐷,𝑢_𝑁∈𝐻¹(Ω) be weak solutions of the 𝐵𝑉𝑃𝑠 (4) and (6), respectively. Then 𝑢_𝐷 and 𝑢_𝑁 also belong to 𝐻²(Ω). More generally, if Ω is of class 𝐶^𝑘+1,1, where 𝑘 is a nonnegative integer, then 𝑢_𝐷 and 𝑢_𝑁 are elements of 𝐻^𝑘+2(Ω).

Proof. We first consider the solution 𝑢_𝐷∈𝐻¹(Ω) to the Dirichlet problem (4). We use Theorem 25 to show that 𝑢_𝐷 is an element of 𝐻²(Ω). Here, (110) is applied with the following settings.
We consider 𝑛=2. The domain Ω is of class 𝐶^1,1. 𝐿=−Δ, and hence 𝑎_𝑖𝑗=𝑎_𝑗𝑖=−1 for 𝑖=𝑗 and 𝑎_𝑖𝑗=𝑎_𝑗𝑖=0 for 𝑖≠𝑗, with 𝑖,𝑗=1,2. We also observe that ∑²_𝑖,𝑗=1𝑎_𝑖𝑗(𝑥)𝜁_𝑖𝜁_𝑗=−(𝜁₁)²−(𝜁₂)²=−1|𝜁|², for all 𝜁=(𝜁₁,𝜁₂)^𝑇. Thus 𝛼=1. Furthermore, we have the following data: 𝑓=0∈𝐿²(Ω),  𝑔=1∈𝐻^3/2(Γ), 𝑔=0∈𝐻^3/2(Σ). Therefore, by using Theorem 24, there exists a unique 𝑢=𝑢_𝐷∈𝐻²(Ω), which is a solution to (4).
For higher regularity of 𝑢_𝐷 we apply Theorem 25. At first we consider 𝐶^2,1-domains. In this case, 𝑘=1. We have 𝑎_𝑖𝑗=𝑎_𝑗𝑖=−1 for 𝑖=𝑗 and 𝑎_𝑖𝑗=𝑎_𝑗𝑖=0 for 𝑖≠𝑗, 𝑖,𝑗=1,2. The operator 𝐵 is the identity operator, thus of order 𝑑=0. From the first consequence, it is known that 𝑢∈𝐻²(Ω) satisfies −Δ𝑢=0 and 𝑢=𝑔 on 𝜕Ω, where 𝑔=1∈𝐻^5/2(Γ), 𝑔=0∈𝐻^5/2(Σ). Therefore, by applying Theorem 25, we have 𝑢=𝑢_𝐷∈𝐻³(Ω). In general, for smoother domains with 𝐶^𝑘+1,1 boundaries, solutions to (4) are elements of 𝐻^𝑘+2.
Next, we recall that, for 𝐶^1,1 domain, there is a weak solution 𝑢_𝑁∈𝐻¹(Ω) to the boundary value problem (6). We also show that the solution actually lies in 𝐻²(Ω) and if the domain is more regular, then so is the solution. More precisely, we want to show that if Ω is a domain whose boundary is of class 𝐶^𝑘+1,1, then 𝑢_𝑁 is in 𝐻^𝑘+2(Ω), where 𝑘 is a nonnegative integer. For this purpose we need Theorem 26 which implies 𝑢_𝑁∈𝐶^∞(Ω).
Choose a bounded connected domain 𝐺 with 𝐶^∞ boundary Γ₁ such that 𝐴⊂𝐺 and 𝐺⊂𝐵, where 𝐴 and 𝐵 are the domains described in Section 2. Let Ω₁ be the annulus having boundaries Γ and Γ₁, and let Ω₂ be the other annulus with boundaries Γ₁ and Σ. First, we consider the following elliptic problem on Ω₁: −Δ𝑤=0inΩ₁,𝑤=1onΓ,𝑤=𝑢_𝑁onΓ₁. (117) Since Ω₁ is bounded with compact boundaries, we have 𝐶^∞(Γ∪Γ₁)⊂𝐻^𝑘+3/2(Γ∪Γ₁). So by applying Theorem 27, we get 𝑤∈𝐻^𝑘+2(Ω₁). Since 𝑢_𝑁|_Ω1 also solves (117), then by uniqueness we have 𝑤=𝑢_𝑁|_Ω1. If 𝐴 is a domain with 𝐶^1,1 boundary Γ, then by Theorem 27 we have 𝑢_𝑁|_Ω1∈𝐻²(Ω₁). Moreover, if 𝐴 is a 𝐶^𝑘+1,1-domain, then 𝑢_𝑁|_Ω1 is in 𝐻^𝑘+2(Ω₁).
Second, we consider the following boundary value problem: −Δ𝑣+𝑣=𝑢_𝑁inΩ₂,𝜕𝑣=𝜕𝐧𝜕𝑢_𝑁𝜕𝐧onΓ₁,𝜕𝑣𝜕𝐧=𝜆onΣ. (118) Because 𝑢_𝑁|_Ω2∈𝐻¹(Ω₂), it follows that 𝑢_𝑁|_Ω2∈𝐿²(Ω₂). We have also shown that 𝑢_𝑁∈𝐶^∞(Ω). This implies that 𝑢_𝑁∈𝐶^∞(Γ₁). Hence, (𝜕𝑢_𝑁/𝜕𝐧)∈𝐻^1/2(Γ₁). Also, 𝜆∈𝐻^1/2(Σ). Since 𝐵 is a domain of class 𝐶^1,1, then by Theorem 28, we infer that (118) has a unique solution 𝑣∈𝐻²(Ω₂). Note, however, that 𝑢_𝑁|_Ω2 also solves (118). So by uniqueness of the solution, we get 𝑢_𝑁|_Ω2=𝑣. Therefore, 𝑢_𝑁|_Ω2∈𝐻²(Ω₂). Now, if domain 𝐵 is of class 𝐶^𝑘+1,1 and 𝑢_𝑁|_Ω2∈𝐻²(Ω₂)⊂𝐻¹(Ω₂), we get 𝑣∈𝐻³(Ω₂) by applying Theorem 25 and so is 𝑢_𝑁|_Ω2. Doing this recursively, we end up with 𝑢_𝑁|_Ω2∈𝐻^𝑘+2(Ω).
Hence, for 𝐶^𝑘+1,1-domains Ω(∶=𝐵⧵𝐴), if we combine 𝑢_𝑁|_Ω1∈𝐻^𝑘+2(Ω₁) and 𝑢_𝑁|_Ω2∈𝐻^𝑘+2(Ω₂), we get 𝑢_𝑁∈𝐻^𝑘+2(Ω₁∪Ω₂). Moreover, 𝑢_𝑁 is 𝐶^∞ in a neighborhood of Γ₁ because 𝑢_𝑁∈𝐶^∞(Ω). Therefore, 𝑢_𝑁∈𝐻^𝑘+2(Ω).

Remark 30. In the computation of the first-order shape derivative, since we are dealing with 𝐶^1,1-domains, we may consider 𝐻²-regularity for the solutions 𝑢_𝐷 and 𝑢_𝑁, as justified by Theorem 29.

4.5. The Shape Derivative of 𝐽

First, we state and prove the following lemma.

Lemma 31. Let Ω⊂ℝ^𝑛 be a bounded Lipschitz domain. Then the following equation: −_Ω(div𝐅)𝑔𝑑𝑥=_Ω𝐅⋅∇𝑔𝑑𝑥−_𝜕Ω𝑔(𝐅⋅𝐧)𝑑𝑠 (119) is valid for vector field 𝐅 and scalar function 𝑔 having the following regularity: (i)𝐅∈𝐻¹(Ω;ℝ^𝑛) and 𝑔∈𝐻¹(Ω);(ii)𝐅∈𝐶¹(Ω;ℝ^𝑛) and 𝑔∈𝑊^1,1(Ω).

Proof. First we recall the Gauss’ divergence theorem in ℝ^𝑛 saying that if a domain Ω⊂ℝ^𝑛 is a bounded Lipschitz domain, then we have _Ωdiv𝐅𝑑𝑥=_𝜕Ω𝐅⋅𝐧𝑑𝑠 (120) for a vector field 𝐅∈𝐶¹(Ω,ℝ^𝑛). Second we take the divergence of the product of a scalar function 𝑔 and the vector field 𝐅 to get div𝑔𝐅=𝑔div𝐅+𝐅⋅∇𝑔. (121) Then, integrating both sides over Ω and applying the divergence theorem to the vector field 𝑔𝐅 we obtain (119).
(i) If 𝐅∈𝐻¹(Ω;ℝ^𝑛) and 𝑔∈𝐻¹(Ω), then (div𝐅)𝑔∈𝐿¹(Ω),  𝐅⋅∇𝑔∈𝐿¹(Ω), and the integral ∫_𝜕Ω𝑔(𝐅⋅𝐧) is bounded. Hence, (119) is well defined. Note that the formula −_Ω(div𝐅)𝑔𝑑𝑥=_Ω𝐅⋅∇𝑔𝑑𝑥−_𝜕Ω𝑔(𝐅⋅𝐧)𝑑𝑠 (122) holds for 𝐅∈𝐶¹(Ω;ℝ^𝑛) and 𝑔∈𝐶¹(Ω). We write −_Ω(div𝐅)𝑔𝑑𝑥=_Ω𝐅⋅∇𝑔𝑑𝑥−_𝜕Ω𝛾₀𝑔𝛾₀𝐅⋅𝐧𝑑𝑠, (123) where 𝛾₀∶𝐻¹(Ω)→𝐿²(𝜕Ω) is a trace operator. Let 𝐅∈𝐻¹(Ω;ℝ^𝑛) and 𝑔∈𝐻¹(Ω). By density, we pick {𝐅_𝑘}⊂𝐶¹(Ω and {𝑔_𝑘}⊂𝐶¹(Ω) such that 𝐅_𝑘→𝐅 in 𝐻¹(Ω;ℝ^𝑛) and 𝑔_𝑘→𝑔 in 𝐻¹(Ω). By (122), −_Ωdiv𝐅_𝑘𝑔_𝑘𝑑𝑥=_Ω𝐅_𝑘⋅∇𝑔_𝑘𝑑𝑥−_𝜕Ω𝑔_𝑘𝐅_𝑘⋅𝐧𝑑𝑠. (124) Note that 𝛾₀𝐅_𝑘=𝐅_𝑘|_𝜕Ω and 𝛾₀𝑔_𝑘=𝑔_𝑘|_𝜕Ω. Also 𝐅_𝑘→𝐅 in 𝐻¹(Ω;ℝ^𝑛) implies that div𝐅_𝑘→div𝐅 in 𝐿²(Ω). Moreover, 𝑔_𝑘→𝑔 in 𝐻¹(Ω) implies ∇𝑔_𝑘→∇𝑔 in 𝐿²(Ω). Furthermore, since 𝛾₀∈ℒ(𝐻¹(Ω),𝐿²(𝜕Ω)), 𝛾₀𝐅_𝑘→𝛾₀𝐅 and 𝛾₀𝑔_𝑘→𝛾₀𝑔 in 𝐿²(𝜕Ω). Therefore, (119) holds for 𝐅∈𝐻¹(Ω;ℝ^𝑛) and 𝑔∈𝐻¹(Ω).
(ii) If 𝐅∈𝐶¹(Ω;ℝ^𝑛) and 𝑔∈𝑊^1,1(Ω), then (div𝐅)𝑔∈𝐿¹(Ω) and 𝐅⋅∇𝑔∈𝐿¹(Ω). Note that 𝛾𝑔∈𝐿¹(𝜕Ω) (cf. [25, page 316]), where 𝛾𝑔 refers to the trace of 𝑔 on 𝜕Ω, hence 𝑔(𝐅⋅𝐧)∈𝐿¹(𝜕Ω). Therefore, (119) is also well defined for this case. Using similar arguments as above and using the density of 𝐶¹(Ω) in 𝑊^1,1(Ω) we can show that (119) is also valid for this case.

Now we apply Lemma 31 to prove the next lemma.

Lemma 32 (see [1]). Let 𝑢_𝐷 and 𝑢_𝑁 belonging to 𝐻²(Ω) satisfy the Dirichlet problem (4) and the Neumann problem (6), respectively. Then _Ω𝐴∇𝑢_𝐷⋅∇𝑢_𝐷=−_Σ∇𝑢_𝐷⋅𝐧²𝐕⋅𝐧,_Ω𝐴∇𝑢_𝑁⋅∇𝑢_𝑁=_Σ∇𝑢_𝑁²𝐕⋅𝐧−2𝜆_Σ𝐕⋅∇𝑢_𝑁, (125)(126) where 𝐴 is given by property (12) of Lemma 11.

Proof. From Lemma 11, we recall the expression for 𝐴, which is given by 𝐴=(div𝐕)𝐼−((𝐷𝐕)+(𝐷𝐕)^𝑇). Our first goal is to derive an expression for ∫_Ω𝐴∇𝑢⋅∇𝑣 for 𝑢,𝑣∈𝐻²(Ω). We begin by writing ∫_Ω𝐴∇𝑢⋅∇𝑣 as follows: _Ω𝐴∇𝑢⋅∇𝑣=_Ω(div𝐕)∇𝑢⋅∇𝑣−_Ω()−𝐷𝐕∇𝑢⋅∇𝑣_Ω(𝐷𝐕)^𝑇∇𝑢⋅∇𝑣. (127) We manipulate each term on the right-hand side of (127). First, because 𝑢,𝑣∈𝐻²(Ω), we have ∇𝑢⋅∇𝑣∈𝑊^1,1(Ω). Hence we can use (119) by taking 𝑔=∇𝑢⋅∇𝑣 and by choosing 𝐅=𝐕. In addition, we take into account that 𝐕 vanishes on the fixed boundary Γ. This leads to −_Ω(div𝐕)∇𝑢⋅∇𝑣=_Ω∇()−∇𝑢⋅∇𝑣⋅𝐕_Σ(∇𝑢⋅∇𝑣)𝐕⋅𝐧. (128) The other two terms on the right hand side of (127) are manipulated as follows. The term ∇(∇𝑢⋅∇𝑣)⋅𝐕 is written as ∇(∇𝑢⋅∇𝑣)⋅𝐕=(∇𝑣)^𝑇∇²𝑢𝐕+𝐕^𝑇∇²𝑣∇𝑢, (129) where ∇²𝑢 represents the Hessian of 𝑢. Because Hessian is symmetric, we obtain ∇(∇𝑢⋅∇𝑣)⋅𝐕=∇²𝑢∇𝑣+∇²𝑣∇𝑢⋅𝐕. (130) Substituting (130) into (128) we get _Ω(div𝐕)∇𝑢⋅∇𝑣=−_Ω∇²𝑢∇𝑣+∇²𝑣∇𝑢⋅𝐕+_Σ(∇𝑢⋅∇𝑣)𝐕⋅𝐧. (131) Next, we expand the expression div((𝐕⋅∇𝑢)∇𝑣) using (121) as ))=()()=())div((𝐕⋅∇𝑢∇𝑣𝐕⋅∇𝑢Δ𝑣+∇𝑣⋅𝐕⋅∇𝑢𝐕⋅∇𝑢Δ𝑣+∇𝑣⋅(𝐷𝐕^𝑇∇𝑢+∇²=()+()𝑢𝐕𝐷𝐕∇𝑣⋅∇𝑢∇𝑣^𝑇∇²𝑢𝐕+(𝐕⋅∇𝑢)Δ𝑣. (132)
Integrating both sides of (132) over Ω, applying Stoke’s theorem, and considering 𝐕=0 on Γ we end up with _Ω(𝐷𝐕)∇𝑣⋅∇𝑢+_Ω(∇𝑣)^𝑇∇²𝑢𝐕+_Ω()=𝐕⋅∇𝑢Δ𝑣_Σ(𝐕⋅∇𝑢)∇𝑣⋅𝐧, (133) or equivalently _Ω(𝐷𝐕)∇𝑣⋅∇𝑢=−_Ω(Δ𝑣)(𝐕⋅∇𝑢)−_Ω∇²𝑢∇𝑣⋅𝐕+_Σ𝜕𝑣().𝜕𝐧𝐕⋅∇𝑢 (134) Interchanging 𝑣 and 𝑢 we get _Ω(𝐷𝐕)∇𝑢⋅∇𝑣=−_Ω(Δ𝑢)𝐕⋅∇𝑣−_Ω∇²𝑣∇𝑢⋅𝐕+_Σ𝜕𝑢().𝜕𝐧𝐕⋅∇𝑣 (135) Also, because (𝐷𝐕)^𝑇∇𝑣⋅∇𝑢=(𝐷𝐕)∇𝑢⋅∇𝑣, we obtain _Ω(𝐷𝐕)^𝑇∇𝑣⋅∇𝑢=−_Ω(Δ𝑢)𝐕⋅∇𝑣−_Ω∇²𝑣∇𝑢⋅𝐕+_Σ𝜕𝑢().𝜕𝐧𝐕⋅∇𝑣 (136) Thus, _Ω(𝐷𝐕)^𝑇∇𝑢⋅∇𝑣=−_Ω(Δ𝑣)𝐕⋅∇𝑢−_Ω∇²𝑢∇𝑣⋅𝐕+_Σ𝜕𝑣().𝜕𝐧𝐕⋅∇𝑢 (137) Adding (131), (135), and (137) altogether, we express (127) as _Ω𝐴∇𝑢⋅∇𝑣=_Ω()+Δ𝑢𝐕⋅∇𝑣_ΩΔ𝑣(𝐕⋅∇𝑢)−_Σ𝜕𝑣()−𝜕𝐧𝐕⋅∇𝑢_Σ𝜕𝑢()+𝜕𝐧𝐕⋅∇𝑣_Σ(∇𝑢⋅∇𝑣)𝐕⋅𝐧. (138) Set 𝑢=𝑣=𝑢_𝐷 in (138). The first two integrals on the right hand side of (138) vanish because −Δ𝑢_𝐷=0 in Ω. Moreover, since 𝑢_𝐷=0 on Σ we have ∇𝑢_𝐷=(𝜕𝑢_𝐷/𝜕𝐧)𝐧. Thus, we can write (138) as follows: _Ω𝐴∇𝑢_𝐷⋅∇𝑢_𝐷=−2_Σ∇𝑢_𝐷⋅𝐧𝐕⋅∇𝑢_𝐷+_Σ∇𝑢_𝐷⋅∇𝑢_𝐷𝐕⋅𝐧=−2_Σ𝜕𝑢_𝐷𝜕𝐧𝜕𝑢_𝐷+𝜕𝐧𝐕⋅𝐧_Σ𝜕𝑢_𝐷𝜕𝐧²𝐕⋅𝐧=−_Σ∇𝑢_𝐷⋅𝐧²𝐕⋅𝐧. (139) Therefore, (125) is satisfied.
On the other hand, by replacing both 𝑢 and 𝑣 by 𝑢_𝑁 and by considering that −Δ𝑢_𝑁=0 in Ω and 𝜕𝑢_𝑁/𝜕𝐧=𝜆 on Σ, we derive (126) as _Ω𝐴∇𝑢_𝑁⋅∇𝑢_𝑁=−2_Σ𝜕𝑢_𝑁𝜕𝐧𝐕⋅∇𝑢_𝑁+_Σ∇𝑢_𝑁⋅∇𝑢_𝑁=𝐕⋅𝐧_Σ∇𝑢_𝑁²𝐕⋅𝐧−2𝜆_Σ𝐕⋅∇𝑢_𝑁. (140)

Now, we derive the explicit form of the first-order shape derivative of 𝐽.

Theorem 33. For 𝐶^1,1 bounded domain Ω, the first-order shape derivative of the Kohn-Vogelius cost functional 𝐽(Ω)=12_Ω∇𝑢_𝐷−𝑢_𝑁²𝑑𝑥 (141) in the direction of a perturbation field 𝐕∈Θ, where Θ is defined by (13) and the state functions 𝑢_𝐷 and 𝑢_𝑁 satisfy the Dirichlet problem (4) and the Neumann problem (6), respectively, is given by ()=1𝑑𝐽Ω;𝐕2_Σ𝜆²−∇𝑢_𝐷⋅𝐧²+2𝜆𝜅𝑢_𝑁−∇𝑢_𝑁⋅𝜏²𝐕⋅𝐧𝑑𝑠, (142) where 𝐧 is the unit exterior normal vector to Σ, 𝜏 is a unit tangent vector to Σ, and 𝜅 is the mean curvature of Σ.

Proof. First we consider the functionals defined on the reference domain and perturbed domains 𝐽(Ω)=12_Ω∇𝑢_𝐷−𝑢_𝑁²𝐽Ω𝑑𝑥,_𝑡=12_Ω𝑡∇𝑢_𝐷,𝑡−𝑢_𝑁,𝑡²𝑑𝑥_𝑡. (143) Let 𝑧_𝑡=𝑢_𝐷,𝑡−𝑢_𝑁,𝑡 and 𝑧=𝑢_𝐷−𝑢_𝑁. Note that if 𝑧^𝑡=𝑧_𝑡∘𝑇_𝑡, then we have 𝐷𝑧^𝑡=(𝐷𝑧_𝑡∘𝑇_𝑡)(𝐷𝑇_𝑡). Hence 𝑀_𝑡∇𝑧^𝑡=∇𝑧_𝑡∘𝑇_𝑡. (144) By this together with Lemma 16, we can write 𝐽(Ω_𝑡) as follows: 𝐽Ω_𝑡=12_Ω𝑡∇𝑧_𝑡²𝑑𝑥_𝑡=12_Ω𝐼_𝑡∇𝑧_𝑡∘𝑇_𝑡²=1𝑑𝑥2_Ω𝐼_𝑡𝑀_𝑡∇𝑧^𝑡21𝑑𝑥=2_Ω∇𝑧^𝑡,𝐴_𝑡∇𝑧^𝑡. (145) Then we write 𝐽(Ω_𝑡)−𝐽(Ω) as 𝐽Ω_𝑡(Ω)=1−𝐽2_Ω∇𝑧^𝑡,𝐴_𝑡∇𝑧^𝑡−⟨=1∇𝑧,∇𝑧⟩𝑑𝑥2_Ω∇𝑧^𝑡,𝐴_𝑡−𝐼∇𝑧^𝑡+∇𝑧^𝑡2−|∇𝑧|²𝑑𝑥=∶𝐽₁(𝑡)+𝐽₂(𝑡). (146) Using Lemma 11, the symmetry of 𝐴, and noting that 𝑧^𝑡→𝑧 as 𝑡→0 we have 𝐽₁(𝑡)=12_Ω∇𝑧^𝑡,𝐴_𝑡−𝐼∇𝑧^𝑡𝑑𝑥,𝐽₁(0)=0. (147) Hence, 𝐽̇₁(0)=lim_𝑡→0𝐽₁(𝑡)−𝐽₁(0)𝑡=12lim_𝑡→0Ω∇𝑧^𝑡,𝐴_𝑡−𝐼𝑡∇𝑧^𝑡=1𝑑𝑥2_Ω⟨⟩=1∇𝑧,𝐴∇𝑧𝑑𝑥2_Ω∇𝑢_𝐷−𝑢_𝑁,𝐴∇𝑢_𝐷−𝑢_𝑁=1𝑑𝑥2_Ω𝐴∇𝑢_𝐷⋅∇𝑢_𝐷1𝑑𝑥+2_Ω𝐴∇𝑢_𝑁⋅∇𝑢_𝑁−𝑑𝑥_Ω𝐴∇𝑢_𝐷⋅∇𝑢_𝑁𝑑𝑥. (148) To manipulate 𝐽₂(𝑡) we use the identity 𝑎²−𝑏²=(𝑎−𝑏)²+2𝑏(𝑎−𝑏) as 𝐽₂(𝑡)=12_Ω∇𝑧^𝑡2−|∇𝑧|²=1𝑑𝑥2_Ω∇𝑧^𝑡−𝑧²𝑑𝑥+_Ω∇𝑧,∇𝑧^𝑡−𝑧𝑑𝑥=∶𝐽₂₁(𝑡)+𝐽₂₂(𝑡). (149)𝐽₂₁(𝑡) is manipulated as follows: 𝐽₂₁(𝑡)=12_Ω∇𝑧^𝑡−𝑧²=12_Ω∇𝑢^𝑡_𝐷−𝑢^𝑡_𝑁−𝑢_𝐷−𝑢_𝑁²≤_Ω∇𝑢^𝑡_𝐷−𝑢_𝐷²+_Ω∇𝑢^𝑡_𝑁−𝑢_𝑁². (150) Applying Theorems 22 and 23 yields 𝐽̇₂₁=0. 𝐽₂₂(𝑡) is treated as follows: 𝐽₂₂(𝑡)=_Ω∇𝑧,∇𝑧^𝑡=−𝑧𝑑𝑥_Ω∇𝑢_𝐷−𝑢_𝑁,∇𝑢^𝑡_𝐷−𝑢_𝐷−𝑢^𝑡_𝑁−𝑢_𝑁=_Ω∇𝑢_𝐷,∇𝑢^𝑡_𝐷−𝑢_𝐷−_Ω∇𝑢_𝑁,∇𝑢^𝑡_𝐷−𝑢_𝐷−_Ω∇𝑢_𝐷−𝑢_𝑁,∇𝑢^𝑡_𝑁−𝑢_𝑁𝑑𝑥. (151) Since 𝑢^𝑡_𝐷−𝑢_𝐷∈𝐻¹₀(Ω), the variational equation for 𝑢_𝐷 and 𝑢_𝑁 implies 𝐽₂₂(𝑡)=−_Ω∇𝑢_𝐷−𝑢_𝑁,∇𝑢^𝑡_𝑁−𝑢_𝑁=𝑑𝑥_Ω∇𝑢_𝑁,∇𝑢_𝐷−𝑢_𝑁−_Ω∇𝑢^𝑡_𝑁,∇𝑢_𝐷−𝑢_𝑁𝑑𝑥. (152) Choosing 𝑢_𝐷−𝑢_𝑁∈𝐻¹_Γ,0(Ω) as test function and applying (65) and (84), 𝐽₂₂(𝑡) can be written as 𝐽₂₂(𝑡)=𝜆_Σ𝑢_𝐷−𝑢_𝑁−_Ω𝐴_𝑡∇𝑢^𝑡_𝑁,∇𝑢_𝐷−𝑢_𝑁+_Ω𝐴_𝑡−𝐼∇𝑢^𝑡_𝑁,∇𝑢_𝐷−𝑢_𝑁=𝜆_Σ1−𝑤_𝑡𝑢_𝐷−𝑢_𝑁+𝑑𝑠_Ω𝐴_𝑡−𝐼∇𝑢^𝑡_𝑁,∇𝑢_𝐷−𝑢_𝑁. (153) From this we obtain 𝐽̇₂₂(0)=𝜆_Σ𝑢_𝑁div_Σ𝐕+_Ω𝐴∇𝑢_𝑁,∇𝑢_𝐷−𝑢_𝑁, (154) and therefore, 𝐽̇₂(0)=𝜆_Σ𝑢_𝑁div_Σ𝐕+_Ω𝐴∇𝑢_𝑁,∇𝑢_𝐷−𝑢_𝑁. (155) Combining 𝐽̇₁(0) and 𝐽̇₂(0) we get 𝑑𝐽(Ω;𝐕)∶=lim_𝑡→0+𝐽Ω_𝑡−𝐽(Ω)𝑡=𝐽̇₁(0)+𝐽̇₂(0)=12_Ω𝐴∇𝑢_𝐷⋅∇𝑢_𝐷+12_Ω𝐴∇𝑢_𝑁⋅∇𝑢_𝑁−_Ω𝐴∇𝑢_𝐷⋅∇𝑢_𝑁+𝜆_Σdiv_Σ𝐕𝑢_𝑁+_Ω𝐴∇𝑢_𝑁,∇𝑢_𝐷−𝑢_𝑁=12_Ω𝐴∇𝑢_𝐷⋅∇𝑢_𝐷−12_Ω𝐴∇𝑢_𝑁⋅∇𝑢_𝑁+𝜆_Σdiv_Σ𝐕𝑢_𝑁. (156)
We know from the previous section that 𝑢_𝐷 and 𝑢_𝑁 exist in 𝐻²(Ω) since Ω is of class 𝐶^1,1. Using this smoothness we can now apply Lemma 32 and write (156) as follows: ()1𝑑𝐽Ω;𝐕=−2_Σ∇𝑢_𝐷⋅𝐧²−1𝐕⋅𝐧2_Σ∇𝑢_𝑁²𝐕⋅𝐧−2𝜆_Σ𝐕⋅∇𝑢_𝑁+𝜆_Σdiv_Σ𝐕𝑢_𝑁1=−2_Σ∇𝑢_𝐷⋅𝐧²𝐕⋅𝐧+𝜆_Σ𝐕⋅∇𝑢_𝑁+𝑢_𝑁div_Σ𝐕−12_Σ∇𝑢_𝑁²𝐕⋅𝐧. (157) Since _Σ∇𝑢_𝑁⋅𝐕+𝑢_𝑁div_Σ𝐕=_Σ𝜆+𝑢_𝑁div_Σ𝐧𝐕⋅𝐧 (158) and |∇𝑢_𝑁|²=𝜆²+(∇𝑢_𝑁⋅𝜏)², we obtain the first-order shape derivative of 𝐽: ()=1𝑑𝐽Ω;𝐕2_Σ−𝜆²−∇𝑢_𝑁⋅𝜏²+2𝜆²+2𝜆𝑢_𝑁𝜅−∇𝑢_𝐷⋅𝐧²𝐕⋅𝐧. (159)

5. Conclusion

In this paper we derived the explicit form of the first-order Eulerian shape derivative of the Kohn-Vogelius cost functional given by (7) in a rigorous manner. As seen in the presentation, we can avoid working on the shape derivatives of the states and apply their Hölder continuity instead. We employed techniques similar to [9, 13] but it was not necessary to introduce adjoint variables. For the shape derivative of the cost functional to be well defined we observe that we can consider domains with 𝐶^1,1 boundaries and we need 𝐻² regularity for the state variables.

Rewriting the first-order shape derivative as 𝑑𝐽(Ω;𝐕)=∫_Σ𝐹𝐕⋅𝐧, where 1𝐹=2−∇𝑢_𝑁⋅𝜏²−∇𝑢_𝐷⋅𝐧²+𝜆²+2𝜆𝑢_𝑁𝜅, (160) we conclude that 𝐽 is shape differentiable at Ω. This is because 𝑑𝐽(Ω;𝐕) exists for all 𝐕∈Θ and the mapping 𝐕↦𝑑𝐽(Ω;𝐕) is linear and continuous with respect to 𝐕∈Θ since |𝑑𝐽(Ω;𝐕)|≤|𝐹|_𝐿¹_(Σ)|𝐕|_𝐶(Σ)≤|𝐹|_𝐿¹_(Σ)|𝐕|_𝐶^1,1_(𝑈). (161)

We also observe that the shape derivative of 𝐽 which is formulated from the Bernoulli free boundary problem depends on the normal component of the deformation field 𝐕 at the free boundary Σ; that is, there exists a function 𝑔_Ω defined on the free boundary Σ such that 𝑑𝐽(Ω;𝐕)=_Σ𝑔_Ω𝐕⋅𝐧𝑑𝑠. (162) This agrees with the Hadamard structure theorem [26, 27].

Theorem 34 (see [26, page 318]). Let Ω be a domain with 𝐶^𝑘+1 boundary Γ for some integer 𝑘≥0. Assuming that at Ω a shape gradient of 𝐽(Ω) exists. Then there exists a scalar distribution 𝑔 in 𝒟^𝑘(Γ)^′ such that 𝑑𝐽(Ω;𝐕)=𝑔(Ω),𝐕⋅𝐧_{𝒟𝑘(Γ)′×𝒟𝑘(Γ)}, (163) where 𝒟^𝑘(Γ) is the space of 𝐶^𝑘 functions from Γ to ℝ^𝑛, and 𝐕⋅𝐧 is the normal component of 𝐕 on Γ.

Proof. See [26].

Acknowledgments

The paper is partially supported by the ÖAD—Austrian Agency for International Cooperation in Education and Research for the Technologiestipendien Südostasien (Doktorat) scholarship in the frame of the ASEA-Uninet; the SFB Research Center ‘‘Mathematical Optimization and Applications in Biomedical Sciences” SFB F32; the University of the Philippines Baguio. The first author is thankful to the Institute for Mathematics and Scientific Computing, University of Graz, for funding a presentation of partial results during the 7th European Conference on Elliptic and Parabolic Problems held on May 21–25, 2012 at Gaeta, Italy. He is likewise grateful to Professor Gilbert Peralta for his helpful insights and suggestions. Moreover, he is thanking UP Baguio for giving a Ph.D. incentive grant for writing this paper. Last but not least, the authors would like to thank the referees for devoting time to review this material.