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

Authors: Wenjing Yan; Axia Wang; Yichen Ma
Journal: Mathematical Problems in Engineering
 (2013)

---

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


---

## Body

## 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 Section2, 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 bodyS with the boundary Γs is located in an external flow. Let Γ denote the boundary of Ω, and suppose that an incompressible viscous flow occupies Ω. The boundary Γ:=Γw∪Γs∪Γin∪Γout, Γin is the inflow boundary, Γout denotes the outflow boundary, Γw represents the boundary corresponding to the fluid wall, and Γs is the boundary which is to be optimized. The unknowns of the full model are the fluid velocity u, the pressure p, and the temperature T.Find(u,p,T) satisfying
(1)-νΔu+∇p=0inΩdivu=0inΩu=0onΓw∪Γsu=gonΓinσ(u,p)·n=0onΓout,(2)-βΔT+u·∇T=0inΩ∂T∂n=0onΓout∪ΓsT=T1onΓwT=T2onΓin,
where ν stands for the inverse of Reynolds number, the stress tensor is defined by σ(u,p):=-pI+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. LetL2(Ω) be the space of square integrable real-valued functions on Ω with the usual norm. The space Hk(Ω), where k=1,2,…,n denotes the standard Sobolev space on Ω (see [9]), that is, the space of functions with generalized derivatives of order up to k in L2(Ω). Consider(3)V0(Ω):={u∈(H2(Ω))2:u=0onΓw∪Γs∪Γin},Vg(Ω):={u∈(H2(Ω))2:u=0onΓw∪Γs,u=gonΓin{(H2(Ω))2}},W(Ω):={S∈H1(Ω):S=0onΓw∪Γin},Q(Ω):={p∈H1(Ω):∫Ωpdx=0}.Our aim is to optimize the shape of the boundaryΓs that minimizes a given cost functional J 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:
(4)minΩ∈𝒪J(Ω)=2ν∫Ω‍|ε(u)|2dx+12∫Ω‍|∇T|2dx,
where the boundary Γw∪Γin∪Γout is fixed. An example of the admissible set 𝒪 is
(5)𝒪:={Ω⊂ℝ2:Γw∪Γin∪Γoutis  fixed,∫Ωdx=constant}.
Let Ω be of piecewise C1, 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 ℝ2 with the boundary ∂Ω which belongs to piecewise Ck, and a velocity space V∈Ek:={V∈C([0,α];𝒟k(Ω-,ℝ2))}, where α is a small positive real number and 𝒟k(Ω-,ℝ2) denotes the space of all k-times continuous differentiable functions with compact support contained in ℝ2. The velocity field
(6)V(t)(x)=V(t,x),x∈Ω,t≥0
belongs to 𝒟k(Ω-,ℝ2) for each t. It can generate transformations
(7)Ft(V)X=x(t,X),t≥0,X∈Ω
through the following dynamical system
(8)dxdt(t,X)=V(t,x(t)),x(0,X)=X
with the initial value X. We denote the transformed domain Ft(V)(Ω) by Ωt(V) at t≥0, and also set Γt:=Ft(Γ). There exists an interval I=[0,β), 0<β≤α, and a one-to-one map Ft from Ω- onto Ω- such that (see [6]) (1)
F0=I;(2)
(t,x)↦Ft(x) belongs to C1(I;Ck(Ω)) with Ft(∂Ω)=∂Ω;(3)
(t,x)↦Ft-1(x) belongs to C(I;Ck(Ω)). In addition, for sufficiently small t>0, the Jacobian Jt of Ft is strictly positive:
(9)Jt(x):=det|DFt(x)|=detDFt(x)>0,
where DFt(x) denotes the Jacobian matrix of the transformation Ft evaluated at a point x∈Ω associated with the velocity field V. We also introduce the following notations: DFt-1(x) represents the inverse of the matrix DFt(x), and DFt-1(x)* is the transpose of the matrix DFt-1(x).Definition 1.
LetJ(Ω) 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:
(10)limt↘0J(Ωt)-J(Ω)t:=dJ(Ω;V).
Furthermore, if the mapV↦dJ(Ω;V):Ek→ℝ is linear and continuous, we say that J is shape differentiable at Ω. In the distributional sense, we have
(11)dJ(Ω;V)=〈∇J,V〉𝒟k(Ω-,ℝ2)′×𝒟k(Ω-,ℝ2),
where ∇J is the shape gradient of J 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 functionalJ(Ω).LetΩ be of class C2, and we derive the variational formulation of the state system (1) and (2) in appropriate Sobolev spaces:
(12)seek(u,p)∈Vg(Ω)×Q(Ω),such  that∫Ω‍[2νε(u):ε(v)-pdivv]dx=0,∀v∈V0(Ω),∫Ω‍divuqdx=0,∀q∈Q(Ω),findT∈H1(Ω),such  that∫Ω(β∇T·∇S+u·∇TS)dx=0,∀S∈W(Ω).
Now, we will prove the main theoretical results of the paper.Theorem 2.
LetΩ be of class C2 and V∈E2, the cost functional J(Ω) possesses the shape gradient ∇J which can be expressed as
(13)∇J=[2ν(ε(u):ε(v)+|ε(u)|2)+12|∇T|2+β∇T·∇S]n,
where the adjoint states v and S satisfy the following adjoint system:
(14)-νΔv+∇q+S∇T=-2νΔuinΩdivv=0inΩσ(v,q)·n-4νε(u)·n=0,onΓoutv=0onΓin∪Γw∪Γs,(15)βΔS+u·∇S=ΔTinΩβ∂S∂n+u·S·n=∇T·nonΓout∪ΓsS=0onΓin∪Γw.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)G(Ω,u,p,T,v,q,S)=J(Ω)-L(Ω,u,p,T,v,q,S),
where
(17)L(Ω,u,p,T,v,q,S)=∫Ω‍[2νε(u):ε(v)-pdivv-divuq]dx+∫Ω‍(β∇T·∇S+u·∇TS)dx.
Problem (16) can be put in the following form: (18)minΩ∈𝒪min(u,p,T)∈Vg(Ω)×Q(Ω)×H1(Ω)max(v,q,S)∈V0(Ω)×Q(Ω)×W(Ω)G(Ω,u,p,T,v,q,S).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 J(Ω) by using the adjoint system.The first optimality condition for the problem can be established as follows:(19)min(u,p,T)∈Vg(Ω)×Q(Ω)×H1(Ω)max(v,q,S)∈V0(Ω)×Q(Ω)×W(Ω)G(Ω,u,p,T,v,q,S).Formally the adjoint equations are defined from the Euler-Lagrange equations of the Lagrange functional G. Clearly, the variation of G with respect to (v,q,S) can recover the state system and its mixed weak formulation (12). In order to find the adjoint state system, we differentiate G with respect to p in the direction δp:
(20)∂G∂p(Ω,u,p,T,v,q,S)·δp=∫Ω‍δpdivvdx=0.
Taking δp with compact support in Ω gives
(21)divv=0.Then, we differentiate G with respect to u in the direction δu and apply Green’s formula:(22)∂G∂u(Ω,u,p,T,v,q,S)·δu=∫Ω‍(-2νΔu+νΔv-∇q)·δudx-∫Ω‍S∇T·δudx-2ν∫∂Ω‍ε(v)·n·δuds+4ν∫∂Ωε(u)·n·δuds+∫∂Ωqδu·nds=∫Ω‍(-2νΔu+νΔv-∇q)·δudx-∫Ω‍S∇T·δudx-∫∂Ω‍σ(v,q)·n·δuds+4ν∫∂Ω‍ε(u)·n·δuds.Considering δu with compact support in Ω, we obtain (23)-νΔv+∇q+S∇T=-2νΔu.Then, varyingδu on Γout gives
(24)σ(v,q)·n-4νε(u)·n=0onΓout.Similarly, we differentiateG with respect to T in the direction δT:
(25)∂G∂T(Ω,u,p,T,v,q,S)·δT=∫Ω(βΔS+u·∇S)·δTdx-∫Γout∪Γs‍(β∂S∂n+u·S·n)δTds+∫Ω‍∇T·∇(δT)dx=∫Ω‍(βΔS+u·∇S)·δTdx-∫Γout∪Γs‍(β∂S∂n+u·S·n)δTds-∫Ω‍ΔT·δTdx+∫∂Ω‍∇T·δT·nds.According to the compact support ofδu in Ω, we can get
(26)βΔS+u·∇S=ΔT.In addition, varyingδu on Γout and Γs leads to (27)β∂S∂n+u·S·n=∇T·n.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Γs and consider the transformation Ft(V); then the flow of the velocity field can be expressed as follows:(28)V∈Vad:={V∈C0(0,α;(C2(ℝ2))2):&V∈Vad:=V=0in  the  neighborhood  ofΓs{(0,α;(C2(ℝ2))2)}}.The perturbed domain can be defined by Ωt=Ft(V)(Ω). Now, we define the cost functional: (29)j(t):=min(ut,pt,Tt)∈Vg(Ωt)×Q(Ωt)×H1(Ωt)max(vt,qt,St)∈V0(Ωt)×Q(Ωt)×W(Ωt)G(Ωt,ut,pt,Tt,vt,qt,St),where (ut,pt,Tt) and (vt,qt,St) satisfy (14) and (15) on the perturbed domain Ωt, respectively.Our object is to derive the derivative with respect to the variablet. Unfortunately, the Sobolev spaces Vg(Ωt), V0(Ωt), Q(Ωt), and W(Ωt) depend on the parameter t, so we need to introduce the so-called function space parametrization technique which consists in transporting the different quantities defined on the variable domain Ωt back into the reference domain Ω which does not depend on the perturbation parameter t. Since the functionals above mentioned are defined in a fixed domain Ω with respect to the parameter t, we can apply the differential calculus. Hence, we define the following functional spaces:
(30)Vg(Ωt)={u∘Ft-1:u∈Vg(Ω)},V0(Ωt)={v∘Ft-1:v∈V0(Ω)},Q(Ωt)={p∘Ft-1:p∈Q(Ω)},W(Ωt)={S∘Ft-1:S∈W(Ω)},H1(Ωt)={T∘Ft-1:T∈H1(Ω)},where “∘” denotes the composition of the two maps.SinceFt and Ft-1 are diffeomorphisms, these parametrizations cannot change the value of the saddle point. We can rewrite (29) as (31)j(t)=min(u,p,T)max(v,q,S)G(Ωt,u∘Ft-1,p∘Ft-1,T∘Ft-1,j(t)=min(u,p,T)max(v,q,S)Gv∘Ft-1,q∘Ft-1,S∘Ft-1),where the Lagrangian functional (32)G(Ωt,u∘Ft-1,p∘Ft-1,T∘Ft-1,v∘Ft-1,q∘Ft-1,S∘Ft-1)=I1(t)+I2(t)+I3(t)with (33)I1(t):=2ν∫Ωt‍|ε(u∘Ft-1)|2dx+12∫Ωt‍|∇(T∘Ft-1)|2dx,I2(t):=-∫Ωt‍[2νε(v∘Ft-1):ε(u∘Ft-1)I2(t)≔-∫Ωt‍-(p∘Ft-1)div(v∘Ft-1)I2(t)≔-∫Ωt‍-div(u∘Ft-1)(q∘Ft-1)]dx,I3(t)≔-∫Ωt‍[β∇(T∘Ft-1)·∇(S∘Ft-1)I3(t)≔-∫Ωt‍+(u∘Ft-1)·∇(T∘Ft-1)(S∘Ft-1)]dx.
To perform the differentiation of the perturbed Lagrangian functional G(Ωt,u∘Ft-1,p∘Ft-1,T∘Ft-1,v∘Ft-1,q∘Ft-1,S∘Ft-1), we introduce the following Hadamard formula:
(34)ddt∫Ωt‍F(t,x)dx=∫Ωt‍∂F∂t(t,x)dx+∫∂Ωt‍F(t,x)V·ntds
for a sufficiently smooth functional Ft:[0,α]×ℝ2→ℝ.By Hadamard formula (34), we obtain
(35)∂tG(Ωt,u∘Ft-1,p∘Ft-1,T∘Ft-1,v∘Ft-1,q∘Ft-1,S∘Ft-1)|t=0=I1′(0)+I2′(0)+I3′(0),
where
(36)I1′(0)=4ν∫Ωɛ(u):ɛ(-Du·V)dx+2ν∫Γs|ɛ(u)|2Vnds+∫Ω∇T·∇(-DT·V)dx+12∫∂Ω‍|∇T|2Vnds,(37)I2′(0)=-∫Ω[2νɛ(-Du·V)·ɛ(v)-∫‍Ω+2νɛ(u)·ɛ(-Dv·V)+Du·u·(-Dv·V)-∫‍Ω+D(-Du·V)·u·v+Du·(-Du·V)·v-∫‍Ω-(-∇p·V)divv-pdiv(-Dv·V)-∫‍Ω-(-∇q·V)divu-qdiv(-Du·V)]dx+∫Γs‍(-2ν(ɛ(u):ɛ(v))-Du·u·v+∫Γs‍+pdivv+divuq)Vnds,(38)I3′(0)=-∫Ω‍β∇(-DS·V)·∇Tdx-∫Ω‍β∇S·∇(-DT·V)dx-∫Ω‍(-Du·V)·∇TSdx-∫Ω‍u·∇(-DT·V)Sdx-∫Ω‍u·∇T(-DS·V)dx-∫Γs‍(β∇T·∇S+u·∇TS)Vnds.To simplify (36)–(38), we introduce the following lemma.Lemma 3 (see [6]).
If two vector functionsu and v vanish on the boundary Γs, the following identities
(39)Du·V·n=divuVn,ɛ(u):ɛ(v)=(ɛ(u)·n)·(ɛ(v)·n),(ɛ(u)·n)·(Dv·V)=(ɛ(u)·n)·(ɛ(v)·n)Vn
hold on the boundary Γs.Notingu and v vanish on the boundary Γs, V|Γin∪Γw∪Γout=0, and employing Lemma 3, (36) can be rewritten as
(40)I1′(0)=-2ν∫ΩΔu(-Du·V)dx+2ν∫Γs‍|ɛ(u)|2Vnds-∫ΩΔT·(-DT·V)dx+12∫Γs‍|∇T|2Vnds.
Since (u,p) and (v,q) satisfy (1) and (14), respectively, (37) reduces to
(41)I2′(0)=∫Ων(Δv+∇divv)(-Du·V)dx-2∫Γs‍(ɛ(v)·n)(-Du·V)ds+∫Ων(Δu+∇divu)(-Dv·V)dx-2∫Γs‍(ɛ(u)·n)(-Dv·V)ds-∫Ω∇p(-Dv·V)dx+∫Γs‍p(-Dv·V)·nds-∫Ω∇q(-Du·V)dx+∫Γs‍q(-Du·V)·nds=∫Ω(2νΔu+∇TS)(-Du·V)dx+∫Γs‍(2νɛ(u):ɛ(v))Vnds.
Similarly, according to Lemma 3, the state system (2), and the adjoint state system (15), (38) becomes
(42)I3′(0)=∫Ω‍β(-DS·V)·ΔTdx-∫∂Ω‍β(-DS·V)∇T·nds+∫Ω‍βΔS·(-DT·V)dx-∫∂Ω‍β(-DT·V)∇S·nds+∫Ω‍divu·(-DT·V)Sdx+∫Ω‍u·∇S·(-DT·V)dx-∫Γs‍u·(-DT·V)S·ndx-∫Ω‍(-Du·V)·∇TSdx-∫Ω‍u·∇T(-DS·V)Sdx-∫Γs‍(β∇T·∇S+u·∇TS)Vnds=∫Ω‍(βΔT-u·∇T)(-DS·V)dx+∫Ω‍(βΔS+u·∇S)(-DT·V)dx-∫Ω‍(-Du·V)·∇TSdx-∫Γs‍(β∇S+uS)·n·(-DT·V)ds.
Adding (40)–(42) together, we finally obtain the boundary expression for the Eulerian derivative of J(Ω):
(43)dJ(Ω;V)=2ν∫Γs‍[ɛ(u):ɛ(v)+|ɛ(u)|2]Vnds+12∫Γs‍|∇T|2Vnds+∫Γs‍β∇T·∇S·Vnds.
Since the mapping V↦dJ(Ω;V) is linear and continuous, we get the expression for the shape gradient
(44)∇J=[2ν(ɛ(u):ɛ(v)+|ɛ(u)|2)+12|∇T|2+β∇T·∇S]n
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 ℝ2 and only consider the conforming finite-element approximations. Let Xh⊂(H2(Ω))N and Sh⊂L2(Ω) be two families of finite dimensional subspaces parameterized by h which tends to zero. We also define the following functional spaces:
(45)V0h:={uh∈Xh:uh=0onΓw∪Γu∪Γin},Vgh:={uh∈Xh:uh=0onΓw∪Γin,uh=gonΓin},Qh:={ph∈Sh:∫Ωphdx=0},Wh:={Sh∈H1(Ω):Sh=0onΓw∪Γin}.
Besides, the following assumptions are supposed to hold.  (H1)
There existsC>0 such that for 0≤m≤l,
(46)infvh∈Vgh∥vh-v∥1≤Chm∥v∥m+1,∀v∈(Hm+1(Ω))N∩Vg(Ω);(H2)
there existsC>0 such that for 0≤m≤l′,
(47)infqh∈Qh∥qh-q∥0≤Chm∥q∥m,∀q∈Hm(Ω)∩Q(Ω);(H3)
the Ladyzhenskaya-Brezzi-Babuska inf-sup condition is verified; that is, there existsC>0 such that
(48)inf0≠qh∈Qhsup0≠vh∈Vh∫Ωqhdivvhdx∥vh∥1∥qh∥0≥C,Vh=VghorV0h.The Galerkin finite-element approximation of the state system (1) and (2) in mixed form are as follows:
(49)seek(uh,ph)∈Vgh×Qh,suchthat(vh,qh)∈V0h×Qh∫Ω‍[2νɛ(uh):ɛ(vh)-phdivvh]dx=0∫Ωdivuhqhdx=0,findTh∈Hh1,suchthat∫Ω(β∇Th·∇Sh+uh·∇ThSh)dx=0,∀Sh∈Wh.
Similarly, we can obtain the finite element approximation for adjoint state system,
(50)seek(vh,qh,Sh)∈V0h×Qh×Wh,suchthat(φh,ξh,γh)∈V0h×Qh×Wh∫Ω‍[2νɛ(vh):ɛ(φh)-qhdivφh+Sh∇Th]=4ν∫Ωɛ(uh):ɛ(φh)dx∫Ω‍divvhξhdx=0∫Ω‍β∇Sh·∇γh+uh·∇Shγhdx=∫Ω‍∇Th·∇γhdx.
We also have the discrete cost functional
(51)Jh(Ω)=2ν∫Ω‍|ɛ(uh)|2dx+12∫Ω‍|∇Th|2dx,
and the discrete shape gradient
(52)∇Jh=[{12|∇Th|2}2ν(ɛ(uh):ɛ(vh)+|ɛ(uh)|2)+12|∇Th|2+β∇Th·∇Sh]n.


### 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 functionalJ(Ω). For the minimization problem (4), we rather work with the unconstrained minimization problem
(53)minΩ∈ℝ2G(Ω)=J(Ω)+lV(Ω),
where V(Ω):=∫Ωdx and l is a positive Lagrange multiplier. The Eulerian derivative of G(Ω) is
(54)dG(Ω;V)=∫Γs‍∇G·Vds,
and the shape gradient ∇G:=[2ν(ɛ(uh):ɛ(vh)+|ɛ(uh)|2)+(1/2)|∇Th|2+β∇Th·∇Sh+l]n. Ignoring regularization, a descent direction is found by defining V=-hk∇G, and then we can update the shape Ω as Ωk=(I+hkV)Ω, where hk is a descent step at kth 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, (H1(Ω))2. In this case, the descent direction is the unique element d∈(H1(Ω))2 of the problem
(55)-Δd+d=0inΩ,d=0,onΓin∪Γout∪Γw,Dd·n=-∇GonΓs.The resulting algorithm consists of the following parts:(1)
give an initial shapeΩ0, an initial step h0, and a Lagrange multiplier l0;(2)
solve the state system and adjoint state system, then we can evaluate the descent directiondk by (55) with Ω=Ωk and l=lk;(3)
setΩk+1=(I-hkdk)Ωk, where hk 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 multiplierl in the optimization problem (53). The value of l 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 l, 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:
(56)dG(Ω;V)=dJ(Ω;V)+ldV(Ω;V)=0,
at least in the average sense on the boundary Γs; that is,
(57)l=-∫Γs‍∇Jds∫Γsds.
Finally in Step (3) of our algorithm, we refresh the Lagrange multiplier by
(58)lk+1=lk+l2+m|V(Ωk)-Vtarget(Ω)|Vtarget(Ω),
where m is a small positive parameter and Vtarget(Ω) denotes the target volume of the shape.


#### 5.2.2. Step Size

The choice of the descent step sizehk 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 h0 so that
(59)h0∇G(Ωk)·dk=hk-1∇J(Ωk-1)·dk-1.


#### 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 theP1 bubble—P1 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 yh can be approximated by using a recovery method, such as the Zienkiewicz-Zhu recovery procedure [14], the simple linear fitting [15], or the double L2 projection:
(60)ℋh=IL2(∇(IL2(∇uh))),
where IL2 denotes the L2 projection on the P1 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 Figure1, corresponding to an external flow around a solid body S. We reduce the problem to a bounded domain D by introducing an artificial boundary ∂D:=Γin∪Γw∪Γout which has to be taken sufficiently far from S so that the corresponding flow is a good approximation of the unbounded external flow around S and Ω:=D∖S- is the effective domain. In addition, the boundary Γs:=∂S is to be optimized.Figure 1
Schematic illustration for the boundary of domainΩ.We chooseD to be a rectangle (0,9)×(0,3) and S 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
(61)𝒪:={Ω⊂ℝ2:∂Disfixed,theareaVtarget(Ω)=2.1}.
We set the initial shape of the body S to be a circle of center (4.5,1.5) with radius r=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 atRe=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Ω.Figures3 and 4 demonstrate the comparison between the initial shape and optimal shape for the distribution of the velocity u=(u1,u2)T, the pressure p, and the temperature T with the Reynolds numbers Re = 100.Comparison of the initial shape and optimal shape (Re=100).
(a)
u
1 for initial shape
(b)
u
2 for initial shape
(c)
u
1 for optimal shape
(d)
u
2 for optimal shapeComparison of the initial shape and optimal shape (Re=100).
(a)
p for initial shape
(b)
T for initial shape
(c)
p for optimal shape
(d)
T for optimal shapeWe 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 Figures5 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).


## 5.1. Discretization of the Optimization Problem

We suppose thatΩ is a bounded polygonal domain of ℝ2 and only consider the conforming finite-element approximations. Let Xh⊂(H2(Ω))N and Sh⊂L2(Ω) be two families of finite dimensional subspaces parameterized by h which tends to zero. We also define the following functional spaces:
(45)V0h:={uh∈Xh:uh=0onΓw∪Γu∪Γin},Vgh:={uh∈Xh:uh=0onΓw∪Γin,uh=gonΓin},Qh:={ph∈Sh:∫Ωphdx=0},Wh:={Sh∈H1(Ω):Sh=0onΓw∪Γin}.
Besides, the following assumptions are supposed to hold.  (H1)
There existsC>0 such that for 0≤m≤l,
(46)infvh∈Vgh∥vh-v∥1≤Chm∥v∥m+1,∀v∈(Hm+1(Ω))N∩Vg(Ω);(H2)
there existsC>0 such that for 0≤m≤l′,
(47)infqh∈Qh∥qh-q∥0≤Chm∥q∥m,∀q∈Hm(Ω)∩Q(Ω);(H3)
the Ladyzhenskaya-Brezzi-Babuska inf-sup condition is verified; that is, there existsC>0 such that
(48)inf0≠qh∈Qhsup0≠vh∈Vh∫Ωqhdivvhdx∥vh∥1∥qh∥0≥C,Vh=VghorV0h.The Galerkin finite-element approximation of the state system (1) and (2) in mixed form are as follows:
(49)seek(uh,ph)∈Vgh×Qh,suchthat(vh,qh)∈V0h×Qh∫Ω‍[2νɛ(uh):ɛ(vh)-phdivvh]dx=0∫Ωdivuhqhdx=0,findTh∈Hh1,suchthat∫Ω(β∇Th·∇Sh+uh·∇ThSh)dx=0,∀Sh∈Wh.
Similarly, we can obtain the finite element approximation for adjoint state system,
(50)seek(vh,qh,Sh)∈V0h×Qh×Wh,suchthat(φh,ξh,γh)∈V0h×Qh×Wh∫Ω‍[2νɛ(vh):ɛ(φh)-qhdivφh+Sh∇Th]=4ν∫Ωɛ(uh):ɛ(φh)dx∫Ω‍divvhξhdx=0∫Ω‍β∇Sh·∇γh+uh·∇Shγhdx=∫Ω‍∇Th·∇γhdx.
We also have the discrete cost functional
(51)Jh(Ω)=2ν∫Ω‍|ɛ(uh)|2dx+12∫Ω‍|∇Th|2dx,
and the discrete shape gradient
(52)∇Jh=[{12|∇Th|2}2ν(ɛ(uh):ɛ(vh)+|ɛ(uh)|2)+12|∇Th|2+β∇Th·∇Sh]n.


## 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 functionalJ(Ω). For the minimization problem (4), we rather work with the unconstrained minimization problem
(53)minΩ∈ℝ2G(Ω)=J(Ω)+lV(Ω),
where V(Ω):=∫Ωdx and l is a positive Lagrange multiplier. The Eulerian derivative of G(Ω) is
(54)dG(Ω;V)=∫Γs‍∇G·Vds,
and the shape gradient ∇G:=[2ν(ɛ(uh):ɛ(vh)+|ɛ(uh)|2)+(1/2)|∇Th|2+β∇Th·∇Sh+l]n. Ignoring regularization, a descent direction is found by defining V=-hk∇G, and then we can update the shape Ω as Ωk=(I+hkV)Ω, where hk is a descent step at kth 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, (H1(Ω))2. In this case, the descent direction is the unique element d∈(H1(Ω))2 of the problem
(55)-Δd+d=0inΩ,d=0,onΓin∪Γout∪Γw,Dd·n=-∇GonΓs.The resulting algorithm consists of the following parts:(1)
give an initial shapeΩ0, an initial step h0, and a Lagrange multiplier l0;(2)
solve the state system and adjoint state system, then we can evaluate the descent directiondk by (55) with Ω=Ωk and l=lk;(3)
setΩk+1=(I-hkdk)Ωk, where hk 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 multiplierl in the optimization problem (53). The value of l 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 l, 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:
(56)dG(Ω;V)=dJ(Ω;V)+ldV(Ω;V)=0,
at least in the average sense on the boundary Γs; that is,
(57)l=-∫Γs‍∇Jds∫Γsds.
Finally in Step (3) of our algorithm, we refresh the Lagrange multiplier by
(58)lk+1=lk+l2+m|V(Ωk)-Vtarget(Ω)|Vtarget(Ω),
where m is a small positive parameter and Vtarget(Ω) denotes the target volume of the shape.


### 5.2.2. Step Size

The choice of the descent step sizehk 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 h0 so that
(59)h0∇G(Ωk)·dk=hk-1∇J(Ωk-1)·dk-1.


### 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.2.1. Lagrange Multiplier

We now discuss the choice of the Lagrange multiplierl in the optimization problem (53). The value of l 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 l, 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:
(56)dG(Ω;V)=dJ(Ω;V)+ldV(Ω;V)=0,
at least in the average sense on the boundary Γs; that is,
(57)l=-∫Γs‍∇Jds∫Γsds.
Finally in Step (3) of our algorithm, we refresh the Lagrange multiplier by
(58)lk+1=lk+l2+m|V(Ωk)-Vtarget(Ω)|Vtarget(Ω),
where m is a small positive parameter and Vtarget(Ω) denotes the target volume of the shape.


## 5.2.2. Step Size

The choice of the descent step sizehk 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 h0 so that
(59)h0∇G(Ωk)·dk=hk-1∇J(Ωk-1)·dk-1.


## 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 theP1 bubble—P1 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 yh can be approximated by using a recovery method, such as the Zienkiewicz-Zhu recovery procedure [14], the simple linear fitting [15], or the double L2 projection:
(60)ℋh=IL2(∇(IL2(∇uh))),
where IL2 denotes the L2 projection on the P1 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 Figure1, corresponding to an external flow around a solid body S. We reduce the problem to a bounded domain D by introducing an artificial boundary ∂D:=Γin∪Γw∪Γout which has to be taken sufficiently far from S so that the corresponding flow is a good approximation of the unbounded external flow around S and Ω:=D∖S- is the effective domain. In addition, the boundary Γs:=∂S is to be optimized.Figure 1
Schematic illustration for the boundary of domainΩ.We chooseD to be a rectangle (0,9)×(0,3) and S 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
(61)𝒪:={Ω⊂ℝ2:∂Disfixed,theareaVtarget(Ω)=2.1}.
We set the initial shape of the body S to be a circle of center (4.5,1.5) with radius r=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 atRe=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Ω.Figures3 and 4 demonstrate the comparison between the initial shape and optimal shape for the distribution of the velocity u=(u1,u2)T, the pressure p, and the temperature T with the Reynolds numbers Re = 100.Comparison of the initial shape and optimal shape (Re=100).
(a)
u
1 for initial shape
(b)
u
2 for initial shape
(c)
u
1 for optimal shape
(d)
u
2 for optimal shapeComparison of the initial shape and optimal shape (Re=100).
(a)
p for initial shape
(b)
T for initial shape
(c)
p for optimal shape
(d)
T for optimal shapeWe 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 Figures5 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.



---

*Source: 728534-2013-04-22.xml*