## Abstract

Rarefied gas flows are highly nonequilibrium flows whose flow physics cannot be discerned accurately within the framework of the Navier–Stokes equations. The Burnett equations and the Grad moment equations, which form a super-set of the Navier–Stokes equations, have been proposed in the literature to model such flows but not much success has been achieved because of some inherent limitations of these equations. In this review article, we mainly focus on the recently proposed Onsager-Burnett equations (Singh et al., 2017, “Derivation of stable Burnett equations for rarefied gas flows,” Phys. Rev. E **96**, p. 013106) for rarefied gas flows, and the progress achieved so far by solving these equations for some benchmark flow problems. Like Burnett and Grad equations, the OBurnett equations form a super-set of the Navier–Stokes equations and belong to the class of higher order continuum transport equations. However, there are two fundamental aspects where the significance of the OBurnett equations is clearly visible. First, the OBurnett equations are unconditionally stable as well as thermodynamically consistent unlike the conventional Burnett and Grad moment equations. Second, the OBurnett constitutive relations for the stress tensor and the heat flux vector do not have any higher order derivatives of velocity, pressure, or temperature. This is quite significant since now the equations need the same number of boundary conditions as that of the Navier–Stokes equations. As such, the OBurnett equations form a complete theory, which cannot be said for the conventional Burnett equations. These two important aspects help to set the OBurnett equations apart from the rest of the higher order continuum theories. The results of the OBurnett equations are compiled for two benchmark rarefied flow problems: force-driven compressible Poiseuille flow and the normal shock wave flow problem. For force-driven compressible Poiseuille flow, the OBurnett equations successfully capture the nonequilibrium effects such as nonuniform pressure profile and presence of normal stresses and tangential heat flux in the flow. The accurate description of highly nonequilibrium internal structure of normal shocks has always been the stringent test for the higher order continuum theories. The results of the OBurnett equations for normal shocks show that there is no theoretical upper Mach number limit for the equations. Further, the equations predict smooth shock structures at all Mach numbers, existence of heteroclinic trajectory, positive entropy generation throughout the shock, and significant improvement over the results of the Navier–Stokes equations. Finally, the recently proposed Grad's second problem, which has the potential to become a benchmark problem, is discussed. The solution of Grad's second problem for different interaction potentials (Maxwell and hard-sphere molecules) within the Burnett hydrodynamics is also presented at length and some important remarks are made in this context.

## 1 Introduction

In the fluid dynamics community, the Navier–Stokes equations (referred to herein as the mass, momentum, and energy equations taken together) based on linear constitutive laws (Newton's law of viscosity and Fourier's law of heat conduction) are considered sacrosanct. However, there are several flow phenomena occurring regularly in nature such as shock waves, flow around space reentry vehicles, air flow in soil, lungs, and other small passages, etc. where these celebrated equations fail to accurately model these flows. Having accurate models is important for precisely calculating the drag force on aerosols, wave propagation in the Martian atmosphere, heat and mass transfer with gas flow in a microchannel, thrust from a microthruster, separation of isotopes in a gas centrifuge, and many other engineering problems. Such flows exhibit strong nonequilibrium behavior and are generally characterized by two dimensionless numbers, namely, Knudsen number (Kn) and Mach number (Ma) while the third nondimensional number, Reynolds number (Re) can be expressed in terms of other two as $Re=\pi \gamma 2MaKn$ with *γ* being the ratio of specific heats. The Mach number, as we know, is a ratio of the speed of the perturbation to the speed of sound. The Knudsen number, typically defined as the ratio of mean free path (*λ*) and characteristic length scale (*L*), measures the degree of rarefaction. This parameter becomes large in microscale flows (small characteristic length scale) or flows at higher altitudes (large mean free path) and generally marks the breakdown of the local equilibrium hypothesis and small gradient approximation that are inherent in the Navier–Stokes equations.

Based on an order of magnitude of Knudsen number, the flow can be demarcated into four different regimes [1–4] as shown in Table 1. For vanishing Knudsen numbers ($Kn\u21920$), the flow is close to equilibrium and can be modeled using the conventional Navier–Stokes equations with regular no-slip boundary condition. However, as the Knudsen number increases (Kn* *>* *0.001), there is a departure from the equilibrium condition and nonequilibrium effects such as velocity slip and temperature jump at the wall boundaries are observed. These nonequilibriated flows can still be modeled within the continuum approach using the Navier–Stokes equations but now the wall boundary conditions need to be accounted for through the velocity slip and temperature jump conditions [5,6]. Although it is not possible to capture all the nonequilibrium effects in entirety through the modification of the boundary conditions alone [2,7–9], the Navier–Stokes equations with slip and jump boundary conditions do provide a reasonable approximation for weak nonequilibrium flows and help to expand the Knudsen number envelope to 0.01. As we go in to the transition regime ($0.1<Kn<10$), we have a strong departure from the equilibrium condition, and this is perhaps the most difficult flow regime to model. The flows with Knudsen number greater than 10 are generally referred to as free molecular flows where intermolecular collisions are less significant as compared to molecule-wall collisions [10]. Such flows represent the highest degree of rarefaction and can be modeled using the collision-less Boltzmann equation.

Knudsen number | Regime | Fluid model |
---|---|---|

Kn < 0.001 | Continuum | Navier–Stokes equations |

Wall boundary conditions: no slip | ||

$0.001<Kn<0.1$ | Slip flow | Navier–Stokes equations |

Wall boundary conditions: Velocity slip and temperature jump conditions | ||

$0.1<Kn<10$ | Transition | Higher order continuum transport equations |

(e.g., Burnett, Grad, OBurnett, O13), DSMC | ||

Kn > 10 | Free molecular flow | Collisionless Boltzmann equation, DSMC |

Knudsen number | Regime | Fluid model |
---|---|---|

Kn < 0.001 | Continuum | Navier–Stokes equations |

Wall boundary conditions: no slip | ||

$0.001<Kn<0.1$ | Slip flow | Navier–Stokes equations |

Wall boundary conditions: Velocity slip and temperature jump conditions | ||

$0.1<Kn<10$ | Transition | Higher order continuum transport equations |

(e.g., Burnett, Grad, OBurnett, O13), DSMC | ||

Kn > 10 | Free molecular flow | Collisionless Boltzmann equation, DSMC |

The present review mainly focuses on the continuum theories that are formulated for modeling flows especially in the slip and the transition regimes. In these flow regimes, since the conventional continuum approach of flow modeling no longer works, it becomes imperative to start from the fundamental Boltzmann kinetic equation to model such strong nonequilibrium flows. Although direct simulation Monte Carlo (DSMC) [11,12] technique provides the solution of the Boltzmann equation, these simulations are often time-consuming or expensive. An alternate way is to build *higher order continuum transport equations* starting from the Boltzmann equation, which form a super-set of the Navier–Stokes equations.

An interesting question can be raised: If DSMC technique provides an accurate solution of the Boltzmann equation over the range of Knudsen number of interest, why to put extra effort in formulating higher order continuum theories? There are two primary reasons for favoring the theoretical equations over the DSMC technique. The first reason is the time-consuming and computationally expensive nature of the DSMC simulations. As we know, the DSMC technique, invented by Bird [11,12], is a mesoscopic, particle-based technique, which simulates the movement and the collision of the particles. As this technique operates at the molecular level, they prove to be expensive, especially in the slip and transition regimes where a large number of molecules need to be simulated to obtain meaningful results. The second reason is that the final results from the DSMC simulations provide little insights into the physical understanding of the process. The human mind has been trained to rely on the mathematical equations and their solutions for a better understanding of the physical phenomena. For example, in the Navier–Stokes equations, the physical significance of each term is well-known and how these terms manifests itself in the final analytical solutions for some of the model flow problems can be appreciated more readily than having the final solution without an analytical basis. The nonequilibrium flows are abundantly rich in flow physics with several nonintuitive flow phenomena [13–19] that are difficult to understand from the final results of the DSMC simulations. When we model such flows using higher order continuum transport equations, we get an idea about the physical effects and interesting interplay happening between different forces/effects that are typical in nonequilibrium flows. Hence, we feel that one cannot abandon the theory at the expense of numerical technique and a great amount of research still needs to be carried out in the formulation of higher order continuum theories.

The Chapman–Enskog-based Burnett equations [20,21] and Grad moment method-based 13 moment equations [22,23] form the two most important sets of higher order continuum transport equations. The fundamental Boltzmann kinetic equation of the kinetic theory of gases forms the basis in both these approaches. The solution of the Boltzmann equation is a single particle distribution function (*f*), which is a function of seven independent variables. This high dimensionality and nonlinear integro-differential nature of the equation makes it extremely difficult to solve the Boltzmann equation and obtain the exact form of the distribution function. The equation can be solved only for the simplest of the cases, for example, for a homogeneous system of gas at rest with no gradients of thermodynamic variables, one can solve the equation to obtain the well-known Maxwell–Boltzmann distribution. However, for the systems out of equilibrium, it is virtually impossible to solve the Boltzmann equation and obtain the exact form of the distribution function.

In the Chapman–Enskog method [20,21], the central idea is to represent the distribution function in an infinite series in terms of Knudsen number around a local equilibrium state (Fig. 1). The series needs to be truncated at some point to determine the form of the distribution function. At the zeroth and first approximation, inviscid Euler and Navier–Stokes equations are obtained with explicit expressions for absolute viscosity and thermal conductivity. This systematic derivation of Euler and Navier–Stokes equations along with the expressions for transport coefficients can be considered as a major success of the Chapman–Enskog method. As we go one step further, we obtain Burnett equations at the second approximation. Essentially, the linear constitutive relations obtained at the Navier–Stokes level are appended with several nonlinear terms involving products of derivatives of velocity, pressure, and temperature. However, the success that is observed at the zeroth and the first approximation is not replicated in the case of Burnett equations. The primary reasons are the unstable nature of the equations (see Bobylev's instability [24,25]) and the possible violation of the second law of thermodynamics [26,27]. Furthermore, since the equations contain third-order derivative terms, they require additional boundary conditions for their complete solution, which are generally unknown. In order to overcome the limitations of Burnett equations, several variants of Burnett equations such as augmented Burnett equations [28], regularized Burnett equations [29], and generalized Burnett equations [30], have been proposed in the literature.

The moment method [22,23] is a simpler yet powerful method proposed by Grad. In this method, the solution of the Boltzmann equation is sought by expanding the single particle distribution function in terms of tensorial orthogonal Hermite polynomials around a local equilibrium state (Fig. 1). Grad hypothesized that by including the higher order moments, the stress tensor *p _{ij}* and heat flux vector

*q*in the primary variables set, it is possible to achieve accurate description of nonequilibrium flows. Hence, we do not stop at the conservation equations of mass, momentum, and energy like we do in Burnett equations; rather we generate separate transport equations for the stress tensor and the heat flux vector by taking moments of the Boltzmann equation. As a result, in the 13 moment method, we have thirteen primary variables:

_{i}*ρ*,

*u*,

_{i}*T*,

*p*and

_{ij,}*q*(1 + 3 + 1 + 5 + 3 = 13) and the problem of closure that is encountered in the momentum and energy equations in the Burnett framework is now shifted in the transport equations for

_{i}*p*and

_{ij}*q*. Once the form of the distribution function is determined by solving the Boltzmann equation, the unknown higher order moments can be evaluated. However, the closure is still one step away, since the production terms arising from the collision integral in the transport equations for

_{i}*p*and

_{ij}*q*still need to be evaluated. These terms signify the net production of stresses and heat fluxes due to collisions and evaluation of these terms is complex with the closed-form expressions known only for the Maxwell molecules [3,31]. Hence, we can safely say that the moment equations are known in closed form only for Maxwell molecules. These types of production terms do not arise in Burnett hydrodynamics since these terms are exactly zero for the first five moments. The vanishing of these terms signify that there is no generation of mass, momentum, and energy during the intermolecular collisions. Similar to Burnett equations, the Grad 13 moment equations also have their own set of limitations. The boundary conditions for full three-dimensional Grad 13 moment equations are unknown; hence, the moment theory is also far from complete. Further,

_{i}*H*-theorem, an equivalent of second law of thermodynamics, is not proven for the equations. The equations form a hyperbolic set of equations; as a result, the problem of subshocks arises in shock profiles [32]. To overcome these limitations, regularized 13 [3,33–35] (R13) and regularized 26 [36] (R26) moment equations are some of the variants proposed in the literature.

The Burnett and the moment equations have a mathematical basis but the thermodynamic consistency of the equations is not yet established. This can be the potential source of the limitations that plagues these equations. To have a sound and meaningful theory, it is necessary to formulate a theory that adheres to the thermodynamic principles. In this context, the Onsager-consistent approach that has been shown as the third branch emanating from the Boltzmann equation in Fig. 1 appears to be a promising approach. Completely independent of the Chapman–Enskog and the moment method, this approach has strong roots in the nonequilibrium thermodynamics. A carefully derived distribution function satisfying the Onsager's reciprocity principle is employed to obtain thermodynamically consistent Burnett-*like* and Grad-*like* equations known as the OBurnett (OBurnett) and Onsager 13 (O13) moment equations.

In this work, we mainly focus on the Burnett hydrodynamics and the results for some of the model flow problems. Since the research area of higher order continuum transport equations is highly involved mathematically, we include a brief overview of the important steps involved in the formulation of such theories. For this purpose, we first introduce the Boltzmann equation followed by the derivation of the basic conservation laws in Sec. 2. The conventional form of Burnett equations is then introduced in Sec. 3 followed by the Bhatnagar–Gross–Krook (BGK)-Burnett equations in Sec. 4. The relative merits and demerits of both these equations are touched upon. This is followed by the introduction of the OBurnett equations in Sec. 5 where a brief description about the underlying derivation is included. The recent results obtained with the OBurnett equations for force-driven plane Poiseuille flow and normal shocks are discussed in Secs. 6 and 7, respectively. Comparison with the conventional Burnett equations is made wherever possible in order to bring forward the advantages of the OBurnett equations. The recently proposed Grad's second problem, which has the potential to serve as an ideal test case, is introduced in Sec. 8. Finally, important remarks are made in Sec. 9 followed by concluding remarks in Sec. 10.

## 2 Boltzmann Equation

In the derivation of higher order hydrodynamic theories, one has to start from the fundamental Boltzmann equation for the kinetic theory of gases. The equation describes the evolution of the thermodynamic system, which is not in a state of equilibrium. It does not take into account the position (**x**) and velocities (**c**) of the individual molecules, but considers the probability of finding a molecule in a very small region of phase space ($dxdc$). Here, phase space represents a six dimensional space formed by three space coordinates (*x*, *y*, *z*) and three components of molecular velocity ($cx,cy,cz$).

Here, *F _{k}* is the external force per unit mass acting on the molecules and the term $J(f,f1)$ represents the collision integral. We have intentionally not opened up the collision integral because of its inherent complications, however interested readers can refer to standard texts on the kinetic theory of gases [21,37,38] for explicit form of the collision integral.

*ρ*,

**u**, and

*T*at the macroscopic level. The information embedded in the distribution function is far too detailed to be of any real use. Further, it is particularly difficult to assign the initial and boundary conditions if we are solving problems at the level of the distribution function. Hence, the macroscopic model with five independent variables (

*ρ*,

**u**, and

*T*) in three-dimensional physical space is a far better choice than the microscopic model with a single variable $f(x,c,t)$, which is a function of seven independent variables. This shift from microscopic to macroscopic can be achieved by taking moments of the distribution function as

*m*is the mass of each molecule,

*k*is the Boltzmann constant,

_{B}**c**is the molecular velocity, and

**C**is the peculiar velocity defined as

*P*can be obtained by taking second-order moment of the distribution function. Please note that the pressure tensor

_{ij}*P*and stress tensor,

_{ij}*τ*as defined in standard fluid dynamics textbooks, have opposite sign, $Pij=\u2212\tau ij$

_{ij}*q*can be obtained by taking third-order moment of the distribution function and then contracting it as

_{i}*P*as defined in Eq. (5) can be decomposed as

_{ij}*p*is traceless, i.e., sum of the diagonal elements is zero (

_{ij,}*p*= 0). Contracting the pressure tensor equation (Eq. (5)), an expression for the thermodynamic pressure is obtained as

_{ii}where we have used $CiCj\u221213C2\delta ij=C\u2329iCj\u232a$. The angular bracket notation denotes the symmetric and trace-less part of the tensor.

Equation 12 represents the balance of internal energy and for monatomic gas, we have $\u03f5=cvT=3RT/2$ During the intermolecular collisions, there is no generation of mass, momentum, and energy so that the production terms involving the collision integral become zero. It is important to note that the form of the distribution function is not required while arriving at the conservation equations from the Boltzmann equation. However, the system of equations is not closed with unknowns, stress tensor (*p _{ij}*), and heat flux vector (

*q*), appearing in the momentum and energy equations. To evaluate these unknowns, we need to solve the Boltzmann equation, obtain the form of the distribution function, and then using Eqs. (9) and (6), the stress tensor and heat flux vector can be obtained in terms of primary variables resulting in a closed set of equations. This forms the standard procedure for building hydrodynamic theories starting from the Boltzmann equation. This approach does not employ phenomenological laws such as Newton's law of viscosity and Fourier's law of heat conduction to obtain expressions for the stress tensor and heat flux vector, as is done in the traditional approach. Furthermore, this kinetic theory approach also gives expressions for transport coefficients, absolute viscosity, and thermal conductivity, in terms of molecular parameters. Hence, we can say that everything is derived from the kinetic theory of gases, starting from the conservation laws to the constitutive relations along with the explicit expressions for transport coefficients.

_{i}## 3 Burnett Equations

These conditions strictly imply that only the zeroth approximation ($f=f(0)$) should reproduce the primary variables (mass, momentum, and energy densities) with contribution from the higher order approximations ($f(1),f(2),$ etc.) to be strictly zero.

*μ*is absolute viscosity and

*k*is thermal conductivity. The tensor

*S*is the symmetric and trace-free part of the velocity gradient tensor defined as

_{ij}In the above equation, *δ _{ij}* is a second-order isotropic Kronecker delta tensor. The numerical values of Burnett coefficients ($\omega \u2032s$ and $\theta \u2032s$) appearing in Eqs. (21) and (22) are given in Table 2. They depend on the interaction potential between the molecules and are different for Maxwell and hard-sphere molecules.

Coefficients | Maxwell molecules | Hard-sphere molecules |
---|---|---|

ω_{1} | $4/3[7/2\u2212(T/\mu )(d\mu /dT)]$ | $1.014\xd74/3[7/2\u2212(T/\mu )(d\mu /dT)]$ |

ω_{2} | 2 | 2.028 |

ω_{3} | 3 | 2.418 |

ω_{4} | 0 | 0.681 |

ω_{5} | $3(T/\mu )(d\mu /dT)$ | $0.806(3T/\mu )(d\mu /dT)\u22120.990$ |

ω_{6} | 8 | 7.424 |

θ_{1} | $(15/4)[7/2\u2212(T/\mu )(d\mu /dT)]$ | 11.644 |

θ_{2} | 45/8 | 5.822 |

θ_{3} | –3 | –3.090 |

θ_{4} | 3 | 2.418 |

θ_{5} | $35/4+(T/\mu )(d\mu /dT)$ | $1.627\xd73[35/4+(T/\mu )(d\mu /dT)]$ |

Coefficients | Maxwell molecules | Hard-sphere molecules |
---|---|---|

ω_{1} | $4/3[7/2\u2212(T/\mu )(d\mu /dT)]$ | $1.014\xd74/3[7/2\u2212(T/\mu )(d\mu /dT)]$ |

ω_{2} | 2 | 2.028 |

ω_{3} | 3 | 2.418 |

ω_{4} | 0 | 0.681 |

ω_{5} | $3(T/\mu )(d\mu /dT)$ | $0.806(3T/\mu )(d\mu /dT)\u22120.990$ |

ω_{6} | 8 | 7.424 |

θ_{1} | $(15/4)[7/2\u2212(T/\mu )(d\mu /dT)]$ | 11.644 |

θ_{2} | 45/8 | 5.822 |

θ_{3} | –3 | –3.090 |

θ_{4} | 3 | 2.418 |

θ_{5} | $35/4+(T/\mu )(d\mu /dT)$ | $1.627\xd73[35/4+(T/\mu )(d\mu /dT)]$ |

The Burnett order terms are second-order accurate in Knudsen number whereas Navier–Stokes order terms are first-order accurate in Knudsen number. Hence, it is expected that these terms will play a dominant role for the flows in the slip and the transition regime. As is evident from Eqs. (21) and (22), the Burnett order stress and heat flux terms are much more complex as compared to that of Navier–Stokes order terms (Eqs. (19) and (20)). Specifically, these terms are nonlinear, involves product of velocity gradients and are also dependent on the temperature gradient terms. An interesting point to note is that the temperature gradient terms (*ω*_{3} and *ω*_{5} terms) can give rise to stresses even in the absence of velocity and pressure gradients in the flow. Further, we note the presence of second order derivatives in the Burnett order stress (Eq. (21)) and heat flux terms (Eq. (22)), which upon substituting in the conservation equations for momentum and energy (Eqs. (11)–(12)) become third-order derivatives, thereby requiring additional boundary conditions for their solution.

## 4 BGK-Burnett Equations

*f*) at a point in physical space to a Maxwellian equilibrium distribution function (

*f*). The rate at which this transition happens can be considered to be proportional to the molecular collision frequency (

_{M}*ν*). This is the primary assumption in the BGK approximation and with this assumption, the collision integral is replaced as

The BGK-kinetic model inherits the basic features of the collision integral. The model satisfies the basic conservation equations, has *H*-theorem, and the distribution function reduces to Maxwellian in equilibrium; as such, the basic features of the original collision integral are retained in the BGK-kinetic model. However, being a *model*, it does not describe the complete physics of the binary collisions, and we are bound to lose some kinetic information with this simplification which we will elaborate in Grad's second problem in Sec. 8. Further, the model incorrectly predicts the Prandtl number as unity for monatomic gases. To overcome this limitation, another BGK-type model known as ES-BGK (ellipsoidal statistical) model [41] was proposed.

The Burnett equations evaluated by employing BGK-kinetic model using Chapman–Enskog method are termed as BGK-Burnett equations [42–45]. The structure of the BGK-Burnett equations is same as that obtained in conventional Burnett equations (Eqs. (17)–(22)); however, the Burnett coefficients (*ω*'s and *θ*'s) are different and they are given in Table 3 for Maxwell and hard-sphere molecules.

## 5 OBurnett Equations

In the Chapman–Enskog method, the first-order correction to the distribution function that gives Navier–Stokes equations satisfies the Onsager's symmetry principle [46,47]. However, the symmetry principle is not satisfied by the second-order correction to the distribution function, which gives Burnett equations [48,49]. This can be the source for the unstable nature of the Burnett equations and also underscores the importance of symmetry principle in the formulation of higher order continuum theories.

*X*and fluxes, $\u03d2i$ [4,50–52] as

_{i}where $f(0)$ represents Maxwell–Boltzmann distribution. The subscript *O* in $f|O$ signifies that the distribution function is consistent with the Onsager's symmetry principle. Truncating the series at the first-order correction, this approach faithfully reproduces the Navier–Stokes equations. Projecting the distribution function with second-order correction to the macroscopic level, we obtain the OBurnett equations [4,52]. The derived form of the distribution function at the second-order correction satisfies the linearized Boltzmann equation, compatibility conditions, and *H*-theorem while preserving the Onsager's reciprocity symmetry principle. The complete form of the Onsager-consistent distribution function is given in our earlier works [51,52].

*p*and heat flux vector

_{ij}*q*are obtained by adding the corresponding Navier–Stokes and Burnett terms as (underlined terms represent the Navier–Stokes order contribution)

_{i}*β*and

*g*are given as

Variable | Base equation | Change of variables | |
---|---|---|---|

$pyyOB$ | $pxxOB$ | $x\u2192y;y\u2192x$ | $u\u2192v;v\u2192u$ |

$pzzOB$ | $pxxOB$ | $x\u2192z;z\u2192x$ | $u\u2192w;w\u2192u$ |

$pyzOB$ | $pxyOB$ | $x\u2192y;y\u2192z;z\u2192x$ | $u\u2192v;v\u2192w;w\u2192u$ |

$pzxOB$ | $pxyOB$ | $x\u2192z;y\u2192x;z\u2192y$ | $u\u2192w;v\u2192u;w\u2192v$ |

$qyOB$ | $qxOB$ | $x\u2192y;y\u2192x$ | $u\u2192v;v\u2192u$ |

$qzOB$ | $qxOB$ | $x\u2192z;z\u2192x$ | $u\u2192w;w\u2192u$ |

Variable | Base equation | Change of variables | |
---|---|---|---|

$pyyOB$ | $pxxOB$ | $x\u2192y;y\u2192x$ | $u\u2192v;v\u2192u$ |

$pzzOB$ | $pxxOB$ | $x\u2192z;z\u2192x$ | $u\u2192w;w\u2192u$ |

$pyzOB$ | $pxyOB$ | $x\u2192y;y\u2192z;z\u2192x$ | $u\u2192v;v\u2192w;w\u2192u$ |

$pzxOB$ | $pxyOB$ | $x\u2192z;y\u2192x;z\u2192y$ | $u\u2192w;v\u2192u;w\u2192v$ |

$qyOB$ | $qxOB$ | $x\u2192y;y\u2192x$ | $u\u2192v;v\u2192u$ |

$qzOB$ | $qxOB$ | $x\u2192z;z\u2192x$ | $u\u2192w;w\u2192u$ |

*Fewer boundary conditions:*Noticeable absence of higher order derivatives of velocity and temperature eliminates the need of prescribing additional boundary conditions. Hence, the well-known Maxwell's velocity slip and temperature jump boundary conditions can serve as the boundary conditions for the equations and now, we have the complete theory in the form of OBurnett equations, unlike the conventional Burnett equations.*Unconditionally stable:*When linear stability analysis is performed for the OBurnett equations, the equations reduce to perturbed form of the Navier–Stokes equations. This implies that the OBurnett equations are unconditionally stable, similar to the Navier–Stokes equations. This is not true for the conventional Burnett equations, which are shown to be unstable for small wavelengths.*Prandtl number value:*The correct value of the Prandtl number is ensured by selecting two different relaxation times for momentum and energy transport in the formulation of the distribution function.*Validity of the equations in the transition regime:*The OBurnett theory does not involve Chapman–Enskog like expansion in terms of Knudsen number; hence, we expect the equations to expand the Knudsen number envelope into the transition regime.

Taking into account these salient features, we have a sound, complete theory in the form of OBurnett equations having a great potential in describing high Knudsen number flows. In Secs. 6 and 7, we compile the results of OBurnett equations for some of the benchmark flow problems which will help to establish the validity of the equations.

## 6 Force-Driven Plane Poiseuille Flow

Force-driven plane Poiseuille flow is one of the few flow problems for which explicit solution of the Navier–Stokes equations is known. Assuming constant transport coefficients, the equations basically reduce to a linear form, which can then be solved analytically to obtain a parabolic velocity profile. However, when the characteristic length scale (distance between the lower and upper plates) becomes comparable to the mean free path, the Navier–Stokes equations are inadequate in capturing certain nonequilibrium effects peculiar to this flow problem as demonstrated in the literature [7,8,54–56]. Particularly, the pressure and temperature profiles as predicted by the Navier–Stokes equations are qualitatively different when compared with the DSMC results.

As in Poiseuille flows, the plates are considered to be infinitely long as shown in Fig. 2 and the flow is assumed to be steady. The origin is selected such that the lower plate is located at $y=\u2212L/2$ while the upper plate is located at $y=L/2$. Both the plates are maintained at constant temperature *T _{R}*. The flow is driven by an external body force per unit mass, $Gi={a,0,0}$,

*a*being constant acceleration acting in the

*x*-direction.

**u**is assumed to have only

*x*-component (

*u*) while the other two components,

*v*and

*w*are taken as zero. Further, all the flow variables are functions of

*y*-direction only. With these assumptions, the continuity equation is satisfied and the momentum and energy conservation equations can be written as

Equations (34*b*) and (34*c*) can be readily integrated to obtain $p+pyy=C1$ and $pyz=C2$ where *C*_{1} and *C*_{2} are the integration constants. Note that these expressions are obtained without substituting the constitutive relations for stress tensor, and hence they remain same irrespective of the equations set. Before reporting the results of the Navier–Stokes and other higher order transport equations, we analyze the problem using the inviscid Euler equations. For Euler equations, dissipative effects are absent, i.e., *p _{ij}* = 0 and

*q*= 0 giving

_{i}*p*=

*C*

_{1}and $C2=0$. Substituting

*p*= 0 in the

_{xy}*x*-momentum Eq. (34

*a*) gives

*a*=

*0, which is a contradiction since the flow is driven by some nonzero external force. This implies that the Euler equations do not provide a meaningful solution for the flow problem driven by an external force.*

*Navier–Stokes equations:*For the Navier–Stokes equations, the constitutive relations for the stress tensor and heat flux vector for the problem considered are

On solving the *y*-momentum equation, it is readily seen that the pressure is a constant. Treating density and transport coefficients (*μ*, *k*) as constants, the equations are reduced to a linear form, which can then be solved to obtain the well-known parabolic velocity and quartic temperature profiles. In the dilute gas regime, the treatment of density and transport coefficients as constants is no longer valid. Hence, density needs to be replaced in the *x*-momentum equation using the ideal gas equation ($p=\rho RT$) and the temperature dependence of the transport coefficients is accounted for by considering the molecules as hard-sphere molecules ($\mu ,k\u221dT1/2$).

Without actually solving the equations, some important remarks can be made from the constitutive relations and the reduced form of the equations. For example, all the normal stresses ($pxx,pyy,pzz$) and heat fluxes in the *x* and *z* directions (*q _{x}*,

*q*) are zero and pressure is uniform throughout the flow. Further, the nature of the equations reveals a quartic temperature profile with a maximum at the center. Comparison with the DSMC results for Knudsen numbers falling in the slip and early transition regime [7] show the limitations of the Navier–Stokes equations. These limitations can be summarized as:

_{z}The DSMC results predict a noticeable variation in the pressure profile in the transverse (

*y*) direction [7,55,57,58], whereas Navier–Stokes equations give a uniform pressure throughout the flow field.The temperature profile as observed in DSMC simulations has two maxima with a characteristic dip at the center [7,8,55,58–60]. This nonmonotonic temperature profile cannot be captured by the Navier–Stokes equations.

The DSMC results show that there is a non-zero, non-negligible tangential heat flux component (

*q*) in the opposite direction to the flow in the bulk region [7,54–56,58,60,61]. Further, the presence of normal stresses ($pxx,pyy,pzz$) in the flow is also confirmed by the DSMC results. The Navier–Stokes equations predict all these fluxes to be zero which are at odds with the DSMC results._{x}

With these observations, it is evident that for rarefied gas flows in the slip and transition regimes, there are several nonequilibrium effects observed in the flow which cannot be captured with the conventional Navier–Stokes equations. We now report the results of the conventional Burnett [7] and OBurnett [9] equations for the same problem.

*Conventional Burnett equations:*While performing analysis using the conventional Burnett equations, the reduced form of the

*x*-momentum and the energy equations remain unchanged. However, the Burnett order contribution to the normal stress

*p*is not zero implying the variation of pressure in the transverse direction. The reduced form of the conservation equations read as

_{yy}*OBurnett equations:*For the OBurnett equations, the constitutive relationships for stresses and heat fluxes in terms of regular variables (

*ρ*,

*u*,

*p,*and

*T*) are

For the conventional Burnett as well as OBurnett equations, it is important to note that the higher order Burnett contribution to shear stress *p _{xy}* and normal heat flux

*q*is zero. On the contrary, normal stresses ($pxx,pyy,pzz$) and tangential heat flux (

_{y}*q*) have contribution only from the Burnett order. The results for conserved (Fig. 3) and nonconserved variables (Fig. 5) at Knudsen number equal to 0.1 according to the Navier–Stokes equations, OBurnett equations, [9] conventional Burnett equations [7], and Nonlinear Coupled Constitutive Relation (NCCR) theory [60] are compared against the DSMC [7] and Molecular dynamics (MD) results [58]. The problem was solved as an initial value problem with the same initial values for

_{x}*p*,

*u,*and

*T*for Navier–Stokes, conventional Burnett [7], and OBurnett equations [9]. In this way, the issue of prescribing accurate boundary conditions is bypassed and the conclusions can then be considered to be independent of the wall boundary conditions.

Out of all the results, the results for pressure profile (Fig. 3(c)) are interesting and warrant more attention. As discussed previously, the Navier–Stokes equations essentially give uniform pressure field. On the other hand, higher order continuum theories are able to resolve the pressure field accurately in the bulk region. However, some discrepancy is observed in the near wall region where conventional Burnett equations predict a peculiar bi-modal profile; also reported by DSMC results of Uribe and Garcia [7]. This result is interesting in a sense that profiles for other thermodynamic variables (density and temperature) do not exhibit such kind of bi-modal behavior in the near wall region. It should be noted that the differential equation for pressure as obtained in conventional Burnett equations has one free parameter (*b*_{2} in Eq. (16) of Uribe and Garcia [7]) for which we need the value of second derivative of pressure at the center. Typical DSMC data, especially the pressure data, has a lot of scatter and hence, accurate determination of $d2pdy2|s=0$ is not possible. The value of this free parameter was fine-tuned in Uribe and Garcia [7] so that the Burnett results for pressure follow the DSMC results. Later, Jadhav et al. [9] reported that the pressure profile is extremely sensitive to *b*_{2} where a slight variation from 0.75 to 0.81 changes the qualitative behavior of pressure as shown in Fig. 4.

For OBurnett equations, the first-order differential equation for pressure has no free parameters. Hence, rather than constructing the pressure profile so as to have better agreement with the benchmark results, it was allowed to take its natural course. The numerical results for pressure show that the profile is monotonic and in good agreement with the MD results of Rana et al. [58] This contrasting result for pressure in the near wall region: bi-modal behavior as per DSMC versus the monotonic profile as per MD is difficult to settle and needs further studies to resolve this issue. Rana et al. [58] conjectured that the DSMC solution in the near wall region is contaminated by the computational errors, thereby resulting in inaccurate resolution of the gas-wall interaction. This observation is supported by Shah et al. [62] where the authors report largest computational error in the cells closest to the wall. As per these two studies, it looks like the stochastic, mesoscopic DSMC technique is somehow unable to resolve the gas–wall interaction accurately and there is a need to carry out microscopic, deterministic MD simulations. For temperature (Fig. 3(d)), both DSMC and MD techniques predict a noticeable dip at the center, which is not captured by conventional Burnett as well as OBurnett equations. As shown by Xu [59], this dip is captured by super-Burnett equations, which are third-order accurate in Knudsen number. Since both conventional Burnett and OBurnett equations are second-order accurate in Knudsen number, it is not surprising that both these theories are unable to capture this dip.

For stresses and heat fluxes (Fig. 5), it is observed that the results of conventional Burnett as well as OBurnett equations are in good quantitative agreement with the benchmark results in the bulk region, whereas some discrepancy is observed in the near wall region. We would like to highlight some of the subtle points here. The DSMC results predict a nonzero value for normal stresses and tangential heat flux at the center. These results are somewhat surprising since they question some of the fundamental principles of nonequilibrium thermodynamics. In a sense, the DSMC results are predicting fluxes even in the absence of thermodynamic forces which is at odds from the thermodynamical perspective. For the problem considered, there is an inherent symmetry in the problem so that the profiles of conserved variables ($u\u2217,\u2009\rho \u2217,\u2009p\u2217$, and $T\u2217$) are symmetric about the centerline axis (*s *=* *0) implying that their first order derivatives at *s *=* *0 to be zero. With no gradients at the center, this condition suggests that the thermodynamic forces vanish at the center and therefore, thermodynamic fluxes should also be essentially zero, i.e., $pij(0)=0$ and $qi(0)=0$. However, the predictions of the DSMC simulations seem to suggest otherwise though these values can be considered to be negligible. This is a difficult issue to settle and we feel that while interpreting the computational results, an emphasis should also be placed on the first principles of thermodynamics.

To conclude, the results of different continuum theories for force-driven plane Poiseuille flow were analyzed. This problem demonstrates the emergence of the nonequilibrium effects when the flow is in the transition regime even for a seemingly simple one-dimensional flow problem, which cannot be captured within the Navier–Stokes framework. Some important remarks regarding the monotonocity of the pressure and temperature profiles were made. We also emphasized that the thermodynamical principles should be taken into account while interpreting the computational results.

## 7 Normal Shock Wave Flow Problem

The accurate description of shock waves has remained an elusive problem in theoretical fluid dynamics. Till date, the internal structure of the shock has not been resolved accurately within the continuum framework. Inside the shock, thermodynamic properties change rapidly over a very short distance (of the order of few mean free paths). Because of these steep gradients, the Navier–Stokes equations, which are based on linear constitutive laws, are unable to discern the flow physics inside this strong nonequilibrium region. This inability of the equations to describe the shock profiles accurately was perhaps the primary driving force in the quest of formulation of better higher order continuum theories. So when any new theoretical formulation comes up, the theory is put to test for the shock wave flow problem to establish its validity. As there are no physical wall boundaries, the intricacies of gas–wall interactions do not arise. Hence, the problem of prescribing wall boundary conditions is resolved and the well-known Rankine–Hugoniot conditions serve as the boundary conditions.

To establish the accuracy of any new theory, the theoretical predictions needs to be compared and benchmarked against the experimental/DSMC/MD data available in the literature and for shock waves, there is no dearth of this data. Specifically, experimental measurements of density across the shock performed using sophisticated electron beam technique (see Alsmeyer [63] and references therein) are recorded in the literature. Although experimental data for other thermodynamic variables are not available owing to the difficulty involved in the measurements, the DSMC and MD simulation techniques furnish accurate data for all the flow field variables that can be considered as reliable. Going one step further, at the microscopic level, experimental measurements of the particle distribution function at different points inside the shock are also available in the literature [64–66]. With all these aspects, it suffices to say that the normal shock wave flow problem serves as an ideal test case for gauging the accuracy of different higher order continuum theories. As we will see, this test problem pushes the continuum theory to an extreme limit and exposes the limitations of the theory, if any.

In a typical shock wave, we have two equilibrium states: the cold region before the shock ($x\u2192\u2212\u221e$) is known as upstream state ($S0$) while the hot region after the shock ($x\u2192\u221e$) is known as downstream state ($S1$) as shown in Fig. 6. The upstream region characterized by density *ρ*_{0}, velocity *u*_{0,} and temperature *T*_{0} is supersonic and is a more rarefied region than the downstream region. The downstream region is subsonic and characterized by *ρ*_{1}, *u*_{1,} and *T*_{1}. Across the shock, the velocity transforms from supersonic to subsonic, whereas density and temperature increases rapidly across this narrow region of the shock. The Mach number (Ma) and the Knudsen number (Kn) are the two important nondimensional numbers that characterize the shock profiles. The Mach number (ratio of velocity of the shock and adiabatic sound velocity) is defined based on the upstream quantities and is greater than unity. In a typical shock wave, the Knudsen number is of the order of unity and falls well within the transition regime.

*x*-component of the velocity vector, stress tensor, and heat flux vector, represented as

*u*,

*σ*, and

*q*, respectively. The time dependence of these quantities can be removed if we solve the problem in the frame of reference moving with the shock. Accordingly, the mass, momentum, and energy conservation Eqs. (10)–(12) are reduced to

*ρ*

_{0}and velocity

*u*

_{0}[67] as

At this point, we substitute the constitutive relationships as per the Euler, Navier–Stokes, Burnett, and OBurnett equations and the analysis differ accordingly. It is easy to see that the Euler equations with no dissipative effects (*σ* = 0 and *q *=* *0) does not give a continuous curve joining the two equilibrium states. This implies that the inviscid Euler model does not provide a good theoretical framework, which we also observed in the force-driven Poiseuille flow problem.

*a*) and 44(

*b*)) with that of Burnett relations (Eq. (24) of Uribe et al. [68]), we observe a relatively simpler structure and absence of problematic second-order derivative terms in the case of OBurnett equations. Substituting Eqs. 44(

*a*) and 44(

*b*) into the reduced form of the integral conservation equations for momentum (43

*b*) and energy (43

*c*), an implicit and nonlinear dynamical system of two ordinary differential equations of order one is obtained. Once this system is solved for velocity and temperature, shock profiles for density, heat flux, and normal stress can be easily obtained. The final system of equations with Rankine–Hugoniot conditions serving as boundary conditions is

The coupled ordinary differential equations that we get after substituting constitutive relations for normal stress and heat flux are highly nonlinear in nature. In addition, it is necessary to account for the temperature dependence of the transport coefficients in the governing equations, further increasing the complexity of the equations. Hence, obtaining shock profiles using analytical tools is highly improbable; although some analytical works using Navier–Stokes equations can be found in the literature [69–73] either for the case of constant transport coefficients or for a particular value of Prandtl number. Linearizing the governing equations is one way to reduce the complexity of the equations, which may aid in obtaining analytical solutions for some flow problems. However, for shock waves, the region inside the shock being strongly nonequilibrium, any attempt to linearize the shock problem will kill the essence of the problem. Therefore, it becomes imperative to attack the shock problem using the numerical tools. In this context, the work of Gilbarg and Paolucci [74] is quite significant where they obtained the shock profiles numerically using the Navier–Stokes equations by accounting the correct value of Prandtl number along with the temperature dependence of transport coefficients. The nature of the upstream and downstream equilibrium states was clearly identified as unstable and saddle nodes, respectively; mathematical analysis of these two states is discussed at length in Garcia-Colin et al. [27] It was suggested to perform numerical integration from the downstream state ($S1$) to the upstream state ($S0$) when the problem is treated as an initial value problem.

The solution curves of the Navier–Stokes equations in the phase space (*u*, *T*) for Ma =* *2 are shown in Fig. 7. The curve connecting the downstream point $S1$ and the upstream point $S0$ is known as heteroclinic trajectory and represents the exact solution. It is interesting to note that many curves are passing through the upstream point $S0$, whereas only two pass through the downstream point $S1$. This implies that no approximate solution from $S0$ will ever reach $S1$ while all solutions from $S1$ will reach $S0$. Hence, if one has to obtain an approximate solution, it becomes necessary to carry out the integration from the downstream point $S1$ to the upstream point $S0$. This method of numerical integration in a negative sense (*x* goes from positive to negative values) proposed by Gilbarg and Paolucci [74] has been followed in several works [53,67, 68,75–77] when the problem is solved as an initial value problem.

The treatment of normal shock within the Navier–Stokes framework has been well documented in the literature [69,72–74,78]. With stability and thermodynamic consistency of the Navier–Stokes equations well established, the theoretical results—smooth shock structures and existence of a heteroclinic trajectory at all Mach numbers [74], was not surprising. In these aspects, it seems that we have a sound mathematical theory for shock waves in the form of Navier–Stokes equations. However, the problem arises when its predictions for hydrodynamic field variables and shock parameters like shock asymmetry, inverse shock thickness are compared against the benchmark results. Even at moderate Mach numbers (Ma* *>* *1.3), the quantitative agreement of the Navier–Stokes results with the benchmark results cannot be said up to the mark.

The shock wave flow problem has also been solved using conventional Burnett equations and interested readers can refer to literature [27,68,79–81]. We would like to highlight some of the important fundamental drawbacks of conventional Burnett equations, which come to the fore when applied to shock waves. First, given the complicated nature of the equations, it is difficult to establish the thermodynamic consistency of the Burnett equations theoretically, like we do for the Navier–Stokes equations. In fact, the works of Comeaux et al. [26] and others [43,53] for normal shocks show that the Burnett equations give negative entropy production in the upstream region of the shock, thereby violating the second law of thermodynamics, as shown in Fig. 8. With regard to the heteroclinic trajectory, Uribe et al. [82] showed that such a trajectory does not exist for the conventional Burnett equations above a certain critical Mach number (Ma* *>* *2.69). Interestingly, this is the critical Mach number above which the equations show negative entropy production in the upstream region of the shock. For Ma =* *3, the orbital structures in the phase space (*u*, *T*) near the upstream point $S0$ are shown in Fig. 9. For conventional Burnett equations, the derivatives of velocity and temperature become absurdly high leading to oscillations in the orbital structures. As a result, the phase trajectory never reaches $S0$ indicating nonexistence of heteroclinic trajectory. These oscillations are greatly amplified for higher Mach numbers and it can be safely concluded that the Burnett equations lose validity for Ma* *>* *2.69. The doubts about the conventional Burnett equations are further cascaded when one considers the unstable nature of the Burnett equations to small wavelengths. This instability, generally referred to as the Bobylev's instability [24,25], restricts to obtain numerical solutions on finer meshes.

The results of the OBurnett equations for normal shock have been recently obtained [53,83] and the theory resolves all the drawbacks associated with the conventional Burnett equations. For the shock waves, the OBurnett equations form a thermodynamic consistent set of equations as shown mathematically in Jadhav and Agrawal [53]. This point is illustrated in Fig. 8 where OBurnett equations predict positive entropy generation throughout the shock for two Mach numbers ($Ma=2.69,5$), thereby satisfying the second law of thermodynamics, similar to the Navier–Stokes equations. The mathematical proof of the existence of heteroclinic trajectory for the Navier–Stokes equations as given by Gilbarg and Paolucci [74] can be extended to the OBurnett equations to show that such a trajectory exists for the OBurnett equations as well. This point becomes clear from Figs. 9 and 10, which confirms the existence of the heteroclinic trajectory for the OBurnett and Navier–Stokes equations for Ma = 3, 7. Further, the stability analysis shows that the OBurnett equations are unconditionally stable, unlike the conventional Burnett equations [4,52]. With these results, we can conclusively say that the shock wave flow problem, while exposing the limitations of the conventional Burnett equations, help to establish the fundamental aspects of the OBurnett equations and they seem to be a superior superset of the Navier–Stokes equations.

Having discussed about the fundamental aspects of OBurnett equations, more pertinent question is how the OBurnett results compare against the benchmark results. This can be quantified by comparing the shock profiles for hydrodynamic field variables as well as the orbital structures in the phase space (*u*, *T*) with the benchmark data. Although the first option of comparing shock profiles seems more appropriate, we would like to bring forward one major caveat associated with it. While comparing the shock profiles, the choice of origin is of paramount importance. Sometimes, the origin is selected such that the density at the origin gives the average of upstream and downstream values [63,84–86], whereas in some cases, the selection of origin can be based on velocity parameter [53,67,68]. Based on these two different selection parameters, the other shock profiles are translated accordingly and the final comparison with the DSMC/MD benchmark data based on these two criteria will be different. In order to overcome this ambiguity, it is advisable to compare the underlying orbital structures in the phase space (*u*, *T*) with the benchmark results. When working in the phase space, the dependence on the space coordinate is removed and a fair comparison with the benchmark results is then possible. In our view, comparison of orbital structures can be a better yardstick for benchmarking the results since other quantities ($\rho ,q,\sigma $) are basically derived quantities from the velocity and the temperature data. Comparison of orbital structures at Ma* *=* *7 (Fig. 10) and shock profiles at Ma = 3 (Fig. 11) with the DSMC results imply that the OBurnett equations perform much better than the Navier–Stokes equations, especially in the upstream region of the shock.

*t*. Integrating over the physical and velocity space, we obtain total number of molecules (

*N*) as

In the second step, we have normalized the distribution function by the total number of molecules. It can be easily recognized that this normalized quantity *f*/*N* denotes the probability of finding a molecule in a phase space volume $dcdx$. As we know, the probability of finding a molecule inside the volume *V* cannot be negative; any violation of this aspect can question the underlying derivation of the particle distribution function. This drawback seems to arise in Grad 13 moment method [22,23] where the Grad distribution function assumes negative values for Ma* *>* *4 as shown in [3,85]; thereby invalidating the 13 moment equations at higher Mach numbers.

As seen from the results, it is safe to say that the shock wave flow problem underscores the superiority of the OBurnett equations over the conventional Burnett equations. This success of the OBurnett equations can be attributed to the underlying derivation method, which essentially complies with the principles of nonequilibrium thermodynamics. As a result, we have much more sound and physically consistent Burnett order constitutive relationships for the stress tensor and the heat flux vector. We can also say that the problematic terms that are present in the conventional Burnett equations do not arise in the OBurnett constitutive relations. To elaborate more on this point, it is important to know that the original Burnett equations and the conventional Burnett equations are the two sets of Burnett equations proposed in the literature. The equations as derived by Burnett [20,87] involved material derivatives terms of velocity gradient tensor and temperature gradient in the Burnett order constitutive relations and are known as original Burnett equations. Chapman and Cowling [21] used an order of magnitude analysis and replaced the material derivatives terms with the spatial gradients using the inviscid Euler equations. This form of equations being relatively simpler, gained more adoption and were called as conventional Burnett equations. These particular material derivative terms were shown to generate negative entropy generation in the upstream region of the shock, as shown in Comeaux et al. [26] Further, some works [42,88] also suspect entropic inconsistency due to the second derivative term of density and suggested to drop this term. An order of magnitude analysis as performed in Jadhav et al. [83] also points that these terms can be neglected at higher Mach numbers. In OBurnett constitutive relations, none of these problematic material derivative term and second-order derivative terms are present and the issues that surface in Burnett theory are not encountered in the OBurnett theory.

A case of very strong shock with Mach number equal to 134 has been studied extensively [53,67,68,84,89–91] and can be considered as an ultimate test for any higher order continuum theory. Although specific reason for selecting this particular Mach number is not known to us, the availability of accurate MD and DSMC results made this case more prevalent in the literature. Since the Burnett equations lose validity for Ma* *>* *2.69 and the results of moment equations and their variants are also questionable for Ma* *>* *4, none of these theories can be applied for this case. In our recent work [53], with no tweaking of the OBurnett equations in any way, it was shown that the equations give smooth shock structures, a clear heteroclinic trajectory exists (Fig. 12) and the entropy generation is always positive across the shock (Fig. 13) for Ma* *=* *134. These results imply that there is theoretically no upper Mach number limit for the OBurnett equations unlike the other higher-order continuum theories. From the variation of hydrodynamic fields across the shock (Fig. 14), it is seen that the equations provide substantial improvement over the results of the Navier–Stokes equations. The agreement with the MD results is quantitative in the cold part of the shock whereas further improvement is desirable in the hot part of the shock. With respect to the shock width, the wider shock wave as predicted by MD is described sufficiently well by the OBurnett equations.

There are other works attempted to describe the shocks, below we make a list of some of these works:

Shock profiles obtained by modification of the Navier–Stokes equations [77,96,97]

Nonlinear coupled constitutive relation theory [98–100] proposed by Myong

Shock profiles using the BGK Burnett equations [42,43] based on the BGK-Boltzmann equation

Mott–Smith theory [91,104] where the distribution function is assumed to be a sum of two Maxwellians

The above list is definitely not exhaustive and not sorted in any order of preference. However, we would like to make some remarks here. In some of the works mentioned above, the basic governing equations were modified or some of the parameters like viscosity index were tweaked to have better agreement with the benchmark results. Some theories built from first principles like Burnett equations or Grad equations look sound but they fail in shock wave flow problem. In these aspects, physically meaningful results of OBurnett equations for shock waves enhance belief in the underlying Onsager consistent approach. The theory has a sound fundamental basis; there are no free parameters in the equations with all the OBurnett coefficients coming from the kinetic theory. With no tweaking of the equations in any way, it was observed that the OBurnett equations provided significant improvement over the results of Navier–Stokes equations for normal shocks. As such, we believe that the equations form an accurate set of higher order transport equations and have great potential in describing strong nonequilibrium flows.

## 8 Grad's Second Problem

The derivation of the higher order hydrodynamics theories, as we know now, is a challenging task and the benchmark problems available in the literature to check the accuracy of the equations are also limited. Since the structure of the higher order continuum transport equations is very complex, the focus should be on obtaining accurate solutions for simple flow problems rather than obtaining crude approximations for the complicated problems. Hence, the problems selected as benchmark cases should have a relatively simple definition and be amenable to analytical treatment. This way, greater insights can be obtained about the nature of the equations and the underlying derivation. In this context, we recently proposed a relatively simple and novel problem which we termed it as Grad's second problem [107,108]. The genesis of this problem can be traced to the Grad's work [22] where he proposed two simple problems to check the validity of the moment equations. These problems are: (i) the solution for temporal variation of stress and heat flux in a flow in absence of density, velocity and temperature gradients, and (ii) steady one-dimensional flow of heat in a static gas. Both these problems are considered in an infinite domain with no physical wall boundaries and we term these two problems as *Grad's first and second problem*, respectively.

*x*

_{1}-direction as shown in Fig. 15. According to these conditions, we have

In the absence of physical wall boundaries, we do not encounter the complex gas-wall interactions and the need to prescribe the wall boundary conditions does not arise. Further, since the bulk velocity is zero, many of the terms involving the velocity and their gradients become zero and the governing equations are greatly simplified. This simplistic setup of the problem may eventually lead to analytical solutions for some of the quantities. With gas essentially at rest, the primary interest is in obtaining the steady-state solution for the pressure and the temperature fields. The solution of pressure for different interaction potentials makes this problem particularly interesting. For Maxwell molecules, the exact normal solution of the Boltzmann equation [109] predicts uniform pressure field with no stresses in the flow domain. On the other hand, when we consider hard-sphere molecules, the DSMC solution to the Boltzmann equation predicts pressure gradient in the *x*_{1}-direction along with presence of normal stress in the flow domain [110]. This appears to be an interesting problem to solve within the Burnett hydrodynamics and we report the results of conventional Burnett equations as obtained in Jadhav and Agrawal [107].

where $pij(2)$ is given from Eq. (53). From the structure of the equations, it is seen that there is one-side coupling between the momentum and energy equations. That is, energy equation can be solved first to obtain analytical solution for temperature which can then be substituted in the momentum equation to obtain solution for pressure. We now discuss the results for hard-sphere and Maxwell molecules.

### 8.1 For Hard-Sphere Molecules.

*A*and

*B*are proportionality constants. The energy Eq. (55) can be solved analytically for temperature as

where *C*_{1} is the integration constant.

*x*

_{1}-direction. Substituting these stresses in momentum Eq. (54), the

*x*

_{2}and

*x*

_{3}-momentum equations are identically zero while the

*x*

_{1}-momentum equation gives

*x*

_{1}-momentum Eq. (60), we obtain

Substituting the solution for temperature (Eq. (57)) in Eq. (61), we obtain a third-order ordinary differential equation for pressure, which is solved as an initial value problem. The three necessary initial conditions for pressure: $[p,\u2009dp/dx1,\u2009d2p/dx12]|x1=x\u2217$ and one initial condition for temperature: [*T*]$|x1=x2$ are obtained from the DSMC solution of the Boltzmann equation [110]. The results for pressure and temperature as per the conventional Burnett equations are compared against that of DSMC results obtained by Montanero et al. [110] as shown in Fig. 16. It is evident that the pressure gradient exists in the *x*_{1}-direction, which is captured reasonably well by the Burnett equations. The direction of this pressure gradient is opposite to that of temperature gradient and the magnitude also being smaller. This weak, nonuniform pressure field as predicted by the DSMC simulations for hard-sphere molecules can be considered as a nonlinear effect.

### 8.2 For Maxwell Molecules.

The similar solution was obtained by Santos [111] within the moment framework for Maxwell molecules starting from the one-dimensional form of the Boltzmann equation.

which are in agreement with the exact normal solution (Eq. (62)) of the Boltzmann equation.

When Grad's second problem is attempted within the BGK-Burnett or OBurnett equations, which are based on BGK-kinetic model, some interesting conclusions can be drawn out. For Maxwell molecules, both the BGK-Burnett and OBurnett equations predict results that are consistent with the solution of Boltzmann equation (Eq. (63)). However, the uniform pressure field as suggested by these equations for hard-sphere molecules does not agree with the DSMC results. It seems that the BGK-kinetic model-based equations do not distinguish between the Maxwell and hard-sphere molecules, at least for the Grad's second problem. Even the solution of the Grad's second problem with the BGK-Boltzmann equation (Eq. (28)) predicts similar results for Maxwell and hard-sphere molecules as shown in Santos et al. [112]. We mentioned in Sec. 4 that the BGK-kinetic model has all the basic features of the original collision integral. However, it is demonstrated through this problem that some kinetic information is lost by replacing the original collision integral with the BGK-kinetic model. The solution of the Grad's second problem within the moment framework is limited to Maxwell molecules since the equations are known in exact form only for the Maxwell molecules. As for the solution, Grad 13 and R13 moment equations give results that are consistent with the solution of Boltzmann equation [108]. For hard sphere molecules, the production terms appearing in the transport equations for the stress tensor and the heat flux vector cannot be evaluated in closed form. As such, it is difficult to comment on how the moment equations will fare for the hard-sphere molecules.

Some final remarks regarding the Grad's second problem can be made now. As we have seen, the problem definition is relatively simple with no macroscopic bulk velocity and the main focus is on obtaining the pressure field for different interaction potentials. The transition from uniform to nonuniform pressure field as we switch from Maxwell to hard-sphere molecules makes this problem particularly interesting. The solution of the Boltzmann equation: exact normal solution for Maxwell molecules and DSMC solution for hard-sphere molecules, can serve as the benchmark results with which the theoretical results of the equations can be compared with. As there are no physical wall boundaries in the flow domain, the issue of prescribing boundary conditions does not arise. With all these features, we can say that the Grad's second problem has the potential to serve as an ideal test case for checking the accuracy of the higher-order continuum transport equations. Out of all the theories, we observed that only the conventional Burnett equations were able to predict results that are consistent with the solution of the Boltzmann equation for hard-sphere as well as Maxwell molecules. Especially, the weak, nonuniform pressure field as observed in hard-sphere molecules is captured well by the conventional Burnett equations. On the other hand, the theories based on BGK-kinetic model were not able to capture the nonuniform pressure field for hard-sphere molecules, though they provided accurate results for Maxwell molecules. It was conjectured that there is some loss of kinetic information when original collision integral is replaced by the BGK-kinetic model.

## 9 Discussion

In this section, important comments regarding the boundary conditions, thermodynamic consistency of the equations and the solutions for standard benchmark problems are made with respect to conventional Burnett and OBurnett equations.

### 9.1 Comment on Thermodynamic Consistency and Boundary Conditions

#### 9.1.1 The OBurnett Equations Along With Slip Boundary Conditions Form a Complete Theory.

Prescribing wall boundary conditions is an important next step in the formulation of higher order continuum theories, and unfortunately, it is at this step where most of the higher order hydrodynamic theories fall short. As we know, for a given set of partial differential equations, it becomes necessary to supply the appropriate initial and boundary conditions for the complete solution. When we speak of the Navier–Stokes equations at the continuum level ($Kn\u21920$), the no-slip boundary conditions are appropriate and no extra effort is required to model these fairly obvious no-slip boundary conditions. This is not the case in dilute gas regime (Kn* *>* *0.001) where there is slipping of the molecules at the wall surface. Furthermore, there cannot be a universal set of wall boundary conditions that can be applied for all the higher-order continuum theories, i.e., each set of equations will have its own specific boundary conditions. For conventional Burnett equations, the appearance of higher order derivatives of velocity, temperature, and pressure in the constitutive relations makes it necessary to prescribe additional boundary conditions which are still unknown. This is perhaps the strongest criticism that can be made against the conventional Burnett equations and the future outlook looks grim since very little efforts are seen in the literature to resolve this issue.

The moment equations and their variants have their own set of problems. By treating stress tensor and heat flux vector as the primary variables, we now have to specify the boundary conditions for these variables as well. Some progress is visible in this front where Grad himself presented his own set of boundary conditions for a two-dimensional case [22]. In the case of R13 equations, we can see several iterative attempts [113–116] to arrive at the correct form of the boundary conditions.

For the OBurnett equations, we saw that there are no higher order derivatives in the constitutive relations of the stress tensor and the heat flux vector. This structure of the OBurnett closure relations can be considered remarkable since now we need the same number of boundary conditions as that of the Navier–Stokes equations. In this respect, the well-known Maxwell velocity slip and temperature jump boundary conditions are sufficient and helps to complete the OBurnett theory. Furthermore, the linearized stability analysis of the OBurnett equations proves the unconditional stable nature of the equations, similar to the Navier–Stokes equations. This makes it possible to obtain accurate solutions on finer grids, which is not possible for the unstable conventional Burnett equations. We believe that with these two important aspects, namely, the stable nature of the equations and the known boundary conditions, the idea of carrying out numerical simulations for different boundary value problems using the OBurnett equations within the computational fluids dynamics framework does not look far stretched.

#### 9.1.2 The OBurnett Equations Form a Thermodynamical-Consistent Set of Equations.

The thermodynamic consistency of the higher order continuum transport equations has been difficult to prove given the complicated nature of the equations. For the Navier–Stokes equations, the proof for positive definite entropy production has been well established. However, we do not have a definite answer regarding the thermodynamic consistency of the conventional Burnett equations. In fact, one instance of the equations violating the second law of thermodynamics like the normal shock problem is sufficient to conclude that they are thermodynamically inconsistent. In the moment framework, the *H*-theorem, which is equivalent of second law of thermodynamics is not yet proven for full three-dimensional nonlinear Grad 13 and regularized 13 moment equations. In the derivation of OBurnett equations, the Onsager principle consistent distribution function is constructed such that it satisfies the *H*-theorem. In a way, this strict compliance of the second law of thermodynamics is fulfilled at the level of the distribution function itself and what we eventually get is the thermodynamically consistent set of equations. This point becomes evident in the normal shock problem where we observed that the equations give positive entropy production across the shock at all Mach numbers including the case of infinite shock strength (Ma* *=* *134).

#### 9.1.3 The Results of OBurnett Equations for Benchmark Problems Look Promising.

Having discussed about the stability, thermodynamic consistency, and boundary conditions part for the OBurnett equations, more pertinent question is how these equations perform for the rarefied gas flows in the slip and transition regime. To quantify the performance of the equations, the numerical results of the OBurnett equations for two important bench mark problems: force driven compressible plane Poiseuille flow and normal shock wave flow problem were reported.

The OBurnett equations successfully captured the several nonequilibrium phenomena peculiar to the force driven Poiseuille flow. For example, monotonic pressure variation in the transverse direction, variation of normal stresses and tangential heat flux were successfully resolved by the equations with the possible exception of temperature dip at the center. While interpreting the results, we mainly focused on the variation of conserved (*ρ*, *u _{i}*,

*T*, and

*p*) and nonconserved variables (

*p*and

_{ij}*q*) in the transverse direction. We believe that once these quantities are properly resolved, the prediction of the integral quantities such as mass flowrate, skin friction coefficient, etc., should be rather straightforward. Another interesting point that we raised is the doubtful presence of normal stresses and tangential heat flux at the center of the channel where there are no gradients of velocity, temperature or pressure. The DSMC results predict the presence of these dissipative fluxes, although these values can be considered to be negligible. This behavior seems odd with the first principles of nonequilibrium thermodynamics which strictly says that the dissipative fluxes should be zero in absence of thermodynamic driving forces. It is difficult to come up with justification for this apparent contradiction of the DSMC results with the thermodynamical principles. We believe that great care should be exercised while interpreting the numerical results and more importantly, a simple sanity check should be performed to check for the possible violation of the first principles of thermodynamics.

_{i}The normal shock wave flow problem can be considered as an ultimate test for higher-order continuum theories. We say this because this problem has been instrumental in bringing forward the fundamental limitations of the theory, if any. For example, nonexistence of heteroclinic trajectory, and violation of second law of thermodynamics in the upstream region of the shock were some of the major drawbacks that were exposed by shock problem for conventional Burnett equations. The Grad 13 moment equations were not spared either by the shock problem. The subshocks appearing in the shock profiles and negative values of the Grad distribution function for Ma* *>* *4 put serious limitations on the moment equations. As a result, these theories always have an upper Mach number limit above which they become invalid. This Mach number limit for conventional Burnett, Grad 13 moment and R13 moment equations are 2.69, 1.65, and 4, respectively. It is interesting to note that none of these issues surfaced for the Navier–Stokes equations with practically no upper Mach number limit for the equations. The only thing that troubled the researchers was the inability of the Navier–Stokes equations to resolve the shock profiles accurately even for moderate Mach numbers. This perhaps was the trigger for the quest in search of the accurate higher order continuum theories. However, the dismal performance for the shock waves came as a major setback for the conventional Burnett and moment equations.

In this context, the theoretical and numerical results of the OBurnett equations assume significance. With no tweaking of equations in any way, mathematical evidences were put forward for the following three fundamental aspects:

Smooth shock structures at all Mach numbers

Existence of heteroclinic trajectory at all Mach numbers

Positive definite entropy production across the shock at all Mach numbers

It is important to note that none of the other higher order continuum theories satisfy all these claims. With respect to the hydrodynamic fields across the shock, the results of the OBurnett equations are shown to be in quantitative agreement with the DSMC results in the upstream region of the shock which is supposedly more difficult to capture. Hence, it can be safely said that the OBurnett equations not only improved upon the theoretical aspects of other higher order continuum theories but also provided significant improvement over the results of Navier–Stokes equations at all Mach numbers.

### 9.2 Possible Ways of Solving the OBurnett Equations for Other Problems.

The structure of the higher order continuum transport equations is complex; hence, finding analytical solutions even for simple one-dimensional flow problems is quite a challenging task. There are two methods, the “Agrawal-Singh iterative approach” [4] and the perturbation method, with the help of which analytical solutions of the Burnett equations can be obtained for some of the model flow problems.

The novel *iterative* approach proposed by Agrawal and Singh [117,118] can be applied to flow problems for which analytical solutions of the Navier–Stokes equations is available. As we know, the OBurnett and conventional Burnett equations form a superset of the Navier–Stokes equations, i.e., all the terms of the Navier–Stokes equations are present in the OBurnett equations. Agrawal and Singh used this fact judiciously to “build-up” analytical solutions of the Burnett equations for some of the model flow problems. The important steps in this approach are shown in Fig. 17. The solution of the lower order model, i.e., Navier–Stokes equations is taken as the base solution of the Burnett equations. Since we have explicit closed-form expressions for the primary variables (velocity, density, pressure, etc.), the Burnett order stresses can be readily obtained for the physical problem in hand. These solutions are substituted in the Burnett equations and the order of magnitude of each term is determined and compared with the highest order term. This allows us to identify the significance of each of the additional term from the Burnett equations. The additional terms, if significant, are appended to the governing equations and the augmented governing equations are needed to be solved further analytically. This process is repeated till there are no additional terms in the Burnett equations or the remaining terms are insignificant. If either of these criteria is met, we can say that the convergence is achieved. It is recommended that this method be applied at several points in the flow field, especially in the near-wall region where the solution of the Burnett equations is expected to be substantially different than that of the Navier–Stokes equations. This approach was demonstrated successfully by obtaining analytical solutions of the conventional Burnett equations for microchannel [118] and microtube [117] flows and can be extended to OBurnett equations and possibly other higher order transport equations as well. Another recent analytical study [119] for the microchannel flow used the closed-form expressions for velocity and pressure obtained at the Navier–Stokes level to evaluate the stresses at the Burnett level for Knudsen numbers falling in the slip and transition regimes.

The perturbation method is another promising way of obtaining analytical solutions of the higher order continuum transport equations. Arkilic [120] employed this approach for isothermal gaseous flow in microchannel to obtain analytical solutions of the Navier–Stokes equations for the case when both the Reynolds and Mach numbers of the flow are small. When both these nondimensional numbers are small, the continuity and momentum equations get uncoupled making it possible to obtain closed form solutions for streamwise and transverse components of velocity vector and the pressure distribution across the length of the channel. The similar approach was employed for conventional Burnett equations by Rath et al. [121] to obtain closed form expression for pressure distribution across the entire flow domain. The pressure profile obtained using this approach was found to be in good agreement with the analytical result of Rath et al. [122].

## 10 Conclusions

We presented a concise review of the higher order continuum transport equations belonging to the family of Burnett hydrodynamics with a special focus on the OBurnett equations. Since this particular research area is mathematically involved, we started from the first principles with the introduction of the fundamental Boltzmann kinetic equation. A brief overview of the standard procedure involved in the formulation of higher order transport equations starting from the Boltzmann equation was presented. This kinetic approach of deriving hydrodynamics theories has several advantages over the traditional approach that we regularly see in the standard fluid dynamics reference books. For example, basic conservation laws and closure relations for the stress tensor and the heat flux vector such as Newton's law of viscosity and Fourier's law can be derived systematically with explicit expressions for the transport coefficients. Within the Burnett hydrodynamics, we discussed about the conventional Burnett, BGK-Burnett and OBurnett equations and also shed some light on the underlying derivation. The remarkable features of the OBurnett equations were put forward which help establish the superiority of the OBurnett equations over the conventional Burnett equations. Next, the results of the OBurnett equations for standard benchmark flow problems in rarefied gas dynamics show that they form a reliable set of higher order continuum transport equations.

As a final remark, we would like to say that developing a continuum theory for rarefied gas flows is a formidable task that has stimulated great interest among the researchers. The task becomes even more challenging knowing that the equations need to be supplied with appropriate wall boundary conditions where we need to take into account the intricacies of the gas-wall interactions. As we know, one negative result is sufficient to put serious question marks over the accuracy of the theory. The early results of the OBurnett equations are promising and more importantly, does not indicate violation of the thermodynamic principles. However, they need to be thoroughly tested for more complex problems to firmly establish its validity. It will be interesting to see whether the OBurnett equations withstand the test of time, like the Navier–Stokes equations did.

## Funding Data

Department of Atomic Energy Science Research Council Outstanding Investigator Award (No. 21/01/2015-BRNS/35038).