1. Introduction
Numerical methods for gas dynamics problems play an important role in computational mathematics, and a vast literature is devoted to them, see, in particular, the monographs [
1,
2,
3,
4,
5] and references therein. Stability conditions for such methods, including those explicit in time, are of great practical and theoretical importance. For the explicit schemes, computational costs are inversely proportional to the time step
but the maximal value of
is restricted by a Courant-type stability condition involving steps of the spatial mesh and the sound speed.
Among the mentioned methods there exist a family of mesh methods based on preliminary viscous regularization of the equations involving so-called quasi-gas dynamics (QGD) and quasi-hydrodynamics (QHD) regularizations [
6,
7,
8,
9]. These methods are widely used for numerical solution of various applied problems. In the barotropic case, gas dynamics systems of equations with such regularizations were introduced and studied in [
10,
11,
12,
13], and their numerous applications to computer simulation were given for various 1D and 2D shallow water models [
14,
15,
16,
17,
18,
19], some 2D astrophysical problems [
20] and the 2D and 3D compressible Navier–Stokes–Cahn–Hilliard models [
21,
22,
23], etc. Despite of a lot of applications, almost nothing was known until recent years on rigorous theoretical conditions for stability of schemes with the QGD and QHD regularizations thus leading to additional time-consuming preliminary numerical experiments in order to choose the adequate parameters of schemes allowing to use larger values of
. Notice that there exist also other regularizations for gas dynamics equations, see in [
24,
25,
26], etc., however, they have not yet undergone such detailed practical testing.
The present paper is a contribution to the theoretical basis of schemes with the QGD and QHD regularizations. We study an explicit two-level in time and symmetric three-point in each spatial direction finite-difference scheme for the 2D and 3D barotropic gas dynamics system of equations with such regularizations linearized on a constant solution (with any velocity, though rather often previously much simpler case of the zero velocity was considered only). The possibility of applying symmetric approximation in space is related namely to application of the regularized equations. For this scheme, we give a matrix criterion as well as simpler and rather close to each other sufficient conditions and necessary conditions for
-dissipativity of solutions to the Cauchy problem, for any Mach number and the uniform rectangular mesh. To this end, we apply the spectral method [
4,
27,
28] and analyze matrix inequalities containing symbols of symmetric matrices of convective and regularizing terms. Notice that we apply the version of the spectral method based on the Fourier series [
4,
28] rather than the integral Fourier transform [
27]. Moreover, the analysis of namely the
-dissipativity is natural since the corresponding results are known for the linearized QGD and QHD systems of equations [
10,
11], though rather often only the von Neumann-type necessary stability conditions are studied that is simpler but does not ensure stability in any norm. Furthermore, the advantage of the spectral method versus an alternative energy approach is the possibility to derive not only sufficient conditions for stability but necessary conditions as well. In this paper, for the first time this analysis is generally performed in the unified manner for the both regularizations, but for the simpler QHD regularization it turns out to be more complicated since its regularizing terms are weaker with respect to the convective ones than in the QGD regularization.
Explicit formulas for the scheme regularizing parameter that guarantee taking the maximal time step are given for the derived conditions. We also present the choice of the “average” spatial step in the Courant-type conditions on the time step and in the regularizing parameter which ensures that the sufficient conditions and necessary conditions are uniform with respect to the rectangular spatial mesh and, for the QGD regularization, with respect to the Mach number as well (that can be valuable for simulation of super- and hypersonic flows). This choice depends not only on the steps of the rectangular spatial mesh, but also the Mach numbers in the respective directions. These results are valuable for practical applications helping one to choose adequately the scheme parameters allowing to use the maximal time step.
For the QHD regularization, these 2D and 3D results are derived for the first time and, for the QGD regularization, they improve those that have recently been obtained in [
29]. Previously, the
-dissipativity analysis of similar schemes in the much simpler 1D barotropic case was accomplished for zero and any Mach number in [
30,
31,
32], respectively. In this paper, we base on the papers in [
29,
32] and aim to develop further their technique. We also include the complete proof of the formula for the norm of the level-to-level transition operator in terms of the eigenvalues of its symmetrized symbol that is essential for the spectral technique; the main items of its proof have recently been presented in [
33]. Importantly, the technique is general enough and applicable for various schemes, other regularizations (for example, see in [
34]) and in more general or different statements of the equations, and some such studies are planned for the near future.
2. The Barotropic Gas Dynamics System of Equations with Two Regularizations, Its Linearization and the Corresponding Difference Scheme
We write the barotropic gas dynamics system of equations with the QGD and QHD regularizations in the form from [
12] in a unified manner, setting for them
and
, respectively. They consist of the mass and momentum balance equations which in the absence of external forces have the form
in
,
, for
. The sought functions are the density
and velocity
of a gas which depend on
and
t, also
is the pressure with
. The operators div and
are taken in
x, also
and
. The symbol ⊗ denote the tensor product of vectors. The divergence of a tensor is taken with respect to its first index.
The regularized mass flux
and the viscous stress tensor
are as follows:
where
, and hereafter the summation from 1 to
n over the repeated indices
(and only over them) is assumed. Furthermore,
is the Navier–Stokes tensor with
and the artificial viscosity coefficients
and
with the parameters
for
or
for
(the Schmidt number) and
,
is the unit tensor,
and
are the regularizing momentum and viscosity tensor, and
is the regularizing parameter.
For and given and , this system becomes the barotropic compressible Navier–Stokes system of equations. For , it is simplified as the barotropic Euler system of equations.
We recall the linearization of system (
1)–(
3) on a constant solution
and
[
10,
11]. Let
be the background sound speed and values of the parameters
,
and
. We substitute the solution in the form
and
into Equations (
1)–(
3), where
and
are the dimensionless small perturbations. We omit the terms of higher than the first order of smallness with respect to them and get the linearized system of equations for
. Its vector symmetrized form is as follows:
in
for
, where
and
are matrices of the order
(we omit the derivation details).
In order to write down these matrices, we introduce the column vectors
of the canonical basis in
, then
for any
k,
i and
j from 1 to
n and
. Hereafter
also
is the Mach number,
is the unit matrix of the order
k, the matrix
has the order
, where
is the diagonal matrix of the order
k with diagonal elements listed sequentially. Clearly
,
and
are symmetric matrices and
that is essential below.
For the solution to systems of equations like (
4) supplemented with the initial condition
, the following uniform in
bound is known [
10,
11]:
We define the uniform mesh
in
with the nodes
,
, and the step
,
, and the mesh
in
t with the nodes
,
, and the step
. We introduce the symmetric difference operators in
and the forward difference operator in
t
where
and
.
We define the rectangular mesh
in
with
. Let
H be the Hilbert space of vector-functions
:
defined and square summable on
, equipped with the inner product
(where, for example, ) and the corresponding norm .
We approximate the system of Equation (
4) using the defined difference operators and get the explicit in
t and symmetric three-point in each direction
difference scheme
on
, where
is the Kronecker symbol. A similar difference scheme arises after the linearization of schemes for the original Equations (
1)–(
3) given, in particular, in [
35].
We pose the question on conditions for the validity of the mesh counterpart of bound (
7), namely,
We define the level-to-level transition operator acting in
H:
where
I is the unit operator. Bound (
9) is equivalent to the properties
and the
H–dissipativity (i.e.,
–dissipativity) of the scheme
Notice that in the case of more general than (
8) non-homogeneous scheme
on
, with any given
and
, under the above mentioned property
, it is easy to derive the following more general stability bound:
3. An Analysis of the –Dissipativity of the Scheme
Let
and
be given by the formulas of the typical form [
6,
7,
8]
with the parameters
(like the Courant number) and
. Here, the “average” spatial step
is arbitrary during the whole analysis and can depend also on
and other parameters. Its adequate choice arises as a result of the analysis, and it will be given at the end of the paper. Below we derive conditions on
in dependence on
related to the validity of bound (
9).
According to the spectral method [
4,
28], we consider particular solutions of scheme (
8) in the form
where
is the imaginary unit and
is the vector parameter. We substitute them into (
8) taking into account the formulas
together with Formula (
10) and derive the explicit recurrent formula
In it,
is the matrix-symbol of the operator
having the form
where the matrices
and
are proportional to the symbols of convective and viscous (regularizing) terms as well as
and
Hereafter, mainly it seems more convenient to take as a parameter instead of ; obviously, .
We define the row vector and number
Lemma 1. The matrices and can be written in the -block symmetric form Proof. For
, the result was proven in ([
29], Lemma 1).
For
, by virtue of Formulas (
5) and (
6) for the matrices
together with
, we can write down
These formulas imply form (
11) with
of the matrix
. □
In ([
29], Lemma 1), the important matrix inequality was proven:
Recall that it follows from the formula
where
is the trace of
Q, which can be straightforwardly verified.
Now, we derive the corresponding more complicated inequality between the matrices and .
Lemma 2. 1. The following inequality holds: where the given constant is maximal in this inequality, with for . Herewith there holds the inequality 2. The following inequality holds: Proof. 1. The following formula holds:
It is not difficult to check that the inequality
is equivalent to the corresponding inequality for the
-blocks of the matrices
and
:
In the proof of ([
32], Theorem 3), it was shown that the maximal constant in this inequality equals
, where
is given in Item 1 of this Lemma. Herewith inequality (
17) implies that
and therefore knowingly
.
2. Let
. Due to the formula
, we can write down
where
and
.
The inequalities
and
mean that
for any
and
. The first of these inequalities is a consequence of the inequality
that follows from the estimates (already used in [
29])
The second inequality (
20) for
and due to the lower estimate
follows from the matrix inequality
It coincides with (
17) for
and, therefore, is valid for
. Consequently, we have
. On the other hand, for
, inequality (
20) implies the inequality with any
and
w that is equivalent to inequality (
21), thus the found constant
b is maximal in (
20).
Finally, we get
with
, since
runs over
when
runs over
S.
Herewith, in the inequality for any , the specified constant is maximal. This additional result is essential below. □
To apply Lemma 2, it is required to study behavior of the function .
Lemma 3. 1. For , the function increases in .
2. Let . Then decreases in and increases in and, moreover, , where 3. As a consequence of Items 1–2, the following formula holds: Herewith, and .
Proof. Let
. The function
is the larger root of the quadratic equation
Notice that due to Formula (
15), we have
excluding the particular case
. From (
14) and (
15), it also follows that
increases in
.
Except for the specified case, there exists
, and the differentiation of (
22) gives
As
, the property
is equivalent to
; moreover,
(otherwise
and
). As
, here
. Inserting this expression for
into (
22), after a series of simplifications we come to the quadratic equation
For
and
it has the real roots
Clearly and . It is not difficult to check that the property is equivalent to .
Therefore, for
, there exists a unique root
of the equation
for
, and as
increases in
, also
decreases in
and increases in
. Moreover, in this case, solving of the equation
for
after simple algebraic transformations leads to the quadratic equation
One of its roots is , therefore . In this case, we set and , and get formulas given in Item 2 of this Lemma.
In the opposite case , the derivative does not vanish for , thus increases for . Notice also that here . Thus, the results of Lemma are valid. □
Notice that in the case
(in particular, for
), we get
in the entire subsonic region
. On the other hand, in the case
(for example, for
and
), the quantity
is minimal in the entire subsonic region
and
.
Moreover, the following properties hold:
thus,
decreases rapidly as
M grows.
We denote by the maximal eigenvalue of the Hermitian matrix A. In the following theorem, we give the formula for the norm of the level-to-level transition operator in terms of the eigenvalues of its symmetrized symbol that is essential for the applied spectral technique. It will be more convenient to consider that . Let be the Hilbert space of vector functions : , equipped with the norm , where .
Theorem 1. The following chain of equalities hold: Proof. The complex Hilbert spaces
H and
are isomorphic by means of the multiple complex Fourier series: to any mesh function
corresponds the function
such that
, and vice versa, for the Fourier coefficients of this function, the following formula
is valid. Herewith, the Parseval equality holds:
Due to Formula (
24) and definition of
, we obtain
i.e.,
are the Fourier coefficients of the function
. Thus, due to the Parseval equality and the equality
we derive
Further, clearly
, and thus
here we take into account the continuity of
on
D.
It remains to prove the inequality of the opposite type
as the formula
with
is well known. To this end, we first write down the formula
for some
. Let
be an eigenvector corresponding to the eigenvalue
:
We construct the function
, where
is the characteristic function of the ball
with
. By definition of
, we can write
where
and
is the measure of
. Next, we transform and bound from below the arisen integral as follows:
Therefore, using the right Formula (
26), we derive the lower bound
Due to continuity of
, the matrix
is continuous on
D as well, thus we get
Therefore,
, and due to the left Formula (
26) we get inequality (
25). □
The main items of the presented full proof have recently been described briefly in [
33] (Theorem 1). Clearly, neither the specific form of
nor the specific dimension
of the involved vectors are essential in it (though the continuity of
on
D has been exploited).
Now, we present a criterion, sufficient conditions and necessary conditions for the validity of bound (
9). Recall that the quantity
has been introduced in (
12), and we set
to unify the form of inequalities (
12) and (
16).
Theorem 2. Let .
1. For the validity of bound (9), the matrix inequality serves as a criterion. Here, is the commutator of the matrices and , and the matrix is Hermitian.
2. For the validity of bound (9), the matrix inequality with any , is a sufficient condition.
Consequently, for , the same is valid concerning the inequality 3. For the validity of bound (9), the inequalities with are necessary conditions.
Proof. 1. Justifications of criterion (
27) and the sufficient condition (
28) (as well as (
29) for
) were given in ([
29], Theorem 2). Herewith, the derivation of the criterion is based on Theorem 1 and the equivalence of the properties
and
.
By virtue of (
12) and (
16) inequality (
28) follows from the inequality
The last inequality means the following inequality for the eigenvalues
of the matrix
:
As
, it is equivalent to
For
, choosing
we get the sufficient condition (
29).
2. For
, we find
and
thus
and
Now criterion (
27) leads to the necessary condition
, i.e.,
. We write down the lower bound by the maximal diagonal element:
by means of which the latter necessary condition implies (
30). Note that, for
, the matrix
is diagonal, and bound (
33) turns into the equality
.
To derive condition (
31), we set
with
and
similarly to the proof of ([
36], Theorem 4). Then, we obtain
with
,
and therefore
and
with
where the asymptotic relations for vectors and matrices are understood componentwise.
Herewith,
and
, therefore here criterion (
27), after division by
and passing to the limit as
, leads to the necessary condition
After division of the both sides by
and the formal replacement of
by
, it takes the form
with the symmetric matrix
Recall that the matrices
and
for
have been introduced in (
18) and (
19).
In the inequality
the constant
is maximal. For
, this has been established in the proof of Lemma 2, whereas for
this is a consequence of the relations
which follow from (
13). Therefore, condition (
35) implies the necessary condition (
31). □
Corollary 1. For the validity of bound (9), the following inequality is a necessary condition.
Notice that
and
as
or
. The maximum of the right-hand side of the sufficient condition (
29) is attained at
and equals
The maximum of the right-hand side of the necessary condition (
36) is attained at
and equals
For comparison, notice that for the validity of bound (
9), the matrix inequalities
serve as necessary conditions as well. For
, this was proved in ([
29], Theorem 2) by reducing to the 1D case and, for
, this is proved in the same way (recall that the 1D case was previously studied in [
32]). The form of matrices
,
, is inessential in this proof.
Due to Lemma 2, Item 1, for
, and due to ([
32], Theorem 1 and Remark 2), for
, these inequalities are equivalent to the number inequalities
where
and
. As
, for
, condition (
31) is more sharp than the second condition (
37), whereas for
they coincide.
It is not difficult to calculate eigenvalues of the matrix
, in particular, the maximal one, which is the maximal eigenvalue of its
block:
(recall that
first appeared in Lemma 2), namely,
see in [
29] and ([
32], Theorem 3). For them, the following two-sided bounds hold:
Furthermore, the following lower bounds hold:
therefore, the necessary condition (
30) is qualitatively stronger than the first condition (
37). It is not difficult to ensure this property on the quantitative level as well. To this end, it is required to strengthen bound (
33) by means of using the
blocks of matrix (
32):
for
, since then
This has not been implemented above in order not to complicate bound (
30) essentially.
To apply the sufficient condition (
29), we present the uniform in
upper bound for
.
Theorem 3. For and , for the eigenvalues of the matrix , see (11), the following bound holds:where , and also for or is arbitrary for . In the particular case , it holds that .
Proof. The following quadratic form corresponds to the matrix
:
with any
and
. As
, we have
and therefore the following bound holds:
As
and
, due to the Cauchy inequality we obtain
using the formula
with
introduced in the Lemma, similarly to the proof of ([
29], Theorem 4). Replacing
with
here, we also obtain
Therefore, using bound (
39) we get
with
, and due to the classical Rayleigh formula for
, we derive
It is not difficult to calculate that
Further, for
, we write down the estimates
On the other hand, for
, the right-hand side of inequality (
41) can be estimated as follows:
Consequently, due to the classical Rayleigh formula for , we get . □
Of course, on the right-hand side of bound (
38), we could strengthen
to the value
given in the proof; above this has not been done in order to avoid too cumbersome result.
In accordance with the derived bounds, the natural choice of
, depending on
and
only, is as follows:
Such a formula is suggested for the first time for schemes based on the QGD or QHD regularizations. For it, the following two-sided bound holds:
together with the equality
Due to this equality, we can estimate the constant on the right-hand side of bound (
38) as follows:
and, furthermore,
For and the latter estimate is obvious, whereas for it arises after taking the minimum in and simplifying the result a little.
Furthermore, for
in the necessary condition (
30), the following lower bounds hold:
Importantly, the given bounds lead to the sufficient condition and necessary condition independent of . Moreover, in the case , i.e., for the QGD regularization, they are also uniform in the Mach number that can be valuable for computing super- and hypersonic gas flows.
In addition, Formula (
10) can be rewritten in the form
For some other schemes, a similar formula for
but without powers 2 and
is contained in ([
1], Chapter 2).