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 Reference Kim and Monroe2014) 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. Reference Pahlani, Schwartzentruber and James2023). 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 Reference Boyd and Schwartzentruber2017)

(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 Reference Dayal and James2012) 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 Reference Rivlin and Ericksen1955), 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 Reference Coleman, Markovitz and Noll2012). 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 Reference Green and Rivlin1957) 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. Reference Coleman, Markovitz and Noll2012). (Generally, OMD flows and viscometric flows are different.) Rivlin used this reduction to find some exact solutions for simple shearing and torsional flow (Rivlin Reference Rivlin1956). 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 (Reference Bird and Huilgol1999) 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. (Reference Coleman, Markovitz and Noll2012). 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. Reference Pahlani, Schwartzentruber and James2022)

(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 Reference Boyd, Chen and Candler1995; Boyd Reference Boyd2003) 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 (Reference Ou and Chen2020) 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 (Reference Singh and Kroells2021) 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 Reference Truesdell and Muncaster1980)

(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. Reference Pahlani, Schwartzentruber and James2022). 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 Reference Truesdell and Muncaster1980; Ordóez, Brey & Santos Reference Ordóez, Brey and Santos1989; Garzó & Santos Reference Garzó and Santos2003).

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. Reference Pahlani, Schwartzentruber and James2022).

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 Reference Bird1970; Boyd Reference Boyd2003):

(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 Reference Truesdell and Muncaster1980). 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 (Reference Myong and Park2012) 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 Reference Fosdick and Rajagopal1979; Dunn Reference Dunn1982). It is shown that instability and unboundedness are unavoidable when this inconsistency occurs (Coleman, Duffin & Mizel Reference Coleman, Duffin and Mizel1965; Coleman & Mizel Reference Coleman and Mizel1966). 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 Reference Dunn and Fosdick1974):

(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 Reference Dunn and Fosdick1974). 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 Reference Boyd and Schwartzentruber2017),

(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. Reference Comeaux, Chapman and MacCormack1995)

(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. Reference Comeaux, Chapman and MacCormack1995). 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 Reference Müller1972) 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 Reference Müller1972; Söderholm Reference Söderholm1976). It was remarked by Truesdell (Reference Truesdell1976) 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 Reference Liu and Müller1986), 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 (Reference Müller1972) 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.