Angelo Carnarius
email: angelo.carnarius@cfd.tu-berlin.de

Frank Thiele
email: frank.thiele@cfd-berlin.com

Emre Özkaya
email: ozkaya@mathcces.rwth-aachen.de

Anil Nemili
email: nemili@mathcces.rwth-aachen.de

Nicolas R Gauger
email: gauger@mathcces.rwth-aachen.de

Optimal Control of Unsteady Flows Using a Discrete and a Continuous Adjoint Approach

Keywords: optimal control, active flow control, discrete adjoint, continuous adjoint, unsteady turbulent flows, URANS

While active flow control is an established method for controlling flow separation on vehicles and airfoils, the design of the actuation is often done by trial and error. In this paper, the development of a discrete and a continuous adjoint flow solver for the optimal control of unsteady turbulent flows governed by the incompressible Reynoldsaveraged Navier-Stokes equations is presented. Both approaches are applied to testcases featuring active flow control of the blowing and suction type and are compared in terms of accuracy of the computed gradient.

Introduction

For many aerodynamic applications in aviation and automotive industry, flow separation has to be taken into account. The lift of an airfoil at a high angle of attack, for instance, decreases drastically, if the flow separates on the suction side.

Many studies in the past decades have shown that the aerodynamic behaviour of a body can be improved by using active flow control [START_REF] Gad-El-Hak | The Taming of the Shrew: Why Is It so Difficult to Control Turbulence[END_REF]. However, the choice of the control parameters is very case-specific and not trivial. An efficient method of finding the optimal set of actuation parameters is the gradient-based optimisation, which requires the calculation of the gradient of the cost function with respect to the control parameters. The control variables are then updated in an iterative manner according to a descent direction, which can be obtained from the gradient vector.

A very efficient way of computing the gradient is by using adjoint methods. Compared to simpler approaches such as Finite Differences or the Complex Taylor Series Expansion (CTSE) [START_REF] Squire | Using complex variables to estimate derivatives of real functions[END_REF], adjoint-based methods compute the gradient vector at a fixed expense independent of the number of actuation parameters. Adjoint methods are commonly divided into the continuous and the discrete approach.

In the continuous adjoint method [START_REF] Jameson | Aerodynamic Design via Control Theory[END_REF], first the optimality system for a given objective function is derived and the resulting PDEs are then discretised and solved numerically. This procedure is called first optimise then discretise. The continuous approach is numerically efficient but it is known to suffer from consistency problems. The gradient can become inaccurate for insufficient time steps and grid spacing, which can be disadvantageous for complex configurations. Furthermore, most statistical turbulence models required for the unsteady Reynoldsaveraged Navier-Stokes equations (URANS) are non-differentiable. The common approach is to use the so-called constant eddy viscosity or frozen turbulence assumption, i.e. the eddy viscosity is treated as independent of the control parameters and therefore taken from the primal solution. This assumption can lead to significant errors in the computed gradient [START_REF] Carnarius | Adjoint approaches for optimal flow control[END_REF].

The concept of the discrete adjoint method [START_REF] Elliot | Practical 3D aerodynamic design and optimization using unstructured meshes[END_REF][START_REF] Giles | Algorithm Developments for Discrete Adjoint Methods[END_REF][START_REF] Nielsen | Aerodynamic design optimization on unstructured meshes using the Navier-Stokes equations[END_REF] is to first discretise then optimise, i.e. the discretised governing equations are used to derive the optimality system. This approach allows the generation of a fully consistent optimality system independent of the grid size, time step and turbulence model, as it does not require analytical differentiability [START_REF] Nemili | Optimal Control of Unsteady Flows Using Discrete Adjoints[END_REF]. Furthermore, Automatic Differentiation (AD) techniques [START_REF] Griewank | Evaluating Derivatives: Principles and Techniques of Algorithmic Differentiation 2nd edition[END_REF] can be used to develop the discrete adjoint solver for a given simulation code in a semi-automatic fashion.

In this paper, we present the development of a continuous and a discrete adjoint solver for the optimal control of unsteady turbulent flows governed by the incompressible URANS equations. Both approaches, which are presented in more detail in sections 3 and 4, are based on the same flow solver ELAN [START_REF] Xue | Entwicklung eines effizienten parallelen Lösungsalgorithmus zur dreidimensionalen Simulation komplexer turbulenter Strömungen[END_REF]. For the current study, the adjoint solvers are applied to testcases which feature active flow control of the blowing and suction type.

Flow Model

Governing equations For this study, the unsteady, incompressible, turbulent flow in the domain Ω is described by the Reynolds-averaged Navier-Stokes equations 4∂u i ∂x i = 0 (1)

∂̺u i ∂t + ∂̺u i u j ∂x j + ∂p ∂x i - ∂ ∂x j (µ + µ t ) ∂u i ∂x j + ∂u j ∂x i = 0 , (2)
where u i and p are the Reynolds-averaged velocity and pressure, respectively. The density ̺ and the dynamic viscosity µ are constant for the cases shown here.

The eddy viscosity µ t is obtained from the Wilcox-k-ω-model [START_REF] Wilcox | Reassesment of the scale determing equation for advanced turbulence models[END_REF], which consists of transport equations for the turbulent kinetic energy k and the turbulent frequency ω.

Boundary Conditions At the farfield boundaries (Γ i ), u i , k and ω were prescribed. On the body surface (Γ b ), u i and ∂k/∂n were set to zero, whereas a high-Re boundary condition [START_REF] Xue | Entwicklung eines effizienten parallelen Lösungsalgorithmus zur dreidimensionalen Simulation komplexer turbulenter Strömungen[END_REF] was used for the turbulent frequency. At the outflow (Γ o ), the gradient of the turbulent quantities normal to the boundary was set to zero. For the Navier-Stokes equations, we want the sum of normal and friction forces to vanish at the outlet, i.e.

-pn i + (µ + µ t )

∂u i ∂x j + ∂u j ∂x i n j = 0 . (3)
At the control segment (Γ c ), the actuation velocity c i (b j ) as well as the turbulent quantities were prescribed. The dependence of c i (b j ) on the vector of the actuation parameters, b j , is case-specific and is given in the testcase descriptions.

Continuous Adjoint Approach

Let J be the objective function to be minimised. Then the optimisation problem can be stated as

J(u i , p, c i (b j )) min over (u i , p, c i (b j )) subject to R(u i , p, c i (b j )) = 0 , (4)
where R represents the state equations including the boundary conditions. In the cases presented here, the objective function can be written as5

J = - 1 T T 0 Γ b,c
(µ + µ t ) ∂u i ∂x j + ∂u j ∂x i n j -pn i -̺u i u j n j e i dA dt

+ γ u ∞ 1 T T 0 Γc ̺u 2 (u i n i ) 2 + ǫ dA dt . (5)
If the unity vector e i is parallel to the mean flow, eq. 5 is the time-averaged drag. If e i is oriented normal to the mean flow, eq. 5 represents the time-averaged downforce. The second integral is a penalty term which accounts for the energy consumption of the actuation and can be scaled by the factor γ. The parameter ǫ is only required for the differentiability of the penalty term, i.e. 1 ≫ ǫ > 0.

To solve the minimisation problem, one first introduces the Lagrange function

L(u i , p, c i (b j ), v i , q) = C J(u i , p, c i (b j ))+ T 0 Ω qR ̺ dV dt+ T 0 Ω v i R u dV dt , (6)
where the Lagrange multipliers v i and q are the adjoint velocity and pressure, respectively, and C is a scaling factor to fix the units. By setting the variation of L with respect to the state variables, ∂L ∂u k δu k and ∂L ∂p δp, to zero, one can obtain the adjoint equations. First the variations δu k and δp have to be separated from other terms by using integration by parts and the boundary conditions have to be applied to the boundary integrals. As the resulting equations have to be fulfilled for any variation δu k and δp, all integrals have to vanish individually, which gives the adjoint equations and boundary conditions. Setting the variation of L with respect to the control to zero and using the same procedure gives the equation for the gradient calculation. Due to the page limitation, a detailed derivation has to be omitted and the adjoint system can only be summarised. The adjoint PDEs read

∂v i ∂x i = 0 Ω - ∂̺v i ∂t + ̺v j ∂u j ∂x i - ∂̺u j v i ∂x j + ∂q ∂x i - ∂ ∂x j (µ + µ t ) ∂v i ∂x j + ∂v j ∂x i = 0 Ω v i = 0 Γ i ̺v i u j n j -qn i + (µ + µ t ) ∂v i ∂x j + ∂v j ∂x i n j = 0 Γ o v i + C T e i = 0 Γ b,c , (7)
with the initial condition v i = 0 at t = T . The gradient w.r.t. the actuation parameters can be evaluated from

dJ db n = T 0 Γc -̺u j v j n m -qn m + (µ + µ t ) ∂v m ∂x j + ∂v j ∂x m n j ∂c m ∂b n dA dt + γ u ∞ C T T 0 Γc ̺c k c k √ c i n i c j n j + ǫ c l n l n m + 2̺c m c i n i c j n j + ǫ ∂c m ∂b n dA dt .
(8) Note, that the frozen turbulence assumption has been used, i.e. an adjoint turbulence model is not required.

Discrete Adjoint Approach

If we consider the discrete implementations of the objective function J and the state equations R, the discrete optimisation problem can be stated as:

J d (y, b i ) min over (y, b i ) subject to R d (y, b i ) = 0 , (9)
where y = (u i , p) is the discrete state vector and J d , R d denote the discrete implementations of J and R. Note, that in the discrete realisation the actuation variables b i are chosen as independent variables. The gradient of J d with respect to the actuation parameters b i can be computed from

dJ d db i = ∂J d ∂b i -ψ ⊤ ∂R d ∂b i , (10)
where the adjoint vector ψ can be determined by solving the adjoint system

∂R d ∂y ⊤ ψ = ∂J d ∂y . ( 11
)
One way of constructing the adjoint system is by computing ∂R d /∂y and ∂J d /∂y using Finite Differences. The linear system of equations is then hand-coded and solved by an iterative method (e.g. GMRES). The resulting adjoint variables are used to calculate the gradient vector in eq. 10. A more promising way of developing the adjoint system is by employing the reverse mode of AD, which has the major advantage that it constructs the adjoint system consistently and computes the gradient vector dJ d /db i accurate to machine precision. In the present work, the discrete adjoint solver is developed by employing the AD tool TAPENADE [START_REF] Hascoët | TAPENADE 2.1 user's guide[END_REF] in reverse mode of differentiation.

If the reverse mode of AD is applied in a black-box fashion, the resulting adjoint code will have tremendous memory requirements. In order to reduce the excessive memory demands, we apply the reverse accumulation and checkpointing strategies for the Automatic Differentiation of the underlying flow solver. The solution strategy for the incompressible URANS equations mainly consists of two iterative loops: the time evolution step and the iterations for the velocitypressure coupling scheme. Inside each time step, the velocity-pressure coupling iterations are performed, which are commonly known as outer iterations in the CFD community. It may be noted, that the outer iterations for the velocitypressure coupling scheme converge to a fixed point in each time step. If AD is applied to these outer iterations in a black-box fashion, the flow solutions at each outer iteration of the primal solver must be saved for the adjoint part. However, the adjoint iterations require only the converged primal solution. Therefore, a lot of memory and run time can be saved, if we make use of the iterative structure and store only the converged flow solution in each physical time step. This can be achieved by employing the reverse accumulation approach [START_REF] Christianson | Reverse accumulation and attractive fixed points[END_REF][START_REF] Gauger | Automatic Differentiation of an Entire Design Chain for Aerodynamic Shape Optimization[END_REF], the details of which are presented in the following.

Consider the total derivative of a discrete objective function J d with respect to the control b i at the converged state solution y * for any time step:

dJ d (y * , b i ) db i = ∂J d (y * , b i ) ∂b i + ∂J d (y * , b i ) ∂y * dy * db i . ( 12
)
On the other hand, if we have a fixed point for the state solution

y * = G(y * , b i ) ⇔ R d (y * , b i ) = 0, we get dy * db i = ∂G(y * , b i ) ∂b i + ∂G(y * , b i ) ∂y * dy * db i = I - ∂G(y * , b i ) ∂y * -1 ∂G(y * , b i ) ∂b i . (13)
Multiplying on both sides with ∂J d (y * ,bi) ∂y *

⊤

, we obtain

∂J d (y * , b i ) ∂y * ⊤ dy * db i = ∂J d (y * , b i ) ∂y * ⊤ I - ∂G(y * , b i ) ∂y * -1 :=y * ⊤ ∂G(y * , b i ) ∂b i . (14)
From the definition of y * ⊤ in equation ( 14) and making use of equation ( 13), the adjoint fixed point iteration can be written as

y * ⊤ = y * ⊤ ∂G(y * , b i ) ∂y * + ∂J d (y * , b i ) ∂y * ⊤ . ( 15
)
The first term on the right hand side of the above equation is the adjoint of a single outer iteration. This can be generated by applying the reverse mode of AD to the wrapper subroutine G, which combines all the steps done within one outer iteration of the flow solver. The gradient vectors ∂J d /∂y * and ∂J d /∂b i come from the adjoint of the post-processor, which is computed only once for each time iteration.

We now focus our attention on adjoining the time iterations. In general, the computation of the unsteady adjoint solution over the time interval [0, T ] with N time steps requires the storage of flow solutions at time steps T 0 to T N -1 . The stored solutions are then used in solving the adjoint equations from T N to T 0 . For many practical aerodynamic configurations with millions of grid points and a large number of unsteady time steps, the storage costs may become prohibitively expensive.

One way of circumventing the excessive storage cost is by employing a checkpointing strategy [START_REF] Griewank | Evaluating Derivatives: Principles and Techniques of Algorithmic Differentiation 2nd edition[END_REF], where the flow solutions are stored only at selective time steps known as checkpoints. These are then used to recompute the intermediate states that have not been stored. In the present example, we chose r (r ≪ N ) checkpoints. We then have 0

= T 0 = T C1 < T C2 < • • • < T Cr-1 < T Cr < T N = T .
Here, T Cr represents the time step at r th checkpoint. During the adjoint computation over the subinterval [T Cr , T N ], required flow solutions at intermediate time steps are recomputed by using the stored solution at T Cr as the initial condition. The above procedure is then repeated over other subintervals T Cr-1 , T Cr until all adjoints are computed. It may be noted, that the checkpoints can be reused when they become free. We designate them as intermediate checkpoints.

Various checkpointing strategies have been proposed based on the storage criteria. If all the checkpoints are stored in main memory, it is called singlestage checkpointing. In yet another approach called multi-stage checkpointing [START_REF] Stumm | Multistage approaches for optimal offline checkpointing[END_REF], the checkpoints are stored both in main memory and on hard-disk, thus reducing the number of flow recomputations. In the present work, we have used the single-stage binomial checkpointing strategy, which is implemented in the algorithm revolve [START_REF] Griewank | Algorithm 799:Revolve: An implementation of checkpointing for the reverse or adjoint mode of computational differentiation[END_REF] and generates the checkpointing schedules in a binomial fashion, so that the number of flow recomputations is proven to be optimal.

Numerical Results

Cylinder with Pulsed Blowing and Suction

The first application is the unsteady laminar flow around a circular cylinder at a Reynolds-number of Re = 100, based on the cylinder diameter D and the freestream velocity u ∞ . The objective is to reduce the drag by applying pulsed blowing or suction according to

c n = u a sin [2πf (t -t 0 )] -u a ( 16
)
on 15 slits, which are equidistantly distributed in 75% of the cylinder surface, see fig. 1(a). In eq. 16, c n is the actuation velocity normal to the slit surface and u a , f and t 0 are the amplitude, frequency and phase shift, respectively. The actuation mode, i.e. blowing or suction, is set by the sign of the amplitude. For the case studied here, the actuation amplitudes at all slits are the parameters to be optimised, while the frequency and phase shift were fixed to f = 1 u ∞ /D and t 0 = 0 D/u ∞ , respectively. Only the continuous adjoint flow solver was applied to this testcase in order to test its accuracy in the unsteady laminar mode, i.e. without the influence of the frozen turbulence assumption. A numerical mesh consisting of about 25000 control volumes (CV) and a time step of ∆t = 0.04 D/u ∞ was used for the computations. In every iteration of the optimisation, which was performed with the steepest descent method, the primal solution was integrated over 15000 time steps. For the calculation of the objective function and the gradient, the first 5000 time steps were neglected to remove the initial transient. The optimisation was terminated when all sensitivities had dropped by two orders of magnitude.

As can be seen from fig. 2(a), the drag coefficient of the cylinder decreases from c d = 1.336 to c d = 0.899 when actuated with the optimal control parameters, which is a reduction of more than 30%. The comparison of the sensitivities at the first optimisation step shows a good agreement of the adjoint-based gradient with Finite Differences, see fig. 2(b). There are only small deviations, which can be attributed to the insufficient grid spacing and time step. Note, that only the slits on the upper half of the cylinder are shown, as the optimisation leads to a symmetric actuation. The same holds for the optimal amplitude distribution, which is presented in fig. 2(c). The blowing and suction at slits one to four generates a symmetric vortex pair which is pushed away from the rear of the cylinder by the blowing at slits five to eight. As is obvious in fig. 1(b), this increases the pressure level behind the cylinder, thus reducing the pressure drag.

NACA0015 Airfoil with Synthetic Jet Actuation

The second application is the lift maximisation for the unsteady turbulent flow around a NACA0015 airfoil at Re = 10 6 , based on the cord length c and the freestream velocity u ∞ . The angle of attack (AoA) is α = 20 o , leading to a massive separation on the suction side. In this case, sinusoidal blowing and suction (also called synthetic jet), which can be modelled according to

c i = u a r i sin [2πf (t -t 0 )] , r i = cos (β -θ) sin (β -θ) , (17)
is applied at four slits on the suction side of the airfoil with a constant frequency of f = 1.28 u ∞ /c. Compared to the pulsed actuation (eq. 16) the blowing angle β can now be varied. The angle of the slit surface, θ, is fixed by the geometry of the airfoil. Computations with the discrete and the continuous adjoint solver were performed on a coarse mesh with 9500 CV and ∆t = 0.005 c/u ∞ over 100 time steps, including the initial transient. The comparison of the sensitivity gradients, summarised in tab. 1, reveals an excellent agreement between the forward and reverse mode AD, giving only very small differences of approx. 1 × 10 -5 . Compared to the AD-based solver, the results of the continuous adjoint code are significantly less accurate. One reason for this is the insufficient grid spacing, which is known to cause consistency problems with the continuous adjoint approach [START_REF] Nemili | Optimal Control of Unsteady Flows Using Discrete Adjoints[END_REF]. Furthermore, this can also be attributed to the frozen turbulence assumption. The active flow control modifies the separation on the suction side of the airfoil considerably, which has a 1. Comparison of the sensitivities for the NACA0015 testcase strong impact on the turbulence field. This is completely neglected by the frozen turbulence assumption.

Summary and Outlook

In this paper, the development of a continuous and a discrete adjoint flow solver for the optimal control of unsteady, turbulent flows governed by the incompressible URANS equations was presented. For the continuous adjoint approach, the wide-spread frozen turbulence assumption was used, while the AD-based discrete approach is fully consistent independent of the grid size, time step and turbulence model, as it does not require analytical differentiability. The numerical efficiency of the discrete solver has been improved by employing the reverse accumulation technique and the binomial checkpointing, which allows the application of the discrete adjoint solver to practical configurations.

The numerical results of the drag reduction of the cylinder flow showed, that the continuous adjoint method works well for unsteady laminar flows. However, it gives fairly inaccurate sensitivity gradients when applied to the turbulent flow around a NACA0015 airfoil at a high Re-number due to the frozen turbulence assumption and insufficient grid spacing. In contrast to this, the sensitivities obtained from the AD-based adjoint solver are of excellent accuracy and match the forward mode AD nearly perfectly.

In future studies, the different approaches will be applied to more complex geometries such as multi-element high-lift configurations or simplified car models, aiming at a more detailed comparison of the adjoint methods in terms of accuracy and numerical efficiency.

Fig. 1 .

 1 Fig. 1. Contour plots for the cylinder flow

Fig. 2 .

 2 Fig. 2. Optimisation results for the cylinder flow

In the following, the Einstein summation convention is used, which implies summation from 1 to 3 over indices which appear twice in a single term. Indices, which appear only once take the value 1, 2 and 3 individually.

Note, that the negative sign of the first integral is a result of the convention that the normal vector ni is directed out of the wall-adjacent control volume, i.e. into the body surface.

Acknowledgments This research was partly funded by the German Science Foundation (DFG) within the scope of the project Instationäre Optimale Strömungskontrolle aerodynamischer Konfigurationen.


Marius Cornel
email: cornel.murea@uha.fr

Murea

Andrei Halanay
email: halanay@mathem.pub.ro

Embedding domain technique for a fluid-structure interaction problem

Keywords:

We present a weak formulation for a steady fluid-structure interaction problem using an embedding domain technique with penalization. Except of the penalizing term, the coefficients of the fluid problem are constant and independent of the deformation of the structure, which represents an advantage of this approach. A second advantage of this model is the fact that the continuity of the stress at the fluid-structure interface does not appear explicitly. Numerical results are presented.

A steady fluid-structure interaction problem

The present paper is devoted to the study of the numerical behavior of an elastic structure immersed in a viscous incompressible fluid. We use Stokes equation to model the flow motion. The displacement of the structure under the flow motion will be modeled by linear elasticity equations, under the small deformations assumption. In this paper, we study the steady case.

Let D ⊂ R 2 be a bounded open domain with boundary ∂D. Let Ω S 0 be the undeformed structure domain, and suppose that its boundary admits the decomposition ∂Ω S 0 = Γ D ∪ Γ 0 , where Γ 0 is a relatively open subset of the boundary. On Γ D we impose zero displacement for the structure. We assume that Ω S 0 ⊂ D. Suppose that the structure is elastic and denote by u = (u 1 , u 2 ) : Ω S 0 → R 2 its displacement. A particle of the structure with initial position at the point X will occupy the position x = ϕ (X) = X + u (X) in the deformed domain Ω S u = ϕ Ω S 0 . We assume that Ω S u ⊂ D and the fluid occupies Ω F u = D \ Ω S u . We set Γ u = ϕ (Γ 0 ) and we suppose that Γ u does not touch the container wall, i.e. ∂D ∩ Γ u = ∅. We recall that Γ 0 is a relatively open subset. The boundary Γ u represents the moving fluid-structure interface. The boundary of the deformed structure is ∂Ω S u = Γ D ∪ Γ u . In the case when Ω S u ⊂ D, the fluid-structure geometrical configuration is represented in Figure 1 and the boundary of the fluid domain admits the decomposition

∂Ω F u = ∂D ∪ Γ D ∪ Γ u . Γ Γ 0 u Γ D Ω F Ω S u u D Fig. 1. Geometrical configuration.
Also, for the first numerical test, we consider the case when Γ D ⊂ ∂D as in Figure 2 and the boundary of the fluid domain admits the decomposition

∂Ω F u = (∂D \ Γ D ) ∪ Γ u .

Weak formulation using an embedding domain technique with penalization

We introduce the tensor ǫ (w) = 1 2 ∇w + (∇w)

T and we assume that the fluid is Newtonian and the Cauchy stress tensor is given by σ F (v, p) = -p I + 2µ F ǫ (v), where µ F > 0 is the viscosity of the fluid and I is the unit matrix. We assume that the structure verifies the linear elasticity equation, under the assumption of small deformations. The stress tensor of the structure written in the Lagrangian framework is σ S (u) = λ S (∇ • u) I + 2µ S ǫ (u), where λ S , µ S > 0 are the Lamé coefficients. We present in an informal and intuitive manner the ideas behind our approximation approach using embedding domain technique with penalization. In the fluid domain, Stokes equations are solved:

-∇ • σ F (v, p) = f F , in Ω F u (1) ∇ • v = 0, in Ω F u (2)
We introduce two more equations concerning the fluid fields, but written on the deformed structure domain:

-∇ • σ F (v, p) + 1 ε P (v) = f F , in Ω S u (3) ∇ • v = 0, in Ω S u (4)
where ε > 0 is a penalization parameter,

P (v) = |v 1 | α-1 sgn (v 1 ) , |v 2 | α-1 sgn (v 2 ) (5)
where v = (v 1 , v 2 ) and 1 < α < 2 is a real number.

Remark 1. This choice of the penalization term is justified in [START_REF] Halanay | Existence and approximation for a steady fluidstructure interaction problem using fictitious domain approach with penalization[END_REF], in order to obtain existence of the fluid-structure interaction problem. For the steady case, the role of the penalization term is to obtain very small values of the fluid velocity in the structure domain. If we take other values for α or for a H 1 penalization, we do not get more regularity of the solution.

Let χ S u be the characteristic function of Ω S u . Combining (1) and (3), it follows that

-∇ • σ F (v, p) + 1 ε χ S u P (v) = f F , in D. (6)
Similarly, we have from ( 2) and ( 4)

∇ • v = 0, in D. (7)
In view of the equation (3), the "fictitious" fluid velocity and pressure defined on the structure domain Ω S u depend on ǫ. In the following, we denote by v ε and p ε the fluid velocity and pressure defined all over the domain D.

Let us introduce the bi-linear forms

a S u, w S = Ω S 0 λ S (∇ • u) ∇ • w S + 2µ S ǫ (u) : ǫ w S dX a F (v, w) = D 2µ F ǫ (v) : ǫ (w) dx b F (w, p) = - D (∇ • w) p dx
and the Hilbert spaces

W S = w S ∈ H 1 Ω S 0 2 ; w S = 0 on Γ D , W = H 1 0 (D) 2 , Q = L 2 0 (D) = {q ∈ L 2 (D) ; D q dx = 0}.
We assume for the moment that

f F ∈ L 2 (D) 2 , f S ∈ L 2 (Ω S 0 ) 2 and g ∈ H 1/2 (∂D) 2 , such that ∂D g • n F ds = 0. For a given u ∈ W 1,∞ (Ω S 0 ) 2 , such that u 1,∞,Ω S 0 < 1 and u = 0 on Γ D , we define: -fluid velocity v ε ∈ H 1 (D) 2 , v ε = g on ∂D, -fluid pressure p ε ∈ Q, -structure displacement u ε ∈ W S ,
as the solution of the following weakly coupled system of PDE's:

a F (v ε , w) + b F (w, p ε ) + 1 ε D H u P (v ε ) • w dx = D f F • wdx, ∀w ∈ W (8) b F (v ε , q) = 0, ∀q ∈ Q (9) a S u ε , w S = Ω S 0 f S • w S dX + Ω S 0 J σ F (v ε , p ε ) • ϕ F -T : ∇ X w S dX + 1 ε Ω S 0 J H u P (v ε • ϕ) • w S dX - Ω S 0 J f F • ϕ • w S dX, ∀w S ∈ W S (10) where ϕ(X) = X + u(X), F(X) = I + ∇ X u(X), J(X) = det F(X).
The equations ( 8) and ( 9) are obtained from ( 6) and ( 7). The coefficient H u in (8) is a regularization of the characteristic function of Ω S u , which is necessary in order to prove the continuity of the solution with respect to the structure displacement.

Remark 2. From the structure equation -∇ • σ S (u ε ) = f S , in Ω S 0 using Green's formula, we obtain for all w S = 0 on Γ D that

a S u ε , w S = Ω S 0 f S • w S dX + Γ0 σ S (u ε ) n S • w S dS.
We can prove (see [START_REF] Halanay | Existence and approximation for a steady fluidstructure interaction problem using fictitious domain approach with penalization[END_REF]) that the sum of the last three terms in (10) is equal to the fluid forces acting on the structure which is also equals to Γ0 σ S (u ε ) n S • w S dS. In fact, from (10) and the above weak formulation of the structure, we can get that the boundary condition at the interface concerning the continuity of the stress is verified in a weak sense (see [START_REF] Halanay | Existence and approximation for a steady fluidstructure interaction problem using fictitious domain approach with penalization[END_REF]). The second boundary condition at the interface is the continuity of the velocity, i.e. v = 0 on Γ u in the steady case. This is obtained by using the penalization term in the structure domain.

Define the nonlinear operator

T ε : u ∈ W 1,∞ (Ω S 0 ) 2 ; u 1,∞,Ω S 0 < 1, u = 0 on Γ D → W 1,∞ (Ω S 0 ) 2 by T ε (u) = u ε .
A solution of the penalized fluid-structure interaction problem will be, by definition, a fixed point of T ε . In [START_REF] Halanay | Existence and approximation for a steady fluidstructure interaction problem using fictitious domain approach with penalization[END_REF], we discuss the existence of a solution of the penalized fluid-structure interaction problem. The convergence of u ε , v ε , p ε when ε goes to 0 is also analyzed.

Partitioned procedures based on fixed point iterations

The penalized term P (v) is non-linear in v for α = 2. But for α = 2, we have P (v) = v. Now, the fluid problem, at the Step 2 of the algorithm below, becomes linear and, for a given u k ǫ , it has a unique solution. For 1 < α < 2, we can prove the existence of a fixed point for the nonlinear operator T ε defined at the end of the previous section, but not for α = 2. We can also replace H u in (8) by χ S u the characteristic function of Ω S u in order to simplify the computation. The regularization of the the characteristic function was necessary in order to prove the continuity of the solution with respect to the structure displacement.

Under the assumption of small displacements for the structure, we can approach the Jacobian determinant J by 1 and the gradient of the deformation F by the identity matrix I. The structure problem at the Step 3 is linear and, for given v k ǫ and p k ǫ , it has a unique solution. Algorithm

Step 1. Given the initial displacement of the structure u 0 ∈ W S , compute the characteristic function χ S u 0 , put k := 0.

Step 2. Find the velocity v ε ∈ H 1 (D) 2 , v ε = g on ∂D and the pressure

p k ε ∈ Q by solving the fluid problem a F v k ε , w + b F w, p k ε + 1 ε D χ S u k ε v k ε • w dx = D f F • w dx, ∀w ∈ W b F v k ε , q = 0, ∀q ∈ Q. Step 3. Find the new displacement of the structure u k+1 ε ∈ W S by solving a S u k+1 ε , w S = Ω S 0 f S -f F • w S dx + Ω S 0 2µ F ǫ v k ε : ǫ w S dx - Ω S 0 ∇ • w S p k ε dx + 1 ε Ω S 0 v k ε • ϕ k ε • w S dx ∀w S ∈ W S where ϕ k ε (X) = X + u k ε (X). Step 4. Stopping test: if u k ε -u k+1 ε 0,Ω S 0 ≤ tol, then Stop.
Step 5. Compute the characteristic function χ S u k+1 ε , put k := k + 1 and Go to Step 2.

It is possible to consider Navier-Stokes equations for the fluid domain, then at the Step 2 we have to solve a non-linear system. In this case we can use the non-linear penalization term for α = 2. If we use a non-linear model for the structure, too, we have to use more accurate approximations for the the Jacobian determinant J and the gradient of the deformation F.

Numerical tests

The numerical tests have been produced using the software FreeFem++ [3].

Test 1. Shell in steady-state cross flow

First, we have performed numerical simulation using a 2D model adapted from [START_REF] Bathe | Benchmark problems for incompressible fluid flows with structural interactions[END_REF] (see Figure 2) where we have changed the physical parameters of the fluid and of the structure. The dimensions of a rectangular elastic structure are: height ℓ = 3 m, thickness h = 0.125 m. The computational of the fluid D is a rectangle of height H = 5 m and length L = 12 m. The distance between the left side of the fluid and the left side of the structure is 2 m. The lower left corner of Figure 2 is (x 1 = 0, x 2 = 0).

We denote by Σ 1 , Σ 3 the left and the right vertical boundaries of D and by Σ 2 , Σ 4 the bottom and the top boundaries of D, respectively.

The mechanical proprieties of the structure (polybutadiene) are: Young modulus E S = 1.6 × 10 6 N/m 2 , Poisson's ratio ν S = 0.49, the applied volume forces on the structure f S : Ω S 0 → R 2 , f S = (0, 0) N/m 3 . The dynamic viscosity of the fluid (glycerin) is µ F = 1.14 N • s/m 2 . The inflow velocity profile on Σ 1 is

v 1 (x 1 , x 2 ) = V × 1.5 2H x 2 -x 2 2 H 2 m/s, V = 1, v 2 (x 1 , x 2 ) = 0.
The other boundary conditions are:

v = 0 (no-slip) on Σ 2 , v • n F = 0 (slip) on Σ 4 and v × n F = 0, n F • σ F (v, p
) n F = 0 (the tangential velocity and the normal traction are zero) on Σ 3 . We use a fixed mesh for the fluid domain of 13096 triangles and 6719 vertices. The mesh of the structure domain has 188 triangles and 145 vertices. The fluid and structure meshes are not compatible, for example, a vertex on the structure boundary is not necessary a vertex on the fluid mesh. For the approximation of the fluid velocity and pressure we have employed the triangular finite elements P 1 +bubble and P 1 respectively, also called "mini" finite elements. The finite element P 1 was used in order to solve the structure problem. The characteristic function was approached by P 0 finite elements. We have performed the simulation using the Algorithm described in the previous section. For the stopping criterion at the Step 4, we have used the tolerance tol = 0.2 × 10 -4 . The penalization parameter is ε = 10 -3 . The stopping criterion holds after 5 iterations of the fixed point algorithm. The maximal structural displacement is 0.14 m.

The fluid velocity into the fictitious domain is very small

v ε 0,Ω S uε = D χ S uε v ε • v ε dx = 0.0703384.

Test 2. Flexible appendix in a flow

We have adapted the benchmark from [START_REF] Turek | Proposal for numerical benchmarking of fluid-structure interaction between an elastic object and laminar incompressible flow[END_REF]. Originally, the structure was placed horizontally, parallel to the flow, but the displacements in this case are very small. We have placed the structure vertically, transversely to the flow, see Figure 5. The structure is composed by a rectangular flexible appendix attached to a fixed circle. The circle center is positioned at (0.2, 0.2) m measured from the left top corner of the channel. The circle has the radius r = 0.5 m and the rectangular appendix is of length ℓ = 0.35 m, thickness h = 0.02 m. The Young modulus is E S = 1.6 × 10 6 N/m 2 and Poisson's ratio is ν S = 0.49 (polybutadiene).

The channel has the length L = 2.5 m and the width H = 0.75 m. The fluid dynamic viscosity is µ F = 1.420 N • s/m 2 (glycerin).

We have used the following boundary conditions: at the inflow the velocity is

v 1 (x 1 , x 2 ) = V × 1.5 H x 2 -x 2 2 (H/2) 2 m/s, V = 1, v 2 (x 1 , x 2 ) = 0;
at the bottom and the top we have imposed the no-slip boundary condition v = 0 and at the outflow the traction free σ F (v, p) n F = 0. We use a fixed mesh for the fluid domain of 30330 triangles and 15461 vertices and a structure mesh of 128 triangles and 97 vertices. We have employed the same finite elements as for the Test 1. We have treated by the embedding domain technique only the flexible part of the structure.

The penalization parameter is ε = 10 -4 and for the stopping criterion tol = 10 -8 . The fixed point algorithm stops after 8 iterations.

The maximal horizontal displacement of the structure is 0.10886 m. The pressure and the velocity of the fluid are presented in Figures 6 and7

Conclusions

We have presented a fixed point algorithm for solving steady fluid-structure interaction problem. Using the embedding domain technique with penalization, the fluid equations as well as the structure equations are solved in fixed meshes. The fluid and structure meshes could be generated independently. The algorithm can be used for the three dimensions problems.

Fig. 2 .

 2 Fig. 2. Geometrical configuration for the Test 1.

Fig. 3 .

 3 Fig. 3. The fluid pressure [Pa].

Fig. 4 .

 4 Fig. 4. The fluid velocity around the final position of the structure. In each point of the grid, there is an arrow giving the direction of the velocity and the length of the arrow is proportional to the euclidean norm of the velocity. The maximal value for the horizontal component v1 is 2.82 m/s and for the vertical component v2 is 1.18 m/s

Fig. 5 .

 5 Fig. 5. Geometrical configuration.

  . The fluid velocity in the fictitious domain is very smallv ε 0,Ω S uε = D χ S uε v ε • v ε dx = 0.080920.

Fig. 6 .

 6 Fig. 6. The fluid pressure [Pa].

Fig. 7 .

 7 Fig. 7. The fluid velocity [m/s] around the final position of the structure. In each point of the grid, there is an arrow giving the direction of the velocity. The length of the arrow is proportional to the euclidean norm of the velocity which is represented in the color bar.

Acknowledgment

The second author gratefully acknowledges support by Grant CNCS Romania 145/2011.


Andrei Halanay
email: halanay@mathem.pub.ro

Cornel Marius Murea
email: cornel.murea@uha.fr

Fixed domain algorithms in shape optimization for stationary Navier-Stokes equations

Keywords: shape optimization, optimal control, penalization, approximate extension, gradient method

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

Introduction

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

The approximating extension technique makes possible the approximation of the solution to the shape optimization problem by a solution of an optimal control problem. The basic reference will be [START_REF] Halanay | Shape optimization for stationary Navier-Stokes equations[END_REF]. For shape optimization, the general references are [START_REF] Pironneau | Optimal shape design for elliptic systems[END_REF], [START_REF] Delfour | Shapes and Geometrics, Analysis, Differential Calculus and Optimization[END_REF] and for optimal control [START_REF] Lions | Optimal control of systems governed by partial differential equations[END_REF], [START_REF] Neittaanmaki | Optimal control of nonlinear parabolic systems, Theory, algorithms and applications[END_REF]. In particular, for shape optimization for fluids a standard work is [START_REF] Mohammadi | Applied Shape Optimization for Fluids[END_REF]. In optimal design related to optimal control, relevant contributions to the topic of this paper are [START_REF] Borrvall | Topology optimization of fluids in Stokes flow[END_REF], [START_REF] Gao | Optimal shape design for viscous incompressible flow[END_REF], [START_REF] Roubicek | Lipschitz stability of optimal control for steady-state Navier-Stokes equations[END_REF].

In Section 2, the shape optimization problem for steady Navier-Stokes is presented. The directional derivative of the cost function is given in Section 3. A gradient type algorithm is also introduced. In Section 4, numerical results are presented in order to design a nozzle.

Formulation of the shape optimization problem and approximating extensions

Let d be a natural number, d ≤ 4, let D ⊂ R d be a bounded fixed domain and suppose a family O of admissible subdomains Ω ⊂ D is given, satisfying an uniform Lipschitz condition on the boundary ∂Ω.

With the standard notations from [START_REF] Temam | Navier-Stokes equations. Theory and numerical analysis[END_REF],

V(Ω) = {y ∈ D(Ω) d |div y = 0}, V (Ω) = closure of V(Ω) in H 1 0 (Ω) d = {y ∈ H 1 0 (Ω) d |div y = 0}
(since ∂Ω is lipschitzian), we have the following weak formulation of the stationary Navier-Stokes equation with Dirichlet boundary (non-slip) conditions:

Ω   ν d i,j=1 ∂y j ∂x i ∂v j ∂x i + d i,j=1 y i ∂y j ∂x i v j   dx = Ω   d j=1 f j v j   dx, ∀v ∈ V (Ω) (2.1)
or (see [START_REF] Temam | Navier-Stokes equations. Theory and numerical analysis[END_REF]), ν((y, v)

) Ω + b Ω (y, y, v) = Ω f • v dx. Here f = (f 1 , . . . , f d ) ∈ H -1 (D) d ,
E ⊂ Ω is a fixed set and y 0 ∈ L 2 (E) d is given.
This functional is a particular case of a larger class,

J(Ω) = ∧ j(x, y(x), ∇y(x))dx, ∧ = E or ∧ = Ω,
that are studied in [START_REF] Halanay | Shape optimization for stationary Navier-Stokes equations[END_REF].

The uniform Lipschitz assumption turns O into a compact with respect to the Hausdorff-Pompeiu complementary metric (see [START_REF] Neittaanmaki | Optimization of elliptic systems. Theory and applications[END_REF], p. 466). Based on this, it is inferred in [START_REF] Halanay | Shape optimization for stationary Navier-Stokes equations[END_REF] that if there exists an admissible Ω and a corresponding solution of (2.1) for which J( Ω) is finite then there exists at least an optimal pair [Ω * , y * ] ∈ O × V (Ω * ). So, the optimization problem is well defined but its solution is generally nonunique.

If

X(D) ⊂ C( D) is a functional space, define, for g ∈ X(D), Ω = Ω g = int{x ∈ D|g(x) ≥ 0}. If E ⊂ Ω is to be fulfiled one must require g(x) ≥ 0 ∀x ∈ E.
g is called a parametrization of Ω g and Ω g an admissible domain. The solutions of (2.1) in Ω g will be denoted as y g .

If H : R → {0, 1} is the Heaviside function, H • g = χ Ωg , the characteristis function of Ωg . For ε > 0 the following smoothing of the Yosida approximation of the maximal monotone extension of H will be used:

H ε (r) =          1, r ≥ 0 (ε -2r)(r + ε) 2 ε 2 , -ε < r < 0 0, r ≤ -ε
(see also [START_REF] Makinen | On a fixed domain approach for shape optimization problem[END_REF], [START_REF] Neittaanmaki | Fixed domain approaches in shape optimization problems with Dirichlet boundary conditions[END_REF]). It is easy to see that H ε ∈ C 1 (R) and is lipschitzian. The boundary value problem (2.1) has an approximate extension

ν((y ε , v)) D + b D (y ε , y ε , v) + 1 ε D [1 -H ε (g)]y ε • v dx = D f • v dx (2.3) Suppose now that d = 3. It is proved in [6] that, for C 1 = 9m(D) 1/6 , if ν 2 > C 1 ||f || V * (2.4) then (2.
3) has an unique solution y ε (g) that depends continuously on g as a function from (C(D), || || ∞ ) to L 2 (D) 3 . The following theorem that is proved in [START_REF] Halanay | Shape optimization for stationary Navier-Stokes equations[END_REF], §3, allows the approximation of the shape optimization problem (2.1), (2.2) by the optimal control problem (2.2), (2.3). Theorem 2.1. If (2.4) holds then there exists a sequence ε n → 0 such that y εn (g)| Ωg → y g weakly in H 1 (Ω g ) 3 and strongly in L 2 (Ω g ) 3 .

The directional derivative and a gradient type algorithm

In order to solve the optimal control problem (2.2), (2.3) through a gradient type algorithm an important step is the calculation of the directional derivative of the mapping g → J[y ε (g)] in the direction w ∈ X(D). It is proved in [START_REF] Halanay | Shape optimization for stationary Navier-Stokes equations[END_REF], §4, that, under the uniqueness condition (2.4), this derivative in direction w,

∂y ε ∂w (g) = (z 1 , z 2 , z 3 ) ∈ V (D) is the solution of the equation in variations D   ν 3 i,j=1 ∂z j ∂x i ∂v j ∂x i + 3 i,j=1 y ε,i ∂z j ∂x i v j + 3 i,j=1 z i ∂y ε,j ∂x i v j   dx+ + 1 ε D [1 -H ε (g)]z • v dx = 1 ε D ((H ε ) ′ (g)w)y ε • r dx (3.1)
It is also proved in [START_REF] Halanay | Shape optimization for stationary Navier-Stokes equations[END_REF], §4, that under condition (2.4), equation (3.1) has an unique solution.

For the optimal control problem (2.2), (2.3) the co-state equation (see [START_REF] Dede | Optimal flow control for Navier-Stokes equations, drag minimization[END_REF], [START_REF] Gunzburger | Adjoint equation-based methods for control problems in viscous, incompressible flows[END_REF], [START_REF] Posta | Optimal control of Navier-Stokes equations by Oseen approximations[END_REF]) is

D   ν 3 i,j=1 ∂p ε,j ∂x i ∂v j ∂x i - 3 i,j=1 y ε,i ∂p ε,j ∂x i v j + 3 i,j=1 ∂y ε,j ∂x i p ε,j v j   dx+ + 1 ε D [1 -H ε (g)]p ε • v dx = E (y ε -y 0 ) • v dx. (3.2)
Under condition (2.4) the equation (3.2) has a unique solution p ε ∈ V (D). The algorithm will result from the following theorem Theorem 3.1 ([6], §4, Th.5). The direction derivative in g ∈ X(D) of J[y ε (g)] in the direction w ∈ X(D) is given by

∂J ∂w [y ε (g)]w = 1 ε D ((H ε ) ′ (g)w)y ε • p ε dx (3.3)
(p ε is the unique solution of (3.2)).

Algorithm

Step 0 Choose a starting parametrization g 0 and a positive scalar ǫ. Set k = 0.

Step 1 Find y ε the solution of (2.3).

Step 2 Find p ε the solution of (3.2).

Step 3 Set the descent direction

w k = -y ε • p ε . If w k < tol stop. Step 4 Determine g k+1 = g k + θ k w k , θ k > 0 by means of an approximate minimization J(g k+1 ) ≈ min θ≥0 J(g k + θw k ).
Step 5 Update k = k + 1 and go to the Step 1.

For the inaccurate line search at the Step 4, the methods of Goldstein and Armijo were used. If we denote by j : [0, ∞) → R the function j(θ) = J(g k +θw k ), we determine θ k > 0 such that

j(0) + (1 -λ) θ k j ′ (0) ≤ j(θ k ) ≤ j(0) + λθ k j ′ (0) (1)
where λ ∈ (0, 1/2).

Numerical results. Shape optimization of a nozzle

Problem setting

We have adapted the nozzle problem from [START_REF] Pironneau | Optimal shape design for elliptic systems[END_REF]. We assume that the flow in a nozzle is governed by the steady Navier-Stokes equation with prescribed traction at the inflow and outflow. The problem is to design a nozzle that gives a prescribed velocity near the exit. This kind of problem arises in rocket engine industries, in the design of a spray for applying a coating or in the manufacture of high-resolution inkjet printer.

We assume that the polyhedron [A

1 A 2 A 3 A 4 A 5 A 6 A 7 ]
is the fixed computational domain D. The coordinates of its vertices are: A 1 (H, 0), A 2 (0, 0), A 3 (L/2, 0), A 4 (L, Hh), A 5 (L + ℓ, Hh), A 6 (L + ℓ, H), A 7 (L, H), where L = 6, ℓ = 1, H = 3, h = 1. We denote by E the observation zone which is the rectangle

[A 4 A 5 A 6 A 7 ].
We denote by Σ in the boundary [A 1 A 2 ] representing the inflow section and by Σ out the boundary [A 5 A 6 ] representing the outflow section. The desired fluid velocity in the observation zone E is

y 0 = v out 4(H -x 2 )(x 2 -H + h) h 2 , 0 , where v out = 4.
The fluid viscosity is µ = 1 and its density is ρ = 1. Let ω ⊂ D such that ω ∩ E = ∅. We look for a connected domain Ω verifying D \ ω ⊂ Ω ⊂ D and minimizing the cost function

E A6 ω A4 A5 A7 A1 A2 A3
J = 1 2 E (y ǫ -y 0 ) • (y ǫ -y 0 ) dx.
The traction imposed on the inflow is (100, 0) and on the outflow it is (0, 0). We impose no-slip condition on the other boundaries, including the free boundary.

Descent direction

In Figure 2, we show 1-H ǫ (g) which is an approximation of the characteristic function of the domain D \ Ω, for a typical admissible parametrization g.

We remark that the (H ǫ ) ′ (r) vanishes on R, excepting for r ∈ (-ǫ, 0). Consequently, the zone in D, where (H ǫ ) ′ (g) = 0 is very narrow, see Figure 3. When ǫ is very small, this zone could be empty. For this season, we have taken as descent direction not -(H ǫ ) ′ (g)y ǫ • p ǫ which is given by (3.3), but w d = -y ǫ • p ǫ . We recall that (H ǫ ) ′ (g) ≥ 0, consequently w d = -y ǫ • p ǫ is a descent direction in view of (3.3). Table 1. Mesh parameters used in Figure 3 Mesh 1 Mesh 2 Mesh 3 Mesh 4

Fig. 3. The zone where the derivative of the Yosida approximation of the Heaviside function is not vanishes for ǫ = 10 -1 . The mesh parameters are prezented in Table 1 Numerical parameters The mesh of D has 15032 triangles and 7697 vertices. We have used the following finite elements: P 1 + bubble for the velocity, P 1 for the pressure and for the g.

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

We have tested our algorithm for three initial values of g: i) for g(x 1 , x 2 ) = 10 -4 (x 2 -

x 2 1 18 -1.5), the initial value of J is 0.518529 and its final value of is 0.00697268; ii) for g(x 1 , x 2 ) = 10 -4 (x 2 -

x 2 1 18 -0.5), the initial value of J is 0.185802 and its final value of is 0.00466894; iii) for g(x 1 , x 2 ) = 10 -4 (x 2 -

x 2 1 18 -1), the initial value of J is 0.0735282 and its final value is 0.00448260.

Numerical results

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

We remark in Figure 4 case ii), that the zero level set of the initial g partially coincides with the zero level set of the final g. In fact, the initial g vanishes on the free boundary. Since, we impose non-slip boundary condition for y ǫ and p ǫ on the free boundary, the descent direction w d = -y ǫ • p ǫ vanishes on the free boundary, also. Consequently, g k+1 could have the same zero level set as g k .

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

Fig. 1 .

 1 Fig. 1. Computing domain

Fig. 2 .

 2 Fig. 2. The value of 1 -H ǫ (g) on D for ǫ = 10 -4 .

Error(g) = ω ( 1 -

 ω1 H ǫ (g)) y ǫ • y ǫ dx = 0.000446 Aknowledgement This work was supported by ANCS Grant 480/17.03.2011 in the framework of the "Brâncuşi" programm of cooperation between Romania and France.

Fig. 4 .

 4 Fig. 4. Case i): final shape (bottom) obtained from the initial shape (top). Case ii): final shape (top) obtained from the initial shape (bottom). Case iii): final shape (bottom) obtained from the initial shape (top).

Fig. 5 .

 5 Fig.5. Fluid velocity for the optimal shape in the case iii)


Toni Lassila
email: toni.lassila@epfl.ch

Andrea Manzoni
email: andrea.manzoni@epfl.ch

Gianluigi Rozza
email: gianluigi.rozza@epfl.ch

Reduction strategies for shape dependent inverse problems in haemodynamics ⋆

Keywords: Model order reduction, reduced basis methods, free-form deformation, radial basis functions, computational fluid dynamics, shape parametrization, blood flows

This work deals with the development and application of reduction strategies for real-time and many query problems arising in fluid dynamics, such as shape optimization, shape registration (reconstruction), and shape parametrization. The proposed strategy is based on the coupling between reduced basis methods for the reduction of computational complexity and suitable shape parametrizations -such as free-form deformations or radial basis functions -for low-dimensional geometrical description. Our focus is on problems arising in haemodynamics: efficient shape parametrization of cardiovascular geometries (e.g. bypass grafts, carotid artery bifurcation, stenosed artery sections) for the rapid blood flow simulation -and related output evaluation -in domains of variable shape (e.g. vessels in presence of growing stenosis) provide an example of a class of problems which can be recast in the real-time or in the manyquery context.

Introduction and motivation

In last decades more and more powerful computers have allowed to solve numerical problems of very large dimensions and describing very complex phenomena. Nevertheless, a computational reduction is still crucial whenever interested to high performances in rapid -even real-time -simulations and/or repeated output evaluations -seen as many queries evaluations-for different values of some inputs of interest.

A general strategy for reduction in shape dependent flows

Flow control and optimization problems can be formulated as the minimization of a given cost functional (or output) controlling some input parameters which can be physical quantities (e.g. source terms or boundary values) or, alternatively, geometrical quantities; we refer to the latter case as flow control by shape variation, and the optimization of the corresponding flow geometries is thus one possibility to reach that goal; we refer to this case -the most difficult one among flow control problems -such as shape optimization or shape registration/reconstruction problems [START_REF] Haslinger | Introduction to shape optimization: theory, approximation, and computation[END_REF]. Concerning applications arising in fluid mechanics, cost functionals are expressed as functions of flow variables (such as velocity, pressure, temperature), while constraints are usually given in form of PDE systems (Stokes, Navier-Stokes equations, with or without coupling with a structural equation to account for fluid-structure interaction effects) describing the flow, besides topological constraints on the shape of the domain, if necessary. Since (i) optimization procedures require repetitive evaluations of outputs, (ii) PDEs can be hard to solve and (iii) discretization is expensive when geometry keeps changing, computational costs are usually very high; we thus want to address suitable strategies to reduce numerical efforts in many-query problems. Substantial computational saving becomes possible thanks to a reduced order model which relies on two reduction steps: (i) parameterization of the admissible shapes and (ii) substitution of the full-order finite element (FE) solution of flow problems with a reduced solution obtained by the reduced basis (RB) method [START_REF] Rozza | Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations[END_REF]. In fact, once an equivalent parametrized formulation of the flow problemnow embedding the shape as a parametric quantity -can be derived, reduced basis method for parametrized PDEs, enables to evaluate the output very rapidly. In the end, at the outer level a suitable iterative procedure for the optimization is performed. A brief presentation of the whole framework can be found in [START_REF] Manzoni | Model order reduction by reduced basis for optimal control and shape optimization[END_REF], while a more detailed analysis has been recently addressed in [START_REF] Lassila | A reduced computational and geometrical framework for inverse problems in haemodynamics[END_REF][START_REF] Lassila | Boundary control and shape optimization for the robust design of bypass anastomoses under uncertainty[END_REF].

Abstract setting

From an abstract point of view, a shape optimization/identification can be seen as an optimal control problem for which the control variable is the shape of the domain Ω itself. This entails the minimization of a cost functional J (•) over a set of admissible shapes O ad , by finding the optimal shape of the domain where the PDE is defined:

find Ω = arg min

Ω∈O ad J (Y (Ω)) (1)
where J (Y (Ω)) depends on the solution Y = Y (Ω) of a PDE state problemdefined on Ω -which can be written in an abstract form as

Y ∈ Y(Ω) : A(Y, W ; Ω) = F(W ; Ω), ∀ W ∈ Y(Ω). (2)
Here A(•, •; Ω) is a continuous, uniformly inf-sup stable bilinear form and F(•; Ω) is a bounded linear form, both defined on the original domain Ω; Y(Ω) denotes a suitable functional space defined over Ω. Let us assume that the shape Ω = Ω(µ) depends on a set of input parameters µ = (µ 1 , . . . , µ p ) ∈ D ⊂ R p ; in this way, problem (1)-( 2) can be reduced to the following parametric optimization inverse problem: find μ = arg min

µ∈D ad J (Y (µ)) (3)
where D ad ⊆ D and Y (µ) solves

Y (µ) ∈ Y(Ω(µ)) : A(Y (µ), W ; µ) = F(W ; µ), ∀ W ∈ Y(Ω(µ)). ( 4
)
For a more general setting and overview, see e.g. [START_REF] Lassila | A reduced computational and geometrical framework for inverse problems in haemodynamics[END_REF].

2 Reduced basis method for computational reduction Our approach to shape dependent flow problems takes advantage of reduced basis (RB) methods for rapid and reliable prediction of engineering outputs associated with parametric PDEs [START_REF] Rozza | Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations[END_REF][START_REF] Patera | Reduced Basis Approximation and A Posteriori Error Estimation for Parametrized Partial Differential Equations[END_REF][START_REF] Quarteroni | Certified reduced basis approximation for parametrized partial differential equations in industrial applications[END_REF]; see e.g. [START_REF] Rozza | Reduced basis methods for Stokes equations in domains with non-affine parameter dependence[END_REF][START_REF] Rozza | On the stability of the reduced basis method for Stokes equations in parametrized domains[END_REF][START_REF] Rozza | Reduced basis approximation and error bounds for Stokes flows in parametrized geometries: roles of the inf-sup stability constants[END_REF] for applications to the Stokes problem and [START_REF] Quarteroni | Numerical solution of parametrized Navier-Stokes equations by reduced basis methods[END_REF][START_REF] Veroy | Certified real-time solution of the parametrized steady incompressible Navier-Stokes equations: rigorous reduced-basis a posteriori error bounds[END_REF][START_REF] Deparis | Reduced basis method for multi-parameter-dependent steady Navier-Stokes equations: Applications to natural convection in a cavity[END_REF] for the Navier-Stokes case. The method is built upon a classical finite element (FE) "truth" approximation space Y N of (typically very large) dimension N and is based on the use of "snapshot" FE solutions of the PDEs, corresponding to certain parameter values, as global approximation basis functions previously computed and stored. The RB framework requires a reference (µ-independent) domain Ω in order to compare, and combine, FE solutions that would be otherwise computed on different domains and grids; moreover, this procedure enables to avoid shape deformation and remeshing that normally occur at each step of an iterative optimization procedure [START_REF] Rozza | Model order reduction by geometrical parametrization for shape optimization in computational fluid dynamics[END_REF]. In Sect. 3 two possible techniques for the construction of such a mapping will be briefly recalled. We thus consider Ω as reference domain related to the parameter-dependent "original" domain of interest Ω(µ) through a parametric mapping T (•; µ), s.t. Ω(µ) = T ( Ω; µ). By mapping the problem (3) back to the reference domain Ω, we obtain the following problem in its abstract form:

find μ = arg min µ∈D ad s(µ) = J (Y (µ)) s.t. Y (µ) ∈ Y( Ω) : Ã(Y (µ), W ; µ) = F(W ; µ), ∀ W ∈ Y( Ω).
(5) Focusing on shape optimization and/or registration problems, and following the so-called discretize than optimize approach, the standard Galerkin FE approximation of (5) reads as follows:

find μ = arg min µ∈D ad s N (µ) = J (Y N (µ)) s.t. Y N (µ) ∈ Y N : Ã(Y N (µ), W ; µ) = F(W ; µ), ∀ W ∈ Y N .
The reduced basis method provides an efficient way to compute an approximation Y N (µ) of Y N (µ) (and related output) by using a Galerkin projection on a reduced subspace made up of well-chosen FE solutions, corresponding to a specific choice S N = {µ1 , . . . , µ N } of parameter values. Denoting

Y N N = span{Y N (µ n ), n = 1, . . . , N }, (6)
the RB space, the RB formulation of ( 5) is as follows:

find μ = arg min

µ∈D ad s N (µ) = J (Y N (µ)) s.t. Y N (µ) ∈ Y N N : Ã(Y N (µ), W ; µ) = F(W ; µ), ∀ W ∈ Y N N .
Thanks to the (considerably) reduced dimension O(N ) ≪ O(N ) of the systems obtained from RB approximation, we can provide both reliable results and rapid response in the real-time and multi-query contexts. In particular:

-Reliability is ensured by rigorous a posteriori estimations for the error in the RB approximation w.r.t. truth FE discretization (see e.g. [START_REF] Rozza | Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations[END_REF][START_REF] Rozza | Reduced basis approximation and error bounds for Stokes flows in parametrized geometries: roles of the inf-sup stability constants[END_REF]); -Rapid response is achieved by an Offline-Online computational strategy and a rapidly convergent RB space assembling, based on a greedy algorithm. To achieve this goal, RB methods rely on the assumption of affine parametric dependence 1 in A(•, •; µ) and F(•; µ).

Hence, in an expensive Offline stage we prepare a very small RB "database", while in the Online stage, for each new µ ∈ D, we rapidly evaluate both the field and the output (with error bounds) whose computational complexity is independent of FE dimension N .

Efficient shape parametrization techniques for geometrical complexity reduction

In general, shape optimization problems feature more difficulties than optimal control problems, such as shape deformation, shape derivatives and the evaluation of shape-dependent quantities: a crucial aspect of optimal shape design is thus the geometrical treatment of the shapes during the optimization process. Common strategies for shape deformation involve the use of (i) the coordinates of the boundary points as design variables (local boundary variation) or (ii) some families of basis shapes combined by means of a set of control point (polynomial boundary parametrizations). These techniques are not well suited within the RB framework, since a global transformation T (•; µ) is needed, rather than a boundary representation [START_REF] Rozza | Model order reduction by geometrical parametrization for shape optimization in computational fluid dynamics[END_REF]. A more versatile parametrization can be introduced by exploiting the free-form deformation (FFD) techniques, in which the deformations of an initial design, rather than the geometry itself, are parametrized [START_REF] Lassila | Parametric free-form shape design with PDE models and reduced basis method[END_REF]. In this case, the shape parametrization is constructed on a regular lattice of control points, by combining the deformations acting on a subset of active control points through a basis of (tensor products of) Berstein polynomials. Input parameters are given by the deformations of the active control points, which have to be properly chosen, following some problem-dependent criteria [START_REF] Manzoni | Shape optimization of cardiovascular geometries by reduced basis methods and free-form deformation techniques[END_REF].

Despite its flexibility, the FFD techniques do not satisfy any interpolation property and control points must reside on a regular lattice. In order to overcome these possible limitations, other different techniques based on interpolation properties may be recovered. In particular, we have been focusing on the radial basis functions (RBF) techniques [START_REF] Manzoni | Model reduction techniques for fast blood flow simulation in parametrized geometries[END_REF], which are traditionally used for nonlinear multidimensional interpolation on scattered data (for example in image registration). With respect to FFD techniques, RBF techniques allow a better local boundary control and a free choice of the position of the control points (also on the boundary of the shape domain).

Application in haemodynamics: real-time blood flow simulations in parametrized cardiovascular geometries

The framework based on the coupling between FFD or RBF techniques (or other low-dimensional shape parametrizations) and RB methods has turned out to be useful also for a real-time simulation of blood flows in arterial vessels which might show a deep variation in geometrical configuration, as for example carotid artery bifurcations. Our goal is twofold:

spanning a variety of carotid configurations through low-dimensional shape parametrizations [START_REF] Bressloff | Parametric geometry exploration of the human carotid artery bifurcation[END_REF], and shape registration of parametrized carotid shapes from patient data measured in the form of flow velocities;

real-time simulation of blood flows in reconstructed geometries and computing indices related to arterial occlusion risk and highly dependent on geometrical configurations, possibly for predictive surgery applications.

In the first approach we might minimize some discrepancy functional between the simulated velocity and the observed velocity in an atlas-based variational data assimilation method (see e.g. [START_REF] Mcleod | Atlas-based reduced models of blood flows for fast patient-specific simulations[END_REF]); in the latter we minimize a cost function such as the viscous energy dissipation

J (Y (µ)) = ν 2 Ω |∇u(µ)| 2 dΩ.
to obtain carotid shapes exhibiting the least disturbance to the blood flow, being Y (µ) = (u(µ), p(µ)) the velocity and the pressure of the fluid, respectively.

Validation of the reduced basis methodology

A first numerical test has been performed exploiting a coupled FFD+RB framework on a simple geometrical configuration (see Fig. 1), given by a stenosed carotid artery parametrized with respect to the displacement of two control points (p = 2) located close to the bifurcation (see e.g. [START_REF] Manzoni | Model reduction techniques for fast blood flow simulation in parametrized geometries[END_REF] for further details about representation of carotid bifurcations). Flow simulations through a steady Stokes model show a remarkable dependence of the flow even on small variation of the shape configuration. In particular, our interest has been focused on the evaluation of an output related both with the flow and the shape, given by the viscous energy dissipation. Some details concerning the reduced basis spaces are listed in Tab. 1; we remark the strong reduction in the system dimensions and a large computational speedup, concerning performances for each new geometrical configuration, of about two orders of magnitude. We provide a certification of the accuracy of the methodology: in Fig. 2 the true errors between the FE and the RB approximation are reported, the related error bounds (see [START_REF] Rozza | Reduced basis approximation and error bounds for Stokes flows in parametrized geometries: roles of the inf-sup stability constants[END_REF] for error bound expression and derivation), as well as the error between the FE and the RB output. We observe fast, nearly exponential convergence in N . Furthermore, the a posteriori error bounds are both reliable and reasonably effective.

A comparison between FFD and RBF parametrizations

Next, we report here some preliminary results on the comparison between a FFD and a RBF setting defined on the carotid configuration already introduced. Also in this case we are interested in the evaluation of the viscous energy dissipation; we just compare the two parametrizations techniques by considering a Navier-Stokes model for the fluid flow. Both the settings deal with p = 6 parameters, given by the vertical displacements of some selected control points; in the FFD case we introduce a 6 × 8 lattice of control points, while in the RBF case we introduce in total 12 control points close to the bifurcation and at the extrema (see Fig. 3), using the thin-plate spline (TPS) and the Gaussian shape functions [START_REF] Buhmann | Radial Basis Functions[END_REF]. In this last case, we deal with the displacement of the six control points located at the center of the configuration.

We compare the shapes obtained by minimizing the viscous energy dissipation: in Tab. 2 are reported the results for the two cases, while the configurations corresponding to the minimum values of the viscous energy dissipation are represented in Fig. 4. We can remark that the three shapes are quite similar, as well as the output reduction. The number of iterations taken by the optimization procedure is comparable among the three options. Regarding the main qualities of these two shape parametrization tools, the RBF technique proves to be more versatile and accurate for this kind of applications -it enables to choose freely the location of control points rather than selecting the most relevant control points on the regular FFD lattice, as well as to impose interpolation constraints -even if its construction and the computation of the related parametrized tensors is much more difficult. Not only, by considering the same amount, location and available displacements of control points, the Gaussian RBF is found to be more suitable for describing local deformations; the TPS option allows to get a more global and regular deformation, where an enhanced shape smoothness is ensured by a minimization of the bending deformation energy property, fullfilled by this kind of RBF.

Table 2. Results for the minimization of the viscous energy dissipation obtained by using the FFD and the RBF settings introduced above.

Conclusion and perspectives

The capability of the reduced basis method to solve shape registration and optimization problems involving incompressible flows in real-time looks promising if coupled with an efficient and versatile geometrical parametrization. The integration of the RBF parametrization technique within the reduced basis framework, as well as its application to blood flow simulation on geometries reconstructed from patient data, looks promising in its flexibility and ability to express a variety of shape deformations. Further elements that may be explored deal with the uncertainty quantification [START_REF] Lassila | A reduced computational and geometrical framework for inverse problems in haemodynamics[END_REF] and/or robust optimization and control problems [START_REF] Lassila | Boundary control and shape optimization for the robust design of bypass anastomoses under uncertainty[END_REF] for patient-specific scenarios.

Fig. 1 .

 1 Fig. 1. Schematic diagram of the FFD setting; bold control points can be freely moved in vertical direction and used as parameters representing small deformations.

Fig. 2 .

 2 Fig. 2. Left: • error estimation (natural norm) and • true error between RB and FE approximation; right: true error between FE and RB output (vorticity).

  Fig.3. FFD setting (top) and RBF setting (bottom); for each case, parameters are given by the displacements of the selected (blue) control points.

Fig. 4 .

 4 Fig. 4. Optimal configuration obtained by minimizing the viscous energy dissipation for the FFD case (top), the RBF case with the thin-plate spline option (middle) and the RBF case with the Gaussian option (bottom).

Table 1 .

 1 RB + FFD for the carotid artery bifurcation: numerical details.

	Number of FE dof Nv + Np	24046
	Number of RB functions N	16
	Number of design variables P	2
	Linear system dimension reduction 500:1
	FE evaluation t online F E	(s)	2.8039
	RB evaluation t online RB	(s)	0.0231

  [START_REF] Deparis | Reduced basis method for multi-parameter-dependent steady Navier-Stokes equations: Applications to natural convection in a cavity[END_REF]. FFD setting (top) and RBF setting (bottom); for each case, parameters are given by the displacements of the selected (blue) control points.

		FFD	RBF (thin-plate) RBF (Gaussian)
	output reduction	39, 1%	45, 9%	36, 7%
	iterations	84	117	91
	parameters	6 (48)	6 (12)	6

If this assumption does not hold, it could be recovered in through an intermediate empirical interpolation process.

⋆ Work supported by Swiss National Science Foundation (SNSF) under grant 200021-122136 and European Research Council under Project Mathcard ERC-2008-AdG-227058.


Jan Pustelnik

Numerical shape optimization via dynamic programming

Keywords: dynamic programming, numerical approximation, contact problem, shape optimization, sufficiency optimality condition, structural optimization, topological derivative, stationary Navier-Stokes equations

In this paper we describe a novel framework for finding numerical solutions to a wide range of shape optimization problems. It is based on classical dynamic programming approach augmented with discretization of the space of trajectories and controls. This allows for straightforward algorithmic implementation. This method has been used to solve a well known problem called the "dividing tube problem", a state problem related to fluid mechanics, that requires simultaneous topology and shape optimization in case of elastic contact problems and involves solving the Navier-Stokes equations for viscous incompressible fluids.

Introduction

In the paper, as a model problem, we consider state problems related to fluid mechanics, namely the Navier-Stokes equations for viscous incompressible fluids. The main problem is to search for optimal shape of a given objective. For an incompressible fluid, conservation laws for momentum and mass are assumed to be in force. The displacement field of the body is governed by the Reynoldsaveraged Navier-Stokes (RANS) equations with an algebraic mixing length turbulence. The volume of the body is assumed to be bounded.

The results pertaining the existence, regularity and finite-dimensional approximation of solutions to mentioned problems are given in [4], [START_REF] Hlavaček | Solving of variational Inequalities in Mechancs[END_REF]. The primaldual algorithms for numerical solving of contact problems were developed in [START_REF] Hüber | A Primal-Dual Active Set Algorithm for Three Dimensional Contact Problems with Coulomb Fristion[END_REF], [START_REF] Stadler | Semismooth Newton and Augmented Lagrangian methods for a Simplified Friction Problem[END_REF]. In the course of solution the necessary optimality condition for simultaneous shape and topology optimization is formulated, while the shape and topological derivatives are employed, what stays close to classical optimization problems and gives sufficient optimality conditions. It is a different approach than the one applied e.g. in [START_REF] Fulmański | Soko lowski A Level Set Method in Shape and Topology Optimization for Variational Inequalities[END_REF] where the notion of topological derivative and results concerning its application in optimization of elastic structures are reported. We describe a new numerical algorithm for that optimization problem.

General shape optimization problem

We consider the following shape optimization problem, which is being analyzed and subsequently solved in [START_REF] Fulmański | Pustelnik Dynamic programming approach to structural optimization problem -numerical algorithm[END_REF]:

minimize J(Ω) = Ω L(x, u(x), ∇u(x))dx (1)
subject to

Ω ∈ Θ, Au(x) = f (x, u(x)), Bu(x) = φ(x) on ∂Ω ( 2
)
where Θ is a certain family of bounded with C 0,1 boundary domains of D ⊂ R n which will be defined precisely in subsection 2.1 and A is a differential operator e.g. defining Navier-Stokes equations and B an operator acting on the boundary. We assume that L : R n × R × R n → R is Lipschitz continuous with respect to all variables, f : R n × R m → R m is continuous and Lipschitz continuous with respect to last variable, φ(•) is continuous on ∂Ω.

Reduction of shape optimization problem to classical control problem

We will summarize in this subchapter the results of research published in [START_REF] Fulmański | Pustelnik Dynamic programming approach to structural optimization problem -numerical algorithm[END_REF] in order to introduce relevant objects on which the presented numerical method operates.

Let U be a given nonempty, compact set in C 0,1 of surfaces defined on E ⊂ R n-1 . We assume that each supremum of each subfamily of U also belongs to U as well as finite concatenation of element of U belongs to U . Let U ∋ v → Ω(v) be a given family of simply connected domains in D ⊂ R n with C 0,1 boundary such that some fixed part of their boundary is changing and is a surface v from U . We assume that Ω(v) depends in a smooth way on v and that there exists a v max ∈ U such that Ω(v) ⊂ Ω(v max ), for all v ∈ U and there exists a v min ∈ U such that Ω(v min ) ⊂ Ω(v), for all v ∈ U . Let us denote the part of the boundary ∂Ω(v min ) corresponding to the surface v min as Γ 0 while that corresponding the surface v as Γ (v). We have the following BVP:

Find

z max ∈ C 1,k (Ω(v max )) such that ∆z(x) = 0 in Ω(v max )\ Ω(v min ), z(x) = 0 on Γ 0 , z(x) = 1 on Γ (v max ). Next for v ∈ U , v = v min , find z ∈ C 1,k (Ω(v)\ Ω(v min )) such that: ∆z(x) = 0 in Ω(v max )\ Ω(v min ), z(x) = 0 on Γ 0 , z(x) = z max (x) on Γ (v). Solutions to the BVP above belong to C 1,k (Ω(v)\ Ω(v min )
)and in fact they are restrictions of z max to Ω(v)\ Ω(v min ). Because z(x) depends (in a continuous way) on v, we will use the notation z(x, v).We define the family Θ of sets over which the problem (1)-( 2) is considered as: Θ = {Ω(v) : v ∈ U }. The sets from Θ are called admissible sets for problem (1)- [START_REF] Fulmański | Soko lowski A Level Set Method in Shape and Topology Optimization for Variational Inequalities[END_REF]. For a given Ω(v)\ Ω(v min ), we introduce the field and the deformation:

V (x, v) = ||∇z(x, v)|| -2 ∇z(x, v), T (w, v) = x(s, w, v), s ∈ [0, 1], where x(•, w, v) is a solution to d ds x(s, w, v) = V (x(s, w, v), v), s ∈ [0, 1], x(0, w, v) = w,
w ∈ Γ 0 . Notice, that for a given fixed w ∈ Γ 0 , the point x(1, w, v) belongs to Γ (v). For a given control v ∈ U we can write:

d ds x(s, w, v) = V (x(s, w, v), v), s ∈ [0, 1], x(0, w, v) = w. (3)
Then the boundary Γ (v) is the image of Γ 0 by the map x(1, •, v). Thus, for a given v = v min , we have an alternative definition of the Ω(v)\ Ω(v min ):

Ω(v)\ Ω(v min ) = {x : x = x(s, w, v), 0 < s < 1, w ∈ Γ 0 }.
This means that we can construct and study some objects over the set Ω(v) with the help of the family

F (v): F (v) = {x(s, w, v) : 0 < s < 1, w ∈ Γ 0 }. The original functional J(Ω) in terms of the family F (v) can be rewritten as J(F (v)) = I(v), I(v) = Ω(vmin) L(y, u(y), ∇u(y))dy + Ω(v)\ Ω(vmin) L(x, u(x), ∇u(x))dx = Ω(vmin) L(y, u(y), ∇u(y))dy + 1 0 Γ0 L(x(s, w, v))dwds, where L(x(s, w, v)) = L(x(s, w, v), u(x(s, w, v)), ∇u(x(s, w, v))) ∂ ∂s x ∂ ∂w
x . Therefore we are able to reduce the original shape optimal control problem to classical optimal control problem (P): minimize

I(v) subject to d ds x(s, w, v) = V (x(s, w, v), v), s ∈ [0, 1], x(0, w, v) = w, w ∈ Γ 0 , v ∈ U , Ω(v) ∈ Θ, Au(x) = f (x, u(x)), x ∈ Ω(v). (4)
In order to formulate any sufficient optimality conditions for this problem we apply classical dynamic programming scheme.

Dynamic programming approach as a method to solution of (P)

Let us take any x ∈ Ω(v max )\ Ω(v min ) and denote by U x a subfamily of U such that x ∈ v for each v ∈ U x . Next denote by v x = max U x , where the maximum over U x means that Ω(v) ⊂ Ω(v x ) for all v ∈ U x . By our assumption on U , v x exists and v x ∈ U . Put Ūx = {v ∈ U : Ω(v) ⊂ Ω(v x )}. By (3) for each v ∈ U x there is a trajectory x(•, w, v) such that x = x(1, w v , v), for some w v ∈ Γ 0 . The problem (P) falls into the category of Lagrange control problems treated in many books (e.g. [START_REF] Fleming | Deterministic and Stochastic Optimal Control[END_REF]). Following Chapter IV of this book we define a value function for (P), for x ∈ Ω(v max ):

S(x) = inf Ω(vmin)
L(y, u(y), ∇u(y))dy

+ 1 0 Γ0 L(x(s, w, v))dwds , ( 5
)
where infimum in ( 5) is taken over all pairs (x(•, w, v), v)

satisfying d ds x(s, w, v) = V (x(s, w, v), v), s ∈ [0, 1], v ∈ Ūx , w ∈ Γ 0 and for v ∈ U x , x(1, w v , v) = x, for some w v ∈ Γ 0 . Each pair (x(•, w, v), v
) satisfying these equations will be called admissible for the point x ∈ Ω(v max )\ Ω(v min ). However, in practice, we cannot expect that S(•) is of C 1 in Ω(v max )\ Ω(v min ), this is why we are interested in numerical approximation of S(•). Therefore, we shall look for ε-value function:

S ε (•). For given ε > 0 we call S ε : Ω(v max )\ Ω(v min ) → R, ε-value function if S(x) ≤ S ε (x) ≤ S(x) + ε, x ∈ Ω(v max )\ Ω(v min ). ( 6
)
It is clear that there exists infinitely many ε-value functions S ε (•).

Numerical approximation of the value function

This section is an adaptation of the method developed by Pustelnik in his Ph.D. thesis [START_REF] Pustelnik | Approximation of optimal value for Bolza problem[END_REF] for numerical approximation of value function for Bolza problem from optimal control theory. Let us define the following set T = x : x ∈ Ω(v max )\ Ω(v min ) . Since Ω(v max )\ Ω(v min ) is bounded, the set T is compact. Let T ∋ x → g(x) be an arbitrary function of class C 1 in T such that g(x) = c, x ∈ Γ 0 , where c is some constant which will be determined later. For a given function g, we define

(x, v) → G g (x, v) as G g (x, v) = g x (x)V (x, v) - Γ0 L(x(1, w, v))dw, (7)
v ∈ Ūx , where x(•, w, v), u are defined as previously. Next, we define the function x → F g (x) as

F g (x) = max G g (x, v) : v ∈ Ūx . (8)
Note that by the assumptions on L and V , the function F g is continuous in T . By the continuity of F g and compactness of T , there exist k d and k g such that k d ≤ F g (x) ≤ k g for x ∈ Ω(v max )\ Ω(v min ).

Definition of covering of T

Let η > 0 be fixed and {q η j } j∈Z be a sequence of real numbers such that q η j = jη, j ∈ Z (Z -set of integers). Denote J = {j ∈ Z : there is x ∈ T , jη < F g (x) ≤ (j + 1)η}, Next, let us divide the set T into the sets P η,g j , j ∈ J, as follows

P η,g j := x ∈ T : q η j < F g (x) ≤ q η j+1 , j ∈ J.

Discretization of F g

Define in T a function

h η,g (x) = -q η j+1 , x ∈ P η,g j , j ∈ J. (9)
Then, by the construction of the covering of T , we have

0 ≤ F g (x) + h η,g (x) ≤ η, x ∈ T. (10)
Let (x(•, w, v), v) be any admissible pair with the trajectory defined in [0, 1], starting at the point x(0, w, v)), w ∈ Γ 0 fixed. We show that there exists an increasing sequence of m points {τ i } i=1,...,m ,

τ 1 = 0, τ m = 1, such that for τ ∈ [τ i , τ i+1 ] |F g (x(τ i , w, v)) -F g (x(τ, w, v))| ≤ η 2 , i = 2, . . . , m -2, ( 11
)
|F g (x(τ 2 , w, v)) -F g (x(τ, w, v))| ≤ η 2 , τ ∈ (τ 1 , τ 2 ], |F g (x(τ m-1 , w, v)) -F g (x(τ, w, v))| ≤ η 2 , τ ∈ [τ m-1 , τ m ).
Indeed, it is a direct consequence of two facts: Lipschitz continuity of x(•, w, v) with a common Lipschitz constant and continuity of F g . From (11) we infer that for each i ∈ {1, . . . , m -1} if x(τ i , w, v) ∈ P η,g j for a certain j ∈ J, then we have for τ ∈ [τ i , τ i+1 )

x(τ, w, v) ∈ P η,g j-1 ∪ P η,g j ∪ P η,g j+1 . Define

h η,g (x(τ 1 , w, v)) = h η,g (x(τ, w, v)) for some τ near τ 1 , h η,g (x(τ m , w, v)) = h η,g (x(τ, w, v)) for some τ near τ m . Thus for τ ∈ [τ i , τ i+1 ] h η,g (x(τ i , w, v)) -η ≤ h η,g (x(τ, w, v)) ≤ h η,g (x(τ i , w, v)) + η, (12)
and so, for i ∈ {2, . . . , m -1}

h η,g (x(τ i , w, v)) -h η,g (x(τ i-1 , w, v)) = η i x(•,w,v) , (13)
where η i x(•,w,v) is equal to -η or 0 or η. Integrating (12) we get

| 1 0 h η,g (x(τ, w, v))dτ - i∈{1,...,m-1} [h η,g (x(τ i , w, v))(τ i+1 -τ i )]| ≤ η.
By using the formula above and the following simple arithmetic transformations i∈{2,...,m-1}

[h η,g (x(τ i , w, v)) -h η,g (x(τ i-1 , w, v))](τ m -τ i ) = i∈{1,...,m-1} [h η,g (x(τ i , w, v))(τ i+1 -τ i )] -h η,g (x(τ 1 , w, v))(τ 1 -τ m ), we obtain i∈{2,...,m-1} [h η,g (x(τ i , w, v)) -h η,g (x(τ i-1 , w, v))](τ m -τ i ) + h η,g (x(τ 1 , w, v))(τ m -τ 1 ) -η(τ m -τ 1 ) ≤ τm τ1 h η,g (x(τ, w, v))dτ ≤ i∈{2,...,m-1} [h η,g (x(τ i , w, v)) -h η,g (x(τ i-1 , w, v))](τ m -τ i ) + h η,g (x(τ 1 , w, v))(τ m -τ 1 ) + η(τ m -τ 1 ).
and, taking into account (13), we infer that i∈{2,...,m-1}

η i x(•,w,v) (τ m -τ i ) + h η,g (x(τ 1 , w, v))(τ m -τ 1 ) -η(τ m -τ 1 ) ≤ τm τ1 h η,g (x(τ, w, v))dτ (14) ≤ i∈{2,...,m-1} η i x(•,w,v) (τ m -τ i ) + h η,g (x(τ 1 , w, v))(τ m -τ 1 ) + η(τ m -τ 1 ).
We would like to stress that ( 14) is very useful from numerical point of view: we can estimate the integral h η,g (•, •) along any trajectory x(•, w, v) as a sum of finite number of values, where each value consists of a number from the set {-η, 0, η} multiplied by τ mτ i . Moreover, for two different trajectories:

x(•, w 1 , v 1 ), x(•, w 2 , v 2 ), the expressions i∈{2,...,m-1} η i x(•,w 1 ,v 1 ) (τ m -τ i ) + h η,g (x(τ 1 , w 1 , v 1 ))(τ m -τ 1 )
and i∈{2,...,m-1}

η i x(•,w 2 ,v 2 ) (τ m -τ i ) + h η,g (x(τ 1 , w 2 , v 2 ))(τ m -τ 1 ) are identical if h η,g (x(τ 1 , w 1 , v 1 )) = h η,g (x(τ 1 , w 2 , v 2 )) (15)
and

η i x(•,w 1 ,v 1 ) = η i x(•,w 2 ,v 2 ) for all i ∈ {2, . . . , m -1}. ( 16
)
The last one means that in the set B of all trajectories x(•, w, v), w ∈ Γ 0 , v ∈ U , we can introduce an equivalence relation r: we say that two trajectories x(•, w 1 , v 1 ) and x(•, w 2 , v 2 ), w 1 , w 2 ∈ Γ 0 , v 1 , v 2 ∈ U are equivalent if they satisfy (15) and ( 16). We denote the set of all disjoint equivalence classes by B r . The cardinality of B r , denoted by ||B r ||, is finite and bounded from above by 3 m+1 . Define

X = x = (x 1 , . . . , x m-1 ) : x 1 = 0, x i = η i x j , i = 2, . . . , m -1, x j ∈ B r , j = 1, . . . , ||B r || .
It is easy to see that the cardinality of X is finite.

The considerations above allow us to estimate the approximation of the value function.

Theorem 1. We have the following estimation

min x∈Br,w0∈Γ0 - τm τ1 h η,g (x(τ, w 0 , v))dτ -g(x(τ m , w 0 , v)) ≤ max x∈Br τm τ1 - Γ0 L(x(s, w, v))dw ds -g(x(τ 1 , w 0 , v)) ≤ max x∈Br,w0∈Γ0 - τm τ1 h η,g (x(τ, w 0 , v))dτ -g(x(τ m , w 0 , v)) + η(τ m -τ 1 ),
where u is a solution to (4) for Ω(v).

Proof. By inequality (10) 0 ≤ F g (x) + h η,g (x) ≤ η we have -h η,g (x) ≤ F g (x) ≤ -h η,g (x) + η. Integrating the last inequality along any x(•, w 0 , ṽ) in the interval [τ 1 , τ m ] we get

- τm τ1 h η,g (x(τ, w 0 , ṽ))dτ ≤ max v∈U τm τ1 g x (x(τ, w 0 , ṽ))V (x(τ, w 0 , ṽ), v) - Γ0 L(x(τ, w, v))dw dτ ≤ - τm τ1
h η,g (x(τ, w 0 , ṽ))dτ + η(τ mτ 1 ).

Hence, we get two inequalities min x∈Br,w0∈Γ0

τm τ1 h η,g (x(τ, w 0 , ṽ))dτg(x(τ m , w 0 , ṽ))

≤ min x∈Br,w0∈Γ0 max v∈U τm τ1 (-g(x(τ m , w 0 , ṽ)) +g x (x(τ, w 0 , ṽ))V (x(τ, w 0 , ṽ), v) - (-g(x(τ m , w 0 , ṽ))

+g x (x(τ, w 0 , ṽ))V (x(τ, w 0 , ṽ), v) - Γ0 L(x(τ, w, v))dw dτ ≤ max x∈Br,w0∈Γ0 - τm τ1
h η,g (x(τ, w 0 , ṽ))dτg(x(τ m , w 0 , ṽ)) + η(τ mτ 1 ).

As a consequence of the above we get min x∈Br,w0∈Γ0

τm τ1 h η,g (x(τ, w 0 , ṽ))dτg(x(τ m , w 0 , ṽ))

≤ max x∈Br τm τ1 - Γ0 L(x(τ, w, v))dw dτ -g(x(τ 1 , w 0 , v)) ≤ max x∈Br,w0∈Γ0 - τm τ1
h η,g (x(τ, w 0 , v))dτg(x(τ m , w 0 , ṽ)) + η(τ mτ 1 )

and thus the assertion of the theorem follows.

Let us now define four following symbols:

F (η,1) (x) := -i=2,...,m-1 η i x (τ m - τ 1 ), F 1 (x) := -i∈{1,...,m-1} x i (τ m -τ 1 ), F (η,i) (x) := -i=2,...,m-1 η i x (τ m -τ i ), F i (x) := -i∈{1,...,m-1} x i (τ m -τ i ). Now,
) := -h η,g (x(τ 1 , w 0 , v))(τ m -τ 1 ) -g(x(τ m , w 0 , v)). Taking into account (14) we get min x∈X {F i (x)} + min x∈Br,w0∈Γ0 {H(x, w 0 )} -η(τ m -τ 1 ) ≤ min x∈Br - τm τ1 h η,g (x(τ, w 0 , v))dτ -g(x(τ m , w 0 , v)) ≤ min x∈X {F i (x)} + max x∈Br,w0∈Γ0 {H(x, w 0 )} + η(τ m -τ 1 )
and a similar formula for supremum. Applying that to the result of the theorem above, we obtain the following estimation

min x∈X {F i (x)} + min x∈Br,w0∈Γ0 {H(x, w 0 )} -2η(τ m -τ 1 ) ≤ max x∈Br τm τ1 - Γ0 L(x(τ, w, v))dw dτ -g(x(τ 1 , w 0 , v)) (17) ≤ max x∈X {F i (x)} + max x∈Br,w0∈Γ0 {H(x, w 0 )} + η(τ m -τ 1 ).
Thus, we come to the main theorem of this section, which allows us to reduce an infinite dimensional problem to the finite dimensional one.

Theorem 2. Let η > 0. Assume that there is θ > 0 and v such that

max x∈X {F i (x)} (18) + max x∈Br,w0∈Γ0 {H(x, w 0 )} ≤ min x∈X {F i (x)} + min x∈Br,w0∈Γ0 {H(x, w 0 )} + θ(τ m -τ 1 ), min x∈Br,w0∈Γ0 {H(x, w 0 )} = min w0∈Γ0 {-h η,g (x(τ 1 , w 0 , v))(τ m -τ 1 ) -g(x(τ m , w 0 , v))} Then (η + θ)(τ m -τ 1 ) + min x∈Br,w0∈Γ0 {H(x, w 0 )} + min x∈X {F i (x)} (19) is ε-optimal value at (τ 1 , w 0 ) for ε = 2η + θ with g(w) = Ω(vmin) L(y, ū(y), ∇ū(y))dy, w ∈ Γ 0 ,
where ū is a solution to (4) for Ω(v).

Proof. From the formulae (17), (18) we infer {H(x, w 0 )} + η(τ mτ 1 ) + θ(τ mτ 1 ).

Next, using the definition of value function ( 5), we get (19).

3.3

The algorithm for numerical solution of (P)

In the previous section the last theorem allows us to estimate an ε-optimal value of function (see [START_REF] Hüber | A Primal-Dual Active Set Algorithm for Three Dimensional Contact Problems with Coulomb Fristion[END_REF]) for problem (P). As can be seen from the formulas (18) and (19) the essence of the approximation is to be able to calculate the value of the following expressions: x i (τ mτ i )    .

To achieve this aim we construct a particular directed weighted graph G, in which the weight of every edge is the value of the expression x i (τ mτ i ). This graph has following properties 1. Every path has length of m -1 edges. 2. Every two vertices connected by an edge correspond to points (τ, x 1 ) and (τ + ∆τ, x 2 ) such, that the point x 2 is reachable from the point x 1 in the next unit of time τ + ∆τ .

Therefore by identifying in the graph G the path with lowest (greatest) cost we find the value of the expression inf

  min x∈X {F i (x)} + min x∈Br,w0∈Γ0 {H(x, w 0 )} -2η(τ mτ 1 ) (τ, w, v))dw dτ -Ω(vmin)L(y, ū(y), ∇ū(y))dy≤ max x∈X {F i (x)} + max x∈Br,w0∈Γ0

  ..,m-1}x i (τ mτ i )

  we use the definition of equivalence class to reformulate the theorem above in a way that is more useful in practice. To this effect let us note that, by definition of equivalence relation r, we have

		min x∈Br	F (η,1) (x) = min x∈X	{F 1 (x)} , max x∈Br	F (η,1) (x) = max x∈X	{F 1 (x)} .
	Let	us	now	also	define	the	following	auxiliary	symbol
	H(x, w 0

The algorithm for generation of the graph G 1. Let B -a set of trajectories -be a finite set Bezier curves. 2. Let P be a set of points. At the beginning, P contains only one point: p = w 0 , where w 0 is any but fixed point such that w 0 ∈ Ω 0 3. Create in graph G node which corresponds to point p. 4. Calculate F g (x 0 ) (equation ( 8)) where x 0 = (0, w 0 ) , i.e. find a Bezier curve α ∈ B which minimize value of G g (x 0 , •).