2. Velocity gradient dynamics in the presence of vibrational non-equilibrium

In this section, we derive the evolution equation of the velocity gradient ($\boldsymbol {A}$) and the pressure-Hessian ($\boldsymbol {P}$) tensors under vibrational non-equilibrium conditions. Mathematically, these tensors are defined as

(2.1a,b)\begin{equation} A_{ij}=\frac{\partial u_i}{\partial x_j} \quad \mathrm{and} \quad P_{ij}=\frac{\partial}{\partial x_j}\left(\frac{1}{\rho} \frac{\partial p}{\partial x_i}\right), \end{equation}

where $u_i$ and $x_i$ represent the $i$th velocity component and $i$th spatial coordinate in the Eulerian description of the flow field. The symbol $\rho$ represents density, and $p$ is pressure.

We start with the governing equations of mass and momentum for a viscous compressible medium

(2.2)$$\begin{gather} \frac{\partial\rho}{\partial t}+u_j\frac{\partial\rho}{\partial x_j}={-}\rho\frac{\partial u_j}{\partial x_j}, \end{gather}$$

(2.3)$$\begin{gather}\frac{\partial u_i}{\partial t}+u_j\frac{\partial u_i}{\partial x_j}={-}\frac{1}{\rho}\frac{\partial p}{\partial x_i}+\frac{1}{\rho}\frac{\partial\tau_{ji}}{\partial x_j}. \end{gather}$$

The viscous stress $(\tau _{ji})$ is expressed using the following relationship:

(2.4)\begin{equation} \tau_{ji}= 2 \mu\left[\frac{1}{2}\left(\frac{\partial u_i}{\partial x_j}+ \frac{\partial u_j}{\partial x_i}\right)\right]+\delta_{ji}\zeta\frac{\partial u_k}{\partial x_k}, \end{equation}

where, $\delta _{ji}$ is the Kronecker delta symbol, $\mu$ denotes dynamic viscosity and $\zeta$ is the bulk viscosity of the fluid medium. Bulk viscosity ($\zeta$) is taken to be equal to $-(2/3)\mu$ (White Reference White1979).

The governing equation of the total energy per unit mass $(e)$ is

(2.5)\begin{equation} \frac{\partial e}{\partial t}+ u_i\frac{\partial e}{\partial x_i} ={-}\frac{1}{\rho} \frac{\partial}{\partial x_i}(u_i p)+\frac{1}{\rho} \frac{\partial}{\partial x_j}(u_i\tau_{ij})-\frac{1}{\rho} \frac{\partial q_i}{\partial x_i}-\frac{1}{\rho}\frac{\partial q_{v,i}}{\partial x_i}, \end{equation}

where $e$ represents the total energy (per unit mass), which equals the sum of the kinetic energy, translational–rotational energy and vibrational energy ($e_v$) per unit mass. It is expressed as

(2.6)\begin{equation} e=\underbrace{\frac{u_iu_i}{2}}_{Kinetic\ energy}+\underbrace{\frac{5}{2}RT}_{Translational\text{-}rotational\ energy}+\underbrace{e_v}_{Vibrational\ energy}. \end{equation}

The symbol $T$ is the temperature of the gas, and $R$ is the gas constant (J kg$^{-1}$ K$^{-1}$). The symbols $q_i$ and $q_{v,i}$ are the thermal and vibrational heat fluxes (J s$^{-1}$ m$^{-2}$) in the $i$th direction.

In order to understand the influence of vibrational non-equilibrium on compressible flows, an additional equation of the evolution of vibrational energy per unit mass is introduced

(2.7)\begin{equation} \frac{\partial e_v}{\partial t}+u_i\frac{\partial e_v}{\partial x_i} ={-}\frac{1}{\rho} \frac{\partial q_{v,i}}{\partial x_i}+\frac{w_v}{\rho}. \end{equation}

This equation uses a source term ($w_v$) that governs the energy transfer between translational–rotational and vibrational modes. We employ the Landau–Teller model for the source term $(w_v)$ expressed as (Vincenti & Kruger Reference Vincenti and Kruger1982)

(2.8)\begin{equation} w_v ={-} \rho \frac{(e_{v}-e_{v}^*)}{\tau}, \end{equation}

where $e_{v}^*$ is the vibrational energy per unit mass in equilibrium with local temperature $(T)$

(2.9)\begin{equation} e_{v}^* = \frac{\dfrac{h\nu}{k_BT}}{\exp\left({\dfrac{h\nu}{k_B T}}\right)-1}RT. \end{equation}

The symbol $h$ represents Plank's constant, $\nu$ is the vibrational frequency of molecules of the gas and $k_B$ is the Boltzmann constant. The quantity $\tau$ appearing in (2.8) is the vibrational relaxation time (Park Reference Park1993)

(2.10)\begin{equation} \tau = \frac{1}{p} \exp[A (T^{{-}1/3} - B) - 18.42], \end{equation}

where $p$ is in atm. Coefficients $A$ and $B$ are calculated using (2.11) and (2.12), respectively,

(2.11)$$\begin{gather} A = 1.16 \times 10^{{-}3} \sqrt{M} \theta^{4/3}, \end{gather}$$

(2.12)$$\begin{gather}B = 0.015 M^{1/4}. \end{gather}$$

The symbol $M$ is the molecular weight of the gas in gram/mol, and $\theta$ is the characteristic vibrational temperature of the gas in $K$. The values of $A$ and $B$ for nitrogen $(N_2)$ are 221.46 and 0.029, respectively (Park Reference Park1993). In line with (2.9), we can express $e_{v}$ also as a function of new variable $T_v$, the so-called local vibrational temperature

(2.13)\begin{equation} e_{v} = \frac{\dfrac{h\nu}{k_BT_v}}{\exp\left({\dfrac{h\nu}{k_B T_v}}\right)-1}RT_v. \end{equation}

Finally, the thermodynamic variables like pressure ($p$), density ($\rho$) and temperature ($T$) are related algebraically through the state equation of a perfect gas

(2.14)\begin{equation} p =\rho RT. \end{equation}

Equation set (2.2)–(2.14) describes a compressible flow field under the vibrational non-equilibrium condition with density ($\rho$), pressure ($p$), temperature ($T$), velocity ($u_i$) and vibrational energy ($e_v$) as seven dependent variables.

To derive the exact evolution equation of the velocity gradient tensor $\boldsymbol {A}$ under a vibrational non-equilibrium condition, we take the gradient of (2.3) to arrive at the following equation:

(2.15)\begin{equation} \frac{{\rm D} A_{ij}}{{\rm D} t}={-A_{ik}A_{kj}}-{P_{ij}}+{\varGamma_{ij}}. \end{equation}

The symbol, ${{\rm D}}/{{\rm D} t}={\partial }/{\partial t}+u_{k}({\partial }/{\partial x_{k}})$ in (2.15) is the material derivative representing the change in $A_{ij}$ following the motion of the local small fluid element. The first term on the right-hand side is the self-stretching term, representing the inertial action on the local state of the velocity gradient tensor. The second term is the pressure-Hessian tensor ($\boldsymbol {P}$). The pressure-Hessian tensor represents the process by which density and pressure gradients directly influence the evolution of velocity gradients in compressible flows. The third term is the viscous transport term

(2.16)\begin{equation} \varGamma_{ij}=\frac{\partial}{\partial x_j}\left\{\frac{1}{\rho} \frac{\partial}{\partial x_k}\left[\mu\left(A_{ik} +A_{ki}- \frac{2}{3}A_{pp}\delta_{ki}\right)\right]\right\}. \end{equation}

Our interest is in having a simple dynamical equation set to represent the dynamics of velocity gradients. While (2.15) describes the exact evolution of the velocity gradient tensor ($\boldsymbol {A}$), it is mathematically unclosed. The unclosed terms are $P_{ij}$ and $\varGamma _{ij}$. In incompressible flows, pressure is completely governed by Poisson's pressure equation. However, in compressible flows, the acoustic time scale responsible for pressure evolution may not be significantly smaller than the fluid time scale, and thus pressure behaves like a $bona fide$ thermodynamic variable. Therefore, in compressible flows, it is desirable to have an evolution equation of the pressure-Hessian tensor ($\boldsymbol {P}$) as well. Suman & Girimaji (Reference Suman and Girimaji2011) derived the exact evolution equation of the pressure-Hessian tensor using the continuity, energy and the state equations for compressible flows with no reference to the vibrational excitation of the constituent molecules. Consequently, their evolution equation of the $\boldsymbol {P}$ tensor does not account for the effects of vibrational non-equilibrium in a flow field. In contrast, in this work, we aim to derive the exact evolution equation of the pressure Hessian for a compressible flow under the influence of the vibrational non-equilibrium process.

To derive such an evolution equation of the pressure Hessian, we first express the material derivative of the pressure Hessian as

(2.17)\begin{align} \frac{{\rm D} P_{ij}}{{\rm D} t} &= \frac{{\rm D}}{{\rm D} t} \left\{\frac{\partial}{\partial x_j}\left(\frac{1}{\rho}\frac{\partial p}{\partial x_{i}}\right)\right\} \nonumber\\ &= \frac{2}{\rho^3}\frac{{\rm D}\rho}{{\rm D} t} \frac{\partial\rho}{\partial x_j} \frac{\partial p}{\partial x_i} - \frac{1}{\rho^2} \frac{\partial p}{\partial x_i}\frac{{\rm D}}{{\rm D} t} \left(\frac{\partial\rho}{\partial x_j}\right)-\frac{1}{\rho^2} \frac{\partial\rho}{\partial x_j}\frac{{\rm D}}{{\rm D} t}\left(\frac{\partial p}{\partial x_i}\right) \nonumber\\ &\quad -\frac{1}{\rho^2} \frac{{\rm D}\rho}{{\rm D} t}\left(\frac{\partial^2p}{\partial x_j\partial x_i}\right) + \frac{1}{\rho}\frac{{\rm D}}{{\rm D} t}\left(\frac{\partial^2p}{\partial x_{j}\partial {x_i}}\right). \end{align}

Using (2.14) we write

(2.18)\begin{equation} \frac{{\rm D} p}{{\rm D} t}=\frac{{\rm D}\rho}{{\rm D} t}RT+\rho R\frac{{\rm D} T}{{\rm D} t}. \end{equation}

Next, we derive the material derivative of temperature. Equation (2.5) is the evolution equation of total energy $[e={u_iu_i}/{2}+\frac {5}{2}RT+e_v]$. The governing equation for kinetic energy is obtained by multiplying (2.3) with $u_i$. Subsequently, subtracting the evolution equation of kinetic energy and that of vibrational energy (2.7) from the total energy equation (2.5) leads us to the evolution equation of the temperature $(T)$

(2.19)\begin{equation} \frac{\partial T}{\partial t}+u_i \frac{\partial T}{\partial x_i}={-}\frac{2}{5}T \frac{\partial u_i}{\partial x_i}+\frac{2}{5}\frac{\tau_{ij}}{\rho R} \frac{\partial u_i}{\partial x_j}+\frac{2}{5}\frac{k}{\rho R} \frac{\partial}{\partial x_k}\left(\frac{\partial T}{\partial x_k}\right)-\frac{2}{5}\frac{w_v}{\rho R}, \end{equation}

where $k$ is the thermal conductivity of the medium (J s$^{-1}$ m$^{-1}$ K$^{-1}$), and it appears when the thermal flux $q_i$ is represented using Fourier's law of heat conduction: ($q_i=-k({\partial {T}}/{\partial {x}_i})$). Now, using (2.19) in (2.18) along with (2.2) leads to the evolution equation of pressure under the vibrational non-equilibrium conditions

(2.20)\begin{equation} \frac{\partial p}{\partial t}+u_i \frac{\partial p}{\partial x_i}={-}\frac{7}{5}p \frac{\partial u_i}{\partial x_i}+\frac{2}{5}\tau_{ij}\frac{\partial u_i}{\partial x_j}+ \frac{2}{5}\frac{k}{R}\frac{\partial}{\partial x_k}\left\{\frac{\partial}{\partial x_k} \left(\frac{p}{\rho}\right)\right\}-\frac{2}{5} w_v. \end{equation}

We subject (2.20) to the gradient operator successively to arrive at the equation $({{\rm D}}/{{\rm D} t})({\partial p}/{\partial x_i})$ and $({{\rm D}}/{{\rm D} t})({\partial ^2p}/{\partial x_i\partial x_j})$. Subsequently, employing (2.17), we finally arrive at the exact evolution equation of the pressure-Hessian tensor in the presence of vibrational non-equilibrium

(2.21)\begin{align} \frac{{\rm D} P_{ij}}{{\rm D} t} &= {-A_{kj}P_{ik}-A_{ki}P_{kj}-\frac{2}{5}A_{kk}P_{ij}} \nonumber\\ &\quad -\frac{7}{5\rho}\frac{\partial A_{kk}}{\partial x_j}\frac{\partial p}{\partial x_i}- \frac{7}{5\rho}\frac{\partial A_{kk}}{\partial x_i}\frac{\partial p}{\partial x_j}-\frac{1}{\rho} \frac{\partial A_{ki}}{\partial x_j}\frac{\partial p}{\partial x_k}+\frac{1}{\rho} \frac{\partial A_{kk}}{\partial x_j}\frac{\partial p}{\partial x_i} \nonumber\\ &\quad -\frac{7}{5}\frac{p}{\rho}\left(\frac{\partial^2A_{kk}}{\partial x_i\partial x_j}-\frac{1}{\rho} \frac{\partial A_{kk}}{\partial x_i}\frac{\partial\rho}{\partial x_j}\right) +\left(\frac{2}{5R}\right)\frac{\partial}{\partial x_j}\left\{\frac{1}{\rho} \frac{\partial^2}{\partial x_i\partial x_k}\left[k \frac{\partial}{\partial x_k} \left(\frac{p}{\rho}\right)\right]\right\} \nonumber\\ &\quad +\frac{2}{5}\frac{\partial}{\partial x_j}\left[\frac{1}{\rho}\frac{\partial\tau_{km}}{\partial x_i}A_{km}+\frac{1}{\rho}\tau_{km}\frac{\partial A_{km}}{\partial x_i}\right] - \frac{2}{5} \frac{\partial}{\partial{x}_j}\left(\frac{1}{\rho}\frac{\partial{w}_v}{\partial{x}_i}\right). \end{align}

To make further progress, we must have a closure model for $w_v$. Employing the Landau–Teller model in (2.21), the pressure-Hessian equation modifies to

(2.22)\begin{align} \frac{{\rm D} P_{ij}}{{\rm D} t} &= \underbrace{-A_{kj}P_{ik}-A_{ki}P_{kj}-\frac{2}{5}A_{kk}P_{ij}}_{(I)} \nonumber\\ &\quad \underbrace{-\frac{7}{5\rho}\frac{\partial A_{kk}}{\partial x_j} \frac{\partial p}{\partial x_i}-\frac{7}{5\rho}\frac{\partial A_{kk}}{\partial x_i} \frac{\partial p}{\partial x_j}-\frac{1}{\rho}\frac{\partial A_{ki}}{\partial x_j} \frac{\partial p}{\partial x_k}+\frac{1}{\rho}\frac{\partial A_{kk}}{\partial x_j} \frac{\partial p}{\partial x_i}}_{(II)} \nonumber\\ &\quad \underbrace{-\frac{7}{5}\frac{p}{\rho}\left(\frac{\partial^2A_{kk}}{\partial x_i\partial x_j}- \frac{1}{\rho}\frac{\partial A_{kk}}{\partial x_i}\frac{\partial\rho}{\partial x_j}\right)}_{(III)} \underbrace{+\left(\frac{2}{5R}\right)\frac{\partial}{\partial x_j}\left\{\frac{1}{\rho} \frac{\partial^2}{\partial x_i\partial x_k}\left[k\frac{\partial}{\partial x_k} \left(\frac{p}{\rho}\right)\right]\right\}}_{(IV)} \nonumber\\ &\quad +\underbrace{\frac{2}{5}\frac{\partial}{\partial x_j}\left[\frac{1}{\rho} \frac{\partial\tau_{km}}{\partial x_i}A_{km}+\frac{1}{\rho}\tau_{km} \frac{\partial A_{km}}{\partial x_i}\right]}_{(V)} \nonumber\\ &\quad +\underbrace{\frac{2}{5}\frac{\partial}{\partial x_j} \left[\frac{1}{\rho}\frac{\partial}{\partial x_i}\left(\frac{\rho{e}_v}{\tau}\right)\right]}_{(VI)} - \underbrace{\frac{2}{5}\frac{\partial}{\partial x_j}\left[\frac{1}{\rho} \frac{\partial}{\partial x_i}\left(\frac{\rho{e}_v^{*}}{\tau}\right)\right]}_{(VII)}, \end{align}

where $e_v^*$ and $e_v$ are defined in (2.9) and (2.13), respectively.

To make this equation mathematically as well as physically more amenable to our understanding, we have gathered various terms as groups $(I-VII)$. We refer to these groups as mechanisms that influence the evolution of the pressure-Hessian tensor (in line with Suman & Girimaji Reference Suman and Girimaji2011).

The first mechanism on the right-hand side involves interactions between our primary unknowns, the $\boldsymbol {A}$ and $\boldsymbol {P}$ tensors only. In line with the nomenclature for the $-A_{ik}A_{kj}$ term in (2.15), this mechanism is called the local stretching of the pressure Hessian field by the local velocity gradient field (Suman & Girimaji Reference Suman and Girimaji2011).

The second mechanism involves higher-order gradients of the velocity gradient components and the first gradient of the pressure field. In contrast, mechanism $III$ involves terms in which higher-order gradients of $\boldsymbol {A}$ interact with the raw thermodynamic field $p/\rho$ itself. Suman & Girimaji (Reference Suman and Girimaji2011) argued that the presence of the raw thermodynamic variable $p/\rho$ (rather than its gradients) makes the underlying physics dependent on an acoustic time scale (and, accordingly, a Mach number). Indeed, the authors subsequently modelled the essential physics incumbent in this mechanism to recover the well-known Mach number effects from their enhanced homogenized Euler equation model of the compressible velocity gradient tensor. Mechanism $IV$ involves thermal conductivity $(k)$. This clearly represents the influence of thermal conduction on the pressure-Hessian evolution. On the other hand, mechanism $V$ involves the effect of viscous heating on the evolution of the $\boldsymbol {P}$ tensor. Mechanisms $I$ to $V$ were discussed by Suman & Girimaji (Reference Suman and Girimaji2011) as well. However, the authors examined a compressible flow field in the absence of any vibrational excitation. Their total energy equation involved changes in the kinetic energy and translational–rotational energy only. In contrast to mechanisms $I$ to $V$, mechanisms $VI$ and $VII$ are the ‘new’ mechanisms appearing in the evolution equation in the presence of vibrational non-equilibrium. These two ‘new’ mechanisms show the explicit influence of vibrational energy on the evolution of the dynamics of the pressure-Hessian field and, in turn, that on the velocity gradient field, as well. These mechanisms are attributable to the employment of the Landau–Teller model (2.8) for the interconversion of translational–rotational and vibrational energy modes. In subsequent sections, we wish to understand how these mechanisms affect the pressure-Hessian and the velocity gradient tensors.

3. The closure problem

We now examine (2.22) from the point of view of mathematical closure. This examination is motivated by the eventual goal of having a computationally viable set of ordinary differential equations for the $\boldsymbol {A}$ and $\boldsymbol {P}$ tensors. Among mechanisms $I$ to $V$, all are mathematically unclosed except for mechanism $I$. Suman & Girimaji (Reference Suman and Girimaji2011) and Danish et al. (Reference Danish, Suman and Srinivasan2014) have already addressed this closure problem and proposed some models. In this paper, we focus on the ‘new’ mechanisms $VI$ and $VII$. These mechanisms are also mathematically unclosed for multiple reasons: the algebraic involvement of higher-order gradients of the velocity gradients, density, pressure, temperature and the gradients of vibrational temperature.

Danish et al. (Reference Danish, Suman and Srinivasan2014) has earlier employed a DNS database of isotropic compressible turbulence to examine the relative importance of mechanisms $I$ to $V$ of the evolution equation of the $\boldsymbol {P}$ tensor for a compressible medium without any vibrational excitation. To identify the dominant mechanisms, they contracted the evolution equation of $P_{ij}$ and arrived at the evolution equation of $({{\rm D}}/{{\rm D} t})(P_{ij}P_{ij})$. All the original mechanisms were still present on the right-hand side of that equation; however, they appeared in their contracted scalar forms. The relative significance of these scalar forms was then compared by plotting the PDFs of the relative magnitudes of these mechanisms using a DNS database. In line with this approach, we, too, first contract equation (2.22) with ${Pij}$ and arrive at the scalar version of (2.22)

(3.1) \begin{align} \frac{1}{2} \frac{{\rm D}}{{\rm D} t}(P_{ij}P_{ij}) &= \underbrace{\left({-}A_{kj}P_{ik}-A_{ki}P_{kj}-\frac{2}{5}A_{kk}P_{ij}\right)P_{ij}}_{(I')} \nonumber\\ &\quad -\underbrace{\left(\frac{7}{5\rho}\frac{\partial A_{kk}}{\partial x_j} \frac{\partial p}{\partial x_i}+\frac{7}{5\rho}\frac{\partial A_{kk}}{\partial x_i} \frac{\partial p}{\partial x_j}+\frac{1}{\rho}\frac{\partial A_{ki}}{\partial x_j} \frac{\partial p}{\partial x_k}-\frac{1}{\rho}\frac{\partial A_{kk}}{\partial x_j} \frac{\partial p}{\partial x_i}\right)P_{ij}}_{(II')} \nonumber\\ &\quad -\underbrace{\frac{7}{5}\frac{p}{\rho}\left(\frac{\partial^2A_{kk}}{\partial x_i\partial x_j}- \frac{1}{\rho}\frac{\partial A_{kk}}{\partial x_i} \frac{\partial\rho}{\partial x_j}\right)P_{ij}}_{(III')}\nonumber\\ &\quad + \underbrace{\left(\frac{2}{5R}\right)\ \frac{\partial}{\partial x_j} \left[\frac{1}{\rho}\frac{\partial^2}{\partial x_i\partial x_k} \left\{k\frac{\partial}{\partial x_k}\left(\frac{p}{\rho}\right)\right\}\right]P_{ij}}_{(IV')} \nonumber\\ &\quad + \underbrace{\left(\frac{2}{5}\ \frac{\partial}{\partial x_j}\left[\frac{1}{\rho} \frac{\partial\tau_{km}}{\partial x_i}A_{km}+\frac{1}{\rho}\tau_{km} \frac{\partial A_{km}}{\partial x_i}\right]\right)P_{ij}}_{(V')} \nonumber\\ &\quad + \underbrace{\frac{2}{5}\frac{\partial}{\partial x_j}\left[\frac{1}{\rho} \frac{\partial}{\partial x_i}\left(\frac{\rho{e}_v}{\tau}\right)\right]P_{ij}}_{(VI')} - \underbrace{\frac{2}{5}\frac{\partial}{\partial x_j}\left[\frac{1}{\rho} \frac{\partial}{\partial x_i}\left(\frac{\rho{e}_v^{*}}{\tau}\right)\right]P_{ij}}_{(VII')}. \end{align}

In this equation, mechanisms $VI'$ and $VII'$ are the scalar counterparts of mechanisms $VI$ and $VII$ of (2.22). On the other hand, mechanisms $I'$ to $V'$ are the scalar counterparts of the non-vibration mechanisms $I$ to $V$ of (2.22). We intend to examine the relative importance of the collective action of the vibrational mechanisms $VI'$–$VII'$ to that of the inertial mechanisms $I', II', III'$. At this point, the physics related to the viscous processes ($IV'$ and $V'$) are not of interest. Thus, for a direct and unambiguous study, we define a fraction $f$

(3.2)\begin{equation} f = \frac{|VI' - VII'|}{| I' + II' + III'| + | VI' - VII'|}. \end{equation}

Examining $f$ can be useful to quantify the extent to which the evolution of the pressure Hessian is affected by the vibrational non-equilibrium process, if at all. The natural algebraic range of $f$ is $[0,1]$. A small (near-zero value) of $f$ would imply that the contribution of the vibrational mechanism is negligible in the governing equation of the pressure Hessian. On the other hand, a high magnitude of $f$ would imply a substantial contribution of the vibrational mechanism in governing and shaping the pressure-Hessian tensor in a compressible flow field.

4. Direct numerical simulations of isotropic turbulence with vibrational non-equilibrium

We have employed the hy2Foam solver of Casseau et al. (Reference Casseau, Palharini, Scanlon and Brown2016) and Casseau (Reference Casseau2017) to perform DNS for our study. hy2foam is an open-source (Git 2022) solver that has been developed on the OpenFOAM platform and has the capability to handle both vibrational and chemical non-equilibrium. Even though the solver offers various turbulence models to perform Favre-averaged simulations, no turbulence model has been used to generate a well-resolved direct numerical simulation database. In this case, the solver solves the instantaneous continuity equation (2.2), velocity equation (2.3) and total energy equation (2.5). The state equation of a perfect gas is used as given by (2.14). Further, the evolution equation of vibrational energy ($e_v$) (2.7) is solved using the Landau–Teller model as presented in (2.8). The vibrational relaxation time scale appearing in (2.8) is calculated using the Millikan–White equation, as detailed in (2.10).

The computational domain is a cubical box with each side of length $2{\rm \pi}$. The computational grid is uniform, with the number of nodes being $192^3$/$256^3$/$512^3$. Periodic boundary conditions are imposed on the opposite sides of the cubic domain. A random initial velocity field is generated with zero mean and zero divergence. The initial spectrum of turbulent kinetic energy $E(\kappa )$ is

(4.1)\begin{equation} E(\kappa,t=0) = A_0\kappa^4\exp({-}2\kappa^2/\kappa^2_0), \end{equation}

where $A_0$ represents the spectrum constants. The symbol $\kappa$ is the magnitude of the wavenumber vector, and $\kappa _0$ is the $\kappa$ where $E(k)$ is maximum at $t=0$.

In homogeneous isotropic decaying turbulence, the Reynolds number based on the Taylor microscale $(Re_{\lambda })$ and initial turbulent Mach number ($M_t$) are defined as

(4.2a,b)\begin{equation} Re_{\lambda}=\sqrt{\frac{20}{3\eta\epsilon}}K \quad \text{and}\quad M_t=\sqrt{\frac{ 2K }{n R\langle T \rangle}}, \end{equation}

where $\eta$, $\epsilon$ and $K$ represent kinematic viscosity, dissipation rate and turbulent kinetic energy, respectively, and $\langle T \rangle$ is the volume-averaged temperature. The symbol $n$ represents the ratio of the specific heat capacity at constant pressure to that at constant volume.

In addition to $M_t$ and $Re_{\lambda }$, two more non-dimensional numbers characterize a simulation in the presence of vibrational non-equilibrium: (i) the initial Damköhler number $(D)$, and the initial ratio of the mean vibrational temperature to mean temperature $(\xi )$. These numbers are defined as in (4.3) and (4.4), respectively,

(4.3)$$\begin{gather} D = \frac{\tau_{LE}}{\tau}, \end{gather}$$

(4.4)$$\begin{gather}\xi=\frac{\langle{T_v}\rangle}{\langle{T}\rangle}, \end{gather}$$

where $\tau _{LE}$ is the large eddy-turnover time defined as $\tau _{LE} = \lambda /u^{\prime }$. The root mean squared value of the fluctuating velocity, denoted as $u'$, is defined as

(4.5)\begin{equation} u' = \sqrt{\left\langle \frac{u_1^2 + u_2^2 + u_3^2}{3} \right\rangle}, \end{equation}

where $\langle \text { }\rangle$ represents the volume-averaged value of the quantity of interest (equivalently, the mean of the quantity), and $\lambda$ is the Taylor microscale, equal to (Honein & Moin Reference Honein and Moin2004)

(4.6)\begin{equation} \lambda = \sqrt{\frac{u'^2}{\langle(\partial u_1/\partial x_1)^2 \rangle}}. \end{equation}

The definition of the Damköhler number ($D$) used in this study is defined in line with the similar dimensionless number used in chemically reacting flows to relate the characteristic time scale of chemical reactions to that of the advective transport phenomenon (or the flow time scale) rate occurring therein (Fogler Reference Fogler2006). As the initial Damköhler number becomes large, the disparity between the fluid time scale and the vibrational relaxation time scale increases, and an initial vibrational non-equilibrium state will more rapidly tend to approach its equilibrium state. In the limit of the Damköhler number being zero, the vibrational state of the flow will remain effectively frozen even though the flow itself evolves. In this study, our objective is to examine the influence of this parameter having a value in the moderate range of 0.25 to 4, wherein both turbulence and vibrational relaxation processes are expected to evolve substantially.

The range of $\xi$ used in our simulation set spans from 0.125 to 0.75. In hypersonic aerospace applications, a shock upstream of a solid obstruction typically causes a substantial disparity in the vibrational temperature values upstream and downstream of the shock front. The translational energy rapidly attains its Maxwellian equilibrium state immediately downstream of the shock. This increases the temperature immediately downstream of the shock front. However, the vibrational temperature after the shock might remain at its pre-shock value. Downstream of such a shock front, the ratio of vibrational temperature to the local temperature is less than unity. Indeed, our chosen range for $\xi$ reflects such typical and realistic scenarios that initiate a vibrational relaxation process alongside the evolution of the background flow field.

In table 1, we present the list of simulation cases employed in our study. Simulation A has identical initial conditions as employed by Honein & Moin (Reference Honein and Moin2004) and is used to validate the hy2FOAM solver against the results of Honein & Moin (Reference Honein and Moin2004). Like Honein & Moin (Reference Honein and Moin2004), in simulation A, the viscosity dependence on temperature is modelled using a power law with the exponent being $0.76$. This simulation does not involve any vibrational non-equilibrium process. On the other hand, simulation cases B1–B4 and C–F all have the vibrational relaxation process included in the simulation. Simulation cases B1–B4 have been performed for the most ‘severe’ combination of the initial parameters employed in our study. For all these simulations, the turbulent Mach number, the Reynolds number, the Damköhler number and $\xi$ have been set as 0.5, 60, 4 and 0.125, respectively. While B1–B3 differ only in terms of the grid size ($192^3, 256^3$ and $512^3$), case B4 differs from B2 only in terms of the Courant–Friedrich–Lewy (CFL) number.

Table 1. The DNS employed in this study.

Cases B2, C and D are employed to evaluate the influence of the initial Damköhler number on the physics of interest since these sets of simulations differ only in terms of that parameter. The initial turbulent Mach number and the Taylor microscale-based Reynolds number are 0.5 and 60, respectively, for these simulations. The initial value of $\xi$ equals 0.5. On the other hand, cases B2, E and F differ in terms of $\xi$. Simulation case G has its settings identical to B2, but has the vibrational relaxation processes switched off. This simulation is used as the baseline case with which the results from the simulation cases involving vibrational relaxation processes are compared.

In all these simulations (B1–B4 and C–G) the gaseous medium is chosen to be nitrogen $(N_2)$. In these simulations, we have kept the coefficient of viscosity ($\mu$) a constant (as determined by the initial Reynolds number). This is a deliberate choice because we wish to isolate and identify the direct influence of the variation in the Damköhler number and $\xi$ on the velocity gradient and the pressure-Hessian dynamics. We wish to mention that our simulations are different from those of Khurshid & Donzis (Reference Khurshid and Donzis2019), who not only use a temperature-dependent viscosity coefficient, but also allow the flow field to decay from an initial state which has been extracted from an already fully developed state of forced turbulence. We plan to include a variable viscosity coefficient in our future work.

Subject to the availability of the library of discretization schemes in the hy2Foam package, we have employed the second-order accurate central-upwind scheme developed by Kurganov, Noelle and Petrova (Kurganov, Noelle & Petrova Reference Kurganov, Noelle and Petrova2001; Greenshields et al. Reference Greenshields, Weller, Gasparini and Reese2010) for the advection term, whereas the diffusive fluxes utilize the central difference method. For temporal discretization, a second-order accurate, backward implicit formulation is used. Uniform time stepping is used in all the simulations. The interpolation schemes employed are of second-order and unbounded. The value of the time step is chosen to be such that, initially, the CFL number is in the range of 0.00625–0.025 (case-specific values are listed in table 1).

Before discussing the results related to the vibrational non-equilibrium phenomenon, we first assess the capability of the hy2Foam solver in simulating standard compressible isotropic turbulence (case A). We compare the results from this simulation against the results of Honein & Moin (Reference Honein and Moin2004). In figure 1(a–c), we present comparisons in terms of the temporal evolution of appropriately normalized turbulence kinetic energy ($K$), root-mean-square (r.m.s.) pressure fluctuations ($p'$) and the r.m.s. of the temperature ($T'$). The normalizing quantities in figure 1(a–c) are the initial value of the turbulent kinetic energy, the initial values of $(np_oM_t^2)$ and the initial values of $(n-1)(T_oM_t^2)$, respectively. The subscript $o$ denotes the initial value of the respective quantity. In each case, the time axis is normalized by the eddy-turnover time ($\tau _{LE} = {\lambda }/{u'}$) calculated based on the initial flow field. The initial condition set for these simulations corresponds to Honein & Moin (Reference Honein and Moin2004), wherein all thermodynamic fluctuations are zero initially. This choice explains the steep variations in the thermodynamic quantities during the first few time steps of the simulations. In figure 1(d), we present the evolution of skewness ($S_u$) from case A

(4.7)\begin{equation} S_u=\frac{\left\langle \left[ \left( \dfrac{\partial u_1}{\partial x_1} \right)^3 + \left(\dfrac{\partial u_2}{\partial{x}_2} \right)^3 + \left(\dfrac{\partial u_3}{\partial{x}_3} \right)^3 \right] / 3 \right\rangle}{\left\langle \left[ \left( \dfrac{\partial u_1}{\partial{x}_1} \right)^2 + \left(\dfrac{\partial u_2}{\partial{x}_2} \right)^2 + \left(\dfrac{\partial u_3}{\partial{x}_3} \right)^2 \right] / 3 \right\rangle^{3/2}}. \end{equation}

The included reference results in figure 1(a–d) are the digitized versions of the data available in Honein & Moin (Reference Honein and Moin2004) and Chen et al. (Reference Chen, Wen, Wang and Guo2020) (for $S_u$). Evidently, the results of the case A simulation agree quite well with the reference DNS results.

Figure 1. Simulation case A. Temporal evolution of (a) turbulent kinetic energy, (b) pressure fluctuations, (c) temperature fluctuations and (d) skewness. Symbols represent the results of Honein & Moin (Reference Honein and Moin2004) and Chen et al. (Reference Chen, Wen, Wang and Guo2020).

In figures 2(a–c) and 3(a–c) we present results from simulation cases B1–B3. Figure 2(a–c) shows the temporal evolution of the turbulent kinetic energy and the thermodynamic fluctuations, whereas figure 3(a–c) shows the temporal evolution involving different aspects of the velocity gradient field: skewness ($S_u$), r.m.s. of vorticity ($\omega ^{'}$) and r.m.s. of dilatation rate ($\theta ^{'}$)

(4.8)$$\begin{gather} \omega^{'}={\langle \omega_i \omega_i \rangle}^{1/2}, \end{gather}$$

(4.9)$$\begin{gather}\theta'=\left\langle \left( \frac{\partial u_i}{\partial x_i} \right)^2 \right\rangle^{1/2}, \end{gather}$$

where $\omega _i$ represents the $i$th component of the instantaneous vorticity vector. All these time-dependent statistics have been normalized by their respective initial values. The difference between the results obtained on the $256^3$ grid and those obtained on the $512^3$ grid is small. The maximum error percentage of $\omega ^{'}$ between the $256^3$ grid and the $512^3$ grid is less than $3.5\,\%$ over the entire duration of $4$ eddy-turnover times. Similarly, the maximum error percentage of $\theta ^{'}$ between the $256^3$ grid and the $512^3$ grid is less than $3.0\,\%$ over the entire duration of $4$ eddy-turnover times. Based on these results, we deem the results of simulation B2 ($256^3$) to be adequately grid independent. Thus, other simulations (cases C–G) have also been performed on a $256^3$ grid. We have performed a similar grid-independence study corresponding to the initial conditions of the simulation case G, as well. Even though those results are not included in the paper, the $256^3$ grid is found to be adequate there too.

Figure 2. Grid independence. Temporal evolution of (a) turbulent kinetic energy, (b) r.m.s. of density, (c) r.m.s. of temperature fluctuations in simulation cases B1, B2 and B3.

Figure 3. Grid independence: (a) skewness, (b) vorticity fluctuations and (c) dilatational fluctuations in simulation cases B1, B2 and B3.

We have performed a time-step-independence study using the results of the simulation cases B2 and B4. In simulation B4, the employed time step is twice that of the time step used in B2. The corresponding difference in the results (turbulent kinetic energy, $p^{'}$, $T^{'}$, $S_u$, $\omega ^{'}$ and $\theta ^{'}$) from these two simulations is found to be insignificant (results not included here to avoid repetition). Thus, the time step employed in all other simulations of table 1 is chosen to be the same as employed in the simulation case B2.

In figure 4, we present the evolution of skewness for simulations B2, C, D, E, F and G (4.7). We observe that, beyond $t/\tau _{LE}>1.5$, in every case, the skewness stabilizes to a value close to $-$0.45, indicating that, beyond that time, the nonlinear turbulent processes are in full effect. In figure 5 we present the decay of $Re_{\lambda }$ and turbulent Mach number ($M_t$, 4.2a,b) in various simulations that differ in terms of initial Damköhler number and the initial value of $\xi$. Evidently, there is no influence of these initial parameters on the decay of $Re_\lambda$ and $M_t$.

Figure 4. Temporal evolution of skewness for simulation cases (a) B2, C, D and G and (b) B2, E, F and G.

Figure 5. Temporal evolution of (a) Taylor micro-scale Reynolds number and (b) turbulent Mach number in simulation cases B2, C, D, E and F.

5. The vibrational relaxation process and plan of study

To visualize the vibrational relaxation process in our simulations, we first present the temporal evolution of the statistics of temperature and the vibrational temperature. In figure 6(a) we present the temporal evolutions of mean temperatures ($\langle T \rangle$ and $\langle T_v \rangle$) in the sets of simulations B2, C and D, which differ in terms of the initial Damköhler number. In figure 6(b) we present the temporal evolutions of the r.m.s. of the two temperature quantities ($\sqrt {\langle T^{'}T^{'} \rangle } \text { and } \sqrt {\langle T_v^{'}T_v^{'} \rangle }$) in the same set of simulations. In figure 7(a) we present the temporal evolutions of mean temperatures ($\langle T \rangle$ and $\langle T_v \rangle$) in the sets of simulations B2, E and F, which differ in terms of the initial $\xi$. In figure 7(b) we present the temporal evolutions of r.m.s. of fluctuations in the two temperature quantities ($\sqrt {\langle T^{'}T^{'} \rangle } \text { and } \sqrt {\langle T_v^{'}T_v^{'} \rangle }$) in the same set of simulations. The mean temperature profiles have been non-dimensionalized with $T_o$, whereas the r.m.s. values of fluctuations have been normalized with the factor $(n-1)T_o{M}_t^2$.

Figure 6. (a) Evolution of mean temperature (with symbols) and mean vibrational temperature (without symbols) in simulation cases B2, C and D. (b) Evolution of r.m.s. of temperature (with symbols) and that of the vibrational temperature (without symbols) in simulation cases B2, C and D.

Figure 7. (a) Evolution of mean temperature (with symbols) and mean vibrational temperature (without symbols) in simulation cases B2, E and F. (b) Evolution of r.m.s. of temperature (with symbols) and that of the vibrational temperature (without symbols) in simulation cases B2, E and F.

In all these simulations (figures 6(a) and 7(a)), it is evident that the mean translational–rotational energy (represented by $\langle T \rangle$) of the flow is diverted to the mean vibrational energy (represented by $\langle T_v \rangle$). Consequently, the mean vibrational temperature increases with a simultaneous decrease in the mean temperature. In simulation B2, ($D=4$), by the end of 1.5 eddy-turnover times, $\langle T \rangle$ and $\langle T_v \rangle$ become almost equal, indicating the achievement of a vibrational equilibrium state for the mean flow field. However, this does not necessarily mean that the fluctuating flow field has also attained vibrational non-equilibrium. In figures 6(b) and 7(b), we observe that the r.m.s. values of the two temperatures tend to differ substantially from each other even at later times (up to 5 eddy-turnover times).

The new vibrational mechanisms appearing in (2.22) can be expressed in terms of the instantaneous and fluctuating energies. To understand the significance of contributions to these mechanisms from the mean and the fluctuating flow fields separately, we perform the Reynolds decomposition of $e_v$ and $e_v^{*}$

(5.1a,b)\begin{equation} e_v=\langle e_v \rangle + e_v^{'} \quad \text{and}\quad e_v^{*}=\langle e_v^{*}\rangle + e_v^{*'}. \end{equation}

Substituting these in the vibrational mechanism $VI^{'}$–$VII^{'}$ in (2.22) gives

(5.2)\begin{align} & \frac{2}{5}\frac{\partial}{\partial x_j}\left[\frac{1}{\rho}\frac{\partial}{\partial x_i} \left( \frac{\rho}{\tau}e_v\right)\right]-\frac{2}{5}\frac{\partial}{\partial x_j} \left[\frac{1}{\rho}\frac{\partial}{\partial x_i}\left( \frac{\rho}{\tau}e_v^{*}\right)\right] \nonumber\\ &\quad =\frac{2}{5}\frac{\partial}{\partial x_j}\left[\frac{1}{\rho}\frac{\partial}{\partial x_i} \left( \frac{\rho}{\tau}\langle{e_v}\rangle+\frac{\rho}{\tau}e_v^{'}\right)\right] -\frac{2}{5}\frac{\partial}{\partial x_j}\left[\frac{1}{\rho}\frac{\partial}{\partial x_i} \left(\frac{\rho}{\tau}\langle{e_v^{*}}\rangle+\frac{\rho}{\tau}e_v^{*'}\right)\right], \end{align}

which, upon regrouping, leads to

(5.3)\begin{align} & \frac{2}{5}\frac{\partial}{\partial x_j}\left[\frac{1}{\rho}\frac{\partial}{\partial x_i} \left(\frac{\rho}{\tau}e_v\right)\right]-\frac{2}{5}\frac{\partial}{\partial x_j} \left[\frac{1}{\rho}\frac{\partial}{\partial x_i}\left( \frac{\rho}{\tau}e_v^{*}\right)\right] \nonumber\\ &\quad =\frac{2}{5}\frac{\partial}{\partial x_j}\left[\frac{1}{\rho}\frac{\partial}{\partial x_i} \left( \frac{\rho}{\tau}\langle{e_v}\rangle-\frac{\rho}{\tau} \langle{e}_v^{*}\rangle\right)\right]+ \frac{2}{5} \frac{\partial}{\partial x_j}\left[\frac{1}{\rho}\frac{\partial}{\partial x_i} \left(\frac{\rho}{\tau}{e_v^{'}}-\frac{\rho}{\tau}e_v^{*'}\right)\right]. \end{align}

Equation (5.3) shows that the vibrational mechanisms take two kinds of contributions: (i) from the disparity of $\langle {e_v}\rangle -\langle {e_v^{*}}\rangle$, and (ii) from the disparity of ${e_v^{'}}-{e_v^{*'}}$. To gauge the temporal variation in the relative importance of the two contributions, we define, with reference to (5.3), the following fraction ($\chi$):

(5.4)\begin{equation} \chi = \frac{| F_{ij}F_{ij}|}{|F_{ij}F_{ij}| + | M_{ij}M_{ij}|}, \end{equation}

where

(5.5)\begin{equation} M_{ij}=\frac{2}{5}\frac{\partial}{\partial x_j}\left[\frac{1}{\rho} \frac{\partial}{\partial x_i}\left(\frac{\rho}{\tau}\langle{e_v}\rangle- \frac{\rho}{\tau}\langle{e}_v^{*}\rangle\right)\right], \end{equation}

and

(5.6)\begin{equation} F_{ij}=\frac{2}{5}\frac{\partial}{\partial x_j}\left[\frac{1}{\rho} \frac{\partial}{\partial x_i}\left( \frac{\rho}{\tau}{e_v^{'}}-\frac{\rho}{\tau}e_v^{*'}\right)\right]. \end{equation}

Clearly, at a given location in the flow field, the quantity $M_{ij}M_{ij}$ represents a scalar measure of vibrational disparity in the mean flow, whereas $F_{ij}F_{ij}$ represents a scalar measure of vibrational disparity in the fluctuating flow in the context of the vibrational mechanisms $VI$–$VII$. With the initial thermodynamic fluctuations in the flow field being zero in our simulations, the value of $\chi$ is not well defined at that instant. However, as the flow field evolves, $\chi$ starts assuming meaningful values. An increasing value of $\chi$ at a location in the flow field would indicate the reducing significance of the disparity in the mean flow. In figures 8 we present the evolution of the volume-averaged (or the mean) value of $\chi$. In figure 8(a) we present the results from simulations B2, C and D, which differ only in terms of the initial Damköhler number. In figure 8(b) we present the results from simulations B2, E and F, which differ only in terms of the initial value of $\xi$. We observe that, in both sets of simulations, $\langle \chi \rangle$ grows and tends to reach the asymptotic value of unity. In the first set of simulations (B2, C and D), the approach towards unity is most rapid for the case with $D=4$ and is the slowest for the case with $D=0.25$. For the $D=4$ case, $\langle \chi \rangle$ reaches its asymptotic state within one eddy-turnover time, which coincides with the mean vibrational temperature equilibrating with the mean temperature (figure 6a). Similarly, for other simulations, the time instant when $\langle \chi \rangle$ reaches its asymptotic state coincides with the mean vibrational temperature equilibrating with the mean temperature. Furthermore, we observe that, in figure 8(a), $\langle \chi \rangle$ crosses the value of $0.5$ quite early (within two eddy-turnover times) in all the simulations. This is clear evidence that, beyond two eddy-turnover times, in all three simulations (B2, C and D), the disparity in the fluctuating flow field (represented by $F_{ij}F_{ij}$) is actually the major contributor to the vibrational mechanisms ($VI$–$VII$). The same conclusion emerges from figure 8(b), wherein we have included the results of the simulations B2, E and F, which differ only in terms of the initial value of $\xi$. In all these cases, $\langle \chi \rangle$ reaches its asymptotic state at the same time (in all cases, the initial Damköhler number is identical). Evidently, beyond one eddy-turnover time, the disparity in the fluctuating flow field is almost the sole contributor to the vibrational mechanisms ($VI$–$VII$) in the pressure Hessian equation (2.22). These results demonstrate that the disparity in the fluctuating field persists for a long time in decaying turbulence (long after the mean vibrational temperature has equilibrated with the mean temperature in decaying turbulence). This, in turn, means that the vibrational mechanism would remain active throughout the simulation time even when the kinetic energy has decayed to a much lower level compared with its initial state.

Figure 8. Temporal evolution of $\chi$. (a) Influence of Damköhler number. (b) Influence of $\xi$.

Our goal is to examine if the presence of the vibrational non-equilibrium phenomenon affects the velocity gradients and the pressure-Hessian tensors at all in a compressible flow field. Accordingly, in §§ 6 and 7, we examine the influence of the Damköhler number on these tensors. Specifically, we examine (i) the principal strain-rate components (eigenvalues of the strain-rate tensor), (ii) the magnitude of the velocity gradient tensor and (iii) the magnitude of the pressure-Hessian tensor relative to that of the velocity gradient tensor with varying Damköhler number. Flow fields from simulation cases B2, C, D and G are employed for this study.

Subsequently, in § 7, to gain deeper insight into the exact mechanisms by which the vibrational non-equilibrium phenomenon influences the evolution of the pressure Hessian tensor (and in turn the velocity gradient tensor itself), we perform a detailed parametric study of the ‘new’ mechanisms $VI$–$VII$ appearing in the governing equation of the pressure Hessian tensor (2.22). Here, we examine the influence of the initial Damköhler number $\xi$ on these mechanisms. Further, we also discuss the influence of the sign of the local dilatational state on the vibrational mechanism.

6. Influence of vibrational non-equilibrium on the velocity gradient and pressure-Hessian tensors

In compressible turbulence, the scalar aspects of the fluctuating velocity gradient tensor are often examined in terms of the so-called solenoidal dissipation ($\epsilon _s$) and the dilatational dissipation ($\epsilon _d$) of the turbulent kinetic energy

(6.1a,b)\begin{equation} \epsilon_s=\langle\mu\omega_i\omega_i\rangle,\quad \epsilon_d = \left\langle \frac{4\mu}{3} \left(\frac{\partial u_i}{\partial x_i} \right)^2 \right\rangle, \end{equation}

where $\omega _i$ represents the $i$th component of the instantaneous vorticity vector and $({\partial u_i}/{\partial x_i} )$ is the velocity divergence. While, in general, the variation in $\mu$ could also alter the temporal variations of $\epsilon _s$ and $\epsilon _d$, for our simulations (B2, C, D, E, F and G) it is held constant. Thus, the temporal changes in $\epsilon _s$ and $\epsilon _d$ must be directly attributable to the changes in $\omega ^{'2}$ and $\theta ^{'2}$. Accordingly, we directly examine the temporal evolutions of $\omega ^{'}$ and $\theta ^{'}$ in these simulation cases to identify and understand the influence of the vibrational parameters $D$ and $\xi$.

In figure 9(a,b), we present the temporal evolutions of $\theta ^{\prime ^2}$ and $\omega ^{\prime 2}$ from the simulation cases B2, C and D (which differ in terms of the initial Damköhler number alone). In figure 9 (c,d), we present the temporal evolutions of $\theta ^{\prime ^2}$ and $\omega ^{\prime 2}$, from the simulation cases B2, E and F (which differ in terms of the initial ratio of vibrational temperature to local temperature ($\xi$) alone). For reference, we also present the results from simulation case G, which has the same initial Mach number and Reynolds number but no vibrational relaxation process present therein. For each curve, the instantaneous variable has been normalized by the initial value of $\omega ^{\prime 2}$.

Figure 9. Temporal evolutions in simulation cases B2, C, D and G; (a) $\theta ^{'2}$ and (b) $\omega ^{'2}$. Temporal evolutions in simulation cases B2, E, F and G;(c) $\theta ^{'2}$ and (d) $\omega ^{'2}$.

We observe that the influence of the Damköhler number (over the range of this parameter considered in this work) on $\omega ^{\prime 2}$ is negligible. In contrast, the initial Damköhler number does affect the evolution of $\theta ^{\prime 2}$. We observe that, as the initial Damköhler number increases, the magnitude of $\theta ^{\prime ^2}$ decreases, and this difference tends to become larger at later times. In contrast, figure 9(c,d) shows that changes in the initial value of $\xi$ influence neither the statistics of dilatational fluctuations ($\theta ^{\prime ^2}$) nor the vortical fluctuations ($\omega ^{\prime ^2}$).

It is well known that, in an isotropic turbulent flow field, the probability density function (PDF) of $A_{ii}$ shows a single peak distribution centred at zero. However, the distribution becomes wider in a simulation with a higher initial turbulent Mach number (Suman & Girimaji Reference Suman and Girimaji2013). This tendency indicates the existence of a significant number of fluid elements having non-zero $A_{ii}$. While ${A_{ii}}<0$ indicates contracting fluid elements, $A_{ii}>0$ indicates expanding fluid elements. We wish to examine if such a distribution of $A_{ii}$ is affected by the presence of vibrational non-equilibrium processes. In figure 10, we present the PDFs of $A_{ii}$ extracted from simulations B2, C, D and G at a representative time of $t/\tau _{LE}=3.5$. As the vibrational process becomes more rapid (represented by higher $D$), the distribution of $A_{ii}$ is indeed affected. The distributions become narrower with an increase in $D$. This shrinking of the distribution function explains the reduction in $\theta ^{\prime 2}$ evident in figure 9(a).

Figure 10. The PDF of $A_{ii}$. Influence of initial Damköhler number ($D$). Cases B2, C, D and G.

To further examine the influence of the Damköhler number on the velocity gradient field, we turn our attention to the strain-rate tensor. We specifically examine the locally normalized strain-rate tensor ($s_{ij}$), where

(6.2a,b)\begin{equation} s_{ij} = \frac{1}{2} (a_{ij} + a_{ji})\quad \text{and}\quad a_{ij} = \frac{A_{ij}}{\sqrt{A_{mn} A_{mn}}}. \end{equation}

The eigenvalues of the locally normalized strain-rate tensor may not give direct information about the absolute magnitudes. Still, they are an excellent indicator of the relative stretching/compression of a small fluid element in the principal coordinate system of $s_{ij}$. We express these self-normalized eigenvalues of the $\boldsymbol {s}$ tensor ((6.2a,b), Suman & Girimaji Reference Suman and Girimaji2009)

(6.3a–c)\begin{equation} \alpha^* = \frac{\alpha}{\sqrt{\alpha^2 + \beta^2 + \gamma^2}}, \quad \beta^* = \frac{\beta}{\sqrt{\alpha^2 + \beta^2 + \gamma^2}}, \quad \gamma^* = \frac{\gamma}{\sqrt{\alpha^2 + \beta^2 + \gamma^2}}, \end{equation}

where the eigenvalues $(\alpha, \beta, \gamma )$ satisfy the condition

(6.4)\begin{equation} \alpha \geq \beta \geq \gamma. \end{equation}

In incompressible turbulence, it is a well-known behaviour that the ratio $\alpha ^*:\beta ^*:\gamma ^*$ attains a most-probable state of 3:1:–4 (Ashurst et al. Reference Ashurst, Kerstein, Kerr and Gibson1987). This presents the physical picture of a small fluid element being severely stretched and severely compressed in two mutually perpendicular directions. In contrast, the third direction still has a moderate stretching to ensure that the volume of the fluid element is preserved ($\alpha ^*+\beta ^*+\gamma ^*=0$). In compressible turbulence, however, there would be an obvious departure from these proportions because the velocity field is no longer divergence free $(\alpha ^*+\beta ^*+\gamma ^*\neq 0)$. Our intent here is to identify whether the presence of vibrational non-equilibrium in a compressible flow field further influences the proportionality of the three eigenvalues or not. In line with this objective, in figure 11(a–c) we present the PDFs of $\alpha ^{*}$, $\beta ^{*}$ and $\gamma ^{*}$ extracted from simulations which differ in terms of the initial Damköhler number at a representative time $t/\tau _{LE}=3.5$. Figure 11(a) includes PDFs of $\alpha ^{*}$, figure 11(b) includes PDFs of $\beta ^{*}$ and figure 11(c) includes PDFs of $\gamma ^{*}$. All the simulations represented here have identical Reynolds number $(Re_{\lambda }=60)$ and Mach number ($M_t = 0.5$). Considering the no-vibration case (simulation case G) to be our datum for comparison at the specified Mach number and Reynolds number, we observe that the Damköhler number does have an influence on the distributions of normalized eigenvalues. The overall trend is that, at increased Damköhler number, the extreme eigenvalues ($\alpha$ and $\gamma$) tend to decrease in terms of their intensity, but they still retain predominantly positive and negative signs, respectively. On the other hand, the intermediate eigenvalue PDFs seem to be mostly unaffected as the Damköhler number increases.

Figure 11. The PDFs of normalized strain-rate eigenvalues ($\alpha ^* \geq \beta ^* \geq \gamma ^*$). Cases B2, C, D and G.

To further understand how a change in the Damköhler number results in the changes observed in figure 9, we examine the evolution equation of the instantaneous rotation-rate component $W_{ij}$ and the trace of the velocity gradient tensor $A_{ii}$

(6.5)$$\begin{gather} \frac{{\rm D} W_{ij}}{{\rm D} t} ={-}S_{ik}W_{kj} - W_{ik}S_{kj} - 0.5(P_{ij} - P_{ji})+ 0.5(\varGamma_{ij} - \varGamma_{ji}), \end{gather}$$

(6.6)$$\begin{gather}\frac{{\rm D} A_{ii}}{{\rm D} t}={-A_{ik}A_{ki}}-{P_{ii}}+{\varGamma_{ii}}, \end{gather}$$

where

(6.7a,b)\begin{equation} S_{ij} = \tfrac{1}{2} (A_{ij} + A_{ji})\quad \text{and} \quad W_{ij} = \tfrac{1}{2} (A_{ij} - A_{ji}). \end{equation}

In (6.5), $S_{ij}$ is the strain-rate tensor, $W_{ij}$ is the rotation-rate tensor, $P_{ij}$ is the pressure-Hessian tensor and $\varGamma _{ij}$ is the viscous transport term. While some aspects of the dynamics of $\theta ^{\prime 2}$ can be understood with reference to the evolution equation of $A_{ii}$, the evolution equation of the components of the rotation-rate tensor ($W_{ij}$) can prove to be helpful in understanding the behaviour of $\omega ^{\prime 2}$. Using (6.5) and (6.6), we write the evolution equation for the magnitude of the vorticity vector (Wilczek & Meneveau Reference Wilczek and Meneveau2014) and that for $A_{ii}A_{mm}$ (neglecting the viscous processes on the right-hand sides)

(6.8)$$\begin{gather} \frac{{\rm D}}{{\rm D} t} (\omega_i \omega_i) = S_{ik} \omega_k \omega_i - 2 \varepsilon_{ijk} \tilde{P}_{jk} \omega_i, \end{gather}$$

(6.9)$$\begin{gather}\frac{{\rm D}}{{\rm D} t} (A_{ii}A_{mm})={-A_{ik}A_{kj}A_{mm}}-{P_{ii}A_{mm}}, \end{gather}$$

where $\tilde {P}_{ij} = \frac {1}{2} (P_{ij} - P_{ji}),\ \omega _i = -\varepsilon _{ijk} W_{jk}, \text {and}\ \varepsilon _{ijk}$ is the Levi-Civita tensor.

We have neglected the viscous processes in (6.8) and (6.9) as our immediate focus in this study to examine the roles of the inertial and pressure-related processes only. On the right-hand side of these equations, we identify the presence of two types of mechanisms: inertial (the respective first terms in the two equations) and pressure related (the respective second terms in the two equations). To better understand how the initial Damköhler number influences the evolution of $\theta ^{^\prime 2}$ in figure 9(a) but does not influence the evolution of $\omega ^{\prime 2}$ in figure 9(b), we examine the behaviour of these inertial and pressure mechanisms in the simulation cases B2, C, D and G.

In figure 12(a,b), we plot the temporal variations in the mean values of the inertial ($\langle |S_{ik}\omega _k\omega _i|\rangle$) and pressure mechanisms ($\langle |-2\varepsilon _{ijk} \tilde {P}_{jk} \omega _i|\rangle$) in (6.8). The symbol $|\text { }|$ denotes the absolute value of the scalar quantity in context. Every plot in these figures has been normalized with the initial value of $\omega ^{^\prime 3}$. In figure 12(a), the inertial mechanism of (6.8) is not much affected by the changed Damköhler number. On the other hand, in figure 12(b), even though the pressure mechanism of (6.8) seems to be affected by the initial value of the Damköhler number, its strength is evidently too small to have any consequence compared with the inertial mechanism of ((6.8), figure 12a). The strength of the pressure mechanism is small because it takes a contribution from the antisymmetric part of the pressure-Hessian tensor (6.5). Danish et al. (Reference Danish, Suman and Srinivasan2014) demonstrate that, in compressible decaying turbulence, the generation of the antisymmetric portion of the pressure-Hessian tensor is directly influenced by the occurrence of shocklets. With the Mach number of our current simulations being small, the antisymmetric part of the pressure Hessian must be quite small, and thus (even though it might be influenced by the initial Damköhler number) is not large enough to have any consequence for the evolution of the strength of the vorticity vector (figure 9b).

Figure 12. Temporal variation in (a) the inertial mechanism and (b) the pressure mechanism in the evolution equation of (6.8). Simulation cases B2, C, D and G.

In figure 13(a,b) we plot the temporal variations in the mean values of the corresponding inertial $(\langle |-A_{ik}A_{kj}A_{mm}|\rangle )$ and pressure mechanisms $(\langle |-P_{ii}A_{mm}|\rangle )$ of (6.9). We observe that a change in the initial Damköhler number has negligible influence on the inertial mechanism of (6.9), as well (like figure 12a). However, the pressure mechanism does show substantial influence of the initial Damköhler number. Further, unlike the behaviour observed in figure 12(a,b), the order of magnitude of the pressure mechanism in figure 13(b) is indeed comparable to that of the inertial mechanism in figure 13(a). Based on these observations, we conclude that the perceivable influence of the initial Damköhler number on the magnitude of the dilatation fluctuations in figure 9(a) is attributable to the changes in the pressure-Hessian tensor induced by the initial Damköhler number.

Figure 13. Temporal variation in (a) the inertial mechanism and (b) the pressure mechanism in the evolution equation of (6.9). Simulation cases B2, C, D and G.

Suman & Girimaji (Reference Suman and Girimaji2011) have examined the statistics of an aptly normalized form of the algebraic sum of the inertial and pressure mechanisms in the evolution equation of $A_{ii}$ (6.6). This algebraic sum does represent the rate of change in the dilatation rate ($A_{ii}$) following a local small fluid element (while neglecting the role of viscous processes). The authors reported that, as the initial turbulent Mach number increases in their simulation (even with no vibrational physics included therein), the PDF of this quantity spreads, implying the presence of higher rates of changes in the dilatation rates of fluid elements. Following a similar line of investigation as followed by those authors, in this work, we define a quantity, $\delta$

(6.10)\begin{equation} \delta={-}A_{ik}A_{ki}-P_{ii}. \end{equation}

In figure 14(a), we present the temporal evolution of the mean of $\delta ^2$ from the simulation cases B2, C, D and G. All these curves are normalized by the initial value of $\omega ^{\prime 4}$. We observe that, in figure 14(a), as the initial Damköhler number increases, there is indeed a suppression in the magnitude of $\delta$. This suppression does provide an explanation for the suppression in $\theta ^{\prime 2}$ observed earlier in figure 9(a). In figure 14(b) we present the temporal evolution of the mean of $\delta ^2$ in simulation cases B2, E, F (which differ in terms of initial $\xi$ alone) and G. Evidently, the initial value of $\xi$ has no consequence for the evolution of the mean of $\delta ^2$. This observation is consistent with the observation made in figure 9(c) wherein no influence of $\xi$ has been observed on the evolution of $\theta ^{\prime 2}$.

Figure 14. Temporal evolution of $\langle \delta ^2\rangle$ in simulation cases (a) B2, C, D and G and (b) B2, E, F and G.

Suman & Girimaji (Reference Suman and Girimaji2011) have demonstrated that, in the context of the velocity gradient dynamics (2.15), normalizing the pressure-Hessian tensor ($P_{ij}$) with the square of the velocity gradient tensor's magnitude is appropriate. This normalized form delineates the relative importance of the pressure-Hessian process in comparison with the self-stretching terms in (2.15). This relative magnitude is presented as $\phi =\sqrt {P_{mn}P_{mn}}/(A_{ij}A_{ij})$, which is a dimensionless quantity. In figure 15(a), we present the PDF of $\ln \phi$ extracted from simulations B2, C, D and G at the time instant $t/\tau _{LE}=3.5$. All these PDFs are single-peaked distributions. However, we do observe some perceivable effect on these PDFs as the Damköhler number changes. With an increase in the value of the initial Damköhler number, the peak of the distribution shifts to the left, indicating a general decrease in the strength of the pressure-Hessian tensor compared with that of the velocity gradient tensor. In figure 15(b), we present the PDF of $\ln \phi$ extracted from simulations B2, E, F and G at the time instant $t/\tau _{LE}=3.5$. Evidently, the influence of $\xi$ on these PDFs is negligible.

Figure 15. The PDF of $\ln {{\sqrt {P_{mn}P_{mn}}}/{A_{ij}A_{ij}}}$: (a) simulation cases B2, C, D and G and (b) simulation cases B2, E, F and G.

Based on the results of this section, we conclude that, even within the moderate range of the Damköhler number considered in this work, there is a considerable influence of the vibrational non-equilibrium process (represented by varying $D$) on the pressure-Hessian tensor, which in turn, affects the velocity gradient statistics. Overall, an increased Damköhler number suppresses the dilatational fluctuations, whereas the initial value of $\xi$ does not seem to have any significant effect. To gain further insights, in the next section, we directly examine the temporal behaviour of the ‘new’ vibrational mechanisms instrumental in the evolution equation of the pressure-Hessian tensor (3.1).

7. Behaviour of vibrational mechanisms $VI^{'}$ and $VII^{'}$

In this section, we focus on the vibrational mechanisms that have emerged in the evolution equation of the pressure-Hessian tensor (3.1). Accordingly, we examine the behaviour of $f$ (3.2) and attempt to understand how the initial Damköhler number affects the vibrational mechanisms.

To understand the influence of the initial Damköhler number on $f$, we present the temporal evolution of $\langle |\,f|\rangle$ from simulation cases B2, C and D in figure 16(a). We observe that, in each case, $\langle |\,f| \rangle$ undergoes a two-phase evolution. In the first phase of evolution, $\langle |\,f| \rangle$ increases rapidly, reaches a peak value and then decays monotonically. The peak value of $\langle |\,f| \rangle$ is higher in a simulation with a higher initial Damköhler number. Indeed, over the entire duration of the simulations (up to 5 eddy-turnover times), there is an apparent incremental shift in $\langle |\,f| \rangle$ with an increase in the initial Damköhler number at all times. We observe the mean value of $f$ is predominantly negative in all simulations at all times. This implies that $VI < VII$ in (2.22) at most locations of the flow domain. This bias in the distribution is plausibly due to our choice of the initial state of the flow field wherein $T_v< T$ (or equivalently, $e_v< e$), which means $\xi < 1$. Figure 16(a) shows that such bias continues to exist in the flow field at later times too.

Figure 16. Temporal evolution of $|\,f|$. (a) Simulation cases B2, C and D and (b) simulation cases B2, E and F.

We have examined the influence of the initial ratio $\xi = T_v/T$, as well, on the vibrational mechanisms $VI'$ and $VII'$. The temporal evolution of $|\,f|$ from simulations B2, E and F are presented in figure 16(b). Evidently, the influence of initial $\xi$ on $|\,f|$ is significantly smaller than the influence of the initial Damköhler number (figure 16a). This contrasting behaviour does provide an explanation for the observations made earlier in the previous section wherein the statistics of the pressure-Hessian tensor and the velocity gradient tensor seem to be affected by the initial Damköhler number but not much by the initial value of $\xi$.

In incompressible fluid flows, the term $a_{ii}$ is consistently zero for all fluid elements. However, in compressible flows, $a_{ii}$ can assume positive, negative or zero values. A value of $a_{ii} > 0$ signifies an expanding fluid element, $a_{ii} < 0$ indicates contraction and $a_{ii} = 0$ implies an instantaneously volume-preserving fluid element.

Suman & Girimaji (Reference Suman and Girimaji2011) demonstrated that some of the flow statistics (related to velocity gradients) in a compressible flow field may not show any variation by changing a global compressibility parameter like $M_t$. However, when examining the statistics of the same quantities conditioned upon discrete levels of $a_{ii}$ or even upon the sign of $a_{ii}$, significant variations are observed. In our study, it is pertinent to examine whether $f$ exhibits any sensitivity to the local compressibility parameter $a_{ii}$. To examine this, we focus on the simulation case B2 ($M_t=0.5$). In figure 17, we present the PDF of the mean of $f$ conditioned upon the sign of instantaneous $a_{ii}$ at the time instant $t/\tau _{LE}=2.5$. We include two PDFs, separately conditioned upon (a) the contracting and (b) the expanding fluid elements. Clearly, these evolutions do not show much difference compared with the unconditional evolution (included in figure 17 for easy comparison). We have examined the conditioned PDFs of $a_{ii}$ at two more time instants ($t/\tau _{LE}=1.5$, and $t/\tau _{LE}=3.5$) (results not included to avoid repetition). However, the PDFs of $f$ remain identical for all time instants. Based on these results, we conclude that the importance of the vibrational mechanism in (3.1) seems to be similar for the contracting and the expanding fluid elements. It is plausible to attribute these observations to the fact that there is no explicit appearance of $A_{ii}$ and, for that matter, any component of $A_{ij}$, in vibrational mechanisms $VI$ and $VII$ in (2.22).

Figure 17. Conditional PDF of $f$ (3.2). Influence of normalized dilatation rate ($a_{ii}$) (case B2).