Geometric Aspects of Shape Optimization

We present a review of known results in shape optimization from the point of view of Geometric Analysis. This paper is devoted to the mathematical aspects of the shape optimization theory. We focus on the theory of gradient flows of objective functions and their regularizations. Shape optimization is a part of calculus of variations which uses the geometry. Shape optimization is also related to the free boundary problems in the theory of Partial Differential Equations. We consider smooth perturbations of geometrical domains in order to develop the shape calculus for the analysis of shape optimization problems. There are many applications of such a framework, in solid and fluid mechanics as well as in the solution of inverse problems. For the sake of simplicity we consider model problems, in principle in two spatial dimensions. However, the methods presented are used as well in three spatial dimensions. We present a result on the convergence of the shape gradient method for a model problem. To our best knowledge it is the first result of convergence in shape optimization. The complete proofs of some results are presented in report (Plotnikov and Sokolowski, Gradient flow for Kohn–Vogelius functional).


Introduction
The theory of shape optimization is a mathematical discipline that is located at the intersection of the calculus of variations and the theory of free boundary value problems. Historically, its beginning is attributed to the time of the appearance of Newton's studies (1685, Principia Mathematica) of the problem of finding the shape of a body which moves in a fluid with minimal resistance to motion. It seems that the paper [72] was the first publication devoted to shape optimization problems in the mechanics of solids. The other early direction included problems of optimizing the eigenvalues of elliptic operators. See monograph [55] for references and the historical remarks. We should also mention pioneering works [45,46] on the application of variational methods to problems of ideal fluid flows with free boundaries.
The beginning of the modern mathematical theory of shape optimization was laid in monographs [34,93,104]. In these monographs, it was first singled out as an independent scientific discipline. At present, the theory of shape optimization includes a large number of various applied problems. A number of different approaches have been developed to solve shape optimization problems. The purpose of this paper is to give the reader an idea of the main problems of the theory and methods for their solution. We will focus on the geometric aspects of the theory.
Typically, the shape optimization problem admits the following general formulation. First, a fixed bounded of the Euclidean space R d , d = 2, 3, is specified. It is supposed to contain the inclusion i such that i ⊂ . The shape of the inclusion is unknown and must be determined together with the solution of the boundary value problem. It is also assumed that the regions i and e = \ i are filled with some physical substances. Among such substances can be solid elastic materials, liquids, physical fields (electric field, gravitation field), or simply void. The state of each substance is described by solutions to the system of governing partial differential equations equipped with appropriate boundary conditions. These solutions are completely determined by the inclusion i . In this framework, the compact set = ∂ i defines the interface between domains occupied by materials with different physical properties.
Finally, an objective function J is specified. Usually it is considered as a function of i and solutions to the governing equations. Its value is completely determined by the inclusion i . Therefore, we will denote the objective function as J ( i ) or equivalently J ( ). The shape optimization problem is to find i that minimizes the objective function, Here the minimum is taken over the admissible set of inclusions. Let us give some basic typical examples of applied shape optimization problems. The simplest are problems of shape identification in electric tomography and geophysics. Electrical impedance tomography is used in medical imaging to reconstruct the electric conductivity of a part of the body from measurements of currents and voltages at the surface [26]. The same technique is also used in geophysical explorations. An important special case consists in reconstructing the shape of an unknown inclusion or void assuming (piecewise) constant conductivities.
Transmission single measurement identification problem Let us assume that a material occupy the bounded region in the space of points x ∈ R d , d = 2, 3. Without loss of generality, we may assume that the boundary of is infinitely differentiable. Furthermore, we assume there are two disjoint open arcs N , D ⊂ ∂ such that cl N ∪ D = ∂ . The inclusion, which is unknown and must be determined together with the solution, occupies the subdomain i with the boundary . The equilibrium equations for the electric field potential u : → R in the simplest case can be written as div (a∇u) = 0 in , (

1.2)
Here n is the outward normal vector to ∂ , h n is a given voltage, and h d is a given distribution of electric potential. We assume that h n and h d are extended to ∂ and h n ∈ L 2 (∂ ), h d ∈ W 1/2,2 (∂ ). (1. 3) The conductivity a is defined by the equalities a = 1 in e , a = a 0 in i , (1.4) where a 0 is a given positive constant, If D = ∅, then for every h n and h d satisfying condition (1.3), problem (1.2) admits a unique solution u ∈ W 1,2 ( ). If in addition, the arcs N , D belong to different connected components of ∂ , h n , h d ∈ C ∞ (∂ ), and ∂ , ∂ i belong to the class C ∞ , then u ∈ C ∞ ( ). The problem on the identification of the inclusion i is formulated as follows. For a given function g : D → R it is necessary to find an inclusion i such that the solution to problem (1.2) satisfies the extra boundary condition a∇u · n = g on D . (1.5) It is assumed that g satisfies the orthogonality condition D g ds + N h n ds = 0.
More generally, the problem of identification is to determine the shape of the inclusion by the additional boundary condition. This inverse problem is ill-posed and in general case has no solution. In practice, its approximate solution can be found by solving the variational problem min i ∈A J ( i ), (1.6) where the objective function J ( i ) is a positive function that vanishes if and only if a solution to problem (1.2) satisfies the condition (1.5), A is some class of admissible inclusions. Notice that the mapping i → u, where u is a weak solution to problem (1.2), determines a nonlinear operator, which takes the set of admissible shapes A into W 1,2 ( ). The most successful choice of the objective function is the Kohn-Vogelius energy functional, which is defined as follows: [61] Here v, w : → R satisfy the equations and boundary conditions div a∇v = 0div a∇w = 0 in , a∇v · n = g w = h d on D , a∇v · n = h n a∇w · n = h n on N . (1.8)

Single measurement identification problem with void
The other example is the electrical impedance tomography problem that can be formulated as follows. For a domain ⊂ R 2 and functions h ∈ W 1/2,2 (∂ ), g ∈ W −1/2,2 (∂ ), to find subdomain i and function u ∈ W 1,2 ( ) such that (1.9) These equations define an overdetermined boundary value problem which has a solution only for the true inclusion i . Following Roche and Sokolowski, [97] we can replace boundary value problem (1.9) by the variational problem for Kohn-Vogelius type functional. To this end, denote by v and w solutions to the boundary value problems (1. 10) In this case the Kohn-Vogelius functional reads Shape optimization problems in mechanics of solids Again consider the standard geometric configuration that consists of bounded domain ⊂ R d , d = 2, 3, and the inclusion i . The state of linear elastic solid is completely characterized by the displacement field u : → R d satisfying the equilibrium equation − div (A e(u)) = F, (1.12) where F is a given mass force, the strain tensor e, and the Hooke's law matrix A are defined by the equality 2e(u) = ∇u + ∇u , A = A i χ i (x) + A e (1 − χ i (x)). (1.13) Here the characteristic function χ i : → {0, 1} of the domain i is defined by the equality The constant matrices A β , β = i, e, with entries A β lmpq characterize the properties of the elastic material and satisfies the symmetry and positivity conditions: c −1 |ξ | 2 ≤ c Aξ : ξ ≤ c|ξ | 2 for all symmetric matrices ξ. (1.14) Note that the stress tensor σ is defined by the equality σ (u) = A e(u). Equation (1.12) should be endowed with boundary conditions. For example, we can take the Neumann and Dirichlet boundary conditions in the form where h n and h d are given tractions and displacements, N and D are open disjoint subsets of ∂ such that N ∪ D = ∂ . There are various formulations of the shape optimization problems in the solid mechanics corresponding to different objective functions. The typical choice of an objective function is We also can consider the single measurement identification problem for elastic material similar to transmission single measurement identification problem formulated above.

Shape optimization problems in fluid mechanics
The considerations of hydrodynamical forces acting on the object traveling within fluid is fundamental to the design of aircrafts, cars, and in many other practical problems. The design of optimal shapes with minimal (maximal) drag is one of the most important problems of applied hydrodynamics. It can be regarded as the shape optimization problem for equations of fluid dynamics. This problem was widely discussed in the literature. We refer the reader to review [75] and to monograph [94]. Again, assume that ⊂ R d , d = 2, 3, is a hold all bounded domain with the smooth boundary ∂ . It is supposed that contains a nonpermeable body i with the boundary . A viscous incompressible fluid occupies the flow domain e = \ i . The state of the fluid is completely characterized by the velocity field u : e → R d and the pressure function p : e → R, which satisfy the Navier-Stokes equations and the boundary conditions (1.16) where the constant vector u ∞ is the flow direction. If ∂ is Lipschitz, then this boundary value problem admits at least one solution u ∈ W 1,2 ( e ), p ∈ L 2 ( e ). If in addition ∂ and belong to the class C l+α with l ≥ 2 and α ∈ (0, 1), then u ∈ C l+α ( e ) and p ∈ C l−1+α ( e ).
The drag F D is the projection of the hydrodynamics force, acting on the body, onto the direction u ∞ , i.e., It was proved in the seminal paper [11] that the expression for the drag can be equivalently rewritten in the form of the volume integral It should be noted that the absolute minimum drag is achieved with an empty set i . Hence the drag minimization problem only makes sense if there are additional constraints on the geometry of i , which guarantee the nontriviality of solution. As such constraints, we can choose the area (length) L of the boundary = ∂ i ds = fixed positive constant or the volume of the body i dx = fixed positive constant.
In addition, in the two-dimensional case, we can define the lift (lifting force) F L as the projection of the hydrodynamic force onto the direction orthogonal to u ∞ , The problem of minimizing the drag for a given lifting force, as well as the problem of the ratio F D /F L optimization, are natural optimum design problems for the Navier-Stokes equations, see e.g., [47,63]. We listed the main applications of optimization theory to problems in solid and fluid mechanics. In fact, the theory of shape optimization finds applications in various fields of science, for example, in biology, [5], and photonics, [64]. Methods Unfortunately, shape optimization problems as stated with no additional geometric constraints are usually ill-posed, see [60,79,107] for examples. The reason is that microstructures tend to form, which are associated with a weak convergence of the characteristic functions χ m i along a minimizing sequence m i , m ≥ 1. Indeed, in the absence of strong compactness of the minimizing sequences of designs, the optimal state should be attained by a fine mixture of different phases. There are two different ways to cope with these difficulties.
First, the well-posed problems can be generated by a relaxation (homogenization) procedure. The homogenization of the material properties lead to the formation of microstructures. In such a way, the set of admissible shapes is extended and includes the microstructures. The quasi-convexification of the integrand in J is performed by taking the infimum over all possible microstructures, therefore, the existence of minimizers is ensured. The relaxation procedure usually yields continuous design variables χ i over the reference domain . In such a case, it is impossible to define any shape from the homogenized solution for solids, liquids, or voids. Hence the relaxed optimal solutions may not lead directly to practical designs. The analysis of the relaxation method is beyond the scope of this paper. We refer the reader e.g., to monographs [19,25], and paper [4] for a description of the relaxation method.
The second approach is the regularization of the objective function with the geometric energy functionals. The first-order penalization of J reads: (1.19) where L is the perimeter of i , p > 0 is the regularization parameter. If = ∂ i is a regular manifold, then L is the area of in 3D case and the length of in 2D case. We refer to monograph [42] for the theory of sets with finite perimeter (Caccioppoli sets). This penalization was proposed in [8] by analogy with the Mumford-Shah functional, [78], in the theory of image segmentation processes. Note that the appearance of the perimeter regularization is motivated by the difficulties regarding the mathematical treatment of shape optimization. If the shape optimization problem is additionally supplemented with a perimeter penalization, then positive results concerning existence of optimal shapes have been obtained (see for instance [105]). However, sets with finite perimeter may be irregular in general case. Hence penalization (1.19) can be regarded as a weak regularization of shape optimization problems. The stronger regularization may be obtained if we impose constraints on the curvatures of . This approach also was motivated by the theory of image processing, [77]. The only possible conformally and geometrically invariant penalization functional depending on curvatures is the Willmore functional defined by the equality (1.20) where H is the mean curvature of . We refer the reader to monographs [53,115], for the basic theory of surfaces with finite Willmore energy. In 2D case E e coincides with the famous Euler elastica functional. Therefore, we can define the strong regularization of an objective function as follows: Here j , j = e, p, are some positive constants. Note that the penalization term can be interpreted as the cost of structure manufacturing. Hence j are not necessary supposed to be small.

Remark 1.1
On the other hand, the influence of the geometric energy penalization on the optimal design should be further studied both from theoretical and numerical points of view. It is known that in the case of level set method with the topological derivatives adding the perimeter requires special construction of the numerical method of solution to obtain useful optimal designs, see [9].
The most important question of the theory is the construction of a robust algorithm for the numerical study of shape optimization problems. The standard approach is to use the steepest descent method based on the shape calculus developed by Sokolowski and Zolesio [104]. See also Delfour and Zolesio [34], and references therein. The shape calculus works for inclusions i with the regular boundary = ∂ i . In this setting, the objective function J is considered as a functional defined on the totality of smooth curves . This assumption is natural from the practical point of view. Without loss of generality we may restrict our considerations to the class of twice differentiable immersions (parametrized surfaces, curves) f : In this framework, we will use the notation J ( f ) along with the notation J ( ). The main goal of the shape calculus is to develop the method of differentiation of objective functions with respect to shapes of geometrical objects.
Following the general method of the shape calculus, we define the shape derivative of an objective function. To this end, choose an arbitrary vector field X : S d−1 → R d and consider the immersion The manifolds t = f t (S d−1 ), t ∈ (−1, 1), define the one-parametric family of perturbations of . The shape derivativeJ of J in the direction X is defined by the equalityJ If it admits the Hadamard representatioṅ where n is the inward normal to = ∂ i , then the vector field is said to be the gradient of J at the point f . The same definition holds for the geometric energy functional E.

The Steepest Descent Method and the Gradient Flow
It follows from the definition that the shape gradient d J can be regarded as a normal vector field on . If f is sufficiently smooth, for example f ∈ C 2+α , then the mapping f + δ d J ( f ) defines an immersion of S d−1 into R d for all sufficiently small δ > 0. In the steepest descent method, the optimal immersion f and the corresponding shape = f (S d−1 ) are determined as a limit of the sequence of immersions 25) and the corresponding sequence of surfaces n = f n (S d−1 ). Here the energy E is defined (1.21), δ is a fixed positive number, usually small, f 0 is an arbitrary admissible initial shape. Relation (1.25) can be considered as the time discretization of the Cauchy problem is a decreasing function of t, a solution to problem (1.26) can be considered as approximate solution to the penalized variational problem Hence the existence of a solution to Cauchy problem (1.26) guarantees the wellposedness of the steepest descent method. In its turn, the existence of the limit lim t→∞ f (t) guarantees the convergence of the method. This paper is devoted to the mathematical aspects of the shape optimization theory. We focus on the theory of gradient flows of objective functions and their regularization. However, a number of important ideas and methods are left out of the scope of this article. For example: 1. Topological optimization, which is based on the concept of a topological derivative used in the level set type method, [30-32, 80, 82-85, 103]. 2. The theory of homogenization method developed in [3,4]. 3. Application of direct methods of the calculus of variations using the theory of capacity, [34]. 4. Shape optimization problems with uncertainty conditions and random data, [27,28,[50][51][52]. 5. The optimal layout theory in optimum design, [14,15,65,72,87].
The paper is organized as follows. Shape sensitivity analysis is one of the main tools of the theory of shape optimization. In Sect. 2, we present the outline of main ideas of the shape calculus. In order to be clear, we restrict the considerations to the relatively simple example of the single measurement identification problem. We give the derivation of the basic formulas for the material derivatives of the solutions to this problem and derive the representations of shape derivatives of objective functions. The formulations are given both in the distributed form and in the form of a contour integral in the Hadamard form. In the general case, shape optimization problems can be attributed to the class of problems with free boundaries of mathematical physics. Such problems are difficult for mathematical analysis. Their numerical solution encounters significant difficulties. There are several approaches that help simplify the problem. In the next two sections, we will cover two of the most popular approaches: the phase field method and the level set method. Sect. 3 is devoted to modeling shape optimization problems using the phase field (diffusive surface) method. This method allows us to reduce the original problem with a free boundary to a boundary value problem for a weakly nonlinear system of parabolic-elliptic equations. In this section, we give the construction of a phase field approximation for shape optimization problems in rigid body mechanics and viscous fluid dynamics. We will also proceed with the derivation of the phase field equations for the corresponding gradient flows. Sect. 4 contains a description of the level set method, which is one of the most common methods for studying shape optimization problems. This method is a special algorithm for the numerical solution of optimization problems. It is based on the representation of the moving surface of the gradient flow of an objective function in the form of a solution to the Hamilton-Jacobi equation. A rigorous mathematical justification of the level set method is hardly possible, but it allows constructing efficient numerical algorithms.
In the last Section A, we consider the question of the correctness of the theory of gradient flows for shape optimization problems. For the model problem of identification of the inclusion form we establish the existence of a smooth solution of the equations.

Shape Calculus
In this section, we give the outline of the main ideas of the shape calculus theory.
This theory traces its origins to Hadamard's pioneering paper [49]. Now the shape calculus is one of the main mathematical tools of the general shape optimization theory. We refer the reader to monographs [10,34,37,54,104,109] and papers [2,86,102] for details and references.
In order to make the explanation more clear we restrict our considerations by the 2D single measurement identification problem and the simplest scalar version of the compliance problem. We start with the analysis of the shape derivative of the transmission problem for the Laplace equation.
Assume as before that an electric field occupies one-connected bounded domain ⊂ R 2 with smooth Jordan ∂ . Furthermore assume that an inclusion occupies a bounded domain i with i ⊂ . Denote by the boundary of i and set e = \ i . Suppose also that the conductivity a : → R satisfies the condition The problem is to find the electric field u : → R satisfying the following equations and boundary conditions: Here u − and u + are restrictions of u on e and i , ν is the outward normal vector to ∂ , n is the inward normal vector to ∂ i = , g is a given voltage. Further we will assume that the given function g satisfies the solvability condition These equations can be rewritten in the equivalent form div (a∇u ) = 0 in , a∇u = g on ∂ . (2.

Material Shape Derivative of Solution to Problem (2.1.4)
The definition of the shape derivatives of solutions to problem (2.1.4) is based on the following construction. Choose an arbitrary mapping ϕ : → R 2 of the class C ∞ ( ) and consider the family of C ∞ mappings y t : → defined by the equality (2.2.1) By the contraction mapping principle, the mapping y t takes diffeomorphically the domain onto itself for all t from the small interval (−t * , t * ). Here the small positive t * depends only on ϕ. Obviously, y t coincides with the identical mapping outside of the support of ϕ. Moreover, it is an analytic function of t in the interval (−t * , t * ).
The diffeomorphism y t defines the one-parametric families of the sets t i , t , and the functions a t , They can be regarded as the perturbations of the inclusion i , the interface , and the conductivity coefficient a. The perturbed electric potential u t (y) serves as a weak solution to the elliptic boundary value problem It is clear that u ≡ u 0 is a weak solution to problem (2.1.4), (2.1.6). In other words, u t defines the perturbation of the original solution u. The calculation of the derivative u t with respect to t (the Eulerian derivative) is difficult, since the set of discontinuity of the coefficient a t strongly depends on t. Hence the derivative ∂ t u t at t = 0 can be defined only outside of . In order to cope with this difficulty, the theory of shape calculus deals with the so-called material derivative which is defined as follows. Introduce the one-parametric family of functions v t : → R given by the equality The material derivativeu : → R is defined by the relatioṅ Here u is a solution to problem (2.1.4), (2.1.6). The limit is taken in some suitable Banach space. In our case, an appropriate space is W 1,2 ( ). Now our task is to obtain the effective representation for the derivativesu andü. Recall that u t is a solution to boundary problem (2.2.3). The change of the independent variable in (2.2.3) leads to the following equations for the Here the symmetric positive matrix N is defined by the equalities where the notation dϕ stands for the Jacobi matrix of the mapping ϕ. By virtue of the Neumann theorem, we have It follows that the matrix N admits the decomposition Note that for a suitable choice of t * , the series in the right-hand side converges in any Calculations give the following representations for the first two terms in decomposition (2.2.8).
The following lemma shows that the formula for S 2 can be essentially simplified.

Lemma 2.1 Under the above assumptions, we have
Proof Introduce the temporary notation It is necessary to prove that We begin with the observation that which yields On the other hand, we have We thus get (2.2.12) Next, we have Combining this result with (2.2.12) we arrive at the identity which obviously yields desired equality (2.2.11) Let us turn to the derivation of the representations foru andü. Notice that a weak solution to problem (2.2.6) satisfies the integral identity for all test functions ζ ∈ W 1,2 ( ). The first integral in the left-hand side of this integral identity defines the positive continuous sesquilinear form in the Hilbert space of all functions v ∈ W 1,2 ( ) with zero average over ∂ . Moreover, this form is analytic function of the parameter t on the interval (−t * , t * ). It follows from the analytic theory of perturbations of self-adjoint operators, [59], ch. 7, that the weak solution v t to problem (2.2.13) is an analytic function of the parameter t ∈ (−t * , t * ). Moreover, v t admits the representation The series in the right-hand side converges strongly in the space W 1,2 ( ). It is clear that the material derivativeu[ϕ] and the second material derivativeü [ϕ, ϕ] in the direction of the vector field ϕ are defined by the equalitieṡ In order to complete the derivation of the material derivatives of u, it remains to obtain the equations for v 1 and v 2 . Substituting decompositions (2.2.8) and (2.2.14) into (2.2.13) and retaining the first three terms in the obtained equality we conclude that the integral identities hold for all test functions ζ ∈ W 1,2 ( ). It follows that v 1 and v 2 serve as weak solutions to the boundary value problems div (a∇v 1 ) = −div (aS 1 ∇u) in ,

Distributed Shape Derivatives of the Kohn-Vogelius Functional
We define the objective function J by the equality where u is the solution to problem (2.1.4), (2.1.6). The diffeomorphism y t (x) = (I + tϕ)(x), ϕ ∈ C ∞ 0 ( ), takes the curve , the inclusion i , and the coefficient a to the perturbed t , t i , and a t given by relations (2.2.2). Let u t be a solution to perturbed problem (2.2.3). Thus we get the one-parametric family J (t) perturbation of J given by the equalities where N and v t are defined by (2.2.6) and (2.2.7). The shape derivative of the functional J in the direction ϕ are given by the equalitieṡ .
In order to obtain the representation for the shape derivatives of J , we substitute the decompositions (2.2.8) and (2.2.14) into (2.3.1) to obtain It is clear that the power series at the right-hand side converges on the interval (−t * , t * ). Thus we getJ The direct calculations show that On the other hand, identities (2.2.16) and (2.2.17) with ζ replaced by u and v 1 imply Substituting these equalities into (2.3.3) and recalling relations (2.2.15) we finally obtainJ give the desired representation for the shape derivatives of J .

Hadamard Representation of the First-Order Derivative of the Kohn-Vogelius Functional
The representation of the shape derivatives of the Kohn-Vogelius functional given by formulae (2.3.4) depends on the vector field ϕ defined in the whole domain . However, the perturbation of J must depend on the perturbation of the interface . Therefore, we may expect that the shape derivatives depends only on the restriction of the vector field ϕ on . In other words, the integrals in (2.3.4) should be independent of an extension ϕ| to . This means that the area integrals in (2.3.4) can be reduced to the integrals over the interface . This leads to the so-called Hadamard formulae for the shape derivatives. In this subsection we obtain the Hadamard representation for the first-order derivativeJ . Our considerations are based on the following auxiliary lemma. Assume that is a C 2 Jordan curve. Let us consider two vector fields p, q : → R 2 satisfying the following conditions: where p − , q − are restrictions of p, q on e and p + , q + are restrictions of p, q on i . We will denote by · the jumps across the interface . For example, we have Next set The similar notation holds for q.

Lemma 2.2 Under the above assumptions we have
Proof Note that From this and the equalities Integrating by parts we obtain Next, we have Substituting this result into (2.4.6) we finally arrive at desired identity (2.4.4).
We are now in a position to derive the representation for the first derivative of the objective function J . The result is given by the following proposition.

Proposition 2.3
Let a weak solution u to problem (2.2.6) and the interface satisfy the condition Then we haveJ The proof is based on Lemma 2.2. Set p = a∇u, q = ∇u.
Since a = const in domains e and i , u is a harmonic function in these domains. It follows that Applying Lemma 2.2 we obtain Since u and a∇u are continuous in , we have Substituting this result into (2.4.9) and recalling formula (2.3.4) forJ we finally obtaiṅ and the proposition follows.
Recall that the gradient d J of an arbitrary objective function is defined by the equalities Thus we get the following:

Corollary 2.4 Under the assumptions of Proposition 2.3, we have
Remark 2.5 Let us consider the "complementary" problem with the Dirichlet boundary condition for the functional Arguing as in the proof of Corollary 2.4 we obtain the following equality:

The Second-Order Shape Derivative of Kohn-Vogelius Functional
In this subsection we derive the Hadamard representation for the second-order deriva-tiveJ . The result is given by the following proposition.
Then we have We split the proof into a sequence of lemmas.

Lemma 2.7
Under the assumptions of Proposition 2.6, we have where D 1 is given by (2.5.3).

Lemma 2.8
Under the assumptions of Proposition 2.6, we have It follows that It remains to note that and the lemma follows.

Lemma 2.9
Under the assumptions of Proposition 2.6, we have We thus get It remains to note that and the lemma follows.

Lemma 2.10
Under the assumptions of Proposition 2.6, we have Proof Since u is harmonic function in \ , it follows from (2.5.6) and (2.5.9) that

Next, we have
Repeating these arguments we conclude that equalities (2.5.10) hold for A 22 + B 22 and A 21 + B 21 .

Lemma 2.12
Under the assumption of Proposition 2.6, we have We begin with the observation that −div (S 1 ∇u) = div (dϕ + dϕ − div ϕ I )∇u .
We have Since u = 0 in \ , straightforward calculations lead the representation − div (S 1 ∇u) = R 1 + R 2 , (2.5.14) where Substituting the expression for R i into (2.5.14) and noting that we finally obtain equality (2.5.13).

Lemma 2.13
Under the assumptions of Proposition 2.6, we have where D i are given by (2.5.7) and (2.5.9).
Proof It follows from Lemma 2.11 that On the other hand, Lemma 2.12 yields Combining these equalities we arrive at (2.5.15).
We are now in a position to complete the proof of Proposition 2.6. Equality (2.5.4) reads Since a is a constant in each component of \ , equality (2.5.15) implies It follows from (2.5.7) and (2.5.9) that Noting that and This leads to the equality Substituting this relation into (2.5.16) we finally arrive at the desired representation

Examples
For many objective functions, there exist the explicit expressions for their gradients.
Below we list some of them.

Transmission Single Measurement Identification Problem
In this case the gradient d J of the Kohn-Vogelius objective function (1.7) is defined as follows, see [2] d J = 2 a∂ n v ∂ n v − a∂ n w ∂ n w ) n − a∇v · ∇v − a∇w · ∇w n, (2.5.17) where · denotes the jump across , n is the inward normal to ∂ i = , v, and w are solutions to equations (1.8).

Single Measurement Identification Problem with Void
In this case the gradient of the objective function is defined by the equality, [97], dJ = (∂ n v 2 − ∂ n w 2 ) n, (2.5.18) where v and w are solutions to problem (1.10).

Drag Minimization Problem for Navier-Stokes Equations
In the drag minimization problem, the objective function J = F D is defined by formulae (1.17) and (1.18). The analysis of the shape derivative and the gradient of J have been conducted in [11,92]. In particular, [92] gives the following expression of d J , The state (u, p) and the costate (w, q) are given by the solutions to the boundary value problems Hence the shape gradient d J , in contrast to the shape gradients of Kohn-Vogelius type functionals, depends on shape derivatives of solutions to governing equations. It is important to note that formula (2.5.19) makes sense if and only if solutions to problems (2.5.20) and (2.5.21) are unique. The uniqueness criterium was discussed in [92] and [47]. We only note that solutions are unique if the viscosity ν is sufficiently large with respect to |u ∞ | and the diameter of .

Preliminaries
In the previous sections, we consider shape optimization problems in a fixed hold-all domain ∈ R d , d = 2, 3, containing an unknown inclusion i . The goal was to find the shape of i in order to minimize the objective function J , depending on i and some physical field u defined in or in i or in e = \ i . The unknown inclusion is completely characterized by the design variable ρ : → {−1, 1} defined by the equality is the characteristic function of i . Hence ρ = 1 in i and ρ = −1 in e . The main idea of the phase field method (diffusive surface method) is to replace the discontinuous design variable ρ by a continuous phase field function (an order parameter) ϕ : → R. The interface between the material components in the phase field approximation can be roughly represented as a small neighborhood of the level set {ϕ = 0} named diffusive interface. The domains i and e occupied by different materials approximately correspond to the sets {ϕ > 0} and {ϕ < 0}. Hence no restrictions are imposed on the topology of the diffusive interface. In addition, the boundaries of the sets occupied by different components may intersect with the boundary of the hold all domain . The latter is important for the analysis of the compliance problem in solid mechanics, see [12,13] for the statement of the compliance problem. Thus, the phase field theory is one of the most appropriate mathematical tools for solving topological optimization problems.
The phase field method first was developed as a way to represent the surface dynamics of phase transition phenomena. It dates back to historical papers [7,20]. This method has been used in many interface dynamic studies as a general interface tracking method. The first application of the phase field method to the structural optimization problem was given in [16]. The idea of applying of the phase field method to shape optimization in hydrodynamics first was proposed in [47]. In this section, we shortly describe the possible applications of the phase field theory to the basic shape optimization problems. It is important to note that the phase field method works only for penalized problems with the perimeter or (and) elastic penalization. Hence, first we consider the phase field approximation for geometric problems. Throughout this section the notation W stands for the Ginzburg-Landau potential Phase field approximation of geometric problems In this paragraph, we consider the phase field approximation for the area functional and the Willmore-Helfrich functional. The approach is based on the general theory of -limit. We refer the reader to monographs [17,25] for the state of art in the domain. Following [17], the -limit of the sequence of functionals is defined as follows.
Definition 3.1 Let us consider a family of functionals F ε : X → [−∞, ∞], ε ∈ (0, 1), defined on a topological space X . In that case we say that F ε -converges to where N (x) denotes the family of all neighborhoods of x in X . In this case we say that F(x) is the -limit of F ε at x and we write If this equality holds for all x ∈ X then we say that F ε -converges to F (on the whole X ).
For the first time, the phase field approximation of the area functional was considered in paper [74], see also [73]. In these papers, it was shown that for every bounded Lipschitz domain ⊂ R d , d = 2, 3, the -limit in L 1 ( ) of the energy functional and F ε (ϕ) = ∞ otherwise admits the representation.
Here H d−1 is (d − 1)-dimensional Hausdorff measure, ∂ * is the essential De Giorgi boundary, see [17], the constant c W is defined by the equality The question on the phase field approximation of the Willmore-Helfrich functional was raised by De Giorgi, [29], Conjecture 4. Later De Giorgi's conjecture was modified in the context of the phase transition theory, see [98] and references therein. The modified De Giorgi's problem can be formulated as follows. Let ⊂ R d be a bounded domain with Lipschitz boundary, let W be the double-well potential defined by (3.2). For every > 0 and γ > 0 define the functional F ε : L 1 ( ) → R by the equalities if ϕ ∈ W 2,2 ( ) and F ε (ϕ) = ∞ if ϕ ∈ L 1 ( )\W 2,2 ( ). Now introduce a set i ⊂ such that the ∂ i ∩ belongs to the class C 2 . Denote by ρ the design variable given by (3.1). In [98] it was proved that Here H is the mean curvature of ∂ i ∩ . The requirement that the interface ∩∂ i be smooth is essential. There are numerous examples of two-dimensional Euler's elastica with singular points. The phase field approximation does not work in neighborhoods of such points. It should also be noted that the phase field approximation (3.4) is neither the only one and nor the most accurate one. For the discussion of the problem, see for instance [18] and references therein.
It follows from what was mentioned above that the functionals (3.3) and (3.4) define the phase field approximations of the area and Willmore-Helfrich functionals. For sufficiently smooth vector fields ϕ and smooth compactly supported perturbations δϕ it follows that L 2 -gradients d F and dF of the functionals F ε and F are defined by the equalities Calculations show that where the chemical potential µ is defined by

5)
Gradient flow There is a massive body of literature devoted to L 2 -gradient flow (mean curvature flow) of the area phase field approximation F given by (3.3). We refer the reader to papers [57,58] for details. The gradient flow of the Willmore functionals was considered in the papers [69], [38], and [44]. It follows from these papers that L 2gradient flow equation for the functional F ε reads where µ is given by (3.5). For the sake of simplicity, we take = 1 and rewrite this equation in the form The natural boundary and initial conditions for equation (3.6) can be taken in the form Equation (3.6) along with boundary and initial conditions (3.7) determine the wellposed boundary value problem for weakly nonlinear fourth-order parabolic equation. The global existence and uniqueness of strong solutions for this and more general phase field models are proved in [23].

Phase Field Approximations of Objective Functions
Two-components problems The main goal of the shape optimization theory is minimizing an objective function J as well as minimizing its penalization E + J . Therefore, in order to develop a theory of the phase field approximation for shape optimization problems, it is necessary to determine the phase field approximations of the objective function J . It is important to note that the theory of the phase field models can be developed only for penalized objective functions E + J or L + J . In the previous paragraph, we considered the phase field approximations of geometric functionals L and E. Now In this setting, the functions v and w satisfy the equations and boundary conditions div a∇v = 0 d i va∇w = 0 i n , where h ∈ W 1/2,2 (∂ ), g ∈ L 2 (∂ ) are given functions, n is the outward normal to ∂ . We also add the following conditions: which guarantee the solvability and uniqueness of solutions to problem (3.8). In the theory of the phase field models, the discontinuous design variable ρ is replaced by a continuous phase function ϕ. Therefore, it is most natural to choose the phase approximation of the conductivity a(ρ) in the form a = a(ϕ) such that a monotone function a ∈ C ∞ (R) satisfies the conditions Thus we arrive at the following expression for the phase field approximation of the Kohn-Vogelius functional J = a(ϕ) |∇v − ∇w| 2 dx, (3.10) where the functions v and w satisfy the equations and boundary conditions div (a(ϕ)∇v) = 0 d i v(a(ϕ)∇w) = 0 i n , and additional conditions (3.9). It should be noted that the solutions to problems (3.11) are uniquely determined by the phase field ϕ. Hence J = J (ϕ) is a function of the phase variable. Now, our task is to derive the expression for the gradient of the Kohn-Vogelius functional (3.10). Recall that the gradient d J of an objective function J is defined by the equality where ψ is an arbitrary smooth function. The gradient of the phase field approximation of the Kohn-Vogelius functional is given by the following lemma. where v and w are solutions to boundary value problems (3.11)-(3.9).
Proof Without loss of generality we may assume that ϕ ∈ C ∞ ( ). Fix an arbitrary ψ ∈ C ∞ ( ) and consider the one-parametric family of coefficients a s = a(ϕ + sψ), s ∈ (−1, 1). (3.13) Denote by v s , w s ∈ W 1,2 ( ) the solutions to the boundary value problems Since a s is an infinitely differentiable function of s, it follows from the general theory of elliptic equations, see [59], that v s and w s are infinitely differentiable functions of s with values in W 1,2 ( ). Set Differentiation of equations (3.14) with respect to s at s = 0 gives the following system of equations and boundary conditions for ω and ς : which are understood in the sense of distributions. This means that ω ∈ W 1,2 ( ), ς ∈ W 1,2 0 ( ), and the integral identities hold for all ξ ∈ W 1,2 ( ) and η ∈ W 1,2 0 ( ). Next notice that for s = 0, Let us evaluate the integral in the right-hand side. Since ς ∈ W 1,2 0 ( ), it follows from (3.11) that a∇v · ∇ς dx = a∇w · ∇ς dx = 0.
It follows from this and the first integral identity in (3.15) Combining the obtained results we arrive at the identity Substituting this equality into (3.16) we obtain the desired relation (3.12).
The following consequence of Lemma 3.2 may be useful. Let us consider the functionals where v and w are solutions to boundary value problems (3.11).

Corollary 3.3 Under the assumptions of Lemma 3.2, we have
Proof Notice that if h = 0 (g = 0), then w = 0 (v = 0). Hence the corollary is the straightforward consequence of Lemma 3.2.
Compliance minimization problem. Assume that a two-component elastic material occupies the hold all bounded domain ⊂ R d , d = 2, 3, with the C ∞ boundary ∂ . Denote by i and e = \ i the domains occupied different components. Usually i is considered as an elastic inclusion. It is important to note that in the compliance minimization problem, it is not assumed that ∂ i ∩ ∂ = ∅. The state of the material is characterized by the displacement vector field u : → R d . Introduce the strain tensor e and the Hooke's law stiffness matrix A defined by the equalities Here χ i is the characteristic function of domain i , the constant matrices A i and A e characterize the properties of the elastic material and satisfies the symmetry and positivity conditions (1.14). With this notation the stress tensor σ is defined by the equality σ = Ae(u). In the absence of material forces, the equilibrium equation reads div σ ≡ div (A e(u)) = 0 in . This equation should be supplemented with boundary conditions. We take the traction boundary conditions where g is a given force acting on the material. The compliance J equals the work done by the applied load (3.17) By virtue of the equilibrium equation, we can rewrite this expression in the equivalent form The functional J can be regarded as an objective function for the compliance minimization problem. Recall the notation ρ = 2χ i − 1 for the design variable. It is easily seen that In the theory of phase field approximation, the design variable ρ should be replaced by a continuous phase field function ϕ. We choose the phase field approximation of the tensor A in the form where b ∈ C ∞ (R) is a monotone function such that b(ϕ) = 0 for ϕ ≤ −1 and b(ϕ) = 1 for ϕ ≥ 1. Obviously A(ϕ) satisfies the symmetry and positivity conditions (1.14). Thus we get the following phase field approximation of the objective function: Arguing as in the proof of Lemma 3.2 and Corollary 3.3 we arrive at the following formula for the gradient of J : : e(u). (3.18) Phase field approximation of objective functions. One-component problems The important characteristic of the phase field models is that the phase function should be defined in the hold all domain . The peculiarity of one-component problems is the presence of void where the phase function is not defined. To overcome this difficulty, we make the classical assumption that the void region say e is filled with some fictitious material. Below we give two examples of the application of this approach. One-component compliance problem. In one-component compliance problem, it is assumed that the domain e is a void region. The means that A e = 0. Next, we suppose that the domain e is filled with some fictitious material with the Hooke's law matrix A e = δ A * , δ ∈ (0, 1). We also assume that a constant matrix A * satisfies the symmetry and positivity condition (1.14). Now take the phase field matrix A(ϕ) in the form  (3.20) where the constant vector u ∞ is the flow direction. The objective function J is the projection of the hydrodynamics force, acting on the body, onto the direction u ∞ . Calculations show that in the case of the viscous incompressible fluid it is defined by the equality, see [11,92], Note that the drag minimization problem makes sense if there are additional constraints on the geometry of i , which guarantee the nontriviality of solution. As such constraints, we choose the area (length) L of the boundary = ∂ i , ds = given positive constant.
In order to derive the expression for the phase field approximation of the objective function and its gradient, we have to define the equations for fictitious fluid in the whole domain . The following approach was proposed in [47]. Notice that the flow domain e approximately coincides with set {ϕ < 0}. On the other hand, the streamlined body i approximately coincides with the set {ϕ > 0}. Now introduce a one-parametric family of functions a δ (ϕ) with the following properties: Next we introduce the fictitious resistance force −a δ (ϕ)u and define the phase field approximation of the Navier-Stokes equations as follows: We take the phase field approximation of the objective function in the form Our task is to derive the expression for the gradient of J by using the shape calculus. On this way we meet the following critical difficulty, which is typical for nonlinear problems. Notice that for every ϕ ∈ L ∞ ( ), problem (3.23) has a weak solution u ∈ W 1,2 ( ), p ∈ L 2 ( ). In particular, for ν > ν 0 > 0, this solution admits the estimate It is worthy noting that c u is independent of a δ and ϕ. Moreover, if ϕ ∈ C l+α ( ), then u ∈ C l+2+α ( ). However, this solution is not unique and the number of solutions may depend on ϕ. If it is the case, then we cannot define the shape gradient of the objective function. In order to cope with this difficulty we impose additional restriction on the data and the viscosity coefficient in order to provide the uniqueness of solutions to problem (3.22). It is known, see [56], that a solution u to problem (3.22) with a δ = 0 is unique, if it satisfies the inequality where In particular, u is unique if ν > ν 0 > 0 with arbitrary fixed ν 0 and ν > √ dC D c u (u ∞ , , ν 0 ). (3.25) This result also holds true for an arbitrary positive a. Further we will assume that the viscosity coefficient satisfies inequality (3.25). In order to derive the formula for d J we proceed as in the proof of Lemma 3.2. We restrict our considerations by the case of fixed small parameter δ and will write simply a(ϕ) instead of a δ (ϕ).
Denote by (u s , p s ) ∈ C ∞ ( ) the solution to the boundary value problem −ν u s + div (u s ⊗ u s ) + ∇ p s + a s u s = 0, div u s = 0 in , Since a s is an infinitely differentiable function of s, it follows that u s and p s are infinitely differentiable functions of s. Set Differentiation of equations (3.26) with respect to s at s = 0 gives the equations and boundary conditions for υ and q Note that for s = 0, we have Integrating by parts and noting that div u = 0, we obtain Combining this relation with equation (3.26) we arrive at the equality Here we use the identity The analysis of the uniqueness proof in [56] shows that condition (3.25) guarantees the existence and boundedness of the linear operator Next notice that In other words, ω is a solution to the adjoint boundary value problem Thus we arrive to the following formula for the gradient of the phase field approximation of the drag functional (3.30)

Gradient Flows
It follows from the analysis in Sects. 3.1 and 3.2 that in the framework of the phase field theory, the shape optimization problem is reduced to problem of minimization of the functionals Here F ε and F ε are given by equalities ( In all these cases the functionals are weakly lower semicontinuous in suitable Sobolev spaces. In particular, the existence of minimizers can be proved by the application of the, direct methods of the calculus of variations. The justification of the steepest descent method requires study of the gradient flows of the functionals (3.31). The correctness of the evolutionary gradient flow equation guarantees the well-posedness of this method. We restrict our consideration be the case of strong regularization with the Willmore-Helfrich penalty functional. It follows from (3.6) that can take the gradient flow equations in the form The problems with constraints. In many applications the additional constraints on admissible shapes are imposed. We consider in details the case of the perimeter constraints.
The key observation is that in the phase field theory the lengths element ds of the unknown interface admits the approximation Hence for = 1, the corresponding constraints condition can be written in the form

35)
L 0 is a given positive constant such that In this case parameter γ in equation (3.32) becomes a function of the temporal variable and can be regarded as the Lagrange multiplier. It is easily seen that In this case the function C is independent of t and equals C(0). At the end of this section we give two examples of gradient flow equations for the problems listed in Sect. 3.2. Two-components compliance minimization problem. In this case we have the following elliptic-parabolic problem for the displacement field u and the phase function ϕ.
Here the matrix A(ϕ) is defined by (3.19). Drag minimization problem. In this case we have the parabolic equation for ϕ, the Navier-Stokes equation with fictitious resistance force for the velocity u and the pressure p, and the adjoin linear equation for the adjoint variables ω and π .

Level Set Method
The level set method first introduced and devised in [91] is a simple and efficient method for computing the motion of an interface in two or three dimensions. The level set method has a wide range of applications, including problems in fluid mechanics, solids mechanics, and image processing over the years [5,90,91,100]. The level set methods for structural shape and topology optimization problems are proposed and developed in papers [5,101,111,112]. There is now a growing massive of literature devoted to the application of the level set method to numerous applied shape optimization problems. We refer the reader to a review [6] and references therein for the state of the art in this domain. However, it is important to note that in general this method does not have a rigorous mathematical justification. It remains unclear whether the differential equations underlying this method are mathematically correct. Nevertheless, the level set method remains a powerful and efficient method for the numerical analysis of applied problems. In this section, we give a brief outline of the main ideas of the level set method. In order to illustrate the main features of the method, we restrict our considerations by a simple geometric configuration. Assume that a two-component material occupies a bounded domain ⊂ R 2 with the smooth boundary. Assume that the components occupy two disjoint subdomains i and e of the hold all domain separated by an interface , i.e., i ∪ e ∪ = . For simplicity we assume that the inclusion i is a simply connected domain with regular boundary . The main idea is to define a one-parametric family of moving surfaces (curves) (t) such that the objective function J ( (t)) decreases as the quasi-time variable t decreases. It separates the moving domains i (t) and e (t).
The level set method occupies an intermediate position between the phase field method and the gradient flow method. Just as in the phase field model, we introduce a phase function ϕ(x, t) satisfying the conditions (4.1) In this setting the interface (t) coincides with the level set {ϕ(·, t) = 0}. Hence the task is to define the evolution of the phase function ϕ(x, t).
In contrast to the phase field model, the phase function ϕ is determined directly from the kinematic equation, which describes the motion of the interface points along the trajectories of some velocity field V (x, t). Denote by V (x, t), x ∈ (t), the velocity of the motions of points x ∈ (t). It assumes that V is extended to such that the extended field V (x, t) is Lipschitz in the variable x and the normal component of V vanishes at ∂ . The evolution of the phase function along the field V is defined by a solution to the Cauchy problem for linear transport equation Here ϕ 0 is the initial distribution of the phase function. A moving interface is defined as a level set {ϕ(·, t) = 0}. The function ϕ can be defined by the characteristic method by the equality ϕ(x, t) = ϕ 0 (X (t, x, 0)), (4.3) where X (t, x, s) is a solution to the Cauchy problem d ds The tangent component of the field V generates the sliding of the interface along itself. Therefore, it is always assumed that the vector field V is directed along the normal to . Note that the normal to the level surface is oriented to the side increasing phase function. We thus get where v is some scalar field. As a result, we arrive at the Hamilton-Jacobi equation for the phase function This equation is widely used in wave optics and acoustics as a mathematical model of the motion of wave fronts. Almost all results on the global existence and uniqueness of solutions to the Hamilton-Jacobi equations have been obtained using the method of viscosity solutions. We refer the reader to the basic monograph [68] and review article [24] to get acquainted with the general theory of viscosity solutions. Note that the viscosity solutions method is applicable to autonomous equations with a Lipschitz velocity field and Lipschitz initial data. In this case, the solution is also Lipschitz.
To obtain the equation of the level set method for optimization problems, the velocity v must be specified. The main requirement is that the objective function J ( (t)) decreases as t increases. This leads to the following expression for the restriction v of the scalar field v to the interface (t): Here, H is an arbitrary smooth function satisfying the conditions where v * is an extension of the vector field v given by (4.5) to the hold all domain . The question of constructing such an extension is nontrivial. We refer the reader to paper [1] for a discussion of this issue.
In particular, v * strongly depends on the choice of an extension operator. Note that v * and v depend on in some implicit and complicated manner. In fact, equation The rigorous treatment of this specific case was developed in [21,43]. See also important monograph [48] and the references therein.
In the conclusion of the section, following [6], we describe the recurrent process for finding approximate solutions to problem (4.7). The process is divided into three steps.
1. The first step is to discretize the process with respect to time. For this purpose, the interval (0, T ) is divided into subintervals 2. The second step is to determine the velocity field v * n . Assume that the approximate velocity v * a and the approximate phase function ϕ a are well defined. Set Following [6] define the vector field v * n ∈ W 1,2 0 ( ) as a weak solution to the transmission problem, where is a small parameter.
3. Finally, the approximate solution ϕ a on the interval [τ n−1 , τ n ] is defined as a solution to the Cauchy problem If 0 is sufficiently smooth, then the process is well defined. However, the proof of its convergence is problematical.
robust algorithm for numeric calculations of an optimal shape. In this section, we give outline of main ideas of the proofs of existence and smoothness results for the gradient flows in the shape optimization theory. In order to be clear, we restrict our considerations to the relatively simple 2D single measurement identification problem for the Kohn-Vogelius functional. Recall the formulation of this problem given in Sect. 1.

Kohn-Vogelius functional
Let us consider the Kohn-Vogelius energy functional, which is defined as follows, [61]: Here v, w : → R satisfy the equations and boundary conditions div a∇v = 0 d i va∇w = 0 i n , a∇v · n = g w = h on ∂ . ∂ , ∈ C 2+α , h ∈ C 2+α (∂ ), g ∈ C 1+α (∂ ), α ∈ (0, 1). (5.6) Denote by v − , w + the restrictions of v, w to e and by v + , w + the restrictions of v, w to i . It follows from the Schauder estimates for solutions to elliptic equations that v − , w − ∈ C 2+α ( e ) and v + , w + ∈ C 2+α ( i ). For every function with − and + continuous in e and i , the notation stands for the jump of across , For strong solutions to transmission problem (5.4), we have a∂ n v ≡ a∇v · n = 0, a∂ n w ≡ a∇w · n = 0, v = w = 0. (5.7) With this notation the gradient d J of the Kohn-Vogelius objective function (5.3) is defined by (2.5.17),

Geometric Functionals
The standard formulation of the geometric flow equations deals with parametrized curves (surfaces). Further we will assume that the interface admits the representation = f (S 1 ), where the immersion f : S 1 → R is unknown and should be defined along with the solution to the geometric flow problem. Note that f a 2π -periodic function of the angle variable θ ∈ R/2π Z. The element of the length of equals where g is the only nontrivial coefficient of the first fundamental form of the curve . In this setting, the derivative with respect to the arc-length variable s, becomes the nonlinear differential operator depending on f . Hereinafter we assume that the point f (θ ) moves around in the positive counterclockwise direction while the parameter θ increases. The tangent vector and the normal vector form the positive oriented moving frame on . Notice that n is the unit inward normal vector to ∂ i = . The curvature vector k is defined by the equalities Notice that k is orthogonal to τ and is directed along the normal vector n. The Euler elastic energy E e and the perimeter L are defined by the equalities Without loss of generality we can take the penalty function of geometric energy in the form The gradient of E is given by the following lemma.
Lemma A.1 Under the above assumptions, we have Here the connection ∇ s for every vector field : → R 2 , is defined by the equality Identities (5.13) are classic (see for instance [39]). Such identities are very particular case of the 3D Willmore variation formula [115].

Gradient Flow Equations
We are now in the position to specify the gradient flow equation for the penalized Kohn-Vogelius functional. Applying Lemma A.1 we can rewrite equation (5.15) in the form Here the gradient d J is defined by relation (5.8) and can be regarded as nonlinear nonlocal operator acting on . Hence (5.16) is a nonlinear operator equation. It may be considered as a nonlocal perturbation of the elastic flow equation In the literature, this equation is also named as straightening equation and, 1D-Willmore flow equation. Now, there is almost complete theory of this equation, see [26,39,62,71,113] for details and references. We use the methods developed in these papers in our analysis.

A.2.1 Geometric Lemmata
In this subsection, we will consider special class of immersions f : S 1 → R 2 satisfying the conditions where E 0 is some fixed constant. Our consideration are based on the following elementary lemmas on the properties of such immersions. The first gives the two-side estimates for the length L in terms of the energy bound E 0 .

Lemma A.2 The estimate
holds true for every curve satisfying condition. (5.18).

Proof
The proof is given in [95].
The second lemma provides the local graph representation of planar curves with square integrable curvature. Let us consider the following construction.
For every 0 < κ < L/2 denote by κ the arc Next, introduce the Cartesian coordinates (x 1 , x 2 ) with origin at z such the axis of abscissa is directed along the tangent vector τ (θ z ) and the axis of ordinate is directed along the normal vector n(θ z ). The further results do not depend on the choice of z. Now our task is to show that the curve locally can be represented as a graph of C 1+α function in a neighborhood of z. Lemma A.3 Under the above assumptions, there exist positive numbers κ, α, β, and c, depending only on the constant E 0 in (5.18), and the function η ∈ C 1 (−α, β), η(0) = 0, with the following properties: (5.20) Here η (x 1 ) = ∂ x 1 η(x 1 ). Moreover, the mapping x 1 → (x 1 , η(x 1 )) defines C 1parametrization of the arc 3κ and takes diffeomorphically the interval (−α, β) onto this arc.
Proof The proof is given in [95].
Lemma A.3 gives the simple criterium of the absence of self-intersections of curves satisfying the energy condition (5.18).

Corollary A.4
Let an immersion f : S 1 → R 2 meets all requirements of Lemma A.3. Furthermore assume that there is ν > 0 with the property Then has no self-intersections. Conversely, if has no self-intersections, then inequality (5.21) holds for some ν > 0.
Proof The corollary is an obvious consequence of Lemma A.3.
The second corollary extends the previous results to the case of families of immersions with finite elastic energy. Let us consider a family of immersions f (t, ·) : . Every immersion f (t, ·), satisfying condition (5.18), defines L(t)periodic function of the arc-length variable s, Note that the periods L(t) are uniformly bounded from below and above by the constants 2/E 0 and E 0 . Moreover, the functions ∂ 2 s f (t, ·) are uniformly bounded in L 2 (−L(t)/2, L(t)/2). It follows that the set of the mappings f (t, ·), t ∈ [0, T ], satisfying (5.18), is relatively compact in C 1 (R).
Assume that a family of immersions f (t), t ∈ [0, T ], satisfies the following conditions: G. 1 The curves (t) = f (t, S 1 ) have no self-intersections. G. 2 The immersions f (t) satisfy energy condition (5.18) with the constant E 0 independent of t. G. 3 The set of the mappings f (t, ·), t ∈ [0, T ] is compact in the space C(S 1 , R 2 ).
It follows from Lemma A.3 that for every f (t, θ), t ∈ (0, T ), there is κ ∈ (0, 2/E 0 ) which meets all requirements of this lemma and is independent of t.
for all t ∈ [0, T ] and for all arcs 3κ (t) given by Lemma A.3.

Proof
The proof obviously follows from Lemma A.4.

A.2.2 Sobolev Spaces of Periodic Functions
For every integer r ≥ 0, denote by H r , the Sobolev space of all L -periodic mappings with the finite norm For real r ≥ 0, the space H r is defined by the interpolation. Note that the equivalent norm in H r may be defined by the equality where the Fourier coefficients If is a rectifiable Jordan curve of the length L, then the curvature of , the gradient of Kohn-Vogelius functional, tangent, and normal vectors of can be regarded as Lperiodic functions of the arc-length variable s. If depends on the temporal variable t, then we will write H r (t) instead of H r .

A.3 Estimates of the Shape Gradient dJ
In this section, we consider in details the shape gradient d J of the Kohn-Vogelius functional. Our goal is to derive the estimates of d J in the Sobolev spaces H r in terms of the geometric characteristics of the interface . By virtue of representation (5.8), the normal vector field d J : → R 2 is the quadratic form of the derivatives of solutions v, w to boundary value problem (5.4). First we derive the estimates for a general transmission problem. Assume that the interface satisfies the following conditions: H. 1 The Jordan curve ⊂ satisfies the energy condition H.2 There is ν > 0 with the property where κ, depending only on E 0 , is given by Lemma A.3 H. 3 There is ρ > 0 such that dist ( , ∂ ) > ρ. H.1-H.3 is a Jordan curve of the class C 1+α , 0 < α < 1/2. It splits the domain into two parts. The first i (inclusion) is the one-connected domain with boundary . The second is the curvilinear annulus e = \ i bounded by and ∂ . For simplicity, we will assume that ∂ is a Jordan curve of the class C ∞ .

Every curve satisfying Conditions
Next, introduce the piecewise constant function a : → R + (conductivity) defined by the equalities a = 1 in e , a = a 0 in i . If is sufficiently smooth, then w is continuous on . In other words, w − = w + on . However, the normal derivative of w has a jump across . Next set ∂ n w − = ∇w − · n, ∂ n w + = ∇w + · n on . (5.27) Our task is to estimate ∂ n w ± via the curvature of . The following theorem, [95], on the estimates of ∂ n w ± is the first main result of this section. for every q ∈ [1, ∞). In this case the constant c depends in addition on q.

A.4 Existence Theory
In this subsection, we prove the main theorem on the existence of global smooth solution to problem (5.15). Assume that the initial data satisfy the following conditions: Moreover, the Jordan curves (t) = f (t, S 1 ), t ∈ [0, T ), are separated from ∂ . If T < ∞, then there is a sequence f (t j ), t j → T as j → ∞, such that dist ( (t j ), ∂ ) → 0, or (and) f (t j ) converge in C 1 (S 1 ) as j → ∞ to some immersion f ∞ such that the limiting curve ∞ has a self-intersection.
Remark A. 9 The celebrated Li-Yau theorem and its 1D version, see [66,76], constitute that a surface (curve) has no self-intersection if its Willmore energy is not too large. It follows that (t) has no self-intersections for all t, if the initial energy E 0 is less than some critical value which is not small.
Next, the natural question is the existence of limit lim t→∞ f (t) in the case T = ∞. The proof of the existence of such limits for general gradient flows usually is based on the Lojasiewicz-Simon inequality. It is known that this inequality holds true for straightening Eq. (5.17), see [71]. For the shape optimization gradient flow Eq. (5.16) the question is open.
Proof The proof of Theorem A.8 is standard and consists of three steps. The first is the proof of the local solvability of problem (5.16) on the small time intervals. The second step is the proof of the global a priori estimates for smooth solutions to problem (5.16) in Sobolev and Hölder classes. These estimates and the extension method entail the existence of smooth solution on an arbitrary time interval.
A detailed proof of short-time existence is outside of the scope of this paper. Note that Eq. (5.16) is a degenerate parabolic equation with added lower order operator d J . In our 2D case the local existence result can be obtained by writing the flow as a graph over the initial date. In particular, the problem can be reduced to a scalar parabolic equation for the distance function, [22]. See also [39,62] for useful arguments.
Hence our main task is to derive global a priori estimates for solutions to problem (5.15). This part of the proof is technical and lengthy. Our proof is based on the estimates for Kohn-Vogelius functional given by Theorem A.7 and methods developed in [26,39,67]. The results are given by the following two theorems. The first constitutes the Sobolev a priori estimates for the curvature k as a function of the arc-length variable s. Proof The proof is given in [95].
The second theorem gives the a priori estimates for solutions to problem (5.16) in the Hölder classes. Proof The proof is given in [95].
In order to complete the proof of Theorem A.8, we use the extension method. Without loss of generality we may assume that f 0 ∈ C ∞ (S 1 ). Hence problem (5.16) has a C ∞ -solution f defined on some small interval (0, T ). By virtue of Theorem A.11, this solution meets all requirements of Theorem A.8. Let (0, T ) be the maximal interval of existence of such a solution.
Let us prove that T = ∞ in the first case. Assume, in contrary to our claim, that T < ∞. There is > 0 such that quantities ν(t) and ρ(t) are uniformly separated from zero on the interval [T − , T ), i.e., ν(t) > ν > 0 ρ(t) > ρ > 0 for some ν and ρ. Hence f (t) meet all requirements of Theorem A.11 on the interval [T − , T ) with the initial data f (T − ). It follows from this theorem that