QUALITATIVE ANALYSIS OF THE GENERALIZED BURNETT EQUATIONS AND APPLICATIONS TO HALF–SPACE PROBLEMS

. The Generalized Burnett Equations, very recently introduced by Bobylev [3, 4], are tested versus Fluid–Dynamic applications, considering the classical steady evaporation/condensation problem. By means of the methods of the qualitative theory of dynamical systems, comparison is made to other kinetic and hydrodynamic models, and indications on an appropriate choice of the disposable parameters are obtained.


1.
Introduction. The relation between Boltzmann equation and hydrodynamics has been studied in the past by means of the Chapman-Enskog method. This method represents a systematic way of derivation of formal equations of hydrodynamics with any given order of accuracy with respect to the Knudsen number ε (mean free path divided by macroscopic length), and it allows to obtain in cascade closed approximate systems of partial differential equations like the Euler, Navier-Stokes, Burnett, and higher order equations. This way of presenting the Chapman-Enskog procedure [10] was questioned by Cercignani [6] and other authors that were sceptic about higher equations of hydrodynamics, and in particular, for the uncertainty of boundary conditions. Furthermore, in 1982 it was proven [1] that Burnett equations for Maxwell molecules are ill-posed, because of their solution, in particular the constant equilibrium solution, that is unstable with respect to shortwave perturbations. In the last two decades, several authors tried to overcome this problem by combining the Chapman-Enskog expansion with moment methods, and by using some higher order terms for regularization of Burnett equations (see for instance [2] for a complete bibliography). It is worth to cite in particular the attempts to stabilize Burnett equations by relaxation due to Jin and Slemrod [9]. An alternative method recently adopted [2,3,4] does not use any information beyond the Burnett level, but it concentrates on the problem of truncation. In particular, this method, based on simple ideas from the theory of dynamical systems, depends on a suitable choice of coordinates which can influence the result of truncation at Burnett level. A first attempt to regularize Burnett equations has led to Hyperbolic Burnett Equations (HBE) [2] which turned out to be, however, more complicated than the original Burnett equations and did not have the standard form of equations of hydrodynamics. Another choice of coordinates has been recently devised in [3,4], and a simpler regularized version of Burnett equations, called Generalized Burnett Equations (GBE), has been derived in detail. Here, we start from these results and we carry on the analysis by focusing on the stationary version of such equations in one-dimensional space for Maxwell molecules. The governing equations to be investigated are then a set of ordinary differential equations, with the space variable as independent variable. In particular, we concentrate on solutions for gas flows driven by evaporation and condensation at interphase surfaces, and we investigate these solutions in comparison with results from kinetic theory and from other hydrodynamic models.
The paper is organized as follows. In section 2 the Generalized Burnett Equations in one-dimensional space for Maxwell molecules are recalled and commented on, and some of their features are presented for the readers' convenience. Section 3 deals with the stationary version of GBE, considered as a finite dimensional dynamical system, going through equilibria and their linear stability, in order to apply it, as a test, to the standard problem of weak evaporation/condensation. Comparison to kinetic results allows a proper selection of the free parameters which characterize this hydrodynamic model. Quantitative analysis of numerical results is left as future work. However, the nonlinear problem and its normal forms are resumed in section 4, and phase trajectories are deduced numerically in order to show the qualitative structure of the whole phase diagram in an exhaustive way. Finally, section 5 contains some concluding remarks. Further application of the present model to other classical fluid-dynamic problems, like e.g. shock wave structure, is in progress.
2. Generalized Burnett equations. The Generalized Burnett Equations (GBE) are a closed set of partial differential equations for auxiliary unknowns ρ, u, T allowing to determine the actual fundamental hydrodynamic fields ρ tr (number density), u tr (mass velocity), and T tr (temperature) of a rarefied gas. In one space dimension and for Maxwellian molecules (we shall stick to this model in the present work), actual and auxiliary fields are related by where ε is the Knudsen number, with a difference affecting only temperatures and determined solely by the function The functions a(T ) and b(T ) contain two free scalar parameters θ 1 and θ 2 which are typical of the model, and read as where the coefficient η appears in the constitutive equation µ = ηT for viscosity and its expression comes from kinetic theory [1]. The GBE for auxiliary ρ, u, T take the form of standard hydrodynamic equations where D 0 = ∂ t + u ∂ x and the viscous stress Π and the heat flux Q are made up by the usual O(ε) Navier-Stokes and O(ε 2 ) Burnett contributions [6,7], with the addition of further O(ε 2 ) contributions, ε 2 π and ε 2 q respectively, representing the sought generalization, aimed at improving the reliability of the truncation [4]. As a result, the actual approximate constitutive equations closing the exact conservation laws turn out to be The main properties of this hydrodynamic closure, in particular hyperbolicity and stability, can be investigated in detail, and regions of stability in the (θ 1 , θ 2 ) plane can be determined. For our analysis it will be sufficient to consider the most significant one, given by the closed triangle We refer once more the interested reader to [4] for such an interesting derivation and discussion. Here it suffices to comment briefly on some points that show some aspects of this new model, and that will be resumed in the following sections. By proceeding to a standard linearization and adimensionalization of (4) around a uniform equilibrium, we obtain (setting ε = 1, with the mean free path as length scale) the following equations in the small dimensionless perturbations ( ρ, u, T ) By resorting to the usual plane-wave ansatz Φ =Φe i(ωt−kx) , where Φ = ( ρ, u, T ), Φ represents the complex amplitude of the wave, ω its frequency and k the wave number, we are led to a linear homogeneous algebraic system with matrix The dispersion relation det L = 0 relating k and ω reads as (9) In this way of looking at stability, a real k is fixed as wave number, and the corresponding complex frequency ω = ω r (k) + iω i (k) is deduced, yielding the damping α = ω i (k) (which must be nonnegative for stability) and the propagation speed ω r (k)/k. The other way around, a real ω could be kept fixed in order to deduce the corresponding k = k r (ω) + ik i (ω). According to [12], this constitutes an additional criterion for stability, the so called "stability in space". Taking into account travelling directions, this would lead to the requirement that k r and k i have opposite sign. Though this second criterion might be questionable, we test it too in order to get more information on the performance of this hydrodynamic model.
For illustrative purposes, we report some figures, labelled from 1 to 4, that are obtained for different values of the parameters θ 1,2 in the region (6), and show on the left ω(k) in the complex plane with k as a positive parameter, and on the right k(ω) in the complex plane with ω as a positive parameter.      Of course, we observe in all cases Im(ω(k)) ≥ 0 for any k, in agreement with stability, guaranteed by (6) [4]. However, the plots exhibit some differences in the selected points of the parameter space, and indeed, the most appropriate choice of the disposable parameters is clearly one of the first important issues to be solved for this proposed hydrodynamic model. Some steps in this direction are in fact moved in the sequel. Differences are even more striking when plotting k versus ω. There is for instance a degeneracy of the six branches to only four in the cases of Figures 1  and 2, and moreover, in the mid point of the vertical side of triangle (6) (Figure 4), there are even branches in the first and third quadrant, violating thus the so called "stability in space" criterion.
It might be worth making a comparison with other equations of hydrodynamics, in particular with the most popular Navier-Stokes Equations (NSE), achieved by deleting the third derivatives in (7). The corresponding dispersion relation is and the relevant plots are shown in Figure 5.  For both (9) and (10) the three branches of ω(k) for k → 0 are The purely imaginary root ω im (k) tends to finite asymptotic limit as k → +∞ for NSE (10), namely while for GBE the limit can be either finite or infinite depending on the value of the parameters θ 1,2 . In particular, for θ 1 = θ 2 = 0 or θ 1 = 0, θ 2 = 1, an easy calculation yields Otherwise, in the case θ = 1, θ 2 = 0 it can be verified that ω im GB is infinite as k → ∞, and diverges as iα k 2 + O(k), with computable α. The other two branches are complex conjugate, and their real and imaginary parts both diverge as O(k 2 ) for k → ∞, again with computable coefficients (details are omitted). This is different from NSE, where such branches coalesce for k sufficiently high, and from that threshold on the two roots become purely imaginary and eventually diverge when k → ∞. This typical loop diagram is common in literature, and it occurs also for Burnett and super Burnett equations (which give rise to instability), as well as for the so called "augmented Burnett" (which are instead stable) [12].
3. Stationary GBE and weak evaporation/condensation. We shall consider here the stationary version of the previous GBE, having in mind a classical fundamental steady problem of fluid dynamics like evaporation/condensation in a half space. Measuring distances in mean free paths (which amounts to taking ε = 1), the governing equations reduce to the set of ordinary differential equations which may be cast in the frame of the theory of dynamical systems, though the independent variable x is not necessarily confined to R + . Because of expressions (5), these are third order differential equations; however, three first integrals can be easily deduced from them, leaving derivatives only up to second order. More precisely, equation (11a) immediately yields ρu = constant = J 1 . Then, equation (11b) may be rearranged as d dx (ρ(u 2 + T ) + Π) = 0, which provides the conservation J 1 u + ρT + Π = J 2 . Finally, equation (11c) may be rewritten in conservative form as d dx J 1 (u 2 + 5T ) + 2(Πu + Q) = 0, and by substituting the second conservation law multiplied by u we obtain −J 1 u 2 + 2J 2 u + 3J 1 T + 2Q = J 3 , involving second order derivative of u through Q. Thus, system (11) may be replaced by containing the three integration constants J k , k = 1, 2, 3. A first elementary step in any dynamical system is determination of fixed points (i.e. x-independent solutions in this case). It is easy to see that there are for (11) just two such points of physical meaning for fixed values of the J k with J 1 = 0, given by whereas the degenerate case J 1 = 0, of no interest here, would yield infinite points with u * = 0. Since in this paper we shall deal with problems having a stationary state (say A) at infinity, it proves convenient to make reference to its coordinates as scaling parameters, and to call them ρ ∞ , u ∞ , T ∞ with ρ ∞ , T ∞ > 0 and u ∞ = 0. This determines integration constants as The other equilibrium point, say B, has then coordinates and is physical only for u 2 ∞ > 1 3 T ∞ . Now, adimensionalization consists in introducing the scaling which spontaneously singles out a dimensionless parameter where M A is the Mach number in the reference state. Notice that the Mach number in the other fixed point B is and that the two points A and B coincide for S = 5 3 (M A = M B = 1). Otherwise, for S = 5 3 and S > 1 3 , one of the states is supersonic and the other subsonic. The situation in this respect remains then essentially unchanged with respect to the analysis performed in [5] for the simpler Navier-Stokes equations. Omitting all tildes, and denoting derivatives by a prime, the scaled GBE read as and their stationary points are related by with α = A, B.
It is important to stress that (12a) implies necessarily a constant sign for mass velocity u, so that, once u ∞ = 0 is given, the dimensionless velocity u in necessarily positive, which explains the definiteness in sign of (19a) and of the dimensionless equilibrium velocities in (20).
We shall apply the present formalism to the analysis of the half-space evaporation/condensation problem, in which the gas is filling the interval 0 < x < +∞, with proper boundary conditions at x = 0, where it is in contact with its condensed phase, and with the reference state A as asymptotic condition for x → +∞. We have evaporation for u ∞ > 0, or else condensation for u ∞ < 0. As observed in [5], the two problems may be treated in a unified way in dimensionless form by reversing the x-axis in the case of condensation. In other words, the scaled equations (19) should be considered all over the real line −∞ < x < +∞ and with asymptotic state (1, 1, 1) for both problems. The difference is that the domain of interest is R + and the asymptotic state is relevant to the limit x → +∞ for evaporation, whereas the domain of interest is R − and the asymptotic state is relevant to the limit x → −∞ for condensation. Before proceeding in the analysis, we remark that the first of (19) is algebraic and allows to get rid of the unknown ρ as ρ = 1 u , which leads to the reduced system of two second order equations in two unknowns and thus to a four dimensional dynamical system, with variables u, T , Y = u ′ , Z = T ′ . With reference to the classical weak evaporation/condensation problem at the kinetic level, we focus on the analysis described by Sone in [11], where a semiinfinite expanse of a gas bounded by its plane condensed phase was considered at a kinetic level and it was analyzed from a numerical point of view. The problem was faced in terms of BGK equations for the velocity distribution function, while pressure, temperature and flow velocity were recovered as proper moments of it. Stationary solutions were looked for and it was found that no time-independent solutions exist for supersonic evaporation, while in the subsonic case there exist two relationships among the three free dimensionless parameters characterizing the problem, and relating quantities at the wall and at infinity. In particular, both pressure and temperature ratios turn out to depend on the flow velocity at infinity and the problem has one free parameter and a one-dimensional manifold of solutions. For condensation processes, in the subsonic case there is one relationship binding together the parameters, so that the pressure ratio can be considered as a function of temperature ratio and flow velocity and the problem has two free parameters and a two-dimensional manifold of solutions. Finally, for supersonic condensation all parameters are free in order to have a stationary solution, there are no constraints on them and then the problem admits a three-dimensional manifold of solutions. Thus, the dimension of the solution manifold grows of one order of infinity until 3 for supersonic condensation.
Here we want to analyze the weak evaporation/condensation problem for GBE and to determine the dimension of the solution manifolds in both subsonic and supersonic case, in order to compare the results with the kinetic level. The linearized version of equations (22a)-(22b) around the equilibrium A = (1, 1, 0, 0) in the phase space can be obtained by putting into equations (22a)-(22b) and neglecting higher order terms in the perturbations w, τ, y, z; this yields 37 18 The usual exponential ansatz e λx leads to the characteristic equation for the exponent λ Approach and expressions resemble of course those of the case studied in section 2 (dispersion relation) with ω = 0 (time independent problem). However the physical context is different from that, since we are dealing here with a half-space problem, where boundary conditions have to be given at the wall x = 0, and with a well defined constraint to be fulfilled in the asymptotic limit when x tends to infinity. The effects of such constraints on form and structure of possible solutions is just the additional information we are seeking for now.
We notice that when S = 5 3 the characteristic equation (25) admits a null solution and correspondingly a so called "transcritical bifurcation" ( [8], p.159) occurs between equilibrium states A and B, namely equilibria cross each other and exchange stability. When the matrix can be put in diagonal form, the general solution of system (23) can be written as  where C i , i = 1, 2, 3, 4, are arbitrary constants, λ i the eigenvalues of system (23) andê i the corresponding eigenvectors. Analogous forms are in order when the matrix in (24) is not diagonalizable, involving modes of the kind x e λix and so on.
It is important to notice that, for the evaporation problem, the asymptotic state A is reached as x → +∞, and the general solution (26) diverges in correspondence of the eigenvalues with positive real part. Vice versa, for the condensation process, the divergent exponentials in (26) are those corresponding to eigenvalues with negative real part. In both cases, the only way to reach the state A is to set C i = 0 for the constants relevant to the divergent modes.
A necessary condition to reproduce the kinetic scheme [11] for evaporation/condensation processes is that (25) reduces to a third degree polynomial. In fact, in the case of four eigenvalues it is impossible to represent in a right way evaporation and condensation processes simultaneously, and analogously in the case of second degree polynomial, because the stable manifold relevant to equilibrium A may have at most dimension two. Thus, we must require that the coefficient of the fourth degree term in (25) vanishes and also that the coefficient of λ 3 never becomes zero. Indeed, any eigenvalue would bring in, regardless of the dimension of its eigenspace, as many fundamental modes as its multiplicity into the solution of the linearized equations, and these modes should be either accepted or rejected according to whether the limit x → +∞ or x → −∞ is considered, i.e., according to whether evaporation or condensation is dealt with. Since kinetic results previously described indicate that the sum of the dimensions of the solution manifolds relevant to evaporation and condensation must be equal to 3 in both supersonic and subsonic conditions (precisely, 3+0 or 2+1), the same results can be in order for the present model only if the overall number of fundamental modes is equal to 3, leading to above necessity of a degeneracy in the characteristic equation (25). Such a constraint leads to θ 1 = 0 and either θ 2 = 0, or θ 1 + θ 2 = 1, covering two boundary edges of the considered triangular stability region (6). Analogous considerations on Hyperbolic Burnett Equations (HBE), another regularized version of Burnett equations analyzed in detail in [2], and on Navier-Stokes equations, show that both sets of equations cannot reproduce the kinetic pattern. In fact, HBE do not depend on free parameters and their linearized version leads to a fifth degree characteristic polynomial which can never be reduced to a third degree polynomial as for GBE (details are omitted for brevity). On the other hand, by focusing on Navier-Stokes equations, it has been shown in [5] that the linearized system has always only two eigenvalues, which do not allow to describe in a correct way the condensation problem.
In order to see whether the GBE characteristic equation (25) with constraints (27) actually reproduces the kinetic structure, we rewrite it as and notice that, for both options in (27), the phase space collapses to three dimensions, (w, τ, y), since z may be eliminated from (23) and then computed "a posteriori". We realize first that, for S > 5 3 , the sign of coefficients in (28) is such that any real solution is necessarily positive. On the other hand, one root is real "a fortiori", and, if the other two were a complex conjugate pair, their real part would again be positive. In fact, direct substitution and elimination of the imaginary part leads to a cubic equation for the real part that, in force of the inequality has exactly the same pattern as (28). As already observed, one real eigenvalue becomes zero for S = 5 3 , and indeed it is clear that a negative root must exist for (28) when S < 5 3 . No other transition of real roots from positive to negative is possible for S ≤ 5 3 , since S = 5 3 is the only bifurcation value, and in such a case λ = 0 is a simple root. On the other hand, possible crossing of the imaginary axis by a conjugate pair changing real part from positive to negative is also excluded. In fact, purely imaginary roots of (28) cannot exist: assuming the opposite would lead, as it is easy to check, to a compatibility condition that is always violated in force of (29). In conclusion, for S < 5 3 , one eigenvalue is negative and the other two are either positive, or complex conjugate with positive real part. Conjugate pairs do actually occur in some ranges of the two parameters S and θ 1 . For instance, for S = 5 3 , one root is zero, the other two have positive real part, and are real for θ 1 ≤ 9409 15540 , complex otherwise. Double or even triple roots may occur, implying coalescence or splitting of eigenvalues. For instance, at S = 5 53 28 + √ 1261 =:S, a triple root λ = 3 97 (23S − 15) occurs when θ 1 =θ 1 = 18818 1665S (23S−15) ≃ 0.5515 . Some typical trends are shown in Figure 6, where solid line denotes a real eigenvalue and dashed line the real part of a complex conjugate pair. We can observe the change in the coalescence process versus S when θ 1 is increased beyond the critical valueθ 1 . In any case, the main conclusion of the present analysis is the following Proposition 1. On the two segments of the parameter space (θ 1 , θ 2 ) defined by in which the phase space of the dynamical system (23) degenerates to three-dimensional, the equilibrium point A defined by u = 1, T = 1 is unstable, and, more precisely: i) its stable manifold is empty and its unstable manifold is three-dimensional, for S > 5 3 ; ii) its stable manifold is one-dimensional and its unstable manifold is two-dimensional, when S < 5 3 . We recall that the asymptotic state A is supersonic when M A > 1 (S > 5 3 ) and subsonic when M A < 1 (S < 5 3 ). It can be noticed that analogous definitions and considerations would hold if the other equilibrium point corresponding to the same conservation laws, namely point B, were chosen as asymptotic state instead of A. The correspondence is trivially determined by the fact that The evaporation/condensation problem in both subsonic and supersonic case gives rise to four different options and the corresponding solutions when the asymptotic state is represented by the equilibrium A are the following (eigenvalues shall be ordered by decreasing real part).
-Supersonic evaporation (S > 5/3, A asymptotic state as x → +∞) All the eigenvalues have positive real part, then we must require C 1 = C 2 = C 3 = 0 in the reduced three-dimensional version of (26) and there exists only the trivial solution (u, T, Y ) = (1, 1, 0). No other non-trivial solutions for outgoing supersonic flow are possible. -Subsonic evaporation (S < 5/3, A asymptotic state as x → +∞) The eigenvalues λ 1 , λ 2 have positive real part, then we must impose that C 1 = C 2 = 0 and the solutions have the following asymptotic form as -Subsonic condensation (S < 5/3, A asymptotic state as x → −∞) The real eigenvalue λ 3 is negative and we must require C 3 = 0. Thus, the solutions have the following asymptotic form as -Supersonic condensation (S > 5/3, A asymptotic state as x → −∞) All the eigenvalues have positive real part, then the complete general solution is admissible as x → −∞, because it is bounded at infinity for all values of the constants C i , i = 1, 2, 3.
In Table 1 we compare the dimension of the solution manifold relevant to the asymptotic equilibrium A for the stationary weak evaporation/condensation problem in both Burnett (Generalized and Hyperbolic equations) and Navier-Stokes hydrodynamic models of the kinetic equation. We may thus conclude that GBE with this choice of parameters describe adequately, from a qualitative point of view, the nature of steady weak evaporation/ condensation problems, reproducing the structures indicated by kinetic theory, and improving crucially the response with respect to other hydrodynamic-type models, like NSE and HBE. 4. GBE as dynamical system and phase diagram. We go back now to the nonlinear stationary problem. The phase space for the dynamical system generated by GBE (22) is in general four dimensional (it was two dimensional for the same problem tackled by the Navier-Stokes equations), but there are values of the free parameters θ 1,2 which can reduce such dimensionality as already pointed out. It would be desirable to put, at least locally, these governing equations in the normal form x ′ = F(x), where x = (u, T, Y, Z) involves the new unknowns Y = u ′ , Z = T ′ . It is immediately realized however that this task is not trivial at all and often leads to quite severe constraints.
We are interested in particular to GBE corresponding to parameters θ 1,2 on the two edges given by the conditions (30). For θ 2 = 0 (and θ 1 = 0), equations (22) states. Moreover, there exists a two-dimensional unstable manifold W u loc (A) tangent at A to the eigenspace relevant to the two eigenvalues with positive real part. All the points of this surface give rise to trajectories converging to A as x → −∞ and the surface represents the condensation solutions. In Figure 7 we report the phase portrait for S = 1 < 5 3 , with the evaporation curve which represents the one-dimensional stable manifold of A, and one of its two branches is the shock wave curve. Some condensation curves lying on the two-dimensional unstable manifold of A are also shown in the plot. Such a surface represents a kind of separatrix in the phase space, similarly to the Navier-Stokes hydrodynamic model [5], whose two-dimensional phase diagram is reported in Figure 8.
We notice that the case S > 5 3 can be analyzed in the same way, by exchanging the role of fixed points A and B. In particular, A is now supersonic, its unstable manifold W u loc (A) is three-dimensional, and one of its curves is the shock solution leading for x → + ∞ to the other fixed point B which is now subsonic. In addition, in agreement with previous results, the stable manifold W s loc (A) is now empty.
5. Concluding remarks. The question whether GBE are a good model for practical application is a challenging problem. Of course, GBE are characterized by the pro's and con's of any hydrodynamic model, but their stability is certainly a big step forward, and allows a significant improvement with respect to a Navier-Stokes description, which is necessarily poorer by its own nature. The preliminary analysis developed here shows at least that for the specific evaporation/condensation problem there is an optimal choice of parameters leading to a solution structure that fits very well the one from a much more difficult kinetic approach. The question anyway deserves further investigation of other classical problems, like steady shock waves, which is presently in progress.