- Email: [email protected]
Essentially non-uniform meshes for the numerical solution of the Navier-Stokes equations
USSR
Comput. tiths. Mzth. Phys Vol. 22, No. 6, pp. 183-193,
1982.
Printed in Great Britain
0041-5553/82$07.50+.00 01984. Pergamon Press Ltd.
ESSENTIALLY NON-UNIFORM MESHES FOR THE NUMERICAL SOLUTION OF THE NAVIER-STOKES EQUATIONS* V. I. KOPCHENOV, A. N. KRAIKO and M. P. LEVIN Moscow (Received 20 January 1981; revised 27 March 1981)
TAKING the example of the problem of plane viscous incompressible fluid flow into a square cavity with moving upper wall, the possibility of solving the problem mnnerially is demonstrated by using essentially non-uniform meshes, with a degree of non-uniformity that increases with the Reynolds’ number. 1. While the computation of laminar flows in the context of the complete Navier-Stokes equations does not usually involve fundamental difficulties at low and moderate Re numbers, it becomes difficult to obtain reliable results as Re increases (Re > 1000); hence the scope for numerical modelling is seriously restricted in the case of laminar flows for which the boundary layer approximation is unsuitable, and in particular, it becomes impossible to study the loss of stability of the flow. This is linked with the fact that, when it is impossible to use a division of the flow into a domain of non-viscous flow and a viscous boundary layer, thin zones of large parameter gradients appear in the flow at large Re numbers along with domains of relatively smooth parameter variation. These thin zones become thinner as Re increases. The main cause of the difficulties that arise at large Re numbers is linked with so-called scheme dissipative effects, which can, in such situations, partially or completely mask the actual effects of the viscosity or heat conduction etc. On the other hand, when solving the equations of the laminar boundary layer, no such problems arise whatever the Re number. This familiar fact is obvious, since the difference mesh here in practice only covers the boundary layer, and is automatically deformed as the layer thickness varies. This suggests a fairly obvious, but by no means easily realized method of overcoming the difficulties mentioned: namely, to use essentially non-uniform meshes, tracking the zones of sharp variation of the parameters. In spite of the fact that what has been said seems natural enough, and the fact that essentially nonuniform meshes are finding wider applications (see, for example, [l-3] ), further explanations seem necessary in view of the attention paid to scheme dissipative effects. Let us assume that the difference scheme approximates to k-th order the convective terms of the equations of motion and energy. Then the corresponding scheme dissipative terms, proportional to hk, where h is the characteristic dimension of the mesh, will sooner or later, as Re increases, with h fured, exceed the true dissipative terms appearing in the same equations with coefficient Re-1. Needless to say, this will not happen if h decreases in proportion to Re-llk, while the total nwnber of points N increases in proportion to Re2lk in the plane case, or Re3/k in the three-dimensional case. Yet it is obvious that, for large Re and reasonable k = 2, . . . ,6, where the upper value of k is taken with a margin, such an increase in N, and hence of the computing time, will not be realistic. However, the actual situation is not quite so hopeless.
*Zh. VFchisl
Mat. mat. Fiz., 22,6,1457-1467,
1982.
183
184
V.t Kopchenov, A. N. Knaiko and M. P. Levin
It is in fact only usually necessary to take accurate account of the non-ideal gas effects in the above-mentioned zones of large gradients (z.1.g.). Since these zones cannot be extended in ah directions, they will be located, in the two-dimensional case, in the neighbourhoods of certain lines 1 or of points, while in the three-dimensional case they are in the neighbourhoods of surfaces u or of manifolds of lower dimensionality. Let n be the coordinate, measured along the normal to 1 or to u, 6 the thickness of the zone along n, and Y, the projection of the velocity vector Vonto the normal. Assuming that the relevant convective and viscous terms are order-wise equal, e.g. u, (&I,/&z) and Re-’ (B2u,l&zz), where v, is the component of V tangential to 1 or u, it can be shown that u,&i-Re-‘.
(1.1)
On the other hand, the error of the approximation of the convective term has order
(1.2) where iT = n/S. Transforming in the viscous term from r to 5 and using (1. l), we find that
a*v, vn a2U, Re-'_--_ Re-’ a%, -Fm' 62 aft dn*
(1.3)
Since, within the zone considered, E has a variation of order unity, then dh+%~,‘dW’ , a2v,lan* are quantities of the same order. Consequently, the ratio of the scheme to the true dissipative terms, i.e. of the right-hand sides of (1.2) and (1.3), is proportional to (h/6)k. Hence it is clear that, if a sufficiently coarse mesh in the n direction is taken in the z.1.g. (hn = 6/i, where I is a fairly large number, say lo), then the ratio mentioned will be small. Notice, incidentally, that in the absence of in-flow into the wall z.1.g. (though not necessarily into the boundary layers in the generally accepted sense), by virtue of the equation of continuity vn - S and by (1.1) we have the boundary layer estimate 6 - Re-s. Hence the h, here, as in the case of an ordinary boundary layer, has to be reduced in proportion to Re-s as Re increases. Naturally, the coefficients of proportionality are in the general case different for different walls and problems, and need to be determined during the computation (e.g. from preliminary computations performed for two values of Re). Of course the estimates made are extremely crude and can often be refined. All in all, it is desirable to compute flows in the context of the complete system of NavierStokes equations as follows. The difference mesh must be much more sparse with respect to n in the z.1.g. than on other domains; as Re increases, h, needs to decrease in proportion to S, The dependence 6 = S (Re) is established and refined during the computation, possibly drawing on estimates of the type given above, or on the results of much more strict asymptotic analysis, valid in the context of laminar flow conditions as Re + 00.Next, in the situations discussed, the flow outside the z.1.g. in the zones of moderate gradients (z.m.g.). - depends on the vorticity distributions on the parts of the zone boundaries where gas flows into the z.m.g. from z.l.g., while in the z.m.g. themselves the dissipative terms are insignificant. Hence the computation in these zones can be performed on much coarser meshes, independent of the Re value. In order to avoid an unjustified increase in the total number of computational points, the passage from the fine mesh in z.1.g. to the coarse mesh in z.m.g. is best made rapidly. This necessitates using numerical schemes providing approximation on essentially non-uniform meshes. This point is extremely important, since by no means all difference schemes have this property; as the order of approximation increases, stricter demands are made on the mesh uniformity. In such cases the situation is not saved by passage to new independent variables, since similar demands are then made on the coefficients of the equations of the transformations. In the light of these remarks, the spline scheme of [4,5] looks very good; it ensures quite a high order of approximation on essentially non-uniform meshes. This is why the scheme is used in the present
Numericnl solution of the Navier-Stokes equations
185
paper. In [6] another method is described for compensating the approximation errors due to mesh non-uniformities in problems of the present type. Before describing the problem, the numerical scheme, and the results obtained, we recall that the cavity problem is often taken as a test problem when solving the Navier-Stokes equations numerically. The reason is that, in this problem, even at Re = 400, the flow includes quite pronounced z.1.g. close to the walls and the core, which is a z.m.g. As Re increases, the qualitative difference between these zones becomes more marked. This fact, along with the impossibility of using here the classical division into a boundary layer and non-viscous core, in which the flow is to a Arst approximation independent of the boundary layer, explains the limited success achieved in most solutions of the cavity problem on moving into the domain of large Re (Re > 103); this is confirmed particularly by the surveys of [7-9 ] , which give a good picture of the relevant research. Progress into the domain of Re greater than lo3 is also rendered difficult by the fact that computational instability arises here for many schemes. It has been suggested, see [8, lo], that one way of overcoming the instability is to approximate the convective terms of the equations to first order; this leads to effects of scheme viscosity that are even more serious, especially for uniform meshes. Notice in this connection that the analysis in [ 1 l] reveals pronounced distortions of the flow picture, which introduce the scheme viscosity for the typical meshes, such as were used notably in [8]. For similar reasons, most publications containing computations for Re numbers up to 104 (see, for example, [ 12]), have merely revealed the stability of the algorithms (the absence of emergency stops), without giving reasonably accurate numerical data. 2. Consider the plane stationary flow of a viscous incompressible fluid in a square cavity with an upper wall that moves from left to right with constant velocity U. It will be convenient to use dimensionless quantities throughout. As the characteristic parameters we choose the cavity height d, V, and the coefficient of kinematic viscosity u, which is henceforth assumed constant. Then the Navier-Stokes equations can be written as (2.1) Here, x and y are Cartesian coordinates with origin at the bottom left corner of the cavity (x and y are referred to d), J, is the stream function, u=aglay , v=-9$lds are the x andy components of the velocity vector, -%=&~/ax-du/ay - is the vorticity, and Re = Ud/v. On the walls we formulate the adhesion conditions (the subscript w is assigned to quantities on the wall, n is measured along the normal to the wall)
(2.2) The stationary solution of our problem will be found in the process of establishment with respect to time f, using the implicit scheme of the method of alternating directions. For this, we add the term a$/&, having no physical significance, to the first of Eqs. (2.1), and the derivative at/at to the second. Along with 5 and $, we shall also define their first and second derivatives with respect to x and y. As mentioned above, we use the scheme of [4], written in the following form, for the computations. Let 27 be the t step, let I and J be the numbers of the nodes of the difference mesh in the x and y directions, including the nodes on the domain boundary, and let fbe r or J/. The derivatives at node (i, j) will be approximated by one-dimensional cubic M,jf==02fldx2 Lijf=ayay2 will be found by using the same interpolation splmes. The derivatives nt ijf=af/dz , Lj’-aj/ay splines. The scheme has the form
186
K 1 Kopchenov, A N. Kraiko and M R Levi?!
2
i=l,
I
,‘**,
j=2,3
,
)
.
. . ) J-l,
~i;+‘=g,y+” +z [ _lij” ( rnijE)“+i+nij* (lijf)n+‘!’+ &(L”i)
+
&(Mi,i)"+l]
1
i=2,3
I-l,
,...,
2,. . . ,J,
j=l,
I$:,+‘”=$ij”+T[ (Lijv)s+‘i2+(Mij’)‘-f;:+’
1,
n+‘h
i=2,3,.
(2.3)
. , , Z-1,
j=2,3 ,*..,J--1,
$Y,”=lj?,~+“‘+T[
‘+“‘f
(Lij”)
j=2,3 ,...I
‘+i-l;il+i],
(Mij*)
i=2,3,
, . . , I-1,
J-l.
Moreover, for cubic splines we have the relations [ 131 hjM:_,,j+2(hi+hi+*)Mijf+hi+lMi:,,, =6
ft+i,j-fij
(
I11t1
_
fij-f+i,j hi
m,,’ = hiMvj’ + +M;_,,j 3
mi,
f=+&f-_
where hi=Xi-X, -i.
hi+, 6
M:+,
i=2,3 , . . . , z-1,
1’ +
fij-fi-l,j
j +
hi
i=2,3,.
’
fi+itj-fii
’
hi+,
’
(2.4) .
9
I
7
i=l, 2, .*,I-&
Similar relations hold for Liif and lijf.
One time layer is computed as follows. The field 5 on it is first calculated; the derivatives * are taken from a known layer. Then, for the { distribution obtained, we solve the lij” 3 mij equation defining $. An algorithm is used in which the difference equations (2.3), supplemented by relations (2.4), are solved simultaneously at each half-step. The system of algebraic equations for 9 is closed by means of boundary conditions (2.2). The second and third of relations (2.4) enable the first derivative of J, with respect to the coordinate to be expressed in terms of the value of J/ and of the second derivatives of $ on the wall and at the nearest interior point of the flow field. With regard to the difference relations for {, in order to complete the algebraic system we have to complete the definition off and of its second derivatives on the wall, whereas the boundary conditions for f are not present: in the initial statement (2.1) and (2.2). To obtain the missing relations, we shall use corollaries of the initial equations and boundary conditions. For the first Eq. (2.1), noting that J/,.+,= 0, we find that LTx= ( iY$/dn2)
“I.
(2.5)
From the second Eq. (2.1) and the adhesion conditions, we obtain the relation
(2.6)
Numeriml solution of the Navier-Stokes
equations
187
Here, u is measured along the wall, X = 1 on the moving wall, and x = 0 on the fmed wall. The quantity {,n+l is found from (2.5) by using the second derivative of $ on the wall on the known time layer. On a half-integer layer n + 45,the quantity SWis found up to 0 (72) from { on the wall at the (n + I)-th, n-th, and (n - l)th, time layers. In accordance with (2.6), k=l, ,
(,qf,;)
n+‘=_
(&K)
n+‘!z,
k=l,
I,
k=J, j=l,
i=l, 2 ,...,I, 2,. . ., J.
It may be recalled that the second derivatives of the vorticity with respect to the coordinate were found on the wall in [4] by using, instead of (2.6), a similar corollary of the non-stationary equation. At the upper comers of the cavity, u breaks off, and hence I{ I becomes i&rite. In this connection the question arises of the influence of the method of completing the definition off at these points or in their neighbourhoods in the flow as a whole. Two ways of specifying { were tried. In the first U, in the celI adjacent to the comer point was smeared from 0 to 1; in the second, 5 at this point was replaced by the half-sum of 5 at adjacent points of the upper and side walls. In both cases, the difference between the calculated fields u and v was virtually everywhere a fraction of a percent. This is easily appreciated, if we note that, in accordance with the solution obtained in [ 141, which holds for small neighbourhoods of these points, their influence manifests itself at relative distances of the order of Re- 1 when the thickness of the wall layers is of the order of Re-‘A. Mention was earlier made in [9] of the small influence of the choice of the conditions at the upper corners.
FIG. l.-for N=51X51,A-for N=41X41, o-for N = 29X 29, 0 - for N = 40X40, from 1201
The choice of the time integration step is essentially dependent on the statement of the boundary conditions for {, especially for large Re numbers; this was revealed by the analysis made in [ 15-171. The statement of the boundary conditions employed in the present paper ensured stable computation over a wide range of Re numbers with a time integration step exceeding the period of the mesh.
188
V. L Kopchenov, A. N. Kraiko and M P. Levin
FIG. 2 Let us consider in more detail some aspects of the solution of the systems of difference equations. The equation for J/ is reducible to a system of scalar difference equations with a tridiagonal matrix of the coefficients with respect to the unknown nodal values of the stream function; these equations were solved by scalar pivotal condensation. The equation for { reduces to a system of vector three-point difference equations with respect to the unknown nodal values {ii and the derivatives of Mijor Lij. These were solved, either by Seidel’s method, or by vector pivotal condensation [ 18,191. On uniform meshes, cubic splines provided fourth-order approximation for the first derivative of the function with respect to the coordinate, while the second derivatives are approximated to second-order regardless of the mesh uniformity [4]. Though the order of approximation of the first derivatives drops to the third on non-uniform meshes, this is still high enough for the convective terms of the initial system of equations. 3. Using the method described, we wrote a program in FORTRAN for the BESM-6 computer. Some results, obtained on uniform meshes for Re = 100 and 400, are illustrated in Figs. 1 and 2. For these Re numbers, the computations were made with different numbers of computational nodes N, and were compared with the data in [20]. In Fig. 1 we show the profiles u = U(Y) on the centre line x = 0.5 for Re = 400, obtained for different N, and the similar profile in [20]. There is clearly a marked difference between our curves and that in [20] , though for Re = 100 the difference between the profile in [20] and that computed for N = 29X 29 is within the error of graphical representation. The difference with Re = 400 is primarily due to the inaccuracy of the data in [20]. This is revealed by Fig. 2, a and b, referring to Re = 100 and 400, where I$ Im is the maximum modulus of $, and h is the mesh step (1 refers to our data, and 2 to the data in [20]). The right-hand branches of the curves in Fig. 2, a are plotted, using points outside the graph. The method based on spline approximation clearly has faster convergence with respect to h than the method in [20]. This is particularly noticeable for Re = 400. For the case Re = 100 on a 21 X 21 uniform mesh we compared the efficiencies of the two methods mentioned above for solving the equations for t, and we studied the domain of stability with respect to T. It was found that, for Seidel’s method, the maximum tolerance half-step rm * 0.9h is the same as the optimal half-step (in the sense of computing time till establishment). For vector pivotal condensation, the maximum tolerable rm = 1.9/z, while the optimal is coopt= 1.3h. By comparison with the algorithm based on Seidel’s method, the computing time with vector
Numerical solution of the Navier-Stokes equations
189
pivotal condensation was reduced by a factor of 5 for r = 0.9h, and by a factor of 7 for r = I.3h In the latter case the computation took 5 minutes till establishment. These results show that, for the difference schemes employed, even with Re = 400 the uniform meshes enable reasonably high accuracy to be attained only for N > 40X 40. To progress to greater Re without seriously increasing N, non-uniform meshes must be used, in accordance with what has been said above. Let us describe how we constructed the mesh in the present paper, on the basis of special methodic computations and the considerations given above. In the longitudinal and transverse directions, the mesh was constructed according to the same principle, and was characterized by its minimum (hmo and h,) and maximum hm steps; the values of hmo were different for the lower and the other (hm) walls (hm < h, 0). Like the thicknesses 6 and S 0 of the corresponding wall layers, they were found from computations with Re = 400 and were reduced in proportion to Re-s as Re increased. The hm step in the core was taken to be virtually the same as for Re = 400. The transition from hm or hm 0 to hm was made outside the wall layers according to a geometric progression with “ratio” 4 = 1.4 - 2. Since it is made very rapidly, the total number N of cells scarcely increases as Re increases. Computations showed that, to obtain stable computation within the upper and side wall layers, several (in our case four) mesh layers have to be introduced, with cell size (along the normal to the wall) h, + < h, . The value of hm + was found from the expression h,+ = lO/Re; analysis of the simplified equations showed that this device is apparently connected with the good conditionality of the difference boundary value problem. The transition from h, + to hm, is made according to a geometric progression with the same ratio. The introduction of extra layers only involves a slight increase of N as Re increases.
-0.5
0
a.5
FIG. 3 The possibilities of non-uniform meshes are demonstrated by the following example. The flow at Re = 400 was computed on a non-uniform 27X 27 mesh with h,+=h,=0.02, hm0=0.035 and hm = 0.065 with (I = 1.4. It has been seen from Fig. 1 that the uniform mesh with roughly the same number of nodes gives extremely crude results. On the other hand, the velocity profile u = u (Y) on the line x = 0.5, obtained for a non-uniform division, is virtually the same as the profile computed on the much finer uniform mesh with N = 5 1 X 5 1. The practical value of the above principles for constructing the non-uniform mesh was confirmed, not only by computations for Re = 400, but also by extra data obtained for Re = 103. In the latter case, two meshes were used. Mesh 1 was constructed in the way explained above. Mesh 2 had twice as many nodes in the wall layers. Mesh 1 contained 39X 34 nodes with h,+=0.01, h,=O.O12, h,0=0.022 and hm = 0.053. Mesh 2 had more nodes, N = 49X 5 1, though the maximum step in
190
V. I. Kopchenov, A. N. Kraiko and M I! Levin
FIG. 4
L.
-n.i5
n
0.5
FIG. 5 the core remained roughly the same. Doubling of the number of nodes in the wall layer did not produce a serious correction of the results, notably for the velocity I( profile at the cavity centre. The accuracy of the results can be judged from the following data. The values of the maximum speed of reverse current 1r.P I, obtained on meshes 1 and 2, are respectively equal to 0.392 and 0.382. In Fig, 3 the present results (continuous curve) are compared with the data in [21] on the u profile in the cross-section passing through the vortex centre (by the vortex centre we mean the point where I $ I takes its maximum value). The broken curve of Fig. 3 shows the results obtained in [21] using a non-uniform mesh with N = 17 X 17 and a spline scheme of higher order than used in the present paper. The dot-dash curve was computed in [21] by a finite-difference scheme with uniform division and the same N = 17 X 17. The analogous results in [21] for a uniform mesh with N= 65 X 65 agree (within the scale of the graph) in different sections either with the continuous or the broken curve. In Fig. 4 we show the streamline picture with Re = 1000 (the figures on some lines refer to the values of $).
Numerical solution of the Navier-Stokes equations
191
FIG. 6 The data quoted show the advantage of using a non-uniform mesh. As Re increases the advantage becomes more obvious. In principle, it seems that this device could be used to move up to the limit of existence of stationary modes. While the establishment time naturally increases with Re (as the viscosity falls, the “acceleration” of the flow requires more “physical” time), the only serious obstacle to advancing towards higher “laminar” Re would seem to be the non-stationary property of the flow. While admittedly the present algorithm does not describe the actual non-stationary flow, due to the “stationary” boundary conditions for c, and the artificial term a$/& added to the first of Eqs. (2.1), nevertheless, establishment took place in the algorithm at least up to Re = 3.5X 103. Some results for this case are shown in Figs. 5,6, where we give the profile u CJ) for x = 0.5 and the and km = streamline picture for a 47 X 45 mesh with h,+ =2.5. 10m3, h,=0.007, h,“=0.01 0.052. Let us conclude with some remarks on the construction of non-uniform meshes. In the present paper the mesh was fmed and was specified at the start of the computation, on the basis of an analysis of the results for relatively low Reynolds numbers. Since this approach can often be very laborious, it is very important to examine the automatic construction of computational meshes, ensuring satisfactory resolution of the zones of large gradients. Of work in this field, we mention [l] and [22], in which it is proposed to use coordinate transformations dependent on the velocity components. In [23], when solving the stationary problem by the establishment method, the mesh is constructed on each time step by minimizing a certain functional. The method described in [24]’ has become quite popular; here, the node coordinates are found by solving a system of two quasi-linear inhomogeneous equations of elliptic type [25]. In this latter case, however, the disposition of the nodes within the computational domain can only be controlled by compressing them on the boundaries, which usually demands preliminary information about the disposition of the zones of large gradients. In the authors’ view, genuine scope is offered by the use of a dialogue type of
192
V: L Kopchenov, A. N. Kraiko and M P. Levin
interaction with the computer, in which the mesh is corrected during the computing process according to the parameter distributions that are realized. Translated by D. E. Brown
REFERENCES 1. TOLSTYKH, A. I., A method for the numerical solution of the Navier-Stokes equations of a compressible gas over wide range of Reynolds numbers, DokL Akud. Nauk SSSR, 210, No. 1,48-51, 1973. 2. DAIKOVSKII, A. G., POLEZHAEV, V. I., and FEDOSEEV, A. I., Study of the structure of transient and turbulent modes of convection in a vertical layer, Izv. Akad. Nauk SSSR, Mekhan zhidkosti gaza, No. 6, 66-75,1978. 3. DAIKOVSKII, A. G., POLEZHAEV, V. I., and FEDOSEEV, A. I., Numerical modelling of transient and turbulent convection modes using the non-stationary Navier-Stokes equations, Preprint IPMekhan. Akud Nauk SSSR, No. 101, Moscow, 1978. 4. RUBIN, S. G., and GRAVES, R. A., Viscous flow solutions with a cubic spline approximation, Fluids, 3, No. 1, l-36, 1975.
Cbmput.
5. RUBIN, S. G., and KHOSLA, P. K., Higher-order numerical solutions using cubic splines, AIAA J., 14, No. 7, 851-858,1976. 6. LEWIS, E., Steady flow between a rotating circular cylinder and a fixed square cylinder, J, Fluid Mech., 95, No. 3,497-514,1979. 7. KUSKOVA, T. V., Numerical study of two-dimensional flows of a viscous incompressible fluid, in: Applications ofrhe mesh method in gas dynamics (Nekotorye primeneniya metoda setok v gazovoi dinamike), No. 3, 7-136, VTs MGU, Moscow, 1971. 8. BOZEMAN, J. D., and DALTON, C., Numerical study of viscous flow in a cavity, J. Comput. Phys., 12, No. 2, 348-363,1973. 9. GHIA, K. N., HANKY, V. L. Jr., and HODGE, J. K., Solution of Navier-Stokes equations allowing for incompressibility in the usual variables, in: Rocket techniques and cosmonautics (Russian translation, Mir, Moscow, 1979). 10. HOSMEN, A. D., et aL, Numerical methods of studying viscousjluidflows 1972).
(Russiantranslation,hfir,Moscow,
11. DE VAHL, D. G., and MALLINSON, G. D., An evaluation of upwind and central difference approximations by a study of recirculating flow, Comput. Fluids, 4, No. 1,29-43,1976. 12. KURTZ, L. A., et aL, A comparison of the method of lines to finite difference techniques in solving timedependent partial differential equations, Comput. Fluids, 6, No. 1,49-70, 1978. 13. AHLBERG, J. H., et aL, Theory of Splines and Their Applications, Academic Press, New York, 1967. 14. BATCHELOR, G. K., Introduction to Fluid Dynamics, Cambridge U.P., 1967. 15. DAIKOVSKII, A. G., POLEZHAEV, V. I., and FEDOSEEV, A. I, On calculation of the boundary conditions for non-stationary Navier-Stokes equations in “vorticity, stream function” variables, in: Numerical methods of the mechanics of a continuous medium (Chisl. metody mekhan. sploshnoi sredy), Vol. 9, No. 7,49-58, Nauka, Novosibirsk, 1978. 16. TARUNIN, E. L., On the choice of approximation formula for the velocity eddy at a rigid boundary when solving problems of viscous fluid dynamics, in: Numerical methods of mechanics of continuous medium (ChisL metody mekhan. sploshnoi sredy), VoL 9, No. 7,97-111, Nauka, Novosibirsk, 1978. 17. PEPPER and HARRIS, Numerical modelling of free convection in closed cavities by completely implicit method, in: Theoretical foundations of engineering calcuhtions, (Russian translation, Mir, Moscow, 1977). 18. VASIN, V. G., and POLEZHAEV, V. I., Difference schemes for Navier-Stokes equations of compressible gas and calculation of thermal-acoustic waves, Preprint IPMatem. Akad. Nauk SSSR, No. 124, Moscow, 1977.
Instability of a convergent spheriml shock wave
193
19. SAMARSKII, A. A., and NIKOLAEV, E. S., Methods for solving mesh equations (Metody resheniya setochnykh uravnenii), Nauka, Moscow, 1979. 20. BURGCRAF, 0. R., Analytic and numerical study of the structure of stalling flows, in: Collection of surveys and translations of foreign periodical literature, Mechanics, No. 6, (IOO), 51-90, Mir, Moscow,1966. 21. RUBIN, S. G., and KHOSLA, P. K., Polynomial interpolation J. Comput. Phys, 24,No. 3,217-244, 1977.
methods for viscous flow calculations:
22. TOLSTYKH, A. I., On the condensation of nodes of difference methods in the process of solution and on application of higher-accuracy schemes for the numerical study of viscous gas flows, Zh. vychisl Mar. mat. Fiz., 18, No. 1, 139-153, 1978. 23. DANAEV, N. T., LISEIKIN, V. D., and YANENKO, N. N., On the numerical solution of viscous thermally conducting gas flow past a solid of revolution on a moving curvilinear mesh, in: Numerical methods of the mechanics of a continuous medium (ChisL metody mekhan. sploshnoi sredy), Vol. 11, No. 2, 51-61, Nauka, Novosibirsk, 1980. 24. THOMPSON, J. F., THAMES, F. C., and MASTIN, W., Automatic numerical generation of a body-fitted curvilinear system for a field containing any number of arbitrary 2-D bodies, J. Comput. Phys., 15, No. 3,299-319,1974. 25. THOMAS, P. D., and MIDDLECOFF, J. F., Direct control of the grid point distribution in meshes generated by elliptic equations, AIAA J. 18, No. 6,652-656, 1980.
U.S.S.R Cornput. Maths Math. Phys. Vol. 22, No. 6, pp. 193-205.1982. Printed in Great Britain
0041-5553/82$07.50+.00 01984. Pergamon Press Ltd.
INSTABILITY OF A CONVERGENT SPHERICAL SHOCK WAVE* K. V. BRUSHLINSKII Moscow (Received
8 December
1980)
A CONVERGENT spherical shock wave is unstable in the linear approximation. Arbitrary small disturbances of the flow can be expanded in series in generalized spherical functions. The eigenvalues of the differential operators, characterizing the growth rate of the disturbances in the case of the first few spherical harmonics, are found numerically.
Introduction The convergence of a spherical shock wave to its centre of symmetry has been studied by Guderley [ 11, Landau and Stanyukovich [2], and Gel’fand et a~!(see survey [3]), and is a familiar example of a similarity problem of gas dynamics. In this problem the power k in the similarity variable { = tr-k is transcendental and can only be determined during numerical solution of the boundary value problem for an ordinary differential equation. This creates special problems when studying the shock wave stability. Particular mention must be made of [4], where, under simplifying assumptions (channel approximation), an analytical solution is obtained for the convergent shockwave problem, and the instability of the wave is proved. In more complete studies, however, there is necessarily an element of numerical analysis, if only because the similarity solution itself is found numerically. *Zh. vychisl. Mat. mat. Fiz., 22,6, 1468-1479, “bSK.‘.‘:h-:,
1982.
Comput. tiths. Mzth. Phys Vol. 22, No. 6, pp. 183-193,
1982.
Printed in Great Britain
0041-5553/82$07.50+.00 01984. Pergamon Press Ltd.
ESSENTIALLY NON-UNIFORM MESHES FOR THE NUMERICAL SOLUTION OF THE NAVIER-STOKES EQUATIONS* V. I. KOPCHENOV, A. N. KRAIKO and M. P. LEVIN Moscow (Received 20 January 1981; revised 27 March 1981)
TAKING the example of the problem of plane viscous incompressible fluid flow into a square cavity with moving upper wall, the possibility of solving the problem mnnerially is demonstrated by using essentially non-uniform meshes, with a degree of non-uniformity that increases with the Reynolds’ number. 1. While the computation of laminar flows in the context of the complete Navier-Stokes equations does not usually involve fundamental difficulties at low and moderate Re numbers, it becomes difficult to obtain reliable results as Re increases (Re > 1000); hence the scope for numerical modelling is seriously restricted in the case of laminar flows for which the boundary layer approximation is unsuitable, and in particular, it becomes impossible to study the loss of stability of the flow. This is linked with the fact that, when it is impossible to use a division of the flow into a domain of non-viscous flow and a viscous boundary layer, thin zones of large parameter gradients appear in the flow at large Re numbers along with domains of relatively smooth parameter variation. These thin zones become thinner as Re increases. The main cause of the difficulties that arise at large Re numbers is linked with so-called scheme dissipative effects, which can, in such situations, partially or completely mask the actual effects of the viscosity or heat conduction etc. On the other hand, when solving the equations of the laminar boundary layer, no such problems arise whatever the Re number. This familiar fact is obvious, since the difference mesh here in practice only covers the boundary layer, and is automatically deformed as the layer thickness varies. This suggests a fairly obvious, but by no means easily realized method of overcoming the difficulties mentioned: namely, to use essentially non-uniform meshes, tracking the zones of sharp variation of the parameters. In spite of the fact that what has been said seems natural enough, and the fact that essentially nonuniform meshes are finding wider applications (see, for example, [l-3] ), further explanations seem necessary in view of the attention paid to scheme dissipative effects. Let us assume that the difference scheme approximates to k-th order the convective terms of the equations of motion and energy. Then the corresponding scheme dissipative terms, proportional to hk, where h is the characteristic dimension of the mesh, will sooner or later, as Re increases, with h fured, exceed the true dissipative terms appearing in the same equations with coefficient Re-1. Needless to say, this will not happen if h decreases in proportion to Re-llk, while the total nwnber of points N increases in proportion to Re2lk in the plane case, or Re3/k in the three-dimensional case. Yet it is obvious that, for large Re and reasonable k = 2, . . . ,6, where the upper value of k is taken with a margin, such an increase in N, and hence of the computing time, will not be realistic. However, the actual situation is not quite so hopeless.
*Zh. VFchisl
Mat. mat. Fiz., 22,6,1457-1467,
1982.
183
184
V.t Kopchenov, A. N. Knaiko and M. P. Levin
It is in fact only usually necessary to take accurate account of the non-ideal gas effects in the above-mentioned zones of large gradients (z.1.g.). Since these zones cannot be extended in ah directions, they will be located, in the two-dimensional case, in the neighbourhoods of certain lines 1 or of points, while in the three-dimensional case they are in the neighbourhoods of surfaces u or of manifolds of lower dimensionality. Let n be the coordinate, measured along the normal to 1 or to u, 6 the thickness of the zone along n, and Y, the projection of the velocity vector Vonto the normal. Assuming that the relevant convective and viscous terms are order-wise equal, e.g. u, (&I,/&z) and Re-’ (B2u,l&zz), where v, is the component of V tangential to 1 or u, it can be shown that u,&i-Re-‘.
(1.1)
On the other hand, the error of the approximation of the convective term has order
(1.2) where iT = n/S. Transforming in the viscous term from r to 5 and using (1. l), we find that
a*v, vn a2U, Re-'_--_ Re-’ a%, -Fm' 62 aft dn*
(1.3)
Since, within the zone considered, E has a variation of order unity, then dh+%~,‘dW’ , a2v,lan* are quantities of the same order. Consequently, the ratio of the scheme to the true dissipative terms, i.e. of the right-hand sides of (1.2) and (1.3), is proportional to (h/6)k. Hence it is clear that, if a sufficiently coarse mesh in the n direction is taken in the z.1.g. (hn = 6/i, where I is a fairly large number, say lo), then the ratio mentioned will be small. Notice, incidentally, that in the absence of in-flow into the wall z.1.g. (though not necessarily into the boundary layers in the generally accepted sense), by virtue of the equation of continuity vn - S and by (1.1) we have the boundary layer estimate 6 - Re-s. Hence the h, here, as in the case of an ordinary boundary layer, has to be reduced in proportion to Re-s as Re increases. Naturally, the coefficients of proportionality are in the general case different for different walls and problems, and need to be determined during the computation (e.g. from preliminary computations performed for two values of Re). Of course the estimates made are extremely crude and can often be refined. All in all, it is desirable to compute flows in the context of the complete system of NavierStokes equations as follows. The difference mesh must be much more sparse with respect to n in the z.1.g. than on other domains; as Re increases, h, needs to decrease in proportion to S, The dependence 6 = S (Re) is established and refined during the computation, possibly drawing on estimates of the type given above, or on the results of much more strict asymptotic analysis, valid in the context of laminar flow conditions as Re + 00.Next, in the situations discussed, the flow outside the z.1.g. in the zones of moderate gradients (z.m.g.). - depends on the vorticity distributions on the parts of the zone boundaries where gas flows into the z.m.g. from z.l.g., while in the z.m.g. themselves the dissipative terms are insignificant. Hence the computation in these zones can be performed on much coarser meshes, independent of the Re value. In order to avoid an unjustified increase in the total number of computational points, the passage from the fine mesh in z.1.g. to the coarse mesh in z.m.g. is best made rapidly. This necessitates using numerical schemes providing approximation on essentially non-uniform meshes. This point is extremely important, since by no means all difference schemes have this property; as the order of approximation increases, stricter demands are made on the mesh uniformity. In such cases the situation is not saved by passage to new independent variables, since similar demands are then made on the coefficients of the equations of the transformations. In the light of these remarks, the spline scheme of [4,5] looks very good; it ensures quite a high order of approximation on essentially non-uniform meshes. This is why the scheme is used in the present
Numericnl solution of the Navier-Stokes equations
185
paper. In [6] another method is described for compensating the approximation errors due to mesh non-uniformities in problems of the present type. Before describing the problem, the numerical scheme, and the results obtained, we recall that the cavity problem is often taken as a test problem when solving the Navier-Stokes equations numerically. The reason is that, in this problem, even at Re = 400, the flow includes quite pronounced z.1.g. close to the walls and the core, which is a z.m.g. As Re increases, the qualitative difference between these zones becomes more marked. This fact, along with the impossibility of using here the classical division into a boundary layer and non-viscous core, in which the flow is to a Arst approximation independent of the boundary layer, explains the limited success achieved in most solutions of the cavity problem on moving into the domain of large Re (Re > 103); this is confirmed particularly by the surveys of [7-9 ] , which give a good picture of the relevant research. Progress into the domain of Re greater than lo3 is also rendered difficult by the fact that computational instability arises here for many schemes. It has been suggested, see [8, lo], that one way of overcoming the instability is to approximate the convective terms of the equations to first order; this leads to effects of scheme viscosity that are even more serious, especially for uniform meshes. Notice in this connection that the analysis in [ 1 l] reveals pronounced distortions of the flow picture, which introduce the scheme viscosity for the typical meshes, such as were used notably in [8]. For similar reasons, most publications containing computations for Re numbers up to 104 (see, for example, [ 12]), have merely revealed the stability of the algorithms (the absence of emergency stops), without giving reasonably accurate numerical data. 2. Consider the plane stationary flow of a viscous incompressible fluid in a square cavity with an upper wall that moves from left to right with constant velocity U. It will be convenient to use dimensionless quantities throughout. As the characteristic parameters we choose the cavity height d, V, and the coefficient of kinematic viscosity u, which is henceforth assumed constant. Then the Navier-Stokes equations can be written as (2.1) Here, x and y are Cartesian coordinates with origin at the bottom left corner of the cavity (x and y are referred to d), J, is the stream function, u=aglay , v=-9$lds are the x andy components of the velocity vector, -%=&~/ax-du/ay - is the vorticity, and Re = Ud/v. On the walls we formulate the adhesion conditions (the subscript w is assigned to quantities on the wall, n is measured along the normal to the wall)
(2.2) The stationary solution of our problem will be found in the process of establishment with respect to time f, using the implicit scheme of the method of alternating directions. For this, we add the term a$/&, having no physical significance, to the first of Eqs. (2.1), and the derivative at/at to the second. Along with 5 and $, we shall also define their first and second derivatives with respect to x and y. As mentioned above, we use the scheme of [4], written in the following form, for the computations. Let 27 be the t step, let I and J be the numbers of the nodes of the difference mesh in the x and y directions, including the nodes on the domain boundary, and let fbe r or J/. The derivatives at node (i, j) will be approximated by one-dimensional cubic M,jf==02fldx2 Lijf=ayay2 will be found by using the same interpolation splmes. The derivatives nt ijf=af/dz , Lj’-aj/ay splines. The scheme has the form
186
K 1 Kopchenov, A N. Kraiko and M R Levi?!
2
i=l,
I
,‘**,
j=2,3
,
)
.
. . ) J-l,
~i;+‘=g,y+” +z [ _lij” ( rnijE)“+i+nij* (lijf)n+‘!’+ &(L”i)
+
&(Mi,i)"+l]
1
i=2,3
I-l,
,...,
2,. . . ,J,
j=l,
I$:,+‘”=$ij”+T[ (Lijv)s+‘i2+(Mij’)‘-f;:+’
1,
n+‘h
i=2,3,.
(2.3)
. , , Z-1,
j=2,3 ,*..,J--1,
$Y,”=lj?,~+“‘+T[
‘+“‘f
(Lij”)
j=2,3 ,...I
‘+i-l;il+i],
(Mij*)
i=2,3,
, . . , I-1,
J-l.
Moreover, for cubic splines we have the relations [ 131 hjM:_,,j+2(hi+hi+*)Mijf+hi+lMi:,,, =6
ft+i,j-fij
(
I11t1
_
fij-f+i,j hi
m,,’ = hiMvj’ + +M;_,,j 3
mi,
f=+&f-_
where hi=Xi-X, -i.
hi+, 6
M:+,
i=2,3 , . . . , z-1,
1’ +
fij-fi-l,j
j +
hi
i=2,3,.
’
fi+itj-fii
’
hi+,
’
(2.4) .
9
I
7
i=l, 2, .*,I-&
Similar relations hold for Liif and lijf.
One time layer is computed as follows. The field 5 on it is first calculated; the derivatives * are taken from a known layer. Then, for the { distribution obtained, we solve the lij” 3 mij equation defining $. An algorithm is used in which the difference equations (2.3), supplemented by relations (2.4), are solved simultaneously at each half-step. The system of algebraic equations for 9 is closed by means of boundary conditions (2.2). The second and third of relations (2.4) enable the first derivative of J, with respect to the coordinate to be expressed in terms of the value of J/ and of the second derivatives of $ on the wall and at the nearest interior point of the flow field. With regard to the difference relations for {, in order to complete the algebraic system we have to complete the definition off and of its second derivatives on the wall, whereas the boundary conditions for f are not present: in the initial statement (2.1) and (2.2). To obtain the missing relations, we shall use corollaries of the initial equations and boundary conditions. For the first Eq. (2.1), noting that J/,.+,= 0, we find that LTx= ( iY$/dn2)
“I.
(2.5)
From the second Eq. (2.1) and the adhesion conditions, we obtain the relation
(2.6)
Numeriml solution of the Navier-Stokes
equations
187
Here, u is measured along the wall, X = 1 on the moving wall, and x = 0 on the fmed wall. The quantity {,n+l is found from (2.5) by using the second derivative of $ on the wall on the known time layer. On a half-integer layer n + 45,the quantity SWis found up to 0 (72) from { on the wall at the (n + I)-th, n-th, and (n - l)th, time layers. In accordance with (2.6), k=l, ,
(,qf,;)
n+‘=_
(&K)
n+‘!z,
k=l,
I,
k=J, j=l,
i=l, 2 ,...,I, 2,. . ., J.
It may be recalled that the second derivatives of the vorticity with respect to the coordinate were found on the wall in [4] by using, instead of (2.6), a similar corollary of the non-stationary equation. At the upper comers of the cavity, u breaks off, and hence I{ I becomes i&rite. In this connection the question arises of the influence of the method of completing the definition off at these points or in their neighbourhoods in the flow as a whole. Two ways of specifying { were tried. In the first U, in the celI adjacent to the comer point was smeared from 0 to 1; in the second, 5 at this point was replaced by the half-sum of 5 at adjacent points of the upper and side walls. In both cases, the difference between the calculated fields u and v was virtually everywhere a fraction of a percent. This is easily appreciated, if we note that, in accordance with the solution obtained in [ 141, which holds for small neighbourhoods of these points, their influence manifests itself at relative distances of the order of Re- 1 when the thickness of the wall layers is of the order of Re-‘A. Mention was earlier made in [9] of the small influence of the choice of the conditions at the upper corners.
FIG. l.-for N=51X51,A-for N=41X41, o-for N = 29X 29, 0 - for N = 40X40, from 1201
The choice of the time integration step is essentially dependent on the statement of the boundary conditions for {, especially for large Re numbers; this was revealed by the analysis made in [ 15-171. The statement of the boundary conditions employed in the present paper ensured stable computation over a wide range of Re numbers with a time integration step exceeding the period of the mesh.
188
V. L Kopchenov, A. N. Kraiko and M P. Levin
FIG. 2 Let us consider in more detail some aspects of the solution of the systems of difference equations. The equation for J/ is reducible to a system of scalar difference equations with a tridiagonal matrix of the coefficients with respect to the unknown nodal values of the stream function; these equations were solved by scalar pivotal condensation. The equation for { reduces to a system of vector three-point difference equations with respect to the unknown nodal values {ii and the derivatives of Mijor Lij. These were solved, either by Seidel’s method, or by vector pivotal condensation [ 18,191. On uniform meshes, cubic splines provided fourth-order approximation for the first derivative of the function with respect to the coordinate, while the second derivatives are approximated to second-order regardless of the mesh uniformity [4]. Though the order of approximation of the first derivatives drops to the third on non-uniform meshes, this is still high enough for the convective terms of the initial system of equations. 3. Using the method described, we wrote a program in FORTRAN for the BESM-6 computer. Some results, obtained on uniform meshes for Re = 100 and 400, are illustrated in Figs. 1 and 2. For these Re numbers, the computations were made with different numbers of computational nodes N, and were compared with the data in [20]. In Fig. 1 we show the profiles u = U(Y) on the centre line x = 0.5 for Re = 400, obtained for different N, and the similar profile in [20]. There is clearly a marked difference between our curves and that in [20] , though for Re = 100 the difference between the profile in [20] and that computed for N = 29X 29 is within the error of graphical representation. The difference with Re = 400 is primarily due to the inaccuracy of the data in [20]. This is revealed by Fig. 2, a and b, referring to Re = 100 and 400, where I$ Im is the maximum modulus of $, and h is the mesh step (1 refers to our data, and 2 to the data in [20]). The right-hand branches of the curves in Fig. 2, a are plotted, using points outside the graph. The method based on spline approximation clearly has faster convergence with respect to h than the method in [20]. This is particularly noticeable for Re = 400. For the case Re = 100 on a 21 X 21 uniform mesh we compared the efficiencies of the two methods mentioned above for solving the equations for t, and we studied the domain of stability with respect to T. It was found that, for Seidel’s method, the maximum tolerance half-step rm * 0.9h is the same as the optimal half-step (in the sense of computing time till establishment). For vector pivotal condensation, the maximum tolerable rm = 1.9/z, while the optimal is coopt= 1.3h. By comparison with the algorithm based on Seidel’s method, the computing time with vector
Numerical solution of the Navier-Stokes equations
189
pivotal condensation was reduced by a factor of 5 for r = 0.9h, and by a factor of 7 for r = I.3h In the latter case the computation took 5 minutes till establishment. These results show that, for the difference schemes employed, even with Re = 400 the uniform meshes enable reasonably high accuracy to be attained only for N > 40X 40. To progress to greater Re without seriously increasing N, non-uniform meshes must be used, in accordance with what has been said above. Let us describe how we constructed the mesh in the present paper, on the basis of special methodic computations and the considerations given above. In the longitudinal and transverse directions, the mesh was constructed according to the same principle, and was characterized by its minimum (hmo and h,) and maximum hm steps; the values of hmo were different for the lower and the other (hm) walls (hm < h, 0). Like the thicknesses 6 and S 0 of the corresponding wall layers, they were found from computations with Re = 400 and were reduced in proportion to Re-s as Re increased. The hm step in the core was taken to be virtually the same as for Re = 400. The transition from hm or hm 0 to hm was made outside the wall layers according to a geometric progression with “ratio” 4 = 1.4 - 2. Since it is made very rapidly, the total number N of cells scarcely increases as Re increases. Computations showed that, to obtain stable computation within the upper and side wall layers, several (in our case four) mesh layers have to be introduced, with cell size (along the normal to the wall) h, + < h, . The value of hm + was found from the expression h,+ = lO/Re; analysis of the simplified equations showed that this device is apparently connected with the good conditionality of the difference boundary value problem. The transition from h, + to hm, is made according to a geometric progression with the same ratio. The introduction of extra layers only involves a slight increase of N as Re increases.
-0.5
0
a.5
FIG. 3 The possibilities of non-uniform meshes are demonstrated by the following example. The flow at Re = 400 was computed on a non-uniform 27X 27 mesh with h,+=h,=0.02, hm0=0.035 and hm = 0.065 with (I = 1.4. It has been seen from Fig. 1 that the uniform mesh with roughly the same number of nodes gives extremely crude results. On the other hand, the velocity profile u = u (Y) on the line x = 0.5, obtained for a non-uniform division, is virtually the same as the profile computed on the much finer uniform mesh with N = 5 1 X 5 1. The practical value of the above principles for constructing the non-uniform mesh was confirmed, not only by computations for Re = 400, but also by extra data obtained for Re = 103. In the latter case, two meshes were used. Mesh 1 was constructed in the way explained above. Mesh 2 had twice as many nodes in the wall layers. Mesh 1 contained 39X 34 nodes with h,+=0.01, h,=O.O12, h,0=0.022 and hm = 0.053. Mesh 2 had more nodes, N = 49X 5 1, though the maximum step in
190
V. I. Kopchenov, A. N. Kraiko and M I! Levin
FIG. 4
L.
-n.i5
n
0.5
FIG. 5 the core remained roughly the same. Doubling of the number of nodes in the wall layer did not produce a serious correction of the results, notably for the velocity I( profile at the cavity centre. The accuracy of the results can be judged from the following data. The values of the maximum speed of reverse current 1r.P I, obtained on meshes 1 and 2, are respectively equal to 0.392 and 0.382. In Fig, 3 the present results (continuous curve) are compared with the data in [21] on the u profile in the cross-section passing through the vortex centre (by the vortex centre we mean the point where I $ I takes its maximum value). The broken curve of Fig. 3 shows the results obtained in [21] using a non-uniform mesh with N = 17 X 17 and a spline scheme of higher order than used in the present paper. The dot-dash curve was computed in [21] by a finite-difference scheme with uniform division and the same N = 17 X 17. The analogous results in [21] for a uniform mesh with N= 65 X 65 agree (within the scale of the graph) in different sections either with the continuous or the broken curve. In Fig. 4 we show the streamline picture with Re = 1000 (the figures on some lines refer to the values of $).
Numerical solution of the Navier-Stokes equations
191
FIG. 6 The data quoted show the advantage of using a non-uniform mesh. As Re increases the advantage becomes more obvious. In principle, it seems that this device could be used to move up to the limit of existence of stationary modes. While the establishment time naturally increases with Re (as the viscosity falls, the “acceleration” of the flow requires more “physical” time), the only serious obstacle to advancing towards higher “laminar” Re would seem to be the non-stationary property of the flow. While admittedly the present algorithm does not describe the actual non-stationary flow, due to the “stationary” boundary conditions for c, and the artificial term a$/& added to the first of Eqs. (2.1), nevertheless, establishment took place in the algorithm at least up to Re = 3.5X 103. Some results for this case are shown in Figs. 5,6, where we give the profile u CJ) for x = 0.5 and the and km = streamline picture for a 47 X 45 mesh with h,+ =2.5. 10m3, h,=0.007, h,“=0.01 0.052. Let us conclude with some remarks on the construction of non-uniform meshes. In the present paper the mesh was fmed and was specified at the start of the computation, on the basis of an analysis of the results for relatively low Reynolds numbers. Since this approach can often be very laborious, it is very important to examine the automatic construction of computational meshes, ensuring satisfactory resolution of the zones of large gradients. Of work in this field, we mention [l] and [22], in which it is proposed to use coordinate transformations dependent on the velocity components. In [23], when solving the stationary problem by the establishment method, the mesh is constructed on each time step by minimizing a certain functional. The method described in [24]’ has become quite popular; here, the node coordinates are found by solving a system of two quasi-linear inhomogeneous equations of elliptic type [25]. In this latter case, however, the disposition of the nodes within the computational domain can only be controlled by compressing them on the boundaries, which usually demands preliminary information about the disposition of the zones of large gradients. In the authors’ view, genuine scope is offered by the use of a dialogue type of
192
V: L Kopchenov, A. N. Kraiko and M P. Levin
interaction with the computer, in which the mesh is corrected during the computing process according to the parameter distributions that are realized. Translated by D. E. Brown
REFERENCES 1. TOLSTYKH, A. I., A method for the numerical solution of the Navier-Stokes equations of a compressible gas over wide range of Reynolds numbers, DokL Akud. Nauk SSSR, 210, No. 1,48-51, 1973. 2. DAIKOVSKII, A. G., POLEZHAEV, V. I., and FEDOSEEV, A. I., Study of the structure of transient and turbulent modes of convection in a vertical layer, Izv. Akad. Nauk SSSR, Mekhan zhidkosti gaza, No. 6, 66-75,1978. 3. DAIKOVSKII, A. G., POLEZHAEV, V. I., and FEDOSEEV, A. I., Numerical modelling of transient and turbulent convection modes using the non-stationary Navier-Stokes equations, Preprint IPMekhan. Akud Nauk SSSR, No. 101, Moscow, 1978. 4. RUBIN, S. G., and GRAVES, R. A., Viscous flow solutions with a cubic spline approximation, Fluids, 3, No. 1, l-36, 1975.
Cbmput.
5. RUBIN, S. G., and KHOSLA, P. K., Higher-order numerical solutions using cubic splines, AIAA J., 14, No. 7, 851-858,1976. 6. LEWIS, E., Steady flow between a rotating circular cylinder and a fixed square cylinder, J, Fluid Mech., 95, No. 3,497-514,1979. 7. KUSKOVA, T. V., Numerical study of two-dimensional flows of a viscous incompressible fluid, in: Applications ofrhe mesh method in gas dynamics (Nekotorye primeneniya metoda setok v gazovoi dinamike), No. 3, 7-136, VTs MGU, Moscow, 1971. 8. BOZEMAN, J. D., and DALTON, C., Numerical study of viscous flow in a cavity, J. Comput. Phys., 12, No. 2, 348-363,1973. 9. GHIA, K. N., HANKY, V. L. Jr., and HODGE, J. K., Solution of Navier-Stokes equations allowing for incompressibility in the usual variables, in: Rocket techniques and cosmonautics (Russian translation, Mir, Moscow, 1979). 10. HOSMEN, A. D., et aL, Numerical methods of studying viscousjluidflows 1972).
(Russiantranslation,hfir,Moscow,
11. DE VAHL, D. G., and MALLINSON, G. D., An evaluation of upwind and central difference approximations by a study of recirculating flow, Comput. Fluids, 4, No. 1,29-43,1976. 12. KURTZ, L. A., et aL, A comparison of the method of lines to finite difference techniques in solving timedependent partial differential equations, Comput. Fluids, 6, No. 1,49-70, 1978. 13. AHLBERG, J. H., et aL, Theory of Splines and Their Applications, Academic Press, New York, 1967. 14. BATCHELOR, G. K., Introduction to Fluid Dynamics, Cambridge U.P., 1967. 15. DAIKOVSKII, A. G., POLEZHAEV, V. I., and FEDOSEEV, A. I, On calculation of the boundary conditions for non-stationary Navier-Stokes equations in “vorticity, stream function” variables, in: Numerical methods of the mechanics of a continuous medium (Chisl. metody mekhan. sploshnoi sredy), Vol. 9, No. 7,49-58, Nauka, Novosibirsk, 1978. 16. TARUNIN, E. L., On the choice of approximation formula for the velocity eddy at a rigid boundary when solving problems of viscous fluid dynamics, in: Numerical methods of mechanics of continuous medium (ChisL metody mekhan. sploshnoi sredy), VoL 9, No. 7,97-111, Nauka, Novosibirsk, 1978. 17. PEPPER and HARRIS, Numerical modelling of free convection in closed cavities by completely implicit method, in: Theoretical foundations of engineering calcuhtions, (Russian translation, Mir, Moscow, 1977). 18. VASIN, V. G., and POLEZHAEV, V. I., Difference schemes for Navier-Stokes equations of compressible gas and calculation of thermal-acoustic waves, Preprint IPMatem. Akad. Nauk SSSR, No. 124, Moscow, 1977.
Instability of a convergent spheriml shock wave
193
19. SAMARSKII, A. A., and NIKOLAEV, E. S., Methods for solving mesh equations (Metody resheniya setochnykh uravnenii), Nauka, Moscow, 1979. 20. BURGCRAF, 0. R., Analytic and numerical study of the structure of stalling flows, in: Collection of surveys and translations of foreign periodical literature, Mechanics, No. 6, (IOO), 51-90, Mir, Moscow,1966. 21. RUBIN, S. G., and KHOSLA, P. K., Polynomial interpolation J. Comput. Phys, 24,No. 3,217-244, 1977.
methods for viscous flow calculations:
22. TOLSTYKH, A. I., On the condensation of nodes of difference methods in the process of solution and on application of higher-accuracy schemes for the numerical study of viscous gas flows, Zh. vychisl Mar. mat. Fiz., 18, No. 1, 139-153, 1978. 23. DANAEV, N. T., LISEIKIN, V. D., and YANENKO, N. N., On the numerical solution of viscous thermally conducting gas flow past a solid of revolution on a moving curvilinear mesh, in: Numerical methods of the mechanics of a continuous medium (ChisL metody mekhan. sploshnoi sredy), Vol. 11, No. 2, 51-61, Nauka, Novosibirsk, 1980. 24. THOMPSON, J. F., THAMES, F. C., and MASTIN, W., Automatic numerical generation of a body-fitted curvilinear system for a field containing any number of arbitrary 2-D bodies, J. Comput. Phys., 15, No. 3,299-319,1974. 25. THOMAS, P. D., and MIDDLECOFF, J. F., Direct control of the grid point distribution in meshes generated by elliptic equations, AIAA J. 18, No. 6,652-656, 1980.
U.S.S.R Cornput. Maths Math. Phys. Vol. 22, No. 6, pp. 193-205.1982. Printed in Great Britain
0041-5553/82$07.50+.00 01984. Pergamon Press Ltd.
INSTABILITY OF A CONVERGENT SPHERICAL SHOCK WAVE* K. V. BRUSHLINSKII Moscow (Received
8 December
1980)
A CONVERGENT spherical shock wave is unstable in the linear approximation. Arbitrary small disturbances of the flow can be expanded in series in generalized spherical functions. The eigenvalues of the differential operators, characterizing the growth rate of the disturbances in the case of the first few spherical harmonics, are found numerically.
Introduction The convergence of a spherical shock wave to its centre of symmetry has been studied by Guderley [ 11, Landau and Stanyukovich [2], and Gel’fand et a~!(see survey [3]), and is a familiar example of a similarity problem of gas dynamics. In this problem the power k in the similarity variable { = tr-k is transcendental and can only be determined during numerical solution of the boundary value problem for an ordinary differential equation. This creates special problems when studying the shock wave stability. Particular mention must be made of [4], where, under simplifying assumptions (channel approximation), an analytical solution is obtained for the convergent shockwave problem, and the instability of the wave is proved. In more complete studies, however, there is necessarily an element of numerical analysis, if only because the similarity solution itself is found numerically. *Zh. vychisl. Mat. mat. Fiz., 22,6, 1468-1479, “bSK.‘.‘:h-:,
1982.