Analytical models of unstable rarefied gas flow in channels and nozzles: simulation and comparison

Three different analytical models of unstable rarefied gas flows in channels and nozzles are compared using numerical simulations by MonteCarlo method. Numerical results have demonstrated the chaotic behavior of constructed nonlinear dynamic system and the limit properties of corresponding trajectories, attractors, and bifurcations of rarefied gas flows in channels. Flow conditions satisfying the experimental requirements are indicated where the instability of considered type can be detected. The advantages and the drawbacks of the considered approximations are detected and the features of obtained solutions are indicated. Recommendations are given for applying the results in practical applications and in numerical calculations of rarefied gas flows.


INTRODUCTION
The methods of nonlinear dynamics are applied ¦rst by Miroshin to study rar-e¦ed gas §ows in channels [13]. Considering the trajectory of each gas atom as a sequence of the re §ections from the channel walls (according to the assumption that gas §ow is free-molecular), the limit trajectories (i. e., attractors) by increasing number of collisions with the surface have been found analytically. Nonlinear dynamic system is studied corresponding to the multiple scattering of rare¦ed gas particles from the walls of a channel or of a nozzle.
The authors£ previous analytical and numerical investigations have shown that the limit behavior of the considered system is not chaotic everywhere [4 7], but there are regions of instability only in a subset of the phase space. In usual applications of the chaos theory, the most interesting cases arise when the chaotic behavior takes place on a ¤strange attractor.¥ Around this attractor, a large set of initial conditions will lead to orbits that converge to this chaotic 613 Article available at https://www.eucass-proceedings.eu or https://doi.org/10.1051/eucass/201911613 region. However, in a rare¦ed gas §ow, the most interesting thing is not a strange attractor which can be observed in rare¦ed gas §ows but a cascade of bifurcations [18]. The transition from regular to chaotic behavior of the system is here especially interesting, because such transition indicates the values of the parameters where the solution changes substantially with a very small deviation of these parameters.
Rare¦ed gas §ows are discrete dynamic systems (unlike continuous dynamic systems speci¦ed by di¨erential equations), because they are described by the laws of mutual interactions of gas particles and of collisions with the surface. Thus, these systems can exhibit strange attractors whatever their dimensionality, whereas in a continuous dynamical system according to the Poincar‚ eBendixson theorem, a strange attractor can only arise if the system has three or more dimensions. The cascade of bifurcations appearing for the ray model of the scattering function V was considered in [5]. The ray model, as well as the specular model, determines only one velocity of re §ected gas atoms by given incident velocity. However, the angles of incident and of re §ected gas atoms for the ray model could be di¨erent. The advantages of the ray model are various analytical expressions of momentum exchange coe©cients (depending on the angles of incidence and of re §ection) and the variety of the shapes of limit gas atoms trajectories. The di¨use addition (multiplied by the coe©cient σ) to the ray model gives us the more general ray-di¨use model. This addition causes a randomization and changes fundamentally the limit behavior of the studied dynamic system [6,7].
The main purpose of these new calculations is the analytical and numerical investigation of the properties of the nonlinear dynamic system corresponding to a rare¦ed gas §ow not only for di¨erent values of the parameters, but also for di¨erent analytical approximations of momentum exchange coe©cients. From the point of view of nonlinear dynamics, it corresponds to di¨erent iterative equations describing the trajectory of a gas atom. Basic assumption is high Knudsen number Kn > 1 (the number of mutual gas atoms interactions must be much less than the number of gassurface collisions). Applying local interaction approximations allows taking into account that Knudsen number is ¦nite, because the local interaction theory, being exact for free-molecular §ow, is con¦rmed by experimental approximations of momentum and energy exchange coe©cients for transition regime between free-molecular and continuum §ows. Meanwhile, the bifurcations, the attractors, and the corresponding physical values in the §ows have been detected. However, the problem of the empirical con¦rmation of the obtained numerically e¨ect is still di©cult, because the scattering conditions, as well as the regions of the parameters, are quite particular. Therefore, they are hardly reproducible experimentally. The considered bifurcations can essentially a¨ect di¨erent gas §ows applied in practice such as §ows in propulsion systems and in microelectronic vacuum devices.
The parameters of rare¦ed gas §ow in a channel are naturally connected with the parameters of the scattering function V (� u, � u ′ ). This function V depends on the velocities � u of incident upon the channel wall and � u ′ of re §ected from the surface gas atoms. The scattering function V determines entirely the parameters of a rare¦ed gas §ow in a channel (if the geometric shape is known) at high Knudsen number Kn (near free-molecular §ow, Kn → ∞). There is no need to eliminate possible mutual collisions of gas atoms (i. e., to assume free-molecular §ow), because these collisions are taken into account by applying local approximations [4]. However, the solution cannot be extended to the whole region of Kn because according to the locality hypothesis, most of the incident gas atoms must come from channel walls, so that the in §uence of mutual collisions of gas atoms is small [4].
In the previous analytical and numerical investigations [48], the instability of the §ow is detected in long enough channels or nozzles for some transition parameter values of V under some speci¦ed conditions. The basic condition is that the model describing the scattering function V is close to the ray model, i. e., to the generalized specular re §ection de¦ned by the expression [2]: Here, δ is the Dirac delta function and � u * is some speci¦ed velocity of re §ected from the surface gas atoms that is, in general, di¨erent from the specular velocity. The momentum exchange coe©cients for the ray scattering function V are expressed by where u and u ′ are the absolute values of the velocities � u and � u ′ ; and θ and θ ′ are the corresponding angles for incident and re §ected gas atoms [4].
The ray model has a better experimental con¦rmation (especially combined with di¨use scattering) in comparison with other surface interaction models widely applied in practical DSMC calculations, in particular with the specular-di¨use model [4]. The approximation of real gassurface characteristics by the ray model is good enough for di¨erent physical conditions. Therefore, it can be assumed that the ray model is valid also in the considered §ows.

ANALYTICAL MODELS OF DYNAMIC SYSTEM
Simulating successive gas atoms re §ections from channel walls, let us obtain the nonlinear dynamic system; the subscript m indicates the number of re §ections (Fig. 1). Denoting x m = tan θ m , x m+1 = tan θ m+1 (x m and x m+1 are the successive values de¦ning the trajectory of a gas particle in the channel) and expressing θ ′ from the presented formulae of momentum exchange coe©cients p and τ , one obtains the iterative equation: Here, the function ψ and the variable l m are determined by the geometrical shape of the channel, and the momentum exchange coe©cients p and τ are assigned by local approximations. If the solutions become unstable, the nonlinear dynamic system corresponding to Eq. (1) gets many di¨erent limit solutions ¡ attractors [18]. The corresponding parameters of the scattering function V represent the values of singularity. So, small regions of parameter values exist where the Feigenbaum period-doubling cascade [4,5,8] can be observed. Numerical calculations of the trajectories of gas atoms in a channel demonstrate signi¦cant changes of the aerodynamic characteristics of the §ow near these values of system parameters. The investigation of the behavior of the nonlinear dynamic system in these regions is very di©cult because the corresponding iterative Eq. (1) is much more complicated than the well-known equation for logistic map Therefore, additional restrictions are taken into account. First, let us consider (following [48]) two-dimensional §ow, i. e., §at or cylindrical channel (see Fig. 1); then, the function ψ in Eq. (1) is identical, and the connection between the angles θ ′ m and θ m+1 in successive mth and (m + 1)th points of collisions of a gas atom with the surface becomes simple: θ ′ m = θ m+1 . Second, only three of the following most applied local approximations are examined.
1. Depending on three-parameter approximation [47,9] based on the expansion of p and τ/ sin θ in terms of cos n θ, n = 1, 2, . . ., this model has been con¦rmed by experiments in many papers (cited in [4]) from St. Petersburg State University and Central Institute of Aviation Motors (Moscow) and can be described by equations: p = p 1 cos θ + p 2 cos 2 θ; τ = τ 0 sin θ cos θ .
Special cases of this approximation have been examined by Miroshin [13]: p 2 = τ 0 (¤model A¥) and p 2 = 2 (¤model Z¥). The coe©cients p 1 , p 2 , and τ 0 can be expressed in terms of aerodynamic values, such as Kn, Mach number, temperature factor, etc. Substitution into Eq. (1) permits to reduce three coe©cients p 1 , p 2 , and τ 0 to two parameters a and b transforming the iterative equation (1) to 2. Another approximation suggested also by Miroshin [1] (¤model B¥) reduces three coe©cients p 2 , p 4 , and τ 0 (p 2 and τ 0 remain the same and p 4 corresponds to higher power cos 4 θ) to two parameters a 1 and b 1 transforming the iterative equation (1) to 3. A more general approximation containing four coe©cients p 1 , p 2 , τ 0 , and τ 1 and not yet examined for nonlinear dynamic system parameters corresponds to the iterative equation with the same parameters a and b like in Eq. (2) and a new parameter d: For d = 0, Eq. (4) coincides with Eq. (2). Hence, the ¦rst approximation (2) is particular case of the third (4). However, the second model (3) cannot be represented in the form (4), but the main properties are similar to the corresponding values of the parameters. Figure 2 illustrates the behavior of the functions (2)(4) for the values of variables a, b, c, and d (and corresponding momentum exchange coe©cients) observed in real physical conditions of gas surface interaction (and in the experiments). In particular, a = 2 and b = 1.7 for the approximation (2) (dashed line), a = 0.2 and b = 0.9 for (3) (dotted line), and a = 2, b = 1.8, and d = 0.3 for (4) (solid line).
They have similar graphs functions (2)(4) and are all unimodal: they have single maximum and are monotone on both sides of this maximum. This property is necessary for the initiation of Feigenbaum cascade of bifurcations under the following additional condition: the Schwartz derivative must be negative [4].
Investigating the limit behavior of Figure 2 The functions in iterative equations (2)(4) corresponding to successive re- §ections of a gas atom from the walls in a channel nonlinear dynamic system by m → ∞ is the main problem related to studying the stability of rare¦ed gas §ow in a channel [17]. The necessary conditions for the existence of Feigenbaum cascade of bifurcation can be proven analytically for iterative equations (2)(4). All three functions are unimodal (for the mentioned values of the parameters as shown in Fig. 2). This comes from studying the sign of the ¦rst derivatives in a very large region of the variables a, b, and d. Schwartz derivative is negative: for all functions f (x m ) on the right in the iterative equations (2)(4). Miroshin has derived it analytically [1,4] for Eq. (3), i. e., for , on the base of the representation of Schwartz derivative which can be written after the transformation as A similar representation and the proof of the inequality Sf 1 (x) < 0 is obtained for the iterative equation (2), i. e., for , expressing Schwartz derivative as The third model (4), i. e., does not allow a simple formula for Schwartz derivative like (5); it is much more complicated. However, the inequality Sf 3 (x) < 0 can be proved analytically for Negative Schwartz derivative is necessary but not su©cient for the initiation of Feigenbaum cascade of bifurcations; a numerical con¦rmation is needed. Analytical and numerical investigation of limit solutions demonstrates that all the three iterative equations (2)(4) (2) and (3) can also have stationary points 2b), and x 0 = ∞, respectively. The last value has a physical meaning of locking the channel, i. e., its conductivity reduces signi¦cantly.

NUMERICAL CALCULATIONS
The main conclusion from numerical study of the third model (4) is that the attractors and bifurcations in this case are noticeably more various than for the ¦rst two models (2) and (3). In particular, up to three attractors exist corresponding to the roots of the equation The attractors of higher degree could be obtained from the equation f · · · f (x) = x (where f (x) could be iterated many times) and the set of the solutions could be very rich for some values of a, b, and d. A typical trajectory of the dynamic system described by the iterative equation (4) is presented in Fig. 3 in coordinates m and x in the region where it is unstable. The dynamic system is discrete; thus, connecting segments are shown only to make the trajectory visible.
Under certain conditions, iterative equations (2)(4) may have unstable solutions in some regions of the values of gassurface interaction parameters a, b, and d. In these regions (which can be found analytically in the described way), comparative small modi¦cation of the parameters a, b, and d causes signi¦cant di¨erence between corresponding limit values x m = tan θ m . From an aerodynamic point of view, it has the following interpretation: macroscopic parameters of the §ow may §uctuate unsteady while the di¨erence in microscopic values describing gas surface interaction remains negligible. However, the regions in which the §ow becomes unstable are very narrow; therefore, it is di©cult to ¦nd them numerically or experimentally. In the coordinate system (a, b), the regions of instability obtained analytically are concentrated near the line a = b. A real physical setup of the considered rare¦ed gas §ows has not yet a detailed experimental con¦rmation because of technical di©culties in performing such experiments.
However, these analytical evaluations are con¦rmed by numerical calculations.
The best way to study the unstable dynamic systems and their bifurcations is by constructing bifurcation diagrams.
The length l of the channel is assumed to be quite large (from 10 to 100) relative to its width in the Monte-Carlo simulations to obtain the discussed nonlinear dynamic e¨ects. Other parameters (local approximation parameters a, b, d, the part of di¨use scattered gas atoms σ, and initial angle θ 0 ) are changed in accordance with the known limit analytical solution of iterative equations (2) (4). In the initial section of the channel, a uniform §ow is set consisting of N gas atoms having identical velocities with the angle of inclination θ 0 (N changes from 10 4 to 5 · 10 4 ).
The ray-di¨use model of scattering of gas atoms is considered with the parameter a (on x-axis) and with the geometrical parameter x = tan(θ) of the trajectory of a gas atom in the channel on y-axis.
The bifurcation diagrams shown in Fig. 4 demonstrate the possible equilibrium (or long-term) values of a system as a function of the bifurcation parameter a. Therefore, nonrandom solutions are represented as lines (see Figs. 4b and 4c) or as black regions (see Figs. 4a4c). However, random solutions can  (2) be observed on bifurcation diagrams as distributed sets. At the ¦rst sight, the points on these diagrams seem to be distributed randomly but in fact, a sharp structure is well seen in all diagrams. The existence of this structure justi¦es the approach based on the nonlinear dynamics, and corresponding borders between the regions with di¨erent density of the points are the most interesting regions of the instability.
To illustrate the obtained results, comparative graphs of the numerical density of gas atoms in some sections of a channel depending on value a are presented (Fig. 5) for the ¦rst model (2).
The scattering function is assumed to be ray-di¨use with the identical value σ = 0.1. The parameters are changed gradually near the points of the instability to demonstrate the variation of the results. For instance, the variable a changes from 1.47 to 1.60 in four steps (see Fig. 5). The largest di¨erence between curves 1 and 2 corresponding to the nearest values a = 1.47 and 1.48 indicates the region of the instability in the segment 1.47 ≤ a ≤ 1.48.
Monte-Carlo simulation of di¨erent nonlinear dynamic systems corresponding to rare¦ed gas §ows in channels has indicated various developments of instability including cascades of period-doubling bifurcations. Several main conclusions follow from the analysis of the comparison: adding the extra parameter d allows one extending the region where the dynamic system (related to rare¦ed gas §ow) is unstable in the parametric space. Simulated unstable states of the system are close to physical situations observed in practice; numerical calculations con¦rm the analytical evaluations of the values of the parameters corresponding to unstable rare¦ed gas §ows. More sig-ni¦cant in applications is the possibility to use relatively simple analytical representations with the parameters found numerically. Hence, the present study is numerically useful; obtained connection between the parameters of the nonlinear dynamic system and momentum exchange coe©cients (or accommodation coe©cients) makes it possible to express analytic evaluations in terms of aerodynamic characteristics including Knudsen and Mach numbers, temperature factor, etc.; and to verify the bifurcations of simulated type experimentally, all considered physical values in the §ows are to be set exactly to the same values as detected in the present calculations. Therefore, this problem may be complicated enough.