Geometrical shape optimization in fluid mechanics using FreeFem++

Abstract

In this article, we present simple and robust numerical methods for two-dimensional geometrical shape optimization problems, in the context of viscous flows driven by the stationary Navier-Stokes equations at low Reynolds number. The salient features of our algorithm are exposed with an educational purpose; in particular, the numerical resolution of the nonlinear stationary Navier-Stokes system, the Hadamard boundary variation method for calculating the sensitivity of the minimized function of the domain, and the mesh update strategy are carefully described. Several pedagogical examples are discussed. The corresponding program is written in the FreeFem++ environment, and it is freely available. Its chief features—and notably the implementation details of the main steps of our algorithm—are carefully presented, so that it can easily be handled and elaborated upon to deal with different, or more complex physical situations.

Appendix A: Sketch of the proof of Theorem 2

The differentiability of the solution $(\mathbf {u},p)$ to the Navier-Stokes system (1) with respect to the domain is a technical, albeit quite classical matter, and we admit the result, referring for instance to Henrot and Pierre (2005) for the rigorous definition of this notion, and to Henrot and Privat (2010) or de La Sablonière et al. (2011) for this precise calculation. The derivative $(\mathbf {u}^{\prime }, p^{\prime })$ of $(\mathbf {u},p)$ with respect to the domain, in the direction $\boldsymbol {\theta } \in {\Theta }_{ad}$, is solution to the problem:

$$ \left\{ \begin{array}{cl} -\nu{\Delta} \mathbf{u}^{\prime}+(\nabla \mathbf{u} )\mathbf{u}^{\prime} + (\nabla \mathbf{u}^{\prime})\mathbf{u} + \nabla p^{\prime} = 0 & \text{in } {\Omega},\\ \text{div}(\mathbf{u}^{\prime}) = 0 & \text{in } {\Omega},\\ \sigma(\mathbf{u}^{\prime},p^{\prime})\mathbf{n} = 0 & \text{on } {\Gamma}_{\text{out}},\\ \mathbf{u}^{\prime}= 0 & \text{on } {\Gamma}_{\text{in}},\\ \mathbf{u}^{\prime}=-\left( \frac{\partial \mathbf{u}}{\partial n}\right)(\boldsymbol{\theta}.n) & \text{on } {\Gamma}. \end{array} \right. $$

(32)

Also, we only present the calculation of the shape derivative of the functional $D({\Omega })$ given by (5), the calculation being on any point easier in the case of $E({\Omega })$; see de La Sablonière et al. (2011) if need be.

Using the chain rule from the definition (5) of $D({\Omega })$ yields:

$$ D^{\prime}({\Omega})(\boldsymbol{\theta}) ={\int}_{{\Gamma}_{\text{out}}} \mathbf{u}^{\prime}\cdot(\mathbf{u}-\mathbf{u}_{\text{ref}})\,ds. $$

(33)

The main idea of the proof consists in using the adjoint state $(\mathbf {v}_{d},q_{d})$, solution to (12): performing several integrations by parts allows to eliminate the unknown derivatives $(\mathbf {u}^{\prime },p^{\prime })$ from the expression (33). More precisely, multiplying the first equation in (32) by $\mathbf {v}_{d}$ and integrating by parts yields

$$\begin{array}{@{}rcl@{}} 0&=& {\int}_{{\Omega}} (-\nu{\Delta} \mathbf{u}^{\prime}+(\nabla \mathbf{u} )\mathbf{u}^{\prime} + (\nabla \mathbf{u}^{\prime})\mathbf{u} + \nabla p^{\prime}) \cdot \mathbf{v}_{d}\,dx\\ &=&{\int}_{{\Omega}}(2\nu \:e(\mathbf{u}^{\prime}):e(\mathbf{v}_{d}) - \text{div}(\mathbf{v}_{d}) p^{\prime} + (\nabla \mathbf{u})\mathbf{u}^{\prime}\cdot\mathbf{v}_{d}\\ &&+ (\nabla \mathbf{u}^{\prime})\mathbf{u}\cdot\mathbf{v}_{d})\,dx - {\int}_{\partial{\Omega}}\sigma (\mathbf{u}^{\prime},p^{\prime})\mathbf{n}\cdot\mathbf{v}_{d}\,ds\\ &=& {\int}_{\Omega}{(2\nu\: e(\mathbf{u}^{\prime}):e(\mathbf{v}_{d}) + (\nabla \mathbf{u})\mathbf{u}^{\prime}\cdot\mathbf{v}_{d} + (\nabla \mathbf{u}^{\prime})\mathbf{u}\cdot\mathbf{v}_{d})\,dx}\\ \end{array} $$

(34)

where the boundary integral has vanished thanks to the boundary conditions satisfies by $(\mathbf {u}^{\prime },p^{\prime })$ and $(\mathbf {v}_{d},q_{d})$. Likewise, multiplying the first equation in (12) by $\mathbf {u}^{\prime }$ and integrating by parts, we obtain:

$$\begin{array}{@{}rcl@{}} 0&=& {\int}_{{\Omega}} \left( -\nu{\Delta} \mathbf{v}_{d}+(\nabla \mathbf{u})^{T}\mathbf{v} -(\nabla \mathbf{v}_{d} )\mathbf{u} + \nabla q_{d}\right) \cdot \mathbf{u}^{\prime} \, dx \\ &=& {\int}_{\Omega}{ \left( 2\nu\: e(\mathbf{u}^{\prime}): e(\mathbf{v}_{d}) +(\nabla\mathbf{u}) \mathbf{u}^{\prime}\cdot\mathbf{v}_{d} -(\nabla\mathbf{v}_{d})\mathbf{u}\cdot\mathbf{u}^{\prime}\right)\, dx}\\ &&-{\int}_{\partial{\Omega}}\sigma (\mathbf{v}_{d},q_{d})\mathbf{n}\cdot\mathbf{u}^{\prime}\,ds. \end{array} $$

(35)

Combining Eqs. 34 and (35) leads to:

$$ -{\int}_{\Omega}((\nabla\mathbf{u}^{\prime})\mathbf{u}\cdot\mathbf{v}_{d}+ (\nabla\mathbf{v}_{d})\mathbf{u}\cdot\mathbf{u}^{\prime})\,dx = {\int}_{\partial{\Omega}} \sigma (\mathbf{v}_{d},q_{d})\mathbf{n}\,ds. $$

(36)

Now using the identity

$$\begin{array}{@{}rcl@{}} {\int}_{\Omega}{(\nabla\mathbf{v}_{d})\cdot \mathbf{u} \cdot\mathbf{u}^{\prime}\,dx} &=& {\int}_{\partial{\Omega}}{(\mathbf{v}_{d}\cdot\mathbf{u}^{\prime})(\mathbf{u}\cdot\mathbf{n})\,ds}\\ &&- {\int}_{\Omega}{(\nabla\mathbf{u}^{\prime})\mathbf{u}\cdot\mathbf{v}_{d}\,dx}, \end{array} $$

(37)

which again follows from integration by parts, (36) rewrites:

$$ {\int}_{\partial{\Omega}} {(\sigma(\mathbf{v}_{d},q_{d})\mathbf{n}\cdot\mathbf{u}^{\prime} + (\mathbf{u}^{\prime}\cdot\mathbf{v}_{d})(\mathbf{u}\cdot\mathbf{n}))\,ds} = 0. $$

(38)

Eventually, taking into account the boundary conditions in the systems (1), (10) and (32) yields:

$$\begin{array}{@{}rcl@{}} D^{\prime}({\Omega})(\boldsymbol{\theta}) &=&{\int}_{{\Gamma}_{\text{out}}} (\mathbf{u}-\mathbf{u}_{\text{ref}})\cdot\mathbf{u}^{\prime}\,ds\\ &=& {\int}_{{\Gamma}_{\text{out}}}{(\sigma (\mathbf{v}_{d},q_{d})\mathbf{n} + (\mathbf{u} \cdot \mathbf{n})\mathbf{v}_{d}) \cdot \mathbf{u}^{\prime}\,ds}\\ &=&- {\int}_{{\Gamma}}{(\sigma (\mathbf{v}_{d},q_{d})\mathbf{n} + (\mathbf{u} \cdot \mathbf{n})\mathbf{v}_{d}) \cdot \mathbf{u}^{\prime}\,ds}\\ &=& {\int}_{{\Gamma}}{(\sigma (\mathbf{v}_{d},q_{d})\mathbf{n} + (\mathbf{u} \cdot \mathbf{n})\mathbf{v}_{d}) \cdot \frac{\partial\mathbf{u}}{\partial n} \:\boldsymbol{\theta} \cdot n\,ds}.\\ \end{array} $$

(39)

We now use the boundary conditions $\mathbf {u} = 0$ and $\mathbf {v}_{d} = 0$ on ${\Gamma }$ to simplify this last expression. For any tangential vector field $\boldsymbol {\tau }: {\Gamma } \to \mathbb {R}^{d}$ to ${\Gamma }$, they imply that $\frac {\partial \mathbf {u}}{\partial \boldsymbol {\tau }} = 0$, and so, using that $\text {div}(\mathbf {u}) = 0$,

$$ \frac{\partial \mathbf{u}}{\partial n} \cdot \mathbf{n} = 0 $$

(40)

the same relation holds for $\mathbf {v}_{d}$. Hence (39) rewrites:

$$D^{\prime}({\Omega})(\boldsymbol{\theta}) = {\int}_{{\Gamma}}{2\nu \: e(\mathbf{v}_{d})\mathbf{n} \cdot \frac{\partial\mathbf{u}}{\partial n} \:\boldsymbol{\theta} \cdot \mathbf{n}\,ds}.$$

After a few algebraic manipulations based again on (40), we eventually obtain:

$$ D^{\prime}({\Omega})(\boldsymbol{\theta}) = {\int}_{\Gamma}{ 2\nu \: e(\mathbf{u}):e(\mathbf{v}_{d})(\boldsymbol{\theta}\cdot\mathbf{n})\,ds}, $$

(41)

which is the desired result, and terminates the proof of Theorem 2.