1 Introduction

In the field of optimization constrained by partial differential equations (PDEs), there is a large class of problems where not only optimal controls are to be found but also an optimal shape of the experimental domain. Here, the contour of the domain \(\Omega \), where the PDE models the effects of interest, plays the role of the optimization variable. Possible variants are that the outer shape of \(\Omega \) is to be determined, e.g., when \(\Omega \) represents a solid body, or interior interfaces, which separate spatially discontinuous coefficients such as material properties. Shape optimization in general is nowadays an active field of research with applications ranging from magnetostatics [8], interface identification in transmission processes [12, 22, 27], fluid dynamics [2, 9, 25], acoustics [31], image restoration and segmentation [14] and composite material identification [23, 28] to nano-optics [15].

In this article, we focus on shape optimization in fluid dynamics, which is also one of the pioneering applications in this field [11, 17, 20]. In general, the optimization problem can be formulated as

$$\begin{aligned}&\min \limits _{\Omega \in {G}_\text {adm}} \qquad j(y, \Omega ) \nonumber \\&\quad \text {s.t.} \qquad E(y, \Omega ) = 0 \end{aligned}$$
(1)

where j is a shape functional depending on a state variable y and the shape of the domain \(\Omega \). Moreover, y fulfills the PDE constraint E, which itself depends on \(\Omega \). A typical example is an obstacle specimen \({\Omega _{\text {obs}}}\) in a flow tunnel \(\Omega \) as depicted in Fig. 1. One of the main questions is an appropriate choice of the set of admissible shapes \({G}_\text {adm}\), in which optimization takes place. For problems of this type, two prominent approaches can be identified in the literature. On the one hand, the Hadamard–Zolésio structure theorem is applied, which allows to trace back changes in the objective j solely to variations of the boundary \({\Gamma _{\text {obs}}}\) (see for instance [3, 30]). It is thus possible to define directional shape derivatives via variations of \({\Gamma _{\text {obs}}}\) in a direction normal to the boundary. Together with the choice of an appropriate shape and tangent space, this allows to represent the sensitivity for j w.r.t. \({\Gamma _{\text {obs}}}\) as a gradient. This is then interpreted as a deformation to \({\Gamma _{\text {obs}}}\) and a new discretization mesh for the resulting domain can be computed. By this step, the mesh quality of the deformed domain can be ensured as pursued in, e.g., [6, 32]. Alternatively, the definition of shape and tangent space includes the surrounding domain \(\Omega \), which immediately results in deformation information for the entire mesh (e.g., [7, 8, 26]) and makes the additional call to a mesh generator superfluous. Typical approaches consider interpreting the shape sensitivity as a force term in linear elastic models described over \(\Omega \). The resulting displacement field is then applied as a mesh deformation. Especially in recent works (see for instance [4, 5, 13]), linear elastic extension equations are considered and, in particular, a very small or even zero first Lamé constant is favored.

Moreover, a descent method allows to control the mesh quality from one iteration to the next, i.e., for one deformation. Yet, in the limit of the sequence of design updates, quality is typically lost. This effect is described in, e.g.,[23], where variable interfaces must be prevented from overlapping.

In this article, we follow a different approach, which gives a higher level of control on the quality of the mesh around the optimal shape. Based on the method of mappings (cf. to [21]), the question for admissible shapes \({G}_\text {adm} :=\lbrace F(\Omega ):F \in \mathcal {F}_\text {adm}\rbrace \) in (1) is translated to the choice of appropriate function spaces, in which a deformation from reference to the optimal configuration is to be found. Here, \(\mathcal {F}_\text {adm}\) denotes a set of admissible mappings. Starting from a reference configuration \(\Omega \), it is then optimized over the transformations \(F(\Omega )\) yet without explicitly performing mesh deformations. For this purpose, the PDE constraint is transformed to the virtual domain as \(E(y, F(\Omega ))\). The optimization problem then turns into a classical optimal control in the form of

$$\begin{aligned}&\min \limits _{F\in \mathcal {F}_\text {adm}} \quad \qquad j(y, F(\Omega )) \nonumber \\&\text {s.t.}\quad \quad E(y, F(\Omega )) = 0. \end{aligned}$$
(2)

This approach is a recent field of studies and applied in, e.g., [2, 19, 29]. Also based on this approach is the investigation in [13], which is the starting point for the consideration in the present article. Here, the problem in (2) is formulated as

$$\begin{aligned} \min \limits _{{c}\in L^2({\Gamma _{\text {obs}}})}&j(y, F(\Omega )) + \frac{\alpha }{2} \Vert c \Vert ^2_{L^2({\Gamma _{\text {obs}}})} \nonumber \\ \text {s.t.}&E(y, F(\Omega )) = 0 \nonumber \\&F= {{\,\mathrm{id}\,}}+ w \quad \text { in } \Omega \nonumber \\&\det (DF) \ge {\eta _{\text {det}}}\quad \text { in } \Omega \nonumber \\&w = S({c}, \Omega ) \end{aligned}$$
(3)

in terms of a regularization parameter \(\alpha > 0\) and a bound \({\eta _{\text {det}}}> 0\) on the determinant of the derivative of the mapping function \(F\). The focus of the investigations therein is on the extension operator S. It is suggested to choose S to be the composition of mappings \(c \mapsto b \mapsto w\). Here, \(c \mapsto b\) is realized via the solution operator of a Laplace–Beltrami equation on \({\Gamma _{\text {obs}}}\). The mapping to the actual displacement, i.e., \(b \mapsto w\), is then chosen to be the solution operator of a vector-valued elliptic equation, such as a linear elastic model. It is proven that—under certain circumstances—the domain mapping \(F\) is locally a \(C^1(\bar{\Omega }, {\mathbb {R}}^d)\)-diffeomorphism provided that \(\det (DF) \ge {\eta _{\text {det}}}\) is fulfilled.

Alternatively, one could tackle problem (1) with an iterative descent algorithm. This would mean to find a suitable deformation for \(\Omega \) leading to a descent in the objective j. Assume a resulting sequence of domains

$$\begin{aligned} \Omega =:\Omega ^1 \overset{F_1}{\rightarrow } \Omega ^2 \overset{F_2}{\rightarrow } \Omega ^3 \cdots \overset{F_k}{\rightarrow } \Omega ^{k+1} \end{aligned}$$
(4)

Together with the extension operator S and the condition \(\det {DF_i} \ge {\eta _{\text {det}}}\), for all i, it would be possible to guarantee the quality of deformations in each step. Yet, the composition of deformations \(\tilde{F}_k :=F_k \circ \cdots \circ F_1\) would potentially violate the condition \(\det {D\tilde{F}_k} \ge {\eta _{\text {det}}}\) and thus lead to bad mesh quality in \(\Omega ^{k+1}\).

The main focus of our present article is a numerical study of different choices of the extension operator S. It turns out that optimization settings with expected larger deformations are a limiting factor for linear operators S. This limitation is due to the fact that the structure of a shape space, that is as large as possible, can hardly be linear since this would require to explain what scalar multiples or sums of shapes are. Yet, with the method of mappings and a linear extension operator S we approximate the set of admissible shapes locally by a linear function space of admissible deformations to a reference configuration.

We thus suggest a nonlinear extension mapping and present numerical studies on the applicability. It should be mentioned that the theory developed so far is not applicable in this case. It only applies to the linear choice of S, which is a special case of the more general consideration in this article.

The motivation for the choice of S in this present work is the observation that, on the one hand, via the condition \(\det (DF) \ge {\eta _{\text {det}}}\) the local injectivity of \(F\) can be ensured. But on the other hand, this limits significantly the subset of admissible transformations \(\mathcal {F}_\text {adm}\) and thus affects optimal shapes as outlined in the last section of this article. It is thus the task to find an operator S which prevents \(\det (DF) \ge {\eta _{\text {det}}}\) from becoming active even for large deformations. We also discuss cases where the reference domain is not of circular shape and illustrate the performance of the extension and influence on the mesh quality in a deformed domain. The intention of this experiment is to demonstrate that the set of shapes \({G}_\text {adm}\), which is constructed via the mappings from \(\mathcal {F}_\text {adm}\), can be extended significantly and its dependence on the choice of a reference domain \(\Omega \) is reduced. In particular, the studies illuminate whether large deformations in the optimization are possible for general reference configurations, which do not fulfill certain properties like convexity or an injective normal vector field.

This article is structured as follows: In Sect. 2 the shape optimization problem is set up and formulated in terms of the method of mappings. Section 3 is devoted to the nonlinear extension model and, furthermore, the derivation of necessary optimality conditions and the presentation of an optimization algorithm. In Sect. 4, numerical studies are conducted and discussed. The article closes in Sect. 5 with a conclusion of the results.

Fig. 1
figure 1

Sketch of the holdall domain \({G}= \Omega \cup {\Omega _{\text {obs}}}\)

Full size image

2 Optimization Problem

We carry out our considerations based on a classical optimization problem in the field of fluid dynamics described in [20]. In a d-dimensional, bounded domain \(\Omega \) with Lipschitz boundary, as sketched in Fig. 1, we consider the minimization of the following energy dissipation functional

$$\begin{aligned} \min \limits _{\Gamma _{\text {obs}}}\; j(v, {\Gamma _{\text {obs}}}) = \frac{\nu }{2} \int _\Omega \sum _{i,j=1}^{d} \left( \frac{\partial v_i}{\partial x_j} \right) ^2 \, \mathrm{d}x \end{aligned}$$
(5)

where the contour \({\Gamma _{\text {obs}}}\) of the obstacle \({\Omega _{\text {obs}}}\) is assumed to be variable. Here \({\Omega _{\text {obs}}}\) is non-empty, connected set and \({\Gamma _{\text {obs}}}\) is a smooth and compact Riemannian manifold without a boundary. The spatial dimension is chosen as \(d\in \lbrace 2,3 \rbrace \). In (5), the velocity field denoted by v is given in terms of the stationary, incompressible Navier–Stokes equations

$$\begin{aligned} -\nu \Delta v + (v\cdot \nabla )v + \nabla p= & {} 0 \quad \text { in } \quad \Omega \nonumber \\ {{\,\mathrm{div}\,}}v= & {} 0 \quad \text { in } \quad \Omega \nonumber \\ v= & {} v_\infty \quad \text { on } \quad {\Gamma _{\text {in}}}\nonumber \\ v= & {} 0 \quad \text { on } \quad {\Gamma _{\text {obs}}}\cup {\Gamma _{\text {wall}}}\nonumber \\ pn - \nu \frac{\partial v}{\partial n} \quad= & {} 0 \quad \text { on } \quad {\Gamma _{\text {out}}}. \end{aligned}$$
(6)

Together with \({\Gamma _{\text {obs}}}\) the fluid domain \(\Omega \) is allowed to change, but the outer boundaries, i.e., \({\Gamma _{\text {in}}}\), \({\Gamma _{\text {out}}}\) and \({\Gamma _{\text {wall}}}\), of the experiment are fixed. In (6), p denotes the pressure, \(v_\infty \) describes the velocity profile at the inflow boundary, n is the outer normal vector and \(\nu \) the viscosity. Here, \(v_\infty \in H^1(\Omega )\) is assumed to satisfy the compatibility condition

$$\begin{aligned} \int \limits _{{\Gamma _{\text {in}}}} v_{\infty } \cdot n = 0. \end{aligned}$$
(7)

Furthermore, we assume that \({\Gamma _{\text {in}}}\) and \({\Gamma _{\text {obs}}}\) have positive Lebesgue measure, \({\Gamma _{\text {obs}}}\cap ({\Gamma _{\text {in}}}\cup {\Gamma _{\text {wall}}}\cup {\Gamma _{\text {out}}})= \emptyset \) holds during the entire optimization.

For the shape optimization of a specimen \({\Omega _{\text {obs}}}\) with respect to functionals of type (5), it is essential to exclude trivial solutions. Here, shrinking \({\Omega _{\text {obs}}}\) to a point or translations toward \({\Gamma _{\text {wall}}}\) represents undesired descent directions. Thus, the optimization problem has to be additionally constrained to geometrical conditions. Our benchmark problem is to find optimal shapes of a specimen with a given volume, which remains located in the center of the flow tunnel. This is achieved by fixing barycenter and volume of the obstacle \({\Omega _{\text {obs}}}\) with the constraints

$$\begin{aligned} {\text {vol}}({\Omega _{\text {obs}}})&= \int _{\Omega _{\text {obs}}}1\, \mathrm{d}x =\; \text {const}, \end{aligned}$$
(8)
$$\begin{aligned} {\text {bc}}({\Omega _{\text {obs}}})&= \frac{1}{{\text {vol}}({\Omega _{\text {obs}}})}\int _{\Omega _{\text {obs}}}x \, \mathrm{d}x =\; \text {const}. \end{aligned}$$
(9)

Since the computation for the barycenter involves the volume of \({\Omega _{\text {obs}}}\) itself, these conditions are coupled in principle. Yet, if (8) is fulfilled, the term \({\text {vol}}({\Omega _{\text {obs}}})^{-1}\) in (9) is constant and can thus be factored out. By further assuming that the barycenter of the specimen \({\Omega _{\text {obs}}}\) is \(0\in {\mathbb {R}}^d\), it is thus sufficient to require \(\int _{\Omega _{\text {obs}}}x \, \mathrm{d}x =\; 0\).

In the following, for a vector-valued function \(f: {\mathbb {R}}^d \rightarrow {\mathbb {R}}^d\), we denote by Df the Jacobian matrix with the ordering \(Df = \left( \frac{\partial f_i}{\partial x_j}\right) _{i,j=1,\dots , d} \in {\mathbb {R}}^{d\times d}\). Let further

$$\begin{aligned} V:= & {} \lbrace v \in H^1(\Omega , {\mathbb {R}}^d) :{{\,\mathrm{div}\,}}(v)=0, v\vert _{\Gamma _{\text {in}}}=v_\infty , v\vert _{{\Gamma _{\text {wall}}}\cup {\Gamma _{\text {obs}}}}=0\; \text {a.e.} \rbrace ,\nonumber \\ V_0:= & {} \lbrace v \in H^1(\Omega , {\mathbb {R}}^d) :{{\,\mathrm{div}\,}}(v)=0, v\vert _{{\Gamma _{\text {in}}}\cup {\Gamma _{\text {wall}}}\cup {\Gamma _{\text {obs}}}}=0\; \text {a.e.} \rbrace ,\nonumber \\ Q:= & {} \lbrace p \in L^2(\Omega ):\int _\Omega p\, \mathrm{d}x = 0\rbrace \end{aligned}$$
(10)

and consider the weak formulation of the PDE constraint (6):

Find \((v,p) \in V\times Q\) such that

$$\begin{aligned}&\nu \int _{\Omega } D v : D {\delta _{v}} + (D v\, v)\cdot {\delta _{v}} - p {{\,\mathrm{Tr}\,}}(D {\delta _{v}}) \,\mathrm{d}x = 0,\nonumber \\&\qquad -\int _\Omega {\delta _{p}} {{\,\mathrm{Tr}\,}}\left( D v \right) \, \mathrm{d}x = 0 \end{aligned}$$
(11)

for all test functions \(({\delta _{v}}, {\delta _{p}}) \in V_0 \times Q\). Note that within this article we are using the symbol \({\delta _{\cdot }}\) for test functions associated with a given variable.

In order to reformulate the optimization problem (5)–(9) as an optimal control problem in appropriate function spaces, we consider from now on the domain \(\Omega \) as a fixed reference configuration. Let \(F = {{\,\mathrm{id}\,}}+ w\) with \(w \in W^{1, \infty }(\Omega , {\mathbb {R}}^d)\) such that F results in an admissible deformation for \(\Omega \). For the method of mappings, we then consider the state (11), objective (5) and the corresponding state variable v in terms of \(F(\Omega )\).

By means of standard computations, we obtain the weak formulation of the optimization problem pulled back to the reference domain \(\Omega \) by

$$\begin{aligned}&\min \limits _{F \in \mathcal {F}_\text {adm}}\quad j(v,F(\Omega )) = \frac{\nu }{2} \int _{\Omega } \left( D v (DF)^{-1} \right) : \left( D v (DF)^{-1} \right) \det (DF)\, \mathrm{d}x\\&\quad \text {s.t. }\quad \int _{\Omega } \bigl [ \nu \left( D v (DF)^{-1} \right) : \left( D {\delta _{v}} (DF)^{-1} \right) + (D v (DF)^{-1}v)\cdot {\delta _{v}} \nonumber \end{aligned}$$
(12)
$$\begin{aligned}&\quad \qquad - p {{\,\mathrm{Tr}\,}}\left( D {\delta _{v}} (DF)^{-1} \right) \bigr ] \det (DF)\, \mathrm{d}x = 0, \end{aligned}$$
(13)
$$\begin{aligned}&\quad \qquad - \int _\Omega {\delta _{p}} {{\,\mathrm{Tr}\,}}(D v(DF)^{-1}) \det (DF) \,\mathrm{d}x = 0, \end{aligned}$$
(14)
$$\begin{aligned}&\quad \qquad \int _{{\Omega _{\text {obs}}}} \det (DF) - 1 \, \mathrm{d}x = 0, \end{aligned}$$
(15)
$$\begin{aligned}&\quad \qquad \int _{{\Omega _{\text {obs}}}} F \, \det (DF) \, \mathrm{d}x = 0 \end{aligned}$$
(16)

for all test functions \(({\delta _{v}}, {\delta _{p}}) \in V\times Q\). The optimization problem (12)–(16) still leaves open the question for the set of admissible mappings \(\mathcal {F}_\text {adm}\). We thus follow the same approach as in [13] and translate it into the form of (3). By reformulating the constraint \(\det (DF) \ge {\eta _{\text {det}}}\) as a penalty term, we obtain the final optimal control problem

$$\begin{aligned}&\min \limits _{{c}\in L^2({\Gamma _{\text {obs}}})} \;J(v,{c}) :=j(v,F) + \frac{\alpha }{2} \int _{\Gamma _{\text {obs}}}c^2\, \mathrm{d}s + \frac{\beta }{2} \int _\Omega \left( ({\eta _{\text {det}}}- \det (DF))_+\right) ^2\, \mathrm{d}x\nonumber \\&\text {s.t. }\quad {(13)\,\text {to}\, (16)}\nonumber \\&F= {{\,\mathrm{id}\,}}+ w\nonumber \\&w = S(c), \end{aligned}$$
(17)

where \((\cdot )_+\) denotes the positive-part function. The missing piece is now mapping from a scalar-valued boundary control \({c}\) to admissible deformation fields w, which is the subject of the next section.

3 Nonlinear Extension Operators

Consider the optimal control problem (17). The core of the reformulated shape optimization is the choice of the extension operator S, which links a scalar-valued boundary control \({c}\) living on \({\Gamma _{\text {obs}}}\) to a vector-valued displacement field w in \(\Omega \). A domain transformation mapping \(F= {{\,\mathrm{id}\,}}+ w\) is then obtained by the so-called perturbation of identity. In particular, w has to fulfill certain regularity properties as investigated in [13]. It yet turns out in Sect. 4 that for large deformations, i.e., when the reference domain and the optimal configuration differ significantly, linear operators S do not lead to satisfying results. Note that the choice of S significantly influences the set of reachable shapes \({G}_\text {adm}\) determined via \(\mathcal {F}_\text {adm}\). It is thus our intention to find S which allows for large deformations without significantly restricting \({G}_\text {adm}\). Simultaneously, the corresponding mesh deformations \(F(\Omega )\) should exhibit high element qualities for further usage in numerical simulations.

The focus of the present article is thus to propose and study nonlinear extensions S given in terms of the solution operator of the coupled PDEs

$$\begin{aligned} b - \Delta _{\Gamma _{\text {obs}}}b= & {} cn \quad \text { on } \quad {\Gamma _{\text {obs}}}\nonumber \\ -{{\,\mathrm{div}\,}}(\nabla w + \nabla w^\top ) + {\eta _{\text {ext}}}(w \cdot \nabla )w= & {} 0 \quad \text { in } \quad \Omega \nonumber \\ (\nabla w + \nabla w ^\top ) \cdot n= & {} b\quad \text { on } \quad {\Gamma _{\text {obs}}}\nonumber \\ w= & {} 0 \quad \text { on } \quad {\Gamma _{\text {wall}}}\cup {\Gamma _{\text {in}}}\cup {\Gamma _{\text {out}}}\end{aligned}$$
(18)

In the equation above, \(\Delta _{\Gamma _{\text {obs}}}\) denotes the vector-valued Laplace–Beltrami operator. Note that by solving (18) the scalar-valued control \({c}\in L^2({\Gamma _{\text {obs}}})\) is mapped to a vector-valued quantity \(b\in H^2({\Gamma _{\text {obs}}})\). The benefit of this particular extension operator, and especially the nonlinearity \({\eta _{\text {ext}}}(w \cdot \nabla )w\), which is in the focus of this article, becomes particularly visible for experiments with large deformations as pointed out in Sect. 4.2. A popular choice, as discussed in the introduction, is to define the extension only via the linear term \({{\,\mathrm{div}\,}}(\nabla w + \nabla w^\top )\). Yet, this restricts the set \(\mathcal {F}_\text {adm}\) significantly. This is visible especially for problems in fluid dynamics, as pointed out in Sect. 4, where the reference shape is of spherical type, but the optimum to be found is stretched and approximates non-smooth tip and back.

Problems arise due to strong compressions of finite elements in the discretization orthogonal to the main deformation direction. This observation motivates to add the nonlinear advection term \({\eta _{\text {ext}}}(w \cdot \nabla )w\), which—geometrically speaking—promotes displacements w where nodes move along large gradients. This results in a homogeneous distribution of finite elements even around approximately non-smooth regions of \({\Gamma _{\text {obs}}}\).

Firstly, we derive the weak form of Laplace–Beltrami equation, i.e., the first equation of (18):

Find \(b \in H^2({\Gamma _{\text {obs}}}):\)

$$\begin{aligned} \int _{\Gamma _{\text {obs}}}b\cdot {\delta _{b}} + D_{{\Gamma _{\text {obs}}}} b : D_{{\Gamma _{\text {obs}}}} {\delta _{b}} \, \mathrm{d}s = \int _{\Gamma _{\text {obs}}}{c}n \cdot {\delta _{b}} \, \mathrm{d}s \end{aligned}$$
(19)

for all \({\delta _{b}} \in H^2({\Gamma _{\text {obs}}})\).

For the second equation in (18), which specifies the mapping from vector-valued function b to the domain deformation, we follow the argumentation in [13] and obtain the weak formulation in the space

$$\begin{aligned} W :=\left\{ w \in H^{\frac{7}{2}}(\Omega , {\mathbb {R}}^d): w\vert _{{\Gamma _{\text {in}}}\cup {\Gamma _{\text {wall}}}\cup {\Gamma _{\text {out}}}} = 0\, \text {a.e.} \right\} \end{aligned}$$

as follows:

Find \(w \in W\) such that

$$\begin{aligned} \int _{\Omega } (Dw + Dw^\top ):D{\delta _{w}} + {\eta _{\text {ext}}}(Dw\, w)\cdot {\delta _{w}} \, \mathrm{d}x = \int _{\Gamma _{\text {obs}}}b {\delta _{w}} \, \mathrm{d}s \end{aligned}$$
(20)

for all \({\delta _{w}} \in W\) and in terms of \({\eta _{\text {ext}}}\ge 0\).

Here, \(D_{\Gamma _{\text {obs}}}\) denotes the derivative tangential to \({\Gamma _{\text {obs}}}\). In (19) the scalar-valued boundary control \({c}\) is multiplied with the outer normal vector field n to \(\Omega \) at \({\Gamma _{\text {obs}}}\). Then, a vector-valued Laplace–Beltrami equation is solved over \({\Gamma _{\text {obs}}}\). This is coupled with nonlinear (20) where the influence of the advection term is controlled via \({\eta _{\text {ext}}}\). Note that the linear extension operators investigated in [13] arise as a special case of system (19) and (20).

Now that the extension operator is chosen, we can combine the optimization problem (17) with the extension operator (19) and (20) to obtain the Lagrangian

$$\begin{aligned}&{\mathcal {L}}( w, v, p, b, {c}, {\psi _{w}},{\psi _{v}},{\psi _{p}}, {\psi _{b}}, {\psi _{{\text {vol}}}}, {\psi _{{\text {bc}}}}) =\frac{\nu }{2} \int _{\Omega } \left( D v (DF)^{-1} \right) : \nonumber \\&\quad \left( D v (DF)^{-1} \right) \det (DF)\, \mathrm{d}x +\frac{\alpha }{2}\int _{\Gamma _{\text {obs}}}{c}^2 \, \mathrm{d}s\nonumber \\&\qquad + \frac{\beta }{2} \int _{\Omega } (( {\eta _{\text {det}}}- \det (DF) )_+ )^2 \, \mathrm{d}x - \int _{\Omega } \bigl [ \nu \left( D v (DF)^{-1} \right) : \nonumber \\&\quad \left( D {\psi _{v}} (DF)^{-1} \right) + (D v (DF)^{-1}v)\cdot {\psi _{v}} - p {{\,\mathrm{Tr}\,}}\left( D {\psi _{v}} (DF)^{-1} \right) \bigr ] \det (DF)\, \mathrm{d}x\nonumber \\&\qquad + \int _\Omega {\psi _{p}} {{\,\mathrm{Tr}\,}}(D v(DF)^{-1}) \det (DF) \,\mathrm{d}x\nonumber \\&\qquad - \int _\Omega (D w + D w^\top ): D {\psi _{w}} + {\eta _{\text {ext}}}(Dw\, w)\cdot {\psi _{w}} \, \mathrm{d}x + \int _{\Gamma _{\text {obs}}}b\cdot {\psi _{w}} \, \mathrm{d}s\nonumber \\&\qquad - \int _{\Gamma _{\text {obs}}}b\cdot {\psi _{b}} + D_{{\Gamma _{\text {obs}}}} b : D_{{\Gamma _{\text {obs}}}} {\psi _{b}} \, \mathrm{d}s + \int _{\Gamma _{\text {obs}}}{c}n \cdot {\psi _{b}} \, \mathrm{d}s \nonumber \\&\qquad - {\psi _{{\text {bc}}}} \cdot \int _{\Omega } (x + w) \, \det (DF) \, \mathrm{d}x - {\psi _{{\text {vol}}}} \int _{\Omega } \det (DF) - 1 \, \mathrm{d}x, \end{aligned}$$
(21)

where \({\psi _{\cdot }}\) denotes for each variable the associated multiplier. Note that there is no variable corresponding to the multipliers \({\psi _{{\text {vol}}}} \in {\mathbb {R}}\) and \({\psi _{{\text {bc}}}} \in {\mathbb {R}}^d\). These are the finite dimensional multipliers for the barycenter and volume condition (8) and (9).

Lemma 3.1

The first-order optimality system associated with the Lagrangian \({\mathcal {L}}\) in (21) is given by the derivatives \({\mathcal {L}}_{w}, {\mathcal {L}}_{v}, {\mathcal {L}}_{p}, {\mathcal {L}}_{b}, {\mathcal {L}}_{{\psi _{w}}}, {\mathcal {L}}_{{\psi _{v}}}, {\mathcal {L}}_{{\psi _{p}}}, {\mathcal {L}}_{{\psi _{b}}}, {\mathcal {L}}_{{c}}, {\mathcal {L}}_{{\psi _{{\text {vol}}}}}, {\mathcal {L}}_{{\psi _{{\text {bc}}}}}\) as

$$\begin{aligned} {\mathcal {L}}_{w} {\delta _{w}}= & {} - \nu \int _\Omega ( Dv (D F)^{-1}): (D v (D F)^{-1} D {\delta _{w}} (D F)^{-1}) \det (DF) \, \mathrm{d}x \nonumber \\&+ \frac{\nu }{2} \int _\Omega (Dv (D F)^{-1}): (Dv (D F)^{-1}) {{\,\mathrm{Tr}\,}}((DF)^{-1} D {\delta _{w}}) \det (DF) \, \mathrm{d}x \nonumber \\&- \beta \int _\Omega ({\eta _{\text {det}}}- \det (DF))_+ {{\,\mathrm{Tr}\,}}((DF)^{-1} D {\delta _{w}}) \det (DF) \, \mathrm{d}x\nonumber \\&+ \nu \int _\Omega (Dv (D F)^{-1} D {\delta _{w}} (D F)^{-1}) : ( D {\psi _{v}} (D F)^{-1}) \det (DF) \, \mathrm{d}x \nonumber \\\&+ \nu \int _\Omega (Dv (D F)^{-1}): ( D {\psi _{v}} (D F)^{-1} D {\delta _{w}} (D F)^{-1}) \det (DF) \, \mathrm{d}x \nonumber \\&- \nu \int _\Omega (Dv (D F)^{-1}) : ( D {\psi _{v}} (D F)^{-1}){{\,\mathrm{Tr}\,}}((DF)^{-1} D {\delta _{w}}) \det (DF) \, \mathrm{d}x \nonumber \\&+ \int _\Omega (D v (D F)^{-1} D {\delta _{w}} (D F)^{-1}\, v)\cdot {\psi _{v}} \det (DF) \, \mathrm{d}x \nonumber \\&- \int _\Omega (D v (D F)^{-1} \, v)\cdot {\psi _{v}} {{\,\mathrm{Tr}\,}}((DF)^{-1} D {\delta _{w}}) \det (DF) \, \mathrm{d}x \nonumber \\&- \int _\Omega p {{\,\mathrm{Tr}\,}}( D {\psi _{v}} (D F)^{-1} D {\delta _{w}} (D F)^{-1}) \det (DF) \, \mathrm{d}x \nonumber \\&+ \int _\Omega p {{\,\mathrm{Tr}\,}}( D {\psi _{v}} (D F)^{-1}) {{\,\mathrm{Tr}\,}}((DF)^{-1} D {\delta _{w}}) \det (DF) \, \mathrm{d}x \nonumber \\&+ \int _\Omega {\psi _{p}} {{\,\mathrm{Tr}\,}}(D v (D F)^{-1} D {\delta _{w}} (D F)^{-1}) \det (DF) \, \mathrm{d}x \nonumber \\&- \int _\Omega {\psi _{p}} {{\,\mathrm{Tr}\,}}( D v (D F)^{-1}) {{\,\mathrm{Tr}\,}}((DF)^{-1} D {\delta _{w}}) \det (DF) \, \mathrm{d}x \nonumber \\&- \int _\Omega ( D {\delta _{w}} + D {\delta _{w}}^\top ) : D {\psi _{w}} + {\eta _{\text {ext}}}((D{\delta _{w}}\, w)+(Dw\, {\delta _{w}}))\cdot {\psi _{w}}\, \mathrm{d}x \nonumber \\&+\beta \int _\Omega ({\eta _{\text {det}}}- \det (DF))_+ {{\,\mathrm{Tr}\,}}((DF)^{-1} D {\delta _{w}}) \det (DF) \, \mathrm{d}x\nonumber \\&- {\psi _{{\text {bc}}}} \cdot \int _\Omega {\delta _{w}} \det (DF) + (x + w) {{\,\mathrm{Tr}\,}}((DF)^{-1} D {\delta _{w}}) \det (DF) \, \mathrm{d}x\nonumber \\&- {\psi _{{\text {vol}}}} \int _\Omega {{\,\mathrm{Tr}\,}}((DF)^{-1} D {\delta _{w}}) \det (DF) \, \mathrm{d}x = 0, \end{aligned}$$
(22)
$$\begin{aligned} {\mathcal {L}}_{{\psi _{w}}} {\delta _{{\psi _{w}}}}= & {} - \int _\Omega (D w + D w^\top ): D {\delta _{{\psi _{w}}}} + {\eta _{\text {ext}}}(Dw\, w)\cdot {\delta _{{\psi _{w}}}} \, \mathrm{d}x + \int _{\Gamma _{\text {obs}}}b\cdot {\delta _{{\psi _{w}}}} \, \mathrm{d}s = 0,\nonumber \\ \end{aligned}$$
(23)
$$\begin{aligned} {\mathcal {L}}_{v} {\delta _{v}}= & {} \nu \int _{\Omega } \left( D {\delta _{v}} (DF)^{-1} \right) : \left( D v (DF)^{-1} \right) \det (DF)\, \mathrm{d}x\nonumber \\&- \nu \int _{\Omega } \left( D {\delta _{v}} (DF)^{-1} \right) : \left( D {\psi _{v}} (DF)^{-1} \right) \det (DF) \,\mathrm{d}x\nonumber \\&-\int _{\Omega } (D {\delta _{v}} (DF)^{-1}v)\cdot {\psi _{v}} +(D v (DF)^{-1}{\delta _{v}})\cdot {\psi _{v}}\det (DF)\, \mathrm{d}x\nonumber \\&- \int _{\Omega } {\psi _{p}} {{\,\mathrm{Tr}\,}}(D {\delta _{v}}(DF)^{-1}) \det (DF) \,\mathrm{d}x = 0, \end{aligned}$$
(24)
$$\begin{aligned} {\mathcal {L}}_{{\psi _{v}}} {\delta _{{\psi _{v}}}}= & {} -\nu \int _{\Omega } \left( D v (DF)^{-1} \right) : \left( D {\delta _{{\psi _{v}}}} (DF)^{-1} \right) \det (DF) \,\mathrm{d}x \nonumber \\&-\int _{\Omega } (D v (DF)^{-1}v)\cdot {\delta _{{\psi _{v}}}} \det (DF)\, \mathrm{d}x\nonumber \\&+ \int _{\Omega } p {{\,\mathrm{Tr}\,}}(D {\delta _{{\psi _{v}}}}(DF)^{-1}) \det (DF) \,\mathrm{d}x = 0, \end{aligned}$$
(25)
$$\begin{aligned} {\mathcal {L}}_{p} {\delta _{p}}= & {} - \int _\Omega {\delta _{p}} {{\,\mathrm{Tr}\,}}\left( D {\psi _{v}} (DF)^{-1} \right) \det (DF)\, \mathrm{d}x = 0, \end{aligned}$$
(26)
$$\begin{aligned} {\mathcal {L}}_{{\psi _{p}}} {\delta _{{\psi _{p}}}}= & {} \int _\Omega {\delta _{{\psi _{p}}}} {{\,\mathrm{Tr}\,}}\left( D v (DF)^{-1} \right) \det (DF)\, \mathrm{d}x = 0, \end{aligned}$$
(27)
$$\begin{aligned} {\mathcal {L}}_{b} {\delta _{b}}= & {} - \int _{\Gamma _{\text {obs}}}{\delta _{b}}\cdot {\psi _{b}} + D_{{\Gamma _{\text {obs}}}} {\delta _{b}} : D_{{\Gamma _{\text {obs}}}} {\psi _{b}} \, \mathrm{d}s + \int _{\Gamma _{\text {obs}}}{\delta _{b}} \cdot {\psi _{w}} \, \mathrm{d}s = 0, \end{aligned}$$
(28)
$$\begin{aligned} {\mathcal {L}}_{{\psi _{b}}} {\delta _{{\psi _{b}}}}= & {} - \int _{\Gamma _{\text {obs}}}b\cdot {\delta _{{\psi _{b}}}} + D_{{\Gamma _{\text {obs}}}} b : D_{{\Gamma _{\text {obs}}}} {\delta _{{\psi _{b}}}} \, \mathrm{d}s + \int _{\Gamma _{\text {obs}}}{c}n \cdot {\delta _{{\psi _{b}}}} \, \mathrm{d}s = 0, \end{aligned}$$
(29)
$$\begin{aligned} {\mathcal {L}}_{{c}} {\delta _{{c}}}= & {} \alpha \int _{\Gamma _{\text {obs}}}{c}{\delta _{{c}}} \, \mathrm{d}s + \int _{\Gamma _{\text {obs}}}{\delta _{{c}}} n \cdot b \, \mathrm{d}s = 0, \end{aligned}$$
(30)
$$\begin{aligned} {\mathcal {L}}_{{\psi _{{\text {vol}}}}} {\delta _{{\psi _{{\text {vol}}}}}}= & {} - {\delta _{{\psi _{{\text {vol}}}}}} \int _{\Omega } \det (DF) - 1 \, \mathrm{d}x = 0, \end{aligned}$$
(31)
$$\begin{aligned} {\mathcal {L}}_{{\psi _{{\text {bc}}}}} {\delta _{{\psi _{{\text {bc}}}}}}= & {} - {\delta _{{\psi _{{\text {bc}}}}}} \cdot \int _{\Omega } (x + w) \, \det (DF) \, \mathrm{d}x = 0, \end{aligned}$$
(32)

for all test functions \({\delta _{w}} \in W\), \({\delta _{v}} \in V\), \({\delta _{p}} \in Q\), \({\delta _{b}} \in H^2({\Gamma _{\text {obs}}})\), \({\delta _{{\psi _{w}}}}\in W\), \({\delta _{{\psi _{v}}}} \in V_0\), \({\delta _{{\psi _{p}}}} \in Q\), \({\delta _{{\psi _{b}}}} \in H^2({\Gamma _{\text {obs}}})\), \({\delta _{{c}}}\in L^2(\Gamma _{\text {obs}})\), \({\delta _{{\psi _{{\text {vol}}}}}}\in {\mathbb {R}}\), and \({\delta _{{\psi _{{\text {bc}}}}}}\in {\mathbb {R}}^d\).

Proof

The derivatives of \({\mathcal {L}}\) are obtained by utilizing standard rules of differentiation and taking into account the definition of a domain transformation mapping \(F= {{\,\mathrm{id}\,}}+ w\). Note that we particularly use the following identities

$$\begin{aligned} \frac{\partial \det (DF)}{\partial w} {\delta _{w}} = {{\,\mathrm{Tr}\,}}((DF)^{-1} {\delta _{w}})\det (DF) \quad \text { and }\quad \frac{\partial (DF)^{-1}}{\partial w} {\delta _{w}} = - (DF)^{-1} D{\delta _{w}} (DF)^{-1}. \end{aligned}$$

For the derivative of the penalty term, we utilize that

$$\begin{aligned} \frac{\partial }{\partial w}(( {\eta _{\text {det}}}- \det (DF) )_+ )^2 (w)\, {\delta _{w}}= & {} 2 ( {\eta _{\text {det}}}- \det (DF) )_+ \, \chi _{\lbrace {\eta _{\text {det}}}> \det (DF)\rbrace }\, \frac{\partial }{\partial w} \det (DF)(w) {\delta _{w}}\\= & {} 2 ( {\eta _{\text {det}}}- \det (DF) )_+ \, {{\,\mathrm{Tr}\,}}((DF)^{-1} {\delta _{w}})\det (DF). \end{aligned}$$

\(\square \)

Recall that the condition \(\det (DF) \ge {\eta _{\text {det}}}\) in the problem formulated in (3) is realized via the penalty term \(\frac{\beta }{2} \int _{\Omega } (( {\eta _{\text {det}}}- \det (DF) )_+ )^2 \, \mathrm{d}x\). The corresponding term in (22) of the optimality system in (3.1) is non-differentiable due to the positive-part function \((\cdot )_+\). Following the discussions in [13, sec. 3.5], the mapping \(w \mapsto -\beta \int _\Omega ({\eta _{\text {det}}}- \det (DF))_+ {{\,\mathrm{Tr}\,}}((DF)^{-1} D {\delta _{w}}) \det (DF) \, \mathrm{d}x\) is semismooth and one can compute an element from the generalized derivative in direction \(\delta ^\prime \) as

$$\begin{aligned}&({\delta _{w}}, \delta ^\prime ) \mapsto \beta \int _\Omega \chi _{\lbrace {\eta _{\text {det}}}> \det (DF)\rbrace } {{\,\mathrm{Tr}\,}}((DF)^{-1} D \delta ^\prime ) {{\,\mathrm{Tr}\,}}((DF)^{-1} D \delta ^\prime ) \det (DF)^2\\&\quad + ({\eta _{\text {det}}}- \det (DF))_+ {{\,\mathrm{Tr}\,}}((DF)^{-1} D\delta ^\prime (DF)^{-1} D {\delta _{w}}) \det (DF)\\&\quad - ({\eta _{\text {det}}}- \det (DF))_+ {{\,\mathrm{Tr}\,}}((DF)^{-1} D {\delta _{w}}) {{\,\mathrm{Tr}\,}}((DF)^{-1} D\delta ^\prime ) \det (DF)\, \mathrm{d}x. \end{aligned}$$

In the following, we briefly present a solution algorithm for the optimality system (22)–(32). For this purpose, we pursue a similar approach as in [13]. The core of this method is to solve the nonlinear shape optimization problem (17) for a decreasing sequence of regularization parameters \(\alpha _k\), starting from \(\alpha _0 = {\alpha _{\text {init}}}\) until the desired level \({\alpha _{\text {target}}}\) is reached. Because each subsequent optimization problem k is nonlinear, this approach benefits from utilizing the known values \(y_k :=(w, v, p, b, {c}, {\psi _{w}},{\psi _{v}},{\psi _{p}}, {\psi _{b}}, {\psi _{{\text {vol}}}}, {\psi _{{\text {bc}}}})_{k}\) as initial guess in the \((k+1)\)-th iteration. Algorithm 1 summarizes this procedure. Since parts of the optimality system are non-differentiable, we apply a semismooth Newton’s method.

figure a

In Sect. 4, we show how to choose the parameter \({\alpha _{\text {init}}}, {\alpha _{\text {dec}}}\) and \({\alpha _{\text {target}}}\) and illustrate their influence.

4 Numerical Results

This section is devoted to different numerical case studies of stationary, incompressible Navier–Stokes shape optimization problems. The purpose is to illuminate features of the nonlinear extension operator S proposed in Sect. 3. In particular, the benefit for optimization benchmark problems, which involve large deformations from the reference configuration to optimal shapes, is numerically investigated. It is moreover discussed how the local injectivity can be extended to globally injective transformation mappings by adding an artificial volume to the aerodynamic specimen. Furthermore, algorithmic solvability of the optimality system (22)–(32) is addressed in the end of this section.

The experimental settings for the tests are chosen to be comparable in 2d and 3d, respectively. The holdall domain \({G}\in \lbrace {G}_\text {2d}, {G}_\text {3d} \rbrace \), which reflects the flow tunnel in the experiment, is chosen as

$$\begin{aligned} {G}_\text {2d} :=[-7, 7 ]\times [-3, 3 ]\;\text { and }\; {G}_\text {3d} :=\lbrace x \in {\mathbb {R}}^3: -7 \le x_1 \le 7, \sqrt{x_2^2 + x_3^2} \le 3 \rbrace . \end{aligned}$$

Let \(\rho \) denote the diameter of the flow tunnel G. We then fix the velocity at the inflow boundary \({\Gamma _{\text {in}}}\) by \(v_\infty = \left( \cos (\frac{2 \pi \Vert x \Vert _2}{\rho }), 0, \dots , 0\right) ^\top \in {\mathbb {R}}^d\). In all experiments where the specimen \({\Omega _{\text {obs}}}\) is a circle or sphere, the radius is given by \(r = \mathrm{0.5}\) and \({\text {bc}}({\Omega _{\text {obs}}}) = 0\in {\mathbb {R}}^d\).

The discretization of all appearing PDEs is carried out with standard, piecewise linear P1 finite elements. In order to guarantee stability, we follow the pressure stabilized Petrov Galerkin approach (see for instance [16]), which utilizes an additional term for the pressure p and its adjoint variable \({\psi _{p}}\). The system under consideration is thus enriched by the two equations

$$\begin{aligned} g_p&:=\mu \sum _{T \in \mathcal {T}_h} h_T^2 \int _T \left( (DF)^{-1}\nabla p\right) \cdot \left( (DF)^{-1}\nabla {\delta _{p}} \right) \det (DF)\, \mathrm{d}x\\ g_{\psi _{p}}&:=\mu \sum _{T \in \mathcal {T}_h} h_T^2 \int _T \left( (DF)^{-1}\nabla {\psi _{p}}\right) \cdot \left( (DF)^{-1}\nabla {\delta _{{\psi _{p}}}} \right) \det (DF)\, \mathrm{d}x \end{aligned}$$

where \(\mathcal {T}_h\) denotes the set of all finite elements and \(h_T\) measures the longest edge of element T. For each of the subsequent experiments, \(\mu = \mathrm{0.1}\) is chosen.

All computations related to the finite element method are carried out using the GETFEM++ library [24]. We utilize a parallel version of the library, which relies on PARMETIS [18] for mesh partitioning and load balancing. All linear systems are handled via the parallel factorization solver MUMPS [1]. Both 2d and 3d discretization meshes are produced with the GMSH toolbox [10] and the Delaunay algorithms therein. If not stated otherwise, all 2d experiments follow the strategy of Algorithm 1 with the choice \({\alpha _{\text {init}}}=\mathrm{1e-4}\), \({\alpha _{\text {dec}}}=\mathrm{1e-1}\) and \({\alpha _{\text {target}}}=\mathrm{1e-10}\).

4.1 Non-convex Shapes with Large Deformations

Fig. 2
figure 2

Magnitude of velocity v computed on reference domain \(\Omega \) (left) and deformed \(F(\Omega )\) (right)

Full size image
Fig. 3
figure 3

Mesh quality preserving shape optimization experiment with large deformation from reference to optimal configuration

Full size image

Our first numerical study demonstrates the nonlinear extension equation for large deformations of a non-convex shape. The reference domain \(\Omega \) is chosen such that the specimen is described by a B-spline curve \({\Gamma _{\text {obs}}}\) given in terms of 6 control nodes. The situation is depicted in Figs. 2 and 3.

The relevance of this test case is to investigate the performance of the proposed approach for reference domains where the normal vector field n does not homogeneously point in all directions as for a circular shape. The influence of the normal vector is significant since it initially links the scalar-valued control \({c}\) to a vector-valued quantity as can be seen in (18). Section 4.3 is devoted to an experiment where several directions are underrepresented in the discretization of the normal vector n due to the shape of \({\Omega _{\text {obs}}}\).

In Fig. 2, the magnitude of velocity \(\Vert v \Vert _2\), computed in the undeformed state \(\Omega \), i.e., when \({c}=0\), is depicted. The right-hand side shows the velocity according to deformation \(F= {{\,\mathrm{id}\,}}+ w\) in terms of the optimal control \({c}\) after solving the optimality system given by (22)–(32). Furthermore, the optimal mapping \(F\) can be seen in Fig. 3 in the displacement and deformation of discretization elements. The relocation of triangles shows the effect of the nonlinear advection in the extension operator. A deeper look on this effect and the resulting mesh quality is provided in Sect. 4.5.

In this experiment, the viscosity of the fluid is chosen as \(\nu = \mathrm{0.01}\). The holdall domain \({G}\) is as described above. It consists of 382 surface segments on \({\Gamma _{\text {obs}}}\) and 10196 triangles in \(\Omega \). Barycenter and volume of \({\Omega _{\text {obs}}}\) in the reference configuration are given by \({\text {bc}}({\Omega _{\text {obs}}}) = (\mathrm{0.0307784},\mathrm{-0.035759})^\top \) and \( {\text {vol}}({\Omega _{\text {obs}}}) = \mathrm{0.809041}\). Thus, an optimal shape also undergoes a small translation since \({\text {bc}}(F({\Omega _{\text {obs}}})) = 0\) is required.

The essential settings in terms of shape optimization are the parameter \({\eta _{\text {det}}}\) and \({\eta _{\text {ext}}}\). Here, we choose \({\eta _{\text {ext}}}= \mathrm{3.0}\), which leads to the condition \(\det (DF) \ge {\eta _{\text {det}}}\) with \({\eta _{\text {det}}}= \mathrm{5e-2}\) being inactive. We can explain the homogeneously and smoothly deformed mesh due to this fact. In contrast, Sect. 4.2 shows examples where \(\det (DF) \ge {\eta _{\text {det}}}\) is active close to the tip of the optimal shape and how the displacement w and thereby the mesh quality in \(F(\Omega )\) is affected.

4.2 Influence of Factors \({\eta _{\text {det}}}\) and \({\eta _{\text {ext}}}\)

Fig. 4
figure 4

Optimal shapes \(F({\Gamma _{\text {obs}}})\) in terms of different \({\eta _{\text {det}}}\in \lbrace 0.5, 0.25, 0.2, 0.1 \rbrace \). The value \({\eta _{\text {det}}}\) decreases from top to bottom

Full size image
Fig. 5
figure 5

Experiments with \({\eta _{\text {det}}}= 0.05\) constant and variable \({\eta _{\text {ext}}}\) factor

Full size image

In this section, we visualize the influence of the choice of \({\eta _{\text {det}}}\) and \({\eta _{\text {ext}}}\) on the optimization. This illustrates how the set of admissible shapes \(\mathcal {F}_\text {adm}\) is determined thereby. The underlying experiment is a flow in \({G}\) over a circular specimen as described at the beginning of Sect. 4. The viscosity is again chosen to be \(\nu = \mathrm{0.01}\). The domain is discretized with 244 segments on \({\Gamma _{\text {obs}}}\) and 12,640 triangles in \(\Omega \).

First, we observe the influence of \({\eta _{\text {det}}}\) on the set of admissible shapes \(\mathcal {F}_\text {adm}\). Figure 4 visualizes how the condition \(\det (DF) \ge {\eta _{\text {det}}}\) acts on the optimal shape \(F({\Gamma _{\text {obs}}})\). Here, we choose \({\eta _{\text {det}}}\in \lbrace 0.5, 0.25, 0.2, 0.1 \rbrace \) beginning with the largest and then decreasing values. This condition can be interpreted such that the allowed, local change of volume in \(\Omega \) is relaxed from top to bottom in Fig. 4. In this experiment, it turns out that in the last computation with \({\eta _{\text {det}}}= \mathrm{0.1}\) the condition is inactive. Here, \({\eta _{\text {ext}}}= \mathrm{3.0}\) is fixed in all computations. With decreasing \({\eta _{\text {det}}}\), we observe a decrease in the objective function J (cf. (17)) as shown in Table 1. For comparison, J evaluated for the reference shape with \({c}=0\), and thus \(F= {{\,\mathrm{id}\,}}\), is 0.285233.

The next experiment follows the same setup with the only difference that \({\eta _{\text {det}}}= \mathrm{5e-2}\) is now fixed and \({\eta _{\text {ext}}}\in \lbrace \mathrm{0.0}, \mathrm{0.25}, \mathrm{0.5}, \mathrm{1.0}, \mathrm{2.0}, \mathrm{3.0}\rbrace \) takes increasing values. The mesh deformations, resulting from the optimal solution \(F= {{\,\mathrm{id}\,}}+ w\), are visualized in Fig. 5. Each of the subfigures shows a clip of size \(\mathrm{0.1}\times \mathrm{0.08}\) around the tip of the optimal shape. Similarly to the previous experiment, but now with increasing factor \({\eta _{\text {ext}}}\), we observe a decrease in the objective function J as shown in Table 2.

Table 1 Optimal objective function values depending on \({\eta _{\text {det}}}\)
Full size table
Table 2 Optimal objective function values depending on \({\eta _{\text {ext}}}\)
Full size table

Besides the significantly improved mesh qualities and more adequate set \(\mathcal {F}_\text {adm}\), we also observe that the Newton solver benefits from the appropriate choice of \({\eta _{\text {det}}}\) and \({\eta _{\text {ext}}}\). As soon as the condition \(\det (DF) \ge {\eta _{\text {det}}}\) becomes active, the optimality system (22)–(32) is not differentiable any further and the solver switches to semismooth Newton’s method. This effect is already documented in [13] for the case of Stokes flows and linear extension operator S.

Fig. 6
figure 6

Close-up visualization of the discretization around the tip of the ellipse experiment

Full size image
Fig. 7
figure 7

Flow over ellipsoidal reference shape and optimal solution with visualization of auxiliary grid in \({\Omega _{\text {obs}}}\). Color denotes norm of velocity field \(\Vert v \Vert _2\)

Full size image

4.3 Extending the Local-Only Injectivity

This section focuses on an effect that is likely to appear for non-spherical reference domains. In particular, for the aerodynamic experiments considered in this article it might happen that the upper surface of the obstacle overlaps the lower one. Especially for large deformations from reference to optimal shape and for shapes that are stretched parallel to the flow axis, we encounter effects as depicted in Fig. 6a for the experiment shown in Fig. 7. In other words, \(\Omega \) \(\mapsto F(\Omega )\) is not globally injective in this situation. This is due to the fact that the condition \(\det (DF)\) ensures injectivity of \(F\) only locally but not globally.

In the following, we propose a modification of the extension operator S in order to extend the injectivity. Recall that in the setting followed up to here the obstacle domain \({\Omega _{\text {obs}}}\) is treated as void and there is no discretization within. We now consider the operator S on the entire holdall domain \({G}= \Omega \cup {\Omega _{\text {obs}}}\) in contrast to the state equation that remains in \(\Omega \). Moreover, the condition \(\det (DF) \ge {\eta _{\text {det}}}> 0\) is now required on \({G}\). Thus, the displacement field is defined by \(w \in H_0^1({G}, {\mathbb {R}}^d)\). Moreover, we reformulate the weak formulation of the extension operator S given in (20) to

$$\begin{aligned} \int _{{G}} (Dw + Dw^\top ):D{\delta _{w}} + {\eta _{\text {ext}}}(Dw\, w)\cdot {\delta _{w}} \, \mathrm{d}x = \int _{\Gamma _{\text {obs}}}b {\delta _{w}} \, \mathrm{d}s \end{aligned}$$
(33)

for all \({\delta _{w}} \in H_0^1({G}, {\mathbb {R}}^d)\) and appropriate \({\eta _{\text {ext}}}\ge 0\). Simultaneously, the penalty term, which enforces the local injectivity, in (21) changes to

$$\begin{aligned} \frac{\beta }{2} \int _{{G}} (( {\eta _{\text {det}}}- \det (DF) )_+ )^2 \, \mathrm{d}x. \end{aligned}$$
(34)

For this experiment, we choose \(\nu = \mathrm{0.01}\) as in the previous sections. The holdall \({G}\) has the same outer dimension, and the specimen \({\Gamma _{\text {obs}}}\) is an ellipse with semimajor-axis \(r_1=\mathrm{2.7}\), semiminor-axis \(r_2=\mathrm{0.2}\) and barycenter \({\text {bc}}({\Omega _{\text {obs}}}) = (0, 0)^\top \). Its surface is subdivided in 884 segments. Further, \({G}\) is discretized by 36,360 triangles, 29,930 in \(\Omega \) and 6430 in \({\Omega _{\text {obs}}}\).

Figure 6 visualizes the effect of the mapping \(F\) on the discretization grid. In Fig. 6a, the optimal solution for \(\alpha =\mathrm{1e-2}\) is shown. Note that in this particular case we stop the optimization for a larger value, since this already leads to singularities. Figure 6 depicts \(\det (DF)\), which is again bound away from zero by \({\eta _{\text {det}}}=\mathrm{5e-2}\). It can be seen that, although this condition is inactive, the non-injective mapping cannot be prevented.

The same experiment is then conducted with the changes proposed in the beginning of this section, which leads to the values of \(\det (DF)\) shown in Fig. 6c. Now \({\Omega _{\text {obs}}}\) is discretized and S also acts on the interior of the specimen. Here, the optimization is performed with the setting \({\alpha _{\text {init}}}=\mathrm{1e-4}\), \({\alpha _{\text {dec}}}=\mathrm{5e-1}\) and \({\alpha _{\text {target}}}=\mathrm{1e-10}\). The resulting optimal solution is visualized in Fig. 7 where Fig. 7a shows the reference domain and the velocity field computed for this configuration. Figure 7b depicts the domain \(F({G})\) and the velocity field computed on \(F(\Omega )\). The result illustrates that by discretizing the obstacle domain \({\Gamma _{\text {obs}}}\) we have prevented elements from overlapping. Note that the relatively fine grid is chosen at the front and the back of the shape due to the large curvature of \({\Gamma _{\text {obs}}}\) in these regions. This experiment turns out to be more challenging than, e.g., a spherical reference shape since on coarse grids the normal vector field in these areas tends to be underresolved. From a computational point of view, it is attractive to have a coarse grid in \({\Omega _{\text {obs}}}\), as chosen in the center of the specimen, to reduce cost for the solution of the operator S.

4.4 Three-Dimensional Results

In this section, we perform a three-dimensional optimization experiment as an illustration of the proposed methodology. As already observed in [13] for the Stokes experiment, more care has to be taken for the decrease strategy of \(\alpha \) in Algorithm 1. Especially the semismooth Newton solver shows to be challenging w.r.t. to convergence when the condition \(\det (DF) \ge {\eta _{\text {det}}}\) becomes active.

The experiment shown in Fig. 8 is within the framework described at the beginning of Sect. 4. The flow tunnel \(\Omega \) is discretized by 632,093 tetrahedrons and the surface of the spherical specimen in the reference configuration \({\Gamma _{\text {obs}}}\) consist of 8558 triangles. Further, the viscosity is chosen to be \(\nu =\mathrm{0.01}\) and in Algorithm 1 we set \({\alpha _{\text {init}}}= \mathrm{1e-4}\), \({\alpha _{\text {dec}}}= \mathrm{5e-1}\) and \({\alpha _{\text {target}}}=\mathrm{1e-6}\). The results shown here are obtained with an extension factor of \({\eta _{\text {ext}}}=15\) and \({\eta _{\text {det}}}= 0.05\). We visualize the impact of the optimization on the fluid by stream lines of the velocity field in Fig. 8. This figure also shows the effect of the particular operator S on the quality of the surface mesh when it undergoes the optimal deformation \(F\). The combination of Laplace–Beltrami (19) and the nonlinear extension equation (19) leads to a homogeneous distribution of triangles on the surface \({\Gamma _{\text {obs}}}\). It can be observed that this is due to tangential components in \(w\vert _{\Gamma _{\text {obs}}}\). This is a benefit of a vector-valued extension equation over approaches which utilize a static extension of the normal vector field in order to extend the boundary control to the surrounding volume. Furthermore, Fig. 9 shows a zoom-in to the tip of the deformed domain \(F(\Omega )\). Here, we can see a crinkled clip in the \(x_1 x_2\)-plane with \(x_3=0\), which shows the quality of the tetrahedrons.

Fig. 8
figure 8

Velocity stream lines computed on reference and optimal domain together with visualization of surface discretization

Full size image
Fig. 9
figure 9

Crinkled clip of \(x_1 x_2\)-plane with \(x_3=0\) showing deformed mesh \(F(\Omega )\) together with surface elements \(F({\Gamma _{\text {obs}}})\)

Full size image

4.5 Quantification of the Influence of \({\eta _{\text {ext}}}\) on Mesh Quality

This section presents numerical experiments which investigate the influence of the nonlinear extension operator S on the mesh quality in 2d and 3d (cf. Fig. 10). Recall that the discretization mesh is not actually deformed within the optimization. In principle, mesh deformations are not a part of our method. Despite that, in this section, we apply mesh deformation and compute the mesh quality to illustrate the properties of the extension operator. The mesh quality is important in case one would like to further perform the finite element simulations. The 2d experiment is conducted on the same computational domain as before with a circular specimen, 312 surface segments and 6168 triangles in \(\Omega \). The fluid viscosity is chosen to be \(\nu = \mathrm{0.01}\). Figure 10a visualizes the influence of \({\eta _{\text {ext}}}\in [0, 1 ]\) on the mesh quality of \(F(\Omega )\). It is measured by the ratio of radii of largest inscribed and smallest circumscribed circle, where the plot shows the value of the worst triangle.

This experiment quantifies the effect which is already visualized in Fig. 5. For a shape optimization with large deformations from reference to optimal configuration, i.e., \(\Vert w \Vert _{L^2({\Gamma _{\text {obs}}})}\) is relatively large, a pure linear extension operator S does not reliably lead to satisfying mesh qualities. Moreover, it can be seen that in this particular experiment there is a saturation effect of the nonlinearity in S starting at approximately \({\eta _{\text {ext}}}\approx \mathrm{1.5}\). Figure 10b shows the results of a similar experiment in 3d. Here, a mesh is chosen with 6040 surface triangles on \({\Omega _{\text {obs}}}\) and 147 385 tetrahedrons in \(\Omega \). Note that we decrease the viscosity to \(\nu = 0.1\) in this experiment in order to be able to obtain results for \({\eta _{\text {ext}}}< 3.0\). In 3d, quality is measured by the radius ratio of smallest circumscribed sphere to the largest inscribed one. Again the worst element is visualized. Also note that the y-axis is in log-scale. In this setting, it turns out that the effect of compressed cells near the tip and back of the shape, which is stretching due to a decrease in \(\alpha \), is stronger than in 2d. Consequently, the solver fails to converge after reaching a certain value of \(\alpha \). This can be explained by the semismoothness in the optimality system that becomes active in a significant number of finite elements in this situation. However, it can be observed that, starting with approximately \({\eta _{\text {ext}}}\approx 8\), a saturation is possible, where the mesh quality of \(F(\Omega )\) remains adequate for further numerical computations.

Fig. 10
figure 10

Quality of worst element (triangle in 2d or tetrahedron in 3d) after applying optimal deformation for a range of extension factors \({\eta _{\text {ext}}}\)

Full size image

4.6 An Iterative Optimization Algorithm

figure b
Fig. 11
figure 11

Iterative solution strategy according to Algorithm 2 with \({\alpha _{\text {init}}}=1\), \({\alpha _{\text {dec}}}=0.5\) and \({\alpha _{\text {target}}}=\mathrm{2e-7}\)

Full size image

In the previous sections, we solve the nonlinear, non-smooth optimality system with the direct solution strategy given in Algorithm 1. Moreover, a direct solver library is applied to the resulting linear systems within semismooth Newton’s method. This approach is clearly limited due to the high memory requirement. Especially, when the state equation results from a time-dependent problem, this procedure becomes impracticable. Hence, in this section we focus on a numerical study of decoupling system (22)–(32). This approach is summarized in Algorithm 2.

We demonstrate that it is possible to decouple the solution process of state (25) and (27), adjoint (24) and (26) and shape related equations, i.e., (22), (23) and (28), (29), (30), (31), (32), from each other. On the one hand, this allows to reuse existing solvers for the state equation and embed them into the shape optimization framework. On the other, the memory requirement for linear solvers significantly reduces. Moreover, the semismooth part (22) is split from the other equations and a solver can be particularly tailored for this purpose.

Algorithm 2 operates on the nonlinear optimality system as a fixed-point strategy. In an outer loop, it is again iterated over a decreasing regularization parameter \(\alpha \) as in Algorithm 1. Thus, approximate solutions for the optimization problem according to \(\alpha _k\) are utilized as initial guess for the nonlinear solver in iteration \(k+1\). Yet, unlike in the direct approach, the subproblems are only solved approximately by a fixed-point iteration, which solves the decoupled equations of the optimality system in turns. The termination criterion for this inner loop is the relative change in the control variable \({c}\) measured in the \(L^2({\Gamma _{\text {obs}}})\)-norm.

In Fig. 11, the results of one run of Algorithm 2 are shown. The underlying optimization experiment is a 2d computation on the same grid as in Sect. 4.5 with 312 surface segments and 6168 triangles in \(\Omega \), \({\alpha _{\text {dec}}}= 0.5\), \({\alpha _{\text {init}}}= 1.0\), \({\alpha _{\text {target}}}=\mathrm{2e-7}\), \(\nu = 0.1\) and \({\eta _{\text {ext}}}= 1.5\). Note that the initial value of \(\alpha \) is significantly larger than the choices made for Algorithm 1. Figure 11 shows the required inner iterations until the condition

$$\begin{aligned} \frac{\Vert {c}_{\ell +1} - {c}_\ell \Vert _{L^2({\Gamma _{\text {obs}}})}}{\Vert {c}_{\ell +1}\Vert _{L^2({\Gamma _{\text {obs}}})}} < \epsilon \end{aligned}$$

is fulfilled for \(\epsilon = \mathrm{1e-2}\). Furthermore, the value of the objective J (cf. (17)) is visualized. It is computed in Algorithm 2 in line 10 at the end of one inner loop. Notice the jumps in the objective function between iteration 5 and 20. In our experiments, it turns out that this is an effect that both influences the minimal possible \({\alpha _{\text {dec}}}\) and \({\alpha _{\text {init}}}\).

In this setting, a total number of 53 inner iterations, i.e., solutions of the state equation, are required to reach the optimal shape. This numerical study can thus be seen as an illustration of how to reduce the computational costs of the large, coupled, nonlinear system (22)–(32). Thus, the proposed method is applicable to more complex problems, such as non-stationary Navier–Stokes flows.

Furthermore, despite the high computational cost of the system with the proposed nonlinear extension equation, the overall complexity it not increased significantly. This is due to the fact that the model PDE in general is also nonlinear and time-dependent.

5 Conclusion

In this article, we have proposed and numerically illustrated choices of nonlinear extension operators within the method of mappings for aerodynamic shape optimization. These operators are based on the idea that an additional, nonlinear advection term leads to a rearrangement of discretization cells along the major direction of deformations.

The main goal we have achieved is to circumvent mesh degeneracy effects that appear under large deformations when the extension of the boundary control is chosen according to linear elastic models. Especially in the underlying aerodynamic drag minimization, where optimal shapes tend to become stretched in flow direction and compressed in the orthogonal directions, we have numerically investigated how mesh quality can be preserved.

We have also demonstrated one possibility to decouple the solution process of the optimality system in order to overcome issues of computational complexity. Moreover, we have studied how the set of admissible shapes depends on the nonlinearity of the operator and how the local injectivity of mappings can be extended to large deformations. Since the proposed methodology is formulated in function spaces without taking a specific discretization into account, another benefit of this approach is that it naturally allows to introduce concepts like adaptivity. An important field for future investigations is a detailed description of properties of the set \(\mathcal {F}_\text {adm}\), which is constructed in terms of the nonlinear extension operator S.