SATURATION-FREE NUMERICAL SCHEME FOR COMPUTING THE FLOW PAST A LATTICE OF AIRFOILS WITH A SHARP EDGE Текст научной статьи по специальности «Математика»

Russian Journal of Nonlinear Dynamics, 2019, vol. 15, no. 2, pp. 135-143. Full-texts are available at http://nd.ics.org.ru DOI: 10.20537/nd190203

NONLINEAR PHYSICS AND MECHANICS

MSC 2010: 76B10

Saturation-Free Numerical Scheme for Computing the Flow Past a Lattice of Airfoils with a Sharp Edge

A. G. Petrov

The Joukowski - Chaplygin condition, which allows us to determine the circulation of the flow past contour with a sharp edge, is one of the most important achievements of S. A. Chaplygin, whose 150 birthday is celebrated this year. This research is devoted to this problem.

We consider the flow past of a lattice of airfoils by a potential fluid flow. The profile line is determined parametrically by an equation in the form of two dependences of the Cartesian coordinates on the parameter. For a smooth closed loop the Cartesian coordinates are periodic analytical functions. Their Fourier series coefficients decrease exponentially depending on the harmonic number. In the meantime, for a sharp edge loop they decrease much slower — inversely to the square of the harmonic number. Using the symmetric continuation of the profile with a sharp edge, a method of presenting it as a Fourier series with exponentially decreasing coefficients is proposed. Based on this idea, a quickly converging numerical scheme for the computation of the flow past airfoils lattice with a sharp edge by a potential fluid flow has been developed.

The problem is reduced to a linear integro-differential equation on the lattice contour, and then, using specially developed quadrature formulas, is approximated by a linear system of equations. The quadrature formulas converge exponentially with respect to the number of points on the profile and can be rather simply expressed analytically.

Thanks to its quick convergence and high accuracy, this method allows one to optimize profiles by using a direct method by any given integral characteristics. We can find the distribution of the shear stress and the breakaway point on the calculated velocity distribution on the profile from the solution of the boundary layer equation. In the method suggested we do not need to perform the difficult work of constructing the lattice. Also, the problem of scheme viscosity at high Reynolds numbers is omitted.

Keywords: lattice of airfoils, potential fluid flow, sharp edge loop

Received February 25, 2019 Accepted June 04, 2019

This research was carried out within the state assignment of FASO of Russia (state registration No. AAAA-A17-117021310382-5), supported in part by RFBR (project No. 17-01-00901).

Alexander G. Petrov petrovipmech@gmail.com

Ishlinsky Institute for Problems in Mechanics of the Russian Academy of Sciences, prosp. Vernadskogo 101-1, Moscow, 119526 Russia

Introduction

The basics of hydrodynamic lattice theory were founded in the classical works of N. E. Jou-kowski and S. A. Chaplygin. Huge steps in the development of this theory were made in Weinig's [1] and N. E. Kochin's [2] monographs. The basics of the lattice theory are described in the monographs of L.I. Sedov [3], N.E.Kochin et al. [4], and in more detail in the monograph of G. Yu. Stepanov [5]. To construct a solution, the classical theory uses conformal mappings. Using them, we may construct a solution for simple lattices — a plate and a circle.

At present, the lattice theory is developed using numerical methods. In Gostelow's monograph [6] the lattice theory is explained based on rich experimental material. Also, the calculation results are based on difference methods and the finite element method.

The topic of lattices of airfoils with sharp edge flow around them is very popular at present both in Russia [9, 10] and abroad [11, 12].

1. Potential flow past a lattice of airfoils

In [7] a numerical scheme for forces acting on the lattice profiles is proposed. The periodic lattice of airfoils is an infinite sequence of airfoils Lj, j = ... — n, ..., —1, 0, 1,2, ..., n, ... The airfoil Lj+l is obtained by translating the airfoil Lj along the axis x through the lattice spacing h (Fig. 1).

Fig. 1. Lattice of airfoils.

The inlet velocities are = U sin a, = —U cos a. The outlet velocities are V^x =

= Vqqx + r/h, V^y = v^y. The arithmetic mean velocities are

Ux (v<xx + V^x)/2 v<xx +

r

2Ti

Uy = (v^y + V^y)/2 = V

Here r is the circulation, and the subscripts x and y denote the vector projections on the axes x and y, respectively.

By using these quantities, the force exerted by the fluid on the airfoil can be precisely determined by applying the theorem on momentum variation. The force projections Fx and Fy onto the axes are given by the formulas (see [3])

Fx = rUy, Fy = rUx

(1.1)

Integral equation. Green's function G(x, y, x', y') is expressed in terms of the real part of an analytic function of complex arguments 2 = x + iy, z' = X + iy'

G = Re ln

2 A(z' - z)

ASm 2

- ln 2

[ch(A(y; - y)) - cos(X(x' - a:))]

(1.2)

Through it the integral operator kernel of the system of integral equations for the velocity distribution on the lattice profiles is expressed [7]

1

AV(s) = 2n(C - xUy + yUx), J V(s)ds = -r,

r

TT___ 7) TT — n,

L x 2/7 oo.tj L y ' ooyi

l

AV (s) = -J G(s, s')V (s')ds'.

(1.3)

Given the circulation r and the inlet velocities and we can determine the airfoil velocity distribution V(s), the mean velocities U1 and U2 and the constant C by solving system (1.3). If the airfoil has a sharp trailing edge, the circulation can be determined using the additional Kutta- Joukowski - Chaplygin condition on the rear stagnation.

0

u

2. Parametrization of the integral identities

Regarded as a function of s', the integrand of A defined on the closed contour L by its argument s' has a natural period equal to the total contour length l. It is convenient to specify the coordinate of point M (s) on the contour by the parametric equations x = ), y = ), where ) and ) are functions of ( with a unit period. For numerical computations, the contour is partitioned by a finite number of points M1, M2, MN = M0 such that the point Mi corresponds to Z = Q = i/N, i = 1, 2, ..., N.

The parameter Z and the coordinate s are related by the differential relation

1

ds = la(Z)dZ, j a(Z)dZ = 1, (2.1)

0

where a(Z) is the density of contour points; its values at Z = Zi can be calculated using a quadrature formula for a periodic function (see below):

N

Si = sJidx/dQ* + (dy/dQf, l = °i = J si- (2-2)

i=1

3. Quadrature formulas for contour integrals

As mentioned above, the integrands A in the integral operator have period l and quadrature formulas for periodic functions can be used for them (see [8]). For analytic periodic func-

tions F(x) of period l, a quadrature formula has the form

l N

J F{x)dx = Xj = (3.1)

And for a function with a periodic logarithmic singularity the quadrature formula has the form (see [8])

f n l N

J F(x) In sin j(x - Xi) dx = - j\)F(Xj)>

j=i

«(m) = -iln2 + ^+^icos^J, M = N/2,

m = 0, 1, 2, ..., N. Both (3.1) and (3.2) are exact for trigonometric polynomials of the form

t-, a0 M 2nkx 2nkx , .

F = Y + l^akcos—— +bksm——. (3.3)

k=1

According to the terminology introduced by K. I. Babenko, quadrature formulas with this property are referred to as saturation-free. By using such formulas, the analytical properties of a function can be numerically determined from the order of decrease of its remainder.

If F(x) is an infinitely differentiable periodic function, then the remainder of the Fourier partial sum (4.3) has an order of smallness higher than any power. Therefore, the order of smallness of the remainder of the quadrature formulas is higher than any power 1/Nn.

Function G(Z, (') has a logarithmic singularity at (' = (. The quadrature formulas (3.1) and (3.2) yield

N1 Av«i) = Y.AvVv Aij = + G^lfj,

j=1

Gij = G(xi^ Vu xj, Vj), i = j Gii = ln(l/i/n), .

nm

/3(m) = -In sin— +a(m), m = I, 2, ..., N - 1, /3(0) = a(0).

4. System of linear equations for computing the potential flow past a lattice of airfoils

With the help of (3.4) system (1.3) is approximated by the linear system of equations

Nl E ¿ijVj = MC ~ xiUy + Viux), = "r'

j=i N (4.1)

r

2h

TT___ 7) TT — 7)

L X r\l lOO.TJ L y ' ooy

0

V

N>

This is a system of N+3 linear equations for determining N airfoil velocity values V1, V2, ..., two mean velocities Ux, Uy and, and the constant C.

Consider the problem in which the rear stagnation point on a smooth airfoil is a given quantity, which is used to calculate the circulation. For this purpose, we use the Joukowski-Chaplygin (Kutta) condition. Let the point M(xm,ym) denote the sharp trailing edge of the airfoil. In the computational scheme used, the airfoil is assumed to be smooth and the sharp edge is replaced by an airfoil point with a maximally high curvature. Then, according to the Joukowski-Chaplygin condition, the velocity at this point vanishes (stagnation point): Vm = 0. To calculate the circulation from this condition, the velocity at the airfoil boundary is partitioned into three terms (the possibility of this partition follows from the linearity of the problem):

V=WL = uv{1)_uv{2)+tv{3)

dn x y

(4.2)

where the flows with velocity distributions V(1) and V(2) have zero circulation, and the flow with V (3) has a circulation equal to —1. Thus, functions V(k), k = 1, 2, 3 can be found in a unique way. From (4.1), we obtain the following systems of equations for calculating Vjk) at the airfoil points Xj, yj and the constants C1, C3:

N

N

(V(1))j = 2n(Ci + Vi), i = 1, 2,...,N, (V(1))j = 0,

j=i

N

j=1

N

(V(2))j = 2n(C2 + Xi), i = 1, 2,...,N, (V(2))j = 0,

j=1

N

j=1

l

N

j=1 N j=1 The condition Vm = 0 yields the circulation equation

(4.3)

Vm = Ux (Vi)m - Uy (V2 )m + r^ = 0.

(4.4)

5. Profile with sharp edge expansion in Fourier series

The Fourier series theory implies that the Fourier series coefficients an, bn of a periodic function f (s), whose first derivative is discontinuous, decrease as 1/n2 as n ^ and for the derivative as 1/n.

Thus, the use of Fourier series for functions x(s) and y(s) which have a derivative discontinuity on the sharp edge is not very convenient. A substantial improvement of the convergence of the Fourier series coefficients may be obtained as follows.

Let the digitized profile be given by the Cartesian coordinates xn, yn, n = 1, 2, ..., N so that the sharp edge corresponds to the origin of coordinates (Fig. 2). Then, if we continue by symmetry the profile curve, we will obtain a self-intersecting line at the origin of coordinates in the shape of a figure-of-eight, as shown in Fig. 3. The points of this curve are defined as follows:

xn xn, XN+n —XN-n;

yn = yn, yN +n = —yN-n; (5.1)

n = 1, 2, ..., N.

3 2 1

Fig. 2. Profile with sharp edge.

Fig. 3. Symmetrically continued profile.

The advantage of this curve compared to the initial one consists in the fact that on a double period 0.2l it is determined by smooth infinitely differentiable functions X(s) and y(s). On a double period they can be calculated using quickly converging Fourier series. Furthermore, symmetry considerations imply that they will contain only sinuses

m

X — ^ an

n=1 1 2N

an = ~N ^n

an sin m, y — Y^ bn sin nY,

n=l

2N

2nni I v-^ _

n=1

2N

n=1

2nni

»sinw

(5.2)

To obtain 3 digits accuracy, it is enough to take 20 expansion terms. For y G (0, 2n) we get the whole curve in the form of a figure-of-eight, as shown in Fig. 3. On the half-period y G (0, n) we obtain a half of the figure-of-eight, that is, the analytic form of the initial profile (Fig. 2) with a sharp edge.

A more accurate description of velocity V(s) in the vicinity of the sharp edge may be obtained by accumulating the grid points in the vicinity of the sharp edge. This can be done using the following substitutions Yn — nF(Zn), (n — n/N, n — I, 2, ..., N:

F(i) = i-£sin(27ri).

(5.3)

The closer the parameter p is to one, the more dense are the points placed near the sharp

edge.

6. Dependence of the circulation error on the nodes number

The velocity at the profile boundary may be found from the system of integral equations (4.1)-(4.4). Determining the discrete distribution with the help of the Fourier series (5.1)-(5.2) on the half-period and the accumulation of points using (5.3), by approximation (3.4) we can calculate the matrix with high accuracy with a relatively small number of grid points. Then, solving systems (4.3), using (4.4) we may calculate the circulation with high accuracy.

To check the accuracy of our calculus, we can use the theory of generalized Joukowski profiles. They can be constructed as follows. The equation

2 - z0e^ (Z - Z0e^\z°/Z°

z + z0e^ \Z + Z0e^ implies that the complex variable z can be expressed in terms of Z

jjl + viZ) _ iZe-<*-Z„\'«/Z°

The complex variable Z is expressed in terms of Z"

Z = (Z0 - ae~il3 + Z")e*. (6.2)

The boundary points on the plane z correspond to a circle of radius a on the plane Z". The form and the position of the profile on the plane z are determined by five parameters Z0, z0, 3, a, 0. The angle of the sharp edge equals t = n(2 — z0/Z0).

The following parameters correspond to the profile from Fig. 2:

Z0 = 0.95, t = 0.4, / = 0.1, a = 1, 0 = -0.3, z0 = Z0(2 - t/k).

The complex flow past the potential of the generalized Joukowski profile on the plane z is easily built using the conformal mapping method. In the flow past of the profile with the velocity at infinity Ue1®1 the complex potential W(z) is found as follows. We calculate the complex potential on the plane Z"

upi® a2 r

X(Z") = Ue~»Z" + + In Z", Q = Qx-i>

and substitute into it function Z". The complex potential on the physical plane is w(z) = = X(Z "(z)).

The circulation is found from the condition that the velocity at the sharp edge equals zero. On the plane Z" the point Z" = ae~%^ corresponds to it. Thus, we obtain

r = -4nUa sin(0 + /). (6.3)

It is convenient to use this exact solution for testing the numerical scheme. The circulation by the described numerical method is calculated as follows.

The discrete set of points zn, n = 1, 2, ..., N can be obtained by substituting into the mapping (6.2) and then into (6.1) the points of circle Z^ = e2mn/N, n = 1, 2, ..., N. Thus, by the method described in Section 5, we accumulate the points in the vicinity of the sharp edge. We calculate the matrix Aij by the approximation formulas (3.4) and then, solving three systems

of equations (4.3), we calculate three fundamental solutions k = 1, 2, 3. The circulation is calculated by formula (4.4).

Let us analyze how the error in the calculation of circulation decreases depending on the increase in the number of grid points.

For the first comparison let us take the profile from Fig. 2, that is, the parameters of the profile are Z0 = 0.95, t = 0.4, /3 = 0.1, a = 1, 0 = 0.

For the angle 6 = 0 and the unit flow velocity U = 1 we obtain the exact value of circulation r = -1.2545437, and for 6 = 0.5 we obtain r = -7.09550658. The results of calculus are presented in the array. In the first and second lines, the values of circulation errors for N = 32; 48; 64 and 80 are presented:

N 32 48 64 80

6 = 0 Ar -6 • 10"4 -1 • 10"4 -7 • 10"5 -1 • 10"5 6 = 0.5 Ar 9 • 10"4 7.7 • 10"5 7 • 10"5 6.7 • 10"5.

For the second comparison let us consider a finer profile with parameters Z0 = 0.97, t = 0.2, 3 = 0.05, a = 1, 0 = 0. The profile is presented in Fig. 4. For the angle 6 = 0 and the unit velocity of the flow U = 1 we obtain a new value of circulation r = -0.6280567, and for 6 = 0.5 we get r = -6.352719.

Fig. 4. Second profile with sharp edge.

The results of calculus are presented in the array. In the first and second lines, the values of circulation errors for N = 32; 48; 64 and 80 are presented:

N 32 48 64 80

e = 0 Ar 1.7 • 10"3 -1.4 • 10"3 -3 • 10"4 -1.1 • 10"4 e = 0.5 Ar -3 • 10"3 -6 • 10"5 3 • 10"4 1.6 • 10"5

References

[1] Weinig, F., Die Strömung um die Schaufeln von Turbomaschinen: Beitrag zur Theorie axial durch-strömter Turbomaschinen, Leipzig: Barth, 1935.

[2] Kochin, N.E., Hydrodynamics Lattice Theory, Moscow: Gostekhizdat, 1949 (Russian).

[3] Sedov, L. I., Two Dimensional Problems in Hydrodynamics and Aerodynamics, New York: Wiley, 1965.

[4] Kochin, N.E., Kibel, I.A., and Roze, N. V., Theoretical Hydrodynamics, New York: Wiley, 1964.

[5] Stepanov, G. Yu., Hydrodynamics of Turbomachine Lattices, Moscow: Fizmatlit, 1962 (Russian).

[6] Gostelow, J. P., Cascade Aerodynamics, Oxford: Pergamon, 1984.

[7] Petrov, A. G., Saturation-Free Numerical Scheme for Computing the Flow past a Lattice of Airfoils and the Determination of Separation Points in a Viscous Fluid, Comput. Math. Math. Phys, 2011, vol. 51, no. 7, pp. 1239-1250; see also: Zh. Vychisl. Mat. Mat. Fiz, 2011, vol. 51, no. 7, pp. 1326-1338.

[8] Petrov, A. G., Quadrature Formulas for Periodic Functions and Their Application to the Boundary Element Method, Comput. Math. Math. Phys., 2008, vol.48, no. 8, pp. 1266-1283; see also: Zh. Vychisl. Mat. Mat. Fiz, 2008, vol.48, no. 8, pp. 1344-1361.

[9] Mamaev, I. S. and Vetchanin, E. V., The Self-Propulsion of a Foil with a Sharp Edge in a Viscous Fluid under the Action of a Periodically Oscillating Rotor, Regul. Chaotic Dyn., 2018, vol. 23, nos. 78, pp. 875-886.

[10] Mamaev, I. S., Tenenev, V. A., and Vetchanin, E. V., Dynamics of a Body with a Sharp Edge in a Viscous Fluid, Rus. J. Nonlin. Dyn., 2018, vol. 14, no. 4, pp. 473-494.

[11] Mason, R. J., Fluid Locomotion and Trajectory Planning for Shape-Changing Robots, PhD Dissertation, Pasadena, Calif.: California Institute of Technology, 2003, 264pp.

[12] Tallapragada, P., A Swimming Robot with an Internal Rotor As a Nonholonomic System, in Proc. of the American Control Conf. (ACC, Chicago, Ill., July 2015), pp. 657-662.