Introduction

In this paper we present a numerical scheme of new type which is designed for fast calculation of non-stationary boundary value problems for hyperbolic systems of PDE’s. We call this scheme ‘explicit-implicit’ type scheme. The main feature of it is the ability to efficiently calculate problems with different scales. Namely, while in some regions (or along some directions, etc.) the flow is rather regular large time step can be used and one can use explicit scheme (which is preferable because of higher resolution, better stability properties, etc.), in other regions the flow can exhibit non-regular behavior with complicated fine structure and require fine greed and, accordingly, a small time step for explicit scheme to be used. In this case one uses in our scheme the implicit version and retains a rather large time step. Such methodology makes it possible to efficiently calculate problems with different scales as a whole. Generally speaking, we propose the family of schemes whose construction is based on the following principles:

-         use locally explicit scheme where it is possible;

-         use locally implicit scheme where the Courant number exceeds unity;

-         use upwind principle;

-         the switching between various local scheme types must be continuous;

-         use minimal possible stencil (depending on the specific problem).

Below we describe in detail a realization of the proposed methodology for the calculation of non-stationary processes in swirled flows in the nozzles with non-trivial geometry. We use Euler equations in cylindrical coordinates as a basic model. The swirling coupled with the longitudinal flow and geometry effects generates complicated non-stationary regimes with almost periodic behavior, possible bifurcations and flow restructuring. Besides, such problem has obviously at least two scales inherited from geometry: the longitudinal direction and the radial direction. In the radial direction we need a much finer greed to capture the flow features. Thus, the problem of swirled flow in tubes/nozzles is suitable for testing the ‘explicit-implicit’ technique. 

Our preliminary experience shows that the subsonic regime for swirled flows, generally speaking, is simpler than the supersonic one. Hence, we perform calculations in case of supersonic flow with respect to the longitudinal coordinate. The calculations reveal the oscillating shock wave upstream the diffuser and evolving vortexes in the diffuser. These vortexes can grow, push the shock wave and destroy it. Then shock wave appears again near the nozzle throat, comes closer to the diffuser and the process repeats.

Our paper is organized as follows. In Section 2 we rigorously formulate the problem and in Section 3 we describe in detail the realization of numerical scheme for the model equations according to proposed above methodology. Section 4 is devoted to such a realization for Euler equations while in Section 5 we concentrate on the peculiarities of the case of axisymmetric flows. Finally, in Section 6 we present the results of calculations and relevant discussion including the possibility to calculate Navier-Stokes equations with this technique.

Section 1. Formulation of the Problem

We consider non-stationary axially symmetric flow of swirled non-viscous gas in a nozzle, see Fig. 1.

Fig. 1 General scheme of calculation domain

The corresponding system of differential equations in cylindrical coordinates (taking into account the fact that the derivatives with respect to the angle vanish) can be presented in following two main forms:

1.     Divergent form.

                                                                                          (1)

where

                    (2)

2.     Non-divergent form.

                                                                                     (3)

where

                                                                 (4)

Remark: in some cases the following another divergent form is more convenient

                                                                                         (5)

where  is the vector column .

Here  denotes the density; – the pressure; ,  and  are the axial, radial and tangential velocity components respectively;  is the speed of sound; , ;  and  are the internal energy and enthalpy per unit volume.

We assume that the gas is perfect, i.e. , , .

The geometry of the problem: is the inlet section,  is the outlet section;  and  are the equations of the nozzle wall and the central body. We note that the central body can extend either along all the nozzle () or only along its front part () or can be absent ().

For the sake of definiteness of imposing the inlet boundary conditions, we suppose here that in some neighborhood of the inlet section both boundaries are cylindrical surfaces, i.e.  in the interval .

Boundary conditions:

1.     At the inlet section.

                                                                                            (6)

The first condition is due to the cylindrical geometry of the initial part of the nozzle while the second condition imitates the swirling generated by the swirling device of vane type. Namely, we suppose that the flow after passing each vane turns by the angle. The dependence of  on models the variability of the vanes angle.

As it is known, four boundary conditions are required at the subsonic inlet section (this case is considered in what follows). Thus, besides two kinematical conditions (6) it is necessary to specify two additional conditions. These conditions might be, for example, the determination of any pair of the following functions:,,, where  is the temperature. However, the dependence of these functions on  is unknown and the assumption that they are constant with respect to  is hardly true. Therefore, a more natural way is to determine the state of the gas at rest in the gas holder – the tank of big volume from which the gas goes to the vanes of the swirling device. Such state can be determined through two stagnation parameters, either,  or . Since the process of gas flowing from the gas holder to the vanes can be considered as stationary and adiabatic, we assume that at the inlet section the entropy and the Bernoulli integral  equal their values at the stagnation point, i.e. equal known values  and . Hence to two linear boundary conditions (6) one has to add two non-linear conditions

                                                                                         (7)

                                                                       (8)

(one takes into account (6) in boundary condition (8)).

2.     At the outlet section.

Here two cases are possible. The first case: the outlet section is supersonic. Then the boundary does not require any boundary conditions. The second case: subsonic outlet section. Then for problem to be well posed one needs to impose one (the positivity of  is assumed) boundary condition. Usually the constant pressure is set as such boundary condition:

                                                                                                           (9)

3.     At the nozzle wall and at the surface of central body.

Conventional non-penetration conditions are posed at the nozzle wall and at the surface of central body:

                                                              (10)

Here and below the prime denotes the differentiation.

Let us note that at the axis  (beyond the central body) no additional boundary conditions are required (more precisely, only the boundedness of the solution is required). Because of the axial symmetry one has for

                                                                     (11)

These relations can be used in numerical algorithm to avoid the loss of approximation, which is connected with the coordinate singularity.

The geometry of the problem allows easy transformation to the simple calculation domain. This procedure is done in conventional way – introduction of normalizing coordinate:

                                                                                  (12)

Now the calculation domain is rectangle.

In coordinates  system (1) takes the form

                                                                                         (13)

Here

                                                           (14)

For the non-divergent form (3) one has

                                                                                     (15)

where

                                                                                      (16)

Thus we consider the solution of initial-boundary value problem for the system (13) (or (15)) in the domain  with boundary conditions (6) – (10).

Section 2. Basic Principles of the Difference Scheme

The specific feature of problem under consideration is the combination of swirling and the elongated geometry of the nozzle. With no swirling the flow is practically one-dimensional and there is no need to use a fine grid with respect to  (to). Hence, a difference scheme explicit with respect to - direction provides sufficient accuracy under acceptable boundedness of the time step. As far as longitudinal coordinate  is concerned, the complex and non-stationary structure of the flow dictates the choice of explicit scheme. Thus a fully explicit scheme makes it possible to achieve the required accuracy for reasonable time for the numerical solution of the problem of the gas flow without the swirling.

In the presence of swirling the picture is totally different. In this case the radial gradients are much times greater then longitudinal gradients. So, to preserve the approximation it is necessary to use a very fine mesh with respect toand this fact makes explicit schemes practically inoperative.

On the other hand, the practice of numerical calculation shows that implicit schemes have lower resolution than explicit schemes and hence it is desirable to reduce the scope of their usage.

The acceptable compromise can be achieved following the principle: use the implicit scheme in those and only those situations when the modulus of corresponding Courant numberexceeds unity. The schemes satisfying this condition can be called ‘locally implicit’. It is clear that such schemes are nonlinear and their structure depends on local properties of the solutions. As far as the longitudinal coordinate  is concerned the structure of the flow defines the choice of the explicit scheme.

To summarize, the desired scheme has to have the following properties:

1.     To be explicit with respect to the - direction;

2.     To be locally implicit with respect to - direction, i.e. depending on the value of Courant number to become either explicit or implicit;

3.     The switching of such scheme should be continuous.

The experience of usage of various locally implicit schemes has confirmed their efficiency [1], [2].

Let us start to discuss the main idea of locally implicit schemes. In present study we deal with the schemes of the first order. We first consider the one-dimensional case.

Consider the equation

                                                              (17)

Following the methodology of concept of finite volumes we split the interval (0,1) into subintervals of length. The subintervals’ boundaries are. We will refer the values of mesh functions at upper and lower levels to the centers of these intervals, i.e. to the points . Let us introduce the fluxes (more exactly the fluxes’ densities). Then

                                                                                     (18)

where

This closing stage is common for all methods based on the notion of finite volumes. Thus, the method of calculation of is the critical stage.

Multiplying (17) by, we obtain

                                                                                                  (19)

Precisely this form dictates the construction of both schemes.

Scheme 1 (explicit)

                                                                              (20)

This scheme of fluxes calculation complemented by the divergent closure (18) generates the simplest upwind difference scheme (see Fig. 2)

                                                                  (21)

where .

The condition for its stability and monotonicity is.

Scheme 2 (implicit)

                                                        (22)

Excluding from (18) and (22) we obtain

                                                            (23)

(implicit analogue of the scheme (21), see Fig. 3).

It is easy to see that the scheme 2 is stable (and also monotonic) only when. Let us emphacise that the schemes 1 and 2 coincide for.

It is convenient to represent both schemes for calculation of  as formally three-point difference equation

                                                            (24)

where  depend on  in the following way

                                                       (25)

Let us note that there is an ambiguity forwhen. In this case one can define.

Now let  and, accordingly, . Then the coefficients in (24) are calculated in accordance with (25) taking into account the variability of: is substituted by. The same scheme is valid also for. In this case.

In a number of cases it is technically easier to use the same construction, but with respect to, rather than to. Then in (20), (22), (24) and (25) andare simply substituted byandwith the same indexes. The solution of three-point difference equation gives, after that is calculated. Such a simplification is admittable only for continuous solutions.

It is obvious that for any smooth  the boundary value problem for the system in finite differences (24) is well-posed provided the determination of boundary conditions for (17) corresponds to the signs of  at the boundaries.

The same idea is applicable also for one-dimensional hyperbolic system: now,  is diagonalizable matrix with real eigenvalues.

Suppose is the -th left eigenvector of the matrix and is the corresponding eigenvalue. Denote by the matrix with rows. Then where is the diagonal matrix with elements.

We transform (19) to the diagonal form. Then we have

                                                                                       (26)

It now only remains to apply the already known construction to each component of system (26), i.e. to the equations

                                                                     (27)

Actually it means that the matrixes and vector in (24) are filled with rows. The corresponding table for -th row has the form (see (25)):

                                                (28)

Here.

Now let us consider two-dimensional analog of equation (17):

                                                                           (29)

We start with the case of the constant coefficients  and and try to construct a difference scheme explicit with respect to and locally implicit with respect to  by means of the generalization of the above scheme. Now the difference mesh is two-dimensional –,

. The mesh function is, the fluxes are .

Suppose that. Then the simplest two-dimensional analog of presented locally implicit scheme should look as follows:

                                                                           (30)

The closure stage is standard

                                                                        (31)

Here, , .

It is easy to check that the stability condition for fully explicit variant of the scheme (30), (31) has the form

                                                                                                     (32)

It is clear that if, then this scheme is unstable for any and a smooth switching between the schemes is impossible. Thus the simplest and seemingly most natural construction shows itself invalid. One should transform it in such a way that the stability conditions (for positive) take the form:

Scheme 1 (fully explicit) –;

Scheme 2 (– explicit, – implicit) –.

Then their unification gives us the desired locally implicit scheme that has the approximation property.

The realization of these schemes takes place in two steps.

First step is the calculation of intermediate values of  and which we denote by tilde  and.

Scheme 1:

                                                                                                  (33)

.                                                                                                  (34)

Scheme 2 differs from Scheme 1 by the substitution instead of (34)

                                                                        (35)

The divergent closure completes the step

                                                                       (36)

Second step is the ultimate calculation of  and values. First, weighted value of the solution is calculated:

                                                                                     (37)

Then for Scheme 1

                                                                                        (38)

                                                                                       (39)

For the Scheme 2 (39) is substituted with

                                                                        (40)

Finally, divergent closure (31) gives ultimate result.

Thus, second step differs from the first one by the substitution of with ‘weighted’ value when we calculate (but not in closure equation!).

Let us use conventional technique to obtain necessary stability conditions. Suppose. Then

,                                                            (41)

where according to (33), (34) and (36)

                                                             (42)

It follows from (35) that

                                                                                (43)

where

                                                                                           (44)

                                                            (45)

Second stage gives

                                                            (46)

It is easy to see that

                                                           (47)

And finally

                                                                              (48)

This scheme should be stable for any positive values of  and for  belonging to the segment. In particular, this should be true in case when, i.e. both schemes coincide. Then. The stability condition is  for all. Taking  one has  and hence. So,  if  and the scheme is unstable.

For, according to (48), . Since, then. It is easy to demonstrate that  under the condition. Hence, , and we prove the stability of the scheme in this case.

Other combinations of signs of  and  lead to the same result.

The generalization to the case of hyperbolic system is fulfilled by the analogy with one-dimensional problem because each stage is actually one-dimensional.

The proof of the stability of the scheme for symmetric matrixes  is again performed by the generalization of above method. Since is now a vector, then, where is constant vector. Then

                                                        (49)

                                                                      (50)

Let transform  to diagonal form

                                                                   (51)

where are diagonal matrixes with the elements and are orthonormal matrixes.

The substitution of (49) into the formulas from the first stage gives

                                                                                      (52)

(compare with (43)), where

                                                                       (53)

Here are diagonal matrixes with elements that are determined via (44), (45), and therefore lie within the unit circle.

Performing averaging with and substituting (50) into the formulas of the second stage, one finds that

                                                                                              (54)

Since , it follows from (53) that , which proves the stability.

It is clear that the scheme is stable and for the systems reducible to symmetric ones, in particular, for linearised Euler equations.

We’ve presented above only the basic ideas of the construction of locally implicit schemes. Each specific realization certainly contains a number of peculiar elements, which take into account the essentials of the problem under consideration.

Section 3. The Difference Scheme For Euler Equations

Proposed algorithm is using both forms of main system: non-divergent (15) as well as divergent (13). The basic element of the algorithm is the calculation of fluxes  and  at lateral sides of calculation cell, i.e. .

Fig. 2 Calculation cell

Let us start with -direction, for which the fluxes can be calculated with respect to the explicit scheme. The construction of this scheme uses left eigenvectors  of Jacobi matrix  and its eigenvalues. Suppose  is the left eigenvector of matrix: . Since  and , then .

It is easy to see that

                                                                 (55)

                                                              (56)

The matrixes  and have the form

                                                                       (57)

                               (58)

where .

It follows from these expressions that the following vectors can be taken as

                                                            (59)

Let us determine unknown fluxes, taking into account the analogy to the scalar case, from the following one-dimensional system

                                                                                              (60)

Let us freeze matrix coefficients and after transforming (60) to diagonal form one get for the invariants the following system

                                                                                             (61)

Then on the basis of the above one infers the algorithm for calculation of :

1.     Calculation of  and  at lateral side with the aid of simple averaging with respect to the points  and ;

2.     Calculation of  provided  for any .

3.     Determination of  from the system of linear equations,

                                               (62)

Since in the flows under consideration strong discontinuities (shocks) arise only in z-direction, this situation has also to be taken into account. For the change of the sign of  do not lead to unphysical solutions, the algorithm includes the analysis of  signs in points neighboring (with respect to) to given facet. In this case, if or , then corresponding equation in system (62) has to be changed using

                                                          (63)

For the inlet section four equations of system (62), which correspond to positive, are substituted with boundary conditions. For subsonic outlet last equation from (62) is substituted with boundary condition. In case of supersonic outlet the right hand side in (62) is fully determined (because of positiveness of all  at such boundary).

The eigenvectors, which are necessary for the calculation of boundary fluxes, are calculated not by averaging with respect to neighboring points (it is impossible), but by simple extrapolation or even transfer of  values from neighboring inner point.

Now let us look at -direction. Since in present problem the cross-sectional shocks are absent, the calculation of corresponding fluxes  can be performed using non-divergent form (15) of the main system. This fact essentially simplifies the calculation formulas. According to (16) the matrix  has the form

                                  (64)

Hence for the eigenvalues of the matrix  one has

                                             (65)

where.

Corresponding left eigenvectors  are

                                                                                      (66)

Here.

Since now basic vector is not, but (see (4)), the analog of (61) will be

                                                                              (67)

where

                                                                                                       (68)

Locally implicit (with respect to -direction) scheme is being constructed by natural generalization of scalar case (17). This approach leads to difference system of type (22), where now are  matrixes:

                                                     (69)

where is vector column with the components .

Each equation (row) of this system is the difference approximation of appropriate equation of the system (67) (index ‘’ is omitted for the sake of simplification of writings):

                                    (70)

Here

                             (71)

Courant numbers and invariants are calculated at lower layer. All values that are rendered to lateral sides of calculation cell, i.e. have fractional index with respect to, are calculated by the averaged values.

The system (69) is closed by boundary non-penetration conditions at both boundaries  and:

                                                                                                 (72)

which can be expressed through the invariants and as

                                                                                                     (73)

(Let us remind that.)

Since at both these boundaries, then under reasonable choice of time step the value also at the points neighboring to the boundary points. Hence, at the boundaries and at adjacent to the boundaries facets first three invariants are calculated with respect to the explicit scheme. And so actually the system (69) contains an isolated subsystem of the same type, but for 3-vectors. The solution of this subsystem is easy to find with the aid of conventional sweep method.

As far as the rest of invariants is concerned, these invariants can be found from corresponding subsystem (69) for 2-vectors. In this subsystem at each boundary the difference equation, which corresponds to the coming into the domain characteristics (i.e. for  – forth equation, and for – fifth equation), has to be substituted with the boundary condition in the form (73). Then at the lower boundary the matrix coefficient (more exactly, the part of this coefficient corresponding to 2-vector) turns to zero, and at the left boundary so does the matrix coefficient. Thus resulting subsystem has the standard three diagonal structure and conventional sweep method can be used for its solution.

The solution of both subsystems gives at all facets the complete set of invariants. After that it remains to find at the same facets from the system

                                                                                (74)

where subscript ‘’ denotes that corresponding values are calculated by the aid of averaging with respect to the values at neighboring integer points of lower layer (analogously to calculation of  in (71)).

Further, is calculated knowing, and finally –. Here the procedure of fluxes calculation is over. Let us note that as  and as. Hence, the value of at the boundaries does not depend on the invariants.

Described above method of fluxes (according to explicit scheme) and (according to locally implicit scheme) calculation constitutes the basis of algorithm of the transfer to the next layer: . As in the scalar case (see Section 2) this procedure is fulfilled in two stages.

First stage is completed by the calculation of intermediate values. To do this one has first to find by the aid of divergent closure, i.e. difference form of system (13):

.                                                              (75)

(For divergent closure it is possible to use other form of equations – analog of (5). Numerical experiments show that in some cases this form is preferable.)

Then one finds using the values of:

                                                          (76)

The second stage differs from the first one only by substitution of, which is involved in calculation of and, with weighted value. Taking into account that the stability of scalar scheme is achieved only as, we have

.                                                                                              (77)

Let us emphasize that we speak only about the calculation of and. Under the divergent closure process the values are not being modified. The result of second stage is desired mesh vector function.

Section 4. The Singularities of Swirled Axisymmetric Flows

Since considered flow does not depend on angular coordinate, the forth equation of system (1) gives

                                                                           (78)

where. This equation expresses the conservation of circulation along the streamline. Under the streamline one understands its trace on the plane. (Real streamline is the spiral one.)

Let us look at some consequences of the law of circulation conservation, and start from stationary flow in the nozzle with central body, see Fig. 1. Suppose at inlet section in the point, which is situated on the surface (generatrix) of central body. Moving forward along the particle goes down to the axis, i.e. its -coordinate decreases. But for this streamline. So increases along this line. But full velocity is bounded from above because the Bernoulli integral is constant. Hence,, where is the value of in gas holder. So this streamline can not reach the axis: . Thus some point on should exist where the flow is singular. For example, in such a point the flow detachment – the line of tangential discontinuity – can form. In this case under such line the stagnation zone should arise.

In non-stationary case the matters do not change because unbounded growth of leads to unbounded growth of, i.e. finally to the formation of singularity. In general the law of circulation conservation forbids any streamline with approaching -axis. This fact actually makes ill-posed the problem of inviscid flow of swirled gas in the nozzle with central body of finite length. The location of detachment point and the structure of the flow strictly speaking essentially depend on the initial data.

One can overcome this difficulty in two ways. First, it is possible to set such inlet conditions that the swirling on vanishes. In practice such requirement can be fulfilled using the variability of vanes angles in swirling device. Second, one can introduce thin tube, for example, cylindrical one of radius, which forbids the streamline coming to -axis. Of course the thickness of this tube should be sufficiently small for not having significant influence on the main part of the flow.

To estimate the sizes (with respect to) of singular domain let us consider the following model problem. Suppose the flow is isentropic (), isoenergetic (), stationary and does not depend on longitudinal coordinate. Then and all other functions, which depend only on, satisfy the following system of equations

                                                                                          (79)

                                                              (80)

                                                                                                     (81)

Introducing, rewrite (80) in the form

                                                                           (82)

It follows from (79) and (81) that

                                                                                              (83)

The exclusion of from (82) and (83) gives

                                                                                             (84)

where

                                                                                 (85)

Solving (84) by standard methods one has

                                                                                         (86)

where is an arbitrary constant.

Suppose. Denoting corresponding value of through , one infers

                                                                                            (87)

Thus,

.                                                                                 (88)

For the longitudinal Mach number one has

.                                                                                      (89)

And finally

                                                                               (90)

Hence, considered flow exists only beyond -vicinity of -axis, i.e. only for. As , the values of as well as tend to zero, while the values of and tend to infinity.

It follows from (90) that

                                                                                                    (91)

which is in fact desired estimate of . The variability of  does not lead to principal change in the flow picture, but only complicates the formulas.

Section 5. The Results of Calculations and Discussion

The calculations are done for various forms of the channel and various flow parameters, but the calculations scenario is the same for all cases: first stationary flow without swirling has been calculated (in formula (6) equals zero), then instantly or for short time period the swirling vanes acquire given angle. At this the channel geometry and flow parameters are such that in minimal (critical) section the sonic velocity is achieved for stationary flow and downstream the flow is supersonic at least in the part of the channel. This allows separating up to the proper degree the investigation of the flows in the premix chamber (the domain from the end of central body to critical section) and downstream.

Let us now touch the point about the role of central body. A number of numerical experiments produced with the presented algorithm confirm that stationary flows in the nozzles and pipes without central body, even in the presence of swirling, exist and can be obtained with the aid of time relaxation method without special tricks. Required relaxation time equals approximately to the time interval for which the perturbations travel 3 – 5 times along the calculation domain (upstream and downstream). The same is true for the flows in the nozzles and pipes with central body, but without the swirling.

However, the presence of any free swirling at the nozzle inlet and as a consequence at the surface of central body drastically changes the flow picture. In particular, stationary flow becomes impossible. Indeed, the value Ω=rw is transferred along the streamline (see equation (78)), and following the arguments of Section 4, the streamline must detach from the central body at some point. Moreover, by the very same reason this streamline cannot come to the axisdownstream. Hence, one infers that in the neighborhood of axis the stagnation zone, is formed. This fact contradicts the structure of the stationary flow near the critical section.

It follows from the above that in the presence of swirling after some time the detached flow is formed at some point of central body and then weak reverse flow is also formed. Thus, the contact discontinuity arises. Further evolution of the flow in the neighborhood of -axis is determined by the numerical approximation scheme, in particular by scheme viscosity.

Yet first calculations of the flow in the premix chamber confirm the theoretical derivations that the flow structure depends essentially on the swirling character. In case if one determines the inlet swirling by the formula, where  (here  is the radius of central body, is the outer radius of the channel at the inlet),  is longitudinal and  is circumferential components of the velocity, the singular point does not arise at the surface of central body, and the flow is rapidly stabilized. Typical distribution of circumferential velocity in the middle of the premix chamber is shown in Fig. 3.

Fig. 3 The profile  in the section 140 mm for the calculation with linearly growing swirling; grid 175х40 (see geometry below in Fig. 5)

If one poses at the inlet, where thus imitating the flow passage through the vanes with turning angle , the flow character changes. The flow rather rapidly stabilizes visually, but then slow quasi-stationary regime starts when thin axial jet of non-swirling flow being very slowly swept out. As it is said above, the singular point arises at the descending surface of central body. This point separates the flow in the premix chamber by inner and outer zones that can be discerned by different distributions of  (see next Fig. 4). In this figure the profiles of circumferential velocity component in the middle section of the premix chamber for three grids at the same moment of time are shown. It can be seen that, first, scheme viscosity allows the swirling to penetrate to the domain below the contact surface and, second, scheme viscosity becomes weaker with the growth of grid size and this leads to the increase of maximal value of w. If one does not use special means further under sufficiently fine grid such configuration destroys itself (either pressure becomes negative or oscillations arise, which finally lead also to negative pressure or/and density).

This fact lies in the very nature of the problem. By the formation of contact surface the flow is separated in two immixing parts: upper ‘active’ flow and axial stagnation zone. The evolution of stagnation zone is mute and, seems, is determined by the factors that are external to the problem under consideration.

Fig.4 The comparison of profiles  in the section 140 mm for three grids: 175х20 (…), 350х40 (ооо) and 350х80 () (see geometry below in Fig. 5)

Because of this reason all other calculations are performed in modified domain as is theoretically suggested in Section 4: a thin tube (needle), which separates the flow from the axis r=0, has been located in the flow domain. The needle radius has been determined in two ways. Either, preliminary calculation without the needle has been performed, the detachment point at the surface of central body has been determined and then main calculation is fulfilled with the needle that contains such detachment point; or the needle radius has been estimated by formula (91). Both methods give close results. The intensity of perturbations that are generated by the needle depends on the parameters of the problem under consideration: the greater swirling entails the greater thickness of the needle. If the perturbations become too large the problem has to be posed in another way (to take viscosity into account, to allow the reverse flow at outlet boundary, etc.) In present calculations the area of the section of the needle is approximately 1% of critical section area.

To estimate the accuracy of the method the series of calculations is performed for the channel with central body, see the geometry in Fig. 5. The turning angle of the vanes is 45º, needle thickness is 0.001 (critical section radius is 0.01). Fig. 5 also shows the isolines of pressure that is scaled by the stagnation pressure p₀. One can see sharp pressure drop towards the nozzle axis. Such pressure drop is the characteristic feature of the swirled flows in the nozzles and pipes. Figures 6 – 8 demonstrate the profile of p(r)/p₀ , and the profiles of u(r), w(r) that are scaled by the local speed of sound. Shown data are referred to the section z=0.180 and obtained for the same configuration but for various grids: 55x30, 110x60, 220x120 and 440x240. The non-uniform convergence of the solution with respect to the grid size can be seen: while number of mesh points increases the solution stays practically unchanged near the nozzle wall but converges rather slowly in the zone near the axis. It seems that this fact can be explained by the sharp change of values in the neighborhood of the needle.

Fig. 5 Isolines of pressure under turning angle of the vanes 45º (needle thickness is 0.001)

Fig. 6 Profiles p(r) as z = 0.180 for four grids: circles – 55x30 points; dot-and-dash – 110x60; dashed line – 220x120; solid line – 440x240.

Fig.7 Profiles u(r), scaled with respect to local sound speed, as =0.180 for four grids: circles – 55x30 points; dot-and-dash – 110x60; dashed line – 220x120; solid line – 440x240.

Fig.8 Profiles w(r), scaled with respect to local sound speed, as =0.180 for four grids: circles – 55x30 points; dot-and-dash – 110x60; dashed line – 220x120; solid line – 440x240.

Moreover, two points deserve close attention. First, one can observe the presence of significant high-speed jet in the neighborhood of the axis. The longitudinal component of the velocity in the jet is two times greater than in the main stream. Second, ‘angular’ speed becomes supersonic in some neighborhood of the axis before the flow passes the critical section. Figures 9 – 11 demonstrate analogues profiles in the section =0.420 where the longitudinal component of the velocity is already supersonic. The abovementioned features of the flow take place in this case also.

Fig. 9 Profiles p(r) as z = 0.420 for four grids: circles – 55x30 points; dot-and-dash – 110x60; dashed line – 220x120; solid line – 440x240.

Fig.10 Profiles u(r), scaled with respect to local sound speed, as =0.420 for four grids: circles – 55x30 points; dot-and-dash – 110x60; dashed line – 220x120; solid line – 440x240.

Fig.11 Profiles w(r), scaled with respect to local sound speed, as =0.420 for four grids: circles – 55x30 points; dot-and-dash – 110x60; dashed line – 220x120; solid line – 440x240.

The purpose of the following numerical experiment is to demonstrate the abilities of the algorithm to calculate the flows with internal, including ‘hanged’, shock waves. For calculations the supersonic flow in the tube of variable section has been chosen. The geometry of the tube is shown in Fig. 12. (The scale is distorted for visualization purpose: real cross sectional sizes are 20 times smaller.) At the tube inlet one has supersonic flow with Mach number M = 1.1. Horizontal and inclined walls are conjugated in a smooth way. Central body is absent, hence the artificial needle is also absent.

Fig. 12 The domain geometry and the distribution of v component of the velocity

The calculations are performed for the sequence of grids, but presented results are referred to the grid 960x640. Such fine grid is necessary to get known theoretical result for axially symmetric non-swirled flows that the regular reflection of the shock waves from the axis is impossible – Mach reflection should arise, i.e. one has lambda-shock: the shock wave impinges the axis with respect to the normal direction to the axis. The same figure shows the distribution of radial component v of the velocity. This distribution rather well reflects the structure of shocks system. One can clearly see the system of oblique shock waves. At this first wave is ‘hanged’. It is also clearly seen that when passing the consecutive shock waves the v component changes sign.

Fig. 13 shows the distribution of pressure p at the tube wall. It is distinctly seen from this figure that the first shock wave is ‘hanged’ because the point of shock origination is located rather close to the wall and at the tube wall one sees the break (the discontinuity of the derivative) in the pressure graph, but not the shock.

Fig. 13 The pressure profile along the tube wall

However the secondary wave coming to the wall at reference point z ~ 0.127 can be clearly seen in Fig. 13 as the shock wave, smeared because of numerical approximation. The next shock at the wall is also seen though its smearing increases.

Next calculation is devoted to the case of swirled flow in the channel with subsonic outlet. The geometry of the channel is shown in Figs. 14, 15. The main difficulty in posing of the problem under the condition of subsonic outlet lies in the fact that in swirled flows the pressure sharply changes across the channel, and so it is unnatural to conventionally determine the pressure as constant at the channel outlet. The consequences of such determination would be crucial: at best the flow hardly distorts near the outlet and this disturbance influences all domain of the flow; but as a rule the inflow arises near the axis which requires other boundary conditions and hence other way of the problem determination. To leave out such difficulties we place here the outlet rather far from the axis in such a way that the swirling is moderate and boundary condition p = const can be satisfied. But even in this case the stationary flow is hardly achievable. In most cases the flow comes to the quasi-periodic regime with rather large amplitude of the shock wave movement along the channel. At this in the diffuser zone (in the domain of sharp channel expansion) large-scale intensively evolving vortex emerges. It emerges either in upper part of the diffuser or in lower part, but in a number of calculations from time to time goes from upper part to lower and inverse. We fail to observe any pattern in this process. Figs. 14, 15 demonstrate the stream lines for two moments of time in one of the calculations. These stream lines are in fact the integral curves for the field (u,v) that are obtained by the integration of the system of ordinary differential equations {dz/dt=u,dr/dt=v} with respect to fictive time t under various initial data. At every stream line the points, which correspond to equal interval of fictive time t, are drawn. This allows scoring the fluid speed along the stream line by the comparison of curves length between the points. The same pictures show the position of shock wave drawn by dash line. The notable displacement of the shock wave is observed. Other calculations reveal that this displacement can be rather big and even can lead to the destruction of the shock wave and to the temporal transformation to subsonic regime. The evolution of the shock wave is not directly connected with the sizes and position of the vortex – the flow has essentially non-stationary character. Such character of the flow is also confirmed by the non-closed stream lines in the vortex.

Fig. 14 Streamlines and shock wave. One moment of time

Fig. 15 Streamlines and shock wave. Another moment of time

To summarize, presented numerical method allows effectively performing the calculations for non-stationary as well as for stationary swirled flows in the channels with central body (and without it also). Such flows can have rather complicated structure: near axis jet, sharp pressure drops, supersonic swirling already in the premix chamber, the system of shocks in the diffuser, nontrivial and non-uniform evolution of shock wave position in the central part of the channel, etc. Typical ratio of the mesh size with respect to the radius and the mesh size along the channel attains 1/25, and such ratio does not influence the algorithm robustness. As the result of the calculations one can see the specifics of the problems formulation for the swirled flows in the channels of complicated geometry: the swirling at the inlet is determined by the vanes angle modeling; the nozzle ‘launching’ is gradual; the central body in the swirled flow generates detachment of the flow with the rise of stagnation zone; the outlet of the calculation zone should be far enough from the axis, etc. Also the features of the swirled flows in the channels with central body can be seen: the stationary regime, if exists, is not an attractor; large scale traveling vortexes arise; the shock wave can travel on the large distance, provided enough place for this, etc.

So, proposed calculation methodology works rather well for shown complicated multi-scale flow. It is desirable to increase the accuracy of the algorithm up to the second order provided that the algorithm robustness is conserved. Then this technique could be applied to other complicated Euler flows with the presence of various scales and to the problems, which involve Navier-Stokes equations.

References