# Fixed Domain Algorithms in Shape Optimization for Stationary Navier-Stokes Equations

Andrei Halanay, Cornel Marius Murea

► To cite this version:

Andrei Halanay, Cornel Marius Murea. Fixed Domain Algorithms in Shape Optimization for Stationary Navier-Stokes Equations. 25th System Modeling and Optimization (CSMO), Sep 2011, Berlin, Germany. pp.378-386, 10.1007/978-3-642-36062-6_38. hal-01347559

HAL Id: hal-01347559

https://inria.hal.science/hal-01347559

Submitted on 21 Jul 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Fixed domain algorithms in shape optimization for stationary Navier-Stokes equations

Andrei Halanay¹ and Cornel Marius Murea²

¹ Department of Mathematics and Informatics, University Politehnica of Bucharest, 313 Splaiul Independenței RO-060042 Bucharest, Romania, **halanay@mathem.pub.ro**

² Laboratoire de Mathématiques, Informatique et Applications, Université de Haute Alsace 4-6, rue des Frères Lumière, 68093 Mulhouse, France **cornel.murea@uha.fr** http://www.edp.lmia.uha.fr/murea/

**Abstract.** The paper aims to illustrate the algorithm developed in the paper [6] in some specific problems of shape optimization issued from fluid mechanics. Using the fictitious domain method with penalization, the fluid equations will be solved in a fixed domain. The admissible shapes are parametrized by continuous function defined in the fixed domain, then the shape optimization problem becomes an optimal control problem, where the control is the parametrization of the shape. We get the directional derivative of the cost function by solving co-state equation. Numerical results are obtained using a gradient type algorithm.

**Keywords:** shape optimization; optimal control; penalization; approximate extension; gradient method;

## 1 Introduction

The paper presents some applications of an algorithm developed in [6]. This algorithm is based on a method that uses a penalization of the stationary Navier-Stokes equation that approximates its solution by functions defined on a larger fixed domain. The unknown domains are parametrized by functions in a certain subspace of the space of continuous functions on the larger fixed domain.

The approximating extension technique makes possible the approximation of the solution to the shape optimization problem by a solution of an optimal control problem. The basic reference will be [6]. For shape optimization, the general references are [13], [3] and for optimal control [7], [10]. In particular, for shape optimization for fluids a standard work is [9]. In optimal design related to optimal control, relevant contributions to the topic of this paper are [1], [4], [15].
In Section 2, the shape optimization problem for steady Navier-Stokes is presented. The directional derivative of the cost function is given in Section 3. A gradient type algorithm is also introduced. In Section 4, numerical results are presented in order to design a nozzle.

## 2 Formulation of the shape optimization problem and approximating extensions

Let $d$ be a natural number, $d \le 4$, let $D \subset \mathbb{R}^d$ be a bounded fixed domain and suppose a family $\mathcal{O}$ of admissible subdomains $\Omega \subset D$ is given, satisfying an uniform Lipschitz condition on the boundary $\partial\Omega$.

With the standard notations from [16], $\mathcal{V}(\Omega) = \{y \in \mathcal{D}(\Omega)^d | \text{div } y = 0\}$, $\mathcal{V}(\Omega)$ = closure of $\mathcal{V}(\Omega)$ in $H_0^1(\Omega)^d = \{y \in H_0^1(\Omega)^d | \text{div } y = 0\}$ (since $\partial\Omega$ is lipschitzian), we have the following weak formulation of the stationary Navier-Stokes equation with Dirichlet boundary (non-slip) conditions:

$$ \int_{\Omega} \left( \nu \sum_{i,j=1}^{d} \frac{\partial y_j}{\partial x_i} \frac{\partial v_j}{\partial x_i} + \sum_{i,j=1}^{d} y_i \frac{\partial y_j}{\partial x_i} v_j \right) dx = \int_{\Omega} \left( \sum_{j=1}^{d} f_j v_j \right) dx, \forall v \in \mathcal{V}(\Omega) \quad (2.1) $$

or (see [16]), $\nu((y,v))_{\Omega} + b_{\Omega}(y, y, v) = \int_{\Omega} f \cdot v \, dx$. Here $f = (f_1, \dots, f_d) \in H^{-1}(D)^d$, and $\nu > 0$ is the viscosity.

To this equation we associate the minimization problem

$$ \min\{J(\Omega) = \int_{E} \|y - y_0\|_{e}^{2} dx, \ E \subset \Omega \in \mathcal{O}; y \text{ verifies } (2.1)\} \quad (2.2) $$

$E \subset \Omega$ is a fixed set and $y_0 \in L^2(E)^d$ is given.

This functional is a particular case of a larger class,

$$ J(\Omega) = \int_{\Lambda} j(x, y(x), \nabla y(x)) dx, \quad \Lambda = E \text{ or } \Lambda = \Omega, $$

that are studied in [6].

The uniform Lipschitz assumption turns $\mathcal{O}$ into a compact with respect to the Hausdorff-Pompeiu complementary metric (see [11], p. 466). Based on this, it is inferred in [6] that if there exists an admissible $\hat{\Omega}$ and a corresponding solution of (2.1) for which $J(\hat{\Omega})$ is finite then there exists at least an optimal pair $[\Omega^*, y^*] \in \mathcal{O} \times V(\Omega^*)$. So, the optimization problem is well defined but its solution is generally nonunique.

If $X(D) \subset C(\bar{D})$ is a functional space, define, for $g \in X(D)$, $\Omega = \Omega_g = \text{int}\{x \in D | g(x) \ge 0\}$. If $E \subset \Omega$ is to be fulfilled one must require $g(x) \ge 0$ $\forall x \in E$. $g$ is called a parametrization of $\Omega_g$ and $\Omega_g$ an admissible domain. The solutions of (2.1) in $\Omega_g$ will be denoted as $y_g$.
If $H: \mathbf{R} \to \{0,1\}$ is the Heaviside function, $H \circ g = \chi_{\bar{\Omega}_g}$, the characteristic function of $\bar{\Omega}_g$. For $\varepsilon > 0$ the following smoothing of the Yosida approximation of the maximal monotone extension of $H$ will be used:

$$H^{\varepsilon}(r) = \begin{cases} 1, & r \ge 0 \\ \frac{(\varepsilon - 2r)(r + \varepsilon)^2}{\varepsilon^2}, & -\varepsilon < r < 0 \\ 0, & r \le -\varepsilon \end{cases}$$

(see also [8], [12]). It is easy to see that $H^{\varepsilon} \in C^1(\mathbf{R})$ and is lipschitzian. The boundary value problem (2.1) has an approximate extension

$$\nu((y_{\varepsilon}, v))_D + b_D(y_{\varepsilon}, y_{\varepsilon}, v) + \frac{1}{\varepsilon} \int_D [1 - H^{\varepsilon}(g)] y_{\varepsilon} \cdot v \, dx = \int_D f \cdot v \, dx \quad (2.3)$$

Suppose now that $d = 3$. It is proved in [6] that, for $C_1 = 9m(D)^{1/6}$, if

$$\nu^2 > C_1 \|f\|_{V^*} \quad (2.4)$$

then (2.3) has an unique solution $y_{\varepsilon}(g)$ that depends continuously on $g$ as a function from $(C(D), \|\cdot\|_{\infty})$ to $L^2(D)^3$. The following theorem that is proved in [6], §3, allows the approximation of the shape optimization problem (2.1), (2.2) by the optimal control problem (2.2), (2.3).

**Theorem 2.1.** If (2.4) holds then there exists a sequence $\varepsilon_n \to 0$ such that $y_{\varepsilon_n}(g)|_{\Omega_g} \to y_g$ weakly in $H^1(\Omega_g)^3$ and strongly in $L^2(\Omega_g)^3$.

## 3 The directional derivative and a gradient type algorithm

In order to solve the optimal control problem (2.2), (2.3) through a gradient type algorithm an important step is the calculation of the directional derivative of the mapping $g \mapsto J[y_{\varepsilon}(g)]$ in the direction $w \in X(D)$. It is proved in [6], §4, that, under the uniqueness condition (2.4), this derivative in direction $w$, $\frac{\partial y_{\varepsilon}}{\partial w}(g) = (z_1, z_2, z_3) \in V(D)$ is the solution of the equation in variations

$$\int_D \left( \nu \sum_{i,j=1}^3 \frac{\partial z_j}{\partial x_i} \frac{\partial v_j}{\partial x_i} + \sum_{i,j=1}^3 y_{\varepsilon,i} \frac{\partial z_j}{\partial x_i} v_j + \sum_{i,j=1}^3 z_i \frac{\partial y_{\varepsilon,j}}{\partial x_i} v_j \right) dx + \frac{1}{\varepsilon} \int_D [1 - H^{\varepsilon}(g)] z \cdot v \, dx = \frac{1}{\varepsilon} \int_D ((H^{\varepsilon})'(g)w) y_{\varepsilon} \cdot r \, dx \quad (3.1)$$

It is also proved in [6], §4, that under condition (2.4), equation (3.1) has an unique solution.

For the optimal control problem (2.2), (2.3) the co-state equation (see [2], [5], [14]) is
$$ \int_D \left( \nu \sum_{i,j=1}^{3} \frac{\partial p_{\varepsilon,j}}{\partial x_i} \frac{\partial v_j}{\partial x_i} - \sum_{i,j=1}^{3} y_{\varepsilon,i} \frac{\partial p_{\varepsilon,j}}{\partial x_i} v_j + \sum_{i,j=1}^{3} \frac{\partial y_{\varepsilon,j}}{\partial x_i} p_{\varepsilon,j} v_j \right) dx + \frac{1}{\varepsilon} \int_D [1 - H^{\varepsilon}(g)] p_{\varepsilon} \cdot v \, dx = \int_E (y_{\varepsilon} - y_0) \cdot v \, dx. \quad (3.2) $$

Under condition (2.4) the equation (3.2) has a unique solution $p_{\varepsilon} \in V(D)$. The algorithm will result from the following theorem

**Theorem 3.1 ([6], §4, Th.5).** The direction derivative in $g \in X(D)$ of $J[y_{\varepsilon}(g)]$ in the direction $w \in X(D)$ is given by

$$ \frac{\partial J[y_{\varepsilon}(g)]}{\partial w} = \frac{1}{\varepsilon} \int_D ((H^{\varepsilon})'(g))w y_{\varepsilon} \cdot p_{\varepsilon} dx \quad (3.3) $$

$(p_{\varepsilon}$ is the unique solution of (3.2)).

### Algorithm

**Step 0** Choose a starting parametrization $g_0$ and a positive scalar $\epsilon$. Set $k = 0$.

**Step 1** Find $y_{\varepsilon}$ the solution of (2.3).

**Step 2** Find $p_{\varepsilon}$ the solution of (3.2).

**Step 3** Set the descent direction $w_k = -y_{\varepsilon} \cdot p_{\varepsilon}$. If $\|w_k\| < tol$ stop.

**Step 4** Determine $g_{k+1} = g_k + \theta_k w_k$, $\theta_k > 0$ by means of an approximate minimization

$$ J(g_{k+1}) \approx \min_{\theta \ge 0} J(g_k + \theta w_k). $$

**Step 5** Update $k = k + 1$ and go to the **Step 1**.

For the inaccurate line search at the **Step 4**, the methods of Goldstein and Armijo were used. If we denote by $j: [0, \infty) \to \mathbb{R}$ the function $j(\theta) = J(g_k + \theta w_k)$, we determine $\theta_k > 0$ such that

$$ j(0) + (1 - \lambda) \theta_k j'(0) \le j(\theta_k) \le j(0) + \lambda \theta_k j'(0) \quad (1) $$

where $\lambda \in (0, 1/2)$.

# 4 Numerical results. Shape optimization of a nozzle

## Problem setting

We have adapted the nozzle problem from [13]. We assume that the flow in a nozzle is governed by the steady Navier-Stokes equation with prescribed traction at the inflow and outflow. The problem is to design a nozzle that gives a prescribed velocity near the exit. This kind of problem arises in rocket engine industries, in the design of a spray for applying a coating or in the manufacture of high-resolution inkjet printer.
We assume that the polyhedron [A₁ A₂ A₃ A₄ A₅ A₆ A₇] is the fixed computational domain D. The coordinates of its vertices are: A₁(H, 0), A₂(0, 0), A₃(L/2, 0), A₄(L, H - h), A₅(L + ℓ, H - h), A₆(L + ℓ, H), A₇(L, H), where L = 6, ℓ = 1, H = 3, h = 1. We denote by E the observation zone which is the rectangle [A₄ A₅ A₆ A₇].

We denote by Σ<sub>in</sub> the boundary [A₁ A₂] representing the inflow section and by Σ<sub>out</sub> the boundary [A₅ A₆] representing the outflow section. The desired fluid velocity in the observation zone E is

$$y_0 = \left( v_{out} \frac{4(H - x_2)(x_2 - H + h)}{h^2}, 0 \right), \text{ where } v_{out} = 4.$$

The fluid viscosity is μ = 1 and its density is ρ = 1.

<figure><img src="image_2.png" /><figcaption>Fig. 1. Computing domain</figcaption></figure>

Let ω ⊂ D such that ω ∩ E = ∅. We look for a connected domain Ω verifying D \ ω ⊂ Ω ⊂ D and minimizing the cost function

$$J = \frac{1}{2} \int_{E} (y_{\epsilon} - y_{0}) \cdot (y_{\epsilon} - y_{0}) \, dx.$$

The traction imposed on the inflow is (100, 0) and on the outflow it is (0, 0). We impose no-slip condition on the other boundaries, including the free boundary.

### Descent direction

In Figure 2, we show 1 − H<sup>ε</sup>(g) which is an approximation of the characteristic function of the domain D \ Ω, for a typical admissible parametrization g.

We remark that the (H<sup>ε</sup>)'(r) vanishes on ℝ, excepting for r ∈ (−ε, 0). Consequently, the zone in D, where (H<sup>ε</sup>)'(g) ≠ 0 is very narrow, see Figure 3. When ε is very small, this zone could be empty. For this season, we have taken as descent direction not −(H<sup>ε</sup>)'(g)y<sub>ε</sub> · p<sub>ε</sub> which is given by (3.3), but w<sub>d</sub> = −y<sub>ε</sub> · p<sub>ε</sub>. We recall that (H<sup>ε</sup>)'(g) ≥ 0, consequently w<sub>d</sub> = −y<sub>ε</sub> · p<sub>ε</sub> is a descent direction in view of (3.3).
<figure><img src="image_3.png" /><figcaption>Fig. 2. The value of $1 - H^\epsilon(g)$ on D for $\epsilon = 10^{-4}$.</figcaption></figure>

<table><caption>Table 1. Mesh parameters used in Figure 3</caption><thead><tr><th>mesh</th><th>no. triangles</th><th>no. vertices</th><th>mesh size</th></tr></thead><tbody><tr><td>1</td><td>828</td><td>460</td><td>0.340972</td></tr><tr><td>2</td><td>3316</td><td>1749</td><td>0.171328</td></tr><tr><td>3</td><td>7412</td><td>3842</td><td>0.146829</td></tr><tr><td>4</td><td>13168</td><td>6765</td><td>0.095475</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>Fig. 3. The zone where the derivative of the Yosida approximation of the Heaviside function is not vanishes for $\epsilon = 10^{-1}$. The mesh parameters are presented in Table 1.</figcaption></figure>
## Numerical parameters

The mesh of D has 15032 triangles and 7697 vertices. We have used the following finite elements: $P_1 + bubble$ for the velocity, $P_1$ for the pressure and for the g.

We set the penalization parameter to be $\epsilon = 0.0001$, the number of iterations for the descent algorithm to be 10 and number of iterations for the line search to be 10.

We have tested our algorithm for three initial values of g:

i) for $g(x_1, x_2) = 10^{-4}(x_2 - \frac{x_1^2}{18} - 1.5)$, the initial value of J is 0.518529 and its final value of is 0.00697268;

ii) for $g(x_1, x_2) = 10^{-4}(x_2 - \frac{x_1^2}{18} - 0.5)$, the initial value of J is 0.185802 and its final value of is 0.00466894;

iii) for $g(x_1, x_2) = 10^{-4}(x_2 - \frac{x_1^2}{18} - 1)$, the initial value of J is 0.0735282 and its final value is 0.00448260.

## Numerical results

We have obtained three different optimal shapes, see Figure 4, that means the algorithm find only local optimum. The minimum final value of the cost function among the three tests is obtained in the case iii).

We remark in Figure 4 case ii), that the zero level set of the initial g partially coincides with the zero level set of the final g. In fact, the initial g vanishes on the free boundary. Since, we impose non-slip boundary condition for $y_\epsilon$ and $p_\epsilon$ on the free boundary, the descent direction $w_d = -y_\epsilon \cdot p_\epsilon$ vanishes on the free boundary, also. Consequently, $g_{k+1}$ could have the same zero level set as $g_k$.

The fluid velocity is plotted in Figure 5. We observe that the fluid velocity is very small in the exterior of the optimal shape, more precisely we have

$$ Error(g) = \int_{\omega} (1 - H^{\epsilon}(g)) y_{\epsilon} \cdot y_{\epsilon} dx = 0.000446 $$

**Acknowledgement** This work was supported by ANCS Grant 480/17.03.2011 in the framework of the “Brâncuși” programm of cooperation between Romania and France.
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Case i): final shape (bottom) obtained from the initial shape (top). Case ii): final shape (top) obtained from the initial shape (bottom). Case iii): final shape (bottom) obtained from the initial shape (top).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Fluid velocity for the optimal shape in the case iii)</figcaption></figure>
References

1. T. Borrvall and J. Petersson (2003), *Topology optimization of fluids in Stokes flow*, International Journal for Numerical Methods in Fluids, 41, 77–107.

2. L. Dede (2007), *Optimal flow control for Navier-Stokes equations, drag minimization*, J. Numer. Meth. Fluids, vol. 44, 4, 347–366.

3. M.C. Delfour and J.P. Zolesio (2001), *Shapes and Geometrics*, Analysis, Differential Calculus and Optimization, SIAM, Philadelphia.

4. Z. Gao and Y. Ma (2007), *Optimal shape design for viscous incompressible flow*, arxiv: math.oc/0701470v1.

5. M. Gunzburger (2000), *Adjoint equation-based methods for control problems in viscous, incompressible flows*, Flow, Turbul., Comb. 65, 249–272.

6. A. Halanay, D. Tiba (2009), *Shape optimization for stationary Navier-Stokes equations*, Control and Cybernetics, 38, no.4, 1359–1375.

7. J.L. Lions (1971), *Optimal control of systems governed by partial differential equations*, Springer, Berlin.

8. R. Makinen, P. Neittaanmaki and D. Tiba (1992), *On a fixed domain approach for shape optimization problem*, In W.F. Ames and P.J. van der Houwer, editors, *Computational and Applied Mathematics II: Differential Equations*, North-Holland, Amsterdam, 317–326.

9. B. Mohammadi and O. Pironneau (2001), *Applied Shape Optimization for Fluids*, Oxford University Press, New York.

10. P. Neittaanmaki and D. Tiba (1994), *Optimal control of nonlinear parabolic systems*, Theory, algorithms and applications, vol. 179 of Monographs and Textbooks in Pure and Applied Mathematics, Marcel Dekker, New York.

11. P. Neittaanmaki, J. Sprekels and D. Tiba (2006), *Optimization of elliptic systems. Theory and applications*, Springer, New York.

12. P. Neittaanmaki, A. Pennanen and D. Tiba (2009), *Fixed domain approaches in shape optimization problems with Dirichlet boundary conditions*, J. of Inverse Problems, 25, 1–18.

13. O. Pironneau (1984), *Optimal shape design for elliptic systems*, Springer, Berlin.

14. M. Posta and T. Roubicek (2007), *Optimal control of Navier-Stokes equations by Oseen approximations*, Preprint 2007–013, Necas Center, Prague.

15. T. Roubicek and F. Troltzsch (2003), *Lipschitz stability of optimal control for steady-state Navier-Stokes equations*, Control and Cybernetics, 32, 683–705.

16. R. Temam (1979), *Navier-Stokes equations. Theory and numerical analysis*, North-Holland, Amsterdam.