- Email: [email protected]
Calculation of the steady flow of a viscous fluid past a circular cylinder
176
K. I. Babenko, N. D. Vvedenskayaand M. G. Orlova
Figure 5 for M = 3 shows the temperature graph corresponding to the tenth usual iteration. It is obvious that behind the shock wave there is some “bulge”, which slowly vanishes in subsequent iterations. Figure 6 shows the temperature profile for M= 10. It must be pointed out that in the temperature graph some excess of the equilibrium value behind the shock wave is observed. Calculations have shown that the inclusion of the conservation laws for determining the macroparameters leads to a decrease in the number of iterations required and permits an entry to be made into the range of Knudsen numbers for which the theory of the Navier-Stokes equations with sliding conditions holds. lYanslated by
J. Berry
REFERENCES 1.
WILLIS, D. R., Theoretical solutions of some problems of an almost free-molecular flow. In: Gas dynamics of rarefied gases (Gazodinamika razrezhennykh gazov), Izd-vo in. lit., Moscow, 1963, 385-400.
2.
SHAKHOV, E. M., A method of approximation of Boltzmann’s kinetic equation. In: Numerical methods in the theory of rarefiedgases (Chislennye metody v teorii razrezhennykh gazov), 84-118, VTs Akad. Nauk SSSR, Moscow, 1969.
3.
SHAKHOV, E. M., Couette’s problem for Krook’s generalized equation. The effect of a stress maximum. Zzi? Akad. Nauk SSSR. Mekhan, zhtiikostiigaza, No. 5, 16-24, 1969.
4.
ZHUK, V. I., RYKOV, V. A., and SHAKHOV, E. M., Kinetic models and the structure of a shock wave. Zzv. Akad. NaukSSSR. Mekhan. zhhikostiigaza, No. 4, 135-141, 1913.
5.
SHAKOV, E. M., Structure of a condensation jump in a monatomic gas. Izv. Akad. Nauk SSSR. Mekhan. zhidkosti igaza, No. 2, 69-75, 1969.
6.
KOGAN, M. N., Dynamics of a rarefied gas (Dinamiki razrezhennogo
7.
LYUBARSKII, G. Ya., On the kinetic theory of shock waves. Zh. eksperim. i teor. fiz., 40,4, 1050-1057, 1961.
gaza), “Nauka”, Moscow, 1967.
CALCULATION OF THE STEADY FLOW OF A VISCOUS FLUID PAST A CIRCULAR CYLINDER* K. I. BABENKO, N. D. WEDENSKAYA and M. G. ORLOVA Moscow (Received
6 July 1973)
(Revised version 24 May 1974)
A BOUNDARY value problem in a bounded region, whose solution approximates the solution of the problem is discussed. The effect of boundary conditions, specified on an additional distant boundary, is investigated. A numerical solution algorithm is described and the results of calculations for Re = l/8 to 100 are presented. 1. In this paper we consider the plane problem of the flow of a viscous fluid past an infinite cylinder. This problem has attracted the attention of many authors and has repeatedly been studied numerically, so that it has even become customary to mention the traditionality of such calculations. The majority of papers have been devoted to the circular cylinder, which is in part explained by the apparent simplicity of this case (see, for example, [l-8] ).
*Zh. vychisl.Mat. mat. Fiz., 15, 1, 183-196, 1975.
Steady flow of a viscousfluid past a circukzrcylinder
177
We first discuss the mathematical and computational questions arising here, and we then describe the numerical algorithm and the results obtained. The steady flow of a fluid past an infmite cylinder is described mathematically by the solution of the boundary value problem:
u,+
urm=coscp,
U~-+u,” =-sin
cp,
r-+00.
(1.3)
Here r, cp are polar coordinates, z+, u, are the velocity components in the radial and angular directions, 5 is a vortex, ZJis the viscosity coefficient, and C is the surface of the cylinder. We will assume that the diameter of the cylinder equals 2. Problem (1 .l)-(1.3) has not been completely solved at the present time. The existence of a solution is known (for fairly large values of V)in the class of so-called physically significant solutions [9-lo] forwhich [ (u,-u,m)2+(~P-~~m)21’~=~~-~m~ =O(r-“‘-e), r*m. h this case it is possible to write down several terms of the asymptotic expansion of u and { in the neighbourhood of r=w [11,12],itbeingthecasethatinactualfact Iu-U” I =O(r-‘h),
r-+00.
Moreover, 1U-U” ( =0 (+) everywhere outside the wake behind the body (the wake may provisionally be defined as the region where r sin’ rpO) .However, even for large values of v it is not known whether the solution is unique (see [lo] ). Instead of the dimensionless viscosity coefficient v another parameter, the Reynolds number introduced, this being defined here as 2 / Y. It must be stressed that the Re number so defined may be called the local Reynolds number, relative to the diameter C and characterizing the properties of the solution close to Cat a distance of several diameters. We also point out that the asymptotic behaviour of the solution as r+w is non-uniform with respect to Re.
Re, is often
The calculation of the solution of (1 .l)-( 1.3) is of interest because it enables us to calculate the qualitative nature of the solution of the problem, confirms the existence and stability (in some sense) of a physically significant solution, and makes it possible to obtain quantitative characteristics. For example, the variation with v of the force F with which the fluid acts on the cylinder can be established; for large values of v it is interesting to compare the function F with Lamb’s asymptotic formula (see [ 13, 141). Finally, in the case of a circular cylinder it is comparatively easy to solve the methodical questions connected with the numerical solution of a boundary value problem of the form (1 .l)-(1.3). (Th e examination of these questions is important for carrying out calculations in many more complex cases.) In the calculations the physically significant solution is sought. It is natural to seek the solution of (1 .l)-( 1.3) in some bounded domain G, specifying on an additional sufficiently distant boundary C, certain boundary conditions I&&, u&J ICI-O,
i=l,
2,
(1.4)
K. I. Babenko, N. D. Vvedenskuyaand M. G. Orlova
178
consistent with the beha~our of the solution of (1 .I)-(1.3) at infmity. It may be hoped that everywhere outside the neighbourhood of C1 the solution of such a “truncated” problem is close to the solution of the original problem. But it is first necessary to analyse the problem of the behaviour of the solutions of the truncated problems close to C’, and the effect of the boundary conditions specified on CI , on the whole domain G. This is all the more essential since in this approach to the problem a second scale appears: 2D, the diameter of G, and thereby a second Reynolds number determined by this scale. The second Reynolds number will be D times the original one. For this reason the solution in the ~ei~bourhood of C, may have a sharply expressed boundary layer, which we will call the external boundary layer. Let Cl be the circle r = D. It is comparatively simple to find appropriate boundary conditions in the part of C1 outside the wake, since outside the wake the solution of (l.l)-(1.3) is well approximated by asymptotic formulas, The difficulty is the specification of the boundary conditions on Cl within the wake, since to obtain satisfactory accuracy it is necessary to take a large number of terms of the asymptotic series. Every inaccuracy in the formulation of the boundary con~tions within the wake causes distortion of the solution in a significant part of the domain and the appearance of sharp external boundary layers of a purely parasitic nature which complicate the numerical solution of the problem. Therefore, the choice of the boundary conditions on CI is one of the important factors in developing a numerical algorithm for solving the problem. Everywhere below, the solution is assumed to be symmetric for which au,./
3-L. aq3=aw=0, ip=o,
2. We consider the following three forms of boundary conditions:
(2.1) E*==b (0,
qd--ur(D, cp)=o, cp)=o,
ZF%(D, $4-U,(D,
12-g (D, Fp)-2 (D, cp)=o,
(2.2)
za=g (D, rp)-Z(D, rp)=o.
(2.3)
Let U,, U,, 2 be identical with the first terms of the asymptotic expansion of the physically significant solution: U,(D,rp)=u,“+ U,(D, Z(D,
U~‘=cosrp-
cp) = uf+ cp)= @)=.__
U$)=
-sinqp
+
(2.4)
Fsinq
(2.5) Here F is the value of the force (in the x-direction), acting on unit length of the cylinder; we recall that the value of F depends on v and is not known in advance. We explain why in the boundary conditions the functions (2.4) have been chosen, and not their limiting values r.zYm,u,“, 0. We explain for this the nature of the external boundary layer which arises if we put F = 0 in Eqs. (2.1)-(2.4). Thereby it appears to be possible to compare conditions (2.1)-(2.3).
Steady flow ofa viscous fluid past a circular cylinder
179
Instead of Eqs. (1 .l) we consider the Navier-Stokes equations linearized for a constant flow u,=cos cp, rz,=-sin rp, c-0 (Oseen’s l~earization): 1 au 2E.$+++--=o,
i _&_~=O, au
2$++T
(2.6)
%I y ~+ff+f~)_~*s~a_;~~=o.
(
On the boundary we impose condition (1.2) and one of the conditions (2.1)-(2.3). Let v/D be fairly small, and F = 0 be substituted in (2.1)-(2.4). The principal term, determining the behaviour of the solution in the external boundary layer will be proportional to the solution decreasing as r-+0 which satisfies the corresponding one of the conditions (2.1)-(2.3), in which we have put
u, =
up,
Ii,
=
z
u!‘,
z
p),
1.
F=
(2.7)
It is easy to verify that aaIa (r / 2v) sin acp exp (r cos cp / 2v),
w3)
u.
where I, is the Bessel function with imaginary argument. For fixed (Y
Considering the function (2.7), it is sufficient in (2.8) to restrict ourselves to the first term, then for large values of r avV2
‘r
I
u-
sinqexp
GZ-
(3rry’
Lsinzg v
t
>
r-32.
[l -+-O(v/r)],
(2.9)
In order to understand how ur, h, behave we expand the solution of (2.6) in a Fourier series in tp. (The behaviour of the Fourier coefficients is also important for us, because they occur in the calculations, see below.) Let u,. =
c
u, cos x cp,
u, =
r.
, v, sin x (p,
0% sin Xcp,
c=c
x
X
1(
It is obvious that ux (r) = -
b,rx-1 + +
[--
P-l\
s-xtlti,
(s) ds + r-x-lj
sy+kox(s)dsj,
(2.10) ux (r) = bxrx-1 +
$
[F-l
1 .P+k,
0
(s) ds + F-l
5 sx +%o, (s) ds] .
ro
It follows from (2.7), (29), (2.10) that, for example, for small values of X, x>O, the principal terms (with respect to r) of uX,G and U,, V, are of the form
K. I. Babenko, N. D. Vuedenshyu and M. G. Oriova
180
ux (r) ==: II,
(r)
FS
&p-l
bKrX-1 +
f
2-;-
9
e+,
@vv,
U%P-1 2vx
Vx (r) ==: -
nr2
’
where cx=4av3 / n. In the case of the boundary conditions (2-l), where ur, up, are specified, equating to one another uX(0), U,(D) and G(D), V,(D), we obtain that the slowly decreasing term by+’ f t=D=- l/d. We similarly obtain that in the case of the boundary conditions (2.2) where ur, 5 are specified, the term b,8-* 17=Dwl f nC), and in the case of the boundary conditions (2.3) where u, t, are specified, this term is of the form b,r+’ 1~~~~3~~/~~z. Therefore, the slowly decreasing components of the coefficients h, vx have least value for boundary conditions of the form (2.3). It can be shown that for fmed r the maximum of ox (r) rs a function of x is attained for x -0 ( (r/v) I”). By comparing the values of the coefficients uX(r) , u, (r) , ox (r) for x =0 ( (r/ v)Ih) in the case of different boundary conditions (2.1~(2.3), (2.7), it also follows that the variation of these coefficients in the external boundary layer will be least for boundary conditions of the form (23). For x = 1/ o ( (r / v) “) all the Fourier coefficients are small exponentially with respect to x . The form of the functions (2.7), (29), (2.10) shows that the absolute value of the vacation of ur, u, 5 in the external boundary layer slowly decreases as D increases, and therefore the effect of the boundary conditions is felt far from D,. From a comparison of the three types of boundary conditions it is obvious that for Eqs. (2.6) conditions (2.3), (2.7) correspond to a smaller solution than the other two conditions. By taking into account subsequent terms in the asymptotic expansion of the solution of (2.6), it can be shown that conditions (2.3), (2.4) are preferable to conditions (2.1), (2.4) and conditions (2.2), (2.4). It is natural to suppose that for the solution of (l.l), (1.2), (1.4) conditions (2.3), (2.4), that is, the specification of u,, f of the form (2.4), will also be better than the other two. Numerical experiments have confirmed this supposition. Besides (2.1)-(2.3) we also consider another simple set of boundary conditions, when on part of Cl in the wake we specify not the values of the unknown functions themselves, but the values of differential operators on them. The specification on this part of the boundary, where fluid flows out of the domain, of differential operators of the unknowns is quite often used in numerical calculations (see, for example, [ 151). Such conditions are usually called soft boundary conditions. The specific nature of the case considered is the fact that Cl is a smooth curve. We will explain by a simple example what features may be introduced into the solution by a discontinuous jump of the boundary condition. Consider the problem u,S q,y-aus-buy=O, u=o,
y
=o,
y>O,
(2.11) 00,
au+@.&,;tyu,,=o,
y=o,
Z
~$0 we are givenv,=O, L’=U exp[ --r/z (ax-i-by) ],then It is easy to verify that if for y=O, the solution of (2.11) bounded at zero may be of the form u=exp [ l/z (ax+by) ]I, (r) sh ((p/2), r= (zz+$)‘h, cp=arctg y/s; lt=r”* (1+0 (r)) sin (q/2) as r-+0, that is, the derivatives of u may be unbounded at zero. Aiso, in the case where vyV=O, is specified for y=O, ~(0 the solution can only be of the form
Steady flow of a viscousfluid past a circular cylinder
u = expVh (as + WI
CL
181
Cr.) sinnrp
II
and its derivatives are unbounded at zero. Returning to the Navier-Stokes equations, we pass from I, (p to g=ln r, cp.The functions f and ru, satisfy the equations
On the boundary r=D (g=ln D) let 5 and u,, be specified for 1cp1>Oo and the values of the differential operators of r with respect to 5 and u.. for 1cp! <‘I),, . Evidently, in accordance with the example discussed above, it is not essential to specify for r=D, 1cp1<(DO the values of the first-order operators, but it is possible to specify $(riiQ) b=
&
L r+r$ru,, (D,
-$
(r-“5) 3
r $
r$
r-bg,
@oh
In the calculations on the boundary I =D the simpler boundary conditions
(2.12)
were specified, where
U,, Z
are functions of the form (2.4).
The numerical experiments showed that the specification of boundary conditions of the form (1.4), (2.12) leads to even smaller variation of the solution in the external boundary layer (and to a smaller effect of this layer), than the specification of the conditions (1.4) (2.3). 3. The solution of (1 .l), (1.2) (1.4), where conditions (1.4) are of the form (2.1)-(2.3) or (2.12), is sought by the build-up method, that is, for the corresponding boundary conditions we seek the limit as t+ 00 of the solution of the equations
(3.1)
satisfying some initial conditions (3.2)
182
K. I. Babenko,N. D. Vvedenskaya and M. G. Orlova
In its turn the solution of problem (3.1), (3.2), (1.2), (1.4) is found numerically. We discuss a description of the finite-dimensional approximation of the solution and a numerical algorithm. In the last of equations (3.1) we first pass from derivatives with respect to t to difference ratios. The unknown functions are now considered only at the discrete instants t=t”=nz, n=O, 1, . . . . (For values Re > 1 we also used in the calculations a scheme that is three-layered with respect to t, in which functions for t= (n-t l/z) T, are also introduced, see [5-71.) Denoting the unknowns specified at the instant t = P by u,*, u+,~,c”, we put (n > 0),
(3.3)
Here u,.~, h”, %” satisfy (1.2), (1.4), the value F = F on the n-th layer being calculated from cm, m
Steady flow of
a
circular cylinder
vismus jluid pasta
183
the system of functions sin kq, cos kq, k=l, 2, . . . . Here the first K eigenvalues of the integral equation (3.4) are found numerically with great accuracy and (3.4) is solved, so to say, by a natural method. We put K
U”-t-
K
c
uxncos xcp,
uqn=
x=0
c
u,% sin xv,
6”_t
x=1
wXnsin x’p. x=1
Let Eqs. (3.3) and the boundary conditions (1.2), (1.4) be satisfied on the rays k=O, 1,. . . , K.
cp=cpk=2nk! (2K+l),
Then r.~~, vwn,6~” satisfy the linear system (the subscript n is omitted for convenience) vx’+vx/r+lcuxlr-ch=O,
ux’+uxlr+wxlr=O,
(3.5) mX-rv (o~“+o~‘/F-~~~WX/~) =&, where pX=pX”l are the coefticients of the interpolation trigonometric polynomial which is identical withfh on the rays (p=rpk. The boundary conditions for r = 1 are obvious: u,(l)=v,(l)=O.
(3.6)
In the case of the boundary conditions (2.1)-(2.3) for r = D the pairs lLt=U,,
vx=
v,,
ux=
ox=sL,
ux,
vx=
vx,
ow=Qw,
respectively, are specified, where U,, V,, CL, x=0, 1,. . . , K, are the coefficients of the interpolation trigonometric polynomials for the functions (2.4). Therefore, the system (3.5) together with the boundary conditions decomposes into three equations (with their own boundary conditions) corresponding to different values of 31. For example, for conditions (2.3) we can calculate the solution as follows. We first find the solutions uX, vx, ox of the homogeneous equations (3.5) (with JL=O), corresponding to the boundary conditions &(I) =vx(D) =w,(D)
w,(l) =l.
=o,
The eigenvalues h,=l/v,( 1) of equation (3.4) are remembered. (Of course, the A, are calculated once and stored during the entire time of solving the series of boundary value problems required.) On the next layer t = tn we first find up’, v’,l’, c$ the solutions of (3.5) (with the correct pX), corresponding to the conditions
u’,” (1) = @(xl)(1) = 0,
V’,l’(D) = v,,
0:) (D) = s&.
(3.7)
It is easy to see that the required solutions uX, vx, ox are the solutions of Eqs. (3.5) corresponding to the conditions i&(l) = 0,
0 X(1) = --h XX P(1)
7
UK(II) -= vx,
wx (D) = cl,.
(3.8)
(The value of -u,t*) (1) acts as the Fourier coefficient of the right side of (3.4)). In the case of the boundary conditions (2.12) the calculation of the solution of the boundary value problem is rather more unwieldy. (The system (3.3) together with the boundary conditions does not decompose into more than three equations corresponding to different values of x.
K. I. Babenko, N. D. Vvedenskayaand M. G. Orlova
184
However, it can be carried out in two stages in each of which the problem is solved with con~tions of the form (2.3). Let (pk@D for k,
AD, fpk)=Q,
g,(f), (~bi)=O,
5, (D, qk) = hi+‘,
&p,j (D, qk) = 0,
2, . . +7 k3,
k=O, 4,. . . , L(,
j = 1 + k,, . . ., .&to, k = 0, 1, . . ., K. 2, . . . 2ko:
For the j-th version we calculate the column (bijf , i=l, L1 [h,j(D,
j=l,
cpi)l= $r7jJ-rq,j(r7
cp)Ird, v=qpi, i=i,
2, * . .) k,,
bij =
i=I
Lz ICj(D, ‘Pi-k,,)]=
1 i_ k,, . . ., 2ko.
The matrix II-‘, B= {bij}, of dimension 2k,X2ko is remembered. On the next layer t = tn we first find the solutions zz,(‘), u,(*), wx(” of Eqs. (3.9, corresponding to conditions (3.6) and the conditions
We then calculate the vectors y= {yi), i=4,2,
LI
‘yi=
L,
. . . ,2k,:
Iu’,l’CD,rpi)- u, (D, c-pi)], i = 1, 2, . . ., ko, [@)(D, (p&+)- z (D, (P%+,)], i = 1 + k,, . . ., 2k,,
and the vector q=-B-‘y. The solutions required zzX,ux, ox are the solutions of (3.5) corresponding to the boundary conditions (3.6) and to the conditions
ucp(D,0) =5 (D, 0) =o, F; (0, u,(D, (pk)=%r ~(0,
cpk)
=uv(D,
qk)
7
q%i)
==qk+ko, 5 (0,
&==I, (PA)
=Z
(D,
2,
. . . , ko,
q)k),
tWkG&
The boundary value problem for the ordinary differential equations is solved by a finitedifference method. (The solution scheme of the boundary value problems for the difference equations is an exact repetition of the solution scheme for the boundary value problems for Eqs. (3.3) described above.) Taking into account the fact that in the solution sought the derivatives with respect to I” are extremely large for r-l and may be fairly large for r-D (if LI is fairly large or the boundary conditions on C, have been unfortunately chosen), difference nets with nodes not equidistant with respect to rare used in the calculations. For this in Eqs. (3.5) a transition is made from r to a computational variable s=f (r) such that f’( 1) >>f’( D/2), f’(D) >f’( D/2>. The differential operators with respect to s are approximated by finite-differences, the difference net for s already has equidistant nodes. (The simple but unwieldy computational formulas are given in I% 61).
Steady flow of a viscousfluid past a circular cylinder
185
The following must be pointed out. It is obvious from the outset that for small values of 3t the functions uX, ux, ox, x+1, are comparatively large for r=O (1) and small for r>l. Conversely, for large values of 3c these functions are small for r=O (1) and large for r=O (WC’). In the finite-difference solution of Eqs. (3.5) it would be desirable to use for the different values of 1~ different nets and perform the calculations in different ranges of variation of r. But this would lead to the necessity for a good interpolation of the functions obtained (for calculating the right sides of (3.3)). However we did not study this question. 4. Calculations were performed for the values Re = l/8 to 100. The results are illustrated by diagrams, a commentary on them is given here. TABLE 1
Computed parameters
In the cases presented for ReG4 if r = D we have put conditions (2.3) that is, we have specified u,, %; for Re> 4 we have put conditions (2.12), that is, on one part of Cl we have specified u,, <, and on the other, differential operators of them. The distance to the external boundary D and the degree of the interpolation trigonometric polynomial K in the calculations are shown in Table 1. (We mention that the values of K for ReF9, l/8 are taken too large.) The number of points of s for Re= 100, for example, equals 75, and for Re=t’/s it equals 250. In many cases calculations were performed for different values of D and K, and also for other forms of the boundary conditions. Intermediate values of l?e were also considered. In all cases for Re>4
we took ko=6.
The initial value for Re=4 is c”=O. For Re>4 the initial value to taken was the solution corresponding to the nearest of the smaller Re computed; for Re<4 the initial value was taken as the solution corresponding to the nearest greater Re. The time step T was taken as r=1.5v for Re>2 andr=v/4 for Re<2. The dimensionless computing time after which 4-5 values are established for Re=lOO, for example, was equal to 20 (co corresponds to Re=60) ; for Re=‘/a this time equalled 300 (co corresponds to Re=l). In the transition to the computed variable we usually took s=f(r)
=(r-0.95)‘“+
@+1-r)‘“.
The value of F (and thereby the values of the boundary functions on C,) were recalculated after 40-60 computing steps. The relation between the force F and Re is shown in Fig. 1. For the cases computed the value of F was established to three places. However, the relative accuracy of F is evidently equal to 3-S% for the greatest values of Re and l-l 5% for the smallest values of Re. The behaviour of F is illustrated by the following.
K. I. Babenko, N. D. Vvedenskayaand M. G. Orlova
186
100
Re F
60
0.94
30
1.2
2
9
1.7
2.98
12.55
lia
l/4
‘is
16.83
27.83
47.35
c
F 507
2t 0
-6 b
I
J -2
-7
0
I
i
3
I
4
5
In Re
FIG. 1
FIG. 2
For small values of Re the computed values of F agree well with the asymptotic formula (see, for example, [ 141) 25%
c = 0.577...
Fas =(0.5-c+ln8-lnR_’
.
(4.1)
Judging from the results of the calculations it may be assumed that as Re decreases the difference between Fas and F tends to the constant c1~1.8 (thus, Fa,=18.64,29.67 and 49.26 for Re= r/a,*/dand l/8 respectively). The behaviour of g (1, cp)and the pressure p (rp) on the body for different values of Re is shown in Figs. 2,3. The constant term in p is so chosen that $ P(cp)Q
P
= 0.
Re = %
FIG. 4
Steady flow of a viscousjluid past a circular cylinder
187
FIG. 5
as 02 0
r
-02 -06 -7.0
FIG. 6. 1 is for ur with (p=O_62;2 is for ur with ‘p=O; 3 is for ZL~ with cp=O.ODd; 4 is for ur with q~=1.87;5 is for u,~withcp=0.62;6 is for U,with (~=I.87 Figure 4 shows graphs of the functions uI (r, 0) for rp=O for various values of Re (only up to r=50), and Fig. 5 shows graphs of u, (r,rp)for T--50. For all the values of Re (except Re=‘,/8) the wake behind the body is clearly distinguished. We notice that although the asymptote (2.4) is non-uniform with respect to Re, nevertheless for all the values of Re considered the expression (2.4) gives a good approximation for ‘u r (r, 0) for r> 100. The calculations enable us to judge the qualitative dependence on Re of the solution in the wake behind the body. Thus, for a fixed r-240 the values of zz, (r,0) decrease as Re increases for Re = 2 to 100, and conversely, decrease with Re for Re = 2 to l/8 (Figs. 4,5). The latter obviously indicates the fact that in the range r>40 for values of Re = 2 to l/8 it is possible to use formulas (2.4) and (4.1) simultaneously. A limitation for the calculation of smaller values of Re is the necessity to calculate the flow in a region whose diameter exceeds a thousand. The narrowness of the wake behind the body for Re>SO(see Figs. 4,7) compels the use of a large number of harmonics, which is a limitation for the calculation of larger values of Re. However, the fact that we have succeeded in calculating a narrow wake without especially separating it, must be regarded as an achievement of the method.
188
K. I. Babe&o, N. D, Vvedenshzyaand M. G. Orlova
6?6’ R4
R
P
FIG. 7
FIG. 8
/
r-4.7
FIG. 9 Figures 6-9 relate to the greatest value of Re considered (Re=100) and illustrate the features of the solution encountered here. It is seen (Fig. 6) how sharply u,., u, vary in the boundary layer of the body and how dependent on ‘p is the rate at which u,, u, approach u,.O”,qm. The behaviour of U, for various fixed values of r is shown in detail in Fig. 7. The behaviour of 5 for various values of rp and r is shown in Figs. 8,9. In conclusion we would like to point out some features of this paper. The method used here for the numerical solution of the problem is not the traditional finite-difference method or any of the projection methods. It appeared to us that the problem considered can be referred to the category of difficult computational problems. This was primarily due to the fact that it is required to find a solution in a fairly large domain and correctly select the boundary conditions on the distant boundary; in addition the non-standard boundary condition for the vortex on the face of the body causes additional difficulties.
Steady flow of a viscousfluid past a circular cylinder
189
We succeeded in performing calculations by the proposed method for a fairly wide range of small and moderate Reynolds numbers. Among the results presented we would like especially to point out the relation obtained for the drag as a function of the Reynolds number, and the dependence on r of the velocity component U, in the wake behind the body for large values of r. The calculations were performed on the BESM-6 computer. The graphical processing of the results was carried out by L. N. Evgrafova. At separate stages we were assisted by N. V. Barykova, M. A. Kukarkina, T. N. Piskareva. Yu. B. Radvogin assisted greatly in organizing the work. The authors thank all these colleagues. Danslated by
J. Berry
REFERENCES 1.
KELLER, H. B. and TAKAMI, N., Numerical studies of steady viscous flow about cylinders. In: Numerical solution of non-linear differentkrl equations John Wiley and Sons, Inc., New York, 1966, p. 115-140.
2.
BABENKO, K. I., WEDEMSKAYA, N. D. and ORLOVA, M. G., The steady flow of a viscous fluid past a circular cyIhider. (0 statsionamom obtekanii krugovogo tsilindra vyazkoi zhidkost’yu), Preprint IPM Akad. Nauk SSSR, 1969, No. 41.
3.
LYUL’KA, V. A. and SHCHENNIKOV, V. V., Numerical solution of the Navier-Stokes equations. In: CoIlection of theoretical papers on hydromechanics (Sb. teoretich. rabot po gidromekhan.), 107-149, VTs Akad. Nauk SSSR, Moscow, 1970.
4.
DENNIS, S. C. R. and CHANG, GAU-ZU., Numerical solutions for steady flow past a circular cylinder at Reynolds numbers up to 100. J. Fluid Me&.. 42,3,471-489, 1970.
5.
BABENKO, K. I., WEDENSKAYA, N. D. and ORLOVA, M. G., On stationary flow of viscous fluid past a circular cylinder. In: Fluid dynamics transactiona Part II, No. 5, 3749, Polish Acad. Sci., Warsaw, 1971.
6.
BABENKO, K. I, WEDENSKAYA, N. D. and ORLOVA, M. G., Results of the calculation of the flow of a viscous fluid past an infmite cylinder (Rezul’taty rascheta obtekaniya beskonechnogo tsilindra vyazkoi zhidkost’yu), Preprint IPM No. 38, Akad. Nauk SSSR, 1971.
7.
BABENKO, K. I., VVEDENSKAYA, N. D. and ORLOVA, M. G., A boundary value problem for the NavierStokes equations in the plane flow problem (Kraevaya zadacha dlya uraveneii Nav’e-Stoksa v ploskoi zadache obtekaniya), Preprint IPM No. 39, Akad. Nauk SSSR, 1971.
8.
ZANDBERGEN, P. J., The viscous flow about a circular cylinder. In: Proc Sec. Internat. Conf: Nurner. Methods in Fh&iDynamics, p. 114.119, Springer-Verl., Berlin-Heidelberg-New York, 197 1.
9.
FINN, R. and SMITH, D. R., On the stationary solutions of the Navier-Stokes Arch. Ration. Mech. and Analysis, 25, 1,26-29, 1967.
equations in two dimensions.
10. FINN, R., Stationary solutions of the Navier-Stokes equations (Statsionarnye resheniya uravnenii Nav’eStoksa), Izd-vo SO Akad. Nauk SSSR, Novosibirsk, 1964. 11. IMAI, I., On asymptotic behaviour of viscous fluid flow at a great distance from a cylindrical body, with special reference to Filton’s paradox. Proc. Roy. Sot., 208A, 1095,481.516, 195 1. 12. BABENKO, K. I., On the asymptotic behaviour of a vortex far from a body when a plane flow of a viscous fluid flows past it. PrikL matem. i mekhan., 34,5,91 l-925, 1970. (English translation: AppL Math. Me&I 13. LAMB, H., tiydrodynamics (Gidrodinamika), Gostekhizdat,
Moscow, 1947.
14. LANDAU, L. D. and LIFSHITZ, E. M., Mechanics of continuous media (Mekhanika sploshnykh sred), Gostekhizdat, Moscow, 1953. 15. ROSLYADKOV, G. S. and CHUDOV, L. A., Editors. Some applicationsof the mesh method in gas dynamics (Nekotorye primeneniya metoda setok v gazovoi dinamike), No. 3, VTs MGU, Moscow, 1971. 16. LADYZHENSKAYA, 0. A., Mathematicalproblems of the dynamics of a viwous incompressible fluid (Matematicheskie voprosy dinamiki vyazkoi neszhimaemoi zhidkosti), “Nauka”, Moscow, 1970.
190
0. M. Belotserkovskii, V. A. Gushchin and V. V. Shchennikov
17. BABENKO, K. I. and VVEDENSKAYA, N. D., The numerical solution of a boundary value problem for the Navier-Stokes equations. Zh. vychisl. Mat. mat. Fiz., 12,5, 1343-1349, 1972. 18. PAL’TSEV, B. V., The expansion of solutions of Dirichlet’s problem and a mixed problem for a biharmonic equation in a series of solutions of reducing problems. Zh. v?chisL Mat. mat. Fiz., 6, 1,43-51, 1966.
USE OF THE SPLITTING METHOD TO SOLVE PROBLEMS OF THE DYNAMICS OF A VISCOUS INCOMPRESSIBLE FLUID* 0. M. BELQTSERKOVSKII, V. A. GUSHCBIN and V. V. SHCBENNMOV Moscow (Received
19 Februav
1974)
A NUMERICAL method of solving problems of the dynamics of a viscous incompressible fluid, using the splitting of the non-stationary equations of motion, first proposed by Harlow and Chorin is considered. Unlike methods of the MAC type and its modifications @MAC, MACRL, ALE etc), the difference scheme of the proposed method enables us to calculate flows without using the boundary condition for the vortex or the pressure on the solid surface. Some results of numerical calculations of the flow past bodies of fmite dimensions are given which illustrate the possibilities of the method. 1. At the present time a fairly large number of numerical methods of solving the Navier-Stoke: equations describing the flow of an incompressible fluid are known. Most of these methods were developed for a system of equations applied to the stream function J/ and the vortex function w. A detailed analysis of methods of this type can be found in [l-3] . A common drawback of these methods is the use in some form of a boundary condition for the vortex at the solid surface, which is absent in the physical formulation of the problem. The occurrence of an additional iterative proca due to the boundary condition for the vortex at the solid surface limits the rate of convergence of the numerical algorithms. The following proposition can be stated. A difference scheme enabling the flow of a viscous incompressible fluid to be calculated without the use of a boundary condition for the vortex on the solid surface, all other conditions being equal, possesses greater efficiency. The results of calculations presented in [4] confirm the correctness of this assumption. We mention that the difference scheme of [4] satisfies the requirement indicated above. The obvious limitation of methods of solving the ($, w)-system, connected with the impossibility of extending them to the case of three-dimensional flows of a viscous fluid and to the flows of a compressible gas, explains the revived interest recently in the numerical solution of the Navier-Stokes equations, expressed in the natural variables: av;a2+(vV)~=-Vp+l/AT/, (p is the pressure,
Vis the velocity
vector,
and Y is the kinematic
V v=o
(1.1)
viscosity).
The principal difficulty in the numerical solution of the system of equations (1 .l) is due to the calculation of the pressure field. The first significant success in overcoming this difficulty was achieved by using the idea of “artificial compressibility” [5, 61. The essence of this idea is the introduction into the continuity equation of an additional term of the form (d/at) (p+ Vz/2). AS a result a modified system of equations of the form (a/81) (p+Vz/2) i-V V=O or Gp/dt-t V V-O, dV/dt+(.VV) V=-Vp+vAV, *Zh. vychisL Mat. mat. Fiz., 15, 1, 197-207, 1975.
K. I. Babenko, N. D. Vvedenskayaand M. G. Orlova
Figure 5 for M = 3 shows the temperature graph corresponding to the tenth usual iteration. It is obvious that behind the shock wave there is some “bulge”, which slowly vanishes in subsequent iterations. Figure 6 shows the temperature profile for M= 10. It must be pointed out that in the temperature graph some excess of the equilibrium value behind the shock wave is observed. Calculations have shown that the inclusion of the conservation laws for determining the macroparameters leads to a decrease in the number of iterations required and permits an entry to be made into the range of Knudsen numbers for which the theory of the Navier-Stokes equations with sliding conditions holds. lYanslated by
J. Berry
REFERENCES 1.
WILLIS, D. R., Theoretical solutions of some problems of an almost free-molecular flow. In: Gas dynamics of rarefied gases (Gazodinamika razrezhennykh gazov), Izd-vo in. lit., Moscow, 1963, 385-400.
2.
SHAKHOV, E. M., A method of approximation of Boltzmann’s kinetic equation. In: Numerical methods in the theory of rarefiedgases (Chislennye metody v teorii razrezhennykh gazov), 84-118, VTs Akad. Nauk SSSR, Moscow, 1969.
3.
SHAKHOV, E. M., Couette’s problem for Krook’s generalized equation. The effect of a stress maximum. Zzi? Akad. Nauk SSSR. Mekhan, zhtiikostiigaza, No. 5, 16-24, 1969.
4.
ZHUK, V. I., RYKOV, V. A., and SHAKHOV, E. M., Kinetic models and the structure of a shock wave. Zzv. Akad. NaukSSSR. Mekhan. zhhikostiigaza, No. 4, 135-141, 1913.
5.
SHAKOV, E. M., Structure of a condensation jump in a monatomic gas. Izv. Akad. Nauk SSSR. Mekhan. zhidkosti igaza, No. 2, 69-75, 1969.
6.
KOGAN, M. N., Dynamics of a rarefied gas (Dinamiki razrezhennogo
7.
LYUBARSKII, G. Ya., On the kinetic theory of shock waves. Zh. eksperim. i teor. fiz., 40,4, 1050-1057, 1961.
gaza), “Nauka”, Moscow, 1967.
CALCULATION OF THE STEADY FLOW OF A VISCOUS FLUID PAST A CIRCULAR CYLINDER* K. I. BABENKO, N. D. WEDENSKAYA and M. G. ORLOVA Moscow (Received
6 July 1973)
(Revised version 24 May 1974)
A BOUNDARY value problem in a bounded region, whose solution approximates the solution of the problem is discussed. The effect of boundary conditions, specified on an additional distant boundary, is investigated. A numerical solution algorithm is described and the results of calculations for Re = l/8 to 100 are presented. 1. In this paper we consider the plane problem of the flow of a viscous fluid past an infinite cylinder. This problem has attracted the attention of many authors and has repeatedly been studied numerically, so that it has even become customary to mention the traditionality of such calculations. The majority of papers have been devoted to the circular cylinder, which is in part explained by the apparent simplicity of this case (see, for example, [l-8] ).
*Zh. vychisl.Mat. mat. Fiz., 15, 1, 183-196, 1975.
Steady flow of a viscousfluid past a circukzrcylinder
177
We first discuss the mathematical and computational questions arising here, and we then describe the numerical algorithm and the results obtained. The steady flow of a fluid past an infmite cylinder is described mathematically by the solution of the boundary value problem:
u,+
urm=coscp,
U~-+u,” =-sin
cp,
r-+00.
(1.3)
Here r, cp are polar coordinates, z+, u, are the velocity components in the radial and angular directions, 5 is a vortex, ZJis the viscosity coefficient, and C is the surface of the cylinder. We will assume that the diameter of the cylinder equals 2. Problem (1 .l)-(1.3) has not been completely solved at the present time. The existence of a solution is known (for fairly large values of V)in the class of so-called physically significant solutions [9-lo] forwhich [ (u,-u,m)2+(~P-~~m)21’~=~~-~m~ =O(r-“‘-e), r*m. h this case it is possible to write down several terms of the asymptotic expansion of u and { in the neighbourhood of r=w [11,12],itbeingthecasethatinactualfact Iu-U” I =O(r-‘h),
r-+00.
Moreover, 1U-U” ( =0 (+) everywhere outside the wake behind the body (the wake may provisionally be defined as the region where r sin’ rp
Re, is often
The calculation of the solution of (1 .l)-( 1.3) is of interest because it enables us to calculate the qualitative nature of the solution of the problem, confirms the existence and stability (in some sense) of a physically significant solution, and makes it possible to obtain quantitative characteristics. For example, the variation with v of the force F with which the fluid acts on the cylinder can be established; for large values of v it is interesting to compare the function F with Lamb’s asymptotic formula (see [ 13, 141). Finally, in the case of a circular cylinder it is comparatively easy to solve the methodical questions connected with the numerical solution of a boundary value problem of the form (1 .l)-(1.3). (Th e examination of these questions is important for carrying out calculations in many more complex cases.) In the calculations the physically significant solution is sought. It is natural to seek the solution of (1 .l)-( 1.3) in some bounded domain G, specifying on an additional sufficiently distant boundary C, certain boundary conditions I&&, u&J ICI-O,
i=l,
2,
(1.4)
K. I. Babenko, N. D. Vvedenskuyaand M. G. Orlova
178
consistent with the beha~our of the solution of (1 .I)-(1.3) at infmity. It may be hoped that everywhere outside the neighbourhood of C1 the solution of such a “truncated” problem is close to the solution of the original problem. But it is first necessary to analyse the problem of the behaviour of the solutions of the truncated problems close to C’, and the effect of the boundary conditions specified on CI , on the whole domain G. This is all the more essential since in this approach to the problem a second scale appears: 2D, the diameter of G, and thereby a second Reynolds number determined by this scale. The second Reynolds number will be D times the original one. For this reason the solution in the ~ei~bourhood of C, may have a sharply expressed boundary layer, which we will call the external boundary layer. Let Cl be the circle r = D. It is comparatively simple to find appropriate boundary conditions in the part of C1 outside the wake, since outside the wake the solution of (l.l)-(1.3) is well approximated by asymptotic formulas, The difficulty is the specification of the boundary conditions on Cl within the wake, since to obtain satisfactory accuracy it is necessary to take a large number of terms of the asymptotic series. Every inaccuracy in the formulation of the boundary con~tions within the wake causes distortion of the solution in a significant part of the domain and the appearance of sharp external boundary layers of a purely parasitic nature which complicate the numerical solution of the problem. Therefore, the choice of the boundary conditions on CI is one of the important factors in developing a numerical algorithm for solving the problem. Everywhere below, the solution is assumed to be symmetric for which au,./
3-L. aq3=aw=0, ip=o,
2. We consider the following three forms of boundary conditions:
(2.1) E*==b (0,
qd--ur(D, cp)=o, cp)=o,
ZF%(D, $4-U,(D,
12-g (D, Fp)-2 (D, cp)=o,
(2.2)
za=g (D, rp)-Z(D, rp)=o.
(2.3)
Let U,, U,, 2 be identical with the first terms of the asymptotic expansion of the physically significant solution: U,(D,rp)=u,“+ U,(D, Z(D,
U~‘=cosrp-
cp) = uf+ cp)= @)=.__
U$)=
-sinqp
+
(2.4)
Fsinq
(2.5) Here F is the value of the force (in the x-direction), acting on unit length of the cylinder; we recall that the value of F depends on v and is not known in advance. We explain why in the boundary conditions the functions (2.4) have been chosen, and not their limiting values r.zYm,u,“, 0. We explain for this the nature of the external boundary layer which arises if we put F = 0 in Eqs. (2.1)-(2.4). Thereby it appears to be possible to compare conditions (2.1)-(2.3).
Steady flow ofa viscous fluid past a circular cylinder
179
Instead of Eqs. (1 .l) we consider the Navier-Stokes equations linearized for a constant flow u,=cos cp, rz,=-sin rp, c-0 (Oseen’s l~earization): 1 au 2E.$+++--=o,
i _&_~=O, au
2$++T
(2.6)
%I y ~+ff+f~)_~*s~a_;~~=o.
(
On the boundary we impose condition (1.2) and one of the conditions (2.1)-(2.3). Let v/D be fairly small, and F = 0 be substituted in (2.1)-(2.4). The principal term, determining the behaviour of the solution in the external boundary layer will be proportional to the solution decreasing as r-+0 which satisfies the corresponding one of the conditions (2.1)-(2.3), in which we have put
u, =
up,
Ii,
=
z
u!‘,
z
p),
1.
F=
(2.7)
It is easy to verify that aaIa (r / 2v) sin acp exp (r cos cp / 2v),
w3)
u.
where I, is the Bessel function with imaginary argument. For fixed (Y
Considering the function (2.7), it is sufficient in (2.8) to restrict ourselves to the first term, then for large values of r avV2
‘r
I
u-
sinqexp
GZ-
(3rry’
Lsinzg v
t
>
r-32.
[l -+-O(v/r)],
(2.9)
In order to understand how ur, h, behave we expand the solution of (2.6) in a Fourier series in tp. (The behaviour of the Fourier coefficients is also important for us, because they occur in the calculations, see below.) Let u,. =
c
u, cos x cp,
u, =
r.
, v, sin x (p,
0% sin Xcp,
c=c
x
X
1(
It is obvious that ux (r) = -
b,rx-1 + +
[--
P-l\
s-xtlti,
(s) ds + r-x-lj
sy+kox(s)dsj,
(2.10) ux (r) = bxrx-1 +
$
[F-l
1 .P+k,
0
(s) ds + F-l
5 sx +%o, (s) ds] .
ro
It follows from (2.7), (29), (2.10) that, for example, for small values of X, x>O, the principal terms (with respect to r) of uX,G and U,, V, are of the form
K. I. Babenko, N. D. Vuedenshyu and M. G. Oriova
180
ux (r) ==: II,
(r)
FS
&p-l
bKrX-1 +
f
2-;-
9
e+,
@vv,
U%P-1 2vx
Vx (r) ==: -
nr2
’
where cx=4av3 / n. In the case of the boundary conditions (2-l), where ur, up, are specified, equating to one another uX(0), U,(D) and G(D), V,(D), we obtain that the slowly decreasing term by+’ f t=D=- l/d. We similarly obtain that in the case of the boundary conditions (2.2) where ur, 5 are specified, the term b,8-* 17=Dwl f nC), and in the case of the boundary conditions (2.3) where u, t, are specified, this term is of the form b,r+’ 1~~~~3~~/~~z. Therefore, the slowly decreasing components of the coefficients h, vx have least value for boundary conditions of the form (2.3). It can be shown that for fmed r the maximum of ox (r) rs a function of x is attained for x -0 ( (r/v) I”). By comparing the values of the coefficients uX(r) , u, (r) , ox (r) for x =0 ( (r/ v)Ih) in the case of different boundary conditions (2.1~(2.3), (2.7), it also follows that the variation of these coefficients in the external boundary layer will be least for boundary conditions of the form (23). For x = 1/ o ( (r / v) “) all the Fourier coefficients are small exponentially with respect to x . The form of the functions (2.7), (29), (2.10) shows that the absolute value of the vacation of ur, u, 5 in the external boundary layer slowly decreases as D increases, and therefore the effect of the boundary conditions is felt far from D,. From a comparison of the three types of boundary conditions it is obvious that for Eqs. (2.6) conditions (2.3), (2.7) correspond to a smaller solution than the other two conditions. By taking into account subsequent terms in the asymptotic expansion of the solution of (2.6), it can be shown that conditions (2.3), (2.4) are preferable to conditions (2.1), (2.4) and conditions (2.2), (2.4). It is natural to suppose that for the solution of (l.l), (1.2), (1.4) conditions (2.3), (2.4), that is, the specification of u,, f of the form (2.4), will also be better than the other two. Numerical experiments have confirmed this supposition. Besides (2.1)-(2.3) we also consider another simple set of boundary conditions, when on part of Cl in the wake we specify not the values of the unknown functions themselves, but the values of differential operators on them. The specification on this part of the boundary, where fluid flows out of the domain, of differential operators of the unknowns is quite often used in numerical calculations (see, for example, [ 151). Such conditions are usually called soft boundary conditions. The specific nature of the case considered is the fact that Cl is a smooth curve. We will explain by a simple example what features may be introduced into the solution by a discontinuous jump of the boundary condition. Consider the problem u,S q,y-aus-buy=O, u=o,
y
=o,
y>O,
(2.11) 00,
au+@.&,;tyu,,=o,
y=o,
Z
~$0 we are givenv,=O, L’=U exp[ --r/z (ax-i-by) ],then It is easy to verify that if for y=O, the solution of (2.11) bounded at zero may be of the form u=exp [ l/z (ax+by) ]I, (r) sh ((p/2), r= (zz+$)‘h, cp=arctg y/s; lt=r”* (1+0 (r)) sin (q/2) as r-+0, that is, the derivatives of u may be unbounded at zero. Aiso, in the case where vyV=O, is specified for y=O, ~(0 the solution can only be of the form
Steady flow of a viscousfluid past a circular cylinder
u = expVh (as + WI
CL
181
Cr.) sinnrp
II
and its derivatives are unbounded at zero. Returning to the Navier-Stokes equations, we pass from I, (p to g=ln r, cp.The functions f and ru, satisfy the equations
On the boundary r=D (g=ln D) let 5 and u,, be specified for 1cp1>Oo and the values of the differential operators of r with respect to 5 and u.. for 1cp! <‘I),, . Evidently, in accordance with the example discussed above, it is not essential to specify for r=D, 1cp1<(DO the values of the first-order operators, but it is possible to specify $(riiQ) b=
&
L r+r$ru,, (D,
-$
(r-“5) 3
r $
r$
r-bg,
@oh
In the calculations on the boundary I =D the simpler boundary conditions
(2.12)
were specified, where
U,, Z
are functions of the form (2.4).
The numerical experiments showed that the specification of boundary conditions of the form (1.4), (2.12) leads to even smaller variation of the solution in the external boundary layer (and to a smaller effect of this layer), than the specification of the conditions (1.4) (2.3). 3. The solution of (1 .l), (1.2) (1.4), where conditions (1.4) are of the form (2.1)-(2.3) or (2.12), is sought by the build-up method, that is, for the corresponding boundary conditions we seek the limit as t+ 00 of the solution of the equations
(3.1)
satisfying some initial conditions (3.2)
182
K. I. Babenko,N. D. Vvedenskaya and M. G. Orlova
In its turn the solution of problem (3.1), (3.2), (1.2), (1.4) is found numerically. We discuss a description of the finite-dimensional approximation of the solution and a numerical algorithm. In the last of equations (3.1) we first pass from derivatives with respect to t to difference ratios. The unknown functions are now considered only at the discrete instants t=t”=nz, n=O, 1, . . . . (For values Re > 1 we also used in the calculations a scheme that is three-layered with respect to t, in which functions for t= (n-t l/z) T, are also introduced, see [5-71.) Denoting the unknowns specified at the instant t = P by u,*, u+,~,c”, we put (n > 0),
(3.3)
Here u,.~, h”, %” satisfy (1.2), (1.4), the value F = F on the n-th layer being calculated from cm, m
Steady flow of
a
circular cylinder
vismus jluid pasta
183
the system of functions sin kq, cos kq, k=l, 2, . . . . Here the first K eigenvalues of the integral equation (3.4) are found numerically with great accuracy and (3.4) is solved, so to say, by a natural method. We put K
U”-t-
K
c
uxncos xcp,
uqn=
x=0
c
u,% sin xv,
6”_t
x=1
wXnsin x’p. x=1
Let Eqs. (3.3) and the boundary conditions (1.2), (1.4) be satisfied on the rays k=O, 1,. . . , K.
cp=cpk=2nk! (2K+l),
Then r.~~, vwn,6~” satisfy the linear system (the subscript n is omitted for convenience) vx’+vx/r+lcuxlr-ch=O,
ux’+uxlr+wxlr=O,
(3.5) mX-rv (o~“+o~‘/F-~~~WX/~) =&, where pX=pX”l are the coefticients of the interpolation trigonometric polynomial which is identical withfh on the rays (p=rpk. The boundary conditions for r = 1 are obvious: u,(l)=v,(l)=O.
(3.6)
In the case of the boundary conditions (2.1)-(2.3) for r = D the pairs lLt=U,,
vx=
v,,
ux=
ox=sL,
ux,
vx=
vx,
ow=Qw,
respectively, are specified, where U,, V,, CL, x=0, 1,. . . , K, are the coefficients of the interpolation trigonometric polynomials for the functions (2.4). Therefore, the system (3.5) together with the boundary conditions decomposes into three equations (with their own boundary conditions) corresponding to different values of 31. For example, for conditions (2.3) we can calculate the solution as follows. We first find the solutions uX, vx, ox of the homogeneous equations (3.5) (with JL=O), corresponding to the boundary conditions &(I) =vx(D) =w,(D)
w,(l) =l.
=o,
The eigenvalues h,=l/v,( 1) of equation (3.4) are remembered. (Of course, the A, are calculated once and stored during the entire time of solving the series of boundary value problems required.) On the next layer t = tn we first find up’, v’,l’, c$ the solutions of (3.5) (with the correct pX), corresponding to the conditions
u’,” (1) = @(xl)(1) = 0,
V’,l’(D) = v,,
0:) (D) = s&.
(3.7)
It is easy to see that the required solutions uX, vx, ox are the solutions of Eqs. (3.5) corresponding to the conditions i&(l) = 0,
0 X(1) = --h XX P(1)
7
UK(II) -= vx,
wx (D) = cl,.
(3.8)
(The value of -u,t*) (1) acts as the Fourier coefficient of the right side of (3.4)). In the case of the boundary conditions (2.12) the calculation of the solution of the boundary value problem is rather more unwieldy. (The system (3.3) together with the boundary conditions does not decompose into more than three equations corresponding to different values of x.
K. I. Babenko, N. D. Vvedenskayaand M. G. Orlova
184
However, it can be carried out in two stages in each of which the problem is solved with con~tions of the form (2.3). Let (pk
AD, fpk)=Q,
g,(f), (~bi)=O,
5, (D, qk) = hi+‘,
&p,j (D, qk) = 0,
2, . . +7 k3,
k=O, 4,. . . , L(,
j = 1 + k,, . . ., .&to, k = 0, 1, . . ., K. 2, . . . 2ko:
For the j-th version we calculate the column (bijf , i=l, L1 [h,j(D,
j=l,
cpi)l= $r7jJ-rq,j(r7
cp)Ird, v=qpi, i=i,
2, * . .) k,,
bij =
i=I
Lz ICj(D, ‘Pi-k,,)]=
1 i_ k,, . . ., 2ko.
The matrix II-‘, B= {bij}, of dimension 2k,X2ko is remembered. On the next layer t = tn we first find the solutions zz,(‘), u,(*), wx(” of Eqs. (3.9, corresponding to conditions (3.6) and the conditions
We then calculate the vectors y= {yi), i=4,2,
LI
‘yi=
L,
. . . ,2k,:
Iu’,l’CD,rpi)- u, (D, c-pi)], i = 1, 2, . . ., ko, [@)(D, (p&+)- z (D, (P%+,)], i = 1 + k,, . . ., 2k,,
and the vector q=-B-‘y. The solutions required zzX,ux, ox are the solutions of (3.5) corresponding to the boundary conditions (3.6) and to the conditions
ucp(D,0) =5 (D, 0) =o, F; (0, u,(D, (pk)=%r ~(0,
cpk)
=uv(D,
qk)
7
q%i)
==qk+ko, 5 (0,
&==I, (PA)
=Z
(D,
2,
. . . , ko,
q)k),
tWkG&
The boundary value problem for the ordinary differential equations is solved by a finitedifference method. (The solution scheme of the boundary value problems for the difference equations is an exact repetition of the solution scheme for the boundary value problems for Eqs. (3.3) described above.) Taking into account the fact that in the solution sought the derivatives with respect to I” are extremely large for r-l and may be fairly large for r-D (if LI is fairly large or the boundary conditions on C, have been unfortunately chosen), difference nets with nodes not equidistant with respect to rare used in the calculations. For this in Eqs. (3.5) a transition is made from r to a computational variable s=f (r) such that f’( 1) >>f’( D/2), f’(D) >f’( D/2>. The differential operators with respect to s are approximated by finite-differences, the difference net for s already has equidistant nodes. (The simple but unwieldy computational formulas are given in I% 61).
Steady flow of a viscousfluid past a circular cylinder
185
The following must be pointed out. It is obvious from the outset that for small values of 3t the functions uX, ux, ox, x+1, are comparatively large for r=O (1) and small for r>l. Conversely, for large values of 3c these functions are small for r=O (1) and large for r=O (WC’). In the finite-difference solution of Eqs. (3.5) it would be desirable to use for the different values of 1~ different nets and perform the calculations in different ranges of variation of r. But this would lead to the necessity for a good interpolation of the functions obtained (for calculating the right sides of (3.3)). However we did not study this question. 4. Calculations were performed for the values Re = l/8 to 100. The results are illustrated by diagrams, a commentary on them is given here. TABLE 1
Computed parameters
In the cases presented for ReG4 if r = D we have put conditions (2.3) that is, we have specified u,, %; for Re> 4 we have put conditions (2.12), that is, on one part of Cl we have specified u,, <, and on the other, differential operators of them. The distance to the external boundary D and the degree of the interpolation trigonometric polynomial K in the calculations are shown in Table 1. (We mention that the values of K for ReF9, l/8 are taken too large.) The number of points of s for Re= 100, for example, equals 75, and for Re=t’/s it equals 250. In many cases calculations were performed for different values of D and K, and also for other forms of the boundary conditions. Intermediate values of l?e were also considered. In all cases for Re>4
we took ko=6.
The initial value for Re=4 is c”=O. For Re>4 the initial value to taken was the solution corresponding to the nearest of the smaller Re computed; for Re<4 the initial value was taken as the solution corresponding to the nearest greater Re. The time step T was taken as r=1.5v for Re>2 andr=v/4 for Re<2. The dimensionless computing time after which 4-5 values are established for Re=lOO, for example, was equal to 20 (co corresponds to Re=60) ; for Re=‘/a this time equalled 300 (co corresponds to Re=l). In the transition to the computed variable we usually took s=f(r)
=(r-0.95)‘“+
@+1-r)‘“.
The value of F (and thereby the values of the boundary functions on C,) were recalculated after 40-60 computing steps. The relation between the force F and Re is shown in Fig. 1. For the cases computed the value of F was established to three places. However, the relative accuracy of F is evidently equal to 3-S% for the greatest values of Re and l-l 5% for the smallest values of Re. The behaviour of F is illustrated by the following.
K. I. Babenko, N. D. Vvedenskayaand M. G. Orlova
186
100
Re F
60
0.94
30
1.2
2
9
1.7
2.98
12.55
lia
l/4
‘is
16.83
27.83
47.35
c
F 507
2t 0
-6 b
I
J -2
-7
0
I
i
3
I
4
5
In Re
FIG. 1
FIG. 2
For small values of Re the computed values of F agree well with the asymptotic formula (see, for example, [ 141) 25%
c = 0.577...
Fas =(0.5-c+ln8-lnR_’
.
(4.1)
Judging from the results of the calculations it may be assumed that as Re decreases the difference between Fas and F tends to the constant c1~1.8 (thus, Fa,=18.64,29.67 and 49.26 for Re= r/a,*/dand l/8 respectively). The behaviour of g (1, cp)and the pressure p (rp) on the body for different values of Re is shown in Figs. 2,3. The constant term in p is so chosen that $ P(cp)Q
P
= 0.
Re = %
FIG. 4
Steady flow of a viscousjluid past a circular cylinder
187
FIG. 5
as 02 0
r
-02 -06 -7.0
FIG. 6. 1 is for ur with (p=O_62;2 is for ur with ‘p=O; 3 is for ZL~ with cp=O.ODd; 4 is for ur with q~=1.87;5 is for u,~withcp=0.62;6 is for U,with (~=I.87 Figure 4 shows graphs of the functions uI (r, 0) for rp=O for various values of Re (only up to r=50), and Fig. 5 shows graphs of u, (r,rp)for T--50. For all the values of Re (except Re=‘,/8) the wake behind the body is clearly distinguished. We notice that although the asymptote (2.4) is non-uniform with respect to Re, nevertheless for all the values of Re considered the expression (2.4) gives a good approximation for ‘u r (r, 0) for r> 100. The calculations enable us to judge the qualitative dependence on Re of the solution in the wake behind the body. Thus, for a fixed r-240 the values of zz, (r,0) decrease as Re increases for Re = 2 to 100, and conversely, decrease with Re for Re = 2 to l/8 (Figs. 4,5). The latter obviously indicates the fact that in the range r>40 for values of Re = 2 to l/8 it is possible to use formulas (2.4) and (4.1) simultaneously. A limitation for the calculation of smaller values of Re is the necessity to calculate the flow in a region whose diameter exceeds a thousand. The narrowness of the wake behind the body for Re>SO(see Figs. 4,7) compels the use of a large number of harmonics, which is a limitation for the calculation of larger values of Re. However, the fact that we have succeeded in calculating a narrow wake without especially separating it, must be regarded as an achievement of the method.
188
K. I. Babe&o, N. D, Vvedenshzyaand M. G. Orlova
6?6’ R4
R
P
FIG. 7
FIG. 8
/
r-4.7
FIG. 9 Figures 6-9 relate to the greatest value of Re considered (Re=100) and illustrate the features of the solution encountered here. It is seen (Fig. 6) how sharply u,., u, vary in the boundary layer of the body and how dependent on ‘p is the rate at which u,, u, approach u,.O”,qm. The behaviour of U, for various fixed values of r is shown in detail in Fig. 7. The behaviour of 5 for various values of rp and r is shown in Figs. 8,9. In conclusion we would like to point out some features of this paper. The method used here for the numerical solution of the problem is not the traditional finite-difference method or any of the projection methods. It appeared to us that the problem considered can be referred to the category of difficult computational problems. This was primarily due to the fact that it is required to find a solution in a fairly large domain and correctly select the boundary conditions on the distant boundary; in addition the non-standard boundary condition for the vortex on the face of the body causes additional difficulties.
Steady flow of a viscousfluid past a circular cylinder
189
We succeeded in performing calculations by the proposed method for a fairly wide range of small and moderate Reynolds numbers. Among the results presented we would like especially to point out the relation obtained for the drag as a function of the Reynolds number, and the dependence on r of the velocity component U, in the wake behind the body for large values of r. The calculations were performed on the BESM-6 computer. The graphical processing of the results was carried out by L. N. Evgrafova. At separate stages we were assisted by N. V. Barykova, M. A. Kukarkina, T. N. Piskareva. Yu. B. Radvogin assisted greatly in organizing the work. The authors thank all these colleagues. Danslated by
J. Berry
REFERENCES 1.
KELLER, H. B. and TAKAMI, N., Numerical studies of steady viscous flow about cylinders. In: Numerical solution of non-linear differentkrl equations John Wiley and Sons, Inc., New York, 1966, p. 115-140.
2.
BABENKO, K. I., WEDEMSKAYA, N. D. and ORLOVA, M. G., The steady flow of a viscous fluid past a circular cyIhider. (0 statsionamom obtekanii krugovogo tsilindra vyazkoi zhidkost’yu), Preprint IPM Akad. Nauk SSSR, 1969, No. 41.
3.
LYUL’KA, V. A. and SHCHENNIKOV, V. V., Numerical solution of the Navier-Stokes equations. In: CoIlection of theoretical papers on hydromechanics (Sb. teoretich. rabot po gidromekhan.), 107-149, VTs Akad. Nauk SSSR, Moscow, 1970.
4.
DENNIS, S. C. R. and CHANG, GAU-ZU., Numerical solutions for steady flow past a circular cylinder at Reynolds numbers up to 100. J. Fluid Me&.. 42,3,471-489, 1970.
5.
BABENKO, K. I., WEDENSKAYA, N. D. and ORLOVA, M. G., On stationary flow of viscous fluid past a circular cylinder. In: Fluid dynamics transactiona Part II, No. 5, 3749, Polish Acad. Sci., Warsaw, 1971.
6.
BABENKO, K. I, WEDENSKAYA, N. D. and ORLOVA, M. G., Results of the calculation of the flow of a viscous fluid past an infmite cylinder (Rezul’taty rascheta obtekaniya beskonechnogo tsilindra vyazkoi zhidkost’yu), Preprint IPM No. 38, Akad. Nauk SSSR, 1971.
7.
BABENKO, K. I., VVEDENSKAYA, N. D. and ORLOVA, M. G., A boundary value problem for the NavierStokes equations in the plane flow problem (Kraevaya zadacha dlya uraveneii Nav’e-Stoksa v ploskoi zadache obtekaniya), Preprint IPM No. 39, Akad. Nauk SSSR, 1971.
8.
ZANDBERGEN, P. J., The viscous flow about a circular cylinder. In: Proc Sec. Internat. Conf: Nurner. Methods in Fh&iDynamics, p. 114.119, Springer-Verl., Berlin-Heidelberg-New York, 197 1.
9.
FINN, R. and SMITH, D. R., On the stationary solutions of the Navier-Stokes Arch. Ration. Mech. and Analysis, 25, 1,26-29, 1967.
equations in two dimensions.
10. FINN, R., Stationary solutions of the Navier-Stokes equations (Statsionarnye resheniya uravnenii Nav’eStoksa), Izd-vo SO Akad. Nauk SSSR, Novosibirsk, 1964. 11. IMAI, I., On asymptotic behaviour of viscous fluid flow at a great distance from a cylindrical body, with special reference to Filton’s paradox. Proc. Roy. Sot., 208A, 1095,481.516, 195 1. 12. BABENKO, K. I., On the asymptotic behaviour of a vortex far from a body when a plane flow of a viscous fluid flows past it. PrikL matem. i mekhan., 34,5,91 l-925, 1970. (English translation: AppL Math. Me&I 13. LAMB, H., tiydrodynamics (Gidrodinamika), Gostekhizdat,
Moscow, 1947.
14. LANDAU, L. D. and LIFSHITZ, E. M., Mechanics of continuous media (Mekhanika sploshnykh sred), Gostekhizdat, Moscow, 1953. 15. ROSLYADKOV, G. S. and CHUDOV, L. A., Editors. Some applicationsof the mesh method in gas dynamics (Nekotorye primeneniya metoda setok v gazovoi dinamike), No. 3, VTs MGU, Moscow, 1971. 16. LADYZHENSKAYA, 0. A., Mathematicalproblems of the dynamics of a viwous incompressible fluid (Matematicheskie voprosy dinamiki vyazkoi neszhimaemoi zhidkosti), “Nauka”, Moscow, 1970.
190
0. M. Belotserkovskii, V. A. Gushchin and V. V. Shchennikov
17. BABENKO, K. I. and VVEDENSKAYA, N. D., The numerical solution of a boundary value problem for the Navier-Stokes equations. Zh. vychisl. Mat. mat. Fiz., 12,5, 1343-1349, 1972. 18. PAL’TSEV, B. V., The expansion of solutions of Dirichlet’s problem and a mixed problem for a biharmonic equation in a series of solutions of reducing problems. Zh. v?chisL Mat. mat. Fiz., 6, 1,43-51, 1966.
USE OF THE SPLITTING METHOD TO SOLVE PROBLEMS OF THE DYNAMICS OF A VISCOUS INCOMPRESSIBLE FLUID* 0. M. BELQTSERKOVSKII, V. A. GUSHCBIN and V. V. SHCBENNMOV Moscow (Received
19 Februav
1974)
A NUMERICAL method of solving problems of the dynamics of a viscous incompressible fluid, using the splitting of the non-stationary equations of motion, first proposed by Harlow and Chorin is considered. Unlike methods of the MAC type and its modifications @MAC, MACRL, ALE etc), the difference scheme of the proposed method enables us to calculate flows without using the boundary condition for the vortex or the pressure on the solid surface. Some results of numerical calculations of the flow past bodies of fmite dimensions are given which illustrate the possibilities of the method. 1. At the present time a fairly large number of numerical methods of solving the Navier-Stoke: equations describing the flow of an incompressible fluid are known. Most of these methods were developed for a system of equations applied to the stream function J/ and the vortex function w. A detailed analysis of methods of this type can be found in [l-3] . A common drawback of these methods is the use in some form of a boundary condition for the vortex at the solid surface, which is absent in the physical formulation of the problem. The occurrence of an additional iterative proca due to the boundary condition for the vortex at the solid surface limits the rate of convergence of the numerical algorithms. The following proposition can be stated. A difference scheme enabling the flow of a viscous incompressible fluid to be calculated without the use of a boundary condition for the vortex on the solid surface, all other conditions being equal, possesses greater efficiency. The results of calculations presented in [4] confirm the correctness of this assumption. We mention that the difference scheme of [4] satisfies the requirement indicated above. The obvious limitation of methods of solving the ($, w)-system, connected with the impossibility of extending them to the case of three-dimensional flows of a viscous fluid and to the flows of a compressible gas, explains the revived interest recently in the numerical solution of the Navier-Stokes equations, expressed in the natural variables: av;a2+(vV)~=-Vp+l/AT/, (p is the pressure,
Vis the velocity
vector,
and Y is the kinematic
V v=o
(1.1)
viscosity).
The principal difficulty in the numerical solution of the system of equations (1 .l) is due to the calculation of the pressure field. The first significant success in overcoming this difficulty was achieved by using the idea of “artificial compressibility” [5, 61. The essence of this idea is the introduction into the continuity equation of an additional term of the form (d/at) (p+ Vz/2). AS a result a modified system of equations of the form (a/81) (p+Vz/2) i-V V=O or Gp/dt-t V V-O, dV/dt+(.VV) V=-Vp+vAV, *Zh. vychisL Mat. mat. Fiz., 15, 1, 197-207, 1975.