Navier–Stokes equations with rovibrational state-resolved transport flux closure

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.¹ 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.¹¹ More formal kinetic theory-based arguments proposed by Wang Chang and Uhlenbeck^7,12 and Mason and Monchick^13,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 equations^15–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 models^18–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.²³ 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) framework^10,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, J_sij, 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 $c, d$ represent the vibrational and rotational quantum numbers, respectively, corresponding to rovibrational level κ of the colliding partner. The quantities f_sij and u_sij denote the rovibrational distribution function, and thermal velocity, respectively. The quantity ε is defined as the ratio $τel/θ$⁠, 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 $ϕsijGCE$⁠, 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, $msusij2/2$⁠, gives the zeroth-order or equilibrium distribution,^5,10

where $C1, C2, C3$ 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 $fsij(0)$⁠, take the form

In the above equations, $εsij$ denotes the rovibrational energy, $N$ and $Ls$ 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, $msus$ and total energy, $msus2/2+εsij$ are substituted for ψ_sij in (4), along with the expansion for f_sij [(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 $Rsij(1)$ 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, $Jsijslow(fsij,facd)$ 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, $Vsij$ defined for each species sij, the stress tensor $T$⁠, and the total energy flux q. The first-order component of these transport fluxes in terms of $fsij(1)$ 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, $ϕsijGCE$⁠, which is yet to be determined.

The first-order perturbation, $ϕsijGCE$⁠, is estimated through an integral equation, which is obtained by equating coefficients of the next highest power of ε in the power series expansion of (3) using (5)

where $ϕsijGCE=fsij(1)/fsij(0)$⁠. The form of the linear integral operator, I_sijacd is provided in Eqs. (A7) and (A8). By employing the relation for $fsij(0)$ 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 $cv,tr=3nkB/2ρ$ is the translational specific heat capacity, and $Rsij(0)$ represents the production term from the zeroth-order component of the distribution, and v denotes the hydrodynamic velocity. The quantity $dsij$ on the right side of (16) is the diffusion driving force for species sij given by

where $p$ is the gas pressure.

Based on (16), the general solution for the perturbation, $ϕsijGCE$⁠, can be taken as

The generalized functions $Asij$ and $Dsijacd$ are vector functions, and $Bsij$ is a tensor function of the peculiar velocity, $Csij$⁠, and $Gsij$ is a scalar. Substituting for $ϕsijGCE$ from (18) in the integral operator, I_sijacd, 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.

The first-order components of the transport fluxes in (12)–(14) can be recast by using the expression for $fsij(0)$ from (6) and $ϕsijGCE$ from (18) as

where $S$ is defined as the rate of shear tensor given by

The quantities $Dsijacd, DT,sij, μsij,and ktr,sij$ 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 $Gsij=0$ and hence, relaxation pressure does not contribute to the stress tensor components.

At this point, the unknowns are the various generalized functions. Evaluating these quantities and their bracket integrals will enable the calculation of the various StS transport coefficients through (23). The approach used to estimate these generalized functions is described in Sec. III.

Since no closed-form expression is available for the generalized coefficients, we introduce trial functions $Asij$⁠, $Dsijacd$⁠, and tensor $B¯¯sij$ in place of the generalized functions $Asij, Dsijacd$⁠, and $Bsij$⁠, respectively. These trial functions are then recast as an expansion in Sonine polynomials, $Sν(n)(x)$ in normalized peculiar velocity $Csij2$ as shown below:

with

where r is the number of terms in the expansion, $Asij,r, Dsij,racd, Bsij,r$ represent the expansion coefficients associated with the corresponding trial functions, and $Csij=mCsij/2kBT$⁠. Replacing the generalized functions in (23) with the respective trial functions provided in (24) and retaining only the lowest non-vanishing expansion terms,¹⁰ 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 Kustova¹⁰ 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 $Asij,r$ and $Dsij,racd$⁠, the trial functions $Asij$ and $Dsijacd$ 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 $Asij,r$⁠:

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 $Dsij,racd$ is given by

The term $Λsijacd$ that appears in (28) and (29) is nothing but the bracket integral of $Qsij(r)$ 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, $Bsij,r$⁠. The trial function, $B¯¯sij$⁠, from (24) is substituted into the integral equation Eq. (A11) and then multiplied on both sides by $Qsij(r)$⁠,

and the resulting expression is integrated over velocity space. Thus, the linear equation set that gives $Bsij,r$ is

where the bracket integral of $Qsij(r)$ is given in terms of partial brackets as

In order to evaluate the Sonine polynomial expansion coefficients using (28), (29), and (32), the quantities $Λrr′sijacd$ and $Hrr′sijacd$ 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 reduced⁵ and expressed as a linear combination of collision integrals, $Ωsijacd(w,x)$⁠, defined as

where $γ2=msag2/2kBT=εtr/kBT$ represents the dimensionless relative translational energy between the colliding particles, $σsijacdel$ denotes the differential cross section for an elastic collision, $dΩ$ is the solid angle into which the relative velocity is scattered after the collision, and the reduced mass is $msa=(msma)/(ms+ma)$⁠. 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, $Ωsijacd(w,x)$⁠, is the key step in calculating the StS transport coefficients.^46–50 

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 $ϕsijGCE$⁠, 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 $Asij,0, Asij,1, Dsij,0acd, Bsij,0$ are replaced in terms of the respective transport coefficients using the relation in (26),

Next, the generalized functions $Asij, Dsijacd$⁠, and $Bsij$ in (18) are replaced by the above (35) to obtain

Recognizing that the coefficients of the normalized peculiar velocity, $Csij$⁠, are the first-order component of the rovibrational transport fluxes [(19)–(21)], expression for $ϕsijGCE$ in terms of fluxes [with the superscript $(1)$ dropped hence forth] becomes

where $βs=ms/2kBT$ is the inverse most probable speed, and the diffusion flux, $Jsij$⁠, is related to the diffusion velocity, $Vsij(1)$ as

Now, it can be seen that the coefficients of $Csij$ 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,⁵⁷ 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, $ϕsijGCE$⁠,

By substituting the expressions for $fsij(0)(Csij)$ and $ϕsijGCE(Csij)$ from (6) and (37), the rovibrationally resolved, species perturbation parameter takes the form

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, $Ωsijacd(w,x)$⁠. 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.

Details pertaining to the development of the governing equations and transport coefficients within the rovibrationally resolved framework are provided below.

The form of the collision operator over rapid processes (i.e., elastic collisions) appearing in the Boltzmann equation is given by

where $f′$ denotes the VDF after collision; $g, σsijacdel$ 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 $dΩ$ 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.¹⁰ 

The macroscopic hydrodynamic equations to zeroth-order are given by

Here, $Rsij(0)$ represents the source term evaluated to zeroth-order and arises due to the collision operator over slow processes,

The term containing the linear integral operator, I_sijacd in Eq. (15) can be related to the collision operator over rapid processes, $Jsijrap(fsij,facd)$⁠. For an expansion of the distribution function up to first-order, this relation becomes

Using $f(1)=ϕGCEf(0)$ and the form of the collision operator, $Jsijrap$ provided in (A1), the integral operator I_sijacd can be written as

The general solution for $ϕsijGCE$ shown in Eq. (18) is substituted in the linear integral operator, $Isijacd(ϕGCE)$⁠, 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: $s, h∈N, (i, j)≡l∈Ls, (y, z)≡ξ∈Lh$⁠.

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 $X′$ and $Y′$ 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

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:

When $fsij(1)$ in the above equation is replaced as a product of $ϕsijGCE$ and $fsij(0)$ with the form for $ϕsijGCE$ specified in Eqs. (18), (A16) becomes

where the generalized function, $Gsij$⁠, is a scalar function of the peculiar velocity C_sij, whereas $Asij, Dsijacd$⁠, and $Bsij$ can be expressed as