As the flow becomes rarefied it has been seen that predictions of continuum formulations, such as the Navier–Stokes equations, become inaccurate. These inaccuracies stem from the linear approximations to the stress and heat flux in the viscous flux terms in the Navier–Stokes equations. Hence, it has long been conjectured that the inclusion of higher-order terms in the constitutive relations for the stress and heat flux may improve the predictive capabilities of such continuum formulations. Following this approach, second-order systems of hydrodynamic equations, such as the Burnett and Woods equations, were applied to the shock structure problem. While it was observed that these equations afforded a better description of the shock structure on coarse grids, they were prone to small wavelength instabilities when the grids were refined. The cause of this instability was subsequently traced to the fact that these equations can potentially violate the second law of thermodynamics when the local Knudsen number exceeds a critical limit. This leads to the fundamental question: is entropy consistency achievable in a system of second-order hydrodynamic equations? To answer this question, a novel set of equations, known as the BGK-Burnett equations, is constructed by taking moments of the Boltzmann equation for the second-order distribution function.

The formulation of second-order hydrodynamic equations by moment methods is beset by three hurdles: (i) the highly nonlinear collision integral in the Boltzmann equation needs to be evaluated, (ii) the second-order distribution function does not satisfy the moment closure criterion and (iii) identification of the approximations to the material derivatives in the second-order distribution function that will correctly account for the difference in time scales between the first- and second-order fluxes. The first three terms of the Chapman–Enskog expansion, that defines the second-order distribution function, and the Bhatnagar–Gross–Krook model of the collision integral, form the basis of the BGK-Burnett equations. The entropy-consistent behaviour of the equations depend on the moment closure coefficients and the approximations to the material derivatives in the second-order fluxes. The requirement of moment closure alone, however, results in non-unique closure coefficients and a family of BGK-Burnett equations, from which an entropy-consistent set must be identified. From this family, two sets of BGK-Burnett equations have been considered, and this paper presents the details of the formulation of these two sets of equations, the identification of entropy-consistent approximations to the material derivatives by a novel entropy consistent relaxation technique, and shock structure computations in a monatomic hard sphere gas for a range of Mach numbers.