Moscow, 2004
Ìîñêâà, 2004

Abstract

The paper describes the numerical scheme of new type for the solving of nonlinear hyperbolic
systems of PDEs. The main advantage of this scheme is the ability to efficiently compute
the non-stationary problems with different scales. The tools for its construction include
conventional Godunov approach, TVD principles and specific switching mechanism, which allows
performing the calculations with respect to explicit or implicit pattern depending on the
flow structure. The numerical computation of the swirled flow with the aid of first order
scheme is demonstrated as an example. The effective calculation of the subtle structures of
such flow demonstrates the relevance of presented schemes construction method.

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.

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. .

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.

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 p0.
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)/p0 , 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)

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.

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.