Optimal shape control of airfoil in compressible gas flow governed by Navier-Stokes equations

The flow around a rigid obstacle is governed by compressible Navier-Stokes equations. The nonhomogeneous Dirichlet problem is considered in a bounded domain with a compact obstacle in its interior. The flight of the airflow is characterized by the work shape functional, to be minimized over a family of admissible obstacles. The lift of the airfoil is given in function of time and should be closed to the flight scenario. Therefore, the minimization for a given lift of the work functional with respect to the shape of obstacle in two spatial dimensions is considered. The shape optimization problems for the compressible Navier-Stokes equations with the nonhomogeneous Dirichlet conditions in a bounded domain with an obstacle are considered in the monograph (Plotnikov, Sokolowski, Springer, 2012) in the general case of three spatial dimensions. In the present paper, for the purposes of wellposednes of shape optimization problems, there is no restrictions on the regularity of admissible obstacles, which are simply connected compact subsets of a given, bounded hold-all domain. The presented results are derived in the general framework established in (Plotnikov, Sokolowski, Springer, 2012). It means that in two spatial dimensions we obtain the same optimal shape existence result for the compressible Navier-Stokes equations as it is for the Laplacian (Sverak, JMPA, 1993). However, the complete existence proof is substantially more complex compared to the Laplacian, and it is given with all details in (Plotnikov, Sokolowski, Springer, 2012) for the general case of three spatial dimensions. In the present paper the general theory of shape optimization developed in (Plotnikov, Sokolowski, Springer, 2012) is adapted to the particular case of two spatial dimensions. The shape optimization problem of work minimization over a class of admissible obstacles is introduced. The continuity of the work functional with respect to the obstacle in two spatial dimensions is shown for a wide class of admissible obstacles. The dependence of local solutions to the governing equations with respect to the boundary variations of obstacles is analyzed. The shape derivatives (Sokolowski, Zolesio, Springer, 1992) of solutions to the compressible Navier-Stokes equations are derived. The shape gradient (Sokolowski, Zolesio, Springer, 1992) of the work functional is obtained. The framework for numerical methods of shape optimization (Plotnikov, Sokolowski, Zochowski, MMAR'09, 2009) is established for nonstationary, compressible Navier-Stokes equations.


Introduction
Mathematical modeling of compressible viscous fluids is a new domain of intensive research, we refer the reader to [4] for the first account of modern theory of global weak solutions to compressible Navier-Stokes equations. The theory and some applications to shape optimization are developed by E. Feireisl and his co-workers, we refer to the recent paper [1] in this direction. The boundary value problems for compressible Navier-Stokes equations are considered from the point of view of shape optimization in [6]- [14].
The shape optimization of an obstacle in compressible fluid governed by the Navier-Stokes equations coupled with the transport equation is considered in this paper. We restrict ourselves to the case of two spatial dimensions. The mathematical analysis of such problems in the general case of three spatial dimensions is performed in the recent monograph [6]. The reduction of the dimension simplifies somehow our presentation for the existence theory of an optimal obstacle since in such a case a simple geometrical criterium for Kuratowski-Mosco convergence of Sobolev spaces can be employed. In order to show the existence of an obstacle for the global renormalized weak solutions of the governing equations it is required in two spatial dimensions that the admissible obstacles are simply connected compact sets.
The existence of an optimal obstacle for the work minimization is established in [6] for a wide class of admissible obstacles. In order to show the existence of an obstacle for the global, weak, renormalized solutions of the governing equations it is required only that the admissible obstacles are simply connected compact sets. In addition we need the The shape derivatives of the solutions to the equations and the shape gradient of the work functional are obtained in the framework of boundary variations technique [15,6] The recent monography [6] is devoted to the study of boundary value problems for equations of viscous gas dynamics, named compressible Navier-Stokes equations. The principal significance of the mathematical theory of Navier-Stokes equations lies in the central role they now play in fluid dynamics. In [6] we focus on existence results for the inhomogeneous in/out flow problem, in particular the problem of the flow around a body placed in a finite domain, on the stability of solutions with respect to domain perturbations, on the domain dependence of solutions to compressible Navier-Stokes equations, and on the drag optimization problem.
In the paper we focus on the problem of the compressible viscous gas flow around a moving rigid body and on the question of optimal choice of the shape of this body in order to minimize the work in the nonstationary case. To make the problem useful for applications we assume in addition that the lift is given, or the volume of the obstacle is fixed for the shape optimization problems under considerations.
The state of the fluid or gas is characterized completely by the macroscopic quantities: the density (x, t), the velocity u(x, t), and the temperature ϑ(x, t). These quantities are called state variables in the following. We assume that the temperature is a given constant ϑ.
We notice that explicit values of the Reynolds Re and Mach Ma 2 numbers are not essential for the general mathematical theory of Navier-Stokes equations, provided that these numbers are separated from 0 and bounded from above. The Coriolis force is also immaterial for the mathematical theory [6]. Hence without loss of generality we may assume that (1.1) Thus we come to the following boundary value problem in the flow domain Problem 1. For given T > 0 and for given functions find a velocity field u : Ω × [0, T ) → R 2 and a density : where S(u) = ∇u + ∇u + (λ − 1) div u, div S(u) = ∆u + λ∇ div u, and the inlet Σ in is defined by Recall that n is the unit outward normal to the boundary ∂Ω.
The problem in the stationary case can be formulated as follows: (1.4) find a velocity field u : Ω → R 2 and a density : where the inlet Σ in is defined by Remark 3. For stationary flows, if U · n = 0 on ∂Ω, the inlet is an empty set. In this particular case the total mass M of the gas should be prescribed:

Weak solutions
All known results in the global existence theory for compressible Navier-Stokes equations concern the so-called weak solutions. The notion of a weak solution arises in a natural way from the formulation of the governing equations as a system of conservation laws. In our particular case weak solutions are defined as follows: Let Ω ⊂ R 2 be a bounded domain with C 3 boundary ∂Ω and Q be the cylinder Ω×(0, T ). A couple ( , u) : Q → R + × R 2 is a weak solution to Problem (1) if the following conditions are satisfied: • The state variables satisfy , p, |u| 2 ∈ L ∞ (0, T ; L 1 (Ω)), u ∈ L 2 (0, T ; W 1,2 (Ω)), holds for all vector fields ψ ∈ C ∞ (Q) satisfying which means that the test vector fields vanish at the lateral side and the lid of the cylinder Q.
holds for all functions ζ ∈ C ∞ (Q) satisfying which means that the support of the test function ζ meets the boundary of the cylinder Q at the inlet Σ in and the bottom Ω × {0}.

Renormalized solutions
The mass balance equation (1.3b) is a first order linear partial differential equation for the density . Given a smooth solution to this equation, the composite function ϕ( ) satisfies the first order differential equation (3.1) Thus, a smooth solution to the mass balance equation generates a solution in the form ϕ( ) to the first order quasi-linear differential equation (3.1); the latter solution is associated with the prescribed function ϕ. Equation (3.1) derived for ϕ( ) is called a renormalized equation. The derivation of a renormalized equation is called renormalization. The renormalization procedure is simple for smooth solutions to the mass balance equation, but becomes nontrivial for weak solutions. It turns out that the renormalized equation is convenient for our purposes. A natural way to avoid the renormalization procedure for weak solutions to the mass balance equation is simple replacement of the mass balance equation by its renormalized version. This means that it is required that the density and the velocity satisfy (3.1) for a specific class of functions ϕ. Thus we come to the following definition: holds for all functions ζ ∈ C ∞ (Q) satisfying and for all functions ϕ ∈ C 1 0 (R). Remark 6. Obviously, any C 1 function ϕ : R → R + can be approximated by a sequence of compactly supported C 1 functions ϕ n such that ϕ n (s) → ϕ(s), ϕ n (s) → ϕ (s) as n → ∞ uniformly on compact subsets of R, Hence the integral identity (3.2) is satisfied for all C 1 functions ϕ such that ϕ( ), ϕ ( ) ∈ L 2 (Q).

Work functional for compressible gas flow around airfoil
If the obstacle S moves in atmosphere like a solid body then its physical position S t at time t is defined by S t = U(t)S + a(t).
(4.1) Here the unitary matrix U and the vector field a(t) can be defined by an appropriate flight planning scenario. Assuming that the hold-all domain is moving along with the body, after change of the variables x → U(t)x + a(t), t → t, we obtain, see [6], Ch.2, the following system equation and boundary conditions in the moving frame where the skew-symmetric matrix C = (C ij ) d×d and the vector field f * are defined by By abuse of notation we will write simply f instead of f * . The boundary and initial data are defined by the exterior conditions. In particular, if in the unmovable frame the gas is in the rest at infinity, then U = −W. From the mathematical point of view this connection is not essential and we can consider U and W as independent smooth vector fields.
The following expression for the work W S of the hydrodynamic forces acting on the moving obstacle S is given in [6]: where η is a smooth function with a compact support which contains the obstacle boundary ∂S, and η(x) ≡ 1 in a vicinity of the obstacle boundary, The second functional required for applications in aerodynamics is the lift t → L S (t), which is defined as the component of the force acting on the obstacle in the given direction defined by a unitary constant vector field ν

4)
A natural question is to minimize this work by choosing an appropriate shape of S from some suitable class S. We can briefly introduce some of the main tools developed in [6] for the domain dependence of the weak solutions to compressible Navier-Stokes equations in bounded domains. These are the kinetic equation method and the Kuratowski-Mosco convergence of function spaces. As a result the existence of optimal obstacles for the energy minimization problem can be obtained in two spatial dimensions for the family of simply connected compact sets.

Stability of solutions with respect to nonsmooth data and domain perturbations.
Propagation of rapid oscillations in compressible fluids. In compressible viscous flows, any irregularities in the initial and boundary data are transferred inside the flow domain along fluid particle trajectories. We develop a new method for the study of the propagation of rapid oscillations of the density, which can be regarded as acoustic waves. The main idea is that any rapidly oscillating sequence is associated with a parametrized family µ xt of probability measures on the real line named the Young measure. We establish that the distribution function f (x, t, s) = µ x,t (−∞, s] satisfies a differential relation named a kinetic equation. A remarkable property of compressible Navier-Stokes equations is that in this particular case the kinetic equation can be written in closed form as The kinetic equation being combined with the momentum balance equations gives a closed system of integro-differential equations which describes the propagation of rapid oscillations in a compressible viscous flow. Notice that oscillations can be induced not only by oscillations of initial and boundary data, but also by irregularities of the boundary of the flow domain. We also prove that if the data are deterministic and the function f satisfies some integrability condition, then any solution to the kinetic equation satisfying some integrability conditions is deterministic.

Domain dependance of solutions to compressible
. We show that if a sequence S n of compact obstacles converges to a compact obstacle S in the Hausdorff and the Kuratowski-Mosco sense, then the sequence of corresponding solutions to the in/out flow problem contains a subsequence which converges to a solution to the in/out flow problem in the limiting domain. Moreover, we prove that the typical cost functionals, such as the work of hydrodynamical forces, are continuous with respect to S-convergence. As a conclusion we establish the solvability of the problem of minimization of the work of hydrodynamical forces in the class of obstacles with a given fixed volume.

Optimal control of airfoil shape
Let us suppose that there is given the lift scenario for the flight t → L * (t). Thus, an optimal shape of airfoil is selected under the following conditions for admissible family of compact obstacles S ⊂ B B, where B stands for the hold-all domain in our model: • The volume of the obstacle is given M * := meas 2 (S); • The lift is prescribed in the interval (0, T ), for the fixed time horizon T > 0; • The energy used for the flight is minimized, thus the shape functional W S is considered for the optimization procedure of the airfoil. The typical shape optimization problem which is well posed in our framework can be formulated as follows.
Problem 7. Minimize the shape functional over the family of admissible obstacles S. Here α, β are arbitrary positive constants In two spatial dimensions the general results of [6] lead to the following theorem.
Theorem 8. Let p = γ , γ > 1, and div W = 0. Assume that the family of admissible obstacles S in two spatial dimensions consists of all connected compact subsets of B B. Then there is an optimal obstacle S ∈ S which minimizes J(Ω).
Proof. Let S n ∈ S be a minimizing sequence, i.e., We make use of the result by Sverak [16] which implies that there are a subsequence of the sequence S n , still denoted by S n and an obstacle S ∈ S such that S n converges to S in the Kuratowski-Mosco topology. After passing to a subsequence we may assume that S n → S in the Hausdorff metric. It follows from Theorem 10.2.15 in [6] that W Sn → W S as n → ∞. On the other hand, stability Theorem 9.3.1 in [6], implies that, after passing to a subsequence, we may assume that For any compact set Ω B \ S, The couple (u, ) is a renormalized solution to problem (4.2). Since in formula (4.4) for the lift, ∇η vanishes in a neighborhood of ∂Ω, it follows from (5.2)-(5.4) that the sequence L Sn (t) is bounded in Thus we get

Boundary variations technique for shape sensitivity analysis of work functional
Beside the existence of an optimal obstacle for the work, lift and drag shape optimization problems [6], it is important for applications to provide necessary optimality conditions and to devise a numerical method for solution of the shape optimization problems under considerations. The numerical methods of gradient or steepest descent types require the local information on the behavior of the shape functional to be minimized. The precise information on the shape gradient of cost functional can be obtained as a result from the appropriate shape sensitivity analysis of the functional. The shape sensitivity analysis requires some regularity of solutions to the governing equations like the Lipschitz continuity with respect to the boundary perturbations of the obstacle. The shape sensitivity analysis is performed in [6] for local solutions defined by small perturbations of a class of approximate solutions to the stationary problem.
6.1. Boundary and distributed shape functionals. We recall here, that the following notation is used for the shape functionals under considerations for the shape sensitivity analysis.
We consider the integral shape functionals denoted by J(Ω), with S B a compact obstacle in a hold all domain B, and Ω := B \ S. In the paper the symbols J(S) ≡ J(Ω) are used for the work functional. There is however a difference between J(S) and J(Ω), the work functional J(S) contains an integral over the lateral boundary ∂S × (0, T ), and the same work functional J(Ω) contains an integral over the cylinder Ω × (0, T ). Usually, the functional J(Ω) is obtained from the functional J(S) by an integration by parts formulae.
It is clear that the distributed shape functionals J(Ω) require less regularity from the solutions to the governing equations compared to the boundary shape functionals J(S). On the other hand the distributed shape functionals formally depend on a choice of a function denoted by η which is required in the integration by parts formulae, however in view of the identity J(S) ≡ J(Ω) the values of the shape functional J(Ω) are independent of the choice of η.

Shape sensitivity analysis within boundary variations technique.
Our goal now is to develop the shape sensitivity analysis which results in the shape derivatives of solutions to the governing equations and in the shape gradients of J(Ω) obtained for stationary and nonstationary governing equations by introduction of appropriate adjoint state equations.
Two different types of velocity fields can be employed. The physical field is the state variable u := u(Ω), Ω = B \ S determined from the governing equations for a given obstacle S. This field is in general nonunique, thus the local classical solutions of governing equations are considered for the purposes of the shape sensitivity analysis. The artificial velocity field V := V(ε, x), x ∈ B, is introduced for the purposes of the shape sensitivity analysis with respect to the small perturbations of the obstacle boundary in the normal direction. This field depends on the small shape parameter ε → 0 and it is associated with the domain transformation mapping T ε : S → S ε , Now, the form of the mapping T ε is specified, where the field T(x), x ∈ B, is compactly supported in a small neighborhood of the obstacle S and the support of T is disjoint with the boundary Σ = ∂B. This means that the boundary Σ is invariant under transformation (6.2). In order to evaluate the shape gradient of the functional J(Ω) the method of boundary variations is applied and the Eulerian semiderivative of the shape functional dJ(Ω; V) is obtained in the direction of a vector field V associated with the change of the variables T ε .
This means that for the mapping (6.2) the family of perturbed obstacles is defined by S ε := T ε (S), where ε → 0 stands for the shape parameter. As a result, the differentiability of the real valued function ε → J(Ω ε ), with Ω ε = B \ S ε , is considered at ε = 0, and the existence of the derivative is established.

Shape sensitivity analysis of Navier-Stokes equations
In order to perform the shape sensitivity analysis of the work functional for the nonstationary equations, first, the framework is established. In the governing equations, most of physical constants are posed to be equal to one, therefore the only constant is λ > 0. It is also assumed at this stage of analysis that there is no intersection between the inlet and the outlet on the boundary of the hold all domain B. The boundary variations technique is applied in order to investigate the dependence of the shape functional J(Ω ε ) on the shape of the obstacle S ε in the variable domain Ω ε = B \ S ε for ε → 0.
The tools we are going to discuss in this section include the shape derivatives of the solutions to the nonstationary, compressible Navier-Stokes equations, the shape gradient of the work functional J(Ω) and its decompositions into the geometrical and dynamical parts, and the adjoint state equations associated with the dynamical part of the shape gradient. The proofs of the results are given in [6] in the case of local solutions to the stationary, compressible Navier-Stokes equations. 7.1. Navier-Stokes state equations. Let us consider the general model for the purposes of shpe sensitivity analysis.
The state equation defined in the reference domain Ω × (0, T ) takes the form in Ω .

(7.1c)
Remark 9. For numerical methods [14,3] It is convenient to introduce the effective viscous pressure and rewrite the state equation in the equivalent form, (7.2d)

Linearized and adjoint state equations. Material and shape derivatives of solutions.
Material and shape derivatives of solutions to the governing equations are given by solutions to the appropriate linearized equations. We are going to derive the equations for the shape derivatives of solutions to Navier-Stokes equations.
For the sake of simplicity we assume that the intersection of the inlet Σ in with the outlet Σ out is empty. We assume also that the primal variables π and v for the density and the velocity in the linearized equations vanish for t = 0 and on Σ in and ∂Ω, respectively. The only exception from the homogeneous initial and boundary conditions is the nonhomegeneous Dirichlet condition for the shape derivative u of the velocity field on the obstacle boundary ∂S. The dual variables denoted by ϕ and φ for the density and the velocity vanish on Σ out and on ∂Ω, respectively.
In order to differentiate the solutions of the state equations with respect to the shape the linearized and the adjoint equations are introduced. To this end it is convenient to rewrite equations (7.1a) and (7.1b) in the following form where S(u) = (∇u + ∇u + (λ − 1) div uI) .
We multiply (7.3) and (7.4) by the smooth test functions φ, ϕ, respectively, φ with compact support in Ω and ϕ which vanishes on Σ out , and integrate by parts. It follows that and for the mass balance equation Now, denote by (π, v) a solution to the linearized system at the sufficiently smooth solution ( , u) of the nonlinear system, the same linearized system is derived for the so-called shape derivatives ( , u ), hence we obtain the integral identities satisfied for all test functions φ, It is also useful to perform the integration by parts for the test functions v which vanish on ∂B and are nonnull on the obstacle boundary ∂S, this leads to We denote by the boundary integrals on ∂S which furnish one part of the dynamical shape gradient. It is clear that L dyn,1 (v) depends only on the trace of v on the obstacle boundary ∂S, and the expression is nontrivial when used with the shape derivative v := u . There are the nonhomogeneous Dirichlet conditions on ∂S for the shape derivative of the velocity field u = − ∂u ∂n (T · n), such conditions result from the homogeneous Dirichlet condition for the velocity field u = 0 prescribed on on ∂S.
In the same way the linearized equation is obtained for the mass balance equation, for all test functions ϕ. We introduce the following notation for the bilinear forms defined for linearized operators, acting on the smooth functions such that, v = 0 on ∂Ω, and π = 0 on Σ in , the above expression can be slightly simplified assuming in addition that the initial values for t = 0 also vanish, π(0) = 0 and v(0) = 0.
We denote also Now we take the sum of bilinear form and decompose in the following way in order to identify the adjoint operators Assuming that for t = 0, we have π(0) = 0, v(0) = 0, and that the values of ϕ(T ) and φ(T ) are prescribed, in view of decomposition (7.7) we define the following adjoint operators, first for π := 0 in (7.7), then for v := 0 in (7.7),

Decomposition of shape gradient
The decomposition of the shape gradient into the geometrical and dynamical components seems to be useful for the numerical methods of shape optimization.
where η is a smooth function with a compact support which contains the obstacle boundary ∂S, and η(x) ≡ 1 in a vicinity of the obstacle boundary, furthermore T = ∇u + (∇u) + (λ − 1) div uI − p( ) I. Recall that if an obstacle S moves in atmosphere like a solid body then its physical position S t at time t is defined by S t = U(t)S + a(t). (8.2) Here a unitary matrix U and a vector field a(t) can be defined by an appropriate flight planning scenario. In this case the vector field W is given by formulae and we have Cu = rot W × u. It is convenient for our purposes to introduce the following notation for the shape functional

State equation and shape functional in variable domain.
We are going to apply the boundary variations technique [15], however in the framework adapted to our problem. The reference domain Ω = B \ S is transformed onto the perturbed domain Ω ε = B \ S ε , ε → 0, by means of the change of variables y = x + εT(x), x ∈ Ω, y ∈ Ω ε , (8.7) with the appropriate vector field T supported in a neighborhood of the obstacle S, since we are only interested in the boundary variations of the obstacle S. Therefore, the state equation in Ω ε takes the same form (7.1) with Ω replaced by Ω ε , with the new unknown functions denoted by (¯ ε ,ū ε ), and withq ε := p(¯ ε ) − λ divū ε , the new unkown functions depend on the small parameter ε → 0 which is omitted in the equations, The expression for the shape functional in Ω ε := B \ S ε takes the form where the functions 0 (x), f (x, t), x ∈ Ω ε , t ∈ (0, T ), are given by the restrictions to Ω ε of the functions which are defined for x ∈ R d , the plan of the flight for the deformed obstacle S ε takes the form of the matrix function given by (8.3), now W depends on x ∈ Ω ε . The functions ε , u ε , T ε are given by the solutions to state equations, therefore the functions depend implicitely on ε → 0.
8.3. Shape derivatives of solutions. We are interested in the form of the derivative for the mapping ε → J(Ω ε ).
From the general theory of shape optimization [15] it follows by the Hadamard structure Theorem of shape gradient that the first order derivative of this mapping, under appropriate regularity assumptions on the domain and on the solutions is of the form where the shape gradient G is in general given by a distribution which lives on the boundary. In [6,13] it is shown that for the drag functional the shape gradient is given by a function G S . Therefore, G is the so-called shape gradient of the shape functional J(Ω) in the direction of the vector field T. There are two distinct parts of the shape gradient, the geometrical part, and the dynamical part, it is shown in [13] that the geometrical part of the shape gradient vanishes in the case of the drag functional.
Actually, the shape gradient G can be decomposed into two parts, one which is easy to evaluate numerically which we call the geometrical part, and another which is very difficult to evaluate by numerical methods which is called the dynamical part. The reason is that in order to evaluate the dynamical part, it is required to solve not only the state equation, but also the so-called adjoint state equation which is a linearized variant of the state equations depending on the right hand sides on the derivative of the shape functional. This type of decomposition is easy for the linear problems which are well posed, and extremaly difficult for the nonlinear problems we consider in the monograph. Now, we explain briefly such a decomposition as well as the concepts of the material and shape derivatives for the solutions of our state equation. In the description it is assumed that the solutions are sufficiently smooth which in a specific application should be justified.
The dynamical part of the shape gradient contains the so-called shape derivatives of the fields , u. The remaining part of the shape gradient is called the geometrical part. Roughly speaking, the dynamical part of the shape gradient takes into account only the variations of the state with respect to the boundary variations of the obstacle. In other words, the geometrical part of the shape gradient is obtained for the solutions of the state equation replaced by restrictions to the variable domain of given functions defined e.g., in all the hold all domain, which means that the shape derivatives of the density and of the velocity are set to be null. The decomposition into two parts of the shape gradient is obtained at the final stage of our procedure.
In order to deduce the shape gradient G some methods are available. One possibility is the change of variables Ω ε y(x) = x + εT(x) → x ∈ Ω in the state equation, and derivation in the reference domain Ω with respect to ε → 0. Now, we give some details on the change of variables. It is convenient to change also the unknown velocity field using the Piola transformation in the following way where we denote Remark 10. Notice that the functions ε (x, t) and u ε (x, t) are only defined for x ∈ Ω ε , where ε → 0, therefore the limits (8.17) and (8.18) are well defined in the open set Ω.
Remark 11. Formally, the equations for the shape derivatives ( , u ) are obtained by the linearization of the state equation at the reference domain Ω × (0, T ), and the formal system is of the following form where we assume that the data in the state equation in the reference domain f , U, ∞ admit the shape derivatives C , f , U , ∞ , u 0 , 0 .
8.5. Adjoint state equations. The dynamical part of the shape derivative is further simplified by introduction of an appropriate adjoint state equations. We introduce two linear forms in order to decompose the dynamical part of the shape gradient. The first linear form for the density shape derivative, where the bilinear forms are defined by (7.8) and (7.9). The smooth test functions satisfy the following boundary conditions, π = 0 on Σ in × (0, T ), v = 0 on ∂Ω × (0, T ), v(0) = 0 and π(0) = 0 in Ω. Now, let us note that by the adjoint state equations we have the identity L dens ( ) + L vel (u ) = L π (ϕ, φ), and by the linearized equations for the shape derivatives it follows that we have the second identity L π (ϕ, φ), + L v (ϕ, φ), u =