A method of computing the potential flow round a solid of revolution in an incompressible fluid

A METHOD OF COMPUTING THE POTENTIAL FLOW ROUND A SOLID OF REVOLUTION IN AN INCOMPRESSIBLE FLUID* V. V. VOINOV, 0. V. VOINOV and A. G. PETROV

Moscow (Received 27 September 1972) (Revised version 11 October 1973) THE PROBLEM of the longitudinal revolution

is considered.

contour of the body is constructed calculation

of the potential

The potential problem of Neumann

potential

flow of an incompressible

Fredholm’s integral equation

fluid round a solid of

of the second kind for the potential

on the

and a method of solving it is presented. The results of the

and the associated masses are compared with the exact solution.

flow of an incompressible

fluid round a body corresponds

for the Laplace equation.

reduced to the solution of an integral equation

to the external

It is known [ 1,2] that this problem can be for the density of a simple layer. A number of

papers exist in which this method is realized numerically,

an advantage of it is that the

two-dimensional

problem. This refers, in particular,

problem is reduced to a one-dimensional

to

[3]. The contour of the streamlined body is specified by a single-valued function, connecting the distance from the axis of revolution with the coordinate along the axis. The non-analytic part of the kernel is integrated

by means of regular quadrature

formulas. The numerical

schemes do

not specially take into account the variation of the distance to the axis along the contour, and near the poles of the body of revolution the calculations are constructed in exactly the same way as near the equator. The consequence of this error in the density of the simple layer near the axis exceeds by an order of magnitude paper concentrates on the following:

the error at other points. Unlike these papers the present

1) instead of the equation for the density of the simple layer Fredholm’s integral equation of the second kind is constructed and solved directly for the potential of the velocity field of the fluid; 2) the contour of the solid of revolution

is specified in parametric

form;

3) the singular part of the kernel of the integral equation is explicitly integrals are approximated

by special quadrature

segregated and the

formulas, taking into account the non-analyticity

of the integrands; 4) the scheme of the calculations of the integrals near the axis additional

with the distance from the axis and in the approximation points are introduced.

Unlike existing methods, this method gives the same accuracy in the neighbourhood poles as at other points and is hundreds of times more economical of computer time.

*Zh. vjkhisl. Mat. mat. Fiz., 14, 3, 797-802,

1974.

263

of the

264

V. V. Voinov, 0. V. Vuinov and A. G. Petrov

1. The integral equation Let S be the surface of a body moving in a fluid with velocity equal to unity. The velocity potential Q, must satisfy Laplace’s equation outside the body, tend to zero at infinity and satisfy on the boundary of the body the condition acD/bh-n,,

(1.1)

where a/an is the derivative along the inward normal to the surface of the body, pfr is the projection of the normal in the direction of the velocity of the body. Using the representation of the harmonic function CP in the external domain as the sum of the potentials of a simple and a double layer [ 1,2] on the surface of the body, we can obtain an integral equation for the values of the function 4, on the surface. After solving this equation the velocity potential Q, in the entire flow region is determined by a quadrature. For a problem with axial symmetry the corresponding equation can be reduced to a one-dimensional equation with an integral along a generator of the surface of the body, which is more convenient for numerical calculations. This problem is solved by the fo~owing theorem. Theorem

Let the Lyapunov surface S be the surface of a solid of revolution with axis parallel to the velocity of motion of the body, and the function @ satisfy all the conditions of the external Neumann problem. Then~(Q) satisfies on the boundary S the equation

@i (Q,

m(Q)=

Q')n,(Q') - Q,(0') %

x(Q') dl

(Q,Q') -

L

1 x(Q)

Q’.

(1.2)

which always has a unique solution. Here L is the contour of the generator of the surface of the solid of revolution, dlQ, is the element of length of the contour at the point Q’, x(Q) is the distance from the axis to the point Q, the normal is directed into the contour, Qli(Q, Q’) is the value at the point Q’ of the potential of the homogeneous circle A(Q) which is the directrix of the surface of the sohd of revolution at the point Q:

To solve the external Neumann problem it is sufficient to determine the values of the potential on the surface, after solving the integral equation (I .2). The values of the potential Q,(P) outside the surface are given by the quadrature

(1.3)

The potential 41, is expressed in terms of a complete elliptic integral of the first kind [2] and using functional relations (see [4] , 157-l 58) may be represented in the form % (Q, Q’)=~“~~(‘i,, */zr1, z) t B=2/fr,2+u’(ri~-4xza2)],

(1.4) ri2=x2-l-y2+a2,

z=B2x2.

Short communications

Here y is the coordinate

265

of the point Q’ along the axis of revolution,

plane of the directrix circle at the point Q, F is the hypergeometric (1.4) can also be obtained

from the solution of the Laplace equation

the method of separation

of variables [5] . The parameter z is a function

coordinates.

The parameter z is convenient

z behaves asymptotically

measured from the

function,

and a=x(Q) .Equation

in toroidal coordinates

for numerical calculations

by

of one of the toroidal

because as y-+0 the quantity

of the distance r between the points Q’ and Q,

like a linear function

so that we can integrate with respect to z. The kernel of Eq. (1.2) is representable

in the form

z;=(Bz)‘h[4;;+; (+$nr)

Z+

y), and n, and ny are the components

l.The product (1 - z)G’ is bounded

function

G(z) can be represented

,

G]

(l.5)

G'=dG/dz.

G=F(‘h, ‘/a 1, z), Here r= (x-a,

(-f-f+%)

G’-

while the function

of the normal. As r-co the parameter G has a logarithmic

singularity.

The

as a series:

G(z)= 1

+c[(2;-;;!‘]z

(1.6)

zk.

k=i The radius of convergence

of the, series (1.6) equals 1. In the domain ZNI it is convenient

nG(z)=(In

to

arising from the theory of elliptic integrals [4] :

calculate G(z) by using another representation 16-m

(1-z))

G(1-z)-H(1-z),

(1.7)

k=l

It is obviously most convenient

and accurate to calculate G(z) by the polynomial

approximation

of G and H for ZE [0, 1/2],and by the use of Eq. (1.7) in the rest of the range.

2. Solution of the integral equation The contour L of the solid of revolution

is specified by the coordinates

reference points of Qi, i= 0, 1,. . . , N, arranged sequentially assumed that a differentiable

xi, yi of the

from one pole to the other. It is

curvature exists everywhere on the contour.

A parametric

representation is assume for the description of the curve. The values of the parameter T, equal to 0, 1 , . . . ,N, correspond to the reference points Qor Q1,. . . , Qd”. To determine the values of CD at the reference points the integral equation equations

(1.2) is approximated

by the system of linear

[6] N

@i

+

c

UinQn = $i,

@t=@ (Qi)

9

i=O, 1,. . . , N.

(2.1)

7l-Cl

The unknown reference values @i:

function

(1, is approximated

by interpolation

polynomials

of degree u at the

266

V. V. Vuinov. 0. V. Foinov and A. C. Petrov

where the form of the function V;(r) depends on the method of interpolation and on the degree /.LIn this paper quadratic interpolation was used (,u = 1). Then

vi=s--fT--i)~,

Tr,=‘/*(lz-il-i)

p-4,

i+q,

ze[i--l,

i+l],

i=zj,

TS[i-l,

i+l],

f=O, 1,. . . .

ls

(lz-il-2),

vi=o,

i=2p+l,

The quadrature formula

(2.3) holds for sufficiently smooth functions_@). Here KiZN, where the mesh @’ includes the reference points Qo, Ott.. . , QN and also additional points. The introduction of addjtionai points is necessary mainly because of the geometry of the problem= A measure of the variation of the function 0%(Qr Q’) for PQQ’ *u is the distance to the o-axis. At points of the basic mesh situated close to the axis, for any construction of the mesh its step will be of the order of the scale of the function (Bi, in the integral. Therefore a correct evaluation of the integrab at points of the basic mesh is impossible without the introduction of additional points. Estimates show that as the number of points N of the basic mesh increases, the total number of points necessary for the calculations increases as -N In N. The actual me~od of obta~ing the quadrature formula (2.3) will be presented below. The substitution of (2.2) into (1.2) and the application of (2.3), accurate to in~nitesimals of order ~i-r-i gives the following expression for the matrix oci,:

% C&i, c CijV,(Tj) ", ai c

rjxT(Qj'), SJ=X(Q~'), ai=x(Qt)b

TI?e right sides of (2. I) ar&=;;eterminedsimiIarly by the approximate evaluation of the integrals

The matrix Q, can also be calculated by evaluating the integrals of functions containing the weightsV,(z)as factors. The corresponding formulas can be obtained by substituting (2.1) in (2.2). However, this method is less rational, since it requires the calculation of the integrals of functions which, because of YE,vary extremely qtickiy. The components of the normal, the curvature and the derivative of the length of arc are determined by numerical differentiation with respect to the parameter T at the reference points of the curve. The corresponding parameters of the contour at points of the contour not identical with the reference points, are determined by interpolation in terms of the parameter from the interpolation nodes at the reference points.

Short communications

267

To evaluate the integrals in Eq. (1.2) it is necessary to determine G’(z) in the range o
These

functions

the functions

G(z) and

are calculated by Eqs. (1.6) and (1.7). In the range

0~2~0.5

the analytic functions G, H, G’, H’ are calculated once by series converging rapidly at the points zh=O.Olk, k=O, 1,. . . , 50, and then quadrature interpolation is used with the nodes situated at these points. In the interval 0.5
G(l-z),

H(1-z),

the functions

G’(l-z),

G and G’ are determined

H’(l-z)by

in terms

means of Eqs. (1.7). In this

parts of the integrands in Eq. (1.2) involving the factor In (I-z)are

exactly. The integrals in Eq. (2.3) were evaluated as follows. In the neighbourhood

of z = 1 or for Q’ sufficiently

close to Qi the integration

was performed with respect to z. The

kernel of Eq. (1.2), as follows from Eqs. (1.5) - (1.7), is representable In (1-z) +fz(z)r

j(z)=j1(z)

in the form

z-1,

of f2 to the integral was calculated by

where fr and fi are analytic functions.

The contribution

Simpson’s formula. For the integration

of fl (z) In (1-z) the formula (2.4)

I+,,

J x-h

m

ca

gi = In s-3

(h/x) 2n

E

2n(2n+l)

n-1

gz = 3

(2n+3) ’

Is n-i

(h/x) Zn (2n+l) (2n+3) ’

was used. If f(x) is a second-degree polynomial, (2.4) we have to put gz=O.5, g,=ln (2.4) converge extremely

Eq. (2.4) is exact. At the first step x = h and in Eq.

(2h) -5/6. At the second and subsequent

steps the series in

quickly.

The step h, in the parameter T in the approximation

of the integrals was determined

from

the condition of sufficient smallness of the variation of the argument z at this step. At the singular point z = 1 the step is smallest and it increases with distance from the initial point of integration of Qi. For z<‘/z all the integrals were evaluated by Simpson’s quadrature formula with the variable of integration 7. If z varies little between two adjacent reference points, the integration

is carried out with respect to the reference points with step 1. Close to the axis the

integration

step is automatically

chosen considerably

less than the distance between adjacent

reference points, since in this case, as indicated above, the integrand changes rapidly. Calculations

were performed

on the Minsk-22 and BESM-6 computers.

composed made it possible to calculate the potential for any sufficiently

smooth contours.

The programmes

on the surface and the associated masses

The time for the calculations

with an error of 5.10+

amounted to less than 1 set for the BESM-6 with the program in ALGOL. For comparison we mention that, for example, the method described in [3] for calculations with the same accuracy on the main part of the contour requires - 180 set on the BESM-6, and close to the axis the accuracy is considerably

reduced.

V. V. Voinov, 0. I’. Voinovand A. G. Petrov

268

TABLE t N

Approximate solution

Exact solution

0.4994 0.4771 0.4056 0.2943 0.1548 0.0000

0.5000 0.4755 0.4045 0.2939 0.1545 0.0000

0.4999 0.4909 0.4623 0.4159 0.3537 0.2778 0.1914 0.0976 0.0300

0.4904 0.4619 0.4157 0.3536 0.2778 0.1913 0.0975 0.00~0

II -I N

I

ll

17

Approximate solution 0.49996

0.49590 0.48309

0.46204 0.43309 0.39674 0.35360 0.30443 0.25003 0.19136 0.12942 0.06527 0.00000

25

0.5000

Exact solution 0.50000 0.49572 0.48296 0.46194 0.43301 0 .39668 0.35355 0.30438 0.25000 0.19134 0.12941 0.06526 0.00000

TABLE 2 Approximate solution

lY

N - li, t - 50 sec.

I

N = 17. f - 80 sec.

0.33445

: 2

0.33354 0.28029 2.9782

f

1

N - 25. t = 165 Sec.

$

0.33340 0.28010 2.9141

Exact solution

0.33333 0.28002 2.9735

Table 1 gives the results of a calculation of the potential on a sphere for various numbers of points N, specifying the contour of the solid of revolution. Table 2 gives the values of the associated masses M for various numbers of points for a sphere and for an oblate (1,==2,&=I) and a proiate (1,=1, E,=2) ellipsoid (I,, lY are the semiaxes). The computing times indicated in Table 2 apply to the Minsk-22 computer.

REFERENCES 1.

SOBOLEV, S. L., The equations of ~the~t~~al Moscow, 1966.

physics (Uravneniya matematicheskoi

2,

SRETENSKII, L. N., Theory ofthe Newtonian potential (Teoriya n’yutonovskogo potentsiala), Gostekhizdat, Moscow-Leningrad, 1946.

3.

SHEPELENKO, V. N., Calculation of the potential steady flow around a solid of revolution, in: Numerical methods of the mechanics of a continuous medium (Chislennye metody mekhaniki sploshnoi sredy), Vol. 2, No. 3, 54-63, Izd-vo SO Akad. Nauk SSSR, Novosibirsk, 1971.

4.

CASENAVE, R., integrales et fonctions el~ipti~uesen vue des appI~~ations.Ed. Centre de Documentation PArmement, 157-188, Paris, 1969.

5.

BUCHHOLZ, H., Calculationof electric and magnetic fields (Raschet elektricheskikh i magnitnykh polei), Izd-vo in. lit., Moscow, 1961.

6.

KANTOROVICH, L. V. and KRYLOV, V. I., Approximate methods of higher analysis(Priblizhennye metody vysshego analiza), Fizmatgiz, Moscow-Leningrad, 1962.

fiziki), “Nauka”,

de