A novel p-harmonic descent approach applied to fluid dynamic shape optimization

We introduce a novel method for the implementation of shape optimization for non-parameterized shapes in fluid dynamics applications, where we propose to use the shape derivative to determine deformation fields with the help of the p-\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p-$$\end{document} Laplacian for p>2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p > 2$$\end{document}. This approach is closely related to the computation of steepest descent directions of the shape functional in the W1,∞-\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$W^{1,\infty }-$$\end{document} topology and refers to the recent publication Deckelnick et al. (A novel W1,∞\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$W^{1,\infty}$$\end{document} approach to shape optimisation with Lipschitz domains, 2021), where this idea is proposed. Our approach is demonstrated for shape optimization related to drag-minimal free floating bodies. The method is validated against existing approaches with respect to convergence of the optimization algorithm, the obtained shape, and regarding the quality of the computational grid after large deformations. Our numerical results strongly indicate that shape optimization related to the W1,∞\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$W^{1,\infty }$$\end{document}-topology—though numerically more demanding—seems to be superior over the classical approaches invoking Hilbert space methods, concerning the convergence, the obtained shapes and the mesh quality after large deformations, in particular when the optimal shape features sharp corners.


Introduction
Adjoint-based local optimization has been matured towards an efficient industrially applied strategy, e.g.[27,28].When attention is given to shape optimization, the aim is to find optimal shapes regarding a physical quantity J, e.g. the drag force experienced by an obstacle.Mathematically speaking, a shape functional J is minimized subject to partial differential equation (PDE) constraints.The latter typically govern the physics, e.g. the conservation of mass and momentum.A crucial part of the shape optimization procedure is the choice of the descent direction that provides an update rule for the design variable.In our case the design variable refers to the underlying geometry Ω, the boundary ∂Ω or parts of the boundary Γ.The descent direction is usually employed to converge a sequence of shape updates via a deformation field u.First attempts used the directional derivative J , also referred to as shape derivative in order to formulate a shape update rule [34,19,9].Because shapes are not elements of a vector space, e.g.there is no meaningful definition of the summation of two shapes, the shape derivative is defined by introducing deformations that make shapes variable and thus allows the definition of directional derivatives.With the help of shape calculus, cf.[34,9] for a mathematical perspective or [32,22] for an engineering perspective, a shape derivative can be computed that relates to a scalar field γ defined on the boundary Γ.We want to refer to γ as the local shape sensitivity.Mind that the computation of the shape derivative with shape calculus is a rather involved task and heavily depends on the objective functional J, the PDE constraints that apply as well as the domain Ω and the boundary ∂Ω, respectively.In [29] the shape is updated by using the local shape sensitivity γ to perform the deformation in normal direction n of the boundary Γ.Using the descent direction u = −γn yields a contribution u • n = −γ 2 to reducing the objective functional, and is a popular approach within shape optimization, cf.[41,42].The attempt is limited since it often yields shapes with rough/noisy boundaries [39,20] and distorted near-wall meshes which in turn hamper the preservation of numerical accuracy during the optimization procedure [37,5].In fluid dynamics and neighboring applications, the computational grid is frequently morphed and not renewed after each optimization step.Thereupon, attempts follow in order to gain higher regularity of the deformation.One approach is based on the definition of a shape gradient gradJ by an inner product and the shape derivative J .Here the shape gradient is identified by the Riesz representation of the directional derivative of the shape functional.Even though this leads to smoother deformations, the approach is algorithmic challenging due to solving a PDE on a hyperplane, and also mathematically questionable in the general case, see also [1] and cf.section 2. In this regard, different gradients are associated with different transformations applied to the shape derivative, and several techniques have been proposed to increase the regularity of the shape updates: a) CAD-related shape definitions connect the node-based shape derivatives to the CAD parameterization using the chain rule of differentiation, cf.[23,30].The procedure couples the various local derivatives and thereby ensures smooth shapes.However, the rigid finite dimensional initial CAD parameterization limits the attainable shapes and different CAD models may result in different optimal shapes.b) A coupling of mesh node updates using either local shape functions, e.g.FE-type functions [35,36], or global shape functions, e.g.Hicks-Henne approaches [15].c) A more rigor approach of Jameson and Vassberg [18,41,42] applies an implicit, continuous smoothing operator to either the shape derivative J or the deformation field u, based on an extended definition of the inner product, frequently labeled 'Sobolev-gradient'.Applying the smoothing operation on the surface leads to the Laplace-Beltrami operator and a related surface metric.For computational reasons, the practice is often performed in an explicit manner, cf.[5].The explicitly filtered local shape sensitivity, e.g. by using consistent kernel functions [20], marks a first-order approximation to the implicit Sobolev-gradient [39].
All strategies (a-c) essentially couple node updates and thereby obtain smooth design updates.Having updated the discrete design surface, the subsequent numerical investigation of the updated design also requires an update of the computational mesh, i.e. the shape gradient of the design surface needs to be extended into the domain.The habitat of the shape gradient depends on the surface metric and can be surface as well as volume based.Prominent examples refer to above mentioned Laplace-Beltrami (LB) or the Steklov-Poincaré (SP) metric, c.f. [33].The first approach (LB) exclusively operates in the tangent space of the design surface, the latter (SP) leads to a domain formulation, where results are subsequently projected on the controlled shape.The SP strategy gives the shape update of the design surface and mesh using the shape sensitivity [33,14,1].The volume based SP approach is particularly attractive for optimization procedures which prefer mesh morphing over re-meshing strategies, as it is customary for engineering simulations.While re-meshing can be automated, the lack of fair restart capabilities becomes prohibitively expensive in practical applications.Moreover, the use of standardized, HPC capable solution routines supplied by the flow solver (assembling, solving, etc.) represent another significant benefit of the SP approach.
A different avenue is taken by the phase field method for fluid mechanic shape optimization, where the shape of the sought domain is approximated by the zero level set of a phase field function.This turns the shape optimization problem into a PDE constrained optimization problem where the phase field enters as control in the coefficients of the PDE.This allows to apply the complete algorithmic machinery of PDE constrained optimization methods to this formulation of the fluid dynamic shape optimization problems, and this approach also naturally allows topology changes of the shape.However, numerical methods also for this approach encounter problems in situations where the sought shape needs to develop kinks and/or corners.This approach is proposed in [7], and in a couple of papers investigated for hydrodynamic shape optimization problems [11,12,13].
Although smooth shapes may be desirable for different reasons, they are not necessarily optimal.If the optimum involves a kink or a corner, the above mentioned strategies to obtain smooth shape updates display difficulties to capture such optima from curved initial configurations, if the respective region is (initially) not resolved by very fine grids.Section 4 of the present paper discusses a classical example of a pointed optimal shape.The same is true for counterpart situations, i.e. to transform an initially kinked shape into a curved optimum.Though this might be possible, the convergence is often fairly slow.
The present study aims to convey the merits of an alternative strategy to compute the shape deformation from shape derivatives in the context of CAD-free -aka node-based -shape optimization.The p-Laplace operator is used in a volume-based formulation along the route of the SP metric.The approach is industrially feasible and supports unstructured meshes.Applications refer to 2D and 3D fluid dynamic shape optimization ranging from laminar to turbulent external flows around free floating objects with fixed displacement.
The remainder of the paper is structured as follows: Section 2 outlines the mathematical framework and the rationale that leads to the p-Laplace problem to approximate the steepest descent direction within the W 1,∞ -topology together with a discussion of our fluid dynamic shape optimization problem.Section 3 presents the solution algorithm.Section 4 applies the approach to three test cases and the manuscript closes with conclusions in Sec. 5.

Mathematical Framework
In this section we want to outline the basic idea behind the p-harmonic approach and briefly recall the concept of shape optimization.
For this purpose let J : A → R denote a shape functional, where A denotes the set of admissible domains which has to be specified in the respective application.In our setting the set A is specified through Fig. 1, and the shape functional with (3) is given in (4a).For the algorithmic minimization of J for a given domain Ω ∈ A we intend to specify descent vector fields u * : R d → R d such that J (Ω; u * ) < 0 holds, where J (Ω; u * ) denotes the shape derivative of J at Ω in direction u * .The perturbed domain Ω then has the form Ω = T t (Ω) := (id + tu * )(Ω), where t > 0 is a step size specified in the respective minimization algorithm.However, descent in this context requires to specify an appropriate topology.Moreover, it frequently is required that Lipschitz domains Ω are mapped to Lipschitz domains Ω. Common practical approaches use Hilbert Space methods, i.e. seek descent vector fields u * determined with the help of the shape derivative by a(u * , w) = J (Ω; w) for all w ∈ H, where (H, a(•, •)) denotes an appropriate Hilbert space, see e.g.[1, for an extensive discussion of approaches related to the Hilbert space setting.This requires to compute the Riesz representative of the functional J (Ω; •).A typical choice of H is H m (Ω, R d ), where however m ∈ N has to be chosen large enough to obtain a Lipschitz transformation.A way around this would be to directly choose u * ∈ W 1,∞ (Ω, R d ) as a direction of steepest descent for J at Ω, where W 1,∞ (Ω, R d ) denotes the set of Lipschitz transformations from Ω to R d .This leads to the minimization problem (1) min which is challenging both from the mathematical and the numerical perspective, since it represents a variational problem in the non-reflexive Banach space W 1,∞ (Ω, R d ).
Variational problems of this kind are studied in e.g.[17], where it is proposed to approach solutions u * of problem (1) by a sequence of solutions u * p ∈ W for the mathematical analysis of the limit process p * → ∞.Note that for p = 2 we recover a Hilbert space setting as described above, and refer to problem (2) as p−Laplace relaxation of problem (1).In the next section we adapt problem (2) to our fluid dynamic setting and use its solutions as descent directions in our augmented Lagrange algorithm for the numerical solution of our shape optimization problem.
2.1.Optimization Problem.We now introduce the mathematical setting for our hydrodynamic shape optimization problem, where we refer to Fig. 1 for the geometrical setting and the notation.The governing equations here are either given by the stationary Navier-Stokes equations for incompressible fluids and laminar flows or the Reynolds averaged Navier-Stokes (RANS) equations which we consider for turbulent flow at high Reynolds numbers.Note that in the later case turbulence modeling is required.Our aim is to find the shape of a generic obstacle E ⊂ B with Lipschitz boundary located within the flow channel B ⊂ R d which has minimal drag.The state then is given by the velocity v : Ω → R d and pressure p : Ω → R, which are assumed to be unique on Ω, and thus a mapping Ω → y(Ω) = (v, p) exists.
It is considered that the domain Ω ⊂ B is completely filled with fluid which flows in from the left boundary Γ in with a prescribed velocity v ∞ towards the outflow boundary Γ out on the right.The top and bottom boundaries Γ slip are considered as frictionless slip wall boundaries.For all test cases, the whole boundary Γ of the . Schematic sketch of the computational domain.
obstacle E is considered to be the control of the total force acting on its boundary, viz. (3) where I is the identity tensor.However, additionally to the flow state also geometric constraints need to be considered.On the one hand, the geometry has to be held in place to prevent the body from moving during the optimization process.On the other hand the geometry would shrink to a single point if the volume is not preserved.Thus, in the problem formulation the location of the barycenter and the wetted volume Ω are prescribed to be the vector b ∈ R d and the constant c > 0, respectively.Both geometrical constraints supplement the PDE constraints and the optimization problem reads (4a) min where the local tangent direction is defined by t = t τ / t τ 2 with the tangent projection of the shear force the inflow velocity, and ρ, µ ∈ R + are the density and viscosity of the fluid.The vector containing the geometrical constraints from (4h)-(4i) is given by g = (b, c) T , where b and c refer to barycenter and volume residuals, respectively.To deal with the constraint problem above, a common approach is to introduce Lagrange multipliers.Therefore, we obtain the Lagrangian of problem (4a)-(4i) as where ρ b , ρ c ∈ R + are sufficiently large penalty factors.The variables (v, p) are the Lagrange multipliers for the PDE constraints (4b) and (4c) which also are associated with the adjoint state and λ ∈ R d is a Lagrange multiplier for considering the Dirichlet boundary condition (4d) that holds on the deformed boundary Γ.Note that when considering RANS equations the turbulence model implies the solution of additional state equations and thus corresponding Lagrange multipliers.Here we assume the influence of turbulence effects to be small which justifies the frequently employed frozen turbulence assumption, cf.[10,26,38].The multipliers λ b ∈ R d and λ c ∈ R, belonging to the geometric constraints (4h) and (4i), are determined with an augmented Lagrange method and considered to be constant during the shape optimization process, as described below in detail.Following the standard approach, as described in textbooks like [16,40], the penalized objective function can be expressed by the Lagrangian ( 5) and therewith the constraint problem (4a) -(4i) can be transformed into an unconstrained problem (7) min The first order optimality system for ( 7) is given by the partial derivative of ( 5) with respect to each argument.The Fréchet derivative, w.r.t. the adjoint state (v, p) gives the boundary value problem (4b) -(4g).Differentiation of ( 5) with respect to the state (v, p) and integration by parts leads to the adjoint equation system and with (8g) λ = pn − µ(∇v + ∇v T ) • n where τ = µ(∇v + ∇v T ) − Ip.Together with the PDE constraints of the optimization problem (4b) -(4g) the adjoint equation system (8a) -(8f) characterizes a saddle point of the Lagrangian (5), which is assumed to be a unique stationary point of L. In order to compute the shape derivative of the augmented Lagrangian (5) the domain Ω is made variable by a family of transformations {T t } t≥0 with the parameterized perturbation of identity where V 1,p 0 := {u ∈ W 1,p (Ω, R d ) : u = 0 a.e. on ∂Ω \ Γ}, and t ≥ 0 is a step size, cf.[22].Application of the first order optimality condition leads to the surface formulation of the shape derivative of the objective function J, compare e.g.[6].For its representation we define (11) is the shape derivative of J at Ω into direction of u.While the Lagrange multipliers (v, p) are given by the solution of the boundary value problem (8a) -(8f), the multipliers λ b and λ c for barycenter and volume constraint are not given by the solution of a system of PDEs and rather have to be determined by an augmented Lagrange method like outlined in algorithm 1 of the next section, c.f. [3,1].Hence, with the shape derivative from (13) the minimization problem in (2) with W 1,p (Ω, R d ) replaced by V 1,p 0 (Ω, R d ) is associated with the variational form to the p-Laplacian problem (14) −div (∇u p : ∇u p ) p−2    Note that Γ ⊂ ∂Ω is the part of the boundary which is free for deformation.This approach now can also be interpreted as a generalization by the p-Laplace setting of the approach proposed in [33] where the descent direction is defined regarding a Steklov-Poincaré operator; in other words the Dirichlet-to-Neumann map is applied to γn.Therewith, the descent direction in [2,33] is obtained by solving an elliptic boundary value problem similar to (15) −∆u = 0 in Ω ∂u ∂n = γn on Γ u = 0 on ∂Ω\Γ    which is the Euler-Lagrange equation of the minimization problem (16) min This obviously represents a special case of (2) with p = 2 and serves as a reference for the numerical experiments discussed in section 4. Let us also note that this concept was enhanced in [25] where a nonlinear extension operator has been introduced.The observation of strong distortion of the discrete grid within the main deformation direction motivated adding a nonlinear advection term to the PDE in (15).Therewith, greater deformations are possible and the mesh quality is reasonable even after large deformations of the initial grid.However, as our numerical results show, the p−Laplace relaxation of problem (1) provides a systematic approach to fluid dynamic shape optimization also guaranteeing meshes of high quality after large deformations, in particular when the optimal shape has sharp corners.

Optimization Algorithm
In Algorithm 1 we specify our minimization algorithm.It was used before, in

Algorithm 1 Augmented Lagrange Optimization
Require: end if 16: [33] to determine the multipliers λ b and λ c for the barycenter and volume constraints of a shape optimization problem.Other than in [33] the update for the penalty factors ρ b , ρ c ∈ R + and λ b ∈ R d and λ c ∈ R, respectively, is separated which is motivated by the big difference in the order of magnitude between both penalty factors ( b ∝ 10 8 and c ∝ 10 2 ) for the numerical experiments in this paper.However, the critical part of algorithm 1 refers to step 3 which is therefore farther explored in algorithm 2. To compute the multipliers λ b and λ c the shape optimization algorithm 2 is called in a sequence of convergence tolerances ∈ R + , c.f. step 9 in algorithm 2, in order to keep the computation numerically stable and circumvent unfeasible shapes during the optimization.
A conventional, pressure-based, second-order accurate finite-volume scheme for a cell-centered variable arrangement is employed to discretize the partial differential equations of the primal (4b) -(4g) and adjoint systems (8a) -(8f), cf.[31,38,21].The existing infrastructure and generic subroutines of the fluid solver can be compute γ according to (13) 6: solve p-Laplace problem ( 14) move grid points according to (9) re-used with limited effort for the implementation.In case of p = 2, the initial guess with u p = 0 in Ω leads to convergence of the implementation.For p > 2 a sequence of problems in p is solved in order to obtain an initial guess for the discrete problem with the desired value for p.Within the procedure the solution of the preceding smaller p is used as initial guess.Approximations that offer an initial guess only experience a smaller convergence tolerance in order to save computational time.

Applications
Three fluid dynamic applications are discussed to investigate the performance of the p-Laplace approach (2) in shape optimization.All cases are concerned with drag minimization at steady state and subjected to conserve the wetted volume and its barycenter, c.f. Fig. 1.Starting from the Laplace expression with p = 2, the first case analyses the influence of increasing p for a frequently referenced 2D Stokes flow example that features a pointed oval optimum.Emphasis is given to (a) the convergence of the optimization, (b) the final shape and the attainable drag reduction as well as (c) the quality of the mesh updates.The second case demonstrates the applicability of the approach for an analogue 3D configuration and analyses the same aspects (a-c).Since both initial applications refer to low Reynolds (low Re) number flows, a third example is supplemented to scrutinize the performance in a 2D turbulent flow at high Reynolds number.The initial grid is displayed in Fig. 2 and features 927 evenly distributed cells along the circumference of the cylinder boundary Γ.The CAD-free optimization procedure employs an unstructured, locally refined mesh with approximately 14.5k control volumes which is deformed along with the surface using (9).Mind that the optimal shape is expected to reveal pointy tips at the front and the aft, which shall deliberately not be preempted by the initial discretization of the boundary Γ.The target barycentre and wetted volume are set to b = (0, 0) T and c = 50 • 10 − π/4[m 2 ] respectively.All parameters employed to initialize the algorithms 1 and 2 are displayed in Table 1.The multipliers for the barycenter and volume are initialized with λ b = (0, 0) T and λ c = 0, respectively.The sequence of tolerances applied to the convergence criteria of the shape optimisation problem is given by = 10 −1 , 10 −2 , . . ., 10 −7 .However, it turned out that the multipliers λ b and λ c converge very fast and after four augmented Lagrange steps the multipliers are determined sufficiently exact in each computation with different values for p.The field values for the fixed point iteration for p = 2 can be initialized with u p = 0 for all discrete points.Larger p-values were initialized by solutions obtained from the previous smaller value, each to a suitable tolerance to provide an initial guess for the next following p-Laplace problem.Although the theory suggests to drive p → ∞, lower p-values are of interest for large-scale applications due to the more exhausting computational effort.The numerical effort for solving the p-Laplace problem is investigated in [24] for different values for p which show that it is of polynomial complexity but depends on p and the number of unknowns.The related experience from this investigation reveals an increase of computing time T p for one iteration of the shape optimization algorithm 2 by T p=4 /T p=3 ≈ 2.8 and T p=5 /T p=4 ≈ 1.6.In addition the representable floating point arithmetic of the machine limits the value for p.
4.1.1.Optimal Shapes and Convergence.Figure 3 shows the normalized evolution of the drag objective (6) over a selected number of optimization steps for the five investigated values of p.The baseline solution refers to p = 2 which is also related to the investigations in [33].Table 2 outlines a comparison of performance indicators obtained for the five investigated values of p.The last column refers to the maximum number of design steps that are needed to reach a sufficiently converging objective function (cf.algorithm 2).The table reveals that the convergence   improves and fewer optimization steps are needed for larger values of p.For p = 2 the optimization could not reach convergence, but terminated after 356 steps due to grid quality issues, which are discussed below.At this step, the convergence criterion was about two orders of magnitude above the threshold.However, from a practical point of view the convergence criteria employed in this study might appear rather strict and practical applications would also reach sufficient optima for   Improvements seen for the objective function are attributed to the more extreme deformations obtained from the p-Laplacian problem (2) with p > 2. In order to judge the final shape, the opening angle at the upstream tip may also serve as a measure.The interior opening angles listed in Tab. 2 decrease with greater values of p. Hence, increasing p clearly yields more pointy tips as also indicated by the comparison of tip shapes in Fig. 4. The Stokes flow problem investigated in [29] reported an opening angle of 120 • .The present results rapidly approach the reported opening angle from above.However, for p = 6 the opening angle falls below the reference value, which is attributed to the use of Navier-Stokes instead of the Stokes flow model [29].When attention is directed to the convergence speed, the half axis ratio at the 50th design step, i.e. a/b (it=50) , mentioned in Tab. 2 may be considered as a measure to assess the convergence speed.As all simulations start with a/b = 1, the tabulated data renders the influence of p-values on the ability of the p-Laplace approach to rapidly adjust the shape.It is also observed, that large deformations take place at an early stage of the optimization.For example, the final half axis ratio for p = 4 refers to approximately 2.7 and already 2/3 of this ratio are reached after 50 design steps.A closer inspection of Fig. 3 reveals, that the value of the objective function slightly increases after a quick descent for p = 5 (step 120-160) and p = 6 (step 90-140).This is possible due to the absence of a step size control and occurs because the multipliers that control the barycenter and the displacement strictly speaking only hold for a preceding iteration.Moreover, the p-Laplace problem (2) is not solved exactly in every iteration of algorithm 2 and preceding results are used as the initial guess for a subsequent optimization step.A similar phenomenon can be seen in results reported by Allaire et al. [2], who did use a similar augmented Lagrange procedure.The contours of the final shapes are depicted in Fig. 4, where optimal shapes for two consecutive values of p are compared with each other, i.e. for p = 2 with p = 3, p = 3 with p = 4 and so forth.The shape contours in Fig. 4a show significant differences, in particular at the tips of the resulting geometry.The computation with p = 2 clearly leads to a more rounded shape and a larger vertical extent than the other investigated p values.Shapes returned by 3 ≤ p ≤ 6 are displayed Fig. 4b -4d.Remarkably, a general difference between the respective contours is hard to identify for p ≥ 3. Thus, close-ups are used to assess the tip region.While a rounded tip region is still observed for p = 3, the tip becomes more pointy for p ≥ 4. Figures 4c and 4d only display small differences between the shapes obtained with p = 4, 5, 6.Thus, one can assume that the predicted optimal shapes tend to be converged for p ≥ 4.
4.1.2.Grid Deformation.Besides the convergence, the attainable optimal shapes and the reduction of the objective function value, another major aspect refers to the quality of the mesh updates.Maintaining the grid quality during the optimization process is crucial to the success of the CAD-free optimization procedure.In this study, we focus on (a) the grid orthogonality near the boundaries as well as (b) the cell aspect ratio in the vicinity of the walls.It is seen, that using the p-Laplacian problem (2) with p > 2 significantly improves the quality of the mesh updates in comparison with updates obtained from p = 2. Figs.5a -5e display the final grids in the upstream tip region for the five investigated values of p.A reasonable grid quality is generally observed for large values of p, even after substantial cumulative deformations due to several hundred optimization steps.On the contrary, the aspect ratio of the near wall cells in Fig. 5a increases significantly for p = 2.The cells become stretched and tend to buckle in normal direction, which hampers the iterative convergence of the primal and adjoint flow solver.Therefore the procedure terminated after 354 optimization steps for p = 2.In line with the change of the shape characteristics, a huge change of the mesh characteristics is observed when p is increased from p = 2 to p = 3, cf.Figs.5a and 5b.The grid is less compressed in the vicinity of the upstream tip for p = 3, despite the larger deviations from the initial grid indicated in Fig. 4a.As outlined by Figs.5c -5e, the grid does further improve when p is augmented to p = 4, 5, 6.A detailed comparison of the grids at the upstream tip of the final shape follows from Fig. 6 for p = 4 (green), p = 5 (blue) and for p = 6 (red).The post processed grids are aligned at the tip to support the comparison.The stretching of the grid increases with increasing p, which becomes obvious when observing the spacing along the horizontal center line.converge after 164 design steps and yields 6.4% drag reduction.Similar to the 2D case, the lower quality of the volume grid update restricted the amount design steps for p = 2. Figure 9 depicts the last shape obtained from computations with p = 2 which referred to 94 design steps and 4.93% drag reduction.Fig. 9a reveals that this shape is clearly characterized by round tips at the upstream and downstream ends.The optimal shape returned by p = 4 is shown in Fig. 10.As also indicated by the data listed in Tab. 3, pointy upstream and downstream ends are seen for p = 4 which also results in a significantly larger average half axis ratio.force is reduced by 32.6% when reaching the convergence criterion.Figure 13 compares the final grids in the vicinity of the upstream tip for p = 2 (left) and p = 4 (right).As indicated by Fig. 13a, the aspect ratio deteriorates for p = 2 since the near wall cells stretch in tangential direction.Moreover, cells (again) cluster at the tip.In contrast, the grid for p = 4 depicted in Fig. 13b features evenly distributed cells along the design surface.The predicted shape displays a pointy upstream tip which is a much better approximation to the solution to the optimization problem.A similar conclusion follows from Fig. 14  the situation at the downstream end.It is observed, that the pointy ends develop at later stages of the optimization process, particularly in the rear.Thus, the (almost) pointy rear is hardly reached for p = 2.In combination with p = 4, the cell distribution at the rear remains evenly distributed.Mind that a small round tip can still observed in Fig. 14b.To receive a sharp pointy tip at the downstream end, a step size control would be beneficial.

Conclusions
We presented a novel approach for shape optimization by approximating Lipschitz continuous transformations based on the relaxation of the definition of the steepest descent direction.Examples included were restricted to fluid dynamic applications, but also apply to other simulation areas.The main goal was to improve the shape optimization algorithm by considering a descent direction as the solution to the p-Laplace problem, and investigate the influence of increasing p on the convergence of the shape optimization procedure as well as the obtained shapes and the related updates of the volume grid.An important aspect refers to the behavior of the relative differences when p is increased, since this might guide towards sufficiently high p values associated with appreciated lower computational effort.
Results show that directions obtained from p-harmonic solutions improve the convergence with increasing p.At the same time, the optimal shapes improve regarding the value of the objective function.A remarkable feature is related to the ability of the p-Laplace approach to yield shapes with edges or pointy shapes, even when the initial shape does not contain such features.Furthermore, the quality of the computational grid is virtually preserved even when large deformations of the initial shape occur and no specific grid adjustment is required.Results of the present study suggest that p = 4 seems a sufficiently large p-value to gather the benefits of the p-harmonic approach.
Within future research different solution algorithms for the p-Laplace problem may be considered to improve the approximation of the steepest descent direction and thus reduce the computational efforts.Moreover, applications to large-scale 3D problems may be investigated.

4. 1 .
Drag Optimization in 2D Low Re Flow.The first case studies the drag minimization of a 2D circular cylinder exposed to low Reynolds number flow for p = 2, 3, . . ., 6.The computational domain is illustrated in Fig.1.The case is related to the setting initially described by Pironneau[29].Instead of the Stokes flow considered in[29], we employ a low Re Navier-Stokes formulation of the boundary value problem.The initial cylinder features a unit diameter (D = 1[m]) and is centered in a channel of length 50[m] and height 10[m].The fluid is characterized by a unit density and dynamic viscosity, i.e. ρ = 1[kg/m 3 ], µ = 1[P s • s].Dirichlet conditions are imposed at the inflow Γ in with v ∞ = (1, 0) T [m/s], which yields a unit Reynolds number Re D = 1.Slip wall boundary conditions are applied to the top and bottom boundaries of the channel.Outflow boundary conditions (4g) are used along Γ out .

Figure 2 .
Figure 2. Initial mesh of the 2D low Reynolds number problem.

Figure 4 .
Figure 4. Influence of the p-value on the predicted optimal (final) shapes.

(a) p = 2 6 Figure 5 .
Figure 5. Influence of the p-value on the final grids in the vicinity of the upstream tip.

Figure 7
depicts contour plots of the cell aspect ratio for the final grids obtained with p = 2 (left) and p = 3 (right).These p-values are particularly illustrative, since

( a )
Side view.(b) Perspective view onto the upstream tip.

Figure 9 .
Figure 9. Drag optimization of a sphere at Re = 1; final shape obtained from p = 2 after 94 steps.

Figure 10 .
Figure 10.Drag optimization of a sphere at Re = 1; optimal shape obtained from p = 4 after 164 steps.

4. 3 .Figure 11 .Figure 12 .Table 4 .
Figure 11.Initial mesh for the drag optimization of a 2D ellipses exposed to horizontal approach flow at Re = 3 • 10 6 .the design surface to ensure an adequate resolution of the high curvature region and is depicted in Fig. 11.The design surface Γ is discretised by approximately 2300 cells of equal size and the volume grid features 46k cells.The final shape

4 Figure 13 .
Figure 13.Upstream tip of final shapes for the 2D high Reynolds case.

4 Figure 14 .
Figure 14.Downstream tip of final shapes for the 2D high Reynolds case.

Table 1 .
Initial values for the parameters of the augmented Lagrange procedure.

Table 2 .
Performance indicators obtained with different p-values for 2D low Re case.

Table 2
also displays the final objective function values, which again reveal improvements for increasing p-values.Drag reductions refer to about 7, 6% for the Laplacian approach with p = 2 and increase to approximately 8, 1% for p = 6.

Table 3 .
Performance indicators obtained with different p-valuesfor the 3D low-Re case.
3] and the initial parameters of the augmented Lagrange algorithm 1 are again denoted in Tab. 1.The performance observed with the two investigated p-values is summarized in Tab. 3.For p = 4 the optimization did