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 nonstationary 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.
Àííîòàöèÿ
Â ïðåïðèíòå îïèñàíà ÷èñëåííàÿ ñõåìà íîâîãî êëàññà äëÿ ðåøåíèÿ íåëèíåéíûõ ñèñòåì óðàâíåíèé â
÷àñòíûõ ïðîèçâîäíûõ ãèïåðáîëè÷åñêîãî òèïà. Îñíîâíûì äîñòîèíñòâîì òàêèõ ñõåì ÿâëÿåòñÿ
ñïîñîáíîñòü ýôôåêòèâíî ïðîâîäèòü ðàñ÷åòû íåñòàöèîíàðíûõ çàäà÷, ñîäåðæàùèõ ðàçëè÷íûå ìàñøòàáû.
Êîíñòðóêöèÿ ñõåì òàêîãî òèïà îñíîâàíà íà õîðîøî èçâåñòíîì ïîäõîäå Ãîäóíîâà, ïðèíöèïàõ ÒÂÄ è
ñïåöèôè÷åñêîì ìåõàíèçìå ïåðåêëþ÷åíèÿ, êîòîðûé ïîçâîëÿåò ïðîâîäèòü âû÷èñëåíèÿ â ñîîòâåòñòâèè
ñ ÿâíîé èëè íåÿâíîé âåðñèåé ñõåìû â çàâèñèìîñòè îò ñòðóêòóðû ðàññ÷èòûâàåìîãî ïîòîêà. Â êà÷åñòâå
ïðèìåðà ïðèâîäèòñÿ ðàñ÷åò çàêðó÷åííîãî ïîòîêà ñ ïîìîùüþ ñõåìû ïåðâîãî ïîðÿäêà òî÷íîñòè.
Ýôôåêòèâíûé ðàñ÷åò òîíêîé ñòðóêòóðû òàêîãî ïîòîêà äåìîíñòðèðóåò âûñîêîå êà÷åñòâî îïèñàííîãî
ìåòîäà ïîñòðîåíèÿ ñõåì.
Introduction
In this paper we present a
numerical scheme of new type which is designed for fast calculation of
nonstationary boundary value problems for hyperbolic systems of PDE’s. We call
this scheme ‘explicitimplicit’ 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 nonregular 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 nonstationary
processes in swirled flows in the nozzles with nontrivial geometry. We use
Euler equations in cylindrical coordinates as a basic model. The swirling
coupled with the longitudinal flow and geometry effects generates complicated
nonstationary 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 ‘explicitimplicit’ 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 NavierStokes equations with this technique.
Section 1. Formulation of the
Problem
We consider nonstationary
axially symmetric flow of swirled nonviscous 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.
Nondivergent 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 nonlinear 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 nonpenetration 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 nondivergent form (3)
one has
_{} (15)
where
_{} (16)
Thus we consider the solution
of initialboundary 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 onedimensional
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 nonstationary 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 to_{}and 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 number_{}exceeds 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 onedimensional 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 threepoint
difference equation
_{} (24)
where _{} depend on _{} in the following way
_{} (25)
Let us note that there is an ambiguity for_{}when_{}. 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) _{}and_{}are simply substituted by_{}and_{}with the same indexes. The solution of threepoint 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 wellposed provided the
determination of boundary conditions for (17) corresponds to the signs of _{} at the boundaries.
The same idea is applicable
also for onedimensional 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
twodimensional 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 twodimensional
–_{},
_{}. The mesh function is_{}, the fluxes are _{}.
Suppose that_{}. Then the simplest twodimensional 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 onedimensional problem
because each stage is actually onedimensional.
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: nondivergent (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
onedimensional 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 zdirection, 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 crosssectional
shocks are absent, the calculation of corresponding fluxes _{} can be performed
using nondivergent 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 nonpenetration 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 3vectors_{}. 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 2vectors_{}. 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
2vector_{}) 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 nonstationary 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 illposed 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 axis_{}downstream.
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 quasistationary regime starts when thin axial jet of nonswirling
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_{0}.
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_{0} , 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 nonuniform 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; dotanddash – 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;
dotanddash – 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;
dotanddash – 110x60; dashed line – 220x120; solid line – 440x240.
Moreover, two points deserve
close attention. First, one can observe the presence of significant highspeed
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; dotanddash – 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;
dotanddash – 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;
dotanddash – 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 nonswirled flows that the regular reflection of the shock
waves from the axis is impossible – Mach reflection should arise, i.e. one has
lambdashock: 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 quasiperiodic 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) largescale 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 nonstationary character. Such character of the flow is also
confirmed by the nonclosed 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
nonstationary 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
nonuniform 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 multiscale 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 NavierStokes equations.
References
[1]

Radvogin Yu. B., Zaitsev N.A. “Multidimensional
minimal stencil supported second order accurate upwind schemes for solving
hyperbolic and Euler systems.” — Keldysh Institute of
Applied Mathematics Preprint, 1996, No. 22, 1 – 28.

[2]

Radvogin Yu. B., Zaitsev N.A. “A locally implicit
second order accurate difference schemes for solving 2D timedependent
hyperbolic systems and Euler equations.” — Proc. Of Intern. Conf.
on Spectral and High order Methods, Herzliya, 1998, in Appl. Num. Math., 33,
(2000), 525 – 532.

