2. Kinetic model

The Boltzmann equation was extended to molecular gases by Wang-Chang & Uhlenbeck (Reference Wang-Chang and Uhlenbeck1951), considering quantum mechanics where each internal energy level is assigned an individual velocity distribution function; this yields a complicated operator for particle collisions that is prohibitively expensive for numerical simulations. Several simplified models using the relaxation-time approach have been developed (Morse Reference Morse1964; Holway & Lowell Reference Holway and Lowell1966; Rykov Reference Rykov1975; Gorji & Jenny Reference Gorji and Jenny2013; Wu et al. Reference Wu, White, Scanlon, Reese and Zhang2015; Wang et al. Reference Wang, Yan, Li and Xu2017).

Generally, rarefied gas damping and resonance problems can be investigated by seeking solutions of kinetic equations via the direct simulation Monte Carlo (DSMC) (Hadjiconstantinou Reference Hadjiconstantinou2002; Emerson et al. Reference Emerson, Gu, Stefanov, Yuhong and Barber2007) method or the discrete velocity method (DVM) (Kalempa & Sharipov Reference Kalempa and Sharipov2009; Wu et al. Reference Wu, Reese and Zhang2014). The DSMC method uses a collection of particles to mimic the behaviour of gas molecules: particles move through the spatial space in a realistic manner, while intermolecular collisions and gas interactions are calculated stochastically according to some collision models and pair-selection schemes (Ivanov & Rogasinskii Reference Ivanov and Rogasinskii1991; Bird Reference Bird1994; Roohi et al. Reference Roohi, Stefanov, Shoja-Sani and Ejraei2018). It was rigorously proved that DSMC and the Boltzmann equation are equivalent for monatomic gases in the dilute limit (Wagner Reference Wagner1992). For rarefied molecular gases, the Borgnakke & Larsen (Reference Borgnakke and Larsen1975) phenomenological collision model was developed to reproduce the energy exchange rate, where one continuous variable was introduced to represent the internal energies of a molecular gas. Although DSMC is widely used for rarefied gas flows, it is very time-consuming for oscillating problems (Park et al. Reference Park, Bahukudumbi and Beskok2004). The time step has to be significantly small (compared with both the molecular mean free time and the time period of oscillations) and the total simulation time should be sufficiently long to ensure that the time-periodic state is achieved. In addition, in an unsteady DSMC algorithm, ensemble averaging over thousands of different simulations for each time step is necessary to yield noise-reduced solutions.

The DVM, on the other hand, falls into the category of deterministic approaches. It relies on direct discretisation of the governing equation over computational grids and so can produce noise-free solutions. Furthermore, by assuming small variations in flow properties, the time-periodic flow can be converted into a quasi-steady-state problem and the computational cost will be greatly reduced (see § 3). Therefore, in this work, we use the DVM and a deterministic-based model (Li et al. Reference Li, Zeng, Su and Wu2021; Su et al. Reference Su, Li, Zhang and Wu2022) to investigate wave propagation in rarefied molecular gases, which is modified from the Rykov model (Rykov Reference Rykov1975). The model is able to recover the general temperature and thermal relaxation rates that are predicted by the Wang-Chang–Uhlenbeck equation and can freely adjust the relevant relaxation rates. Therefore, it can simultaneously obtain experimentally measured values of the bulk viscosity and thermal conductivity for a given molecular gas (Wu et al. Reference Wu, Li, Liu and Ubachs2020; Li et al. Reference Li, Zeng, Su and Wu2021; Su et al. Reference Su, Li, Zhang and Wu2022), and their effects on the wave propagation can be separately investigated. To avoid having too many parameters in the analysis, we assume that vibrations of gas molecules are not activated. This assumption has limitations for some polyatomic gases such as carbon dioxide and methane, where vibrational relaxation plays a crucial role in wave attenuation (Ejakov et al. Reference Ejakov, Phillips, Dain, Lueptow and Visser2003; Kustova et al. Reference Kustova, Mekhonoshina, Bechina, Lagutin and Voroshilova2023). In such cases, a model taking into account both rotational and vibrational relaxations (Li et al. Reference Li, Zeng, Huang and Wu2023) is necessary, where additional relaxation rates are required. The rotational mode is considered through the classical mechanics approach, with a continuous variable representing the rotational energy. A brief description of the present kinetic model is given in the following.

In the gas kinetic theory, the state of a molecular gas with excited rotational mode is described by a one-particle velocity–energy distribution function $f (t, \boldsymbol {x}, \boldsymbol {v}, I)$, which is a function of time $t$, spatial coordinate $\boldsymbol {x}=(x,y,z)$, molecular translational velocity $\boldsymbol {v}=(v_x,v_y,v_z)$ and rotational energy $I$. In the absence of external force, the evolution of $f$ is governed by (Rykov Reference Rykov1975; Su et al. Reference Su, Li, Zhang and Wu2022)

(2.1)\begin{equation} \frac{\partial f}{\partial t} + \boldsymbol{v} \boldsymbol{\cdot} \frac{\partial f}{\partial \boldsymbol x} = \underbrace{\frac{p_t\delta_0}{\mu}\left(g_t - f\right)}_{elastic} + \underbrace{\frac{p_t\delta_0}{\mu Z}\left(g_r - g_t\right)}_{inelastic}. \end{equation}

It can be seen that the collision operator (the term on the right-hand side) that describes the change of $f$ due to particle collisions is split into two parts: elastic collisions that preserve translational energy and inelastic ones that exchange translational and rotational energies. Both elastic and inelastic parts are expressed as a simple relaxation term, where the relaxation time related to elastic collisions is $\mu /(p_t\delta _0)$ with $\mu$, $p_t$ and $\delta _0$ being the gas shear viscosity, the pressure related to translational motions and a constant rarefaction parameter, respectively. Here $Z$ is the rotational collision number such that a gas molecule would roughly experience one inelastic collision in every $Z$ collisions, and $g_t$ and $g_r$ are the translational and rotational reference distribution functions, respectively. The reference distributions are series of orthogonal polynomials, expanding at the local equilibrium distribution function $E = E_t(T) \times E_r(T)$, and can be expressed as

(2.2)$$\begin{gather} g_t = E_t(T_t) \times E_r(T_r)\left[1+\frac{4\boldsymbol{q}_t \boldsymbol{\cdot} \boldsymbol{c}}{15T_tp_t} \left(\frac{|\boldsymbol{c}|^2}{T_t}-\frac{5}{2}\right)+\frac{4\boldsymbol{q}_r \boldsymbol{\cdot} \boldsymbol{c}}{dT_tp_r} \left(\frac{I}{T_r}-\frac{d}{2}\right)\right], \end{gather}$$

(2.3)$$\begin{gather}g_r = E_t(T) \times E_r(T)\left[1+\frac{4\boldsymbol{q}' \boldsymbol{\cdot} \boldsymbol{c}}{15Tp}\left(\frac{|\boldsymbol{c}|^2}{T}-\frac{5}{2}\right)+\frac{4\boldsymbol{q}'' \boldsymbol{\cdot} \boldsymbol{c}}{dTp}\left(\frac{I}{T}-\frac{d}{2}\right)\right], \end{gather}$$

with $E_t$ and $E_r$ being the local Maxwellian distributions

(2.4)$$\begin{gather} E_t(T) = \frac{n}{({\rm \pi} T)^{3/2}}\exp\left(-\frac{|\boldsymbol{c}|^2}{T}\right), \end{gather}$$

(2.5)$$\begin{gather}E_r(T) = \frac{I^{d/2-1}}{\varGamma(d/2) T^{d/2}}\exp\left(-\frac{I}{T}\right). \end{gather}$$

In the above formulae: $d$ is the number of rotational degrees of freedom; $\varGamma$(${\cdot }$) is the gamma function; $T_t$, $T_r$, $\boldsymbol {q}_t=(q_{t}^{x},q_{t}^{y},q_{t}^{z})$ and $\boldsymbol {q}_r=(q_{r}^{x},q_{r}^{y},q_{r}^{z})$ are the translational and rotational temperatures and the related heat fluxes, respectively; and $\boldsymbol {c} = \boldsymbol {v} - \boldsymbol {u}$ is the peculiar velocity with $\boldsymbol {u}=(u_x,u_y,u_z)$ being the gas bulk velocity. The overall temperature $T$ is a weighted sum of the translational and rotational temperatures: $T=(3T_t+dT_r)/(3+d)$. The pressures are defined as $p_t = nT_t$ and $p = nT$ in terms of the translational and overall temperatures, respectively, where $n$ is the gas number density. Here $\boldsymbol {q}'$ and $\boldsymbol {q}''$ are two auxiliary heat fluxes, defined as linear combinations of the translational and rotational heat fluxes (Li et al. Reference Li, Zeng, Su and Wu2021):

(2.6)\begin{equation} \begin{bmatrix} \boldsymbol{q}'\\ \boldsymbol{q}'' \end{bmatrix}= \begin{bmatrix} (2-3A_{tt})Z+1 & -3A_{tr}Z \\ -A_{rt}Z & 1-A_{rr}Z \end{bmatrix} \begin{bmatrix} \boldsymbol{q}_t \\ \boldsymbol{q}_r \end{bmatrix}, \end{equation}

where $A_{ij}$ ($i,j = t$ or $r$) are the thermal relaxation rates.

Here, the density and temperatures are normalised by the reference density $n_0$ and temperature $T_0$, respectively; the shear viscosity by its value at the reference temperature $\mu _0$; velocities by the most probable speed $v_m = \sqrt {2k_BT_0/m}$, where $k_{B}$ is the Boltzmann constant and $m$ is the molecular mass; spatial coordinates by the characteristic flow length $H$; time by $H/v_m$; the internal energy by $k_BT_0$; heat fluxes by $n_0k_BT_0v_m$; pressures by $n_0k_BT_0$; and the distribution functions by $n_0v_m^{-3}(k_BT_0)^{-1}$. The rarefaction parameter is therefore defined as

(2.7)\begin{equation} \delta_0 = \frac{p_0H}{\mu_0 v_{m}}, \end{equation}

where $p_0=n_0k_BT_0$. It is inversely proportional to the unconfined length-based Knudsen number:

(2.8)\begin{equation} Kn_l=\frac{\sqrt{\rm \pi}}{2}\frac{1}{\delta_0}. \end{equation}

It is straightforward to show that the kinetic model recovers the Jeans–Landau temperature relaxation model, i.e.

(2.9)\begin{equation} \frac{\partial T_r}{\partial t}=\frac{p_t\delta_0}{\mu Z}\left(T-T_r\right). \end{equation}

In addition, the heat flux relaxation of the present model is consistent with that derived from the Wang-Chang–Uhlenbeck equation by Mason & Monchick (Reference Mason and Monchick1962). The heat flux relaxation directly determines the thermal conductivity, as well as the translational and rotational contributions (Mason & Monchick Reference Mason and Monchick1962; McCormack Reference McCormack1968), i.e.

(2.10) \begin{equation} \left[\begin{array}{@{}c@{}} \partial\boldsymbol{q}_t/\partial t\\ \partial\boldsymbol{q}_r/\partial t \end{array}\right]={-}\frac{p_t\delta_0}{\mu}\left[\begin{array}{@{}cc@{}} A_{tt} & A_{tr}\\ A_{rt} & A_{rr} \end{array}\right]\left[\begin{array}{@{}c@{}} \boldsymbol{q}_t\\ \boldsymbol{q}_r \end{array}\right]. \end{equation}

Applying the Chapman–Enskog multiscale expansion to the kinetic model equation and retaining the terms up to the order of $O(1/\delta )$, the Navier–Stokes equations can be derived, and the transport coefficients are obtained immediately. The dimensionless bulk viscosity $\mu _b$ (normalised by $\mu _0$) is given by

(2.11)\begin{equation} \mu_b = \frac{2dZ}{3(d+3)}\mu, \end{equation}

whereas the dimensionless translational conductivity $\kappa _t$ and rotational conductivity $\kappa _r$ (normalised by $2k_B\mu _0/m$) are determined as

(2.12)\begin{equation} \begin{bmatrix} \kappa _{t}\\ \kappa _{r} \end{bmatrix}= \frac{\mu}{4} \begin{bmatrix} A_{tt} & A_{tr} \\ A_{rt} & A_{rr} \end{bmatrix}^{{-}1} \begin{bmatrix} 5 \\ d \end{bmatrix}. \end{equation}

In practice, it is more convenient to express the thermal conductivity in terms of the Eucken factors (Eucken Reference Eucken1913):

(2.13a–c)\begin{equation} \frac{2\kappa_t}{\mu}=\frac{3}{2}f_{t},\quad \frac{2\kappa_r}{\mu}=\frac{d}{2}f_{r},\quad \frac{2\kappa}{\mu}=\frac{2\left(\kappa_t+\kappa_r\right)}{\mu} = \frac{3+d}{2}f_{eu}, \end{equation}

where $\kappa =\kappa _t+\kappa _r$ is the overall thermal conductivity and $f_{eu}$, $f_t$ and $f_r$ are the total, translational and rotational Eucken factors, respectively.

For the numerical solution of our kinetic equation, a boundary condition is required to determine the value of $f$ at the boundary of a computational domain. Considering a non-absorbing wall with velocity $\boldsymbol {u}_w$ at temperature $T_w$, all the gas molecules $(\boldsymbol {v}',\boldsymbol {I}')$ hitting the wall will return to the flow field with a new state $(\boldsymbol {v},I)$. Given the unit normal vector, $\boldsymbol {n}$, of the wall pointing towards the gas, the velocity–energy distribution function of the molecules in the vicinity of the wall is

(2.14)\begin{equation} f_w =\left\{\begin{array}{@{}cc} f^{-}, & (\boldsymbol{v}'-\boldsymbol{u}_w) \boldsymbol{\cdot} \boldsymbol{n} \le 0, \\ f^+, & (\boldsymbol{v}'-\boldsymbol{u}_w) \boldsymbol{\cdot} \boldsymbol{n} > 0, \end{array}\right. \end{equation}

where $f^-$ and $f^+$ are the distributions of the incident and reflected molecules, respectively. The correlation between the incident and reflected distribution functions is defined through the reflection kernel $R(\boldsymbol {v}' \rightarrow \boldsymbol {v}, I' \rightarrow I)$ as

(2.15)\begin{align} & (\boldsymbol{v}-\boldsymbol{u}_w) \boldsymbol{\cdot} \boldsymbol{n} f^+(\boldsymbol{v} ,I) \nonumber\\ &\quad ={-}\iint _{({\boldsymbol{v}'-\boldsymbol{u}_w)}\ \boldsymbol{\cdot}\ \boldsymbol{n} \le 0} (\boldsymbol{v}'-\boldsymbol{u}_w) \boldsymbol{\cdot} \boldsymbol{n}f^{-}(\boldsymbol{v}',I')R(\boldsymbol{v}' \rightarrow \boldsymbol{v}, I' \rightarrow I)\,{\rm d}\boldsymbol{v}'\,{\rm d}I'\end{align}

for all $(\boldsymbol {v}'-\boldsymbol {u}_w) \boldsymbol {\cdot } \boldsymbol {n} > 0$. For a fully diffuse wall, the reflection kernel is

(2.16)\begin{equation} R(\boldsymbol{v}' \rightarrow \boldsymbol{v}, I' \rightarrow I) = \frac{2(\boldsymbol{v}-\boldsymbol{u}_{w}) \boldsymbol{\cdot} \boldsymbol{n}I^{d/2-1}}{{\rm \pi} \varGamma (d/2)T_{w}^{2+d/2}}\times \exp \left(-\frac{(\boldsymbol{v}-\boldsymbol{u}_w)^2}{T_w}-\frac{I}{T_w}\right). \end{equation}

3. Formulation of wave propagation

The schematic diagram of the simulation set-up is shown in figure 1. Here, we specify the formulation of wave propagation in rarefied molecular gases confined between two infinite, parallel and impermeable plates located at $x = 0$ and $x = 1$, which can be treated as a simple prototype of wave propagation in confined space. The plate at $x=0$ may oscillate harmonically at a constant frequency $\omega$ in its normal direction so that the velocity varies as

(3.1)\begin{equation} u_w(t) = U_{m}{\rm Re}\{\exp\left({\rm i}{St}\,t\right)\}, \end{equation}

and a sound wave is induced in the gas. This plate may attain unsteady heating so that its temperature varies as

(3.2)\begin{equation} T_w(t) = 1+\tau_{m}{\rm Re}\{\exp\left({\rm i}{St}\,t\right)\}, \end{equation}

and a thermal–acoustic wave is consequently induced. In (3.1) and (3.2), ‘${\rm i}$’ is the imaginary unit, ${\rm Re}\{{\cdot }\}$ denotes the real part of a complex expression, and $U_m$ and $\tau _m$ are the maximum speed and perturbed temperature of the plate, respectively. The constant ${St}$ is the normalised oscillation frequency, defined as the Strouhal number

(3.3)\begin{equation} {St} = \frac{\omega H}{v_{m}}. \end{equation}

Note that the two driving forces have the same frequency. The time-based Knudsen number correlates with ${Kn}_l$ and ${St}$ through

(3.4)\begin{equation} {Kn}_t=\frac{\mu_0\omega}{p_0}=\frac{2{Kn}_l{St}}{\sqrt{\rm \pi}}, \end{equation}

where $p_0/\mu _0$ is the reference collision frequency of gas molecules.

Figure 1.

Schematic diagram of sound and thermal–acoustic wave propagation in rarefied molecular gases between two infinite, parallel and impermeable plates. The oscillating plate is positioned at $x=0$ with harmonically varying velocity (along the $x$ direction) or temperature, and the stationary plate is at $x=1$. The stationary plate has a fixed temperature of $T_{w}=1$. Both plates are fully diffuse walls.

We are interested in the flow state when the oscillation has been fully established. As a consequence, the time-dependent behaviour of the gas is periodic with the same frequency as the oscillation stimuli. We further assume that the velocity and perturbed temperature of the plate are sufficiently small quantities, i.e. $u_m\ll 1$ and $\tau _m\ll 1$; hence, the gas system deviates slightly from its reference equilibrium state, which is described by the global Maxwellian

(3.5)\begin{equation} E_0 = \frac{I^{d/2-1}}{{\rm \pi} ^{3/2} \varGamma (d/2)} \exp\left(-|\boldsymbol{v}|^2-I\right), \end{equation}

and the kinetic equations can be linearised by representing the distribution function as

(3.6)\begin{equation} f(t,x,\boldsymbol{v},I) = E_0[1+{\rm Re}\{h(x,\boldsymbol{v},I)\exp({\rm i}{St}\,t)\}], \end{equation}

where $|h| \ll 1$ is the amplitude of the perturbed distribution function. The macroscopic flow properties of interest are also expressed by the following complex functions:

(3.7) \begin{equation} \left.\begin{gathered} n(t, x) = 1 + {\rm Re}\left\{\rho(x)\exp\left({\rm i}{St}\,t\right)\right\},\\ u_x(t,x) = {\rm Re}\left\{U(x)\exp\left({\rm i}{St}\,t\right)\right\},\\ T_t(t,x) = 1 + {\rm Re}\left\{\tau_t(x)\exp\left({\rm i}{St}\,t\right)\right\}, \\ T_r(t,x) = 1 +{\rm Re}\left\{\tau_r(x)\exp\left({\rm i}{St}\,t\right)\right\},\\ q^x_t(t,x) ={\rm Re}\left\{Q_t(x)\exp\left({\rm i}{St}\,t\right)\right\},\\ q^x_r(t,x) ={\rm Re}\left\{Q_r(x)\exp\left({\rm i}{St}\,t\right)\right\},\\ p_{xx}(t,x) =1+{\rm Re}\left\{P(x)\exp\left({\rm i}{St}\,t\right)\right\}, \end{gathered}\right\} \end{equation}

where $p_{xx}$ is the gas pressure in the direction of wave propagation. Here $[\rho,U,\tau _t,\tau _t,Q_t,Q_r,P]$ are complex quantities, which can be represented in exponential form as

(3.8a,b)\begin{equation} M\left(x\right)=M_{am}\left(x\right)\exp\left[{\rm i}\phi_M\left(x\right)\right],\quad M=\rho,U,\tau_t,\tau_r,Q_t,Q_r,P, \end{equation}

where $M_{am}$ and $\phi$ are the real functions corresponding to the amplitude and phase of a periodic function. Note that for the spatial coordinates, all the variables depend only on $x$, since the propagation of the wave is restricted to this direction.

Substituting (3.6) and (3.7) into the kinetic equation (2.1) and neglecting all the nonlinear terms, we obtain the equation for $h$, written as

(3.9) \begin{align} {\rm i}St\,h + v_x \frac{\partial h}{\partial x} &= \delta_0 \left[\rho + 2Uv_x + \tau_t \left( |\boldsymbol{v}|^{2}-\frac{3}{2}\right)+ \tau_r \left(I-\frac{d}{2}\right)\right. \nonumber\\ &\quad +\left.\frac{4}{15}Q_tv_x \left(|\boldsymbol{v}|^{2}-\frac{5}{2}\right)+ \frac{4}{d}Q_rv_x \left(I-\frac{d}{2}\right)-h\right] \nonumber\\ &\quad + \frac{\delta}{Z} \left[(\tau - \tau_t)\left(|\boldsymbol{v}|^2-\frac{3}{2}\right)+(\tau - \tau_r) \left(I-\frac{d}{2}\right) \right. \nonumber\\ &\left.\quad +\frac{4(Q'-Q_t) \boldsymbol{v}}{15}\left(|\boldsymbol{v}|^2-\frac{5}{2}\right)+\frac{4(Q''-Q_r) \boldsymbol{v}}{d} \left(I-\frac{d}{2}\right)\right], \end{align}

where $\tau =(3\tau _t+d\tau _r)/(3+d)$, and $Q'$ and $Q''$ are related to $Q_t$ and $Q_r$ according to (2.6). The boundary condition for $h$ is

(3.10) \begin{equation} \left.\begin{gathered} h^+ = \sqrt{\rm \pi}U_m - 2\sqrt{\rm \pi} \iint _ {v'_x \le 0} v'_xE_{0}h^{-}\,{\rm d}\boldsymbol{v}'\,{\rm d}I' \\ \quad\quad\quad\quad\quad\quad\qquad\quad +\, 2U_mv_x + \tau_m \left(|\boldsymbol{v}|^{2}-\frac{5}{2}\right)+ \tau_m\left(I-\frac{d}{2}\right)\quad \text{at}\ x=0,\\ h^+ = 2\sqrt{\rm \pi} \iint _ {v'_x > 0} v'_xE_{0}h^{-}\,{\rm d}\boldsymbol{v}'\,{\rm d}I'\quad \text{at}\ x=1, \end{gathered}\right\} \end{equation}

where $h^+$ and $h^-$ denote the reflected and incident perturbed distribution functions. Once $h$ is known, the macro-quantities are calculated from the velocity moments as

(3.11)\begin{align} & \left[\rho, U, \tau_t, \tau_r, Q_t, Q_r, P\right] \nonumber\\ &\quad =\iint \left[1, v_x, \frac{2}{3}|\boldsymbol{v}|^{2}-1, \frac{2}{d}I-1, v_x\left(|\boldsymbol{v}|^{2}-\frac{5}{2}\right),v_x \left(I-\frac{d}{2}\right),v^2_x\right]E_{0}h\,{\rm d}\boldsymbol{v}\,{\rm d}I. \end{align}

By introducing complex expressions, the time variable in the governing equation system is eliminated. The technique of computational fluid dynamics for time-independent problems is adopted to find deterministic solutions, which can greatly reduce the computational cost. The details about the computational issues are given in the following section and Appendix A. Note that this method cannot be used to simulate the propagation of a strong wave with large amplitude. In such a circumstance, the assumption that variant flow properties have the same oscillating frequency as the external stimuli may break down, and a weak shock wave may be generated and propagate in space (Cox, Mortell & Reck Reference Cox, Mortell and Reck2001; Tang, Cheng & Xu Reference Tang, Cheng and Xu2001). Then, time-dependent nonlinear kinetic equations should be applied to investigate such flows. Numerical methods with efficient shock-capture schemes, e.g. the unified gas kinetic schemes, should be used (Wang & Xu Reference Wang and Xu2012; Wang et al. Reference Wang, Zhu, Su, Wu and Zhang2018).

5. Onsager–Casimir reciprocal relationship

In this section, we examine whether the OCRR also holds for molecular gases. If we consider a weakly non-equilibrium system where some irreversible processes occur, then the corresponding physical law can be described in a general linear form as

(5.1)\begin{equation} J_{k} = \sum _{n} \varLambda _{kn}X_{n}, \end{equation}

where $X_{n}$ are driven forces, $J_{k}$ are conjugated thermodynamic fluxes and $\varLambda _{kn}$ is the matrix of kinetic coefficients. The non-diagonal elements $\varLambda _{kn}$ ($k\neq n$) contain solutions corresponding to the forces $X_{k}$ and $X_{n}$. If the form of the thermodynamic flow $J_{k}$ is chosen so that the entropy production $s$ in the statistical system is expressed as

(5.2)\begin{equation} s = \sum _{k} J_{k}X_{k}, \end{equation}

then the kinetic coefficients satisfy

(5.3)\begin{equation} \varLambda_{kn}^{t} = \varLambda_{nk}^{t}, \end{equation}

where the superscript $t$ means that each coefficient $\varLambda _{kn}$, composed of two solutions $k$ and $n$, is calculated from the time-reverted molecular state of one of these solutions.

From the Boltzmann equation, the thermodynamic fluxes and forces can be identified for gas flows, and the kinetic coefficients are derived as (Sharipov Reference Sharipov2006)

(5.4) $$\begin{gather} \varLambda_{kn} = \int_{\varSigma}\left(\boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{n} h_{w}^{(k)},h^{(n)}- \frac{1}{2}h_{w}^{(n)}\right){\rm d} \varSigma - {\rm i}St\left(\left(h^{(k)}, h^{(n)}\right)\right), \end{gather}$$

(5.5) $$\begin{gather} \varLambda_{kn}^{t} = \int_{\varSigma}\left(\mathcal{T} \boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{n} h_{w}^{(k)},h^{(n)}-\frac{1}{2}h_{w}^{(n)}\right){\rm d} \varSigma - {\rm i} St\left(\left(\mathcal{T}h^{(k)}, h^{(n)}\right)\right), \end{gather}$$

where $\varSigma$ represents solid walls that enclose the flow field $\varOmega$, $h_w$ is given by the wall velocity and temperature as

(5.6)\begin{equation} h_w=2\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{u}_w+\tau_w\left(|\boldsymbol{v}|^2-\frac{5}{2}\right)+\tau_w\left(I-\frac{d}{2}\right), \end{equation}

and the binary operations $({\cdot },{\cdot })$ and $(({\cdot },{\cdot }))$ are defined for the two arbitrary functions as

(5.7)$$\begin{gather} \left(\phi, \psi \right) = \iint E_{0}\phi(\boldsymbol{x}, \boldsymbol{v}, I) \psi(\boldsymbol{x}, \boldsymbol{v}, I)\,{\rm d}\boldsymbol{v}\,{\rm d}I, \end{gather}$$

(5.8)$$\begin{gather}((\phi, \psi )) = \int_{\varOmega}\left(\phi, \psi \right){\rm d}\boldsymbol{x}. \end{gather}$$

Finally, $\mathcal {T}$ is the operator reversing the state of molecules in time,

(5.9)\begin{equation} \mathcal{T}\phi(\boldsymbol{x}, \boldsymbol{v}, I) = \phi(\boldsymbol{x}, \boldsymbol{v}^{t}, I^{t}), \end{equation}

which is self-conjugate, i.e.

(5.10)\begin{equation} \left(\mathcal{T}\phi, \psi\right) = \left(\mathcal{T}\psi, \phi\right). \end{equation}

For the sound wave stimulated by oscillation of the transducer and the thermoacoustic wave induced by a temperature oscillation, the thermodynamic forces are defined as

(5.11)\begin{equation} X_{U} = u_{m},\quad X_{T} = \tau_{m}, \end{equation}

and then the perturbed distribution can be expressed as

(5.12)\begin{equation} h(x, \boldsymbol{v}, I) = h^{(U)}(x, \boldsymbol{v}, I)X_{U}+h^{(T)}(x, \boldsymbol{v}, I)X_{T}. \end{equation}

In addition, $h_{w}$ that defines the boundary condition can also be decomposed as

(5.13)\begin{equation} h_{w}(x, \boldsymbol{v}, I) = h_{w}^{(U)}(x, \boldsymbol{v}, I)X_{U}+h_{w}^{(T)}(x, \boldsymbol{v}, I)X_{T}, \end{equation}

where

(5.14)$$\begin{gather} h_{w}^{(U)} =\left\{\begin{array}{@{}lc} 2v_{x}, & x = 0, \\ 0, & x = 1, \end{array}\right. \end{gather}$$

(5.15)$$\begin{gather}h_{w}^{(T)} = \left\{\begin{array}{@{}lc} \left(|\boldsymbol{v}|^{2}-\dfrac{5}{2}\right)+\left(I-\dfrac{d}{2}\right), & x = 0, \\ 0, & x = 1. \end{array}\right. \end{gather}$$

Using the reciprocal relation, i.e. (5.1), the kinetic coefficients given by (5.4) and (5.5), and the self-conjugate property described by (5.10), we obtain

(5.16)\begin{align} & \left(\mathcal{T}v_{x}h_{w}^{(U)},h^{(T)}\right)_{x = 0} - \left(\mathcal{T}v_{x}h_{w}^{(U)},h^{(T)}\right)_{x = 1} \nonumber\\ &\quad = \left(\mathcal{T}v_{x}h_{w}^{(T)},h^{(U)}\right)_{x = 0} - \left(\mathcal{T}v_{x}h_{w}^{(T)},h^{(U)}\right)_{x = 1} \nonumber\\ &\qquad +\left(\mathcal{T}v_{x}h_{w}^{(U)},h_{w}^{(T)}\right)_{x = 0} - \left(\mathcal{T}v_{x}h_{w}^{(U)},h_{w}^{(T)}\right)_{x = 1}. \end{align}

Inserting $h_{w}^{(U)}$ and $h_{w}^{\textrm {T}}$, we have

(5.17)\begin{align} \left(\mathcal{T}2|\boldsymbol{v}|^{2}, h^{(T)}\right)_{x = 0} &= \left(\mathcal{T}v_{x}\left(|\boldsymbol{v}|^{2}-\frac{5}{2}\right),h^{(U)}\right)_{x = 0} + \left(\mathcal{T}v_{x}\left(I-\frac{d}{2}\right),h^{(U)}\right)_{x = 0} \nonumber\\ &\quad + \left(\mathcal{T}2|\boldsymbol{v}|^{2}, \left(|\boldsymbol{v}|^{2}-\frac{5}{2}\right)\right)_{x = 0} + \left(\mathcal{T}2|\boldsymbol{v}|^{2}, \left(I-\frac{d}{2}\right)\right)_{x = 0}, \end{align}

and thus

(5.18)\begin{equation} P^{(T)} ={-}Q_{t}^{(U)}-Q_{r}^{(U)} ={-}Q^{(U)},\quad x = 0. \end{equation}

If we use $h^{(T)}(1-x,\boldsymbol {v},I)$ instead of $h^{(T)}(x,\boldsymbol {v},I)$, the boundary condition is then transformed into

(5.19)\begin{equation} h_{w}^{(T)}(1-x) =\left\{\begin{array}{@{}lc} 0, & x = 0,\\ \left(|\boldsymbol{v}|^{2}-\dfrac{5}{2}\right)+\left(I-\dfrac{d}{2}\right), & x = 1. \end{array}\right. \end{equation}

The reciprocal relation gives

(5.20)\begin{align} \left(\mathcal{T}2|\boldsymbol{v}|^{2}, h^{(T)}(1-x)\right)_{x = 0} &={-} \left(\mathcal{T}v_{X}\left(|\boldsymbol{v}|^{2}-\frac{5}{2}\right),h^{(U)}(x)\right)_{x = 1} \nonumber\\ &\quad - \left(\mathcal{T}v_{X}\left(I-\frac{d}{2}\right),h^{(U)}(x)\right)_{x = 1}, \end{align}

and thus

(5.21)\begin{equation} P^{(T)} = Q_{t}^{(U)}+Q_{r}^{(U)} = Q^{(U)},\quad x = 1. \end{equation}

Therefore, the OCRR for wave propagation in molecular gases is confirmed by (5.18) and (5.21), which show that the magnitude of the deviation of gas pressure at the walls induced by the thermoacoustic wave is equal to the magnitude of the total heat flux (sum of the translational and internal heat fluxes) induced by the sound wave.

We further show the numerical results for the pressure amplitude and the total heat flux under a wide range of ${St}$ and ${Kn}_{l}$. The results on the transducer and receptor are shown in figure 9. In all the flow regimes, the pressure amplitude in the thermoacoustic wave agrees with the heat flux in the sound wave, suggesting that the OCRR also holds for molecular gases.

Figure 9.

The variation of pressure amplitude and heat flux under different $St$: (a) $x = 0$ and (b) $x = 1$. The solid lines denote the pressure amplitude while the hollow points represent the heat flux. Three colours denote different flow regimes: red represents the slip regime ($Kn_{l} = 0.01$), blue represents the transition regime ($Kn_{l} = 0.1$) and yellow represents the free molecular flow regime ($Kn_{l} = 10$).

The comparison between the pressure amplitude of the thermoacoustic wave and the total heat flux of the sound wave under different values of the bulk viscosity and thermal conductivity is shown in figure 10. The simulation is conducted in the slip regime with $Kn = 0.01$. It is evident that the OCRR is valid under different values of bulk viscosity and thermal conductivity. From the left column of figure 10, we can see the influence of the bulk viscosity and thermal conductivity on the thermoacoustic wave. In figure 10(a), the resonance frequency is shifted when the bulk viscosity is altered, which is similar to the behaviour of the sound wave. The only difference lies in the value of $P_{am}$. In the case of the thermoacoustic wave, the pressure variation is generated by the oscillation of temperature, resulting in a significantly lower pressure amplitude than that in the sound wave. Figure 10(c,d) shows that thermal conductivity has little influence on the pressure amplitude of the thermoacoustic wave and the heat flux amplitude of the sound wave.

Figure 10.

Comparison between the pressure amplitude of the thermoacoustic wave and the total heat flux of the sound wave under different values of the bulk viscosity and thermal conductivity. The pressure amplitude of the thermoacoustic wave is shown in (a,c) while the total heat flux of the sound wave is shown in (b,d).