A constitutive relation generalizing the Navier–Stokes theory to high-rate regimes

Abstract

We propose a constitutive equation for flows of gases in high-rate regimes where the Navier–Stokes theory breaks down. The model generalizes the Navier–Stokes relation and agrees well with that model in all lower rate flows examined. Our proposed constitutive relation is calibrated with the method of objective molecular dynamics (OMD) using families of compressible and incompressible flows of Lennard-Jones argon. The constitutive relation makes use of the higher-order objective strain rates due to Rivlin and Ericksen (J. Rat. Mech. Anal., vol. 4, 1955, pp. 323–425). The constitutive relation is fully frame-indifferent, and the macroscopic flows corresponding to the OMD simulations are exact solutions for the proposed model. The model is shown to agree with atomistic results much better than the Navier–Stokes equations in the transition regime. The success of our model indicates that it is not higher gradients that become important in the high-rate regime, but rather higher rates of change of the strain rate tensor. While somewhat more complicated to implement than the Navier–Stokes relation, the proposed model is expected to be compatible with existing methods of computational fluid dynamics and may extend those methods to higher rate regimes, while preserving their ability to handle large spatial scales.

1. Introduction

Hypersonic vehicles during re-entry encounter varying flow regimes during their descent, varying from free-molecular to transition to continuum. The ability to make better predictions of the flows in the entire range determines the feasibility of space missions and development of future hypersonic vehicles. Computational modelling plays an indispensable role due to limitations of the replication of the extreme environment encountered in such missions in ground-based experiment facilities. Computational fluid dynamics (CFD) modelling, which is based on the numerical solution of the partial differential equations of macroscopic conservation equations, works well in the continuum regime where the collision rate is high enough such that near-equilibrium molecular distribution functions are maintained. The CFD methods are less successful in the continuum-transition regime due to the breakdown of the underlying Navier–Stokes–Fourier (NSF) constitutive model. The approximations made there for the stress and heat flux, though applicable in the near-equilibrium regime, fall short in the transition regime.

A viable technique to be used in this regime is Bird's stochastic particle-based method of direct simulation Monte Carlo (DSMC) (Bird 1994; Boyd & Schwartzentruber 2017). The DSMC statistically solves Boltzmann transport equation (BTE) which defines the evolution of the phase density (velocity distribution function (VDF)). It is believed that it can describe dilute flows in all the regimes. Though accurate, DSMC is computationally expensive as a large number of computational cells and molecules, and small time-steps, are required to model near-continuum behaviour. Hence, there is a need to develop a better set of constitutive relations to be incorporated in CFD modelling. In this work we focus on that goal by building a higher-order constitutive model calibrated by a special family of flows in the continuum-transition regime.

Alternative approaches are discussed in detail by Garzó & Santos (2003). One can take the kinetic viewpoint and work with the BTE. However, its closed form analytical solutions are difficult to obtain except for some simple flows. It can be solved approximately via the Chapman–Enskog (CE) perturbative method where the velocity distribution function is assumed to be a perturbation expansion in the Knudsen number (${Kn}$) about the equilibrium Maxwell–Boltzmann distribution. Retaining terms up to second and third order yields the Burnett and super-Burnett equations, respectively (Cercignani 1975; Chapman & Cowling 1990). The Burnett equations were studied extensively but are known to violate frame-indifference and the second law of thermodynamics at high $Kn$ (Müller 1972; Comeaux, Chapman & MacCormack 1995). Moreover, although the CE expansion is believed to be asymptotic, there is no reasoning that it is convergent (Grad 1963). These issues raise questions over possible improvements of the Navier–Stokes constitutive law by including higher-order terms in the CE expansion. In recent years, several authors have proposed augmented forms of the Burnett equation which contains additional terms from super-Burnett order to stabilize the Burnett equations or Bhatnagar–Gross–Krook (BGK)–Burnett equations (Zhong, MacCormack & Chapman 1993; Balakrishnan, Agarwal & Yun 1999).

Another approximate approach is the moment method. In this method, moments of the BTE are taken and evolution equations of the higher-order moments are derived. The number of moments needed to describe a particular process is unsettled but is known to increase with Knudsen number. For some strongly non-equilibrium system, one may need hundreds of moment equations, which soon become intractable (Weiss & Müller 1995; Müller 2008). Furthermore, each moment depends on the divergence of the next higher moment, so closure conditions are needed. Classically, this difficulty is addressed by using physically motivated closure conditions (Grad 1949, 1958; Müller & Ruggeri 2013); Grad's 13 moment method, Grad's 26 moment method (Grad 1949; Levermore 1996). Some works that develop this approach include Eu's evolution equations (Eu 1992) leading to a nonlinear constitutive coupled relation, further developed by Myong (1999), and Singh's moment equations based on a distribution function consistent with Onsager's principle (Singh & Agrawal 2016) leading to Burnett-type stable constitutive relations (Singh, Jadhav & Agrawal 2017). In spite of extensive work in this direction, there does not exist a simple non-classical, well accepted constitutive law like the Navier–Stokes relation for high gradient/rarefied flows.

Significant development of methods to obtain constitutive equations or related macroscopic models comes from the community of researchers studying granular flows (Saha & Alam 2016, 2017; Rongali & Alam 2018a,b). The study of normal stress effects and the derivation of constitutive relations is in the spirit of the present work.

In this work, we take an alternative route and approach the problem using patterns of thought of classical constitutive theory. We examine theories proposed by Reiner, Rivlin & Ericksen, et al. (Reiner 1945; Rivlin & Ericksen 1955; Noll 1974; Rivlin 1997) using the deterministic method of objective molecular dynamics (OMD) (Dumitrică & James 2007; Dayal & James 2010, 2012; Pahlani et al. 2021). In contrast to the kinetic theory approach, transport coefficients do not emerge naturally from the theory itself. We propose and calibrate a frame-indifferent constitutive relation using various flows of Lennard-Jones (LJ) argon gas. The model is termed the Rivlin–Ericksen (RE) constitutive relation. The present work builds on our previous study where a related non-classical constitutive equation was proposed, but limited to the special case of uniform simple-shear (Pahlani, Schwartzentruber & James 2022).

The OMD is derived from the invariance of the equations of molecular dynamics. It provides a (time-dependent) invariant manifold of these equations. It allows for simulation of macroscopic velocity fields of the form

(1.1)\begin{equation} {\boldsymbol{v}}=\boldsymbol{\mathsf{A}}(\boldsymbol{\mathsf{I}}+t\boldsymbol{\mathsf{A}})^{{-}1}{\boldsymbol{x}}, \end{equation}

where $\boldsymbol{\mathsf{A}}$ is an arbitrary assigned $3\times 3$ matrix. This velocity field includes numerous examples of steady and unsteady compressible and incompressible flows. We term this family of OMD flows as ‘universal flows’ (Dayal & James 2012) because the equations of motion of continuum mechanics (no body force) are identically satisfied by (1.1) for every choice of $\boldsymbol{\mathsf{A}}$ and every accepted material model (fluids or solids). In addition, all OMD simulations have a particular statistics that is represented by an ansatz for the molecular density function of the form $f(t, {\boldsymbol {x}}, {\boldsymbol {v}}) = g(t, {\boldsymbol {v}} - \boldsymbol{\mathsf{A}}(\boldsymbol{\mathsf{I}} + t \boldsymbol{\mathsf{A}})^{-1}{\boldsymbol {x}})$. Substitution of this ansatz into the Boltzmann equation gives an exact reduction, see Dayal & James (2012, § 6.1). At atomic level there are no restrictions on the atomic forces or the species of atom. This synergy of atomistic and continuum-level theories is the main motivation behind choosing the method of OMD to test higher-order non-classical constitutive models. An added advantage is that molecular dynamics (MD) equations are only solved for a finite set of atoms denoted by ${\boldsymbol {y}}_k(t),\ k=1,2,\ldots,M$ called ‘simulated atoms’. All the other infinite atoms (the non-simulated atoms) of the system are obtained by applying time-dependent translation group

(1.2)\begin{equation} G_T=\{(\boldsymbol{\mathsf{I}}\,|\,\nu^1\hat{\boldsymbol{ e}}_1+\nu^2\hat{\boldsymbol{ e}}_2+\nu^3\hat{\boldsymbol{ e}}_3): \nu = (\nu^1, \nu^2, \nu^3) \in \mathbb{Z}^3\}, \end{equation}

to the set of simulated atoms. The relationship between the positions and velocities of the simulated and non-simulated atoms given by

(1.3)\begin{equation} \left.\begin{gathered} {\boldsymbol{y}}_{\nu,k}(t) = g_{\nu}({\boldsymbol{y}}_{k}(t)),\quad {\boldsymbol{y}}_{\nu,k}(t) = {\boldsymbol{y}}_{k}(t) + \sum_{i = 1}^3 \nu^i (\boldsymbol{\mathsf{I}} + t \boldsymbol{\mathsf{A}}) {\boldsymbol{e}}_i,\quad g_{\nu} \in G_T,\\ {\boldsymbol{v}}_{\nu,k}(t) = {\boldsymbol{v}}_{k}(t) + \sum_{i = 1}^3 \nu^i \boldsymbol{\mathsf{A}} {\boldsymbol{e}}_i, \end{gathered}\right\} \end{equation}

where $\{{\boldsymbol {y}}_{\nu,k}(t),{\boldsymbol {v}}_{\nu,k}(t)\}$ and $\{{\boldsymbol {y}}_{k}(t),{\boldsymbol {v}}_{k}(t)\}$ are the trajectory of non-simulated and simulated atoms, respectively. Here $\boldsymbol{\mathsf{I}}$ is an identity matrix, $t$ is time and $\boldsymbol{\mathsf{A}}$ is an arbitrary assigned constant second-order tensor that guides the macroscopic flow. This reduces the complicated MD equations to a system of non-autonomous ordinary differential equations for the motions of the simulated atoms in standard form. The main theorem of OMD is that the non-simulated atoms also satisfy the MD equations exactly, even though only the simulated atoms are subjected to MD equations and all the other atoms are obtained by an explicit formula, (1.3). Thus, OMD is an exact atomistic analogue of the flow fields of the form (1.1). This relies on the well-known symmetries of quantum mechanics: the frame-indifference and permutation invariance of atomic forces. The main proof of the theorem and computational details to perform the simulation is detailed elsewhere (Dumitrică & James 2007; Dayal & James 2010; Pahlani et al. 2021, 2022; Pahlani, Schwartzentruber & James 2023), and the connections with the Boltzmann equation can be found in Dayal & James (2010) and James, Nota & Velázquez (2019).

The paper is organized as follows. In § 2 we compare OMD simulations with fluid dynamics. In § 3 we define and calibrate the RE constitutive model and show that it improves the NSF constitutive relation for various incompressible, compressible and unsteady flows. In § 4 we discuss the thermodynamics and stability of the proposed RE model. In § 5 we explain a connection between the RE and Burnett equations. The connection with the exact moment method for Maxwellian gas (molecules interacting by an inverse fifth-power law of force) under simple shear motion was established in our previous work (Pahlani et al. 2022). Finally the conclusions are contained in § 6.

General background for this work can be found in Appendices A (numerical method), B (OMD flows that are viscometric flows) and C (the Burnett equation for OMD flows).

2. Comparison between OMD and fluid dynamics with the NSF model

A heat conducting monoatomic ideal compressible gas is described by the following set of balance laws of mass, momentum and energy:

(2.1)\begin{equation} \left.\begin{gathered} \rho_{t}+\textrm{div}({\rho}{\boldsymbol{v}})=0,\\ \rho({\boldsymbol{v}}_t+\boldsymbol{\nabla}{\boldsymbol{v}}{\boldsymbol{v}})=\textrm{div}\ {{\boldsymbol\sigma}},\\ {\rho}(e_{t}+\boldsymbol{\nabla}e\boldsymbol{\cdot}{\boldsymbol{v}})={{\boldsymbol\sigma}}: \boldsymbol{\nabla}{\boldsymbol{v}}-\textrm{div}\ {{\boldsymbol{q}}}, \end{gathered}\right\} \end{equation}

where $\rho _t$ and $v_t$ are partial derivatives of density and velocity field with respect to time, respectively. The OMD velocity field, (1.1), satisfies the momentum conservation equation identically and OMD simulations generate time-dependent stress and temperature fields which are spatially homogeneous. This reduction substantially simplifies the coupled partial differential equations (2.1) to a system of ordinary differential equations in density $\rho (t)$ and temperature $T(t)$, when substituted with the Navier–Stokes constitutive relation with the ideal gas law, i.e.

(2.2)\begin{equation} \left.\begin{gathered} {\boldsymbol\tau}=\mu_{NSF}\left(\boldsymbol{\nabla}{\boldsymbol{v}}+\boldsymbol{\nabla}{\boldsymbol{v}}^{\rm T}- \frac{2}{3}(\boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{v}})\boldsymbol{\mathsf{I}}\right),\quad e=c_{v}T,\\ \frac{{\rm d}\rho}{{\rm d} t}={-}\rho {\textrm{tr}}\ {\boldsymbol{\mathsf{A}}}({\boldsymbol{I}}+t{\boldsymbol{\mathsf{A}}})^{{-}1},\quad \rho(0)=\rho_0, \end{gathered}\right\} \end{equation}

(2.3)\begin{align} \left.\begin{aligned} \frac{{\rm d} T}{{\rm d} t}&=\frac{-R}{c_{v}}\textrm{tr}({\boldsymbol{\mathsf{A}}}({\boldsymbol{\mathsf{I}}}+t{\boldsymbol{\mathsf{A}}})^{{-}1})T(t)+ \frac{\mu_{NSF}(T(t))}{\rho_{0}c_{v}}(|{\boldsymbol{\mathsf{A}}}({\boldsymbol{\mathsf{I}}}+t{\boldsymbol{\mathsf{A}}})^{{-}1}|^2\\ &\quad +\textrm{tr}(({\boldsymbol{\mathsf{A}}}({\boldsymbol{\mathsf{I}}}+t{\boldsymbol{\mathsf{A}}})^{{-}1})^2) -\frac{2}{3}{\textrm{tr}}({\boldsymbol{\mathsf{A}}}({{\boldsymbol{\mathsf{I}}}}+t{{\boldsymbol{\mathsf{A}}}})^{{-}1}))^2 \\ &\quad\times\exp\left(\left(-\int_0^t \textrm{tr}\ {\boldsymbol{\mathsf{A}}}({\boldsymbol{\mathsf{I}}}+s{\boldsymbol{\mathsf{A}}})^{{-}1}\ {{\rm d}}s\right)\right),\\ T(0)=T_{0}, \end{aligned}\right\} \end{align}

where $R$ is the specific gas constant and $c_v$ is the specific heat at constant volume. This is solved using Runge–Kutta approach for a given viscosity model $\mu _{NSF}(T)$ (Kim & Monroe 2014) which needs to be consistent with the OMD force field for direct comparison.

For OMD simulations, we consider argon gas of initial density $1.25\ \textrm {kg}\ \textrm {m}^{-3}$ and initial temperature $T=400\ \textrm {K}$. The system is initialized by first choosing a set of simulated argon atoms. These can lie anywhere in space but, to reduce computational complexity, they are assumed to lie initially in a fundamental domain, with number density dependent on the initial density of the gas. Within the fundamental domain, the positions are generated randomly and their velocities are sampled from the Maxwell–Boltzmann distribution at preassigned initial temperature $T$. Once the simulation starts, the simulated atoms quickly diffuse into the non-simulated atoms. The complete computational details are provided in the previous study (Pahlani et al. 2023). The LJ potential was utilized to model interaction between argon atoms. Only interactions which are within a certain cutoff are computed since the LJ potential tends to zero rapidly at large distances. The potential is given by

(2.4)\begin{equation} \phi(r_{ij})=4\epsilon_{LJ} \left[\left(\frac{\sigma_{LJ}}{r_{ij}}\right)^{12}- \left(\frac{\sigma_{LJ}}{r_{ij}}\right)^{6}\right], \end{equation}

where $r_{ij}$ is the distance between atoms, $\epsilon _{LJ}= 1.65 \times 10^{-21}$ J and $\sigma _{LJ}=3.4 \times 10^{-10}$ m. Note that the forces acting on simulated atoms are a function of each and every atom of the system. The second-order accurate velocity Verlet scheme is used to integrate the equations of motion. Evolution of OMD trajectories can be used to define time-dependent macroscopic temperature and stress. Temperature is defined by the variance of the kinetic energy and the Boltzmann definition is used for the stress, where there is no contribution to momentum flux from collisions, due to the operational regime of a dilute gas. These are given by (Boyd & Schwartzentruber 2017)

(2.5)\begin{equation} \left.\begin{gathered} T=\left\langle\frac{m}{3k_b N}\left[\sum_{i}^{N}(|{\boldsymbol{v}}^{'}_{i}|^2) \right]\right\rangle,\\ {\boldsymbol\sigma}=\left\langle\frac{\rho}{N}\sum_{i}^{N}\left[{\boldsymbol{v}}^{'}_{i}\otimes{\boldsymbol{v}}^{'}_{i}\right]\right\rangle. \end{gathered}\right\} \end{equation}

In these expressions, $k_b$ denotes Boltzmann constant, $m$ denotes atomic mass, ${\boldsymbol {v}}^{'}_i$ denotes thermal velocity (i.e. the difference between the particle velocity and the mean velocity of flow) of particle $i$, $N$ denotes number of simulated atoms and an angular bracket denotes an ensemble average over multiple simulations initialized at same macroscopic equilibrium conditions with varying random seeds for initial positions and velocities.

In figure 1 we compare the temperature evolution by OMD with that of the NSF theory for an unsteady, compressible velocity field of the form

(2.6)\begin{equation} {{\boldsymbol{v}}}({{\boldsymbol{x}}},t)=\frac{1}{1+({{\boldsymbol{a}}}\boldsymbol{\cdot}{{\boldsymbol{n}}})t}({{\boldsymbol{n}}}\boldsymbol{\cdot}{{\boldsymbol{x}}}){\boldsymbol{a}} . \end{equation}

This is achieved by choosing ${\boldsymbol{\mathsf{A}}}={\boldsymbol {a}} \otimes {\boldsymbol {n}},\ {\boldsymbol {a}} \boldsymbol {\cdot } {\boldsymbol {n}} \ne 0$, where we select a particular case: ${\boldsymbol {a}}=a{\boldsymbol {e}}_1,\ {\boldsymbol {n}}=n_1{\boldsymbol {e}}_1+n_2{\boldsymbol {e}}_2,\ an_{1,2}=\gamma _{1,2}>0$. This flow reveals the competition between two effects: dilatation and shear. When dilatation dominates the thermal energy of the gas, the temperature decreases; and when the shear effect overpowers, the temperature of the system increases, as expected. (The ‘pressure-shear viscometer’ (Dayal & James 2012) is intended to produce such flows.) Different values of $\gamma _{1,2}$ are chosen for comparison with the NSF equation. For $\gamma _1=\gamma _2 = 2.304 \times 10^8$ s$^{-1}$, the temperature of the system drops and NSF and OMD agree very well. For $\gamma _1=\gamma _2 = 2.30 \times 10^9$ s$^{-1}$, the NSF temperature increases. On the other hand, OMD predicts a continual decrease in temperature as shown in figure 1(b). Under high macroscopic flow gradient and rate, the gas exhibits such strong non-equilibrium behaviour that its correct prediction is beyond the range of applicability of NSF. This qualitative failure of NSF motivates the development of a higher-order theory for highly non-equilibrium flows. Note that the comparison between the two theories is performed once the OMD system is out of the transient regime and gradients are fully developed in the simulation.

Figure 1. Comparison of temperature evolution by OMD with that of the Navier–Stokes theory for an unsteady, compressible flow of pressure shear for LJ argon gas: (a) $\gamma _1=\gamma _2 = 2.304 \times 10^8$ s$^{-1}$; (b) $\gamma _1=\gamma _2 = 2.30 \times 10^9$ s$^{-1}$. Here, red line, OMD; blue solid circle, Navier–Stokes.

3. The RE constitutive model

The RE model (Rivlin & Ericksen 1955), which is more general than the Reiner–Rivlin relation (discussed in Supplementary material in detail, available at https://doi.org/10.1017/jfm.2023.727), assumes the stress to be dependent on the deformation gradient, velocity gradient, acceleration gradient, etc. Intuitively, it can capture the dependence of the state of the fluid on past collisions, not captured by the velocity gradient alone. In general, a fluid where stress depends on only a finite number of these time derivatives is called ‘fluid of differential type’. If the fluid is assumed to be isotropic, possessing material symmetry of unimodular group, $(G=SL(3))$, then the principle of frame indifference leads to the conclusion that stress must depend on the density $\rho$ and time derivatives of the deformation gradient through the objective RE tensors $\boldsymbol{\mathsf{A}}_1, \boldsymbol{\mathsf{A}}_2,\ldots$ defined by

(3.1)\begin{equation} \boldsymbol{\mathsf{A}}_i=\partial^i(\boldsymbol{\mathsf{F}}_t(\tau)^{{\rm T}}\boldsymbol{\mathsf{F}}_t(\tau))/\partial\tau^i|_{\tau=t}, \end{equation}

where $\boldsymbol{\mathsf{F}}_t(\tau )$ is the relative deformation gradient (Coleman, Markovitz & Noll 2012). The RE tensors satisfy the recursive relation

(3.2a,b)\begin{equation} \boldsymbol{\mathsf{A}}_1=\boldsymbol{\nabla} {\boldsymbol{v}}+\boldsymbol{\nabla} {\boldsymbol{v}}^\mathrm{T},\quad \boldsymbol{\mathsf{A}}_{n}={\dot{\boldsymbol{\mathsf{A}}}}_{n-1}+(\boldsymbol{\nabla} {\boldsymbol{v}})^\mathrm{T} \boldsymbol{\mathsf{A}}_{n-1}+\boldsymbol{\mathsf{A}}_{n-1} (\boldsymbol{\nabla} {\boldsymbol{v}}), \end{equation}

where the dot denotes the material time derivative. The physical dimension of $\boldsymbol{\mathsf{A}}_n$ are ${(time)}^{-n}$.

Thus, the stress tensor $\sigma$ is expressible as an isotropic tensor function of the RE tensors $\boldsymbol{\mathsf{A}}_1, \boldsymbol{\mathsf{A}}_2,\ldots, \boldsymbol{\mathsf{A}}_n$, as follows:

(3.3)\begin{equation} {\boldsymbol\sigma}=f(\rho, \boldsymbol{\mathsf{A}}_1, \boldsymbol{\mathsf{A}}_2, \boldsymbol{\mathsf{A}}_3,\ldots). \end{equation}

This reduces to Reiner–Rivlin theory by setting $n=1$. There is also a relation between RE theory and constitutive models of materials with memory (Green & Rivlin 1957) which is given by

(3.4)\begin{equation} {\boldsymbol\sigma}=\boldsymbol{\mathsf{F}}(t)\phi\{\boldsymbol{\mathsf{C}}(s)\}\boldsymbol{\mathsf{F}}^{{\rm T}}(s),\quad -\infty < s \le t, \end{equation}

where $\boldsymbol{\mathsf{C}}(s)=\boldsymbol{\mathsf{F}}_t^\textrm {T}(s)\boldsymbol{\mathsf{F}}_t(s)$ is the right (relative) Cauchy–Green tensor. The constitutive (3.4) accounts for the dependence on the full past history of deformation and temperature. Its connection with the RE model comes from the hypothesis that for sufficiently smooth history of deformation, $\boldsymbol{\mathsf{C}}(s)$ can be expanded as a Taylor series about time $t$ at which the stress is measured, and for the fluid with instantaneous memory which does not exhibit gradual stress relaxation (for example a dilute gas considered in this work) (3.4) can be expressed in terms of the instantaneous values of the deformation gradient and its time derivatives and hence the RE tensors. Thus, for a properly invariant constitutive relation that embodies the dependence of the past history of deformation on the present, the RE model is quite general.

An interesting feature of OMD flows is that these satisfy $\boldsymbol{\mathsf{A}}_i=0$ for order three and higher $(i\ge 3)$, a feature shared by some viscometric flows as well (Coleman et al. 2012). (Generally, OMD flows and viscometric flows are different.) Rivlin used this reduction to find some exact solutions for simple shearing and torsional flow (Rivlin 1956). Viscometric flows are widely used in experimental research to find viscometric functions which can characterize complex fluid behaviour and determine its fundamental properties. Unlike OMD flows, not all viscometric flows are exact solutions of the equations of motion. Also there does not exist an exact atomistic analogue for viscometric flows except simple shearing, which is also an OMD flow (the intersection between OMD and viscometric flows is provided in Appendix B). The reduction of differential type constitutive relations for affine velocity fields was noticed earlier by Bird & Huilgol (1999) but its connection with atomic-level theory dealt in this work was previously unnoticed in the field. For rigorous definitions of viscometric flows, the reader is referred to literature by Coleman et al. (2012). In view of these facts we believe that the OMD flows have been underutilized both theoretically and experimentally.

In this work, the stress tensor is assumed to have a representation of the form

(3.5) \begin{equation} \left.\begin{gathered} {\boldsymbol\sigma}=p\boldsymbol{\mathsf{I}}+\tilde{f}(\rho, \boldsymbol{\mathsf{A}}_1, \boldsymbol{\mathsf{A}}_2), \\ \boldsymbol{\sigma}=p\boldsymbol{\mathsf{I}}+ \mu \boldsymbol{\mathsf{A}}_1 + \alpha_1 \boldsymbol{\mathsf{A}}_2 +\alpha_2 \boldsymbol{\mathsf{A}}_1^2 +\alpha_3 \boldsymbol{\mathsf{A}}_2^2+ \alpha_4 (\boldsymbol{\mathsf{A}}_1 \boldsymbol{\mathsf{A}}_2+ \boldsymbol{\mathsf{A}}_2\boldsymbol{\mathsf{A}}_1) \\ +\,\alpha_5 (\boldsymbol{\mathsf{A}}_1^2 \boldsymbol{\mathsf{A}}_2+ \boldsymbol{\mathsf{A}}_2\boldsymbol{\mathsf{A}}_1^2)+\alpha_6 (\boldsymbol{\mathsf{A}}_1 \boldsymbol{\mathsf{A}}_2^2+ \boldsymbol{\mathsf{A}}_2^2\boldsymbol{\mathsf{A}}_1) \\ +\,\alpha_7 (\boldsymbol{\mathsf{A}}_1^2 \boldsymbol{\mathsf{A}}_2^2+ \boldsymbol{\mathsf{A}}_2^2\boldsymbol{\mathsf{A}}_1^2), \end{gathered}\right\} \end{equation}

where $\alpha _i$ is a scalar function of the invariants of $\boldsymbol{\mathsf{A}}_1$ and $\boldsymbol{\mathsf{A}}_2$. These simplifications reduce the viscous stress tensor to

(3.6)\begin{equation} {\boldsymbol\tau}={-}({\boldsymbol\sigma}-({\rm tr}({{\boldsymbol\sigma}})/3)\boldsymbol{\mathsf{I}})={-}({\boldsymbol\sigma}-p\boldsymbol{\mathsf{I}}) \end{equation}

(3.7)\begin{align} {\boldsymbol\tau} &= \mu\left[\boldsymbol{\mathsf{A}}_1-\frac{\mathrm{tr}(\boldsymbol{\mathsf{A}}_1)}{3}\boldsymbol{\mathsf{I}}\right]+ \alpha_1\left[\boldsymbol{\mathsf{A}}_2-\frac{\mathrm{tr}(\boldsymbol{\mathsf{A}}_2)}{3}\boldsymbol{\mathsf{I}}\right] +\alpha_2 \left[\boldsymbol{\mathsf{A}}_1^2-\frac{\mathrm{tr}(\boldsymbol{\mathsf{A}}_1^2)}{3}\boldsymbol{\mathsf{I}}\right] \nonumber\\ &\quad +\alpha_3 \left[\boldsymbol{\mathsf{A}}_2^2-\frac{\mathrm{tr}(\boldsymbol{\mathsf{A}}_2^2)}{3}\boldsymbol{\mathsf{I}}\right] + \alpha_4\left[(\boldsymbol{\mathsf{A}}_1\boldsymbol{\mathsf{A}}_2+\boldsymbol{\mathsf{A}}_2\boldsymbol{\mathsf{A}}_1)-\frac{\mathrm{tr}(\boldsymbol{\mathsf{A}}_1\boldsymbol{\mathsf{A}}_2+\boldsymbol{\mathsf{A}}_2\boldsymbol{\mathsf{A}}_1)}{3}\boldsymbol{\mathsf{I}}\right] \nonumber\\ &\quad +\alpha_5 \left[(\boldsymbol{\mathsf{A}}_1^2 \boldsymbol{\mathsf{A}}_2+\boldsymbol{\mathsf{A}}_2\boldsymbol{\mathsf{A}}_1^2)-\frac{\mathrm{tr}(\boldsymbol{\mathsf{A}}_1^2 \boldsymbol{\mathsf{A}}_2+\boldsymbol{\mathsf{A}}_2\boldsymbol{\mathsf{A}}_1^2)}{3}\boldsymbol{\mathsf{I}}\right] \nonumber\\ &\quad +\alpha_6 \left[(\boldsymbol{\mathsf{A}}_1 \boldsymbol{\mathsf{A}}_2^2+\boldsymbol{\mathsf{A}}_2^2\boldsymbol{\mathsf{A}}_1) -\frac{\mathrm{tr}(\boldsymbol{\mathsf{A}}_1 \boldsymbol{\mathsf{A}}_2^2+\boldsymbol{\mathsf{A}}_2^2\boldsymbol{\mathsf{A}}_1)}{3}\boldsymbol{\mathsf{I}}\right] \nonumber\\ &\quad +\alpha_7 \left[(\boldsymbol{\mathsf{A}}_1^2 \boldsymbol{\mathsf{A}}_2^2+\boldsymbol{\mathsf{A}}_2^2\boldsymbol{\mathsf{A}}_1^2)- \frac{\mathrm{tr}(\boldsymbol{\mathsf{A}}_1^2 \boldsymbol{\mathsf{A}}_2^2+\boldsymbol{\mathsf{A}}_2^2\boldsymbol{\mathsf{A}}_1^2)}{3}\boldsymbol{\mathsf{I}}\right], \end{align}

where $p$ is the pressure given by the equation of state $p=\rho RT$, $R$ is the gas constant and $T$ is the temperature. The condition $\alpha _i=0$ reduces this solution to Newton's viscosity law, where stress is linearly proportional to strain rate, given by

(3.8)\begin{equation} {\boldsymbol\tau}=2\mu \boldsymbol{\mathsf{E}},\quad \boldsymbol{\mathsf{E}} =\frac{1}{2} \left(\boldsymbol{\mathsf{A}}_1-\frac{\textrm{tr}(\boldsymbol{\mathsf{A}}_1)}{3}\boldsymbol{\mathsf{I}}\right)= \frac{1}{2}(\boldsymbol{\nabla} {\boldsymbol{v}}+\boldsymbol{\nabla} {\boldsymbol{v}}^{{\rm T}}) -\frac{1}{3}\boldsymbol{\nabla}\boldsymbol{\cdot} {\boldsymbol{v}} . \end{equation}

Since $\tau$ is trace-free, it reduces to only six independent terms of the form

(3.9)\begin{equation} \tau=\beta_i \boldsymbol{\mathsf{S}}_i, \end{equation}

where the $\boldsymbol{\mathsf{S}}_i$ form an orthogonal basis of the form $\boldsymbol{\mathsf{S}}_1 = {\boldsymbol {e}}_1 \otimes {\boldsymbol {e}}_2 + {\boldsymbol {e}}_2 \otimes {\boldsymbol {e}}_1,\ \boldsymbol{\mathsf{S}}_2 = {\boldsymbol {e}}_1 \otimes {\boldsymbol {e}}_3 + {\boldsymbol {e}}_3 \otimes {\boldsymbol {e}}_1,\ \boldsymbol{\mathsf{S}}_3 = {\boldsymbol {e}}_2 \otimes {\boldsymbol {e}}_3 + {\boldsymbol {e}}_3 \otimes {\boldsymbol {e}}_2,\ \boldsymbol{\mathsf{S}}_4 = {\boldsymbol {e}}_1 \otimes {\boldsymbol {e}}_1,\ \boldsymbol{\mathsf{S}}_5 = {\boldsymbol {e}}_2 \otimes {\boldsymbol {e}}_2, \boldsymbol{\mathsf{S}}_6 = {\boldsymbol {e}}_3 \otimes {\boldsymbol {e}}_3$ and $\boldsymbol{\mathsf{S}}_i\boldsymbol{\mathsf{S}}_j=\delta _{ij}$.

A form of the RE model for simple shear which was shown to work reasonably well under the high shear rate $\kappa$ in our previous work is given by (Pahlani et al. 2022)

(3.10)\begin{gather} {\boldsymbol\tau}=\mu\left[\boldsymbol{\mathsf{A}}_1-\frac{\mathrm{tr}(\boldsymbol{\mathsf{A}}_1)}{3}\boldsymbol{\mathsf{I}}\right]+ \alpha_1\left[\boldsymbol{\mathsf{A}}_2-\frac{\mathrm{tr}(\boldsymbol{\mathsf{A}}_2)}{3}\boldsymbol{\mathsf{I}}\right] +\alpha_2 \left[\boldsymbol{\mathsf{A}}_1^2-\frac{\mathrm{tr}(\boldsymbol{\mathsf{A}}_1^2)}{3}\boldsymbol{\mathsf{I}}\right], \end{gather}

(3.11)\begin{gather} \left.\begin{gathered} \mu^*=\frac{1}{1+c_1(s^*)^{c_2}},\quad \alpha_1^*=\frac{1}{2}(c_3+c_4(s^*)^{c_5})^{c_6},\quad \alpha_2={-}2\alpha_1, \\ \mu^*=\mu/\mu_{NSF},\quad \alpha_1^*=p\alpha_1/\mu_{NSF}^2, \end{gathered}\right\} \end{gather}

where the functional forms of the coefficients are inspired from the connection of RE theory with the kinetic theory of Maxwellian molecules in uniform simple shear. Here, $c_i$ are fitting constants listed in table 1, $\mu _{NSF}$ is LJ viscosity corresponding to NSF equations, $\mu$ can be understood as an effective viscosity which, in contrast to Newtonian theory, shows the opposite trend with the dimensionless breakdown parameter $s^*$, where the functional forms of the coefficients are inspired from the connection of RE theory with the kinetic theory of Maxwellian molecules in uniform simple shear. Here, $c_i$ are fitting constants listed in table 1, $\mu _{NSF}$ is LJ viscosity corresponding to NSF equations, $\mu$ can be understood as an effective viscosity which, in contrast to Newtonian theory, shows the opposite trend with the dimensionless breakdown parameter $s^*$, given by

(3.12)\begin{equation} s^*(t)=\frac{\sqrt{2}\kappa}{\nu(T(t))}=\frac{\sqrt{2}\kappa\mu_{NSF}(t)}{p(t)}= \sqrt{2}\frac{|\tau_{12}({\rm NSF})|}{p} . \end{equation}

Here, $\nu$ is collision frequency of the gas. The quantity $s^*$ is also equivalent to product of gradient-length local (GLL) Knudsen number (Boyd, Chen & Candler 1995; Boyd 2003) and local Mach number ${Ma}_L$ given by

(3.13) \begin{align} s^*(t) &= \frac{|{{\rm v}}_1|}{\sqrt{\gamma RT}}\frac{\mu_{NSF}}{p}\sqrt{\frac{{\rm \pi} RT}{2}} \frac{1}{|{{\rm v}}_1|}\left|{\frac{{\rm d}{{\rm v}}_1}{{\rm d}\kern 0.06em x_2}}\right| 2\sqrt{\frac{\gamma}{\rm \pi}} \nonumber\\ &= {Ma_L}\frac{\lambda}{|{\rm v}_1|}\left|{\frac{{\rm d} v_1}{{\rm d}\kern 0.06em x_2}}\right|2\sqrt{\frac{\gamma}{\rm \pi}} ={Ma}_L{Kn}_{GLL}2\sqrt{\frac{\gamma}{\rm \pi}} \nonumber\\ &= \sqrt{2}\gamma \frac{{Ma}_L^2}{{{{Re}}}_L}. \end{align}

Table 1. Fitted parameters of the RE constitutive law.

Note that the relationship between the coefficients $\alpha _1$ and $\alpha _2$ derived here in the above equation holds true for simple shear flow of LJ molecules. The ansatz $\alpha _2=-2\alpha _1$ is derived from OMD simulation which leads to zero second normal stress difference $(\tau _{33}(t)-\tau _{22}(t)=0)$. It may be true that the difference between $\tau _{22}(t)$ and $\tau _{33}(t)$ for simple shear flow of an LJ gas is too small to be captured by MD simulation, whereas for other potentials it may be significant. The behaviour of gas is a strong function of the interatomic potential used and hence the coefficients might need to be recalibrated for other force fields even for the same flow investigated here.

In these bulk flows the characteristic distance associated with the Knudsen number is the velocity-gradient-based length scale. Both of the dimensionless numbers are required to capture the breakdown of quasi-thermodynamic-equilibrium-based NSF model. Large values of ${Ma}_{L} {Kn}_{GLL}$ or ${{Ma}_{L}^2}/{{Re}_{L}}$ correspond to a regime of transition and free-molecular flow where the NSF relations are not valid. Note that $(|{\textrm {v}}_1|/\kappa )$ is used as an intrinsic length to define a local Reynolds number of the flow. Here $s^*$ is also referred to as the local shear-stress Knudsen number in the literature. Ou & Chen (2020) conducted DSMC analysis of wall-bounded rarefied shear flows and also showed that a shear-stress Knudsen number is the only parameter which is needed to characterize the macroscopic flow profile and nonlinear momentum transport in the bulk. This is in accordance with the finding presented here.

In our previous work we showed that under large flow gradients, the RE equation with coefficients given by (3.11) for simple shear significantly improves NSF predictions. In this work, we generalize this calibration to include a bigger family of multiaxial general incompressible and compressible flows. The aim is to model the OMD family of flows by consistently choosing unified variables which define the non-equilibrium physics in a generic manner. This is closely tied to finding breakdown parameters corresponding to Navier–Stokes equations. Various studies have proposed different strategies to define this breakdown regime. A recent study by Singh & Kroells (2021) provides a short review. In our work, dimensionless frame-indifferent invariants of strain rate tensor $\boldsymbol{\mathsf{A}}_1$ are shown to capture the relevant breakdown physics.

3.1. Incompressible flows

For general OMD flows we define a breakdown parameter that takes into account the presence of multiaxial gradients in the flow. It is defined in terms of the second invariant of $\boldsymbol{\mathsf{A}}_1$ by

(3.14)\begin{equation} s^*=\frac{\sqrt{\left(\left[\boldsymbol{\mathsf{A}}_1-\dfrac{\mathrm{tr}(\boldsymbol{\mathsf{A}}_1)}{3}\boldsymbol{\mathsf{I}}\right]: \left[\boldsymbol{\mathsf{A}}_1-\dfrac{\mathrm{tr}(\boldsymbol{\mathsf{A}}_1)}{3}\boldsymbol{\mathsf{I}}\right]\right)}}{\nu}=\frac{2|\boldsymbol{\mathsf{E}}|}{\nu}, \end{equation}

where $|\boldsymbol{\mathsf{E}}|$ denotes the Frobenius norm of a strain rate tensor $\boldsymbol{\mathsf{E}}$. The parameter $s^*$ characterizes the failure of NSF and it is equivalent to truncation number $Tr$ defined by (Truesdell & Muncaster 1980)

(3.15)\begin{equation} Tr=\sqrt{2}\frac{|\boldsymbol{\mathsf{E}}|}{\nu}=\sqrt{2}\frac{\mu}{p}|\boldsymbol{\mathsf{E}}|=\frac{s^*}{\sqrt{2}}. \end{equation}

Here $s^*$ serves as a dimensionless measure of the rate of dissipation of energy through distortion according to classical fluid mechanics, as given by the second term on the right-hand side of the equation for entropy generation,

(3.16)\begin{align} \rho T\dot{\eta} &= \left(\lambda+\tfrac{2}{3}\mu\right) (\boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{v}})^2-{\bf\sigma}:\boldsymbol{\mathsf{E}}-\text{div}\ {\boldsymbol{q}} \nonumber\\ &=\left(\lambda+\tfrac{2}{3}\mu\right)(\boldsymbol{\nabla} \boldsymbol{\cdot}{\boldsymbol{v}})^2+2\mu|\boldsymbol{\mathsf{E}}|^2+\text{div}(k\boldsymbol{\nabla} T), \end{align}

where $\mu,\ (\lambda +\frac {2}{3}\mu )$ and $\kappa$ are shear viscosity, bulk viscosity and thermal conductivity of the fluid, respectively. Here, the Navier–Stokes and Fourier relations are used to obtain (3.16). The first term on the right-hand side is associated with compressibility and bulk viscosity, which vanishes iff $\lambda = -2\mu /3$ (Stokes’ hypothesis/assumption). The third term vanishes for the homogeneous affine motions of OMD.

Equation (3.16) guides us to the correct invariants of the flow which are associated with positive entropy generation and hence can drive the system far from equilibrium beyond the regime of applicability of NSF. In addition to $s^*$ flow compressibility also contributes towards entropy generation. Its effect is taken into account in § 3.2 when dealing with compressible flows.

The usage of $s^*$ is also motivated by its appearance as a natural measure of departure from equilibrium in the kinetic theory (Pahlani et al. 2022). Moreover, this invariant, dimensionless description helps in generalizing this constitutive framework across a wide range of flows, flow gradients and state points of the dilute gas.

It can be seen from (3.10) that the contribution of the coefficients $\alpha _3$ and $\alpha _5$ are not present in simple shear due to identity

(3.17a,b)\begin{equation} \alpha_5={-}\alpha_3,\quad \boldsymbol{\mathsf{A}}_2^2=\boldsymbol{\mathsf{A}}_1^2\boldsymbol{\mathsf{A}}_2+\boldsymbol{\mathsf{A}}_2\boldsymbol{\mathsf{A}}_1^2. \end{equation}

The contributions of these terms in an OMD motion where (3.17a,b) does not hold is investigated here using biaxial shear flow based on $\boldsymbol {A}=\kappa \boldsymbol {e_{2}\otimes e_{1}}+\kappa \boldsymbol {e_{3}\otimes e_{2}}$. The velocity field and breakdown parameter $s^*$ are then given by

(3.18)\begin{equation} \left.\begin{gathered} {\boldsymbol{v}}({\boldsymbol{x}},t)=\boldsymbol{\mathsf{A}}(\boldsymbol{\mathsf{I}}+t\boldsymbol{\mathsf{A}})^{{-}1}{\boldsymbol{x}}=\kappa x_1 {\boldsymbol{e}}_{2}+ (\kappa x_2-\kappa^2 t x_1) {\boldsymbol{e}}_{3},\\ s^*=\frac{\sqrt{(2\kappa^2+(\kappa^2t)^2)2}}{\nu}. \end{gathered}\right\} \end{equation}

The viscous stress tensor becomes

(3.19)\begin{align} {\boldsymbol\tau} &= (\kappa \mu+2\kappa^5t\alpha_3)\boldsymbol{\mathsf{S}}_1+ (-\kappa^2t \mu-2\kappa^2\alpha_1-2\kappa^4\alpha_3 -4\kappa^6t^2\alpha_3)\boldsymbol{\mathsf{S}}_2 \nonumber\\ &\quad +(\kappa \mu+2\kappa^3t\alpha_1+4\kappa^5t\alpha_3)\boldsymbol{\mathsf{S}}_3 \nonumber\\ &\quad +\left(\frac{4}{3}\kappa^2 \alpha_1 +\frac{2}{3}\kappa^4t^2\alpha_1+\frac{4}{3}\kappa^4\alpha_3\right) \boldsymbol{\mathsf{S}}_4 +\left(-\frac{2}{3}\kappa^2 \alpha_1 +\frac{2}{3}\kappa^4t^2\alpha_1 -\frac{8}{3}\kappa^4\alpha_3\right) \boldsymbol{\mathsf{S}}_5\nonumber\\ &\quad +\left(-\frac{2}{3}\kappa^2 \alpha_1 -\frac{4}{3}\kappa^4t^2\alpha_1+ \frac{4}{3}\kappa^4\alpha_3\right) \boldsymbol{\mathsf{S}}_6, \end{align}

where the dimensionless number $\alpha _3^*=\alpha _3p/\mu _{NSF}^2$ can be computed by taking the inner product of (3.19) with $\boldsymbol{\mathsf{S}}_1$. This yields

(3.20)\begin{equation} \alpha_3^*=\frac{\tau_{12}^*-k^*\mu^*(s^*)}{2(\kappa t)(\kappa^{*4})}, \end{equation}

where $\tau ^*=\tau /p$, $\kappa ^*=\kappa /\nu$ and $\mu ^*$ is given by (3.11).

The OMD solution predicts that $\alpha _3<< \mu,\alpha _{1,2}$. Thus, the effect of $\alpha _3$ can be neglected and the final proposed constitutive model only depends on coefficients $\mu$, $\alpha _1$ and $\alpha _2$ which are functions of $s^*$ according to (3.11). To validate this proposed constitutive model, flow features of biaxial-shear $(\kappa =4.6083\times 10^9\ \textrm {s}^{-1})$ are compared with the predictions of OMD.

It is evident from this comparison shown in figure 2 that RE theory closely agrees with OMD and dramatically improves predictions based on the Navier–Stokes solutions. In particular, figure 2 shows that NSF predicts considerably higher shear stresses ($\tau _{12}, \tau _{23}$) than OMD and RE. This happens because of the inability of the NSF to capture viscosity thinning arising from the flow field dependence on the transport coefficients $\mu, \alpha _1, \alpha _2$. Viscosity thinning inhibits momentum transport which reduces the shear stresses present in the system under large gradients. This is accounted for in the RE model by the factor $({1}/{1+c_1(s^*)^{c_2}})$ in the definition of effective viscosity $\mu$. Also in contrast to RE and OMD, NSF theory does not capture normal stress effects $(\tau _{11} \ne \tau _{22} \ne \tau _{33} \ne 0 )$. This is closely tied to the non-coaxiality of ${\boldsymbol \sigma }$ and $\boldsymbol{\mathsf{d}}$ (symmetric part of velocity gradient) which is taken care of by higher-order terms of the form $\boldsymbol{\mathsf{A}}_2$ and $\boldsymbol{\mathsf{A}}_1^2$ in the RE model which are absent in the Newtonian model. Such non-Newtonian behaviour exhibited by dilute gases has also been reported in the literature for uniform simple shear flows of hard-sphere and Maxwell-type molecules (in the Supplementary material, we also investigate the non-equilibrium behaviour of Maxwellian gas under uniform simple shear flow) using the framework of BTE equations and DSMC (Truesdell & Muncaster 1980; Ordóez, Brey & Santos 1989; Garzó & Santos 2003).

Figure 2. Comparison of the OMD evolution (dashed lines) of (a) shear stress and (b) deviatoric normal stresses with theoretical solutions of Navier–Stokes (thick squares and triangles) and RE model (solid lines) for biaxial shear.

The OMD flow considered here never attains a stationary state. This is because external work of the form $({{\boldsymbol \sigma }}\boldsymbol {\cdot }\boldsymbol {\nabla } {\boldsymbol {v}})$ is continuously being done on the system and there is no balancing source of energy dissipation. This results in a monotonic increase in temperature of the gas due to viscous heating, since we do not introduce any artificial thermostats that might raise other issues of validity. At given $|\boldsymbol{\mathsf{E}}|$, the increase in temperature increases the collision frequency of the gas which results in decrease in departure of the system from equilibrium, as measured by $s^*$. In other words, as the simulation evolves, the flow approaches near-equilibrium conditions. Note that the nonlinear transport guided behaviour of gas is realized when parameter $s^*$ is comparatively higher. This is achieved either by increasing the numerator (introduction of sharp gradients) or decreasing the denominator (low collision frequency of the gas, i.e. making the gas more dilute). In this work, we capture similar non-equilibrium physics in rarefied monoatomic gas by increasing the temporal gradients in the system. As the gas becomes more and more dilute, the strain rate required to capture the NSF breakdown decreases (Pahlani et al. 2022).

Next, we consider a general isochoric case, that of general incompressible unsteady flows. The necessary condition for this to hold true is $\boldsymbol {\nabla }\boldsymbol {\cdot } {\boldsymbol {v}} =0$. This condition of incompressibility is equivalent in Lagrangian form to det$({\boldsymbol{\mathsf{I}}}+t{\boldsymbol{\mathsf{A}}})=1$ for all $t>0$.

Necessary and sufficient conditions for $\textrm {det}({\boldsymbol{\mathsf{I}}}+t{\boldsymbol{\mathsf{A}}})=1$ for all $t>0$ are that $\textrm {det}\ \boldsymbol{\mathsf{A}}= \textrm {tr}\,\boldsymbol{\mathsf{A}}=\textrm {tr}({\boldsymbol{\mathsf{A}}}^{2})=0$. In a suitable orthonormal basis necessary and sufficient conditions are that

(3.21)\begin{equation} \boldsymbol{\mathsf{A}}=\begin{bmatrix} 0 & \kappa & 0 \\ 0 & 0 & 0 \\ \gamma_1 & \gamma_2 & 0 \end{bmatrix} \end{equation}

in this basis. In abstract form, $\boldsymbol{\mathsf{A}}=\kappa {\boldsymbol {e}}_1\otimes {\boldsymbol {e}}_2+{\boldsymbol {e}}_3\otimes {\boldsymbol {g}}$, where ${\boldsymbol {e}}_1, {\boldsymbol {e}}_2, {\boldsymbol {e}}_3$ are orthonormal and ${{\boldsymbol {g}}}=\gamma _1{{\boldsymbol {e}}_1}+\gamma _2{\boldsymbol {e}}_2$.

The corresponding Eulerian description of this motion is given by

(3.22)\begin{equation} {{\boldsymbol{v}}}({{\boldsymbol{x}}},t)={\boldsymbol{\mathsf{A}}}{{\boldsymbol{x}}}-\kappa t \gamma_1 {\boldsymbol{e}}_3. \end{equation}

This flow differs from traditional viscometric flows by the presence of a time-dependent vorticity, $\textrm {curl} \,{\boldsymbol {v}}=(\gamma _2-\kappa \gamma _1t){\boldsymbol {e}}_1-\gamma _1{\boldsymbol {e}}_2-\kappa {\boldsymbol {e}}_3$. Figure 3 compares the OMD viscous stress evolution with that of RE model for $\gamma _1=\kappa =\gamma _2= 2.3041 \times 10^9\ \textrm {s}^{-1}$. The RE theory is shown to agree with OMD predictions reasonably well even under these extreme conditions.

Figure 3. Comparison of the OMD evolution (dashed lines) of (a) shear stress and (b) deviatoric normal stresses with theoretical solutions of RE model (solid lines) for general unsteady incompressible flows.

3.2. Compressible flows

The OMD also provides a method to simulate a nine-parameter family $(\boldsymbol{\mathsf{A}})$ of compressible flows where $\boldsymbol {\nabla }\boldsymbol {\cdot } {\boldsymbol {v}}$ is non-zero. In this case a finite value of the first invariant of $\boldsymbol{\mathsf{A}}_1$ needs to be considered when defining the breakdown regime of gas flows. The evolution of density in these compressible flows is given by

(3.23)\begin{equation} \rho = \rho(0)\exp\left(-\int_0^t E\ {{\rm d}}s\right)=\rho(0) \exp\left(-\int_0^t{{\rm tr}}\ {\boldsymbol{\mathsf{A}}}({\boldsymbol{\mathsf{I}}}+s{\boldsymbol{\mathsf{A}}})^{{-}1}\,{{\rm d}}s\right),\end{equation}

where $\rho (0)$ is the initial density and $E$ is expansion given by $E=\boldsymbol {\nabla }\boldsymbol {\cdot }{\boldsymbol {v}}$. Expansion of the gas can lead to a decrease in the molecular collision rate, which in turn can contribute to a failure of the NSF equations. For expanding flows Bird proposed the following parameter for continuum breakdown (Bird 1970; Boyd 2003):

(3.24)\begin{align} P&=\frac{1}{\nu}\left(\frac{{\rm D}({\ln}\rho)}{{\rm D} t}\right)=\frac{1}{\rho\nu} \left(\frac{{\rm D}\rho}{{\rm D} t}\right), \nonumber\\ &=\frac{1}{\rho\nu}(-\rho\boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{v}})=\frac{-\textrm{tr}({\boldsymbol{A}} ({\boldsymbol{I}}+t{\boldsymbol{A}})^{{-}1})}{\nu}=\frac{-\textrm{tr}({\boldsymbol{A}} ({\boldsymbol{I}}+t{\boldsymbol{A}})^{{-}1})\mu_{NSF}}{p}, \end{align}

where $D$ is the total derivative. We use this criterion to incorporate terms of the form $P=\rho ^*/\rho$ in the earlier computed coefficients $\mu ^*$ and $\alpha _1^*$ to include the contribution of time-dependent density variations. We suggest the following coefficients of the RE model which work well for both compressible and incompressible flows:

(3.25)\begin{equation} \left.\begin{gathered} \mu^*=\frac{1}{1+c_1(s^*)^{c_2}\left(1-c_{7}\dfrac{P}{s^*}\right)^{c_8}},\\ \alpha_1^*=\frac{1}{2}(c_3+c_4(s^*)^{c_5}+c_9({-}P)^{c_{10}})^{c_6}, \end{gathered}\right\} \end{equation}

where $c_i,\ i=1,\ldots,6$, remains the same as that derived from incompressible flows. The $c_i, i=7,\ldots 10$, will now be calibrated using a family of compressible flows.

To extend the calibration of RE equations, we consider the pressure-shear flow field $(\boldsymbol{\mathsf{A}}=\kappa ({\boldsymbol {e}}_1\otimes {\boldsymbol {e}}_1+{\boldsymbol {e}}_1\otimes {\boldsymbol {e}}_2))$ which was compared in § 2 with OMD. The RE model can be used to obtain the viscous stress given by

(3.26) \begin{align} \left.\begin{aligned} {\boldsymbol\tau}&=\left(\frac{\kappa \mu}{(1+\kappa t)}-\frac{2\kappa^2\alpha_1}{(1+\kappa t)^2}\right)\boldsymbol{\mathsf{S}}_1+ \left(\frac{4}{3}\frac{\kappa \mu}{(1+\kappa t)}-\frac{16}{3}\frac{\kappa^2\alpha_1}{(1+\kappa t)^2}\right)\boldsymbol{\mathsf{S}}_4 \\ &\quad+\left(-\frac{2}{3}\frac{\kappa \mu}{(1+\kappa t)}+\frac{8}{3}\frac{\kappa^2\alpha_1}{(1+\kappa t)^2}\right)\boldsymbol{\mathsf{S}}_5+ \left(-\frac{2}{3}\frac{\kappa \mu}{(1+\kappa t)}+\frac{8}{3}\frac{\kappa^2\alpha_1}{(1+\kappa t)^2}\right)\boldsymbol{\mathsf{S}}_6, \end{aligned}\right\} \end{align}

and $\mu ^*$ and $\alpha _1^*$ can be obtained by using identities

(3.27)\begin{equation} \left.\begin{gathered} \mu^*=\left({\boldsymbol\tau}^*:\boldsymbol{\mathsf{S}}_1-\frac{3}{8}{\boldsymbol\tau}^*:\boldsymbol{\mathsf{S}}_4\right)\frac{2(1+\kappa t)}{\kappa^*},\\ \alpha_1^*=\left(\frac{4}{3}{\boldsymbol\tau}^*:\boldsymbol{\mathsf{S}}_1-{\boldsymbol\tau}^*:\boldsymbol{\mathsf{S}}_4\right) \frac{(1+\kappa t)^2}{\kappa^{*^2}}\frac{3}{8}. \end{gathered}\right\} \end{equation}

These coefficients are fitted to obtain the values listed in table 1. In figure 4 we compare the fitted data with the OMD evolution of $\mu ^*$ and $\alpha _1^*$ as a function of $s^*$ and $P$. Corresponding numerical values are given in tables 2 and 3. These are obtained by simulating hundreds of flows with varying $\kappa$ and state points $(\rho, T)$. Note that both $s^*$ and $P$ are invariants of tensor $\boldsymbol{\mathsf{d}}$ and hence are not independent variables. It can be seen from (3.25) that the generalized transport coefficients reduce to those for Navier–Stokes in the limit of zero gradients.

Figure 4. Fitting of the coefficients of RE model for universal flows of LJ argon gas.

Table 2. The OMD computed coefficients for compressible pressure shear flow.

Table 3. The OMD computed coefficients for incompressible simple shear flow $(P=0)$.

3.2.1. One-dimensional expansion and compression

The final fitted model is next validated by simulating one-dimensional (1-D) extensional and compression flow where the gas grows steadily rarer and denser, respectively. The input matrix is chosen to be

(3.28)\begin{equation} \boldsymbol{\mathsf{A}}=\begin{bmatrix} \kappa & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix}. \end{equation}

The velocity field is given by $\textrm {v}_1=({\kappa }/{(\kappa t+1)})\textrm {x}_1, \textrm {v}_2=0, \textrm {v}_3=0$; where $\kappa >0$ results in expansion and $\kappa <0$ results in compression. The RE viscous stress tensor for this flow reduces to

(3.29)\begin{align} {\boldsymbol\tau} &= \left(\frac{4}{3}\frac{\kappa\mu}{(1+\kappa t)}-\frac{4\kappa^2\alpha_1}{(1+\kappa t)^2}\right)\boldsymbol{\mathsf{S}}_4 + \left(-\frac{2}{3}\frac{\kappa \mu}{(1+\kappa t)}+\frac{2\kappa^2\alpha_1}{(1+\kappa t)^2}\right)\boldsymbol{\mathsf{S}}_5 \nonumber\\ &\quad +\left(-\frac{2}{3}\frac{\kappa \mu}{(1+\kappa t)}+\frac{2\kappa^2\alpha_1}{(1+\kappa t)^2}\right)\boldsymbol{\mathsf{S}}_6, \end{align}

where $\tau _{22}=\tau _{33}=-\tau _{11}/2$ and shear stresses are not present in the system. Figure 5 compares the NSF and RE predictions of stresses with OMD results for $\rho (0)=0.1784\ \textrm {kg}\ \textrm {m}^{-3},\ T(0)=10\,000\ \textrm {K},\ \kappa =1.3825\times 10^9\ \textrm {s}^{-1}$. The RE model is shown to agree with OMD results reasonably well. On the other hand, the NSF model overpredicts stresses in the system by a significant amount. Figure 6(a) compares the evolution of temperature in the case of 1-D expansion. The initial equilibrium state of the system is given by $\rho (0)=0.1784\ \textrm {kg}\ \textrm {m}^{-3},\ T(0)=10\,000\ \textrm {K}$. Non-equilibrium conditions in the gas are seen by the development of anisotropy, where the translational temperature of the gas associated with each coordinate direction $T_{e_{\alpha }}=\langle ({m}/{k_b N})[\sum _{i}^{N}({v^{\prime}_{{\alpha },i}})^2]\rangle,\ \alpha =1, 2, 3,$ separates. The temperature in the ${\boldsymbol {e}}_1$ direction $T_{{\boldsymbol {e}}_1}$ is widely different from the other two directions and decreases as the simulation evolves. The OMD evolution of the total temperature agrees with the RE prediction exceedingly well. On the other hand, NS theory contrasts sharply with the observed behaviour. The predicted temperature by RE theory is computed by incorporating a RE constitutive relation in the ordinary differential equation for temperature, (2.3). Similarly, figure 6(b) shows various temperature evolutions for a 1-D compression case where $\kappa$ is non-positive and is given by $\kappa =-2.3041\times 10^8\ \textrm {s}^{-1}$. The initial equilibrium state of the system is: $\rho (0)=0.01784\ \textrm {kg}\ \textrm {m}^{-3},\ T(0)=300\ \textrm {K}$. In contrast to the previous case, the temperature of the system increases with $T_{{\boldsymbol {e}}_1}$ being far higher than $T_{{\boldsymbol {e}}_2}=T_{{\boldsymbol {e}}_3}$. As shown, the RE theory is able to capture the behaviour of gas significantly better than NS. Note that since $\kappa$ is negative, this unsteady compressible flow has a singularity at $t=\tilde {t}={1}/{\kappa }$. Therefore, we stop the simulation before $t$ reaches $\tilde {t}$. Interestingly, from these simulations it can be inferred that a multitemperature theory might not be required for the continuum description of translational non-equilibrium, given the correct nonlinear terms. These effects are automatically incorporated in the underlying RE constitutive model with one temperature.

Figure 5. Comparison between OMD, NSF and RE theory for argon gas under 1-D expansion.

Figure 6. Different temperature evolution comparison between OMD, RE and NSF solution for (a) 1-D expansion, $\kappa = 1.8433 \times 10^9$ s$^{-1}$ and (b) 1-D compression, $\kappa = -2.3041 \times 10^8$ s$^{-1}$.

In figure 7 we examine the 1-D VDF $g(w_1)$ for 1-D compression and expansion, respectively. This is obtained by averaging over all the thermal velocities in the ${\boldsymbol {e}}_2$ and ${\boldsymbol {e}}_3$ directions. To compute $g(w_1)$ from OMD, the fundamental domain is divided into equally spaced bins in the ${\boldsymbol {e}}_1$ direction and the population of atoms having a certain thermal velocity is analysed. For both the computations, we start with a Maxwellian velocity distribution. As the gas evolves under compression, the temperature increases and the distribution gradually flattens and attains a non-Boltzmann character. It is apparent that high velocity tails are significantly overpopulated as compared with a local Maxwellian evaluated at that instant. Similarly, under expansion, the distribution gradually shifts from Maxwellian at $t=0$. Under expansion, the temperature decreases with time and the distribution becomes narrower. Unlike the previous system, the population of low velocity bins is high compared with that of the local equilibrium distribution.

Figure 7. The VDF $g(w_1)$ evolution for 1-D compression. Blue and black dashed lines depict Maxwellians evaluated at local flow conditions.

We conclude that RE theory has the potential to provide a rather complete description of universal flows (OMD flows) in situations far from equilibrium. It makes a significant improvement on the Newtonian transport model. The final form for LJ gas is given by

(3.30)\begin{align} {\boldsymbol\tau}=\mu^*\mu_{NSF}\left[\boldsymbol{\mathsf{A}}_1-\frac{\mathrm{tr}(\boldsymbol{\mathsf{A}}_1)}{3}\boldsymbol{\mathsf{I}}\right]+ \frac{\alpha_1^*\mu_{NSF}^2}{p} \left[\left(\boldsymbol{\mathsf{A}}_2-\frac{\mathrm{tr}(\boldsymbol{\mathsf{A}}_2)}{3}\boldsymbol{\mathsf{I}}\right) -2\left(\boldsymbol{\mathsf{A}}_1^2-\frac{\mathrm{tr}(\boldsymbol{\mathsf{A}}_1^2)}{3}\boldsymbol{\mathsf{I}}\right)\right], \end{align}

where $\mu ^*$ and $\alpha _1^*$ are given by (3.25). This generalizes the one proposed in our previous work focused on uniform simple shear. The coefficients are generalized to take into account multiaxial gradients and compressibility effects. In our OMD simulations, the temperature and density of the gas vary with time due to viscous heating with no artificial thermostats and time-dependent velocity gradients. Thus, the RE model is calibrated using the temperature and density-dependent behaviour of gas and hence is universal across different state points. In addition to $\mu _{NSF}$, further dependence on temperature occurs through breakdown parameters $s^*(T)$ and $P(T)$ on which the coefficients of the RE model depend.

There exists a class of monoatomic flows for which the prediction of kinetic theory is in exact agreement with Navier–Stokes and Euler equations. This is pure dilatation of fluid with null bulk viscosity and sufficiently high density (still within dilute gas regime) to have a valid continuum description. The $\boldsymbol{\mathsf{A}}$ tensor for the flow is given by $k{\boldsymbol {e}}_1\otimes {\boldsymbol {e}}_1+k{\boldsymbol {e}}_2\otimes {\boldsymbol {e}}_2+k{\boldsymbol {e}}_3\otimes {\boldsymbol {e}}_3$. It is a dissipationless flow and thus indistinguishable from the corresponding flow of an inviscid fluid (Truesdell & Muncaster 1980). The model proposed in this work agrees with these results, i.e. $s^*$ vanishes for these flows and ${\boldsymbol \sigma }=p\boldsymbol{\mathsf{I}},\ \tau = 0,\ \forall \ k$.

It is interesting to note that there exists a non-classical relationship between the two coefficients $\mu$ and $\alpha _1$ as shown in figure 8. They closely obey a power law form given by

(3.31)\begin{equation} \alpha_1^*=0.282\mu^*{^{2.771}}. \end{equation}

This can be interpreted as the existence of an unexpected scaling relationship between shear and normal stresses in the flow. The idea of having a non-classical stress constraint was also suggested by Myong & Park (2012) for simple shear using DSMC computations.

Figure 8. Non-classical constitutive relationship between coefficients of RE model of complexity two for dilute gas universal flows.

Next we look at the thermodynamics and stability aspects of the RE model.

4. Thermodynamics and stability of RE fluid of complexity 2

A well-known special example of RE fluids are fluids of second grade which satisfy the following constitutive model:

(4.1) \begin{equation} \left.\begin{gathered} \boldsymbol{\mathsf{M}}=\tfrac{1}{3}\text{tr}(\boldsymbol{\mathsf{M}})\boldsymbol{\mathsf{I}}+\alpha_0\boldsymbol{\mathsf{I}}+\mu(T)\boldsymbol{\mathsf{A}}_1+\alpha_1(T){\boldsymbol{\mathsf{A}}_2+\alpha_2(T)\boldsymbol{\mathsf{A}}_1^2},\\ \tilde{\boldsymbol{\mathsf{M}}}=\boldsymbol{\mathsf{M}}-\tfrac{1}{3}(\text{tr}(\boldsymbol{\mathsf{M}})\boldsymbol{\mathsf{I}})={-}(\sigma-p\boldsymbol{\mathsf{I}}) \\ ={\boldsymbol\tau}=\alpha_0\boldsymbol{\mathsf{I}}+\mu(T)\boldsymbol{\mathsf{A}}_1+\alpha_1(T){\boldsymbol{\mathsf{A}}_2+\alpha_2(T)\boldsymbol{\mathsf{A}}_1^2}, \end{gathered}\right\} \end{equation}

where coefficients $\mu, \alpha _1, \alpha _2$ are constants which may depend on temperature. Note that Cauchy stress tensor $\boldsymbol{\mathsf{M}}$ is negative of the stress tensor ${\boldsymbol \sigma }$ defined in this work.

When using this model as a general constitutive relation, it is known to cause an unphysical response when the coefficients are not compatible with thermodynamics (Fosdick & Rajagopal 1979; Dunn 1982). It is shown that instability and unboundedness are unavoidable when this inconsistency occurs (Coleman, Duffin & Mizel 1965; Coleman & Mizel 1966). The RE-type second-grade model in the literature has typically been questioned due to the discrepancy between thermodynamics and experimental measurements of the coefficients. Dunn & Fosdick studied stability and thermodynamics for second-grade fluid and showed the following conditions on coefficients to be true (Dunn & Fosdick 1974):

(4.2a–c)\begin{equation} \mu(T)\ge 0,\quad \alpha_1(T)\ge 0,\quad \alpha_1(T)+\alpha_2(T)=0 . \end{equation}

A large body of rheological measurements on polymeric fluids suggest otherwise. For example, consider a simple shear flow $(\boldsymbol{\mathsf{A}}=\kappa {\boldsymbol {e}}_1\otimes {\boldsymbol {e}}_2)$. For RE fluid of second grade, we can define first normal stress difference as

(4.3)\begin{equation} N_1={-}(\sigma_1-\sigma_2)=(\tau_1-\tau_2)={-}2\alpha_1\kappa^2. \end{equation}

Experimenters find $N_1>0$, which is inconsistent with the restriction imposed by thermodynamics. This leads to rejection of the idea of the second-grade model to be exact in its own right. On the other hand, our molecular simulation results on dilute monoatomic gas seems to be consistent with this analysis as well as its extension by Dunn for a new class of complexity 2 incompressible fluids considered in this work (Dunn & Fosdick 1974). For these fluids the viscous stress tensor is given by

(4.4)\begin{equation} {\boldsymbol\tau}=\alpha_0\boldsymbol{\mathsf{I}}+\mu(T,\boldsymbol{\mathsf{A}}_1)\boldsymbol{\mathsf{A}}_1+\alpha_1(T,\boldsymbol{\mathsf{A}}_1)\boldsymbol{\mathsf{A}}_1^2+\alpha_2(T,\boldsymbol{\mathsf{A}}_1)\boldsymbol{\mathsf{A}}_2,\end{equation}

where the coefficients are arbitrary isotropic functions of temperature as well as $\boldsymbol{\mathsf{A}}_1$. This family of models is wide enough to include ‘generalized-Newtonian fluids’ $(\alpha _1=0, \alpha _2=0, \mu =\mu (|\boldsymbol{\mathsf{A}}_1^2|,T))$ and ‘second-grade fluids’. It was shown that for second law of thermodynamics to hold and for the Helmholtz free energy to have a minimum in equilibrium, the following conditions are necessary and sufficient:

(4.5)$$\begin{gather} \alpha_1(T,\boldsymbol{\mathsf{A}}_1)=\alpha_1(T, |\boldsymbol{\mathsf{A}}_1|) , \end{gather}$$

(4.6)$$\begin{gather}\int_{0}^{z}\alpha_1(T,\zeta)\,{\rm d}\zeta \ge 0, \quad z\in[0,\infty)\quad \text{for } z\in[0,\epsilon), \epsilon>0 , \end{gather}$$

(4.7)\begin{gather} \left.\begin{gathered} \alpha_1(T,\boldsymbol{\mathsf{A}}_1)+\alpha_2(T,\boldsymbol{\mathsf{A}}_1) \to 0 \\ \text{for path } \boldsymbol{\mathsf{A}}_1=\boldsymbol{\mathsf{A}}_1(\tau),|\boldsymbol{\mathsf{A}}_1(\tau)|\to \infty \quad \text{as } \tau \to \infty \end{gathered}\right\}. \end{gather}

The coefficients calibrated in this work are functions of $s^*$ and $P$, which in turn are functions of $|\boldsymbol{\mathsf{A}}_1|$; thus they obey (4.5).

The condition given by (4.6) is satisfied since the coefficient $\alpha _1$ computed is positive everywhere in the domain; and the last condition given by (4.7) is applicable since

(4.8) \begin{equation} \left.\begin{gathered} \alpha_1(T,\boldsymbol{\mathsf{A}}_1)+\alpha_2(T,\boldsymbol{\mathsf{A}}_1)={-}\alpha_1(T,\boldsymbol{\mathsf{A}}_1) \\ ={-}\frac{\mu_{NSF}^2}{p}\frac{1}{2}(c_3+c_4(s^*)^{c_5}+c_9({-}P)^{c_{10}})^{c_6} \to 0 \\ \textrm{as}\ |\boldsymbol{\mathsf{A}}_1(\tau)|\to \infty, \textrm{for}\ c_6<0. \end{gathered}\right\} \end{equation}

Thus, the final rheological model proposed and calibrated in this work, in addition to being properly invariant, also satisfies necessary restrictions coming from the second law of thermodynamics for an incompressible fluid. This analysis needs to be extended to compressible flows to complete the argument. This investigation is beyond the scope of this work.

5. Connection with Burnett equations

The Burnett equations are suggested as an extension of the Navier–Stokes equations which are derived from CE expansion where the Maxwellian distribution $\widehat {f_0}$ is expanded for small values of $\zeta$ to obtain a VDF of the system,

(5.1)\begin{equation} \hat{f}=\hat{f}_0(1+\zeta\phi_1+\zeta^2\phi_2+\cdots), \end{equation}

where $\zeta$ determines the degree of departure from the distribution. It can be written in a physically meaningful manner (Boyd & Schwartzentruber 2017),

(5.2)\begin{equation} \zeta=\left(\frac{\text{mean-collision-time}}{\text{characteristic flow time}}\right)= \left(\frac{\text{mean free path}}{\text{characteristic length}}\right), \end{equation}

where the relation between mean free path and collision frequency is used. Here $\zeta$ is equivalent to breakdown parameter $s^*$ defined in this work.

Substitution of the expansion of the VDF for $\hat {f}$ into the non-dimensionalized Boltzmann equation yields a general constitutive relation for the stress tensor of the form

(5.3)\begin{equation} {\boldsymbol\tau}_{ij}={\boldsymbol\tau}_{ij}^0+{\boldsymbol\tau}_{ij}^1+{\boldsymbol\tau}_{ij}^2+\cdots+o(\zeta^n). \end{equation}

The zeroth- and first-order approximation yields the Euler and Navier–Stokes equations. Retaining terms up to second order provides the Burnett equation given by (Comeaux et al. 1995)

(5.4) \begin{align} {\boldsymbol\tau}={\boldsymbol\tau}^{1}+{\boldsymbol\tau}^{2} &={-}\left({-}2\mu_{NSF}\overline{\boldsymbol{\nabla}{\boldsymbol{v}}}+\omega_1 \frac{\mu^2}{p}\boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{v}}\overline{\boldsymbol{\nabla}{\boldsymbol{v}}}+ \omega_2\frac{\mu^2}{p}\left(\frac{{\rm D}}{{\rm D} t}\overline{\boldsymbol{\nabla}{\boldsymbol{v}}}-2 \overline{\boldsymbol{\nabla}{\boldsymbol{v}}\overline{\boldsymbol{\nabla}{\boldsymbol{v}}}}\right) \right. \nonumber\\ &\quad +\left.\omega_3\frac{\mu^2}{\rho T}\overline{\boldsymbol{\nabla}\boldsymbol{\nabla} T}+\omega_4 \frac{\mu^2}{\rho p T}\overline{\boldsymbol{\nabla} p \boldsymbol{\nabla} T} +\omega_5 \frac{\mu^2}{\rho T^2}\overline{\boldsymbol{\nabla} T \boldsymbol{\nabla} T}+\omega_6 \frac{\mu^2}{p}\overline{\overline{\boldsymbol{\nabla} {\boldsymbol{v}}}\,\overline{\boldsymbol{\nabla} {\boldsymbol{v}}}}\right), \end{align}

where ${\boldsymbol \tau }^1$ and ${\boldsymbol \tau }^2$ show Navier–Stokes and Burnett contributions, respectively. Time derivatives in (5.4) are eliminated using the conservation relation in the conventional Burnett equations, but we stick to this original form to establish its connection with RE theory. The $\omega _{i}$ are constant coefficients which depend on the interaction force field model. These are derived for some simple intermolecular force fields such as for Maxwellian gas (Comeaux et al. 1995). For incompressible OMD flow (5.4) reduces to

(5.5)\begin{align} {\boldsymbol\tau} &={-}\left(-\mu_{NSF}\boldsymbol{\mathsf{A}}_1+\omega_2\frac{\mu^2}{2p}\boldsymbol{\mathsf{A}}_2+\left(\omega_6 \frac{\mu^2}{2p}\right)\boldsymbol{\mathsf{A}}_1^2+\omega_2\frac{\mu^2}{p}(\boldsymbol{\nabla} {\boldsymbol{v}}^\mathrm{T} \boldsymbol{\mathsf{A}}_{1}-\boldsymbol{\mathsf{A}}_{1}\boldsymbol{\nabla} {\boldsymbol{v}}) \right.\nonumber\\ &\quad +\left.\omega_2\frac{2\mu^2}{3p}\left({\rm tr}\left(\boldsymbol{\nabla}{\boldsymbol{v}}^{{\rm T}} \frac{\boldsymbol{\mathsf{A}}_1}{2}\right)\boldsymbol{\mathsf{I}}\right)-\omega_6 \frac{\mu^2}{3p}({\rm tr}(\boldsymbol{\mathsf{A}}_1^2)\boldsymbol{\mathsf{I}})\right), \end{align}

where the $\boldsymbol{\mathsf{A}}_i$ are RE tensors. For OMD compressible flows, (5.4) reduces to a similar form but includes more terms in the summation which depends on the non-zero divergence of velocity gradient $(\boldsymbol {\nabla }\boldsymbol {\cdot }{\boldsymbol {v}})$.

Writing the Burnett equation in this form shows its connection with the RE constitutive model. We see here that Burnett equations also contain terms which are of the form $\boldsymbol{\mathsf{A}}_1$, $\boldsymbol{\mathsf{A}}_2$ and $\boldsymbol{\mathsf{A}}_1^2$ in accordance with the RE equation but these are not frame-indifferent. The applicability of the principle of material frame invariance (PMFI) for the kinetic theory of gases under extreme deformation rates is a widely studied unsettled subject. Müller (Müller 1972) employed an iterative scheme developed by Ikenberry & Truesdell to derive material equations for a gas of Maxwellian molecules from Boltzmann's equation which is of Burnett's form. The dependence of stresses and heat flux on the spin tensor indicates the limitations of the range of validity of PMFI. It was remarked that frame dependence in the sense mentioned above is due to the action of microscopic Coriolis forces (Müller 1972; Söderholm 1976). It was remarked by Truesdell (1976) that inertial effects of the kind obtained by Müller arise as a consequence of the principle of linear momentum, a notoriously frame-dependent statement. Subsequently, with the work on extended thermodynamics (Liu & Müller 1986), this frame dependence was seen in a new light by Müller. The basic equations of balance or the field equations are believed to reflect the material frame dependence, while the theory is frame independent in the constitutive relations. Altogether, there is no consensus on the general applicability of PMFI and different schools of thought exist. We believe that the iterates adopted in Müller (1972) are approximations of the balance equation itself, therefore may contain a contribution from the Coriolis force due to the rotation of the frame in a non-inertial frame. The constitutive functions should be observer-independent in a certain sense since they characterize the intrinsic properties of the material body itself. We show that with the workhorse of continuum classical theory with PMFI incorporated, our proposed frame independent RE constitutive relation improves on the Navier–Stokes theory for universal flows. This needs to be further extended for other flows with thermal transport to make a more conclusive statement.

In our constitutive framework, coefficients depend on the strength of flow field gradients captured by breakdown parameters $s^*$ and $P$ as opposed to the Burnett model. Similar to the Burnett equation, a non-equilibrium VDF underlies the RE model; see figure 9. Evidently, one needs to resort to a non-perturbative approach to find its full form. The spatial variation of density and temperature may also affect nonlinear momentum transport as postulated by Burnett's theory. Since OMD generates time-dependent stress, density and temperature fields that are spatially homogeneous ($\boldsymbol {\nabla } T=\boldsymbol {\nabla }\rho =0$) in this work, it does not capture departure from equilibrium due to the spatial variation of these fields; hence, OMD captures departure from equilibrium only in the temporal frame.

Figure 9. The VDF $g(w_1)$ evolution for 1-D expansion. Dashed lines depict Maxwellians evaluated at local flow conditions.

6. Conclusion

In this paper we introduce a nonlinear, non-classical constitutive relation for highly non-equilibrium flows of gases based on the use of the RE tensors. We show that the resulting RE theory gives results that agree well with the NSF theory in near-equilibrium regimes but also performs well in diverse non-equilibrium where the NSF theory fails, in some cases severely. We also introduce diagnostic parameters $s^*$ and $P$ that are sufficient to assess the applicability of NSF equations in a variety of flows. The identification of these parameters opens the way to a hybrid RE/NSF method.

The method is applicable to highly non-equilibrium situations where the velocity distribution function departs significantly from a Maxwellian distribution. The success of our method indicates that the paradigm of constitutive relations as closure conditions for the continuum balance laws is a valid approach in this regime. Though the modelling and calibration is limited to LJ argon gas, the framework adopted is quite general. We think that the fitted parameters will have a strong dependence on the particular force field but the (non-dimensionalized) functional forms may be generic across various monatomic gases. The extension of the proposed momentum transport model for diatomic gas will require consideration on bulk viscosity (to account for the deviation of average normal stress from pressure) which is non-zero for diatomic flows in non-equilibrium. Also, the complexity associated with modelling diatomic flows at high-temperatures is further increased due to the presence of active internal energy modes and chemical reactions. The OMD method is equally capable for developing better thermochemistry continuum models in a far from equilibrium regimes.

Essentially, we have used OMD as a method of computational viscometry in this work. This is particularly effective because the macroscopic motion ${\boldsymbol {v}}({\boldsymbol {x}}, t) =\boldsymbol{\mathsf{A}}(\boldsymbol{\mathsf{I}}+t\boldsymbol{\mathsf{A}})^{-1}{\boldsymbol {x}}$ for any $3 \times 3$ matrix $\boldsymbol{\mathsf{A}}$, together with a suitable temperature $T(t)$ satisfying an explicit ordinary differential equation, are exact solutions of the balances of mass, momentum and energy together with general constitutive relations. The atomistic simulations produce exactly this macroscopic motion and a time-dependent temperature. The OMD works for any force field including those generated by adiabatic quantum mechanics, and also works for molecular gases, liquids or solids. It therefore can handle phase change (Pahlani et al. 2023) without ad hoc assumptions.

However, it is important to note that OMD only captures nonlinear momentum transport, since the heat flux is identically zero in OMD. Our future aim is to explore further if our RE constitutive relation for the stress tensor is also suitable for other flows where energy transport is as significant as momentum transport, which may lead to coupling. There is a significant body of literature which suggests the generation of stress due to temperature gradients under highly non-equilibrium regimes (Maxwell 1879; Sharipov & Kremer 1999). However, we do not know of any universally accepted constitutive relation for any material having this dependence.