2nd International Conference on Engineering Optimization September 6-9, 2010, Lisbon, Portugal
Design of an Airfoil from a Given Pressure Distribution Using an Approximate Inverse Operator ˇ ak, Jaroslav Pelant Jan Sim´ Aeronautical Research and Test Institute, Beranov´ ych 130, 199 05 Prague - Letˇ nany, Czech Republic
[email protected],
[email protected]
Abstract This contribution describes a numerical solution of an inverse problem for a flow around an airfoil. Its aim is to design an airfoil shape based on a given pressure distribution on its surface. The pressure distribution affects forces acting on the airfoil, so it is possible through its modification to obtain shapes which satisfy requested parameters. The method is suitable to solve problems in subsonic regimes. The main idea of the method is a combination of a direct and an approximate inverse operator. Using both operators, the resulting shape is corrected in each iteration until the desired precision is achieved. The direct operator means, in this notation, a solution of a 2D flow around the airfoil. There is some degree of freedom when choosing the model. In this case, the Navier-Stokes equations completed with a k-omega turbulence model are used. A stationary solution of this system is obtained by an implicit FVM. The approximate inverse operator is based on the thin-airfoil theory for a potential flow. The airfoil is described by a mean camber line and a thickness function, which ensures the capability to deal not only with thin airfoils. Numerical examples are presented. Keywords: inverse problem, CFD, turbulence, airfoil
1. Introduction This paper concerns a numerical method for a solution of an airfoil design inverse problem. Using this method, it is possible to get an airfoil shape that satisfies a prescribed pressure distribution on its surface under given flow conditions. It can be used as a part of some optimization process. The literature describes various approaches to the solution of this problem. The geometrical inverse design methods use usually some auxiliary equations, based on the actual pressure distribution, to correct the shape. The main idea of the presented method is to construct an inexact inversion, a mapping between a pressure distribution and an airfoil shape. The inversion in its original form is almost exact for a potential incompressible flow. Hence this mapping is used in an iterative process until the desired airfoil is obtained. It is possible to couple this inexact inversion with an arbitrary flow solver. The method is aimed to a subsonic flow, the angle of attack is one of the results of the method. The method presented here is an extension of methods described in previous works [3], [6], [7]. The original method deals with an inviscid flow, later with a laminar viscous flow. The current method deals with a turbulent viscous flow described by the Navier-Stokes equations equipped with the k − ω turbulence model which improves its applicability. In the following text a short description of the method is given.
2. Inverse problem As was said above, the goal of the method is to find an airfoil shape corresponding to a prescribed pressure distribution. The method utilizes the possibility of derivation of an approximate inversion. Thus the problem could be written as: Find a pressure pseudo-distribution p such that P L(p) = f ,
(1)
where f is the given pressure distribution, P is an operator representing a solution of a flow problem and finally L represents an approximate inversion to P. This problem is solved by the method of successive iterations, where the solution of the Eq.(1) is a limit of the sequence ∞ {pk }k=0 , pk+1 = pk + α (f − P Lpk ) . (2) 1
The parameter α is a positive real number chosen such that the sequence converges. According to experiences from numerical results, the choice α = 0.6 is sufficient to achieve convergence in most cases. Lower values ensures better convergence but also increases the number of iterations. The approximate inverse operator L is derived using the thin airfoil theory. The details of this can be found in [3]. In short, the airfoil is composed of a mean camber line and a thickness function, mathematically written ψ1 (x)
=
ψ2 (x)
=
x ± t(x) q
s′ (x) 1 + s′ 2 (x)
s(x) ∓ t(x) q
,
1 1 + s′ 2 (x)
,
x ∈ h0, 1i .
(3)
In this notation the airfoil coordinates are denoted by (ψ1 , ψ2 ), the upper sign is for the upper part of the airfoil and the bottom sign for the bottom part. The functions s(x) (=the mean camber line) and t(x) (=the thickness function) are derived on the idea of artificial distributions of vorticity γ and sources q along the chord. The velocities induced by this distributions are described by Z 1 γ(ξ) 1 dξ, (4) uy (x) = 2π 0 ξ − x Z 1 1 q(ξ) ux (x) = dξ. (5) 2π 0 x − ξ Under assumption that the velocity is tangent to the airfoil, the following result can be obtained Z 1 1 − ξ x dξ − (uup (ξ) − ulo (ξ)) ln s(x) = 2π 0 ξ Z 1 x − ξ 1 dξ , − (uup (ξ) − ulo (ξ)) ln (6) 2π 0 ξ p Z 1 1 uup (ξ) + ulo (ξ) 1 + (ξ − xξ)/(x − xξ) p − 1 ln (7) t(x) = dξ . 1 − (ξ − xξ)/(x − xξ) π 0 2
The variables x and ξ run along the normalized chord line, given by the interval h0, 1i. The symbols uup and ulo represent velocity distributions on the airfoil surface, normalized by the free stream velocity. These functions need not represent a physically relevant distribution, they are related to the sequence (2) instead. Hence, similar to the notation used in Eq.(1), they can be called pseudo-distributions. The transformation between the pressure pseudo-distribution and the velocity pseudo-distribution, in the formula denoted as u(x), is the following: 2 (γ−1)/γ ! 2 u(x) p(x) 2/M∞ +γ−1 1− . (8) = u∞ γ−1 p0
The symbol p0 is the pressure at zero velocity, M∞ is the Mach number in the free stream, u∞ is the free stream velocity and γ is the Poisson adiabatic constant. From the construction of the airfoil coordinates, it is clear that uup and ulo need to be functions. Since the chord line doesn’t need to connect the leftmost point with the rightmost one, simply taking the x-coordinate of a point on the surface doesn’t satisfy this requirement, in general. From that reason the distribution is assumed along the mean camber line, with the x-coordinate as the leading variable x. The mean camber line is evaluated in each iteration, so the only additional expense is the inversion of the mapping Eq.(3). The integrals in Eqs.(6) and (7) are evaluated using a suitable numerical quadrature. Care must be taken, since the integrands have singularities for ξ = x, ξ = 0 and ξ = 1. From this reason a ChebyschevGauss quadrature was chosen, where the interval is divided by the roots of a Chebyschev polynomial. The integrating formula is Z
0
1
N/2 p 2π X f (xi , x2k−1 ) x2k−1 (1 − x2k−1 ), f (xi , ξ) dξ ≈ N k=1
2
i = 0, 2, 4 . . . , N
(9)
with xi = 1/2 (1 + cos (iπ/N )) , i = 0, 1, . . . , N. 3. Flow problem The viscous compressible flow around an airfoil is described by the system of the Navier-Stokes equations. Since most of the flow in real situations is turbulent, the laminar model seems insufficient. To improve the quality of the predicted flow and also the stability of the method, a k − ω model of turbulence is included (see [2], [8]). 3.1. Model of the flow The system of the equations can be written in a vector form 2
2
X ∂R ~ j (w, ~ ∇w) ~ ~ ∂w ~ X ∂ F~j (w) ~ (w, = +S ~ ∇w) ~ , + ∂t ∂x ∂x j j j=1 j=1
(10)
where w ~
=
F~j (w) ~ = ~ j (w, R ~ ∇w) ~ =
~ w, S( ~ ∇w) ~ =
T
(̺, ̺v1 , ̺v2 , E, ̺k, ̺ω) ,
(11) T
(̺vj , ̺v1 vj + δ1j p, ̺v2 vj + δ2j p, (E + p)vj , ̺kvj , ̺ωvj ) , µ ∂e µT ∂k 0, τj1 , τj2 , τj1 v1 + τj2 v2 + γ + , (µ + σk µT ) , Pr PrT ∂xj ∂xj T ∂ω , (µ + σω µT ) ∂xj T 0, 0, 0, 0, Pk − β ∗ ̺ωk, Pω − β̺ω 2 + CD .
(12)
(13) (14)
Using the common notation, ̺ denotes a density, p is a pressure, v1 , v2 are velocity components, E is an energy, k is a turbulent kinetic energy and finally ω is a specific turbulent dissipation. The symbol Pr represents the Prandtl number (the subscript T denotes the turbulence). The viscosity coefficient µ is evaluated using the Sutherland’s formula. The symbol µT denotes an eddy viscosity coefficient, which is given by the formula ̺k µT = . (15) ω The stress tensor in the N.-S. equations is given by relations 2̺k 2 ∂v2 4 ∂v1 − − , τ11 = (µ + µT ) 3 ∂x1 3 ∂x2 3 2 ∂v1 2̺k 4 ∂v2 τ22 = (µ + µT ) − − + , (16) 3 ∂x1 3 ∂x2 3 ∂v2 ∂v1 . + τ12 = τ21 = (µ + µT ) ∂x2 ∂x1 The production of the turbulence Pk and the production of the dissipation Pω are expressed as ∂v2 ∂v1 ∂v1 ∂v2 Pk = τ 11 + τ 22 + τ 12 + , ∂x1 ∂x2 ∂x1 ∂x2 Pk Pω = αω ω , k where τ ij = τij for µ = 0. Finally, the cross-diffusion term CD is given by the relation ∂k ∂ω ̺ ∂k ∂ω + ,0 . CD = σD max ω ∂x1 ∂x1 ∂x2 ∂x2
(17) (18)
(19)
√ The turbulence model is closed by parameters β ∗ = 0.09, β = 5β ∗ /6, αω = β/β ∗ − σω κ2 / β ∗ (where κ = 0.41 is the von K´ arm´an constant), σk = 2/3, σω = 0.5 and σD = 0.5. This choice of parameters resolves the dependence of the k − ω model on the free stream values [2]. If the turbulent kinetic energy k is set to zero, the turbulence model has no influence upon the N.-S. equations and the laminar model is described. 3
3.2. Numerical treatment The system of equations mentioned above is solved by the implicit finite volume method. The variables are normalized using critical values of the density, velocity and pressure. The resulting dimensionless system has the same form as the original one and thus no modification to the system is needed. The computational domain is discretized by a structured quadrilateral C-type mesh. Since the coupling between the equations describing the flow and the equations describing the turbulence is only by the viscous terms, it is possible to solve the problem in two parts [8]. In the first part (continuity equation, momentum equations, energy equation), the variables k and ω are assumed time independent. Similarly, in the second part (k − ω equations), the variables p, ̺, v1 , v2 are held constant in time and the system of two equations is solved with respect to the unknowns k and ω. These systems can be solved independently of each other. This approach reduces computational costs and allows to easily modify a laminar solver into a turbulent one. The linearized system of algebraic equations is solved by the GMRES method (using software SPARSKIT2 [4]). The convective terms F~i are evaluated using the Osher-Solomon numerical flux in the case of the flow part and by the Vijayasundaram numerical flux in the turbulent part. A higher order reconstruction based on the Van Albada limiter is also implemented. The numerical evaluation of a gradient on an edge of two cells is based on the values in the centres of the six neighbouring cells. Since a suitable angle of attack has to be found in order to satisfy the condition on the position of the stagnation point on the leading edge, the airfoil is rotated round a chosen point. 4. Numerical examples The first example shows an asymmetric case. The method was examined on the NACA4412 airfoil. The Reynolds and Mach numbers are Re = 6 · 106 and M∞ = 0.6. Because of the restriction on the stagnation point, the angle of attack is set α∞ = 1.8 ◦ . The obtained airfoil together with the obtained pressure distribution after 40 iterations is in Fig. 1. The next Fig. 2 shows a distribution of the error between the computed and prescribed pressure. A convergence history of the L2 -norm of error kp − f k2 is shown in Fig. 3.
1.4 0.8 1.2 0.6 1 P
0.4 0.01
0.8
0.2
Y
0
p-f
0.005
0.6
0
0.4
0
0.2
0.4
X
0.6
0.8
1
-0.005
0
0.2
0.4 X 0.6
0.8
1
Figure 1: Example 1 - resulting airfoil and pressure Figure 2: Example 1 - difference between the obdistribution (normalized) tained and prescribed pressure along the chord The second example is an airfoil shape resulting from a by-hand prescribed distribution. The Mach number M∞ = 0.7 and Reynolds number Re ≈ 15 · 106 . The pressure distribution after 30 iterations together with the resulting shape are in Fig. 6. The angle of attack is α∞ = 0.82 ◦ . The error along the chord is in Fig. 7. 5. Conclusion An extension of a numerical method for a solution of an inverse problem for a flow around an airfoil was described. This method is suitable for a design based on a prescribed pressure distribution in the case of 4
||p-f||
10
-1
10-2 10
-3
0
5
10
15
20 25 iterations
30
35
40
Figure 3: Example 1 - convergence history of an error kp − f k2
1.4
naca 4412 1st it. 2nd it. 5th it. 20th it.
0.1 0.05
1.2
Y
P
1
0
naca 4412 1st it. 2nd it. 5th it. 20th it.
0.8 0.6
-0.05 0
0.2
0.4 X 0.6
0.8
0.4
1
0
0.2
0.4 X 0.6
0.8
1
Figure 4: Example 1 - airfoil shapes during iteration Figure 5: Example 1 - pressure distributions during process iteration process
1.4 0.8 1.2 0.6 1 P
0.4 0.01
0.8
0.2
Y
0
p-f
0.005 0.6
0 0.4
0
0.2
0.4
X
0.6
0.8
1
-0.005
0
0.2
0.4 X 0.6
0.8
1
Figure 6: Example 2 - Resulting airfoil and pressure Figure 7: Example 2 - Distribution of the error distribution (normalized) along the chord
5
a 2D turbulent viscous flow. It is necessary to prescribe the distribution with care because it is possible to happen that the solution doesn’t exist. Thus the method can be used to modify existing airfoil with some knowledge. Since the angle of attack is one of the results, the method is applicable in the cases, where a specific angle is not demanded. Acknowledgements This work was supported by the Grant MSM 0001066902 of the Ministry of Education, Youth and Sports of the Czech Republic. The authors acknowledge this support. References [1] M. Feistauer, J. Felcman and I. Straˇskraba, Mathematical and Computational Methods for Compressible Flow, Clarendon Press, Oxford, 2003. [2] J. C. Kok, Resolving the Dependence on Freestream Values for the k − ω Turbulence Model, AIAA Journal, 38 (7), 1292-1295, 2000. ´ Report Z–69, Prague, [3] J. Pelant, Inverse Problem for Two-dimensional Flow around a Profile, VZLU 1998. [4] Y. Saad, Iterative Methods for Sparse Linear Systems, 2nd ed., SIAM, 2003. ˇ ak, Solution of 2D Navier-Stokes Equations by Implicit Finite Volume Method and Application [5] J. Sim´ ´ Report R–4003, Prague, 2006. in Inverse Problem, VZLU ˇ ak, J. Pelant, Solution of an Airfoil Design Problem With Respect to a Given Pressure Distri[6] J. Sim´ ´ Report R–4186, Prague, 2007. bution for a Viscous Laminar Flow, VZLU ˇ ak, J. Pelant, A contractive operator solution of an airfoil design inverse problem, PAMM, 7, [7] J. Sim´ 2007 (DOI: 10.1002/pamm.200700136). [8] D. C. Wilcox, Turbulence Modeling for CFD, 2nd ed., DCW Industries Inc., 1998.
6
