About this Journal Submit a Manuscript Table of Contents
Mathematical Problems in Engineering
Volume 2013 (2013), Article ID 728534, 11 pages
Research Article

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

1School of Mathematics and Statistics, Xi'an Jiaotong University, Shaanxi 710049, China
2Center for Computational Geosciences, Xi'an Jiaotong University, Shaanxi 710049, China
3School of Science, Chang'an University, Shaanxi 710064, China

Received 20 November 2012; Revised 4 March 2013; Accepted 7 March 2013

Academic Editor: Jun Jiang

Copyright © 2013 Wenjing Yan et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


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.

1. 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 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.

2. 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 , is the inflow boundary, denotes the outflow boundary, represents the boundary corresponding to the fluid wall, and is the boundary which is to be optimized. The unknowns of the full model are the fluid velocity , the pressure , and the temperature .

Find satisfying where stands for the inverse of Reynolds number, the stress tensor is defined by with the rate of deformation tensor . represents the transpose of the matrix , denotes the identity tensor, is the unit normal vector on the boundary , and denotes the inverse of Peclet number.

We introduce the following functional spaces which will be used throughout this paper. Let be the space of square integrable real-valued functions on with the usual norm. The space , where denotes the standard Sobolev space on (see [9]), that is, the space of functions with generalized derivatives of order up to in . 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 is fixed. An example of the admissible set is Let be of piecewise , the minimization problem (4) has at least one solution with given area in two dimensions [10].

3. 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 gradient-type 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 with the boundary which belongs to piecewise , and a velocity space , where is a small positive real number and denotes the space of all -times continuous differentiable functions with compact support contained in . The velocity field belongs to for each . It can generate transformations through the following dynamical system with the initial value . We denote the transformed domain by at , and also set . There exists an interval , , and a one-to-one map from onto such that (see [6]) (1)(2) belongs to with ; (3) belongs to . In addition, for sufficiently small , the Jacobian of is strictly positive: where denotes the Jacobian matrix of the transformation evaluated at a point associated with the velocity field . We also introduce the following notations: represents the inverse of the matrix , and is the transpose of the matrix .

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 if the limit exists:
Furthermore, if the map is linear and continuous, we say that is shape differentiable at . In the distributional sense, we have where is the shape gradient of in .

4. 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 , 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 and , the cost functional possesses the shape gradient which can be expressed as where the adjoint states 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 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 Then, we differentiate with respect to in the direction and apply Green’s formula:Considering with compact support in , we obtain

Then, varying on gives

Similarly, we differentiate with respect to in the direction :

According to the compact support of in , we can get

In addition, varying on 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 ; then the flow of the velocity field can be expressed as follows:The perturbed domain can be defined by . Now, we define the cost functional: where and 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 , , , 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 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 , we introduce the following Hadamard formula: for a sufficiently smooth functional .

By Hadamard formula (34), we obtain where

To simplify (36)–(38), we introduce the following lemma.

Lemma 3 (see [6]). If two vector functions and vanish on the boundary , the following identities hold on the boundary .

Noting and vanish on the boundary , , and employing Lemma 3, (36) can be rewritten as Since and satisfy (1) and (14), respectively, (37) reduces to Similarly, according to Lemma 3, the state system (2), and the adjoint state system (15), (38) becomes Adding (40)–(42) together, we finally obtain the boundary expression for the Eulerian derivative of : Since the mapping is linear and continuous, we get the expression for the shape gradient by (11). This completes the proof.

5. Finite-Element Approximation and Numerical Simulation

5.1. Discretization of the Optimization Problem

We suppose that is a bounded polygonal domain of and only consider the conforming finite-element approximations. Let and 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 such that for , (H2) there exists such that for , (H3) the Ladyzhenskaya-Brezzi-Babuska inf-sup condition is verified; that is, there exists such that 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

5.2. 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 and is a positive Lagrange multiplier. The Eulerian derivative of is and the shape gradient . Ignoring regularization, a descent direction is found by defining , and then we can update the shape as , 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, . In this case, the descent direction is the unique element of the problem

The resulting algorithm consists of the following parts: (1)give an initial shape , an initial step , and a Lagrange multiplier ; (2)solve the state system and adjoint state system, then we can evaluate the descent direction by (55) with and ; (3)set , 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.

5.2.1. 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 of our algorithm, we refresh the Lagrange multiplier by where is a small positive parameter and 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 so that

5.2.3. 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.

5.3. Numerical Examples

In all computations, the finite-element discretization is effected using the bubble— 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 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 projection: where denotes the projection on the Lagrange finite-element 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 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.

Figure 1: Schematic illustration for the boundary of domain .

We choose to be a rectangle 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 We set the initial shape of the body to be a circle of center with radius .

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 have been shown in Figure 2, and the initial finite-element mesh consists of 1676 elements with 922 vertices.

Figure 2: Initial mesh of the domain .

Figures 3 and 4 demonstrate the comparison between the initial shape and optimal shape for the distribution of the velocity , the pressure , and the temperature with the Reynolds numbers Re = 100.

Figure 3: Comparison of the initial shape and optimal shape ().
Figure 4: Comparison of the initial shape and optimal shape ().

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.

Figure 5: Convergence history of the cost function ().
Figure 6: Convergence history of the cost function ().

6. 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.


This work is supported by National Natural Science Foundation of China (Grants nos.10901127 and 11001216), and the Fundamental Research Funds for the Central Universities.


  1. O. Pironneau, Optimal Shape Design for Elliptic Systems, Springer Series in Computational Physics, Springer, Berlin, Germany, 1984. View at Zentralblatt MATH · View at MathSciNet
  2. B. Mohammadi and O. Pironneau, Applied Shape Optimization for Fluids, The Clarendon Press, New York, NY, USA, 2001. View at Zentralblatt MATH · View at MathSciNet
  3. J. Simon, “Domain variation for drag in Stokes flow,” in Control Theory of Distributed Parameter Systems and Applications, vol. 159 of Lecture Notes in Control and Information Sciences, pp. 28–42, Springer, Berlin, Germany, 1991. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  4. J. A. Bello, E. Fernández Cara, and J. Simon, “The variation of the drag with respect to the domain in Navier-Stokes flow,” in Optimization, Optimal Control, Partial Differential Equations, vol. 107 of International Series of Numerical Mathematics, pp. 287–296, Birkhäuser, Basel, Switzerland, 1992. View at Zentralblatt MATH · View at MathSciNet
  5. J. A. Bello, E. Fernández-Cara, J. Lemoine, and J. Simon, “The differentiability of the drag with respect to the variations of a Lipschitz domain in a Navier-Stokes flow,” SIAM Journal on Control and Optimization, vol. 35, no. 2, pp. 626–640, 1997. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  6. M. C. Delfour and J.-P. Zolésio, Shapes and Geometries: Analysis, Differential Calculus, and Optimization, vol. 4 of Advances in Design and Control, SIAM, Philadelphia, Pa, USA, 2001. View at MathSciNet
  7. W.-J. Yan and Y.-C. Ma, “Shape reconstruction of an inverse Stokes problem,” Journal of Computational and Applied Mathematics, vol. 216, no. 2, pp. 554–562, 2008. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  8. D. Chenais, J. Monnier, and J. P. Vila, “Shape optimal design problem with convective and radiative heat transfer: analysis and implementation,” Journal of Optimization Theory and Applications, vol. 110, no. 1, pp. 75–117, 2001. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  9. R. A. Adams, Sobolev Spaces, vol. 65, Academic Press, New York, NY, USA, 1975. View at MathSciNet
  10. O. Pironneau, “Optimal shape design by local boundary variations,” in Optimal Shape Design, vol. 1740 of Lecture Notes in Mathematics, pp. 343–384, Springer, Berlin, Germany, 2000. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  11. G. Dogan, P. Morin, R. H. Nochetto, and M. Verani, “Discrete gradient flows for shape optimization and applications,” Computer Methods in Applied Mechanics and Engineering, vol. 196, no. 37–40, pp. 3898–3914, 2007. View at Publisher · View at Google Scholar · View at MathSciNet
  12. J. Simon, “Differentiation with respect to the domain in boundary value problems,” Numerical Functional Analysis and Optimization, vol. 2, no. 7-8, pp. 649–687, 1980. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  13. A. Quarteroni, R. Sacco, and F. Saleri, Numerical Mathematics, vol. 37 of Texts in Applied Mathematics, Springer, Berlin, Germany, 2nd edition, 2007. View at MathSciNet
  14. O. C. Zienkiewicz and J. Z. Zhu, “The superconvergent patch recovery and a posteriori error estimates. I. The recovery technique,” International Journal for Numerical Methods in Engineering, vol. 33, no. 7, pp. 1331–1364, 1992. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  15. J. Dompierre, P. Labb, and F. Guibault, “Controlling approximation error,” in Proceedings of the 2nd MIT Conference on Computational Fluid and Solid Mechanics, pp. 1929–1932, Elsevier, Cambridge, Mass, USA, 2003.
  16. F. Hecht, “A few snags in mesh adaptation loops,” in Proceedings of the 14th International Meshing Roundtable, pp. 301–302, Springer, San Diego, Calif, USA, 2005.
  17. R. Temam, Navier-Stokes Equations, Theory and Numerical Analysis, AMS Chelsea, Providence, RI, USA, 2001. View at MathSciNet
  18. V. Girault and P.-A. Raviart, Finite Element Methods for Navier-Stokes Equations: Theory and Algorithms, vol. 5 of Springer Series in Computational Mathematics, Springer, Berlin, Germany, 1986. View at Publisher · View at Google Scholar · View at MathSciNet