Drag minimization for the obstacle in compressible flow using shape derivatives and finite volumes

In the paper the shape optimization problem for the static, compressible Navier-Stokes equations is analyzed. The drag minimizing of an obstacle immersed in the gas stream is considered. The continuous gradient of the drag is obtained by application of the sensitivity formulas derived in the works of one of the co-authors. The numerical approximation scheme uses mixed Finite Volume - Finite Element formulation. The novelty of our numerical method is a particular choice of the regularizing term for the non-homogeneous Stokes boundary value problem, which may be tuned to obtain the best accuracy. The convergence of the discrete solutions for the model under considerations is proved. The non-linearity of the model is treated by means of the fixed point procedure. The numerical example of an optimal shape is given.

1. Introduction. One of the main applications of the theory of compressible viscous flows [27] is the optimal shape design in aerodynamics. The classical example is the minimization problem of the drag of an airfoil travelling in the atmosphere with uniform speed.
In the paper we present numerical computations for a simple geometry in the plane. The theory leading to the formulas for the shape derivatives is based on a series of papers by P.I. Plotnikov and J. Sokolowski [23]- [26], while computational scheme is based on [28]. For the approximation of solutions the method of Finite Volumes is used, what required supplementing the results of [10], [8] with some additional analysis. This analysis constitutes the main novelty of the paper.
The optimal control of compressible flows is still poorly investigated, in contrast to incompressible flows. We can cite here see e.g. [7], [29] or a more recent work [6] where the dynamic version of the problem was considered. 2. Mathematical model. Let B ⊂ R 2 be a doubly connected, bounded hold all domain with a smooth boundary Σ = ∂B. We denote two connected components of ∂B by Σ in and Σ out , Σ = Σ in ∪ Σ out . Suppose that B contains an obstacle S, which is a compact set with a sufficiently smooth boundary Γ. Suppose also that the volume of S is given and small with respect to the volume of B and the obstacle is far from the boundary of B. We define the computational domain as Ω = B \ S. Thus the boundary of Ω consists of three parts: ∂Ω = Σ in ∪ Σ out ∪ Γ. We assume that a viscous, compressible fluid is flowing across Ω. Fluid enters Ω through Σ in and leaves through Σ out . The computational domain is presented in Fig. 1. We restrict our consideration to stationary flows in the framework of the existence theory established in [27], Theorem 11.2.6 on page 311. In fact we need the existence of unique strong solutions to the model under considerations in order to assure the weak differentiability of solutions with respect to small parameter.
The dynamics of the fluid is determined by three functions: the velocity field u, the density and the pressure p( ), where p = p( ) is a smooth, strictly monotone function of the density with p(0) = 0. In computational experiments we set p = p 0 ( 0 ) θ , with p 0 = 1 and θ = 5 3 . The quantities: Ma, Re and λ are physical constants, denoting respectively the Mach number (defined as a proportion of the speed of a fluid to the speed of sound), the Reynolds number (which is a proportion of fluid's inertial forces to its friction forces), and the viscosity coefficient. The quantity 0 , present in the boundary condition for , is a constant, in computations we set 0 = 1.2. Finally, U : R 2 → R 2 is a given smooth vector field, satisfying certain compatibility conditions, see Section 9.2.
2.1. General theory in three spatial dimensions. The stationary boundary value problem for compressible Navier-Stokes equations in three spatial dimensions can be formulated as follows.
Problem. Assume that a gas occupies the flow domain Ω = B \ S, where B ⊂ R n , n = 3 in general case, is a bounded hold-all domain, with an obstacle S in its interior.
Find a couple (u, ) such that ∆u + λ∇ div u = Re u · ∇u + Re Ma 2 ∇p( ) in Ω, (1a) div( u) = 0 in Ω, (1b) where U is a given vector field and Σ = ∂B. The inlet Σ in and the outlet Σ out are defined by Σ in = {x ∈ Σ : U · n < 0}, Σ out = {x ∈ Σ : U · n > 0}, respectively. Here n stands for the outward normal to ∂(B \ S) = Σ ∪ ∂S. For the sake of simplicity we assume that the mass force equals zero, and 0 is a given positive constant.

Change of variables in Navier-Stokes equations.
We consider the general case, the results apply as well in two space dimensions. In order to perform the shape sensitivity analysis in three space dimensions we choose a vector field T ∈ C 2 (R 3 ) 3 vanishing in a neighborhood of Σ, and introduce the mapping Φ ε (x) = x + εT(x), which defines the perturbation S ε = Φ ε (S) of the obstacle S. For small ε, the mapping Φ ε takes diffeomorphically the flow region Ω = B \ S onto Ω ε = B \ S ε . In the perturbed domain Problem 2.1 reads Problem. Find a solution (ū ε ,¯ ε ), to the following boundary value problem posed in the variable domain Ω ε = B \ S ε , for the shape parameter ε ∈ (−δ, δ) with δ > 0: Now, we perform a change of variables in equations (2) in order to reduce Problem 2.2 in the variable domain Ω ε , depending on a small parameter ε, to a problem in the fixed domain Ω, named the reference domain. Denote by M the Jacobi matrix of the mapping Φ ε and by g the determinant of M, i.e., M(x) = I + εD T(x), g(x) = det(I + εD T(x)). (3) Further the notation N stands for the adjugate matrix Definition 2.1.
• The shape gradient of the drag functional for stationary Navier-Stokes equations is given by the derivative of the function Introduce the functions The function u ε is called the Piola transformation ofū ε .
Lemma 2.2. Let (ū ε (y),¯ ε (y)) be a solution to Problem 2.2. Then the couple (u ε (x), ε (x)) defined by (5) satisfies Here, the linear operator A and the trilinear form B are defined by Thus we come to the following problem: Here the operators A and B are defined by (7).
Observe that equations (8) depend only on the matrix N and do not depend on ε directly. In this framework, the fact that N is the adjugate matrix to the Jacobi matrix I + εD T is not so important. Hence our task is to show that Problem 2.2 is well-posed and its solution is differentiable with respect to N. We begin by analysing the well-posedness.
2.3. Solution scheme for stationary problem in reference domain. In [27] we use the simple algebraic scheme proposed in [22]. The basic idea is to introduce the effective viscous pressure q ε and a large parameter σ 0 , and write equations (8) in the form System (10) consists of perturbed Stokes equations (10a)-(10b) for (u ε , q ε ) and a perturbed transport equation (10c) for ε . It is easy to see that equations (10b) and (10c) are equivalent to (9) and the mass balance equation div( ε u ε ) = 0. Therefore, equations (10) are equivalent to (8). It was shown in [22] and [23] that problem (10) has a strong solution if the parameters satisfy the conditions which are more restrictive than the standard conditions Re 1, Ma 2 1 adopted in [2], [16], and [24]. Perturbations. We will look for a local solution to problem (8) which is close to an approximate solution (u , q , ). In order to define such a solution notice that for small Mach and Reynolds numbers equations (10a)-(10b) are close to Stokes equations and the density ε is close to 0 . Thus we take = 0 (11) and choose (u , q ) as the solution to the following boundary problem for the Stokes equation: ∆u − ∇q = 0, div u = 0 in Ω, where The solution to problem (8) is decomposed into the approximate solution (u , q , ) and small perturbations (v ε , π ε , ϕ ε , m ε ), with some unknown functions ϑ = (v ε , π ε , ϕ ε ) and an unknown constant m ε . We point out that our decomposition (14) holds exclusively in the reference domain Ω = B \S. The unknown functions ϑ and the unknown constant m ε contain the necessary and sufficient information on the shape dependence of solutions to the equations (8) upon the shape parameter ε → 0. This information is used in order to determine the material derivatives of the solutions to the model (2) with respect to the shape perturbation. However, the Stokes boundary value problem (1) defined in the reference domain Ω is meaningless for the model (2) defined in the variable domain Ω ε .

Function spaces.
To establish existence results for problem (15), we begin with the definition of an admissible class of solutions to this problem. Notice that we are looking for strong solutions, which means that in equations (15) the unknown functions and their derivatives are integrable, i.e., these equations are satisfied a.e. in Ω. To this end it suffices to assume that v ∈ W 2,2 (Ω) and (π, ϕ, ζ) ∈ W 1,2 (Ω). Next, we are looking for solutions such that (π, ϕ, ζ) are continuous and v is continuously differentiable in Ω. In order to satisfy this demand it suffices to take v ∈ W s+1,r (Ω) and (π, ϕ, ζ) ∈ W s,r (Ω) with sr > 3. To address both issues we introduce the Banach spaces equipped with the norms u X s,r = u W s,r (Ω) + u W 1,2 (Ω) , u Y s,r = u W s+1,r (Ω) + u W 2,2 (Ω) .
We intensively exploit the fact that the spaces X s,r and Y s,r are reflexive. The proof of this fact is nontrivial and follows the line given below. It is well known (see [1,Thm. 3.5]) that for any bounded domain Ω with C 1 boundary, any integer m ≥ 0 and 1 < r < ∞, the Sobolev space W m,r (Ω) is reflexive. The space W 1,r 0 (Ω) is obviously reflexive as a closed subspace of a reflexive Banach space (see [1,Thm. 1.21]).

Lemma 2.5.
Let Ω be a bounded domain with C 1 boundary and 0 < s < 1 < r < ∞. Then the spaces X s,r and Y s,r are reflexive.
While W s,r (Ω) and W 1,2 (Ω) are separable, we cannot conclude from this that X s,r is separable. However, the dual spaces (X s,r ) and (Y s,r ) are separable: Lemma 2.6. Let Ω ⊂ R 2 be a bounded domain with C 1 boundary and 0 < s < 1 < r < ∞. Then (X s,r ) and (Y s,r ) are separable.
Existence theory. Define the Banach space and denote by B τ ⊂ E the closed ball of radius τ centered at 0. Next, note that for sr > 3, elements of the ball B τ satisfy the inequality where the norm in E is defined by Further, c denotes generic constants, which are different in different places and depend only on Ω, U, σ and r, s. We assume that the flow domain and the given data satisfy the following conditions.

Condition 1.
• ∂Ω is a closed surface of class C 3 and the set Γ = cl Σ in ∩Σ\Σ in is a closed C 3 one-dimensional manifold such that Σ = Σ in ∪ Γ ∪ Σ out .
Moreover, we can assume that it is extended to R 3 in such a way that the extension U ∈ C 3 (R 3 ) vanishes in a neighborhood of S. • There is a positive constant c such that U · ∇(U · n) > c > 0 on Γ.
Since the vector field U is tangent to ∂Ω on Γ, the left hand side is well defined.
Let us consider the augmented problem of the following form: Problem. Given (u , q , ), find (v ε , π ε , ϕ ε ) and a constant m ε , depending on ε, such that The constant m ε is determined from the following relations: where the auxiliary function ζ ε is a solution to the adjoint boundary value problem The following existence theorem is the main result of this section.
Theorem 2.7. Assume that the surface Σ and the given vector field U satisfy Condition 1. Furthermore, let numbers r and s satisfy Then there exists σ * > 1, depending only on U, Ω and s, r, with the following property. For every σ > σ * there are positive τ * and c, depending only on U, Ω, r, s, and σ, such that whenever Problem 2.4 has a solution ϑ ∈ B τ , ζ ε ∈ X s,r , m ε ∈ R. The auxiliary function ζ ε and the constants κ ε , m ε admit the estimates 3. Minimization of the shape functional. Now we are going to adopt the general theory to the specific problem in two spatial dimensions. Our goal is to propose the numerical scheme for drag minimization and obtain numerical results. We refer the reader to e.g., [13], [14], [30], [31], [32] for incompressible modeling and numerical methods. The drag minimization is considered also in [12].
3.1. Shape functional. It is known that the theory of compressible fluids has many essential applications in aerodynamics. For example, it is used in aircraft construction. Thus the obstacle S can be identified with an airfoil. We work in the reference frame with the fixed obstacle S. (see Fig.1). Aerodynamical force, acting on the obstacle S can be expressed as where n denotes the unit normal field defined on Γ, pointing outward from S, and I stands for the unit matrix. Similarly as in the case of an airfoil moving in the horizontal direction, the drag functional J can be defined as the horizontal component of the aerodynamical force F a . For the horizontal unit vector e 1 = [1, 0], it reads: Then, applying any smooth function η (in the program η is a solution to the Laplace equation) satisfying suitable boundary conditions, the aerodynamical force can be expressed as an integral over Ω. For η ≡ 1 on Γ and η = 0 on Σ, integrating (23) by parts gives: In consequence, the drag functional is: 3.2. Discretization of continuous shape gradient. The existence of shape gradient of the drag functional is shown for compressible model in [27]. In the numerical procedure we are going to use the material derivatives of the solutions to the governing equations. The obtained continuous shape gradient is discretized for the purposes of numerical method. We point out that the shape derivatives of the continuous compressible model also lead to the shape gradient of the cost [26]. The shape derivatives are obtained by the singular limit procedure in the boundary value problem that is formulated for the material derivatives. Formally, the same result can be derived by an application of the boundary variations technique to the strong formulation of the stationary compressible Navier-Stokes equations. However, it is important to remember that even the existence of the weak solutions to the stationary model with nonhomogenous boundary conditions is a difficult problem in the physical range of parameters. Therefore, the formal considerations for compressible model are not in general justified, and are performed in the framework of experimental mathematics. The results obtained here are proved for almost incompressible model.

3.3.
Deformations of Γ and Ω. We seek a deformation of Γ, such that the value of the drag functional J is minimal relatively to a certain family of admissible shapes. We introduce the family of admissible shapes by defining the set of admissible deformations of the initial shape Γ. The set of admissible deformations of Γ consists of linear combinations of basic transformations of Γ, defined as bell-shaped smooth hat functions (see (109)) depending on the discretization of Γ. It is assumed that h = [h 1 , h 2 ] is a smooth vector field normal to Γ, that is a deformation of Γ in the normal direction. It defines a transformation Φ ε (x) = x + εT(x) of Ω into Ω ε , depending on the vector field T constituting the unique solution to the Laplace equation: Observe that, for h = 0 and consequently T = 0, this transformation reduces to identity. Since both boundaries of Ω and boundary data for h 1 and h 2 are smooth, such T is a vector field on Ω of class C ∞ (Ω), normal on Γ and vanishing on Σ. The derivative of the functional J(Ω) in the direction T is given by the limit We use T instead of h in the definition of this derivative because such a formulation could be used for arbitrary field T as well. By applying the same substitutions as in Section 2.2, the shape functional, defined previously as an integral over the domain Ω , can be written as an integral over Ω: 4. Shape derivative. Differentiating formally the integrand in (28), we obtain the formula for the derivative of J, which can be symbolically written as the sum: where the linear forms L e and L u are given by the formulas [25] L e (T) = e 1 · and (δv, δπ, δϕ, δζ, δm) denote the material derivatives d d (v , π , ϕ , ζ , m ) =0 for the solution to the correction system. Since by definition (1) (u , q , ρ ) do not vary then, from (14), we have (δu, δq, δ ) = (δv, δπ, δϕ). The material derivatives (δv, δπ, δϕ, δζ, δm) can be derived by formal differentiation of the correction system (10), which gives the following system for derivatives: − div(δζu) + σ 1 δζ = div(ζδv) + σ 1 T rDT, with the boundary conditions: where and the coefficients b ij equal: It is easy to see that the left hand sides of (19a) and (33)-(37) are identical, which will simplify further computations.
5. Computational scheme. In order to make an optimization step, we need to calculate the shape derivative, which is given by the formula (29). Thus we have to find the unique solution of the system (19a) and also the material derivatives (δv, δπ, δϕ, δζ, δm), that are the solution to the system (33)-(37). Both these systems can be written in the common form: (Ω)) and σ 2 = σ 1 − div u, with the boundary conditions: We rewrite the above system in the compact form: where the coefficients of (39) depend on u that is computed in the previous step (u, in turn, depends on ϑ and u 0 through (14)). Denote by the inverse of the linear operator ξ −→ L(u(ϑ); ξ) such that L(u(ϑ); ξ) = R(ϑ).
In order to find the approximate solution [5], in which an iteration step is defined by the formula: where ϑ n is a value of solution from the previous step, and We take ϑ 0 = 0. Elements of the sequence α n (in the program we set: α n = 1 2+n , n = 0, 1, ..) have to satisfy the following conditions: In every iteration of the fixed point method we find the solution (42) of the system (40) with fixed right-hand sides f 1 1 , f 1 2 , f 2 , ..., f 5 ∈ L 2 (Ω). Let us notice that if we substitute in the system (40) the null operator for R and consider only the first three equations, we obtain the Stokes system (it has to be solved once in every iteration). To solve (40) numerically in every iteration, we apply the modification of the Finite Element (FE), used in [28], replacing it by the mixture of Finite Volume and Finite Element (which we will denote shortly by FV). The latter one was used in [10] and [8] for homogeneous systems of the type (1).
Similarly as FE, the FV method serves for solving systems of differential equations in the weak form. Let us then multiply the first 5 equations of (40) by the test functions t ∈ C ∞ 0 (Ω) and integrate them over Ω. Next, applying integration by parts, we get: (43) The first three equations of (43) constitute a system which is hardest to solve (the last three equations are independent and can be solved separately). The boundary conditions will be discussed later on (cf. (60)).
6. Approximation method. Assume that there is a given triangulation of the domain Ω, so Ω ≈ T ⊂T T . Denote by h the diameter of the biggest triangle in that mesh (that is the diameter of the circumscribed circle). We will denote the set of all triangle edges from T by E, and the sets of all inner edges and outer edges by E int and E ext , respectively. We will approximate the functions appearing in the system (39) by functions from certain finitely dimensional spaces (the same spaces as in [10] and [8]). The approximations of the functions z 1 1 , z 1 2 , z 3 , z 4 and t, except from the test functions from (43/3), are the elements of the space V h : where t i | t j denotes the common edge of triangles t i and t j , P 1 (t i ) -the space of first order polynomials defined on a triangle t i , andk h | cl ( with the basis of characteristic functions over the triangles from T : where n t denotes the number of triangles in the mesh. Thus, approximations of functions from (43) in the above bases will have the form: with real coefficients z 1 Substituting above sums to the system (43), we obtain a system of algebraic equations in the matrix form: where are the values of the first and the second velocity component on the edges of a given triangle t k . The matrices used above are defined in the following way: while the right-hand side vectors are: The elements of the matrices, as well as the components of the right-hand side vectors can be computed using elementary calculations, like the formula for the gradient of linear function. In order to improve numerical properties of the system (49), we add to its third equation the stabilizing terms T ST AB,1 and T ST AB,2 . The term T ST AB,3 is added to the next two equations. These terms are defined as follows.
Let β > 0, γ > 0, then where |t k | denotes the area of k − th triangle. Next, Here e ij = t i |t j (t i |t j is a common edge of triangles t i and t j ), N (t i ) denotes the set of neighbours of t i , and g ij = (s i h i + s j h j ), where h i , h j are heights from t i and t j which are falling on e ij . Let t k ∈ T . Then: if t k has 1 neighbour, 1 4 if t k has 2 neighbours, 1 6 if t k has 3 neighbours, The factor s i h i | e ij | corresponds to the 1 , 1 2 , or 1 3 part of the area of triangle t i . Thus the penalty for discontinuity between neighbouring triangles t i and t j is proportional to their areas. The last term is and the value of small κ > 0 is in the range 10 −4 -10 −5 ; for the definition of S see formulas (50). Owing to the first term, the system (49) becomes invertible (it corresponds to the condition Π(π) = π in (19). T ST AB,2 is a regularizing term and T ST AB,3 is a term that stabilizes the discretized transport equations, equivalent to adding on the left-hand side the Laplace operator multiplied by κ. Thus the final form of the discretized scheme is the following: 7. Correctness of the approximation method. For the convergence of numerical method we are going to establish: 1. The existence of solutions to the discrete scheme (57) by a fixed point argument. 2. The convergence of the sequence of discrete solutions to the solution (43) of the model, taking into account appropriate a priori estimates for discrete solutions. First of all, the linear Stokes system (1) with the nonhomogeneous Dirichlet boundary conditions is transformed by translation to the equivalent system with homogeneous Dirichlet boundary conditions. Recall that the Stokes system under considerations takes the form: Let u satisfy the system: Substituting w = u − u , we get: Of course, the second equation satisfies the compatibility condition if and only if Σ U · ndS = 0, this condition is imposed in the subsection 8.2. It means that we can apply the scheme (40) with homogeneous Dirichlet boundary conditions to the whole systems which we solve, including (1) . Thus, we replace the approximation space V h (defined in (44)) by the space V 0 h : corresponding to zero boundary conditions. Denote by T h a triangulation T of the diameter h, and by (z 1 h , z 2 h ) -the sequence of solutions for (61) in the space (V 0 h ) 2 × L h corresponding to the triangulation T h . The last three equations of the system (57) are solved, after regularization, in a standard way. Therefore we will concentrate only on its first three equations having the form (59). We rewrite them in the slightly different, but analogical form: where v e,h = |e|z 1 e,h n ti , n ti is a normal vector to the edge e, pointing outwards the triangle t i and Ω,b denotes the following sum of integrals: By the symbol · b we will denote the norm of the space (V 0 h ) 2 , given by the formula: 7.1. Existence of a solution to the discrete scheme. By an application of the Leray-Schauder Fixed Point Theorem (cf. [11]) we can show that the discrete scheme (61) admits a solution. We recall the theorem.
Theorem(Leray-Schauder) 1. Let B be a Banach space, T -compact mapping from B × [0, 1] to B, such that T (X, 0) = 0 for all X ∈ B. Assume that there exists a constant M such that: λ). Then the mapping T 1 : B → B of the form T 1 X = T (X, 1) has a fixed point.
h , z 2 h ) satisfies the following system: where A e,h is defined in the same manner as v e,h . We need to show that the mapping T satisfies all the assumptions. For λ = 0 it has the form: Obviously, the only solution to the above system is the trivial one. Thus T satisfies the condition: The zero solution obviously meets the estimate (62) too. Furthermore, T -as a linear operator, acting in finite dimensional space -is compact. In order to prove the estimate (62) for the solution to the system: we can use estimates derived for the scheme (61) (which is equal to (65) for λ = 1 ), that are presented further in this paper (see 7.2.1). Since, employing the analogical technique as for (61), we get the following inequality: (symbols c 2 , c 1 denote certain positive constants, described later -see 7.2.1).In turn, estimates for z 2 h can be derived from the equation (65/2). Let us multiply it by z 2 h and sum over all triangles from T h : After substituting for λ Ω,b z 2 h div z 1 h the identity from (65/1) (with k h = z 1 h ), we get (from definition of norm · b in the space (V 0 h ) 2 ): Since the first and the third term on the left are non-negative, the following inequality must hold: and after applying definition of the norm in L 2 (Ω) and Hölder's inequality, we have: then: and in consequence: if only Otherwise it holds: which already gives an estimate for z 2 h L 2 (Ω) . From (72) results: In the presence of inequalities (66), (73) and the fact that λ ∈ (0, 1), we have: and The right-hand sides of the above inequalities are given and independent of λ, so also the assumption (62) is satisfied. That ensures the existence of a solution to the system (61).

7.2.
Convergence of discrete solutions to the solution of the model. We will show that the sequence (z 1 h , z 2 h ) has an estimate in V 1 0 (Ω) 2 × L 2 (Ω) independent of the diameter h of the mesh. This makes possible to select a weakly convergent subsequence. Furthermore, we have to know that this subsequence converges to the model's solution when the mesh diameter h → 0. 7.2.1. Estimates for discrete solutions. Multiplying the second equation of the scheme (57) by z 2 h and summing it over all triangles from T h , we get: Next, substituting for Ω,b z 2 h div z 1 h dx the identity (61/1) with k h = z 1 h , we get from definition of the norm · b in the space V 0 h : Removing the last two terms on the left-hand side ,we obtain: Employing Hölder's inequality to the right-hand side gives: wherein the last inequality results from Friedrichs-Poincare inequality for the space (V 0 h ) 2 , having a form: where the constant c 1 depends only on Ω and the regularity of the mesh. In consequence: assuming that f 2 h L 2 (Ω) > 0 (otherwise the inequality (78) simplifies). Now we can use the inf-sup stability of the pair of spaces ((V 0 h ) 2 , L h ) (see [10], [8]): where the constant c 2 > 0 is independent of the size of triangulation. Substituting for Ω,b z 2 h div(k h )dx the expression (61/1), we obtain: Next, the following sequence of estimates holds: Then, substituting the estimate from (84) to inequality (83), we get: Employing (81), we obtain the following relation leading directly to the the estimate for z 1 h : Above inequality gives: Of course, the right-hand side of that inequality can be estimated from below by the expression: if only it is non-negative; otherwise we arrive in the estimate of the form: Thus: so finally we have: Simultaneously, in the presence of (85) and (91), we easily obtain a bound for z 2 h L 2 (Ω) : For further computations it will be also necessary to estimate the expression |z 2 h | T h , defined in the space L h in the following way: From the equality (77) follows that: Furthermore, as s i ≥ 1 6 , then 6(s i h i + s j h j ) ≥ (h i + h j ) and: Thus, an estimate for |z 2 h | T h can be obtained on the basis of the estimate of the right-hand side of (94).
By Theorem 3.3 in [10] it follows that the estimates for discrete solutions (z 1 h , z 2 h ) in V 1 0 (Ω) 2 × L 2 (Ω) which are independent of the mesh size imply the strong convergence of the sequence z 1 h in L 2 (Ω) 2 to the limit z 1 from H 1 0 (Ω) 2 and the weak convergence of z 2 h in L 2 (Ω) (these convergences hold for subsequences only).

7.2.2.
Convergence of discrete schemes to the model. We want to show that the scheme: converges to the system of the form: Since the proof of convergence of the first equation in the scheme (96) to the first one in (97) is considered in the paper [10], we can direct the reader to that source.
Regarding the convergence of the second one, we do not have a ready result but we may apply some ideas from [10]. First, multiplying (96/2) by the approximation l h ∈ L h of a smooth function l ∈ C ∞ 0 (Ω) and summing over all triangles from T h , we obtain: denoting: where l ti,h is the mean value of l on the triangle t i . We will show that the first term converges to Ω div z 1 l dx and the other two tend to zero when h → 0. At the beginning, one can notice that the following equalities hold: since: l h → l in L 2 (Ω) (according to Lemma 2.1(3) from [10]), and div z 1 h is bounded in L 2 (Ω) on every triangle t i , as: Next, integrating the second term by parts, we get because it was already demonstrated that z 1 h → z 1 h→0 in L 2 (Ω) 2 to a function from Applying integration by parts one more time, we obtain the needed result for the first term. Obviously, the second term T 2 tends to zero when h tends to zero since, as it was already shown, the expression z 2 h L 2 (Ω) is bounded. As for the term T 3 , from the boundedness of | z 2 h | T h we get T 3 → 0, owing to: Finally, again from Lemma 2.1(3) from [10], it is obvious that: and from the weak convergence of f 2 ti,h in L 2 (Ω) we have: This, with convergence of the terms T 1 , T 2 and T 3 from the left hand side of (98), gives the result. The novelty of the analysis performed above consists in considering the system with f 2 = 0. The reasoning is based on [10], [8], where a nonlinear system was considered, but had to be supplemented by the elements taking into account this non-homogeneity. 8. Example of obstacle shape optimization.
8.1. Shape and size of the computational domain Ω. We define the computational domain Ω as an ellipse with two small circular holes inside(one of them is the inlet, and the other one -the initial shape of the obstacle, see Fig.1), situated symmetrically with respect to its centre. The radii of the considered ellipse are 30 and 16 units. The centres of the circles lie in the horizontal axis of the ellipse, each in the distance of 3,5 from its centre. The length of the radii of both circles equals 1.

8.2.
Boundary conditions -vector field U. Vector field U, that appearsin the boundary condition (2c) is chosen in such a way, that it satisfies the compatibility condition for the equation (1), which reads: It corresponds to the physical condition that the same amount of the fluid enters and leaves Ω. We assume that the fluid passes through the boundary ∂Ω in normal direction. We take U = a (x−x0) x−x0 2 (a ∈ R), which is the velocity field in the solution to the system (116) defined for x 0 , the source of the gas, taken as the coordinates of the centre of Σ in . 8.3. Optimization step. Using the gradient method, we want to find a transformation of Ω acting in the direction of the steepest descent of J. For a fixed discretization of Γ, the vector field T is a linear combination of basic fields T i related to the boundary points i = 1..n b , namely: where t = [t i ], i ∈ 1..n b are real coefficients. Thus the gradient of the functional J contains the derivatives with respect to the basic fields In order to construct these fields, we approximate the boundary Γ by a closed spline curve ξ of class C 2 , passing through all the points of the discretized boundary (which are the vertices of triangles from T that touch Γ) and parametrized by the arc length s: where L is the length of Γ. Since we want to construct normal movements of the boundary, for every point p k we need a formula for the outer normal vector. It has the form: Then we define the hat function around a given boundary discretization point of Γ by the following formula: where dist(·, s k ) is a minimal distance of s from s k : dist(s, s k ) = min{|s − s k |, L − |s − s k |} (110) and the constant w 0 determines the "width" of the function b k (s). To obtain a normal shift of every boundary point p j , we multiply b k (s j ) by the normal vector n j : Applying the formula (27), all vector fields T k on Ω can be evaluated. During our computations we impose two obvious constraints. The first is the constancy of the volume, since otherwise the optimal obstacle would reduce to the single point. If we assume, that the boundary nodal points of S h move by the vector k∈1..n b ζ k h k , then the requirement of fixed volume reduces to The second constraint concerns the position of the mass centre of the boundary of the obstacle. The speed of the flow decreases with the distance from the inlet, therefore, in order to reduce the drag, the obstacle would move as far to the right as it is possible, see Fig. 1. We want to prevent it, because we are interested in the optimal shape at the specified position.

ANNA KAŹMIERCZAK, JAN SOKOLOWSKI AND ANTONI ZOCHOWSKI
The condition for the mass center to be fixed at the origin reads ∂S h x j ds = 0, j = 1, 2.
This may be transformed as follows: ∂S h Here we have assumed that the lengths of boundary edges l k are approximately the same and equal l 0 . This condition holds at the beginning of each optimization step, see next paragraph. Therefore, in view of (112), x k j = 0, j = 1, 2.
Now if we move each nodal point by n b k=1 ζ k h k (p i ), i = 1, . . . , n b , the above condition takes on the form Because of (113), the condition for keeping mass centre of Γ h in place thus becomes h j k (p i ), j = 1, 2.
As a result we must project the gradient g on the intersection of three hyperplanes It is easily done by taking and finding α 0 , α 1 , α 2 from equations g T 1 d = 0, g T 1h 1 = 0, g T 1h 2 = 0.
As it was mentioned above, after each optimization step we distribute s k approximately uniformly along ξ(s), see (107). This simplifies the condition for keeping the obstacle in place, but also prevents other undesirable behaviour. In places of higher curvature the nodal points on Γ h may get nearer and nearer in the process of optimization. Keeping the distances between these points equal prevents them from overlapping during the consecutive optimization steps. The results of the computations are shown in Fig.2.
with the field U generated by a fundamental solution (Stoklet) to the system: where x 0 / ∈ Ω. Thus, the solution to (115) is given by the formula: In order to compute the rates of convergence of u and p, we make 3 subsequent mesh condensations (every time we take h := 1 2 h). Regularity of the mesh is assured by Matlab's algorithm. Denote by u 1 err, u 2 err and perr errors of the velocity components and the pressure, computed in the norm of L 2 (Ω). We take the total error of the velocity as: uerr = u 1 err 2 + u 2 err 2 .
In both methods the same scheme (57) and the same stabilizing terms (52), (53), (56) have been employed. The values of β and γ are chosen in such a way(β = 1, 2, γ = 0, 6) that they minimize the errors of the Stoklet solution. The same values were used in the computations done during the optimization. The convergence rates of the velocity and pressure obtained from comparison with exact solutions are 3.0714 and 0.9158 for FE compared to 2.5763 and 1.5480 for FV. One can remark that error values themselves are not actually important, only the convergence rate -since for the same triangulation the number of discrete variables for FV is approximately 3 times bigger than for FE. The obtained results are presented on the plot below: It can be concluded that FV approximates the pressure more accurately, whereas its advantage in accuracy for the velocity decreases with subsequent mesh condensation steps. However, the accuracy of the pressure is crucial in this case.