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.
Problems ofthe dynamics ofa viscous in~mpre~~ble j&id
191
is obtained which can be solved by the use of different versions of the splitting method [7,8]. As an example we mention papers [9, lo]. Unlike [9] where the classical scheme of formal splitting into one-dimensional operators is used, in [IO] splitting mainly based on the PIC method was used [l l] . In this scheme we first consider the field from the equation
F---Vv” ‘G
+ (If?)
V” = YAP”
(the field of Y” is assumed to be known). Then the field of Vn+l is determined by the formula ‘6f”+‘=F%vp,
where p is the steady-state solution of the equation i3pidt-k QV=~hpa
A similar method of calculating the pressure is used in the MAC method 1121. As the boundary condition on the solid surface these approaches use the equation of motion projected onto the normal to the surface at the boundary points. The latter fact reduces the efficiency of the numerical methods, since these conditions are not present in the physical formulation of the problem. In the SMAC method [ 131 because of the successive use of the splitting idea [l 1, 121 and the replacement of the pressure p by the potential function q, ensuring that the field of Vn+l is solenoidal and satisfying homogeneous boundary conditions, it has been possible to avoid the difficulties of the MAC method indicated. An original modi~cation of the boundary conditions of the MAC method, enabling homogeneous boundary conditions to be obtained for the pressure, is presented in [14]. It must be pointed out that in the SMAC method and the modified MAC method 1141, because of the difference schemes chosen, the satisfaction of the adhesion condition leads to the need to determine the boundary value of the vortex on the solid surface, satisfying condition [2], The latter is known to be a condition of the first order of accuracy with respect to the step of the spatial net. Moreover, the adhesion condition in the SMAC method does not guarantee the balance of forces on the solid surface. The error here is of order O(V). 2. We will dwell briefly on a desc~ption of the difference scheme of the proposed method and the main stages of the calculation. As initial equations we use the system of equations (1 .l). We consider the following splitting scheme, similar to [ 131:
f-vn = -- (j’“Q z Ap=QBlz,
V” + yA’v”
(2.la)
Vp+‘=p”‘=o 7
(2Sb) (2.lc)
192
0. M. Belotserkovskii, V. A. Gushchin and V. V. Shchennikov
0
ui,-72
%+1/ 7
’
e %-!h,@” A m ,,,,,,,,,,,,,,,,
m
0
V. I
-‘/I
%+
,,,,,,,,,
‘/r, 0 ///,,,I
B
FIG. 1 For the case of a Cartesian coordinate system and a uniform net (Fig. 1) the two-dimensional difference scheme has the following form:
(2.24
ntl Di,j = 0;
(2.2b)
193
Problems of the dynamics of a viscous incompressible jluid
n+1 Uit*/,,j
=
&tljs3j
-
+
(Pit1.j-
Pi,j)~
zi,j+l/r
-
f
(Pi,j+l
P&i)-
(2.24
nt1 Vi,j+l/, =
-
It is easy to see that the scheme (2.2) approximates (2.1) with an error of order 0 (z, h2), where h=max (h,, h,). The essential feature of the proposed method is the choice of the boundary conditions. Since we are interested in solving problems of the flow of a viscous incompressible fluid past a body of finite dimensions, two fundamental types of boundary conditions may be distinguished: conditions on the solid surface and conditions on lines fairly far from the body flowed round. We will discuss each of these conditions. The boundary conditions on the solid surface (see Fig. 1): (the nonflow condition)
$, - ‘/z= 0
(2.3) $+ %/2sy = u:-
%,N (the adhesion condition)
from the latter it follows that n
Uitl,‘,.O =
-
UitV,. 0
2
n
-+ -
Uitl,i. 1
6
+ 0 (hz3).
(2.4)
Condition (2.4) enables us to avoid the necessity to introduce a layer of fictitious cells (within the solid body), which in schemes of the MAC, SMAC and modified MAC type [14] have led to the implicit calculation of the vortex value on the solid surface with first-order accuracy. We notice that within the scope of the proposed approach it is not required to calculate the value of the vortex on the solid surface. The latter may be found from the calculated velocity field by using any of the difference representations of the expression for the vortex
at the boundary points. The boundary conditions on a line distant from the body flowed round (EF, CD, GH), are the conditions in the unperturbed flow, which in the case U, ((OX are of the form (for example, on CD) V&*/, = 0,
u&,,.v
= ucc.
That homogeneous boundary conditions will be obtained for the calculation of the pressure field is guaranteed by the approach of [20], which consists of the following. Putting v?:,; = 0 (for the case of a solid surface) and VT;,,,,, = 0 (in the case of a line distant from the body), from (2.2~) we have F&A,, =
*(Pi,0
-
Pi,-I),
zi, N+*/, =
2
(pi, N-+1-
Pi,N)*
(2.5)
Taking into account (2.5) the difference equation for calculating the pressure in the boundary cells assumes the form
0. M. Belotserkovskii, K A. Gushchin and K V. Shchennikov
194
Pi-Lo~~Pitl.0
+ z
Pi,l-t2Pi,O
_
a;,,)
r
(2.6)
pip0= (2z,/hr2 : 2z/h,2) ( To where ai,,
=
@+li,,
0 ;
gi-a/,,
0
I “;11/, ,
h2 z0=Sy.
1
2
The relation for pi, N is written similarly. The numerical realization of the proposed method includes the following three stages. Stage I. At this stage the values of u”i+l/,, j, v”i,j+j/%are determined from Eqs. (2.2a) by using
conditions (2.3) (after eliminating the boundary values Ej+l,,,O ). By calculating the values of n j at the internal nodes of the domain we find u”i+~l,, ui+,,,, o from (2.4). Stage II. The pressure field is determined from (2.2b). The relations (2.6) are used in the boundary cells. Stage III. The values of &&, v?$,,~ are found from (2.26). As the initial approximations u?j+li, we choose arbitrary functions satisfying the required boundary conditions.
for r&,,j,
The steady-state solution of the system of equations (2.2) is obtained by repeating the stages indicated above until the following build-up criterion is satisfied:
We remark that it is advisable to use the upper relaxation method [15] to solve Eqs. (2.2b) and (2.6). 3. The stability can be investigated stage-by-stage. A stability criterion for the first stage can be obtained from the first differential approximation [ 161. For equations (2.2a) the first differential approximation has the form
The stability criterion of the difference sch;y a<
used follows from (3.1) (h,=h,)
:
u” + Y2 *
Eliminating p from (2.2b) and (2.2c), the unconditional stability of the second and third stages is easily demonstrated by Fourier’s method. 4. Solutions of a whole series of problems of external hydrodynamics were obtained by the method described. The flows round various bodies of finite dimensions (a rectangular beam, a cylinder of finite length with axis parallel to the velocity vector of the incident flow U,, a sphere, a cylinder with axis perpendicular to U,, and also bodies of more complex form) were studied
Problems of the dynamics of a viscous incompressible fluid
195
over a wide range of Reynolds numbers (I< RI++ /103). The flow fields, pictures of the wake, the distribution of pressure and friction over the surface of the body flowed round, the drag coefficients, the coordinates of the points of separation of the flow and the extent of the stagnation zone were determined. Here we present some results, mainly on the methodical scheme, illustrating the accuracy and possibilities of the method. As an example we consider some results of the solution of the problem of the transverse flow past a circular cylinder of radius a, a = 1, whose axis is perpendicular to the velocity vector of the incident flow U,.
5
3y -z-
-
--Irs’
-U
\-L?u3
FIG. 2a
196
0. hf. Belotserkovskii, V..A. Gushchin and V. V. Shchennikov
-55
FIG. 2b Pictures of the flow for the Reynolds numbers 1, lo,30 and 50 (Re=2aU,ly) respectively are shown in Fig. 2,a. Figure 2b shows pictures of the flow for Re = lo3 at the instants T1 = 162.0, T2 = 166.0 and T3 = 170.0 respectively. In the last case (Re = 103) a non-stationary picture of the flow is observed (there is a definite growth of the stagnation zone and at a certain instant a “collapse” and ejection of fluid from the stagnation zone occurs). Possibly the result noticed here requires further study.
I 2
ff, zpaa /y
120
FIG. 3
FIG. 4
60
d
e, 2paa
Problems of the dynamics of a viscousincompressible fluid
197
PO+P, (f/2)$ 3
UC2
A
2
.
X
FIG. 5. 1 is the results of this paper, 2 is from [17], 3 is from [Ml, 4 is from [19]. Figure 3 shows the pressure distribution over the surface of the cylinder (Re = SO). Curve 1 corresponds to the results of this paper, curve 2 to the results obtained by the method of joint (block) relaxation [ 171. The distribution of the vortex function for the same Re number is given in Fig. 4. A detailed comparison with the results of other authors [ 18-201 and with some experimental data [21] is given. As illustrations we show the dependence on the Reynolds number of the pressure at the leading stagnation point (Fig. 5), the angle of separation of the flow measured from the trailing critical point (Fig. 6), and the total drag coefficient (Fig. 7). The agreement between the results is fairly good.
I
I
SO
100 Re
FIG. 6.1 is the results of this paper, 2 is from [17], 3 is from [18], 4 is from [19], 5 is the results of [20].
198
0. hf. Belotserkovskii. V. A. Gushchin and V. K Shchennikov
;9
CB
I.@\t
v, 5 \
B\ \
\
m% A\
0.5 -
8 \ #\ \
if’
$\
n* %
I
0
5
93TE:
7
2
l.g Re
FIG. 7. 1 is the results of this paper, 2 is the results of [21], 3 is from [18], 4 is from [ 191,s is the results of [20]. Figure 8 shows a general picture of the flow past a cube (2a is the linear dimension of the / v=l cube) and the velocity field u in the sections Q (x=- 00, Fig. 8,b) and 2=3u, Re=kU, (Fig. 8,~).
a
FIG. 8
Problems of the dynamics of a viscous incompressible fluid
199
5. In conclusion we point out that the proposed method possesses certain strong aspects of the SMAC method and of the modified MAC method [14] and in our opinion is free from a number of defects inherent in these methods. The difference scheme of this method permits the flow field to be calculated without the use of .the vortex and pressure values on the solid surface. The results of the calculations performed, and also trial calculations of three-dimensional flows carried out using this approach, confirm its efficiency. The difference scheme of the proposed method makes it possible to carry out by a single algorithm calculations of the flow of a viscous incompressible fluid past plane, axisymmetric and three-dimensional bodies of complex configuration, and also of internal flows over a wide range of Reynolds numbers. i?anshted by J. Berry REFERENCES 1.
BRAILOVSKAYA, I. Yu., KUSKOVA, T. V. and CHUDOV, L. A., Difference methods of solving the Navier-Stokes equations (review). In: Computing methods and programming (Vychisl. metody i programmirovanie). No. lI, 3-18, Izd-vo MGU, Moscow, 1968.
2. KUSKOVA, T. V. and CHUDOV, L. A., Approximate boundary conditions for the vortex in calculating the flows of a viscous incompressible fluid. In: Computing methods and programming (VychisL metody i progmmmirovanie). No. 11, 27-31, Izd-vo MGU, Moscow, 1968.
3.
LYUL’KA, V. A. and SHCHENNIKOV, V. V., Numerical solution of the Navier-Stokes equations. In: Theoretical papers on hydrodynamics (Sb. teor. rabot po gidrodinamike), 107-149, VTs Akad. Nauk SSSR, Moscow, 1970.
4.
LOLEZHAEV, V. I. and GRYAZNOV, V. L., A method of approximating the boundary conditions for the Navier-Stokes equations of an incompressible fluid. Proceedings of “The All-Union school for theoretical investigationof numerical methods of the mechanics of continuous medh” (Tr. “Wses. shkoly po teoretich. issL chislennykh metodov mekhan. sploshnoi sredy”), p. 20-26, December 1973, Zvenigorod. Theses of reports (Tezisy dokladov), 39-40, rotaprint IPMekhan. Akad. Nauk SSSR, Moscow 1973.
5.
VLADIMIROV, N. N., KUZNETSOV, B. G. and YANENKO, N. N., Numerical calculations of the symmetric plane flow of a viscous incompressible fluid past a lamina. In: Some questions of computational and applied mathematics (Nekotorye voprosy vychisl. i prikL matem.), 186-192, Izd-vo SO Akad. Nauk SSSR, Novosibirsk, 1966.
6.
CHORIN, A. .I., A numerical method for solving incompressible viscous flow problems. J. Comput. whys., 2, 1, 12-26, 1967.
7.
YANENKO, N. N., The weak approximation 6, 1431-1434, 1964.
8.
SAMARSKII, A. A., Introduction to the theory of difference schemes (Vvedenie v teoriyu raznostnykh skhem), “Nauka”, Moscow, 197 1.
9.
TAYLOR, T. D. and NDEFO, E., Calculation of the flow of a viscous fluid in a channel by the splitting method. In: Numerical methods in the mechanics of fluids (Chislennye metody v mekhan. zhidkostei), 218229, “Mu”, Moscow, 1973.
of systems of differential equations. Sibirskii matem zh., 5,
10. CHORIN, A. J., Numerical solution of Navier-Stokes
equations. Math. Comput., 22,745-762,
1968.
11. HARLOW, F. H., Numerical method of particles in cells for problems of hydrodynamics. In: Computational methods in hydrodynamics (Vychisl. metody v gidrodinamike), 316-342, “Mir”, Moscow, 1967. 12. HARLOW, F. H. and WELCH, J. E., Numerical calculation of time-dependent fluid with free surface. Phys. Fhaiis, 8,2182-2189, 1965.
viscous incompressible flow of
13. AMSDEN, A. A. and HARLOW, F. H., The SMAC method. Los Alamos Scient. Lab. Rept. La-4370, 1970. 14. EASTON, C. R., Homogeneous boundary conditions for pressure in MAC method. J. Comput. Phys., 9,2, 375-379, 1972. 15. MARCHUK, G. N., Methods of computational mathematics (Metody vychislitel’noi matematiki), “Nauka”, Novosibirsk, 1973.
Yu. A. Budyanskii,et al.
200
16. YANENKO, N. N. and SHOKIN, Yu. I. On the correctness of the first differential approximations schemes. Dokl.Akad. NaukSSSR, 182,4,116-178, 1968.
of difference
17. GUSHCHIN, V. A. and SHCHENNIKOV, V. V. A numerical method of solving the Navier-Stokes Zh. vychisL Mat. mat. Fiz., 14,2, 512-520, 1914.
equations.
18. TAKAMI, H. and KELLER, H. B., Steady twodimensional circular cylinder. Phys. Fluids, 12, 12, 51-56, 1969.
viscous flow of an incompressible fluid past a
19. DANNIS, S. C. and CHANG, G., Numerical solution for steady flow past a circular cylinder at Re numbers up to 100. J. FIuidMech., 48,3,471-489, 1970. 20. HAMILEC, A. E. and RAAL, J. D., Numerical studies of viscous flow around a circular cylinder. Phys. Fluids, 12, 1, 11-17, 1969. 21. TRITTON, D. J., Experiments on the flow past a circular cylinder at low Re.J. FWiMech.,
6,547-567,1959.
DETERMINATION OF A SEISMIC SECTION BY THE WAVE REFLECTION METHOD* Yu.
A. BUDYANSKII, A. V. KUDRYA, N. I. OGURTSOV and A. A. CHEKASIN Khar’kov-Novoaleksandrovsk (Received 3 May 1972) (Revised version 1 July 1974)
FOR THE solution of the problem indicated an algorithm is proposed which is applicable to a fairly wide class of layered structures, including the case of continuous variation of the speed of the seismic ray with depth. The reflected wave method is widely used in the solution of an extensive class of geophysical problems associated with the study of the structure of the upper part of the Earth’s crust, The basis of this method was laid in [l] where a theoretical justification was given, and in [2] uniqueness theorems were proved for some classes of sections for the solution of direct and converse problems. For the simplest sections (the plane problem), consisting of two or three horizontal layers, various approximate solution schemes have been proposed. The solution of the direct problem for a plane horizontally-layered homogeneous section reduces to the evaluation of the following integrals:
(1)
(2) where M is the number of strata, Vi(Z) is a function of the velocity in the i-th stratum, Zi, Zi+l are the coordinates of the roof and base of the i-th stratum respectively, &=ain c~, j/vi= sin ai+t, j/vi++ is a parameter, and fIM is the time of passage of the j-th seismic ray from the point of the explosion *Zh. vjkhisL Mat. mat. Fiz., 15, 1, 208-216, 1975.