# Shape Sensitivity Analysis of Incompressible Non-Newtonian Fluids

Jan Sokolowski, Jan Stebel

► To cite this version:

Jan Sokolowski, Jan Stebel. Shape Sensitivity Analysis of Incompressible Non-Newtonian Fluids. 25th System Modeling and Optimization (CSMO), Sep 2011, Berlin, Germany. pp.427-436, 10.1007/978-3-642-36062-6_43. [hal-01347565](https://hal-inriaHal) .

HAL Id: hal-01347565

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

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
# Shape sensitivity analysis of incompressible non-Newtonian fluids

Jan Sokołowski¹ and Jan Stebel²

¹ Institut Élie Cartan, UMR 7502 (Nancy Université, CNRS, INRIA), Laboratoire de Mathématiques, Université Henri Poincaré Nancy 1, B.P. 239, 54506 Vandoeuvre Lès Nancy Cedex, France, Jan.Sokołowski@iecn.u-nancy.fr

² Institute of Mathematics of the Academy of Sciences of the Czech Republic, Žitná 25, 115 67 Praha 1, Czech Republic, stebel@math.cas.cz

**Abstract.** We study the shape differentiability of a cost function for the steady flow of an incompressible viscous fluid of power-law type. The fluid is confined to a bounded planar domain surrounding an obstacle. For smooth perturbations of the shape of the obstacle we express the shape gradient of the cost function which can be subsequently used to improve the initial design.

**Keywords:** shape optimization, shape gradient, incompressible fluid, non-Newtonian fluid, Navier-Stokes equations

## 1 Introduction

Shape optimization for nonlinear partial differential equations is a growing field in the contemporary optimum design of structures. In this field systems of the solid and fluid mechanics as well as e.g., the coupled models of fluid-structure interaction are included for real life problems. The main difficulty associated with the mathematical analysis of nonlinear state equations is the lack of existence of global strong solutions for mathematical models in three spatial dimensions.

In numerical methods of shape optimization the common approach is the discretization of continuous shape gradient. Therefore, the proper derivation and analysis of the regularity properties of the shape gradient is crucial for numerical solution of the shape optimization problem. The shape sensitivity analysis requires, in particular, the proof of the Lipschitz continuity of solutions to the state equations with respect to the boundary variations. This property of solutions can be obtained e.g. by analysis of the state equation transported to the fixed reference domain which is explained in the case of linear elliptic boundary value problems in monograph [11]. For the nonlinear problems the Lipschitz continuity is not obvious and it requires the additional regularity of solutions to the state equation. In addition, for the applications of levelset method of shape optimization it is required that the obtained shape gradient of the cost functional is given by a function while the general theory gives only the existence
of a distribution. In conclusion, it seems that the shape sensitivity analysis in the case of a nonlinear state equation is the main step towards the numerical solution of the shape optimization problems.

In gas dynamics described by the compressible Navier-Stokes there is the existence of weak global solutions. However, the shape sensitivity analysis can be performed only for specific local solutions. The state of art in shape optimization for compressible Navier-Stokes equations is presented in the forthcoming monograph [8], see also [7]. For incompressible Navier-Stokes equations, the sensitivity analysis of shape functionals is performed e.g. in [2] and [6]. In this paper we are concerned with the non-Newtonian model where the stress is a (nonlinear) function of the velocity gradient. Optimal control problem for this model was studied in [9, 13]. Numerical shape optimization was done in [1], see also [3]. We present new results on the existence of the shape gradient.

We consider the steady flow of an incompressible fluid in a bounded domain $\Omega := B \setminus S$ in $\mathbb{R}^2$, where $B$ is a container and $S$ is an obstacle. Motion of the fluid is described by the system of equations

$$ 
\begin{array}{lll}
\operatorname{div}(\mathbf{v} \otimes \mathbf{v}) - \operatorname{div}\mathcal{S}(\mathbb{D}\mathbf{v}) + \nabla p + \mathbb{C}\mathbf{v} = \mathbf{f} & \text{in } \Omega, \\
\operatorname{div}\mathbf{v} = 0 & \text{in } \Omega, & (P(\Omega)) \\
\mathbf{v} = \mathbf{g} & \text{on } \partial\Omega.
\end{array}
$$

Here $\mathbf{v}$, $p$, $\mathbb{C}$, $\mathbf{f}$ stands for the velocity, the pressure, the constant skew-symmetric Coriolis tensor and the body force, respectively. The traceless part $\mathcal{S}$ of the Cauchy stress can depend on the symmetric part $\mathbb{D}\mathbf{v}$ of the velocity gradient in the following way:

$$ \mathcal{S}(\mathbb{D}\mathbf{v}) = \nu(|\mathbb{D}\mathbf{v}|^2)\mathbb{D}\mathbf{v}, \qquad (1) $$

where $\nu$, $|\mathbb{D}\mathbf{v}|^2$ is the viscosity and the shear rate, respectively. In particular, we assume that $\nu$ has a polynomial growth (see Section 2.1 below), which includes e.g. the Carreau and the power-law model.

In the model the term of Coriolis type is present. This term appears e.g. when the change of variables is performed in order to take into account the flight scenario of the obstacle in the fluid.

The aim of this paper is to investigate differentiability of a shape functional depending on the solution to $(P(\Omega))$ with respect to the variations of the shape of the obstacle. We consider a model problem with the drag functional

$$ J(\Omega) := \int_{\partial S} (\mathcal{S}(\mathbb{D}\mathbf{v}) - p\mathbb{I})\mathbf{n} \cdot \mathbf{d}, \qquad (2) $$

with a given constant unit vector $\mathbf{d}$. Instead of $J$ one could take other type of functional, since our method does not rely on its specific form.

Our main interest is the rigorous analysis of the shape differentiability for $(P(\Omega))$ and (2). We follow the general framework developed in [11] using the speed method and the notion of the material derivative. Let us point out that due to (1) the state problem is nonlinear in its nature. We refer the reader to [12] for an introduction to optimization problems for nonlinear partial differential equations.
## 1.1 Shape derivatives

We start by the description of the framework for the shape sensitivity analysis. For this reason, we introduce a vector field $\mathbf{T} \in C^2(\mathbb{R}^2, \mathbb{R}^2)$ vanishing in the vicinity of $\partial B$ and define the mapping

$$ \mathbf{y}(\mathbf{x}) = \mathbf{x} + \varepsilon \mathbf{T}(\mathbf{x}). $$

For small $\varepsilon > 0$ the mapping $\mathbf{x} \mapsto \mathbf{y}(\mathbf{x})$ takes diffeomorphically the region $\Omega$ onto $\Omega_{\varepsilon} = B \setminus S_{\varepsilon}$ where $S_{\varepsilon} = \mathbf{y}(S)$. We consider the counterpart of problem $(P(\Omega))$ in $\Omega_{\varepsilon}$, with the data $\mathbf{f}|_{\Omega_{\varepsilon}}$ and $\mathbf{g}|_{\Omega_{\varepsilon}}$. The new problem will be denoted by $(P(\Omega_{\varepsilon}))$ and its solution by $(\bar{\mathbf{v}}_{\varepsilon}, \bar{\mathbf{p}}_{\varepsilon})$.

For the nonlinear system $(P(\Omega))$ we introduce the shape derivatives of solutions. To this end we need the linearized system of the form:

Find the couple $(\mathbf{u}, \pi)$ such that

$$ 
\begin{array}{lll}
\mathrm{div} \left[ \mathbf{u} \otimes \mathbf{v} + \mathbf{v} \otimes \mathbf{u} - S'(\mathbb{D}\mathbf{v})\mathbb{D}\mathbf{u} \right] + \nabla\pi + \mathbb{C}\mathbf{u} = \mathbf{F} & \text{in } \Omega, \\
\mathrm{div} \mathbf{u} = 0 & \text{in } \Omega, & (P_{\mathrm{lin}}(\Omega)) \\
\mathbf{u} = \mathbf{h} & \text{on } \partial\Omega,
\end{array}
$$

where $\mathbf{F}$ and $\mathbf{h}$ are given elements.

The shape derivative $\mathbf{v}'$ and the material derivative $\dot{\mathbf{v}}$ of solutions are formally introduced by

$$ \mathbf{v}' := \lim_{\varepsilon \to 0} \frac{\bar{\mathbf{v}}_{\varepsilon} - \mathbf{v}}{\varepsilon}, \quad \dot{\mathbf{v}} := \lim_{\varepsilon \to 0} \frac{\bar{\mathbf{v}}_{\varepsilon} \circ \mathbf{y} - \mathbf{v}}{\varepsilon}, $$

where $\bar{\mathbf{v}}_{\varepsilon} \circ \mathbf{y}(\mathbf{x}) := \bar{\mathbf{v}}_{\varepsilon}(\mathbf{y}(\mathbf{x}))$, and are related to each other as follows:

$$ \dot{\mathbf{v}} = \mathbf{v}' + (\nabla \mathbf{v})\mathbf{T}. $$

The standard calculus for differentiating with respect to shape yields that $\mathbf{v}'$ is the solution of $(P_{\mathrm{lin}}(\Omega))$ with the data $\mathbf{F} = \mathbf{0}$ and $\mathbf{h} = -\partial\mathbf{v}/\partial\mathbf{n}(\mathbf{T} \cdot \mathbf{n})$. Using (7) as the definition of $J$ we obtain the expression for the shape gradient:

$$ 
\begin{aligned}
dJ(\Omega; \mathbf{T}) &:= \lim_{\varepsilon \to 0} \frac{J(\Omega_{\varepsilon}) - J(\Omega)}{\varepsilon} \\
&= \int_{\Omega} [ (\mathbb{C}\mathbf{v}') \cdot \boldsymbol{\xi} + (S'(\mathbb{D}\mathbf{v})\mathbb{D}\mathbf{v}' - \mathbf{v}' \otimes \mathbf{v} - \mathbf{v} \otimes \mathbf{v}') : \nabla\boldsymbol{\xi} ] - \int_{\partial S} (\mathbf{f} \cdot \mathbf{d})\mathbf{T} \cdot \mathbf{n}. \quad (3)
\end{aligned}
$$

In the above formula, the part containing $\mathbf{v}'$ depends implicitly on the direction $\mathbf{T}$. This is not convenient for practical use, hence we introduce the adjoint problem for further simplification of (3):

Find the couple $(\mathbf{w}, s)$ such that

$$ 
\begin{array}{lll}
-2(\mathbb{D}\mathbf{w})\mathbf{v} - \mathrm{div} [S'(\mathbb{D}\mathbf{v})^{\top}\mathbb{D}\mathbf{w}] + \nabla s - \mathbb{C}\mathbf{w} = \mathbf{0} & \text{in } \Omega, \\
\mathrm{div} \mathbf{w} = 0 & \text{in } \Omega, \\
\mathbf{w} = \mathbf{d} & \text{on } \partial\Omega.
\end{array}
\quad (P_{\mathrm{adj}}(\Omega))
$$
Consequently, the expression for $dJ$ reduces to

$$dJ(\Omega; \mathbf{T}) = - \int_{\partial S} \left[ (\mathcal{S}'(\mathbb{D}\mathbf{v})^{\top} \mathbb{D}\mathbf{w} - s\mathbb{I}) : \frac{\partial \mathbf{v}}{\partial \mathbf{n}} \otimes \mathbf{n} + \mathbf{f} \cdot \mathbf{d} \right] \mathbf{T} \cdot \mathbf{n}. \quad (4)$$

In order to prove the result given by (3) and (4) we need the material derivatives. In particular, it is sufficient to show that the linear mapping

$$\mathbf{T} \mapsto dJ(\Omega; \mathbf{T})$$

is continuous in an appropriate topology, see the structure Theorem in the book [11] for details.

## 2 Preliminaries

We impose the structural assumptions on the data, state the known results on well-posedness of $(P(\Omega))$ and introduce the elementary notation for shape sensitivity analysis.

### 2.1 Structural assumptions

We require that $\mathcal{S}$ has a potential $\Phi: [0, \infty) \to [0, \infty)$, i.e. $\mathcal{S}_{ij}(\mathbb{D}) = \partial\Phi(|\mathbb{D}|^2)/\partial\mathbb{D}_{ij}$. Further we assume that $\Phi$ is a $\mathcal{C}^3$ function with $\Phi(0) = 0$ and that there exist constants $C_1, C_2, C_3 > 0$ and $r \ge 2$ such that

$$C_1(1 + |\mathbf{A}|^{r-2})|\mathbb{B}|^2 \le \mathcal{S}'(\mathbf{A}) :: (\mathbb{B} \otimes \mathbb{B}) \le C_2(1 + |\mathbf{A}|^{r-2})|\mathbb{B}|^2, \quad (5a)$$

$$|\mathcal{S}''(\mathbf{A})| \le C_3(1 + |\mathbf{A}|^{r-3}) \quad (5b)$$

for any $0 \ne \mathbf{A}, \mathbb{B} \in \mathbb{R}_{sym}^{2 \times 2}$. Here the symbol $:\:$ stands for the usual scalar product in $\mathbb{R}^4$. The above inequalities imply the monotone structure of $\mathcal{S}$, see e.g. [5].

### 2.2 Weak formulation

For the definition of the weak solution we will use the space

$$\mathbf{W}_{0,\mathrm{div}}^{1,r}(\Omega) := \{\phi \in \mathbf{W}_{0}^{1,r}(\Omega); \mathrm{div} \, \phi = 0\}.$$

Let $\mathbf{f} \in (\mathbf{W}_{0,\mathrm{div}}^{1,2}(\Omega))^*$ and $\mathbf{g} \in \mathbf{W}^{1,r}(\Omega)$ with $\mathrm{div} \, \mathbf{g} = 0$. Then a function $\mathbf{v} \in \mathbf{g} + \mathbf{W}_{0,\mathrm{div}}^{1,r}(\Omega)$ is said to be a weak solution to the problem $(P(\Omega))$ if

$$\int_{\Omega} [\mathcal{S}(\mathbb{D}\mathbf{v}) : \mathbb{D}\phi - \mathbf{v} \otimes \mathbf{v} : \nabla\phi + \mathbb{C}\mathbf{v} \cdot \phi] = \int_{\Omega} \mathbf{f} \cdot \phi \quad (6)$$

for every $\phi \in \mathbf{W}_{0,\mathrm{div}}^{1,r}(\Omega)$. Note that the pressure is eliminated since test functions are divergence free.

The following result was shown in [4].
**Theorem 1 (Kaplický et al. [4]).** Let $\Omega \in C^2$, $f \in L^{2+\epsilon_0}(\Omega)$, $\epsilon_0 > 0$ and (5a)-(5b) hold with $r > \frac{3}{2}$. Then there exists a constant $\delta > 0$ such that for every $g$ satisfying

$$ \|g\|_{3,q} \le \delta, \quad (q > 2), $$

problem $(P(\Omega))$ has a weak solution satisfying $v \in W^{2,2+\epsilon}(\Omega)$, $p \in W^{1,2+\epsilon}(\Omega)$, $\epsilon > 0$.

Note that the above result applies only to the unperturbed domain, i.e. $\epsilon = 0$. Assuming smallness of $\|f\|_{2,B}$ and $\|g\|_{3,q,B}$, one can prove that $(P(\Omega))$, $(P(\Omega_\epsilon))$ has a unique weak solution satisfying

$$ \|v\| \le C_E(\|f\|_{2,B}, \|g\|_{3,q,B}) \quad \text{and} \quad \|\bar{v}_\epsilon\| \le C_E(\|f\|_{2,B}, \|g\|_{3,q,B}), $$

respectively, where $C_E$ is independent of $\epsilon$. At this point we summarize the main hypotheses.

**Assumption 1** In what follows, $\Omega \in C^2$ is a bounded planar domain of the form $\Omega = B \setminus S$, $f \in L^{2+\epsilon_0}(B)$, $\epsilon_0 > 0$, $g \in W^{3,q}(B)$ ($q > 2$) is supported in the vicinity of $\partial B$, (5a)-(5b) hold with $r \in [2, 4)$ and $\|f\|_{2,B}$, $\|g\|_{3,q,B}$ are small enough.

Let us point out that equation (2) which defines $J$ is not suitable for weak solutions in general, since the energy inequality does not provide enough information about the trace of $p$ and $Dv$. We therefore introduce an alternative definition that requires less regularity. Let us fix an arbitrary divergence free function $\xi \in C_c^\infty(B, \mathbb{R}^2)$ such that $\xi = d$ in a vicinity of $S$. Then, integrating (2) by parts and using $(P(\Omega))$ yields:

$$ J(\Omega) = \int_{\Omega} [(Cv - f) \cdot \xi + (S(Dv) - v \otimes v) : \nabla \xi]. \quad (7) $$

Note that this identity is finite for any $v \in W^{1,2}(\Omega)$.

## 2.3 Deformation of the shape

Let us introduce the following notation: We will denote by $DT$ the Jacobian matrix whose components are $(DT)_{ij} = (\nabla T)_{ji} = \partial_i T_j$. Further,

$$ N(\mathbf{x}) := g(\mathbf{x})M^{-1}(\mathbf{x}), \quad M(\mathbf{x}) := I + \varepsilon D T(\mathbf{x}), \quad g(\mathbf{x}) := \det M(\mathbf{x}). $$

One can easily check that the matrix $N$ and the determinant $g$ admit the expansions:

$$ g = 1 + \varepsilon \mathrm{div} \, T + O(\varepsilon^2), \quad N = I + \varepsilon N' + O(\varepsilon^2), \quad N' = (\mathrm{div} \, T)I - DT, \quad (8) $$

where the symbol $O(\varepsilon^2)$ denotes a function whose norm in $C^1(\bar{\Omega})$ is bounded by $C\varepsilon^2$, see [11].
The value of the shape functional for $\Omega_{\varepsilon}$ is given by

$$J(\Omega_{\varepsilon}) := \int_{\Omega_{\varepsilon}} [(C\bar{\mathbf{v}}_{\varepsilon} - \mathbf{f}) \cdot \boldsymbol{\xi}_{\varepsilon} + (S(D\bar{\mathbf{v}}_{\varepsilon}) - \bar{\mathbf{v}}_{\varepsilon} \otimes \bar{\mathbf{v}}_{\varepsilon}) : \nabla \boldsymbol{\xi}_{\varepsilon}]$$

where $\boldsymbol{\xi}_{\varepsilon} := (\mathbb{N}^{-\top}\boldsymbol{\xi}) \circ \mathbf{y}^{-1}$. Using the properties of the Piola transform one can check that $\mathrm{div} \boldsymbol{\xi}_{\varepsilon} = 0$. If $\bar{\mathbf{v}}_{\varepsilon}$ and $\bar{p}_{\varepsilon}$ were sufficiently smooth, it would hold that

$$J(\Omega_{\varepsilon}) = \int_{\partial S_{\varepsilon}} (S(D\bar{\mathbf{v}}_{\varepsilon}) - \bar{p}_{\varepsilon}\mathbb{I})\mathbf{n}_{\varepsilon} \cdot \mathbf{d}. \quad (9)$$

Nevertheless, as opposed to $(P(\Omega))$, we do not require any additional regularity of the solution to the perturbed problem $(P(\Omega_{\varepsilon}))$ and hence the expression in (9) need not be well defined.

We introduce the auxiliary function $\tilde{\mathbf{v}}$:

$$\tilde{\mathbf{v}} := \lim_{\varepsilon \to 0} \frac{\mathbb{N}^{\top}\bar{\mathbf{v}}_{\varepsilon} \circ \mathbf{y} - \mathbf{v}}{\varepsilon},$$

which is related to the material derivative $\dot{\mathbf{v}}$ by the identity

$$\tilde{\mathbf{v}} = \mathbb{N}'^{\top}\mathbf{v} + \dot{\mathbf{v}}.$$

For the justification of the results of the paper we will use $\tilde{\mathbf{v}}$ since, unlike the material derivative, it preserves the divergence free condition.

## 3 Main results

The first result is the existence of $\tilde{\mathbf{v}}$ and hence also of the material derivative.

**Theorem 2.** Let Assumption 1 be satisfied. Then the function $\tilde{\mathbf{v}}$ exists and is the unique weak solution of $(P_{\mathrm{lin}}(\Omega))$ with the data

$$\mathbf{F} = \mathbf{A}'_0 := \mathrm{div}(\mathbf{v} \otimes \mathbb{N}'^{\top}\mathbf{v}) + \mathbb{N}'\mathrm{div}(\mathbf{v} \otimes \mathbf{v}) \\ + \mathrm{div} [S'(D\mathbf{v}) ((\mathbb{N}' - \mathbb{I}\mathrm{tr}\mathbb{N}')\nabla\mathbf{v})_{\mathrm{sym}} - D(\mathbb{N}'^{\top}\mathbf{v})) + \mathbb{N}'^{\top}S(D\mathbf{v})] \\ - \mathbb{N}'\mathrm{div} S(D\mathbf{v}) + ((\mathbb{N}' - \mathbb{I}\mathrm{tr}\mathbb{N}')\mathbb{C} + \mathbb{C}\mathbb{N}'^{\top})\mathbf{v} + (\mathbb{I}\mathrm{tr}\mathbb{N}' - \mathbb{N}')\mathbf{f} + (\nabla\mathbf{f})\mathbf{T}, \quad (10a)$$

$$\mathbf{h} = \mathbf{0}. \quad (10b)$$

The following estimate holds:

$$\|\tilde{\mathbf{v}}\|_{1,2,\Omega} \le C \|\mathbf{A}'_0\|_{\mathbf{W}_{0,\mathrm{div}}^{1,2}(\Omega)^*} \le C \|\mathbf{T}\|_{C^2(\bar{\Omega})}. \quad (11)$$

Next we establish the existence of the shape gradient of $J$.
**Theorem 3.** Let Assumption 1 be satisfied and $\mathbf{f} \in \mathbf{W}^{1,2}(\Omega)$. Then the shape gradient of $J$ reads

$$dJ(\Omega, \mathbf{T}) = J_{\mathbf{v}}(\tilde{\mathbf{v}}) + J_e(\mathbf{T}),$$

where the dynamical part $J_{\mathbf{v}}$ and the geometrical part $J_e$ is given by

$$J_{\mathbf{v}}(\tilde{\mathbf{v}}) = \int_{\Omega} [(C\tilde{\mathbf{v}}) \cdot \boldsymbol{\xi} + (S'(D\mathbf{v})D\tilde{\mathbf{v}} - \tilde{\mathbf{v}} \otimes \mathbf{v} - \mathbf{v} \otimes \tilde{\mathbf{v}}) : \nabla \boldsymbol{\xi}],$$

$$J_e(\mathbf{T}) = \int_{\Omega} \left\{ \begin{aligned} & [(\mathrm{I} \operatorname{tr} \mathbf{N}' - \mathbf{N}') C\mathbf{v} - C\mathbf{N}'^{\top}\mathbf{v} - (\mathrm{I} \operatorname{tr} \mathbf{N}' - \mathbf{N}') \mathbf{f} - (\nabla \mathbf{f})\mathbf{T}] \cdot \boldsymbol{\xi} \\ & + [\mathbf{v} \otimes \mathbf{N}'^{\top}\mathbf{v} + S'(D\mathbf{v}) ((\mathbf{N}'\nabla\mathbf{v} - \nabla(\mathbf{N}'^{\top}\mathbf{v}))_{\mathrm{sym}} - (\mathrm{tr} \mathbf{N}')D\mathbf{v}) + \mathbf{N}'^{\top}S(D\mathbf{v})] : \nabla \boldsymbol{\xi} \\ & + [\mathbf{v} \otimes \mathbf{v} - S(D\mathbf{v})] : \nabla(\mathbf{N}'^{\top}\boldsymbol{\xi}) \end{aligned} \right\},$$

respectively. In particular, as $\tilde{\mathbf{v}}$ depends continuously on $\mathbf{T}$, the mapping

$$\mathbf{T} \mapsto dJ(\Omega, \mathbf{T})$$

is a bounded linear functional on $C^2(\mathbb{R}^2, \mathbb{R}^2)$.

Based on the previous result we can deduce that the shape gradient has the form of a distribution supported on the boundary of the obstacle. Since this representation is unique, the formal results derived in Section 1.1 are justified provided that the shape derivatives and adjoints exist and are sufficiently regular.

**Corollary 1.** Let Assumption 1 be satisfied. Then

(i) the shape derivative $\mathbf{v}'$ exists and is the unique weak solution to $(P_{\mathrm{lin}}(\Omega))$ with $\mathbf{F} = \mathbf{0}$, $\mathbf{h} = -\frac{\partial \mathbf{v}}{\partial \mathbf{n}}(\mathbf{T} \cdot \mathbf{n})$;

(ii) the adjoint problem $(P_{\mathrm{adj}}(\Omega))$ has a unique weak solution that satisfies: $\mathbf{w} \in \mathbf{W}^{2,2}(\Omega)$ and $s \in W^{1,2}(\Omega)$.

If in addition $\mathbf{f} \in \mathbf{W}^{1,2}(\Omega)$, then

(iii) the shape gradient of $J$ satisfies (3);

(iv) the representation (4) is satisfied in the following sense:

$$dJ(\Omega; \mathbf{T}) = - \int_{\partial S} \left[ (S'(D\mathbf{v})^{\top}D\mathbf{w} - s\mathbf{I}) : \frac{\partial \mathbf{v}}{\partial \mathbf{n}} \otimes \mathbf{n} + \mathbf{f} \cdot \mathbf{d} \right] \mathbf{T} \cdot \mathbf{n}. \quad (12)$$

In the remaining part we show the main steps of the proof of Theorem 3. Details can be found in [10], where the time-dependent problem is treated.
## 4 Formulation in the fixed domain

In this section we transform the problem $(P(\Omega_\varepsilon))$ to the fixed domain $\Omega$. Let us introduce the following notation:

$$ \mathbf{v}_{\varepsilon}(\mathbf{x}) := \mathbb{N}^{\top}(\mathbf{x})\bar{v}_{\varepsilon}(\mathbf{y}(\mathbf{x})), \mathbf{x} \in \Omega. $$

Note that the definition of $\mathbf{v}_\varepsilon$ implies that $\mathrm{div} \mathbf{v}_\varepsilon = 0$. The new function $\mathbf{v}_\varepsilon \in \mathbf{g} + \mathbf{W}_{0,\mathrm{div}}^{1,r}(\Omega)$ satisfies the equality

$$ \int_{\Omega} [\mathbf{g}\mathcal{S}(\mathcal{D}_{\varepsilon}\mathbf{v}_{\varepsilon}) : \mathcal{D}_{\varepsilon}\phi - \mathbf{v}_{\varepsilon} \otimes \mathbf{v}_{\varepsilon} : \nabla\phi + \mathcal{C}\mathbf{v}_{\varepsilon} \cdot \phi] \\ = \int_{\Omega} \mathbf{f} \cdot \phi + \langle \mathbf{A}_{\varepsilon}^{1}, \phi \rangle_{\mathbf{W}_{0,\mathrm{div}}^{1,2}(\Omega)} \quad \text{for all } \phi \in \mathbf{W}_{0,\mathrm{div}}^{1,r}(\Omega), \quad (13) $$

where the term $\mathbf{A}_\varepsilon^1$ on the right hand side is defined for $\phi \in \mathbf{W}_{0,\mathrm{div}}^{1,2}(\Omega)$ by

$$ \langle \mathbf{A}_{\varepsilon}^{1}, \phi \rangle_{\mathbf{W}_{0,\mathrm{div}}^{1,2}(\Omega)} = \int_{\Omega} \left[ \mathbf{v}_{\varepsilon} \otimes \mathbb{N}^{-\top} \mathbf{v}_{\varepsilon} : \nabla(\mathbb{N}^{-\top} \phi) - \mathbf{v}_{\varepsilon} \otimes \mathbf{v}_{\varepsilon} : \nabla \phi \right. \\ \left. + (\mathcal{C} - \mathbf{g}\mathbb{N}^{-1}\mathcal{C}\mathbb{N}^{-\top})\mathbf{v}_{\varepsilon} \cdot \phi + (\mathbf{g}\mathbb{N}^{-1}\mathbf{f} \circ \mathbf{y} - \mathbf{f}) \cdot \phi \right]. \quad (14) $$

Here $\mathcal{D}_{\varepsilon}\mathbf{v}_{\varepsilon} := \mathbf{g}^{-1}(\mathbb{N}\nabla(\mathbb{N}^{-\top}\mathbf{v}_{\varepsilon}))_{\mathrm{sym}}$.

Applying change of coordinates we further get:

$$ J(\Omega_\varepsilon) = \int_{\Omega} \left[ \mathbf{g} (\mathbb{N}^{-1}\mathcal{C}\mathbb{N}^{-\top}\mathbf{v}_\varepsilon - \mathbb{N}^{-1}\mathbf{f} \circ \mathbf{y}) \cdot \boldsymbol{\xi} \right. \\ \left. + (\mathbb{N}^{\top}\mathcal{S}(\mathcal{D}_\varepsilon\mathbf{v}_\varepsilon) - \mathbf{v}_\varepsilon \otimes (\mathbb{N}^{-\top}\mathbf{v}_\varepsilon)) : \nabla(\mathbb{N}^{-\top}\boldsymbol{\xi}) \right]. \quad (15) $$

Now after all quantities and equations have been transformed to the fixed domain $\Omega$, we can analyze the limit $\varepsilon \to 0$.

**Lemma 1.** The sequence $\{\mathbf{v}_\varepsilon\}_{\varepsilon>0}$ is bounded in $\mathbf{W}_{0,\mathrm{div}}^{1,r}(\Omega)$ and satisfies:

$$ \begin{array}{ll} \mathbf{v}_{\varepsilon} \to \mathbf{v} & \text{weakly in } \mathbf{W}_{0,\mathrm{div}}^{1,r}(\Omega), \\ \mathbb{N}^{\top}\mathcal{S}(\mathcal{D}_{\varepsilon}\mathbf{v}_{\varepsilon}) \to \mathcal{S}(\mathcal{D}\mathbf{v}) & \text{weakly in } L^{r'}(\Omega, \mathbb{R}^{2 \times 2}), \\ \mathbf{A}_{\varepsilon}^{1} \to \mathbf{0} & \text{weakly in } \mathbf{W}_{0,\mathrm{div}}^{1,r}(\Omega)^*. \end{array} $$

In particular, $\mathbf{v}$ is the unique weak solution to $(P(\Omega))$.

## 5 Existence of material derivative

Our next task is to identify $\tilde{\mathbf{v}}$ as the limit of the sequence $\{\mathbf{u}_\varepsilon\}$, where

$$ \mathbf{u}_{\varepsilon} := \frac{\mathbf{v}_{\varepsilon} - \mathbf{v}}{\varepsilon}. $$
First we write down the system for the differences $\mathbf{u}_\varepsilon$. Subtracting (13) and (6) we find that $\mathbf{u}_\varepsilon \in \mathbf{W}_{0,\mathrm{div}}^{1,r}(\Omega)$ satisfies the equality

$$ \int_{\Omega} \left[ \frac{1}{\varepsilon} \mathbf{g}(\mathcal{S}(\mathcal{D}_{\varepsilon} \mathbf{v}_{\varepsilon}) - \mathcal{S}(\mathcal{D}_{\varepsilon} \mathbf{v})) : \mathcal{D}_{\varepsilon} \phi + \mathcal{C} \mathbf{u}_{\varepsilon} \cdot \phi - (\mathbf{v}_{\varepsilon} \otimes \mathbf{u}_{\varepsilon} + \mathbf{u}_{\varepsilon} \otimes \mathbf{v}) : \nabla \phi \right] = \frac{1}{\varepsilon} \langle \mathbf{A}_{\varepsilon}, \phi \rangle_{\mathbf{W}_{0,\mathrm{div}}^{1,2}(\Omega)} \quad (16) $$

for all $\phi \in \mathbf{W}_{0,\mathrm{div}}^{1,r}(\Omega)$. The term $\mathbf{A}_\varepsilon \in \mathbf{W}_{0,\mathrm{div}}^{1,2}(\Omega)^*$ on the right hand side is defined as follows:

$$ \mathbf{A}_{\varepsilon} := \mathbf{A}_{\varepsilon}^{1} + \mathbf{A}_{\varepsilon}^{2}, $$

$\mathbf{A}_{\varepsilon}^{1}$ is given by (14),

$$ \langle \mathbf{A}_{\varepsilon}^{2}, \phi \rangle_{\mathbf{W}_{0,\mathrm{div}}^{1,2}(\Omega)} := \int_{\Omega} [\mathbf{N}^{\top} \mathcal{S}(\mathcal{D}_{\varepsilon} \mathbf{v}) : \nabla(\mathbf{N}^{-\top} \phi) - \mathcal{S}(\mathcal{D} \mathbf{v}) : \mathcal{D} \phi]. $$

Next we state the properties of the sequence $\{\mathbf{u}_\varepsilon\}_{\varepsilon>0}$.

**Lemma 2.** The sequence $\{\mathbf{u}_\varepsilon\}_{\varepsilon>0}$ is bounded in $\mathbf{W}_{0,\mathrm{div}}^{1,2}(\Omega)$. Further it holds:

$$ \frac{\mathbf{A}_{\varepsilon}}{\varepsilon} \to \mathbf{A}_{0}^{\prime} \qquad \text{weakly in } \mathbf{W}_{0,\mathrm{div}}^{1,2}(\Omega)^{*}, $$

$$ \mathbf{u}_{\varepsilon} \to \tilde{\mathbf{v}} \qquad \text{weakly in } \mathbf{W}_{0,\mathrm{div}}^{1,2}(\Omega), $$

$$ \frac{1}{\varepsilon} (\mathbf{g}(\mathcal{S}(\mathcal{D}_{\varepsilon} \mathbf{v}_{\varepsilon}) - \mathcal{S}(\mathcal{D}_{\varepsilon} \mathbf{v})), \mathcal{D}_{\varepsilon} \phi) \to (\mathcal{S}'(\mathcal{D} \mathbf{v}) \mathcal{D} \tilde{\mathbf{v}}, \mathcal{D} \phi) \quad \text{for all } \phi \in \mathbf{W}^{1, \frac{2r}{4-r}}(\Omega), $$

where $\mathbf{A}'_0$ is defined in (10a) and $\tilde{\mathbf{v}}$ is the solution of $(P_{\mathrm{lin}}(\Omega))$ with $\mathbf{F} := \mathbf{A}'_0$ and $\mathbf{h} = \mathbf{0}$.

This completes the proof of Theorem 2.

# 6 Shape gradient of J

To prove Theorem 3, we decompose the fraction

$$ \frac{J(\Omega_{\varepsilon}) - J(\Omega)}{\varepsilon} = J_{1}^{\varepsilon} + J_{2}^{\varepsilon} $$

in a suitable way. Using Lemma 1 and Lemma 2 and the properties of $\mathbf{g}$ and $\mathbf{N}'$, it is then possible to show that

$$ J_{1}^{\varepsilon} \to J_{\mathbf{v}}(\tilde{\mathbf{v}}) \quad \text{and} \quad J_{2}^{\varepsilon} \to J_{e}(\mathbf{T}). $$

The continuity of the map $\mathbf{T} \mapsto dJ(\Omega; \mathbf{T})$ follows from the estimate (11).

# Acknowledgement

The work of J. Stebel was supported by the Czech Science Foundation (GAČR) grant No. 201/09/0917 and by the ESF grant Optimization with PDE Constraints. The research of J. Sokołowski is performed in the framework of the ANR project GAOS.
Bibliography

[1] F. Abraham, M. Behr, and M. Heinkenschloss. Shape optimization in steady blood flow: A numerical study of non-newtonian effects. *Computer Methods in Biomechanics and Biomedical Engineering*, 8(2):127-137, 2005.

[2] L. Consiglieri, Š. Nečasová, and J. Sokolowski. New approach to the incompressible Maxwell-Boussinesq approximation: existence, uniqueness and shape sensitivity. *J. Differential Equations*, 249(12):3052-3080, 2010. ISSN 0022-0396.

[3] J. Haslinger and J. Stebel. Shape optimization for navier-stokes equations with algebraic turbulence model: Numerical analysis and computation. *Applied Mathematics and Optimization*, 63(2):277-308, 2011.

[4] P. Kaplický, J. Málek, and J. Stará. On global existence of smooth two-dimensional steady flows for a class of non-newtonian fluids under various boundary conditions. *Applied Nonlinear Analysis*, pages 213-229, 2002.

[5] J. Málek and K. R. Rajagopal. Mathematical issues concerning the Navier-Stokes equations and some of its generalizations. In *Evolutionary equations. Vol. II*, Handb. Differ. Equ., pages 371-459. Elsevier/North-Holland, Amsterdam, 2005.

[6] M. Moubachir and J.-P. Zolésio. *Moving shape analysis and control*, volume 277 of *Pure and Applied Mathematics (Boca Raton)*. Chapman & Hall/CRC, Boca Raton, FL, 2006. ISBN 978-1-58488-611-2; 1-58488-611-0. Applications to fluid structure interactions.

[7] P. I. Plotnikov and J. Sokolowski. Shape derivative of drag functional. *SIAM J. Control Optim.*, 48(7):4680-4706, 2010. ISSN 0363-0129.

[8] P. I. Plotnikov and J. Sokolowski. *Compressible Navier-Stokes equations. Theory and shape optimization*, volume 73 of *Mathematical Monographs*. Birkhäuser Basel, 2012. ISBN 978-3-0348-0366-3.

[9] T. Slawig. Distributed control for a class of non-newtonian fluids. *Journal of Differential Equations*, 219(1):116-143, 2005.

[10] J. Sokolowski and J. Stebel. Shape sensitivity analysis of time-dependent flows of incompressible non-Newtonian fluids. *Control and Cybernetics*, 40 (4):1077-1097, 2011.

[11] J. Sokolowski and J.-P. Zolésio. *Introduction to shape optimization. Shape sensitivity analysis*, volume 16 of Springer Series in Computational Mathematics. Springer-Verlag, Berlin, 1992.

[12] F. Tröltzsch. *Optimal control of partial differential equations*, volume 112 of *Graduate Studies in Mathematics*. American Mathematical Society, Providence, RI, 2010. ISBN 978-0-8218-4904-0. Theory, methods and applications, Translated from the 2005 German original by Jürgen Sprekels.

[13] D. Wachsmuth and T. Roubíček. Optimal control of planar flow of incompressible non-Newtonian fluids. *Z. Anal. Anwend.*, 29(3):351-376, 2010.