The Application of Shape Gradient for the Incompressible Fluid in Shape Optimization

This paper is concerned with the numerical simulation for shape optimization of the Stokes flow around a solid body. The shape gradient for the shape optimization problem in a viscous incompressible flow is evaluated by the velocity method. The flow is governed by the steady-state Stokes equations coupled with a thermal model. The structure of continuous shape gradient of the cost functional is derived by employing the differentiability of a minimax formulation involving a Lagrange functional with the function space parametrization technique. A gradient-type algorithm is applied to the shape optimization problem. Numerical examples show that our theory is useful for practical purpose, and the proposed algorithm is feasible and effective.


Introduction
The problem of finding the shape design of a system described by the incompressible Stokes equations arises in many design problems such as aerospace, automotive, hydraulic, ocean, and structural and wind engineering.The problem is to optimize the shape of the domain in order to minimize a cost functional that depends on the solutions.In the numerical implements, the shape optimization problem can be transformed into the minimization problem without constraint condition by the Lagrange multiplier and the adjoint equations using adjoint variables corresponding to the state equations.
The optimal shape design of a body subjected to the minimum viscous dissipated energy has been a challenging task for a long time, and it has been investigated by several authors.For instance, Pironneau in [1,2] computed the derivative of the cost functional using normal variation approach; Simon [3] used the formal calculus to deduce an expression for the derivative; Bello et al. in [4,5] considered this problem theoretically in the case of Navier-Stokes flow by the formal calculus.
In the present paper, we will use the so-called function space parametrization technique which was advocated by Delfour and Zolésio to solve Poisson equation with Dirichlet and Neumann condition (see [6]).In our paper [7], we solved a shape reconstruction problem for the inverse Stokes problem and investigated the numerical simulation by the domain derivation and the regularized Gauss-Newton iterative method.However, many authors studied optimal shape design problems in fluid without heat transfer, steady state or not.Chenais et al. studied a shape optimal design problem in a potential flow coupled with a thermal model in [8].
In this paper, we will study the energy minimization problem for Stokes flow with convective heat transfer in spite of its lack of rigorous mathematical justification in case the Lagrange formulation is not convex.We shall show how this theorem allows, at least formally, to bypass the study of material derivative and obtain the expression of shape gradient for the dissipated energy functional.For the numerical solution of the viscous energy minimization problem, we introduce a gradient-type algorithm with mesh adaptation technique, while the partial differential systems are discretized by means of the finite-element method.Finally, we give some numerical examples concerning the optimization of a two-dimensional obstacle located in the viscous flow.
This paper is organized into five parts.In Section 2, we introduce the shape optimization problem of Stokes flow Mathematical Problems in Engineering coupled with conductive heat transfer.The cost functional which we considered is general and depends on the solutions.In Section 3, we briefly recall the velocity method which is used for the characterization of the deformation of the shape.Section 4 is devoted to the computation of the shape gradient of the Lagrange functional due to a minimax principle concerning the differentiability of the minimax formulation by the function space parametrization technique.Finally in Section 5, we propose a gradient-type algorithm for the shape optimization problem, and the examples verify that our method is efficient and useful in the numerical implementations.

Statement of the Shape Optimization Problem
In two dimensions, we consider a typical problem in which a solid body  with the boundary Γ  is located in an external flow.Let Γ denote the boundary of Ω, and suppose that an incompressible viscous flow occupies Ω.The boundary where ] stands for the inverse of Reynolds number, the stress tensor is defined by (u, ) := −I + 2](u) with the rate of deformation tensor (u) := (Du+ * Du)/2.* Du represents the transpose of the matrix Du, I denotes the identity tensor, n is the unit normal vector on the boundary Γ out , and  denotes the inverse of Peclet number.We introduce the following functional spaces which will be used throughout this paper.Let  2 (Ω) be the space of square integrable real-valued functions on Ω with the usual norm.The space   (Ω), where  = 1, 2, . . .,  denotes the standard Sobolev space on Ω (see [9]), that is, the space of functions with generalized derivatives of order up to  in  2 (Ω).Consider Our aim is to optimize the shape of the boundary Γ  that minimizes a given cost functional  which depends on the velocity and the temperature.The cost functional may represent a given objective related to specific characteristic features of the fluid flow.We consider the following minimization problem for an incompressible viscous Stokes flow with convective heat exchanges: where the boundary Let Ω be of piecewise  1 , the minimization problem (4) has at least one solution with given area in two dimensions [10].

The Velocity Method
The mathematical difficulty of this problem is as follows.On the one hand, the set of domains Ω is not a vectorial space; on other hand, we need an expression of the differential of the cost functional in order to use a gradient-type algorithm.
To avoid this difficulty, we recall the definition of admissible domains and define the derivative of a real-valued function with respect to the domain in a classical manner.Then, we are able to give an expression of the differential of the cost functional with the intent to construct a gradienttype algorithm.There are about three types of techniques to perform the domain deformation: Hadamard's normal variation method, the perturbation of the identity method [11], and the velocity method [6].We will use the velocity method which contains the others.At first, we choose an open set Ω in R 2 with the boundary Ω which belongs to piecewise   , and a velocity space V ∈ E  := {V ∈ ([0, ]; D  (Ω, R 2 ))}, where  is a small positive real number and D  (Ω, R 2 ) denotes the space of all -times continuous differentiable functions with compact support contained in R 2 .The velocity field belongs to D  (Ω, R 2 ) for each .It can generate transformations   (V)  =  (, ) ,  ≥ 0,  ∈ Ω (7) through the following dynamical system d d (, ) = V (,  ()) , (0, ) =  (8) with the initial value .We denote the transformed domain   (V)(Ω) by Ω  (V) at  ≥ 0, and also set Γ  :=   (Γ).There exists an interval  = [0, ), 0 <  ≤ , and a one-to-one map   from Ω onto Ω such that (see [6]) (1)  0 = I; (2) (, )  →   () belongs to  1 (;   (Ω)) with   (Ω) = Ω; (3) (, )  →  −1  () belongs to (;   (Ω)).In addition, for sufficiently small  > 0, the Jacobian   of   is strictly positive: where D  () denotes the Jacobian matrix of the transformation   evaluated at a point  ∈ Ω associated with the velocity field V. We also introduce the following notations: D −1  () represents the inverse of the matrix D  (), and * D −1  () is the transpose of the matrix D −1  ().
Definition 1.Let (Ω) be a real-valued functional associated with any regular domain Ω, we call that the functional has an Eulerian derivative at Ω in the direction V if the limit exists: Furthermore, if the map V  → d(Ω; V) : E  → R is linear and continuous, we say that  is shape differentiable at Ω.In the distributional sense, we have where ∇ is the shape gradient of  in Ω.

Function Space Parametrization
In order to compute the exact differential or the shape gradient, few approaches are possible.In the direct differentiation, it requires to derive the state equations with respect to the shape variables.In practice, it implies to solve as many PDE systems as discrete shape variables.In order to avoid this extra computational cost, we use the classical adjoint state method which requires to solve only one extra PDE system.There are two ways for it.The first one is to discretize the equations, using a finite-element method for example, and to derive the discrete equations and obtain the discrete shape gradient.The second one is to calculate the expression of the exact differential of the cost functional and to discretize it.In this paper, we follow the last approach.We will introduce the adjoint state equation and obtain an expression of the exact differential of the cost functional (Ω).
Let Ω be of class  2 , and we derive the variational formulation of the state system (1) and (2) in appropriate Sobolev spaces: Now, we will prove the main theoretical results of the paper.

Theorem 2.
Let Ω be of class  2 and V ∈ E 2 , the cost functional (Ω) possesses the shape gradient ∇ which can be expressed as where the adjoint states k and  satisfy the following adjoint system: Proof.We will utilize the differentiability of a minimax formulation involving a Lagrangian functional with a function space parametrization technique.Now, we introduce the following Lagrangian functional associated with (12): where Problem ( 16) can be put in the following form: The minimax framework can be used to avoid the study of the state derivative with respect to the shape of the domain.
The Karush-Kuhn-Tucker conditions will furnish the shape gradient of the cost functional (Ω) by using the adjoint system.
The first optimality condition for the problem can be established as follows: Formally the adjoint equations are defined from the Euler-Lagrange equations of the Lagrange functional .Clearly, the variation of  with respect to (k, , ) can recover the state system and its mixed weak formulation (12).In order to find the adjoint state system, we differentiate  with respect to  in the direction : Taking  with compact support in Ω gives div k = 0.
Then, we differentiate  with respect to u in the direction u and apply Green's formula: Considering u with compact support in Ω, we obtain Then, varying u on Γ out gives Similarly, we differentiate  with respect to  in the direction : According to the compact support of u in Ω, we can get In addition, varying u on Γ out and Γ  leads to Hence, we obtain the following adjoint state system ( 14) and (15).The velocity method is employed to simulate the domain deformations.We perturb the boundary Γ  and consider the transformation   (V); then the flow of the velocity field can be expressed as follows: V ∈  ad := {V ∈  0 (0, ; ( 2 (R 2 )) 2 ) :  = 0 in the neighborhood of Γ  } . (28) The perturbed domain can be defined by Ω  =   (V)(Ω).Now, we define the cost functional: where (u  ,   ,   ) and (k  ,   ,   ) satisfy ( 14) and ( 15) on the perturbed domain Ω  , respectively.Our object is to derive the derivative with respect to the variable .Unfortunately, the Sobolev spaces   (Ω  ),  0 (Ω  ), (Ω  ), and (Ω  ) depend on the parameter , so we need to introduce the so-called function space parametrization technique which consists in transporting the different quantities defined on the variable domain Ω  back into the reference domain Ω which does not depend on the perturbation parameter .Since the functionals above mentioned are defined in a fixed domain Ω with respect to the parameter , we can apply the differential calculus.Hence, we define the following functional spaces: where "∘" denotes the composition of the two maps.
Since   and  −1  are diffeomorphisms, these parametrizations cannot change the value of the saddle point.We can rewrite (29) as where the Lagrangian functional with To perform the differentiation of the perturbed Lagrangian functional (Ω  , u∘ −1  , ∘ −1  , ∘ −1  , k∘ −1  , ∘ −1  , ∘ −1  ), we introduce the following Hadamard formula: for a sufficiently smooth functional   : [0, ] × R 2 → R. By Hadamard formula (34), we obtain Adding ( 40)-( 42) together, we finally obtain the boundary expression for the Eulerian derivative of (Ω): Since the mapping V  → d(Ω; V) is linear and continuous, we get the expression for the shape gradient by (11).This completes the proof.

Discretization of the Optimization Problem.
We suppose that Ω is a bounded polygonal domain of R 2 and only consider the conforming finite-element approximations.Let  ℎ ⊂ ( 2 (Ω))  and S ℎ ⊂  2 (Ω) be two families of finite dimensional subspaces parameterized by ℎ which tends to zero.We also define the following functional spaces: Besides, the following assumptions are supposed to hold.
(H1) There exists  > 0 such that for 0 ≤  ≤ , inf (47) (H3) the Ladyzhenskaya-Brezzi-Babuska inf-sup condition is verified; that is, there exists  > 0 such that inf The Galerkin finite-element approximation of the state system ( 1) and ( 2) in mixed form are as follows: Similarly, we can obtain the finite element approximation for adjoint state system, We also have the discrete cost functional and the discrete shape gradient

A Gradient-Type Algorithm.
Next, we will present a gradient-type algorithm and numerical examples in two dimensions to verify that our previous methods could be very useful and efficient for the numerical implementation of the shape optimal design problem.We describe the gradient-type algorithm for the minimization of a cost functional (Ω).For the minimization problem (4), we rather work with the unconstrained minimization problem where (Ω) := ∫ Ω d and  is a positive Lagrange multiplier.The Eulerian derivative of (Ω) is and the shape gradient Ignoring regularization, a descent direction is found by defining V = −ℎ  ∇, and then we can update the shape Ω as Ω  = (I + ℎ  V)Ω, where ℎ  is a descent step at th iteration.However, in this paper in order to avoid boundary oscillations (and irregular shapes) and due to the fact that the gradient-type algorithm produces shape variations which have less regularity than the original parametrization, we change the scalar product with respect to which we compute a descent direction [2,12], for instance, ( 1 (Ω)) 2 .In this case,  the descent direction is the unique element d ∈ ( 1 (Ω)) 2 of the problem The resulting algorithm consists of the following parts: (1) give an initial shape Ω 0 , an initial step ℎ 0 , and a Lagrange multiplier  0 ; (2) solve the state system and adjoint state system, then we can evaluate the descent direction   by (55) with Ω = Ω  and  =   ; (3) set Ω +1 = ( − ℎ    )Ω  , where ℎ  is a small positive real number and can be chosen by some rules (see [1]).
Next, we will discuss some details of the gradient-type algorithm and they will make our algorithm truly efficient and effective.

Lagrange Multiplier.
We now discuss the choice of the Lagrange multiplier  in the optimization problem (53).The value of  is updated at each iteration so the shape satisfies the fixed-volume constraint when the algorithm converges.Due to the relatively high cost in moving the mesh, we do not impose exactly the volume constraint before convergence.If the present volume is greater than the target volume, we increase the multiplier , otherwise we decrease it.However, this may lead to the oscillation of the volume.Hence in this  paper, we relax it by assuming that the first order optimality condition is satisfied: at least in the average sense on the boundary Γ  ; that is, Finally in Step (3) of our algorithm, we refresh the Lagrange multiplier by where  is a small positive parameter and  target (Ω) denotes the target volume of the shape.

5.2.2.
Step Size.The choice of the descent step size ℎ  is not an easy task.Too big, the algorithm is unstable; too small, the rate of convergence is insignificant.The classical exact line search method can be very expensive and is unnecessary to guarantee convergence in shape optimization problems.
Here, we use the backtracking approach [13].To limit the number of the required state solutions and to prevent the solver from crashing because of bad shape, it is important to provide the backtracking procedure with a good initial guess.
Here, we choose the initial guess ℎ 0 so that

Stopping Criterion.
In our algorithm, we do not choose any stopping criterion.A classical stopping criterion is to find that whether the shape gradients in some suitable norm are small enough.However, since we use the continuous shape gradients, it is hopeless for us to expect very small gradient norm because of numerical discretization errors.Instead, we fix the number of iterations.If it is too small, we can restart it with the previous final shape as the initial shape.

Numerical Examples.
In all computations, the finiteelement discretization is effected using the  1 bubble- 1 pair of finite element spaces on a triangular mesh.The mesh is performed by a Delaunay-Voronoi mesh generator (see [1]), and during the shape deformation, we utilize a metric-based anisotropic mesh adaptation technique where the metric can be computed automatically from the Hessian of a solution.The Hessian H ℎ of  ℎ can be approximated by using a recovery method, such as the Zienkiewicz-Zhu recovery procedure [14], the simple linear fitting [15], or the double  2 projection: where   2 denotes the  2 projection on the  1 Lagrange finiteelement space (see [16]).Here, we use (60) to get the Hessian.
As it has been said in [16], there is no convergence proof of this method but the result is better.We consider shape optimization of the Stokes flow around a solid body in two dimensions.The schematic geometry of the fluid domain is described in Figure 1, corresponding to an external flow around a solid body .We reduce the problem to a bounded domain  by introducing an artificial boundary  := Γ in ∪Γ  ∪Γ out which has to be taken sufficiently far from  so that the corresponding flow is a good approximation of the unbounded external flow around  and Ω :=  \  is the effective domain.In addition, the boundary Γ  :=  is to be optimized.
We choose  to be a rectangle (0, 9) × (0, 3) and  is to be determined in our simulations.The no-slip boundary conditions are imposed at all the other boundaries.The admissible set is defined as O := {Ω ⊂ R 2 :  is fixed, the area  target (Ω) = 2.1} . (61) We set the initial shape of the body  to be a circle of center (4.5, 1.5) with radius  = 1.
The state systems and the adjoint system are discretized by a mixed finite-element method.Spatial discretization is effected using the Taylor-Hood pair [17,18] of finite-element spaces on a triangular mesh; that is, the finite-element spaces are chosen to be continuous piecewise quadratic polynomials for the velocity and continuous piecewise linear polynomials for the pressure.
The finite-element meshes used for the calculations at Re = 100 have been shown in Figure 2, and the initial finiteelement mesh consists of 1676 elements with 922 vertices.
Figures 3 and 4 demonstrate the comparison between the initial shape and optimal shape for the distribution of the velocity u = ( 1 ,  2 )  , the pressure , and the temperature  with the Reynolds numbers Re = 100.
We run many iterations in order to show the good convergence and stability properties of our algorithm; however, it is clear that it has converged in a much smaller number of iterations (see Figures 5 and 6).We also find that, when the Reynolds number increases, the reduced energy increases with fixed numbers of iterations.However, as the Reynolds number increases, the computational cost associated with the computation of the Stokes system raises; hence, the cost of the fully optimization procedure increases.

Conclusion
In this paper, the optimal shape problem in a Stokes flow coupled with convective heat transfer has been presented.We derived the structure of shape gradient for the cost functional by function space parametrization technique without the usual study of the derivative of the state.Though for the time being this technique lacks a rigorous mathematical framework, a gradient-type algorithm is effectively used for the minimization problem for various Reynolds numbers and we also compare the results for different Reynolds numbers in numerical tests.Further research is necessary on efficient implementations for the Navier-Stokes flow with higher Reynolds number and much more real problems in the industry.

Figure 1 :
Figure 1: Schematic illustration for the boundary of domain Ω.

Figure 2 :
Figure 2: Initial mesh of the domain Ω.

( a )Figure 4 :
Figure 4: Comparison of the initial shape and optimal shape (Re = 100).