In the time since the Navier–Stokes equations were introduced more than two centuries ago, their application to problems involving real gas effects has relied on appropriate closure for the mass, momentum, and energy transport fluxes via the constitutive laws. Determination of the corresponding transport coefficients, most readily obtained through generalized Chapman–Enskog theory, requires knowledge of the intermolecular potentials for rotationally and vibrationally excited molecules. Recent advances in computational chemistry provide extraordinary detail of interactions involving rovibrationally excited molecules, offering a means for transport flux closure with unprecedented accuracy. Here, the bracket integrals for rovibrationally resolved molecular states are developed, and the resulting transport flux closure is presented for the rovibrationally resolved Navier–Stokes equations. The accompanying continuum breakdown parameters are also derived as a rigorous metric to establish the range of applicability of the aforementioned equations in flow conditions approaching the rarefied regime.
One of the great triumphs of kinetic theory concerning the Navier–Stokes equations demonstrates that, for a simple gas, the transport coefficients of viscosity and thermal conductivity are independent of number density. James Clerk Maxwell had first predicted and experimentally verified such behavior based on mean free path arguments in the mid-1800s.1 Later in the early 1900s, independent work by Sydney Chapman and David Enskog demonstrated this result,2–4 manifested in the form of the integral equations (derived through what is now known as Chapman–Enskog theory), which ultimately lead to explicit expressions for viscosity and thermal conductivity.
Much of what is understood about the transport coefficients in simple (single species, monatomic) gases cannot be formally extended to molecular gases.5,6 This is particularly true at high temperatures where the excitation of internal energy states becomes important.5,7–10 Developments from Eucken offer a correction for thermal conductivity to include contributions from internal energy transport, and assume, among other things, that rotational energy and vibrational energy are transported through identity via diffusion.11 More formal kinetic theory-based arguments proposed by Wang Chang and Uhlenbeck7,12 and Mason and Monchick13,14 were introduced later to provide an explicit expression for the heat conductivity of polyatomic gases. These results are often utilized for transport flux closure in a multi-temperature (MT) form of the Navier–Stokes equations15–17 in which a distinct temperature is used to describe each of the translational, rotational, vibrational, and electronic energy states of a gas.
While MT models18–22 allow for a non-equilibrium description between internal energy modes, the population distribution within each internal energy mode is assumed to follow an equilibrium Boltzmann distribution.23 Yet in numerous cases, e.g., hypersonic flight, laser energy deposition, and high-enthalpy experimental facilities, to name a few,24–29 non-Boltzmann distributions persist, contributing to—even dominating—the complex transport processes that drive the evolution of a non-equilibrium gas. To overcome this limitation, improved techniques that leverage a state-based or state-to-state (StS) framework10,30–36 have been widely adopted. The StS approach requires accurate intermolecular potentials capable of describing interactions of molecules at or near the dissociation limit,37–39 and has been developed extensively for studying non-equilibrium flows. While much of this development has focused on establishing the chemical kinetics for these systems,32,40–45 full utility of the StS framework for the Navier–Stokes equations requires an equivalently modernized form of the transport fluxes.
This paper presents the full rovibrational state-resolved transport flux closure for the Navier–Stokes equations. The features of the rovibrational StS model are incorporated into the generalized Chapman–Enskog (GCE) theoretical framework, yielding the rovibrationally resolved Navier–Stokes equations. The procedure for expressing transport coefficients in terms of Sonine polynomial expansion coefficients is outlined, and these transport coefficient definitions are used to obtain forms of the “mechanism-based” GCE parameters as well as the “species perturbation” parameter for the prediction of continuum breakdown.
In the StS framework, each rovibrational energy level of a molecule is treated as a separate “pseudo” species. This allows the rotational and vibrational modes of the molecule to remain coupled, unlike commonly used MT models which assume the separability of these modes and describe them according to distinct temperatures. The GCE method begins by specifying a timescale relation between the various collisional processes that best describes the internal energy modeling paradigm being employed. For the present case, this relation takes the form
where τin denotes the timescales associated with non-reactive inelastic collisions, τreact denotes reaction timescales, and θ corresponds to flow timescales. Thus, elastic processes occurring on timescale τel are assumed to be rapid, and the translational energy mode can, therefore, be described by an equilibrium, or near equilibrium, velocity distribution function governed by temperature T. All inelastic collisions and chemical reactions proceed on slower timescales, and the internal energy distribution can, therefore, take any arbitrary form, as dictated by competing thermochemical processes in the system. Following the method of rapid and slow processes based on the GCE theory, the collision operator of the Boltzmann equation for each pseudo species, sij, is separated as
where s denotes the chemical species, and i and j denote the vibrational and rotational quantum number of rovibrational level, l. Details on the collision operators, Jsij, for the rovibrationally resolved system can be found in the Appendix Eqs. (A1) and (A2) (A = Appendix). Next, a weighting parameter, ε, is introduced in the Boltzmann equation for the rovibrationally resolved system as
where a is the chemical species, and represent the vibrational and rotational quantum numbers, respectively, corresponding to rovibrational level κ of the colliding partner. The quantities fsij and usij denote the rovibrational distribution function, and thermal velocity, respectively. The quantity ε is defined as the ratio , and through it, the timescale variability is introduced into the Boltzmann equation. At this point, it is useful to introduce the equation of transfer within the GCE framework, with the integrals evaluated only for rapid processes,
where ψsij denotes a molecular property of species sij. An additional summation of (4) over all species sij results in the equation of transfer over both rapid and slow processes.
Now, the distribution function in (3) is expanded in a power series in ε up to first order and can also be expressed in terms of a perturbation , acting on the zeroth-order distribution as
Equating coefficients with the lowest like powers of ε, and considering the quantities that are invariant over just the rapid collisions, namely, 1 (i.e., a constant), and the translational energy, , gives the zeroth-order or equilibrium distribution,5,10
where are peculiar velocity components. Thus, it can be seen that the zeroth-order distribution for each species is the equilibrium Maxwellian distribution, and the macroscopic quantities, solely defined by , take the form
In the above equations, denotes the rovibrational energy, and are the sets of chemical species and “pseudo-species,” or rovibrational levels in species s, present in the gas mixture, respectively.
To obtain the first-order flow governing equation for a rovibrationally resolved system, the collisional invariants over rapid processes: 1 (i.e., a constant), and invariants over both rapid and slow processes, namely, momentum, and total energy, are substituted for ψsij in (4), along with the expansion for fsij [(5)]. For the momentum and total energy invariants, (4) is also summed over all sij. Thus, the rovibrationally resolved Navier–Stokes equations take the form
It can be seen that a separate conservation equation for the number density of each species needs to be solved at every time step and spatial location of the flowfield. The source term in (8) denotes the change in the population of rovibrational levels due to inelastic and reactive collisions, and arises due to the collision operator over slow processes, with the distribution function expanded up to first order as per (5)
The set of conservation equations in (8)–(10) also contains the rovibrational transport fluxes up to first order, namely, the diffusion flux with diffusion velocity, defined for each species sij, the stress tensor , and the total energy flux q. The first-order component of these transport fluxes in terms of can be defined as
An additional summation over species, sij has been carried out in the case of the stress tensor and heat flux vector. However, these rovibrational fluxes depend on the first-order perturbation, , which is yet to be determined.
where . The form of the linear integral operator, Isijacd is provided in Eqs. (A7) and (A8). By employing the relation for shown in (6), and replacing all temporal derivatives of macroscopic quantities in terms of their spatial gradients using the zeroth-order governing equations [i.e., Euler equations, see Eqs. (A3)–(A6)], (15) becomes
where is the translational specific heat capacity, and represents the production term from the zeroth-order component of the distribution, and v denotes the hydrodynamic velocity. The quantity on the right side of (16) is the diffusion driving force for species sij given by
where is the gas pressure.
Based on (16), the general solution for the perturbation, , can be taken as
The generalized functions and are vector functions, and is a tensor function of the peculiar velocity, , and is a scalar. Substituting for from (18) in the integral operator, Isijacd, in (16) and equating coefficients of terms having gradients of the same macroscopic parameters produce a set of integral relations provided in Eqs. (A9)–(A12). These integral relations shall be used in evaluating the set of StS rovibrational transport coefficients.
where is defined as the rate of shear tensor given by
The quantities are the various StS transport coefficients, namely, the mass and thermal diffusion coefficients, viscosity, and translational thermal conductivity, respectively. These transport coefficients are related to the generalized coefficients in (18) via bracket integrals, which are defined for the rovibrational system in Appendix-Bracket Integrals. It is clear from the relations in Appendix-Bracket Integrals that the bracket integrals are to be evaluated only over rapid processes, which in this case are just elastic collisions. Additional constraints that guarantee uniqueness (detailed in Appendix-Uniqueness constraints) are also employed, yielding the following definitions of the transport coefficients:
It should be noted that due to the uniqueness constraints given in Eq. (A17), the generalized function and hence, relaxation pressure does not contribute to the stress tensor components.
III. TRANSPORT FORMULATION
Since no closed-form expression is available for the generalized coefficients, we introduce trial functions , , and tensor in place of the generalized functions , and , respectively. These trial functions are then recast as an expansion in Sonine polynomials, in normalized peculiar velocity as shown below:
where r is the number of terms in the expansion, represent the expansion coefficients associated with the corresponding trial functions, and . Replacing the generalized functions in (23) with the respective trial functions provided in (24) and retaining only the lowest non-vanishing expansion terms,10 the StS transport coefficients become
Thus, the problem of estimating the state-based transport coefficients reduces to that of computing the Sonine polynomial expansion coefficients.
The method illustrated in Nagnibeda and Kustova10 to obtain these Sonine polynomial expansion coefficients for a vibrationally resolved system has been adopted in the present work for the rovibrational StS system. To estimate and , the trial functions and are substituted in integral equations Eqs. (A9)–(A12), respectively. The resulting expressions are multiplied on both sides by (27) as follows:
and integrated over rapid processes (i.e., over velocity space). Retaining only the lowest non-zero terms in the expansion in (24), the following set of linear algebraic equations are obtained to compute the expansion coefficients :
It should be noted that the additional equation in (28) arises based on the constraints presented in Appendix-Uniqueness constraints. Similarly, the linear equation set and the additional constraint equation to estimate is given by
The term that appears in (28) and (29) is nothing but the bracket integral of evaluated over rapid processes, and in terms of partial brackets defined in Appendix-Bracket integrals, can be expressed as
A similar approach is adopted to compute the Sonine polynomial expansion coefficients corresponding to shear viscosity, namely, . The trial function, , from (24) is substituted into the integral equation Eq. (A11) and then multiplied on both sides by ,
and the resulting expression is integrated over velocity space. Thus, the linear equation set that gives is
where the bracket integral of is given in terms of partial brackets as
In order to evaluate the Sonine polynomial expansion coefficients using (28), (29), and (32), the quantities and must first be computed. However, the definition of bracket integrals for the rovibrationally resolved system provided in Appendix-Bracket integrals reveals that these quantities (evaluated over rapid processes only) are eight-fold integrals over the thermal velocity space of the colliding particles as well as parameters that influence the collision process, i.e., the interaction potential, impact parameter, etc. For ease of evaluation, the partial brackets (and thereby, the bracket integrals) are systematically reduced5 and expressed as a linear combination of collision integrals, , defined as
where represents the dimensionless relative translational energy between the colliding particles, denotes the differential cross section for an elastic collision, is the solid angle into which the relative velocity is scattered after the collision, and the reduced mass is . For the lowest (first) order approximation of transport coefficients, the necessary relation between the various partial brackets and collision integrals is given in Eqs. (A19) and (A20). Estimating these rovibrationally resolved transport collision integrals, , is the key step in calculating the StS transport coefficients.46–50
IV. CONTINUUM BREAKDOWN PARAMETERS
The limits of the rovibrational Navier–Stokes equations and thereby the onset of continuum breakdown can be identified through quantities known as continuum breakdown parameters.51–56 When the value of the breakdown parameter exceeds a determined threshold, the perturbation of the zeroth-order distribution becomes large, and the continuum approximation underlying the Navier–Stokes equations fails. Here, forms of the “mechanism-based” as well as the “species perturbation” breakdown parameters are derived.
These mechanism-based GCE breakdown parameters are the coefficients of the first-order perturbation , when the perturbation is expressed in terms of the various transport fluxes.52,54–56 Only the lowest non-zero terms are retained within the trial function expansions in (24), and the expansion coefficients are replaced in terms of the respective transport coefficients using the relation in (26),
Recognizing that the coefficients of the normalized peculiar velocity, , are the first-order component of the rovibrational transport fluxes [(19)–(21)], expression for in terms of fluxes [with the superscript dropped hence forth] becomes
where is the inverse most probable speed, and the diffusion flux, , is related to the diffusion velocity, as
Now, it can be seen that the coefficients of in (37) are none other than dimensionless forms of the various transport fluxes. Thus, the mechanism-based GCE continuum breakdown parameters for the rovibrationlly resolved system due to mass and thermal diffusion fluxes, translational heat flux, and components of the stress tensor are given by
The mechanism-based GCE breakdown parameters presented above treat each transport process separately in determining the onset of continuum breakdown. However, all transport mechanisms simultaneously perturb the underlying VDF, leading to continuum breakdown. The species perturbation parameter, originally proposed by Tiwari for a simple gas,57 captures this combined effect. For the rovibrationally resolved system, the species perturbation parameter is evaluated as the 2-norm in Hilbert space of the first-order perturbation, ,
This work presented a rigorous derivation of the rovibrationally resolved Navier–Stokes governing equations. This offers the first complete set of state-resolved hydrodynamic equations with transport flux closure that is consistent with rovibrationally resolved chemical kinetics used in modern computational fluid dynamics techniques. The solution of these governing equations, along with the mechanism-based generalized Chapman–Enskog breakdown parameters, and the species perturbation parameter proposed in this work require evaluation of the rovibrationally resolved transport coefficients. The quantity that incorporates detailed collisional information into the transport coefficients is the collision integral, . Thus, establishing rigorous approaches to evaluate the rovibrational collision integrals through ab initio interaction potentials is of utmost importance in ensuring the accuracy of any transport formulation and continuum breakdown prediction.
This work was supported by the Presidential Early Career Award for Scientists and Engineers (PECASE) Award from the NASA Space Technology Research Grants Program and by the Air Force Office of Scientific Research under Award No. FA9550-17-1-0127. S.S. gratefully acknowledges the support from the Zonta International Amelia Earhart Fellowship Program.
The data that support the findings of this study are available from the corresponding author upon reasonable request.
APPENDIX: Governing Equations Formulations
Details pertaining to the development of the governing equations and transport coefficients within the rovibrationally resolved framework are provided below.
1. Collision Operators
The form of the collision operator over rapid processes (i.e., elastic collisions) appearing in the Boltzmann equation is given by
where denotes the VDF after collision; denote the relative velocity between the colliding species and the differential cross section for an elastic collision between species s and a, with vibrational levels i and c, and rotational levels, j and d, respectively, and represents the solid angle into which the relative velocity is directed after the collision.
The collision operator over slow processes can be separated into collisions involving excitation/de-excitation (E), exchange reactions (EX), and dissociation/recombination (D) as shown below:
The form of each of these operators can be found in Nagnibeda and Kustova.10
2. Euler equations: Zeroth-order governing equations
The macroscopic hydrodynamic equations to zeroth-order are given by
Here, represents the source term evaluated to zeroth-order and arises due to the collision operator over slow processes,
3. Integral relations
The term containing the linear integral operator, Isijacd in Eq. (15) can be related to the collision operator over rapid processes, . For an expansion of the distribution function up to first-order, this relation becomes
Using and the form of the collision operator, provided in (A1), the integral operator Isijacd can be written as
a. Integral relations of the generalized functions
The general solution for shown in Eq. (18) is substituted in the linear integral operator, , on the left-hand side of Eq. (16). When coefficients of terms having like gradients in macroscopic quantities are equated, the following set of integral relations are obtained:
The above relations are evaluated for all: .
4. Bracket integrals
The bracket integral that appears in the derivation of state-based transport coefficients is a bilinear operator and can be expressed, in general, as
where both X and Y are scalars, vectors, or tensors. The partial brackets take the form
where and are post-collision values. It is clear from the above equations that bracket integrals are evaluated only over rapid collisions, which in the present case includes just elastic collisions. The symmetry relations for the brackets defined above are given by
5. Uniqueness constraints
To guarantee uniqueness, the constraints imposed within the GCE method is that macroscopic quantities are solely defined by the zeroth-order component of the distribution function. Thus, moments of the higher-order components of the distribution function, in this case the first-order, are set to zero as shown below:
where the generalized function, , is a scalar function of the peculiar velocity Csij, whereas , and can be expressed as
6. Relation between partial brackets and collision integrals
The various StS transport coefficients depend of the values of the Sonine polynomial expansion coefficients [Eq. (26)], which in turn can be evaluated only if quantities and are known. Now, and , defined as per Eqs. (30) and (33), contain partial brackets of and , respectively. It is clear from (A14) that these partial brackets are eight-fold integrals requiring details of the exact collision process. Hence, for ease of evaluation, these partial brackets are expressed as linear combination collision integrals. Below, the relation between partial brackets of and and collision integrals for a first-order evaluation of StS transport coefficients is provided below: