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

Abstract

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 Γ∶=Γ_𝑤∪Γ_𝑠∪Γ_in∪Γ_out, Γ_in is the inflow boundary, Γ_out 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 −𝜈Δ𝐮+∇𝑝=0inΩdiv𝐮=0inΩ𝐮=0onΓ_𝑤∪Γ_𝑠𝐮=𝑔onΓ_in𝜎𝐮,𝑝⋅𝐧=0onΓ_out,−𝛽Δ𝑇+𝐮⋅∇𝑇=0inΩ𝜕𝑇𝜕𝐧=0onΓ_out∪Γ_𝑠𝑇=𝑇₁onΓ_𝑤𝑇=𝑇₂onΓ_in, (1)(2) where 𝜈 stands for the inverse of Reynolds number, the stress tensor is defined by 𝜎(𝐮,𝑝)∶=−𝑝I+2𝜈𝜀(𝐮) with the rate of deformation tensor 𝜀(𝐮)∶=(D𝐮+^∗D𝐮)/2. ^∗D𝐮 represents the transpose of the matrix D𝐮, I denotes the identity tensor, 𝐧 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 𝐿²(Ω) 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 𝐿²(Ω). Consider𝑉₀(Ω)∶=𝐮∈𝐻²(Ω)²∶𝐮=0onΓ_𝑤∪Γ_𝑠∪Γ_in,𝑉_𝑔(Ω)∶=𝐮∈𝐻²(Ω)²∶𝐮=0onΓ_𝑤∪Γ_𝑠,𝐮=𝐠onΓ_in,𝑊(Ω)∶=𝑆∈𝐻¹(Ω)∶𝑆=0onΓ_𝑤∪Γ_in,𝑄(Ω)∶=𝑝∈𝐻¹(Ω)∶_Ω𝑝d𝑥=0. (3)

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: min_Ω∈𝒪𝐽(Ω)=2𝜈_Ω|𝜀(𝐮)|²1d𝑥+2_Ω|∇𝑇|²d𝑥, (4) where the boundary Γ_𝑤∪Γ_in∪Γ_out is fixed. An example of the admissible set 𝒪 is 𝒪∶=Ω⊂ℝ²∶Γ_𝑤∪Γ_in∪Γ_outisfixed,_Ωd𝑥=constant. (5) 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 𝐕∈E^𝑘∶={𝐕∈𝐶([0,𝛼];𝒟^𝑘(Ω,ℝ²))}, 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 𝐕(𝑡)(𝑥)=𝐕(𝑡,𝑥),𝑥∈Ω,𝑡≥0 (6) belongs to 𝒟^𝑘(Ω,ℝ²) for each 𝑡. It can generate transformations 𝐹_𝑡(𝐕)𝑋=𝑥(𝑡,𝑋),𝑡≥0,𝑋∈Ω (7) through the following dynamical system d𝑥()((𝑡()d𝑡𝑡,𝑋=𝐕𝑡,𝑥)),𝑥0,𝑋=𝑋 (8) with the initial value 𝑋. We denote the transformed domain 𝐹_𝑡(𝐕)(Ω) by Ω_𝑡(𝐕) at 𝑡≥0, and also set Γ_𝑡∶=𝐹_𝑡(Γ). There exists an interval 𝐼=[0,𝛽), 0<𝛽≤𝛼, and a one-to-one map 𝐹_𝑡 from Ω onto Ω such that (see [6]) (1)𝐹₀=I;(2)(𝑡,𝑥)↦𝐹_𝑡(𝑥) belongs to 𝐶¹(𝐼;𝐶^𝑘(Ω)) with 𝐹_𝑡(𝜕Ω)=𝜕Ω; (3)(𝑡,𝑥)↦𝐹_𝑡⁻¹(𝑥) belongs to 𝐶(𝐼;𝐶^𝑘(Ω)). In addition, for sufficiently small 𝑡>0, the Jacobian 𝐽_𝑡 of 𝐹_𝑡 is strictly positive: 𝐽_𝑡(𝑥)∶=detD𝐹_𝑡(𝑥)=detD𝐹_𝑡(𝑥)>0, (9) where D𝐹_𝑡(𝑥) denotes the Jacobian matrix of the transformation 𝐹_𝑡 evaluated at a point 𝑥∈Ω associated with the velocity field 𝐕. We also introduce the following notations: D𝐹_𝑡⁻¹(𝑥) represents the inverse of the matrix D𝐹_𝑡(𝑥), and ^∗D𝐹_𝑡⁻¹(𝑥) is the transpose of the matrix D𝐹_𝑡⁻¹(𝑥).

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: lim_𝑡↘0𝐽Ω_𝑡−𝐽(Ω)𝑡∶=d𝐽(Ω;𝐕). (10)
Furthermore, if the map 𝐕↦d𝐽(Ω;𝐕)∶E^𝑘→ℝ is linear and continuous, we say that 𝐽 is shape differentiable at Ω. In the distributional sense, we have d𝐽(Ω;𝐕)=⟨∇𝐽,𝐕⟩_{𝒟𝑘(Ω,ℝ2)′×𝒟𝑘(Ω,ℝ2)}, (11) 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: seek𝐮,𝑝∈𝑉_𝑔(Ω)×𝑄(Ω),suchthat_Ω2𝜈𝜀(𝐮)∶𝜀(𝐯)−𝑝div𝐯d𝑥=0,∀𝐯∈𝑉₀(Ω),_Ωdiv𝐮𝑞d𝑥=0,∀𝑞∈𝑄(Ω),find𝑇∈𝐻¹(Ω),suchthat_Ω𝛽∇𝑇⋅∇𝑆+𝐮⋅∇𝑇𝑆d𝑥=0,∀𝑆∈𝑊(Ω). (12) Now, we will prove the main theoretical results of the paper.

Theorem 2. Let Ω be of class 𝐶² and 𝐕∈E², the cost functional 𝐽(Ω) possesses the shape gradient ∇𝐽 which can be expressed as ∇𝐽=2𝜈𝜀(𝐮)∶𝜀(𝐯)+|𝜀(𝐮)|²+12|∇𝑇|²+𝛽∇𝑇⋅∇𝑆𝐧, (13) where the adjoint states 𝐯 and 𝑆 satisfy the following adjoint system: 𝜎(𝐮)−𝜈Δ𝐯+∇𝑞+𝑆∇𝑇=−2𝜈Δ𝐮inΩdiv𝐯=0inΩ𝐯,𝑞⋅𝐧−4𝜈𝜀⋅𝐧=0,onΓ_out𝐯=0onΓ_in∪Γ_𝑤∪Γ_𝑠,𝛽𝛽Δ𝑆+𝐮⋅∇𝑆=Δ𝑇inΩ𝜕𝑆𝜕𝐧+𝐮⋅𝑆⋅𝐧=∇𝑇⋅𝐧onΓ_out∪Γ_𝑠𝑆=0onΓ_in∪Γ_𝑤. (14)(15)

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): 𝐺Ω,𝐮,𝑝,𝑇,𝐯,𝑞,𝑆=𝐽(Ω)−𝐿Ω,𝐮,𝑝,𝑇,𝐯,𝑞,𝑆, (16) where 𝐿=Ω,𝐮,𝑝,𝑇,𝐯,𝑞,𝑆_Ω(𝐮)(𝐯)+2𝜈𝜀∶𝜀−𝑝div𝐯−div𝐮𝑞d𝑥_Ω𝛽∇𝑇⋅∇𝑆+𝐮⋅∇𝑇𝑆d𝑥. (17) Problem (16) can be put in the following form: min_Ω∈𝒪min_{𝑢,𝑝,𝑇∈𝑉𝑔(Ω)(Ω)×𝑄×𝐻1(Ω)}max_{𝑣,𝑞,𝑆∈𝑉0(Ω)(Ω)(Ω)×𝑄×𝑊}𝐺Ω,𝐮,𝑝,𝑇,𝐯,𝑞,𝑆. (18)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: min_{𝑢,𝑝,𝑇∈𝑉𝑔(Ω)(Ω)×𝑄×𝐻1(Ω)}max_{𝑣,𝑞,𝑆∈𝑉0(Ω)(Ω)(Ω)×𝑄×𝑊}𝐺Ω,𝐮,𝑝,𝑇,𝐯,𝑞,𝑆. (19)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 𝛿𝑝: 𝜕𝐺𝜕𝑝Ω,𝐮,𝑝,𝑇,𝐯,𝑞,𝑆⋅𝛿𝑝=_Ω𝛿𝑝div𝐯d𝑥=0. (20) Taking 𝛿𝑝 with compact support in Ω gives div𝐯=0. (21)Then, we differentiate 𝐺 with respect to 𝐮 in the direction 𝛿𝐮 and apply Green’s formula:𝜕𝐺=𝜕𝐮Ω,𝐮,𝑝,𝑇,𝐯,𝑞,𝑆⋅𝛿𝐮_Ω−2𝜈Δ𝐮+𝜈Δ𝐯−∇𝑞⋅𝛿𝐮d𝑥−_Ω𝑆∇𝑇⋅𝛿𝐮d𝑥−2𝜈_𝜕Ω𝜀(𝐯)⋅𝐧⋅𝛿𝐮d𝑠+4𝜈_𝜕Ω𝜀(𝐮)+⋅𝐧⋅𝛿𝐮d𝑠_𝜕Ω=𝑞𝛿𝐮⋅𝐧d𝑠_Ω−2𝜈Δ𝐮+𝜈Δ𝐯−∇𝑞⋅𝛿𝐮d𝑥−_Ω−𝑆∇𝑇⋅𝛿𝐮d𝑥_𝜕Ω𝜎𝐯,𝑞⋅𝐧⋅𝛿𝐮d𝑠+4𝜈_𝜕Ω𝜀(𝐮)⋅𝐧⋅𝛿𝐮d𝑠. (22)Considering 𝛿𝐮 with compact support in Ω, we obtain −𝜈Δ𝐯+∇𝑞+𝑆∇𝑇=−2𝜈Δ𝐮. (23)

Then, varying 𝛿𝐮 on Γ_out gives 𝜎𝐯,𝑞⋅𝐧−4𝜈𝜀(𝐮)⋅𝐧=0onΓ_out. (24)

Similarly, we differentiate 𝐺 with respect to 𝑇 in the direction 𝛿𝑇: 𝜕𝐺=𝜕𝑇Ω,𝐮,𝑝,𝑇,𝐯,𝑞,𝑆⋅𝛿𝑇_Ω−𝛽Δ𝑆+𝐮⋅∇𝑆⋅𝛿𝑇d𝑥_{Γout∪Γ𝑠}𝛽𝜕𝑆+𝜕𝐧+𝐮⋅𝑆⋅𝐧𝛿𝑇d𝑠_Ω()=∇𝑇⋅∇𝛿𝑇d𝑥_Ω−𝛽Δ𝑆+𝐮⋅∇𝑆⋅𝛿𝑇d𝑥_{Γout∪Γ𝑠}𝛽𝜕𝑆−𝜕𝐧+𝐮⋅𝑆⋅𝐧𝛿𝑇d𝑠_ΩΔ𝑇⋅𝛿𝑇d𝑥+_𝜕Ω∇𝑇⋅𝛿𝑇⋅𝐧d𝑠. (25)

According to the compact support of 𝛿𝐮 in Ω, we can get 𝛽Δ𝑆+𝐮⋅∇𝑆=Δ𝑇. (26)

In addition, varying 𝛿𝐮 on Γ_out and Γ_𝑠 leads to 𝛽𝜕𝑆𝜕𝐧+𝐮⋅𝑆⋅𝐧=∇𝑇⋅𝐧. (27)

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:𝐕∈𝑉_ad∶=𝐕∈𝐶⁰0,𝛼;𝐶²ℝ²²∶𝑉=0intheneighborhoodofΓ_𝑠. (28)The perturbed domain can be defined by Ω_𝑡=𝐹_𝑡(𝐕)(Ω). Now, we define the cost functional: 𝑗(𝑡)∶=min_{(𝑢𝑡,𝑝𝑡,𝑇𝑡)∈𝑉𝑔(Ω𝑡)×𝑄(Ω𝑡)×𝐻1(Ω𝑡)}max_{(𝑣𝑡,𝑞𝑡,𝑆𝑡)∈𝑉0(Ω𝑡)×𝑄(Ω𝑡)×𝑊(Ω𝑡)}𝐺Ω_𝑡,𝐮_𝑡,𝑝_𝑡,𝑇_𝑡,𝐯_𝑡,𝑞_𝑡,𝑆_𝑡, (29)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: 𝑉_𝑔Ω_𝑡=𝐮∘𝐹_𝑡⁻¹∶𝐮∈𝑉_𝑔(Ω,𝑉₀)Ω_𝑡=𝐯∘𝐹_𝑡⁻¹∶𝐯∈𝑉₀(Ω,𝑄Ω_𝑡)=𝑝∘𝐹_𝑡⁻¹(Ω,𝑊Ω∶𝑝∈𝑄)_𝑡=𝑆∘𝐹_𝑡⁻¹(Ω,𝐻∶𝑆∈𝑊)₁Ω_𝑡=𝑇∘𝐹_𝑡⁻¹∶𝑇∈𝐻¹(Ω), (30)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 𝑗(𝑡)=min_{𝑢,𝑝,𝑇}max_{𝑣,𝑞,𝑆}𝐺Ω_𝑡,𝑢∘𝐹_𝑡⁻¹,𝑝∘𝐹_𝑡⁻¹,𝑇∘𝐹_𝑡⁻¹,𝑣∘𝐹_𝑡⁻¹,𝑞∘𝐹_𝑡⁻¹,𝑆∘𝐹_𝑡⁻¹, (31)where the Lagrangian functional 𝐺Ω_𝑡,𝐮∘𝐹_𝑡⁻¹,𝑝∘𝐹_𝑡⁻¹,𝑇∘𝐹_𝑡⁻¹,𝐯∘𝐹_𝑡⁻¹,𝑞∘𝐹_𝑡⁻¹,𝑆∘𝐹_𝑡⁻¹=𝐼₁(𝑡)+𝐼₂(𝑡)+𝐼₃(𝑡) (32)with 𝐼₁(𝑡)∶=2𝜈_Ω𝑡𝜀𝐮∘𝐹_𝑡⁻¹²1d𝑥+2_Ω𝑡∇𝑇∘𝐹_𝑡⁻¹²𝐼d𝑥,₂(𝑡)∶=−_Ω𝑡2𝜈𝜀𝐯∘𝐹_𝑡⁻¹∶𝜀𝐮∘𝐹_𝑡⁻¹−𝑝∘𝐹_𝑡⁻¹div𝐯∘𝐹_𝑡⁻¹−div𝐮∘𝐹_𝑡⁻¹𝑞∘𝐹_𝑡⁻¹𝐼d𝑥,₃(𝑡)∶=−_Ω𝑡𝛽∇𝑇∘𝐹_𝑡⁻¹⋅∇𝑆∘𝐹_𝑡⁻¹+𝐮∘𝐹_𝑡⁻¹⋅∇𝑇∘𝐹_𝑡⁻¹𝑆∘𝐹_𝑡⁻¹d𝑥. (33) To perform the differentiation of the perturbed Lagrangian functional 𝐺(Ω_𝑡,𝐮∘𝐹_𝑡⁻¹,𝑝∘𝐹_𝑡⁻¹,𝑇∘𝐹_𝑡⁻¹,𝐯∘𝐹_𝑡⁻¹,𝑞∘𝐹_𝑡⁻¹,𝑆∘𝐹_𝑡⁻¹), we introduce the following Hadamard formula: dd𝑡_Ω𝑡𝐹(𝑡,𝑥)d𝑥=_Ω𝑡𝜕𝐹()+𝜕𝑡𝑡,𝑥d𝑥_𝜕Ω𝑡𝐹(𝑡,𝑥)𝐕⋅𝐧_𝑡d𝑠 (34) for a sufficiently smooth functional 𝐹_𝑡∶[0,𝛼]×ℝ²→ℝ.

By Hadamard formula (34), we obtain 𝜕_𝑡𝐺Ω_𝑡,𝐮∘𝐹_𝑡⁻¹,𝑝∘𝐹_𝑡⁻¹,𝑇∘𝐹_𝑡⁻¹,𝐯∘𝐹_𝑡⁻¹,𝑞∘𝐹_𝑡⁻¹,𝑆∘𝐹_𝑡⁻¹_𝑡=0=𝐼^′₁(0)+𝐼^′₂(0)+𝐼^′₃(0), (35) where 𝐼^′₁(0)=4𝜈_Ω𝜀(𝐮)∶𝜀(−D𝐮⋅𝐕)d𝑥+2𝜈_Γ𝑠|𝜀(𝐮)|²𝐕_𝑛+d𝑠_Ω()1∇𝑇⋅∇−D𝑇⋅𝐕d𝑥+2_𝜕Ω|∇𝑇|²𝐕_𝑛𝐼d𝑠,_′2(0)=−_Ω[()(𝐯)(𝐮)()()()()−()−(+2𝜈𝜀−D𝐮⋅𝐕⋅𝜀+2𝜈𝜀⋅𝜀−D𝐯⋅𝐕+D𝐮⋅𝐮⋅−D𝐯⋅𝐕+D−D𝐮⋅𝐕⋅𝐮⋅𝐯+D𝐮⋅−D𝐮⋅𝐕⋅𝐯−∇𝑝⋅𝐕div𝐯−𝑝div−D𝐯⋅𝐕−∇𝑞⋅𝐕div𝐮−𝑞div−D𝐮⋅𝐕)d𝑥_Γ𝑠((𝜀(𝐮)(𝐯𝐕−2𝜈∶𝜀))−D𝐮⋅𝐮⋅𝐯+𝑝div𝐯+div𝐮𝑞_𝑛𝐼d𝑠,_′3(0)=−_Ω()−𝛽∇−D𝑆⋅𝐕⋅∇𝑇d𝑥_Ω()−𝛽∇𝑆⋅∇−D𝑇⋅𝐕d𝑥_Ω()−−D𝐮⋅𝐕⋅∇𝑇𝑆d𝑥_Ω()−𝐮⋅∇−D𝑇⋅𝐕𝑆d𝑥_Ω()−𝐮⋅∇𝑇−D𝑆⋅𝐕d𝑥_Γ𝑠𝛽∇𝑇⋅∇𝑆+𝐮⋅∇𝑇𝑆𝐕_𝑛d𝑠. (36)(37)(38)

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 D𝐮⋅𝐕⋅𝐧=div𝐮𝐕_𝑛,𝜀(𝐮)(𝐯)=(𝜀(𝐮))⋅(𝜀(𝐯)),(𝜀(𝐮))⋅()=(𝜀(𝐮))⋅(𝜀(𝐯))𝐕∶𝜀⋅𝐧⋅𝐧⋅𝐧D𝐯⋅𝐕⋅𝐧⋅𝐧_𝑛 (39) hold on the boundary Γ_𝑠.

Noting 𝐮 and 𝐯 vanish on the boundary Γ_𝑠, 𝐕|_{Γin∪Γ𝑤∪Γout}=0, and employing Lemma 3, (36) can be rewritten as 𝐼^′₁(0)=−2𝜈_ΩΔ𝐮(−D𝐮⋅𝐕)d𝑥+2𝜈_Γ𝑠|𝜀(𝐮)|²𝐕_𝑛−d𝑠_Ω()1Δ𝑇⋅−D𝑇⋅𝐕d𝑥+2_Γ𝑠|∇𝑇|²𝐕_𝑛d𝑠. (40) Since (𝐮,𝑝) and (𝐯,𝑞) satisfy (1) and (14), respectively, (37) reduces to 𝐼^′₂(0)=_Ω𝜈(Δ𝐯+∇div𝐯)(−D𝐮⋅𝐕)d𝑥−2_Γ𝑠(𝜀(𝐯))+⋅𝐧)(−D𝐮⋅𝐕d𝑠_Ω𝜈(Δ𝐮+∇div𝐮)(−D𝐯⋅𝐕)d𝑥−2_Γ𝑠(𝜀(𝐮))−⋅𝐧)(−D𝐯⋅𝐕d𝑠_Ω()+∇𝑝−D𝐯⋅𝐕d𝑥_Γ𝑠𝑝()−−D𝐯⋅𝐕⋅𝐧d𝑠_Ω()+∇𝑞−D𝐮⋅𝐕d𝑥_Γ𝑠𝑞()=−D𝐮⋅𝐕⋅𝐧d𝑠_Ω()+2𝜈Δ𝐮+∇𝑇𝑆)(−D𝐮⋅𝐕d𝑥_Γ𝑠(2𝜈𝜀(𝐮)∶𝜀(𝐯))𝐕_𝑛d𝑠. (41) Similarly, according to Lemma 3, the state system (2), and the adjoint state system (15), (38) becomes 𝐼^′₃(0)=_Ω𝛽()−−D𝑆⋅𝐕⋅Δ𝑇d𝑥_𝜕Ω𝛽()+−D𝑆⋅𝐕∇𝑇⋅𝐧d𝑠_Ω()−𝛽Δ𝑆⋅−D𝑇⋅𝐕d𝑥_𝜕Ω𝛽()+−D𝑇⋅𝐕∇𝑆⋅𝐧d𝑠_Ω()+div𝐮⋅−D𝑇⋅𝐕𝑆d𝑥_Ω()−𝐮⋅∇𝑆⋅−D𝑇⋅𝐕d𝑥_Γ𝑠()−𝐮⋅−D𝑇⋅𝐕𝑆⋅𝐧d𝑥_Ω()−−D𝐮⋅𝐕⋅∇𝑇𝑆d𝑥_Ω()−𝐮⋅∇𝑇−D𝑆⋅𝐕𝑆d𝑥_Γ𝑠𝛽∇𝑇⋅∇𝑆+𝐮⋅∇𝑇𝑆𝐕_𝑛=d𝑠_Ω()+𝛽Δ𝑇−𝐮⋅∇𝑇−D𝑆⋅𝐕d𝑥_Ω()−𝛽Δ𝑆+𝐮⋅∇𝑆−D𝑇⋅𝐕d𝑥_Ω()−−D𝐮⋅𝐕⋅∇𝑇𝑆d𝑥_Γ𝑠𝛽∇𝑆+𝐮𝑆⋅𝐧⋅(−D𝑇⋅𝐕)d𝑠. (42) Adding (40)–(42) together, we finally obtain the boundary expression for the Eulerian derivative of 𝐽(Ω): d𝐽(Ω;𝐕)=2𝜈_Γ𝑠𝜀(𝐮)∶𝜀(𝐯)+|𝜀(𝐮)|²𝐕_𝑛+1d𝑠2_Γ𝑠|∇𝑇|²𝐕_𝑛d𝑠+_Γ𝑠𝛽∇𝑇⋅∇𝑆⋅𝐕_𝑛d𝑠. (43) Since the mapping 𝐕↦d𝐽(Ω;𝐕) is linear and continuous, we get the expression for the shape gradient ∇𝐽=2𝜈𝜀(𝐮)∶𝜀(𝐯)+|𝜀(𝐮)|²+12|∇𝑇|²+𝛽∇𝑇⋅∇𝑆𝐧 (44) 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 S_ℎ⊂𝐿²(Ω) be two families of finite dimensional subspaces parameterized by ℎ which tends to zero. We also define the following functional spaces: 𝑉_0ℎ∶=𝐮_ℎ∈𝑋_ℎ∶𝐮_ℎ=0onΓ_𝑤∪Γ_𝑢∪Γ_in,𝑉_𝑔ℎ∶=𝐮_ℎ∈𝑋_ℎ∶𝐮_ℎ=0onΓ_𝑤∪Γ_in,𝐮_ℎ=𝐠onΓ_in,𝑄_ℎ∶=𝑝_ℎ∈𝑆_ℎ∶_Ω𝑝_ℎ,𝑊d𝑥=0_ℎ∶=𝑆_ℎ∈𝐻¹(Ω)∶𝑆_ℎ=0onΓ_𝑤∪Γ_in. (45) Besides, the following assumptions are supposed to hold. (H1) There exists 𝐶>0 such that for 0≤𝑚≤𝑙, inf_{𝑣ℎ∈𝑉𝑔ℎ}𝐯_ℎ−𝐯₁≤𝐶ℎ^𝑚‖𝐯‖_𝑚+1,∀𝐯∈𝐻^𝑚+1(Ω)^𝑁∩𝑉_𝑔(Ω); (46) (H2) 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_{0≠𝑞ℎ∈𝑄ℎ}sup_{0≠𝑣ℎ∈𝑉ℎ}∫_Ω𝑞_ℎdiv𝐯_ℎd𝑥𝐯_ℎ1𝑞_ℎ0≥𝐶,𝑉_ℎ=𝑉_𝑔ℎor𝑉_0ℎ. (48)The Galerkin finite-element approximation of the state system (1) and (2) in mixed form are as follows: seek𝐮_ℎ,𝑝_ℎ∈𝑉_𝑔ℎ×𝑄_ℎ,suchthat𝐯_ℎ,𝑞_ℎ∈𝑉_0ℎ×𝑄_ℎΩ2𝜈𝜀𝐮_ℎ∶𝜀𝐯_ℎ−𝑝_ℎdiv𝐯_ℎd𝑥=0_Ωdiv𝐮_ℎ𝑞_ℎd𝑥=0,find𝑇_ℎ∈𝐻¹_ℎ,suchthat_Ω𝛽∇𝑇_ℎ⋅∇𝑆_ℎ+𝐮_ℎ⋅∇𝑇_ℎ𝑆_ℎd𝑥=0,∀𝑆_ℎ∈𝑊_ℎ. (49) Similarly, we can obtain the finite element approximation for adjoint state system, seek𝐯_ℎ,𝑞_ℎ,𝑆_ℎ∈𝑉_0ℎ×𝑄_ℎ×𝑊_ℎ,suchthat𝝋_ℎ,𝜉_ℎ,𝛾_ℎ∈𝑉_0ℎ×𝑄_ℎ×𝑊_ℎΩ2𝜈𝜀𝐯_ℎ∶𝜀𝝋_ℎ−𝑞_ℎdiv𝝋_ℎ+𝑆_ℎ∇𝑇_ℎ=4𝜈_Ω𝜀𝐮_ℎ∶𝜀𝝋_ℎd𝑥_Ωdiv𝐯_ℎ𝜉_ℎd𝑥=0_Ω𝛽∇𝑆_ℎ⋅∇𝛾_ℎ+𝐮_ℎ⋅∇𝑆_ℎ𝛾_ℎd𝑥=_Ω∇𝑇_ℎ⋅∇𝛾_ℎd𝑥. (50) We also have the discrete cost functional 𝐽_ℎ(Ω)=2𝜈_Ω𝜀𝐮_ℎ²1d𝑥+2_Ω∇𝑇_ℎ²d𝑥, (51) and the discrete shape gradient ∇𝐽_ℎ=2𝜈𝜀𝐮_ℎ∶𝜀𝐯_ℎ+𝜀𝐮_ℎ²+12∇𝑇_ℎ²+𝛽∇𝑇_ℎ⋅∇𝑆_ℎ𝐧. (52)

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 min_Ω∈ℝ2𝐺(Ω)=𝐽(Ω)+𝑙𝑉(Ω), (53) where 𝑉(Ω)∶=∫_Ωd𝑥 and 𝑙 is a positive Lagrange multiplier. The Eulerian derivative of 𝐺(Ω) is d𝐺(Ω;𝐕)=_Γ𝑠∇𝐺⋅𝐕d𝑠, (54) and the shape gradient ∇𝐺∶=[2𝜈(𝜀(𝐮_ℎ)∶𝜀(𝐯_ℎ)+|𝜀(𝐮_ℎ)|²)+(1/2)|∇𝑇_ℎ|²+𝛽∇𝑇_ℎ⋅∇𝑆_ℎ+𝑙]𝐧. Ignoring regularization, a descent direction is found by defining 𝐕=−ℎ_𝑘∇𝐺, and then we can update the shape Ω as Ω_𝑘=(I+ℎ_𝑘𝐕)Ω, 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 −Δ𝐝+𝐝=0inΩ,𝐝=0,onΓ_in∪Γ_out∪Γ_𝑤,D𝐝⋅𝐧=−∇𝐺onΓ_𝑠. (55)

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 Ω_𝑘+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.

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: d𝐺(Ω;𝐕)=d𝐽(Ω;𝐕)+𝑙d𝑉(Ω;𝐕)=0, (56) at least in the average sense on the boundary Γ_𝑠; that is, ∫𝑙=−^Γ𝑠∇𝐽d𝑠∫_Γ𝑠d𝑠. (57) Finally in Step (3) of our algorithm, we refresh the Lagrange multiplier by 𝑙_𝑘+1=𝑙_𝑘+𝑙2+𝑚𝑉Ω_𝑘−𝑉_target(Ω)𝑉_target(Ω), (58) 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 ℎ₀ so that ℎ₀∇𝐺Ω_𝑘⋅𝐝_𝑘=ℎ_𝑘−1∇𝐽Ω_𝑘−1⋅𝐝_𝑘−1. (59)

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: ℋ_ℎ=𝐼_𝐿²∇𝐼_𝐿²∇𝑢_ℎ, (60) 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 𝜕𝐷∶=Γ_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.

Figure 1 

Schematic illustration for the boundary of domain Ω.

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 𝒪∶=Ω⊂ℝ²∶𝜕𝐷isfixed,thearea𝑉_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 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.

(a) 𝐮1 for initial shape

(b) 𝐮2 for initial shape

(c) 𝐮1 for optimal shape

(d) 𝐮2 for optimal shape

Figure 3 

Comparison of the initial shape and optimal shape (Re=100).

(a) 𝑝 for initial shape

(b) 𝑇 for initial shape

(c) 𝑝 for optimal shape

(d) 𝑇 for optimal shape

Figure 4 

Comparison of the initial shape and optimal shape (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.

Figure 5 

Convergence history of the cost function (Re=100).

Figure 6 

Convergence history of the cost function (Re=400).

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.

Acknowledgment

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.