#  $H$-theorem and boundary conditions for the linear R26 equations:
application to flow past an evaporating droplet

Anirudh S. Rana1 , Vinay Kumar Gupta2<EMAIL_ADDRESS>James E. Sprittles3  Manuel Torrilhon4  1 Department of Mathematics, BITS
Pilani, Pilani Campus, 333031 Rajasthan, India 2 Department of Mathematics,
Indian Institute of Technology Indore, Indore 453552, India 3 Mathematics
Institute, University of Warwick, Coventry CV4 7AL, UK 4 Department of
Mathematics, RWTH Aachen University, D-52062 Aachen, Germany

###### Abstract

Determining physically admissible boundary conditions for higher moments in an
extended continuum model is recognised as a major obstacle. Boundary
conditions for the regularised 26-moment (R26) equations obtained using
Maxwell’s accommodation model do exist in the literature; however, we show in
this article that these boundary conditions violate the second law of
thermodynamics and the Onsager reciprocity relations for certain boundary
value problems, and hence are not physically admissible. We further prove that
the linearised R26 (LR26) equations possess a proper $H$-theorem (second-law
inequality) by determining a quadratic form without cross-product terms for
the entropy density. The establishment of the $H$-theorem for the LR26
equations in turn leads to a complete set of boundary conditions that are
physically admissible for all processes and comply with the Onsager
reciprocity relations. As an application, the problem of a slow rarefied gas
flow past a spherical droplet with and without evaporation is considered and
solved analytically. The results are compared with the numerical solution of
the linearised Boltzmann equation, experimental results from the literature
and/or other macroscopic theories to show that the LR26 theory with the
physically admissible boundary conditions provides an excellent prediction up
to Knudsen number $\lesssim 1$ and, consequently, provides transpicuous
insights into intriguing effects, such as thermal polarisation. In particular,
the analytic results for the drag force obtained in the present work are in an
excellent agreement with experimental results even for very large values of
the Knudsen number.

##  1 Introduction

The behaviour of aerosols (small droplets suspended in vapour/gases) covers a
wide range of phenomena of interest to both scientists and engineers (Sazhin,
2014). In the context of airborne disease transmission (Tang et al., 2006; Xie
et al., 2007), large droplets ($\gtrsim\mu$m in radius) fall rapidly to the
ground due to gravity, and are therefore transmitted over short distances;
smaller droplets ($\lesssim\mu$m in radius), on the other hand, remain
suspended in the air for a significant period of time, and hence are
transmitted over large distances. The immediate intuitive questions are “How
far does a droplet travel?” and “How long does a droplet survive?” The former
is answered by evaluating the momentum of the droplet while the latter by the
heat and mass transfer properties of the droplet, with the two being
inevitably coupled. At small scales, it is difficult to measure the relevant
quantities (e.g. drag forces) on single and multiple droplets; therefore, an
accurate modelling is key to understanding transport phenomena (the trajectory
and thermal state) during evaporation.

When a spherical particle of radius $\hat{R}_{0}$ moves slowly with a steady
translational velocity in a gas, with $\hat{R}_{0}$ being much larger than the
mean free path ($\hat{\lambda}$) of the gas, the flow fields around the
particle can be obtained by solving the Stokes equation, and therefore one can
obtain the drag force acting on the sphere (Lamb, 1932). However, the Stokes
equation as well as the other classical continuum theories of gas dynamics
based on the Euler or Navier-Stokes-Fourier (NSF) equations, which are in
general accurate for describing a flow in the so-called hydrodynamic regime
(Knudsen number $\mathrm{Kn}\lesssim 0.001$), cannot explain the gas flow past
a micro-/nano-sized drop (for which $\mathrm{Kn}=\hat{\lambda}/\hat{R}_{0}\sim
1$) accurately (Cercignani, 2000; Sone, 2002; Struchtrup, 2005). Such flows
can only be predicted by employing tools of kinetic theory.

Within the framework of kinetic theory, the vapour/gas is described by the
Boltzmann equation (Cercignani, 1975; Kremer, 2010), which is the evolution
equation for the (velocity) distribution function of the gas and accurately
describes flows for all Knudsen numbers. Notwithstanding, the exact analytical
solution of the Boltzmann equation is intractable for a general flow problem
while numerical solutions of the Boltzmann equation demand a very high
computational cost in general (particularly, for transition-regime flows).
Although, very recently, some general synthetic iterative schemes have been
proposed for solving the Boltzmann equation (Su et al., 2021; Zhu et al.,
2021) that reduce the computational cost of numerical solutions of the
Boltzmann equation significantly, it is still appreciable to describe a
process with a set of partial differential equations in physical variables (or
moments of the distribution function) as they provide a deeper insight to the
underlying physics of the process and permit the deployment of classical
mathematical techniques. A variety of approximation techniques, aiming to
derive sets of partial differential equations that can describe processes at
moderate Knudsen numbers, have been developed; see e.g. Struchtrup (2005);
Torrilhon (2016) and references therein. Among others, the two well-known
approximation techniques are the Chapman–Enskog (CE) expansion (Chapman &
Cowling, 1970) and method of moments (Grad, 1949).

The CE expansion is suitable for near-equilibrium gas flows and relies on an
asymptotic expansion of the distribution function around the equilibrium
distribution function (which is a Maxwellian) in powers of a small
parameter—typically the Knudsen number. The distribution function in the
Boltzmann equation is replaced with this expanded distribution function, and
coefficients of each power of the Knudsen number on both sides of the
resulting equation are compared. The comparison of coefficients of
$\mathrm{Kn}^{0}$ yields the Euler equations, of $\mathrm{Kn}^{1}$ results
into the NSF equations, of $\mathrm{Kn}^{2}$ leads to the Burnett equations,
of $\mathrm{Kn}^{3}$ gives the so-called super-Burnett equations, and so on
(Struchtrup, 2005; Bobylev, 2018). The Burnett (also the super-Burnett and
beyond) equations are unstable (Bobylev, 1982), and due to lack of any
coherent entropy inequality these equations lead to physically inadmissible
solutions (Struchtrup & Nadler, 2020).

Grad—in a seminal work (Grad, 1949)—proposed an asymptotic solution procedure
for the Boltzmann equation through the moment equations. The method is
referred to as Grad’s moment method and derives the governing equations for
the macroscopic quantities, namely mass density, velocity, temperature,
stress, heat flux and other higher-order moments, through the Boltzmann
equation since these macroscopic quantities are related to the distribution
function via its moments. The well-known Grad 13-moment equations, at least in
the linearised case, are accompanied with a proper entropy inequality (Rana &
Struchtrup, 2016) and are stable.

Over the last few decades, a significant progress has been made towards
improving Grad’s approach while keeping their advantageous features (see e.g.
Karlin et al., 1998; Ohr, 2001; Struchtrup & Torrilhon, 2003; Struchtrup,
2004; Cai et al., 2014; Karlin, 2018; Cai & Wang, 2020). In particular,
Struchtrup & Torrilhon (2003) proposed a regularisation procedure that derives
a closed set of moment equations—referred to as the regularised 13-moment
(R13) equations. The R13 equations, since their derivation, have paved the way
to several developments, including numerical frameworks (see e.g. Koellermeier
et al., 2014; Claydon et al., 2017; Cai & Wang, 2020). It can be shown via the
order of magnitude approach (Struchtrup, 2004) that the R13 equations possess
the accuracy of order $\mathrm{Kn}^{3}$ whereas the Euler, NSF and Grad
13-moment equations have the accuracy of order $\mathrm{Kn}^{0}$,
$\mathrm{Kn}^{1}$ and $\mathrm{Kn}^{2}$, respectively (Struchtrup, 2004, 2005;
Torrilhon, 2016). It is generally accepted that to maintain accuracy, the
number of moments ought to be increased with the Knudsen number. Aiming to
increase the applicability of the regularisation procedure to describe
processes in the moderate transition regime, Gu & Emerson (2009), by following
the approach of Struchtrup & Torrilhon (2003), derived the regularised
26-moment (R26) equations. The R26 equations are anticipated to possess the
accuracy of order $\mathrm{Kn}^{5}$ in the bulk region (i.e. outside of the
Knudsen layer) (Struchtrup, 2004, 2005).

Boundary conditions complementing the Grad or regularised moment equations are
typically obtained by employing the Maxwell accommodation model (Maxwell,
1879) along with Grad’s closure at the boundary yielding the so-called
macroscopic boundary conditions (MBC). The paper by Gu & Emerson (2007) may be
regarded as the pioneering work on wall boundary conditions for the R13
equations using the Maxwell accommodation model. Nevertheless, some
inconsistencies and oscillations were reported in this work. To overcome these
shortcomings, Torrilhon & Struchtrup (2008) proposed a theory of boundary
conditions for the R13 equations using the Maxwell accommodation model but
taking physical and mathematical requirements for the system into account.
Moreover, following the approach of Torrilhon & Struchtrup (2008), the
boundary conditions for the R26 equations—based on the Maxwell accommodation
model—were obtained in Gu & Emerson (2009). The boundary conditions derived
through the Maxwell accommodation model are generically referred to as the
Maxwell boundary conditions. Solutions of the R13 and R26 equations along with
the Maxwell boundary conditions have been obtained for some benchmark problems
(Taheri et al., 2009; Struchtrup & Taheri, 2011; Gu et al., 2010; Gu &
Emerson, 2014; Claydon et al., 2017; Gu et al., 2019) in order to investigate
rarefaction effects in gas at moderate Knudsen numbers. Recently, a hybrid
approach combining the discrete velocity method with the MBC in the near-wall
region, where the moment method is expected to be less accurate due to strong
nonequilibrium effects, has been proposed by Yang et al. (2020). The approach
not only yields more accurate solutions for wall-bounded flows but also
reduces the computational cost in numerical simulations.

Thus far, analytical solutions for evaporation and condensation problems have
been obtained from the linearised R13 (LR13) (Beckmann et al., 2018) and
linearised R26 (LR26) (Rana et al., 2018b) equations and compared with
numerical solutions from kinetic theory. The R26 equations can capture Knudsen
layers—a region of strong non-equilibrium within a few mean free paths of the
boundary—more accurately than the R13 equations giving a good agreement with
kinetic theory for the Knudsen number up to 1. The moment equations are now
being implemented into state-of-the-art open-source softwares (Torrilhon &
Sarna, 2017; Westerkamp & Torrilhon, 2019), providing a numerical platform to
solve these equations for engineering applications.

Despite the success of the moment equations, there have been concerns over the
accuracy of MBC obtained using Grad’s closure at the boundary. Rana &
Struchtrup (2016) showed that the MBC for the LR13 equations violate the
Onsager reciprocity relations (Onsager, 1931a, b), and derived a new set of
phenomenological boundary conditions (PBC) for them. For PBC, one writes the
entropy generation at the boundary and finds the boundary conditions as
phenomenological laws that guarantee positivity of the entropy inequality at
the boundary, which in turn leads to stable numerical schemes and convergence
of moments methods (Ringhofer et al., 2001; Sarna & Torrilhon, 2018a). The PBC
for the LR13 equations are thermodynamically admissible for all processes and
comply with the Onsager reciprocity relations (Rana & Struchtrup, 2016).

In this paper, we first determine a quadratic form without cross-product terms
for the entropy density to prove the $H$-theorem (entropy inequality) for the
LR26 equations, which in turn leads to the PBC for them. We then solve the
LR26 equations along with the associated PBC to investigate a slow flow of a
gas/vapour past an evaporating/non-evaporating spherical droplet analytically,
and compare the results with those obtained from other continuum theories,
existing experimental data and/or numerical results obtained from the
linearised Boltzmann equation. Moreover, we estimate the drag force acting on
the evaporating/non-evaporating spherical droplet with the LR26 equations. The
results not only give an excellent match but—thanks to the analytic
solutions—also provide a deeper understanding of the process involved. Through
the examples considered in this paper, the present approach evidently provides
the best ever results obtained via macroscopic models.

The remainder of the paper is organised as follows. The LR26 equations in the
dimensionless form are presented in § 2. The $H$-theorem for the LR26
equations is proved in § 3. Physically admissible boundary conditions for the
LR26 equations are derived in § 4, where different contributions to the
entropy production rate at the interface are also expounded by exploiting the
Curie principle and Onsager symmetry relations. Canonical boundary value
problems are solved with the LR26 equations with the PBC and MBC in § 5 to
assess the Onsager reciprocity relations and entropy generation at the
interface. The problem of a steady gas flow past a spherical droplet with and
without evaporation is investigated in § 6. Concluding remarks are made in §
7.

##  2 The linearised R26 equations

The treatment developed in this paper is based entirely on the linearised
equations, which result from omitting the nonlinear terms (inertial terms,
viscous heating terms, for instance) from the governing equations. This
simplification is justified since many micro- and nano-fluid devices involve
sufficiently slow flows. For succinctness, we shall present only the LR26
equations (and subsystems), which can be obtained from the original R26
equations derived in Gu & Emerson (2009) in a straightforward way. We shall
consider small deviations from a constant equilibrium state—given by a
constant reference density $\hat{\rho}_{0}$, a constant reference temperature
$\hat{T}_{0}$ and all other fields as zero. Thereby, all the equations can be
linearised with respect to the reference state. Moreover, the equations will
also be made dimensionless by introducing a length scale $\hat{L}_{0}$, a
velocity scale $\sqrt{\hat{R}\hat{T}_{0}}$, a time scale
$\hat{L}_{0}/\sqrt{\hat{R}\hat{T}_{0}}$ and the dimensionless deviations in
field variables from their respective reference states

$\displaystyle\left.\begin{gathered}\rho:=\frac{\hat{\rho}-\hat{\rho}_{0}}{\hat{\rho}_{0}},\quad
v_{i}:=\frac{\hat{v}_{i}}{\sqrt{\hat{R}\hat{T}_{0}}},\quad
T:=\frac{\hat{T}-\hat{T}_{0}}{\hat{T}_{0}},\\\
\sigma_{ij}:=\frac{\hat{\sigma}_{ij}}{\hat{\rho}_{0}\hat{R}\hat{T}_{0}},\quad
q_{i}:=\frac{\hat{q}_{i}}{\hat{\rho}_{0}\big{(}\hat{R}\hat{T}_{0}\big{)}^{3/2}},\\\
m_{ijk}:=\frac{\hat{m}_{ijk}}{\hat{\rho}_{0}\big{(}\hat{R}\hat{T}_{0}\big{)}^{3/2}},\quad
R_{ij}:=\frac{\hat{R}_{ij}}{\hat{\rho}_{0}\big{(}\hat{R}\hat{T}_{0}\big{)}^{2}},\quad\Delta:=\frac{\hat{\Delta}}{\hat{\rho}_{0}\big{(}\hat{R}\hat{T}_{0}\big{)}^{2}},\\\
\Phi_{ijkl}:=\frac{\hat{\Phi}_{ijkl}}{\hat{\rho}_{0}\big{(}\hat{R}\hat{T}_{0}\big{)}^{2}},\quad\Psi_{ijk}:=\frac{\hat{\Psi}_{ijk}}{\hat{\rho}_{0}\big{(}\hat{R}\hat{T}_{0}\big{)}^{5/2}},\quad\Omega_{i}:=\frac{\hat{\Omega}_{i}}{\hat{\rho}_{0}\big{(}\hat{R}\hat{T}_{0}\big{)}^{5/2}},\end{gathered}\right\\}$
(5)

where $\hat{R}$ is the gas constant; $\rho$, $v_{i}$, $T$, $\sigma_{ij}$ and
$q_{i}$ are the dimensionless deviations in the density, macroscopic velocity,
temperature, stress tensor and heat flux from their respective reference
state; $m_{ijk}$, $R_{ij}$, $\Delta$, $\Phi_{ijkl}$, $\Psi_{ijk}$ and
$\Omega_{i}$ are dimensionless deviations in the higher-order moments. Note
that the notation for the field variables are the same as those in Gu &
Emerson (2009) except for the quantities with hats here are with dimensions.
As a consequence of the linearisation, the dimensionless deviation in the
pressure $p$ is given by linearised equation of state, i.e. $p=\rho+T$. The
non-dimensionalisation also introduces a dimensionless parameter—the Knudsen
number

$\displaystyle\mathrm{Kn}=\frac{\hat{\mu}_{0}}{\hat{\rho}_{0}\sqrt{\hat{R}\hat{T}_{0}}L_{0}}$
(6)

where $\hat{\mu}_{0}$ is the viscosity of the gas at the reference state.
Accordingly, the dimensionless and linearised conservation laws for the mass,
momentum and energy read

$\displaystyle\frac{\partial\rho}{\partial t}+\frac{\partial v_{i}}{\partial
x_{i}}$ $\displaystyle=0,$ (7) $\displaystyle\frac{\partial v_{i}}{\partial
t}+\frac{\partial\sigma_{ij}}{\partial x_{j}}+\frac{\partial p}{\partial
x_{i}}$ $\displaystyle=F_{i},$ (8) $\displaystyle\frac{3}{2}\frac{\partial
T}{\partial t}+\frac{\partial q_{i}}{\partial x_{i}}+\frac{\partial
v_{i}}{\partial x_{i}}$ $\displaystyle=0,$ (9)

where $F_{i}$ is the dimensionless external force, $\partial(\cdot)/\partial
t$ denotes the dimensionless time derivative and $\partial(\cdot)/\partial
x_{i}$ denotes the dimensionless spatial derivative. In (7)–(9) and in what
follows, indices denote vector and tensor components in Cartesian coordinates,
and the Einstein summation convention is assumed over the repeated indices
unless stated otherwise.

The higher moment system also includes the governing equations for the stress
and heat flux, which in the linear-dimensionless form read

$\displaystyle\frac{\partial\sigma_{ij}}{\partial t}+\frac{\partial
m_{ijk}}{\partial x_{k}}+\frac{4}{5}\frac{\partial q_{\langle i}}{\partial
x_{j\rangle}}$
$\displaystyle=-\frac{1}{\mathrm{Kn}}\left[\sigma_{ij}+2\mathrm{Kn}\frac{\partial
v_{\langle i}}{\partial x_{j\rangle}}\right],$ (10)
$\displaystyle\frac{\partial q_{i}}{\partial t}+\frac{1}{2}\frac{\partial
R_{ij}}{\partial x_{j}}+\frac{1}{6}\frac{\partial\Delta}{\partial
x_{i}}+\frac{\partial\sigma_{ij}}{\partial x_{j}}$
$\displaystyle=-\frac{\mathrm{Pr}}{\mathrm{Kn}}\left[q_{i}+\frac{5}{2}\frac{\mathrm{Kn}}{\mathrm{Pr}}\frac{\partial
T}{\partial x_{i}}\right],$ (11)

where $\mathrm{Pr}$ is the Prandtl number and the angular brackets around
indices denote the symmetric-tracefree part of a tensor (Struchtrup, 2005).
The right-hand sides of these balance equations contain NSF laws, which can
also be obtained at the first order of a CE-like expansion performed on the
stress and heat flux balance equations ((10) and (11)) in powers of the
Knudsen number, see e.g. Torrilhon & Struchtrup (2004); Struchtrup (2005). The
first order of the CE-like expansion on these equations essentially amounts to
setting their left-hand sides to zero, yielding the NSF laws.

The balance equations (7)–(11) form the so-called 13-moment system. However,
this system is not closed since they contain additional higher moments
$m_{ijk}$, $R_{ij}$ and $\Delta$, which are defined as

$\displaystyle m_{ijk}$ $\displaystyle=\int_{\mathbb{R}^{3}}C_{\langle
i}C_{j}C_{k\rangle}f\,\mathrm{d}\bm{c},$ (12) $\displaystyle R_{ij}$
$\displaystyle=\int_{\mathbb{R}^{3}}C^{2}C_{\langle
i}C_{j\rangle}f\,\mathrm{d}\bm{c}-7\sigma_{ij},$ (13) $\displaystyle\Delta$
$\displaystyle=\int_{\mathbb{R}^{3}}C^{4}f\,\mathrm{d}\bm{c}-15,$ (14)

where $C_{i}=c_{i}-v_{i}$ is the dimensionless peculiar velocity with $c_{i}$
being the dimensionless microscopic velocity. These higher moments are
constructed in such a way that they vanish on being computed with the Grad
13-moment (G13) distribution function, i.e.
$m_{ijk|\textrm{G13}}=R_{ij|\textrm{G13}}=\Delta_{|\textrm{G13}}=0$.

The (linearised) 26-moment theory is comprised of the conservation laws
(7)–(9), the stress balance equation (10), the heat flux balance equation (11)
and the balance equations for $m_{ijk}$, $R_{ij}$ and $\Delta$. The balance
equations for $\hat{m}_{ijk}$, $\hat{R}_{ij}$ and $\hat{\Delta}$ have already
been obtained from the Boltzmann equation by Gu & Emerson (2009); in the
linear-dimensionless form, they read

$\displaystyle\frac{\partial m_{ijk}}{\partial
t}+\frac{\partial\Phi_{ijkl}}{\partial x_{l}}+\frac{3}{7}\frac{\partial
R_{\langle ij}}{\partial x_{k\rangle}}$
$\displaystyle=-\frac{\mathrm{Pr}_{m}}{\mathrm{Kn}}\left[m_{ijk}+3\frac{\mathrm{Kn}}{\mathrm{Pr}_{m}}\frac{\partial\sigma_{\langle
ij}}{\partial x_{k\rangle}}\right]$ (15) $\displaystyle\frac{\partial
R_{ij}}{\partial t}+\frac{\partial\Psi_{ijk}}{\partial x_{k}}+2\frac{\partial
m_{ijk}}{\partial x_{k}}+\frac{2}{5}\frac{\partial\Omega_{\langle i}}{\partial
x_{j\rangle}}$
$\displaystyle=-\frac{\mathrm{Pr}_{R}}{\mathrm{Kn}}\left[R_{ij}+\frac{28}{5}\frac{\mathrm{Kn}}{\mathrm{Pr}_{R}}\frac{\partial
q_{\langle i}}{\partial x_{j\rangle}}\right]$ (16)
$\displaystyle\frac{\partial\Delta}{\partial
t}+\frac{\partial\Omega_{i}}{\partial x_{i}}$
$\displaystyle=-\frac{\mathrm{Pr}_{\Delta}}{\mathrm{Kn}}\left[\Delta+8\frac{\mathrm{Kn}}{\mathrm{Pr}_{\Delta}}\frac{\partial
q_{i}}{\partial x_{i}}\right].$ (17)

The coefficients $\mathrm{Pr}$, $\mathrm{Pr}_{m}$, $\mathrm{Pr}_{R}$ and
$\mathrm{Pr}_{\Delta}$ in the above equations depend on the choice of
intermolecular potential function appearing in the Boltzmann collision
operator. For Maxwell molecules (MM), these parameters read:
$\mathrm{Pr}=2/3$, $\mathrm{Pr}_{m}=3/2$, $\mathrm{Pr}_{R}=7/6$ and
$\mathrm{Pr}_{\Delta}=2/3$ (Gu & Emerson, 2009).

For a third-order accuracy in $\mathrm{Kn}$, it is sufficient to consider only
the terms on the right-hand sides of (15)–(17) and dropping all the terms on
the left-hand sides of these equations; this gives the linear R13 constitutive
relations (Struchtrup & Torrilhon, 2003; Struchtrup, 2005). Formally, the R13
constitutive relations stem from (15)–(17) on performing a CE-like expansion
on these equations in powers of $\mathrm{Kn}$, i.e. by setting the left-hand
sides of (15)–(17) to zero (Struchtrup, 2004). Thus, the (dimensionless) LR13
equations consist of (7)–(11) closed with

$\displaystyle
m_{ijk}=-3\frac{\mathrm{Kn}}{\mathrm{Pr}_{m}}\frac{\partial\sigma_{\langle
ij}}{\partial x_{k\rangle}},\quad
R_{ij}=-\frac{28}{5}\frac{\mathrm{Kn}}{\mathrm{Pr}_{R}}\frac{\partial
q_{\langle i}}{\partial
x_{j\rangle}}\quad\textrm{and}\quad\Delta=-8\frac{\mathrm{Kn}}{\mathrm{Pr}_{\Delta}}\frac{\partial
q_{i}}{\partial x_{i}}.$ (18)

###  2.1 The R26 constitutive relations

The system of 26 moment equations (7)–(11) and (15)–(17) contains the
additional moments $\Phi_{ijkl}$, $\Psi_{ijk}$ and $\Omega_{i}$:

$\displaystyle\Phi_{ijkl}$ $\displaystyle=\int_{\mathbb{R}^{3}}C_{\langle
i}C_{j}C_{k}C_{l\rangle}f\,\mathrm{d}\bm{c},$ (19) $\displaystyle\Psi_{ijk}$
$\displaystyle=\int_{\mathbb{R}^{3}}C^{2}C_{\langle
i}C_{j}C_{k\rangle}f\,\mathrm{d}\bm{c}-9m_{ijk},$ (20)
$\displaystyle\Omega_{i}$
$\displaystyle=\int_{\mathbb{R}^{3}}C^{4}C_{i}f\,\mathrm{d}\bm{c}-28q_{i},$
(21)

which need to be specified in order to close the system of the 26 moment
equations. The R26 constitutive relations, which provides a closure for the 26
moment equations, follow from a CE-like expansion on the evolution equations
for the moments $\Phi_{ijkl}$, $\Psi_{ijk}$ and $\Omega_{i}$ in a 45-moment
theory (see Gu & Emerson, 2009). The R26 constitutive relations read (Gu &
Emerson, 2009)

$\displaystyle\left.\begin{aligned}
\Phi_{ijkl}&=-4\frac{\mathrm{Kn}}{\mathrm{Pr}_{\Phi}}\frac{\partial m_{\langle
ijk}}{\partial x_{l\rangle}},\\\
\Psi_{ijk}&=-\frac{27}{7}\frac{\mathrm{Kn}}{\mathrm{Pr}_{\Psi}}\frac{\partial
R_{\langle ij}}{\partial x_{k\rangle}},\\\
\Omega_{i}&=-\frac{7}{3}\frac{\mathrm{Kn}}{\mathrm{Pr}_{\Omega}}\left(\frac{\partial\Delta}{\partial
x_{i}}+\frac{12}{7}\frac{\partial R_{ij}}{\partial
x_{j}}\right),\end{aligned}\right\\}$ (22)

where the values of the parameters are $\mathrm{Pr}_{\Phi}=2.097$,
$\mathrm{Pr}_{\Psi}=1.698$ and $\mathrm{Pr}_{\Omega}=1$ for MM (Gu & Emerson,
2009). The system of 26 moment equations (7)–(11) and (15)–(17) along with the
constitutive relations (22) form the system of the LR26 equations.

##  3 $H$-theorem for the linearised R26 equations

The $H$-theorem plays a salient role in constructing constitutive relations in
the bulk (de Groot & Mazur, 1962; Struchtrup & Torrilhon, 2007) and in
developing physically admissible boundary conditions (Struchtrup & Torrilhon,
2007; Rana & Struchtrup, 2016; Rana et al., 2018a; Sarna & Torrilhon, 2018a).
Recently, Struchtrup & Nadler (2020) deduced that the Burnett equations, which
do not comply with the $H$-theorem, lead to thermodynamically inadmissible
solutions (producing work without a driving force). In this section, we
shall—for the first time—show that the LR26 equations are accompanied with an
entropy law given as

$\displaystyle\frac{\partial\eta}{\partial
t}+\frac{\partial\Gamma_{i}}{\partial x_{i}}=\Sigma\geqslant 0,$ (23)

where $\eta$ denotes the dimensionless entropy density, $\Gamma_{i}$ the
dimensionless entropy flux and $\Sigma$ the dimensionless non-negative entropy
generation rate.

The Boltzmann entropy density $\hat{\eta}$ and entropy flux $\hat{\Gamma}_{i}$
(in dimensional form) are given by (Müller & Ruggeri, 1998)

$\displaystyle\hat{\eta}$
$\displaystyle=-\hat{k}_{B}\int\hat{f}\left(\ln\frac{\hat{f}}{\hat{y}}-1\right)\mathrm{d}\hat{\bm{c}},$
(24) $\displaystyle\hat{\Gamma}_{i}$
$\displaystyle=-\hat{k}_{B}\int\hat{f}\left(\ln\frac{\hat{f}}{\hat{y}}-1\right)\hat{c}_{i}\,\mathrm{d}\hat{\bm{c}},$
(25)

where $\hat{k}_{B}$ is the Boltzmann constant, $\hat{y}$ is a constant having
dimensions of the distribution function $\hat{f}$ and $\hat{\bm{c}}$ is the
microscopic velocity (with dimensions). Starting from the Boltzmann entropy
density (24), it can be shown that the entropy density of a system of the
linearised moment equations—obtained from the Boltzmann equation—is a
quadratic form without cross-product terms (see e.g., Struchtrup & Torrilhon,
2007; Rana & Struchtrup, 2016; Sarna & Torrilhon, 2018b; Beckmann et al.,
2018; Rana et al., 2018a). It is worthwhile noting that the hyperbolic part of
the LR26 equations is the same as the linearised Grad 26-moment (G26)
equations (in the same way the hyperbolic part of the linearised NSF equations
is the same as the linearised Euler equations); therefore, it is sufficient to
use the G26 distribution function for computing the entropy density for the
LR26 equations. Using the G26 distribution function for $\hat{f}$ in (24) and
(25), we show that entropy density for the LR26 equations is a quadratic form
without cross-product terms and that the entropy flux obtained with the G26
distribution function contains all the hyperbolic terms present in the
expression of the entropy flux for the LR26 equations, although the
calculations are relegated to appendix A for better readability. The
dimensionless entropy density for the LR26 equations is given by (see
appendices A and B)

$\displaystyle\eta=\underline{a_{0}-\frac{1}{2}\rho^{2}-\frac{1}{2}v^{2}-\frac{3}{4}T^{2}-\frac{1}{4}\sigma^{2}-\frac{1}{5}q^{2}}-\frac{1}{12}m^{2}-\frac{1}{56}R^{2}-\frac{1}{240}\Delta^{2},$
(26)

where $a_{0}$ is just a constant and the notation $A^{2}:=A_{i_{1}i_{2}\dots
i_{k}}A_{i_{1}i_{2}\dots i_{k}}$ has been used to denote the contraction of a
$k$-rank tensor $A_{i_{1}i_{2}\dots i_{k}}$ with itself for brevity. The
entropy for the LR26 equations (26) comprises of the contribution from the
13-moment theory (Rana & Struchtrup, 2016)—given by the underlined terms in
(26)—and additional contributions from the higher-order moments.

The dimensionless entropy density (26) can also be obtained in a more
comprehensible way by assuming it to be a quadratic form without cross-product
terms. The approach not only leads to the dimensionless entropy density (26)
but also yields the complete entropy flux and entropy production rate for the
LR26 equations alongside. The approach assumes the following form of
dimensionless entropy density:

$\displaystyle\eta=a_{0}+\frac{a_{1}}{2}\rho^{2}+\frac{a_{2}}{2}v^{2}+\frac{a_{3}}{2}T^{2}+\frac{a_{4}}{2}\sigma^{2}+\frac{a_{5}}{2}q^{2}+\frac{a_{6}}{2}m^{2}+\frac{a_{7}}{2}R^{2}+\frac{a_{8}}{2}\Delta^{2}.$
(27)

Here $a_{0},a_{1},a_{2},\dots,a_{8}$ are unknown coefficients, which are
computed by taking the time derivative of the entropy density (27), then
eliminating the time derivatives of the moments by means of the G26 equations,
and finally comparing the resulting equation with (23). It turns out that one
of the variables from $a_{1},a_{2},\dots,a_{8}$ can be chosen arbitrarily. On
taking $a_{1}=-1$, the other coefficients $a_{2},a_{3},\dots,a_{8}$ turn out
to be such that (27) is the same as (26); see appendix B for details.

For the computation of the entropy flux and entropy production rate for the
LR26 equations, we therefore can start with the entropy density (26). Taking
the time derivative of $\eta$ in (26), and substituting the time derivatives
of the field variables $\rho$, $v_{i}$, $T$ from (7)–(9), the time derivative
of $\sigma_{ij}$ from (10), the time derivative of $q_{i}$ from (11), and the
time derivatives of $m_{ijk}$, $R_{ij}$, $\Delta$ from (15)–(17), one obtains

$\displaystyle\frac{\partial\eta}{\partial t}=$
$\displaystyle\frac{\partial}{\partial
x_{i}}\left(pv_{i}+Tq_{i}+v_{j}\sigma_{ij}+\frac{2}{5}\sigma_{ij}q_{j}+\frac{1}{2}m_{ijk}\sigma_{jk}+\frac{1}{5}R_{ij}q_{j}+\frac{1}{15}q_{i}\Delta\right.$
$\displaystyle\left.+\frac{1}{14}m_{ijk}R_{jk}+\frac{1}{6}\Phi_{ijkl}m_{jkl}+\frac{1}{28}\Psi_{ijk}R_{jk}+\frac{1}{70}R_{ij}\Omega_{j}+\frac{1}{120}\Delta\Omega_{i}\right)$
$\displaystyle-\frac{1}{6}\Phi_{ijkl}\frac{\partial m_{jkl}}{\partial
x_{i}}-\frac{1}{28}\Psi_{ijk}\frac{\partial R_{jk}}{\partial
x_{i}}-\frac{1}{70}\Omega_{j}\frac{\partial R_{ij}}{\partial
x_{i}}-\frac{1}{120}\Omega_{i}\frac{\partial\Delta}{\partial x_{i}}$
$\displaystyle+\frac{1}{2}\frac{1}{\mathrm{Kn}}\sigma^{2}+\frac{2}{5}\frac{\mathrm{Pr}}{\mathrm{Kn}}q^{2}+\frac{1}{6}\frac{\mathrm{Pr}_{m}}{\mathrm{Kn}}m^{2}+\frac{1}{28}\frac{\mathrm{Pr}_{R}}{\mathrm{Kn}}R^{2}+\frac{1}{120}\frac{\mathrm{Pr}_{\Delta}}{\mathrm{Kn}}\Delta^{2}.$
(28)

Comparison of (3) with (23) gives the expression for the entropy flux as

$\displaystyle\Gamma_{i}$
$\displaystyle=-\left(pv_{i}+Tq_{i}+v_{j}\sigma_{ij}+\frac{2}{5}\sigma_{ij}q_{j}+\frac{1}{2}m_{ijk}\sigma_{jk}+\frac{1}{5}R_{ij}q_{j}+\frac{1}{15}q_{i}\Delta\right.$
$\displaystyle\left.\quad+\frac{1}{14}m_{ijk}R_{jk}+\frac{1}{6}\Phi_{ijkl}m_{jkl}+\frac{1}{28}\Psi_{ijk}R_{jk}+\frac{1}{70}R_{ij}\Omega_{j}+\frac{1}{120}\Delta\Omega_{i}\right),$
(29)

and the entropy production rate as

$\displaystyle\Sigma$ $\displaystyle=-\frac{1}{6}\Phi_{ijkl}\frac{\partial
m_{jkl}}{\partial x_{i}}-\frac{1}{28}\Psi_{ijk}\frac{\partial R_{jk}}{\partial
x_{i}}-\frac{1}{70}\Omega_{j}\frac{\partial R_{ij}}{\partial
x_{i}}-\frac{1}{120}\Omega_{i}\frac{\partial\Delta}{\partial x_{i}}$
$\displaystyle\quad+\frac{1}{2}\frac{1}{\mathrm{Kn}}\sigma^{2}+\frac{2}{5}\frac{\mathrm{Pr}}{\mathrm{Kn}}q^{2}+\frac{1}{6}\frac{\mathrm{Pr}_{m}}{\mathrm{Kn}}m^{2}+\frac{1}{28}\frac{\mathrm{Pr}_{R}}{\mathrm{Kn}}R^{2}+\frac{1}{120}\frac{\mathrm{Pr}_{\Delta}}{\mathrm{Kn}}\Delta^{2}.$
(30)

Using the linearised R26 constitutive relations (22) in (3), the entropy
production rate $\Sigma$ can be written as

$\displaystyle\Sigma$
$\displaystyle=\frac{1}{24}\frac{\mathrm{Kn}}{\mathrm{Pr}_{\Phi}}\Phi^{2}+\frac{1}{108}\frac{\mathrm{Kn}}{\mathrm{Pr}_{\Psi}}\Psi^{2}+\frac{1}{280}\frac{\mathrm{Kn}}{\mathrm{Pr}_{\Omega}}\Omega^{2}$
$\displaystyle\quad+\frac{1}{2}\frac{1}{\mathrm{Kn}}\sigma^{2}+\frac{2}{5}\frac{\mathrm{Pr}}{\mathrm{Kn}}q^{2}+\frac{1}{6}\frac{\mathrm{Pr}_{m}}{\mathrm{Kn}}m^{2}+\frac{1}{28}\frac{\mathrm{Pr}_{R}}{\mathrm{Kn}}R^{2}+\frac{1}{120}\frac{\mathrm{Pr}_{\Delta}}{\mathrm{Kn}}\Delta^{2},$
(31)

which is non-negative for $\mathrm{Pr}_{\Psi}$, $\mathrm{Pr}_{\Phi}$,
$\mathrm{Pr}_{\Omega}$, $\mathrm{Pr}_{m}$, $\mathrm{Pr}_{R}$,
$\mathrm{Pr}_{\Delta}$, $\Pr$, $\mathrm{Kn}\geq 0$. This completes the proof
of the $H$-theorem for the linearised R26 equations.

##  4 Phenomenological boundary conditions for the LR26 equations

In this section, we shall demonstrate how the $H$-theorem can be used to
generate a proper and complete set of boundary conditions for the LR26
equations. We shall derive the PBC for the LR26 equations from the entropy
production rate at the interface exploiting the force-flux relationships.

Let there be two phases, namely the liquid (or solid) and the gaseous (or
vapour), connected by a massless interface $I$ of negligible thickness and let
there be no surface tension and surface energy at the interface. The entropy
production rate at such an interface $\Sigma^{I}$ is given by the difference
between the entropy fluxes into and out of an interface (de Groot & Mazur,
1962; Rana & Struchtrup, 2016; Beckmann et al., 2018), i.e.

$\displaystyle\Sigma^{I}=\left[\Gamma_{i}^{\mathrm{(gas)}}-\Gamma_{i}^{\mathrm{(liquid)}}\right]n_{i}\geq
0.$ (32)

where $n_{i}$ is the unit normal to the interface pointing from the liquid
into the gas. Assuming that the gaseous (or vapour) phase is comprised of an
ideal gas, the entropy flux of the gas (or vapour) medium is
$\Gamma_{i}^{\mathrm{(gas)}}=\Gamma_{i}$, given by (3), and assuming that the
liquid phase is comprised of an incompressible liquid and can be described by
the NSF equations, the entropy flux of the liquid phase is (Beckmann et al.,
2018)

$\displaystyle\Gamma_{i}^{\mathrm{(liquid)}}=-p^{\ell}v_{i}^{\ell}-T^{\ell}q_{i}^{\ell}-v_{j}^{\ell}\sigma_{ij}^{\ell},$
(33)

where the superscript ‘$\ell$’ on the variables are used to denote that they
are the properties of the liquid. Substituting
$\Gamma_{i}^{\mathrm{(gas)}}=\Gamma_{i}$ from (3) and
$\Gamma_{i}^{\mathrm{(liquid)}}$ from (33) into (32), the entropy flux at the
interface—after some algebra—turns out to be (see appendix C for details)

$\displaystyle\Sigma^{I}$
$\displaystyle=-\left(\mathcal{P}\mathcal{V}_{i}+\mathcal{T}q_{i}+\mathscr{V}_{j}\sigma_{ij}+\frac{2}{5}\sigma_{ij}q_{j}+\frac{1}{2}m_{ijk}\sigma_{jk}+\frac{1}{5}R_{ij}q_{j}+\frac{1}{15}q_{i}\Delta\right.$
$\displaystyle\left.\quad\,+\frac{1}{14}m_{ijk}R_{jk}+\frac{1}{6}\Phi_{ijkl}m_{jkl}+\frac{1}{28}\Psi_{ijk}R_{jk}+\frac{1}{70}R_{ij}\Omega_{j}+\frac{1}{120}\Omega_{i}\Delta\right)n_{i},$
(34)

where

$\displaystyle\mathcal{P}=p-p_{\mathrm{sat}},\quad\mathcal{T}=T-T^{\ell},\quad\mathcal{V}_{i}=v_{i}-v_{i}^{I}\quad\text{and}\quad\mathscr{V}_{j}=v_{j}-v_{j}^{\ell}$
(35)

with $v_{i}^{I}$ being the velocity of the interface and
$p_{\mathrm{sat}}\equiv p_{\mathrm{sat}}(T^{\ell})$ being the saturation
pressure corresponding to the temperature $T^{\ell}$.

For further simplification, it is imperative to decompose the vectors and
tensors into their components in the normal and tangential directions. Such a
decomposition for symmetric-tracefree tensors of rank up to three is given in
appendix D. This decomposition allows one to write the entropy production at
the wall (4) as a sum of the product of the fluxes with an odd degree in $n$
(i.e. $\mathcal{V}_{n}$, $q_{n}$, $\Omega_{n}$, $m_{nnn}$, $\Psi_{nnn}$,
$\bar{\sigma}_{ni}$, $\bar{R}_{ni}$, $\bar{\Phi}_{nnni}$, $\tilde{m}_{nij}$,
$\tilde{\Psi}_{nij}$, $\check{\Phi}_{nijk}$) and the moments with an even
degree in $n$ as

$\displaystyle\Sigma^{I}$
$\displaystyle=-\mathcal{V}_{n}(\mathcal{P}+\sigma_{nn})-q_{n}\left(\mathcal{T}+\frac{2}{5}\sigma_{nn}+\frac{1}{5}R_{nn}+\frac{1}{15}\Delta\right)$
$\displaystyle\quad-
m_{nnn}\left(\frac{3}{4}\sigma_{nn}+\frac{3}{28}R_{nn}+\frac{5}{12}\Phi_{nnnn}\right)-\Psi_{nnn}\left(\frac{3}{56}R_{nn}\right)$
$\displaystyle\quad-\Omega_{n}\left(\frac{1}{120}\Delta+\frac{1}{70}R_{nn}\right)-\bar{\sigma}_{ni}\left(\bar{\mathscr{V}}_{i}+\frac{2}{5}\bar{q}_{i}+\bar{m}_{nni}\right)$
$\displaystyle\quad-\bar{R}_{ni}\left(\frac{1}{5}\bar{q}_{i}+\frac{1}{7}\bar{m}_{nni}+\frac{1}{14}\bar{\Psi}_{nni}+\frac{1}{70}\bar{\Omega}_{i}\right)-\bar{\Phi}_{nnni}\left(\frac{5}{8}\bar{m}_{nni}\right)$
$\displaystyle\quad-\tilde{m}_{nij}\left(\frac{1}{2}\tilde{\sigma}_{ij}+\frac{1}{14}\tilde{R}_{ij}+\frac{1}{2}\tilde{\Phi}_{nnij}\right)-\tilde{\Psi}_{nij}\left(\frac{1}{28}\tilde{R}_{ij}\right)-\check{\Phi}_{nijk}\left(\frac{1}{6}\check{m}_{ijk}\right).$
(36)

While writing (4), the relation
$\sigma_{nn}\mathscr{V}_{n}\approx\sigma_{nn}\mathcal{V}_{n}$ in light of
(110) and the Curie principle—which states that only forces and fluxes of the
same tensor type (scalars, vectors, 2-tensors, etc.) can be combined (de Groot
& Mazur, 1962)—have been exploited. In the above equation,
$\\{\mathcal{V}_{n}$, $q_{n}$, $m_{nnn}$, $\Psi_{nnn}$, $\Omega_{n}\\}$ are
the odd (in $n$) scalar fluxes, $\\{\bar{\sigma}_{ni}$, $\bar{R}_{ni}$,
$\bar{\Phi}_{nnni}\\}$ are the odd vector fluxes, $\\{\tilde{m}_{nij}$,
$\tilde{\Psi}_{nij}\\}$ are the odd tensor fluxes of rank two and
$\check{\Phi}_{nijk}$ is the odd tensor flux of rank three. The positive
entropy production at the interface can be attained by taking the odd fluxes
proportional to the even (in $n$) driving forces (written in brackets in
(4))—with proportionality constants being positive. This step to obtain the
positive entropy production at the interface yields the required boundary
conditions. Mathematically, the positive entropy production, and hence the
thermodynamically consistent PBC at the interface are achieved by writing a
linear force-flux relationship of the form

$\displaystyle\mathcal{J}_{r}=-\mathcal{L}_{rs}\mathcal{F}_{s}$ (37)

for scalar, vector and tensor fluxes of each rank, successively. Here,
$\mathcal{J}_{r}$ is the $r^{\mathrm{th}}$ component of a vector
$\bm{\mathcal{J}}$ containing the fluxes, $\mathcal{F}_{s}$ is the
$s^{\mathrm{th}}$ component of a vector $\bm{\mathcal{F}}$ containing the
forces, $\mathcal{L}_{rs}$ is the $(r,s)^{\mathrm{th}}$ element of an unknown
matrix $\bm{\mathcal{L}}$ and Einstein summation is assumed over $s$ in (37).
The unknown matrix $\bm{\mathcal{L}}$ consists of phenomenological
coefficients as its elements. The matrix $\bm{\mathcal{L}}$ must be symmetric
in order to comply with the Onsager reciprocity relations (Onsager, 1931a, b)
and must be positive semidefinite to ensure a non-negative entropy production.

###  4.1 Boundary conditions for the scalar fluxes

The scalar fluxes in (4) and their corresponding forces are collected in
$\bm{\mathcal{J}}$ and $\bm{\mathcal{F}}$ as

$\displaystyle\bm{\mathcal{J}}=\begin{bmatrix}\mathcal{V}_{n}\\\ q_{n}\\\
m_{nnn}\\\ \Psi_{nnn}\\\
\Omega_{n}\end{bmatrix}\qquad\textrm{and}\qquad\bm{\mathcal{F}}=\begin{bmatrix}\mathcal{P}+\sigma_{nn}\\\
\mathcal{T}+\frac{2}{5}\sigma_{nn}+\frac{1}{5}R_{nn}+\frac{1}{15}\Delta\\\
\frac{3}{4}\sigma_{nn}+\frac{3}{28}R_{nn}+\frac{5}{12}\Phi_{nnnn}\\\
\frac{3}{56}R_{nn}\\\ \frac{1}{120}\Delta+\frac{1}{70}R_{nn}\end{bmatrix}.$
(38)

Substituting $\bm{\mathcal{J}}$ and $\bm{\mathcal{F}}$ from (38) in (37), one
obtains the PBC for the scalar fluxes. Nevertheless, the matrix
$\bm{\mathcal{L}}$ is still unknown. The elements of the matrix
$\bm{\mathcal{L}}$ can be obtained by comparing the coefficients of
$\mathcal{P}$, $\mathcal{T}$, $\sigma_{nn}$, $R_{nn}$ and $\Delta$ in the
corresponding terms of the PBC and MBC for the scalar fluxes. The MBC for the
scalar fluxes (obtained through the Maxwell accommodation model) read (Rana et
al., 2018b)

$\displaystyle\left.\begin{aligned}
\mathcal{V}_{n}&=-\varsigma\left(\mathcal{P}-\frac{1}{2}\mathcal{T}+\frac{1}{2}\sigma_{nn}-\frac{1}{28}R_{nn}-\frac{1}{120}\Delta+b_{1}\Phi_{nnnn}\right),\\\
q_{n}&=-\frac{\mathcal{V}_{n}}{2}-\varkappa\left(2\mathcal{T}+\frac{1}{2}\sigma_{nn}+\frac{5}{28}R_{nn}+\frac{1}{15}\Delta+b_{2}\Phi_{nnnn}\right),\\\
m_{nnn}&=-\frac{2}{5}\mathcal{V}_{n}+\varkappa\left(\frac{2}{5}\mathcal{T}-\frac{7}{5}\sigma_{nn}-\frac{1}{14}R_{nn}+\frac{1}{75}\Delta+b_{3}\Phi_{nnnn}\right),\\\
\Psi_{nnn}&=\frac{6}{5}\mathcal{V}_{n}+\varkappa\left(\frac{6}{5}\mathcal{T}+\frac{9}{5}\sigma_{nn}-\frac{93}{70}R_{nn}+\frac{1}{5}\Delta+b_{4}\Phi_{nnnn}\right),\\\
\Omega_{n}&=3\mathcal{V}_{n}+\varkappa\left(8\mathcal{T}+2\sigma_{nn}-R_{nn}-\frac{4}{3}\Delta+b_{5}\Phi_{nnnn}\right),\end{aligned}\right\\}$
(39)

where

$\displaystyle b_{1}=-\frac{1}{24},\quad b_{2}=-\frac{1}{12},\quad
b_{3}=-\frac{13}{30},\quad b_{4}=-\frac{1}{2},\quad b_{5}=0,$ (40)
$\displaystyle\varsigma=\frac{\vartheta}{2-\vartheta}\sqrt{\frac{2}{\pi}}\qquad\textrm{and}\qquad\varkappa=\frac{\vartheta+\chi(1-\vartheta)}{2-\vartheta-\chi(1-\vartheta)}\sqrt{\frac{2}{\pi}}$
(41)

with $\vartheta$ being the evaporation/condensation coefficient and $\chi$
being the accommodation coefficient. Clearly, for $\vartheta=0$,
$\mathcal{V}_{n}=0$ from (39), hence one obtains the boundary conditions for a
non-evaporating boundary (Gu & Emerson, 2009). It is worthwhile noting that
the MBC (39) with (40) cannot be expressed in the force-flux formalism (37),
and we shall show later that this leads to violation of the Onsager
reciprocity relations and/or the second law at the interface.

Comparison of the coefficients of $\mathcal{P}$, $\mathcal{T}$, $\sigma_{nn}$,
$R_{nn}$ and $\Delta$ in the PBC for the scalar fluxes (obtained on
substituting $\bm{\mathcal{J}}$ and $\bm{\mathcal{F}}$ from (38) in (37)) and
MBC (39) yields

$\displaystyle\bm{\mathcal{L}}=\varsigma\begin{bmatrix}[r]1&-\frac{1}{2}&-\frac{2}{5}&\frac{6}{5}&3\\\
-\frac{1}{2}&\frac{1}{4}&\frac{1}{5}&-\frac{3}{5}&-\frac{3}{2}\\\
-\frac{2}{5}&\frac{1}{5}&\frac{4}{25}&-\frac{12}{25}&-\frac{6}{5}\\\
\frac{6}{5}&-\frac{3}{5}&-\frac{12}{25}&\frac{36}{25}&\frac{18}{5}\\\
3&-\frac{3}{2}&-\frac{6}{5}&\frac{18}{5}&9\end{bmatrix}+\varkappa\begin{bmatrix}[r]0&0&0&0&0\\\
0&2&-\frac{2}{5}&-\frac{6}{5}&-8\\\
0&-\frac{2}{5}&\frac{52}{25}&-\frac{44}{25}&\frac{8}{5}\\\
0&-\frac{6}{5}&-\frac{44}{25}&\frac{916}{25}&-\frac{72}{5}\\\
0&-8&\frac{8}{5}&-\frac{72}{5}&224\end{bmatrix}.$ (42)

Both the matrices on the right-hand side of (42) are symmetric and positive
semidefinite. Therefore, the matrix $\bm{\mathcal{L}}$ in (42) is also
symmetric and positive semidefinite as $\varsigma$ and $\varkappa$ are non-
negative. Thus, the PBC for the scalar fluxes are obtained from (37) with
$\bm{\mathcal{J}}$, $\bm{\mathcal{F}}$ and $\bm{\mathcal{L}}$ being given by
(38) and (42). On simplification, the PBC for the scalar fluxes are given by
(39) with

$\displaystyle b_{1}=-\frac{1}{6},\quad b_{2}=-\frac{1}{6},\quad
b_{3}=-\frac{13}{15},\quad b_{4}=\frac{11}{15},\quad b_{5}=-\frac{2}{3}.$ (43)

###  4.2 Boundary conditions for the vector fluxes

The vector fluxes in (4) and their corresponding forces are collected in
$\bm{\mathcal{J}}$ and $\bm{\mathcal{F}}$ as

$\displaystyle\bm{\mathcal{J}}=\begin{bmatrix}\bar{\sigma}_{ni}\\\
\bar{R}_{ni}\\\
\bar{\Phi}_{nnni}\end{bmatrix}\qquad\textrm{and}\qquad\bm{\mathcal{F}}=\begin{bmatrix}\bar{\mathscr{V}}_{i}+\frac{2}{5}\bar{q}_{i}+\bar{m}_{nni}\\\
\frac{1}{5}\bar{q}_{i}+\frac{1}{7}\bar{m}_{nni}+\frac{1}{14}\bar{\Psi}_{nni}+\frac{1}{70}\bar{\Omega}_{i}\\\
\frac{5}{8}\bar{m}_{nni}\end{bmatrix}.$ (44)

Substituting $\bm{\mathcal{J}}$ and $\bm{\mathcal{F}}$ from (44) in (37), one
obtains the PBC for the vector fluxes. Nevertheless, the matrix
$\bm{\mathcal{L}}$ is still unknown. The elements of the matrix
$\bm{\mathcal{L}}$ can be obtained by comparing the coefficients of
$\bar{\mathscr{V}}_{i}$, $\bar{q}_{i}$ and $\bar{m}_{nni}$ in the
corresponding terms of the PBC and MBC. The MBC for the vector fluxes
(obtained through the Maxwell accommodation model) read (Gu & Emerson, 2009;
Struchtrup et al., 2017)

$\displaystyle\left.\begin{aligned}
\bar{\sigma}_{ni}&=-\varkappa\left(\bar{\mathscr{V}}_{i}+\frac{1}{5}\bar{q}_{i}+\frac{1}{2}\bar{m}_{nni}+b_{6}\bar{\Psi}_{nni}+b_{7}\bar{\Omega}_{i}\right),\\\
\bar{R}_{ni}&=-\varkappa\left(-\bar{\mathscr{V}}_{i}+\frac{11}{5}\bar{q}_{i}+\frac{1}{2}\bar{m}_{nni}+b_{8}\bar{\Psi}_{nni}+b_{9}\bar{\Omega}_{i}\right),\\\
\bar{\Phi}_{nnni}&=-\varkappa\left(-\frac{4}{7}\bar{\mathscr{V}}_{i}-\frac{12}{35}\bar{q}_{i}+\frac{9}{7}\bar{m}_{nni}+b_{10}\bar{\Psi}_{nni}+b_{11}\bar{\Omega}_{i}\right),\end{aligned}\right\\}$
(45)

where

$\displaystyle b_{6}=-\frac{1}{36},\quad b_{7}=-\frac{1}{280},\quad
b_{8}=\frac{5}{12},\quad b_{9}=\frac{5}{56},\quad b_{10}=\frac{1}{14},\quad
b_{11}=-\frac{3}{490}.$ (46)

Comparison of the coefficients of
$\bar{\mathscr{V}}_{i},\bar{q}_{i},\bar{m}_{nni}$ in the PBC for the vector
fluxes (obtained on substituting $\bm{\mathcal{J}}$ and $\bm{\mathcal{F}}$
from (44) in (37)) and MBC (45) yields

$\displaystyle\bm{\mathcal{L}}=\varkappa\begin{bmatrix}[r]1&-1&-\frac{4}{7}\\\
-1&13&-\frac{4}{7}\\\ -\frac{4}{7}&-\frac{4}{7}&\frac{152}{49}\end{bmatrix}$
(47)

The matrix $\bm{\mathcal{L}}$ in (47) is symmetric and positive semidefinite
as $\varkappa$ is non-negative. Thus, the PBC for the vector fluxes are
obtained from (37) with $\bm{\mathcal{J}}$, $\bm{\mathcal{F}}$ and
$\bm{\mathcal{L}}$ being given by (44) and (47). On simplification, the PBC
for the vector fluxes are given by (45) with

$\displaystyle b_{6}=-\frac{1}{14},\quad b_{7}=-\frac{1}{70},\quad
b_{8}=\frac{13}{14},\quad b_{9}=\frac{13}{70},\quad b_{10}=-\frac{2}{49},\quad
b_{11}=-\frac{2}{245}.$ (48)

###  4.3 Boundary conditions for the tensor fluxes

####  4.3.1 Boundary conditions for the rank-2 tensor fluxes

The rank-2 tensor fluxes in (4) and their corresponding forces are collected
in $\bm{\mathcal{J}}$ and $\bm{\mathcal{F}}$ as

$\displaystyle\bm{\mathcal{J}}=\begin{bmatrix}\tilde{m}_{nij}\\\
\tilde{\Psi}_{nij}\end{bmatrix}\qquad\textrm{and}\qquad\bm{\mathcal{F}}=\begin{bmatrix}\frac{1}{2}\tilde{\sigma}_{ij}+\frac{1}{14}\tilde{R}_{ij}+\frac{1}{2}\tilde{\Phi}_{nnij}\\\
\frac{1}{28}\tilde{R}_{ij}\end{bmatrix}.$ (49)

Substituting $\bm{\mathcal{J}}$ and $\bm{\mathcal{F}}$ from (49) in (37), one
obtains the PBC for the rank-2 tensor fluxes. Nevertheless, the matrix
$\bm{\mathcal{L}}$ is still unknown. The elements of the matrix
$\bm{\mathcal{L}}$ can be obtained by comparing the coefficients of
$\tilde{\sigma}_{ij}$ and $\tilde{R}_{ij}$ in the corresponding terms of the
PBC and MBC for the rank-2 tensor fluxes.

The MBC for the rank-2 tensor fluxes (obtained through the Maxwell
accommodation model) read

$\displaystyle\left.\begin{aligned}
\tilde{m}_{nij}&=-\varkappa\left(\tilde{\sigma}_{ij}+\frac{1}{14}\tilde{R}_{ij}+b_{12}\tilde{\Phi}_{nnij}\right),\\\
\tilde{\Psi}_{nij}&=\varkappa\left(\tilde{\sigma}_{ij}-\frac{15}{14}\tilde{R}_{ij}+b_{13}\tilde{\Phi}_{nnij}\right),\end{aligned}\right\\}$
(50)

where

$\displaystyle b_{12}=\frac{1}{2}\quad\textrm{and}\quad b_{13}=-\frac{1}{2}.$
(51)

Comparison of the coefficients of $\tilde{\sigma}_{ij}$ and $\tilde{R}_{ij}$
in the PBC for the rank-2 tensor fluxes (obtained on substituting
$\bm{\mathcal{J}}$ and $\bm{\mathcal{F}}$ from (49) in (37)) and MBC (50)
yields

$\displaystyle\bm{\mathcal{L}}=\varkappa\begin{bmatrix}[r]2&-2\\\
-2&34\end{bmatrix}$ (52)

The matrix $\bm{\mathcal{L}}$ in (47) is symmetric and positive semidefinite
as $\varkappa$ is non-negative. Thus, the PBC for the rank-2 tensor fluxes are
obtained from (37) with $\bm{\mathcal{J}}$, $\bm{\mathcal{F}}$ and
$\bm{\mathcal{L}}$ being given by (49) and (52). On simplification, the PBC
for the rank-2 tensor fluxes are given by (50) with

$\displaystyle b_{12}=1\quad\textrm{and}\quad b_{13}=1.$ (53)

####  4.3.2 Boundary conditions for the rank-3 tensor fluxes

The rank-3 tensor flux in (4) and its corresponding forces are

$\displaystyle
J_{r}=\check{\Phi}_{nijk}\qquad\textrm{and}\qquad\mathcal{F}_{s}=\frac{1}{6}\check{m}_{ijk}.$
(54)

Substituting $\mathcal{J}_{r}$ and $\mathcal{F}_{s}$ from (54) in (37), one
obtains the PBC for the rank-3 tensor flux. Nevertheless, the coefficient
$\mathcal{L}_{rs}$ in (37) is still unknown but can be obtained by comparing
the coefficients of $\check{m}_{ijk}$ in the PBC and MBC for the rank-3 tensor
flux.

The MBC for rank-3 tensor flux (obtained through the Maxwell accommodation
model) read

$\displaystyle\check{\Phi}_{nijk}=-\varkappa\left(\check{m}_{ijk}+b_{14}\check{\Psi}_{ijk}\right)$
(55)

where

$\displaystyle b_{14}=\frac{1}{18}.$ (56)

Comparison of the coefficients of $\check{m}_{ijk}$ in the PBC for rank-3
tensor fluxes (obtained on substituting $\mathcal{J}_{r}$ and
$\mathcal{F}_{s}$ from (54) in (37)) and MBC (55) yields

$\displaystyle\mathcal{L}_{rs}=6\varkappa$ (57)

Thus, the PBC for the rank-3 tensor flux is obtained from (37) with
$\mathcal{J}_{r}$, $\mathcal{F}_{s}$ and $\mathcal{L}_{rs}$ being given by
(54) and (57). On simplification, the PBC for the rank-3 tensor flux is given
by (55) with

$\displaystyle b_{14}=0.$ (58)

##  5 Onsager reciprocity relations in some canonical problems

In this section, we consider two canonical boundary value problems, namely (i)
the Poiseuille and thermal transpiration gas flows and (ii) evaporation of a
spherical liquid droplet, for which the Onsager reciprocity relations have
been shown to hold for the linearised Boltzmann equation (LBE) (Wu et al.,
2017; Chernyak & Margilevskiy, 1989). We shall employ the LR26 equations with
the PBC and MBC to these problems to demonstrate that the solutions of the
LR26 equations with the MBC for these problems either do not satisfy the
Onsager reciprocity relations or do not respect the second law of
thermodynamics or do not comply with both. On the other hand, the solutions of
the LR26 equations with the PBC satisfy the Onsager reciprocity relations and
always respect the second law of thermodynamics. For comparison purpose, we
shall also include the results obtained through the LR13 equations with the
MBC and PBC.

###  5.1 Poiseuille flow and thermal transpiration flow

Let us first consider the steady state [$\partial_{t}(\cdot)=0$] of a gas
confined between two stationary ($v^{w}_{i}=0$), isothermal ($T^{w}=1$) and
fully diffusive ($\chi=1$) walls of a channel located at $y=\mp 1/2$; see
figure 1. The flow is assumed to be fully developed and driven either by a
uniform external force $F$ (Poiseuille flow) in the positive $x$-direction or
by a constant temperature gradient parallel to the walls $\tau$ (thermal
transpiration flow). The normals for the bottom and top walls are pointing in
$\pm y$-directions, respectively.

Figure 1: Schematic of the Poiseuille flow ($F=1$, $\tau=0$) and thermal
transpiration flow ($F=0$, $\tau=1$) in a channel.

The solution for the velocity and heat flux—for the problem under
consideration—from the LR26 equations read

$\displaystyle v_{x}$
$\displaystyle=C_{1}+\underline{\frac{F}{2\mathrm{Kn}}\left(\frac{1}{4}-y^{2}\right)}-\frac{2}{5}\sum_{i=1}^{4}\left(C_{i}^{K}-\frac{5}{2}\beta_{i}\right)\exp\left(\frac{y\gamma_{i}}{\mathrm{Kn}}\right),$
(59a) $\displaystyle q_{x}$
$\displaystyle=\underline{-\frac{5\mathrm{Kn}}{2\Pr}\tau}-\frac{\mathrm{Kn}}{\Pr}F+\sum_{i=1}^{4}C_{i}^{K}\exp\left(\frac{y\gamma_{i}}{\mathrm{Kn}}\right),$
(59b)

where $C_{1}$ and $C^{K}_{1}$, $C^{K}_{2}$, $C^{K}_{3}$ and $C^{K}_{4}$ are
the constants of integration. Moreover, due to symmetry about the $x$-axis,
$C^{K}_{2}=C^{K}_{1}$, $C^{K}_{4}=C^{K}_{3}$, $\beta_{2}=\beta_{1}$ and
$\beta_{4}=\beta_{3}$. The unspecified three constants of integration, i.e.
$C_{1}$, $C^{K}_{1}$ and $C^{K}_{3}$, are to be evaluated from the three
(vector) boundary conditions at $y=-1/2$.

Figure 2:  Poiseuille flow and thermal transpiration flow: (a) the total mass
flux $M^{P}$ in the Poiseuille flow and (b) violation of the Onsager
reciprocity relation, $M^{T}=Q^{P}$, for the MBC in the case of the LR26
equations (main panel) and in the case of the LR13 equations (inset). The
triangles in the left panel denote the results obtained from the LBE in Wu et
al. (2017).

The underlined terms in (59) represent the classical solution of the NSF
equations without velocity slip for the Poiseuille flow. The constant term
$C_{1}$ in (59) accounts for the velocity slip and terms with exponential
functions describe the Knudsen layer contributions to the velocity and heat
flux. In classical hydrodynamics, the velocity profile for the Poiseuille flow
is parabolic; the LR26 equations add velocity slip $C_{1}$ and Knudsen layer
contributions. Similarly, the Fourier’s contribution to the heat flux is given
by the underlined term in (59b). The coefficients $\gamma_{i}$’s
($i\in\\{1,2,3,4\\}$) in (59) determine the thickness of the Knudsen layer;
hence the LR26 equations describe the Knudsen layer as a superposition of four
exponential functions. The coefficients $\beta_{1}$, $\beta_{3}$ and
$\gamma_{i}$’s ($i\in\\{1,2,3,4\\}$) for MM read $\beta_{1}=-0.24697$,
$\beta_{3}=2.25494$, $\gamma_{1}=-0.510199$, $\gamma_{2}=0.510199$,
$\gamma_{3}=-1.26497$ and $\gamma_{4}=1.26497$ (Gu et al., 2010; Gu & Emerson,
2014).

The total mass flow and heat flow rates of the gas are, respectively, defined
as

$M=\int_{-1/2}^{1/2}v_{x}\,\mathrm{d}y\quad\textrm{and}\quad
Q=\int_{-1/2}^{1/2}q_{x}\,\mathrm{d}y.$ (60)

We shall use the superscripts ‘$P$’ and ‘$T$’ for the Poiseuille ($F=1$,
$\tau=0$) and transpiration ($F=0$, $\tau=1$) flows, respectively.

Figure 2(a) exhibits the mass flow rate, $M^{P}$, in the Poiseuille flow
plotted over the Knudsen number. The results computed from the MBC applied to
the LR26 equations with the MM collision model are compared with the LBE data
from Wu et al. (2017). The Knudsen minimum is observed in the Poiseuille flow,
in which the mass flow rate $M^{P}$ first decreases with the Knudsen number,
attains a minimum value around a critical Knudsen number and then increases;
see figure 2(a). The LR26 equations with PBC closely follow the mass flow rate
profile from the LBE up to $\mathrm{Kn}\lesssim 5$.

Due to the microscopic reversibility, the Onsager reciprocity relations should
hold, which in this case state that $M^{T}=Q^{P}$ (Wu et al., 2017). It is
evident from figure 2(b) that the MBC do not comply with the Onsager
reciprocity relations for the LR26 (and also for the LR13; shown in the inset)
equations, while the PBC do.

Figure 3:  The entropy generation rates at the boundary for (a) the Poiseuille
flow and (b) the thermal transpiration flow. The main panels and insets
exhibit the results computed with the LR26 and LR13 equations, respectively.

Figure 3 illustrates the entropy generation rates at the interface computed
with the LR26 equations (main panels) and the LR13 equations (insets) in the
case of (a) the Poiseuille flow and (b) the transpiration flow. It is clear
from figure 3(a) that the MBC for the LR26 equations (as well as for the LR13
equations) always violate the second law of thermodynamics (i.e. the entropy
production $\Sigma$ becomes negative) for large Knudsen numbers. For the
transpiration flow case (figure 3(b)), the violation of the second law of
thermodynamics by the MBC for the LR26 equations occurs already at a
relatively small Knudsen number ($\mathrm{Kn}<1$), which is within the
transition regime. Surprisingly, the LR13 equations with MBC in the
transpiration case do not contravene with the second law of thermodynamics
whilst they do in the Poiseuille flow case around $\mathrm{Kn}=1$; see the
insets in figure 3. Nevertheless, the PBC always predict a positive entropy
production in both the cases for both models in contrast to the conventional
MBC.

###  5.2 Evaporation of a spherical liquid droplet

We consider a liquid droplet of (fixed) initial radius $\hat{R}_{0}$ with a
given interface temperature $\hat{T}_{l}$ and the corresponding saturation
pressure $\hat{p}_{\mathrm{sat}}$, immersed in its own vapour; see figure 4.

Figure 4:  Schematic of an evaporating liquid droplet surrounded by its own
vapour.

The far-field (dimensional) temperature $\hat{T}_{\infty}$, pressure
$\hat{p}_{\infty}$ and the initial radius of the droplet $\hat{R}_{0}$ are
used for the purpose of non-dimensionalisation.

We consider two cases: (i) the flow is driven by a unit pressure difference
while the temperature of the liquid is equal to the far-field temperature
[i.e. $\Delta p:=(\hat{p}_{\mathrm{sat}}-\hat{p}_{\infty})/\hat{p}_{\infty}=1$
and $\Delta T:=(\hat{T}_{l}-\hat{T}_{\infty})/\hat{T}_{\infty}=0$] and (ii)
the flow is driven by a unit temperature difference while the saturation
pressure is the same as the far-field pressure (i.e. $\Delta p=0$ and $\Delta
T=1$). The problem has been investigated in detail by Rana et al. (2018b). We
skip the details here for the sake of succinctness and focus mainly on looking
at the Onsager reciprocity relations; for computational (and other) details,
the reader is referred to Rana et al. (2018b). The Onsager reciprocity
relations for this problem lead to $Q^{\Delta p}=M^{\Delta T}$ (Chernyak &
Margilevskiy, 1989), where $Q^{\Delta p}$ is the (dimensionless) heat flux in
the pressure-driven process, $M^{\Delta T}$ is the (dimensionless) mass flux
in the temperature-driven process, and the superscripts ‘$\Delta p$’ and
‘$\Delta T$’ are used to denote the pressure- and temperature-driven
processes, respectively.

The mass and heat fluxes in the pressure- and temperature-driven processes
computed from the R13 and R26 theories using the PBC and MBC are illustrated
in figure 5 for different values of the Knudsen number. The top and bottom
rows of figure 5 exhibit the results for the pressure- and temperature-driven
processes, respectively, while the left and right columns display the mass and
heat fluxes, respectively. For comparison purpose, the results obtained by
solving the LBE (Sone et al., 1994) are also included and they are depicted by
symbols in the figure. In general, the results from both the R13 and R26
theories agree with those from the LBE for $\mathrm{Kn}\lesssim 0.5$. However,
for the mass flux in the pressure-driven case [figure 5(a)], the R26 equations
on using the PBC (blue solid lines) give an excellent match with the LBE
results (symbols) for $\mathrm{Kn}\lesssim 0.5$ whereas the results from the
other theories are comparatively higher than the LBE results even for small
Knudsen numbers. Let us focus on panels (b) and (c) of figure 5, which
illustrate the cross-effects, i.e. the heat flux in the pressure-driven case
and mass flux in the temperature-driven case, respectively, and account for
the verification of the Onsager reciprocity relations. The LBE predicts the
same values the heat flux in the pressure-driven case and mass flux in the
temperature-driven case, thanks to the Onsager reciprocity relations (Chernyak
& Margilevskiy, 1989; Sone et al., 1994). Evidently from panels (b) and (c) of
figure 5, the heat flux in the pressure-driven case computed with the
LR26/LR13 theory using the PBC (blue solid line for the LR26 equations and
green solid line for the LR13 equations) seems to be the same as the mass flux
in the temperature-driven case computed with the LR26/LR13 theory using the
PBC (blue solid line for the LR26 equations and green solid line for the LR13
equations) whereas the heat flux in the pressure-driven case computed with the
LR26/LR13 theory using the MBC (red dashed line for the LR26 equations and
magenta dashed line for the LR13 equations) seems to be different from the
mass flux in the temperature-driven case computed with the LR26/LR13 theory
using the MBC (red dashed line for the LR26 equations and magenta dashed line
for the LR13 equations). In order to have a clear picture of the Onsager
reciprocity relations, the difference of $Q^{\Delta p}$ and $M^{\Delta T}$ is
plotted in figure 6, which clearly shows that both the LR13 and LR26 equations
on using the MBC for this problem violate the Onsager reciprocity relations
while both the LR13 and LR26 equations on using the PBC respect the Onsager
reciprocity relations for all Knudsen numbers.

Figure 5: Profiles of the mass and heat fluxes plotted over the Knudsen number
for the problem of evaporation of a spherical liquid droplet: (a) mass flux in
the pressure-driven case, (b) heat flux in the pressure-driven case, (c) mass
flux in the temperature-driven case and (d) heat flux in the temperature-
driven case. The LBE data (symbols) are taken from Sone et al. (1994). Figure
6:  Violation of the Onsager reciprocity relations for an evaporating droplet.

We have also computed the entropy generation rate at the interface for this
problem and found (but not shown here for the sake of conciseness) that it
remains positive for the R13 and R26 equations with the PBC in both pressure-
and temperature-driven cases, also remains positive for the R13 equations with
the MBC in the pressure-driven case but does not remain positive for the R13
equations with the MBC in the temperature-driven case for all Knudsen numbers.
The R26 equations with the MBC yield a positive entropy generation rate for
this problem.

##  6 Gas flow past an evaporating droplet: a new analytic solution

As an application of the above findings, we compute the analytic solution of
the LR26 equations for a steady vapour passing over an evaporating spherical
droplet. Furthermore, we also consider a special case of this problem—when
there is no evaporation; this is essentially the problem of a streaming flow
of a gas past a rigid sphere, which is a classical problem in fluid dynamics.
The results presented in this section are determined using the PBC only, since
we have established that they are the physically admissible boundary
conditions. Nevertheless, we show in the supplementary material that the
difference between the results obtained using the PBC and MBC for this problem
is negligibly small.

Figure 7:  Schematic of a flow past a sphere.

We consider a rarefied gas flow past a spherical droplet as illustrated in
figure 7. It is convenient to study the problem in the spherical coordinate
system $(r,\theta,\phi)$ that relates to the Cartesian coordinate system
$(x,y,z)$ via
$(x,y,z)\equiv(r\sin{\theta}\cos{\phi},r\sin{\theta}\sin{\phi},r\cos{\theta})$,
where $r\in[0,\infty)$, $\theta\in[0,\pi]$ and $\phi\in[0,2\pi)$. The flow is
assumed to be in the $(r,\theta)$-plane coming from the left as shown in the
figure. Owing to spherical symmetry, the flow is independent of the direction
$\phi$. Consequently, all the field variables can depend only on $r$ and
$\theta$. All the vectors and tensors involved in the problem need to be
converted to spherical coordinates. To this end, it is worthwhile noting that
a symmetric-tracefree $k$-rank tensor has $2k+1$ independent components in
three dimensions (Gupta, 2020).

For the considered problem, any vector $a_{i}\in\\{v_{i},q_{i},\Omega_{i}\\}$
in the spherical coordinate system is given by

$\displaystyle a_{i}=(a_{r},a_{\theta},0)^{\mathsf{T}}.$ (61)

A symmetric-tracefree second-rank tensor $a_{ij}\in\\{\sigma_{ij},R_{ij}\\}$
in the spherical coordinate system is given by

$\displaystyle a_{ij}=\begin{bmatrix}a_{rr}&a_{r\theta}&0\\\
a_{r\theta}&-\frac{1}{2}a_{rr}&0\\\ 0&0&-\frac{1}{2}a_{rr}\end{bmatrix}.$ (62)

Note that a symmetric-tracefree second-rank tensor has five independent
components; however, for the problem under consideration, it is clear from
(62) that all components of $a_{ij}$ can be expressed in terms of only two
components, namely $a_{rr}$ and $a_{r\theta}$. A symmetric-tracefree third-
rank tensor has seven independent components; however, for the problem under
consideration, all components of $a_{ijk}$ can again be expressed in terms of
$a_{rrr}$ and $a_{rr\theta}$. The five independent components (other than
$a_{rrr}$ and $a_{rr\theta}$) of $a_{ijk}$, for the problem under
consideration, read

$\displaystyle a_{r\theta\theta}=-\frac{1}{2}a_{rrr},\quad
a_{\theta\theta\theta}=-\frac{3}{4}a_{rr\theta},\quad
a_{rr\phi}=a_{\theta\theta\phi}=a_{r\theta\phi}=0.$ (63)

A symmetric-tracefree fourth-rank tensor $a_{ijkl}$ has nine independent
components; however, for the problem under consideration, all components of
$a_{ijkl}$ can again be expressed in terms of $a_{rrrr}$ and $a_{rrr\theta}$.
The seven independent components (other than $a_{rrrr}$ and $a_{rrr\theta}$)
of $a_{ijkl}$, for the problem under consideration, read

$\displaystyle\left.\begin{gathered}a_{rr\theta\theta}=-\frac{1}{2}a_{rrrr},\quad
a_{\theta\theta\theta\theta}=\frac{3}{8}a_{rrrr},\quad
a_{r\theta\theta\theta}=-\frac{3}{4}a_{rrr\theta},\\\
a_{rrr\phi}=a_{rr\theta\phi}=a_{r\theta\theta\phi}=a_{\theta\theta\theta\phi}=0.\end{gathered}\right\\}$
(66)

Accordingly, the conservation laws (7)–(9) for this problem reads

$\displaystyle\frac{\partial v_{r}}{\partial
r}+\frac{2v_{r}}{r}+\frac{\mathscr{D}v_{\theta}}{r}$ $\displaystyle=0,$ (67a)
$\displaystyle\frac{\partial p}{\partial
r}+\frac{\partial\sigma_{rr}}{\partial
r}+\frac{3\sigma_{rr}}{r}+\frac{\mathscr{D}\sigma_{r\theta}}{r}$
$\displaystyle=0,$ (67b)
$\displaystyle\frac{\partial\sigma_{r\theta}}{\partial
r}+\frac{3\sigma_{r\theta}}{r}-\frac{1}{2r}\frac{\partial\sigma_{rr}}{\partial\theta}+\frac{1}{r}\frac{\partial
p}{\partial\theta}$ $\displaystyle=0,$ (67c) $\displaystyle\frac{\partial
q_{r}}{\partial r}+\frac{2q_{r}}{r}+\frac{\mathscr{D}q_{\theta}}{r}$
$\displaystyle=0,$ (67d)

where $\mathscr{D}\equiv\cot{\theta}+\frac{\partial}{\partial\theta}$. The
stress balance equation (10) reduces to

$\displaystyle\frac{\partial m_{rrr}}{\partial r}$
$\displaystyle+\frac{4m_{rrr}}{r}+\frac{4}{5}\frac{\partial q_{r}}{\partial
r}+2\frac{\partial v_{r}}{\partial
r}+\frac{\mathscr{D}m_{rr\theta}}{r}=-\frac{1}{\mathrm{Kn}}\sigma_{rr},$ (68a)
$\displaystyle\frac{\partial m_{rr\theta}}{\partial r}$
$\displaystyle+\frac{4m_{rr\theta}}{r}+\frac{2}{5}\left(\frac{\partial
q_{\theta}}{\partial r}-\frac{q_{\theta}}{r}\right)+\frac{\partial
v_{\theta}}{\partial r}-\frac{v_{\theta}}{r}$
$\displaystyle-\frac{1}{2r}\frac{\partial
m_{rrr}}{\partial\theta}+\frac{1}{r}\frac{\partial
v_{r}}{\partial\theta}+\frac{2}{5r}\frac{\partial
q_{r}}{\partial\theta}=-\frac{1}{\mathrm{Kn}}\sigma_{r\theta}.$ (68b)

The heat flux balance equation (11) gives

$\displaystyle\frac{1}{2}\left(\frac{\partial R_{rr}}{\partial
r}+\frac{\mathscr{D}R_{r\theta}}{r}+\frac{3R_{rr}}{r}\right)+\frac{1}{6}\frac{\partial\Delta}{\partial
r}-\frac{\partial p}{\partial r}+\frac{5}{2}\frac{\partial T}{\partial r}$
$\displaystyle=-\frac{\mathrm{Pr}}{\mathrm{Kn}}q_{r},$ (69a)
$\displaystyle\frac{1}{2}\left(\frac{\partial R_{r\theta}}{\partial
r}+\frac{3R_{r\theta}}{r}\right)+\frac{1}{6r}\frac{\partial\Delta}{\partial\theta}-\frac{1}{4r}\frac{\partial
R_{rr}}{\partial\theta}-\frac{1}{r}\frac{\partial
p}{\partial\theta}+\frac{5}{2r}\frac{\partial T}{\partial\theta}$
$\displaystyle=-\frac{\mathrm{Pr}}{\mathrm{Kn}}q_{\theta}.$ (69b)

The $m_{ijk}$ balance equation (15) simplifies to

$\displaystyle\frac{1}{r}\mathscr{D}\bigg{[}-\frac{6}{5}\sigma_{r\theta}$
$\displaystyle+\Phi_{rrr\theta}-\frac{6}{35}R_{r\theta}\bigg{]}+\frac{9}{5}\left(\frac{\partial}{\partial
r}-\frac{2}{r}\right)\sigma_{rr}$
$\displaystyle+\left(\frac{\partial}{\partial
r}+\frac{5}{r}\right)\Phi_{rrrr}+\frac{9}{35}\left(\frac{\partial}{\partial
r}-\frac{2}{r}\right)R_{rr}=-\frac{\mathrm{Pr}_{m}}{\mathrm{Kn}}m_{rrr},$
(70a) $\displaystyle\frac{6}{5r}\frac{\partial\sigma_{rr}}{\partial\theta}$
$\displaystyle+\frac{6}{35r}\frac{\partial
R_{rr}}{\partial\theta}-\frac{1}{2r}\frac{\partial\Phi_{rrrr}}{\partial\theta}+\frac{8}{5}\left(\frac{\partial}{\partial
r}-\frac{2}{r}\right)\sigma_{r\theta}$
$\displaystyle+\frac{8}{35}\left(\frac{\partial}{\partial
r}-\frac{2}{r}\right)R_{r\theta}+\left(\frac{\partial}{\partial
r}+\frac{5}{r}\right)\Phi_{rrr\theta}=-\frac{\mathrm{Pr}_{m}}{\mathrm{Kn}}m_{rr\theta}.$
(70b)

The $R_{ij}$ balance equation (16) reduces to

$\displaystyle\frac{1}{r}\mathscr{D}\bigg{[}2m_{rr\theta}$
$\displaystyle-\frac{2}{15}\Omega_{\theta}+\psi_{rr\theta}-\frac{28}{15}q_{\theta}\bigg{]}+\frac{56}{15}\left(\frac{\partial}{\partial
r}-\frac{1}{r}\right)q_{r}+2\left(\frac{\partial}{\partial
r}+\frac{4}{r}\right)m_{rrr}$ $\displaystyle+\left(\frac{\partial}{\partial
r}+\frac{4}{r}\right)\psi_{rrr}+\frac{4}{15}\left(\frac{\partial}{\partial
r}-\frac{1}{r}\right)\Omega_{r}=-\frac{\mathrm{Pr}_{R}}{\mathrm{Kn}}R_{rr},$
(71a) $\displaystyle 2\left(\frac{\partial}{\partial
r}+\frac{4}{r}\right)m_{rr\theta}$
$\displaystyle+\left(\frac{\partial}{\partial
r}+\frac{4}{r}\right)\psi_{rr\theta}+\frac{1}{5}\left(\frac{\partial}{\partial
r}-\frac{1}{r}\right)\Omega_{\theta}+\frac{14}{5}\left(\frac{\partial}{\partial
r}-\frac{1}{r}\right)q_{\theta}$ $\displaystyle-\frac{1}{r}\frac{\partial
m_{rrr}}{\partial\theta}+\frac{14}{5r}\frac{\partial
q_{r}}{\partial\theta}-\frac{1}{2r}\frac{\partial\psi_{rrr}}{\partial\theta}+\frac{2}{5r}\frac{\partial\Omega_{r}}{\partial\theta}=-\frac{\mathrm{Pr}_{R}}{\mathrm{Kn}}R_{r\theta}.$
(71b)

The $\Delta$ balance equation (17) gives

$\frac{1}{r}\mathscr{D}\left(8q_{\theta}+\Omega_{\theta}\right)+\left(\frac{\partial}{\partial
r}+\frac{2}{r}\right)\left(8q_{r}+\Omega_{r}\right)=-\frac{\mathrm{Pr}_{\Delta}}{\mathrm{Kn}}\Delta.$
(72)

Finally, the LR26 closure (22) simplifies to

$\displaystyle\frac{4}{7}\left(\frac{\partial}{\partial
r}-\frac{3}{r}\right)m_{rrr}-\frac{3}{7}\frac{\mathscr{D}m_{rr\theta}}{r}$
$\displaystyle=-\frac{1}{4}\frac{\mathrm{Pr}_{\Phi}}{\mathrm{Kn}}\Phi_{rrrr},$
(73a) $\displaystyle\frac{15}{28}\left(\frac{\partial}{\partial
r}-\frac{3}{r}\right)m_{rr\theta}+\frac{5}{14r}\frac{\partial
m_{rrr}}{\partial\theta}$
$\displaystyle=-\frac{1}{4}\frac{\mathrm{Pr}_{\Phi}}{\mathrm{Kn}}\Phi_{rrr\theta},$
(73b) $\displaystyle\frac{3}{5}\left(\frac{\partial}{\partial
r}-\frac{2}{r}\right)R_{rr}-\frac{2}{5}\frac{\mathscr{D}R_{r\theta}}{r}$
$\displaystyle=-\frac{7}{27}\frac{\mathrm{Pr}_{\psi}}{\mathrm{Kn}}\psi_{rrr},$
(73c) $\displaystyle\frac{8}{15}\left(\frac{\partial}{\partial
r}-\frac{2}{r}\right)\text{$R_{r\theta}$}+\frac{2}{5r}\frac{\partial
R_{rr}}{\partial\theta}$
$\displaystyle=-\frac{7}{27}\frac{\mathrm{Pr}_{\psi}}{\mathrm{Kn}}\psi_{rr\theta},$
(73d)

and

$\displaystyle\frac{\partial\Delta}{\partial
r}+\frac{12}{7r}\mathscr{D}R_{r\theta}+\frac{12}{7}\left(\frac{\partial}{\partial
r}+\frac{3}{r}\right)R_{rr}$
$\displaystyle=-\frac{3}{7}\frac{\mathrm{Pr}_{\Omega}}{\mathrm{Kn}}\Omega_{r},$
(74a)
$\displaystyle\frac{1}{r}\frac{\partial\Delta}{\partial\theta}+\frac{12}{7}\left(\frac{\partial
R_{r\theta}}{\partial
r}+\frac{3R_{r\theta}}{r}\right)-\frac{6}{7r}\frac{\partial
R_{rr}}{\partial\theta}$
$\displaystyle=-\frac{3}{7}\frac{\mathrm{Pr}_{\Omega}}{\mathrm{Kn}}\Omega_{\theta}.$
(74b)

The mathematical structure of the LR26 equation is the same as that of the
LR13 equations, with which Torrilhon (2010) obtained an analytic solution for
the problem of gas flow past a rigid sphere. Following Torrilhon (2010), we
shall adopt the same solution method to solve (67)–(74). In the next
subsection, we shall outline the steps involved in solving the LR26 equations
for the considered problem. For more details, the reader is referred to
Torrilhon (2010).

###  6.1 Analytic solution of the LR26 equations

The analytic solutions of the LR26 equations (67)–(74) for the aforesaid
problem can be obtained by converting these partial differential equations to
ordinary differential equations via assuming an explicit dependency on the
azimuthal angle $\theta$. To this end, the vectorial and tensorial components
with an odd number of indices in $\theta$ are taken to be proportional to
$\sin{\theta}$ while scalars and tensorial components with an even number of
indices in $\theta$ are taken to be proportional to $\cos{\theta}$ (Torrilhon,
2010). The components for vector fields, namely the velocity
$\bm{v}(r,\theta)$ and the heat flux $\bm{q}(r,\theta)$, are taken to be

$\displaystyle\left.\begin{aligned}
v_{r}(r,\theta)&=u_{\infty}[1+a_{1}(r)]\cos{\theta},&\quad
v_{\theta}(r,\theta)&=-u_{\infty}[1+a_{2}(r)]\sin{\theta},\\\
q_{r}(r,\theta)&=\alpha_{1}(r)\cos{\theta},&\quad
q_{\theta}(r,\theta)&=\alpha_{2}(r)\sin{\theta}.\end{aligned}\right\\}$ (75)

The scalars, namely $T(r,\theta)$, $p(r,\theta)$ and $\Delta(r,\theta)$, are
taken to be

$\displaystyle T(r,\theta)=c(r)\cos{\theta},\quad
p(r,\theta)=d(r)\cos{\theta},\quad\Delta(r,\theta)=\delta(r)\cos{\theta}.$
(76)

The components of the second- and third-rank tensors are taken to be

$\displaystyle\left.\begin{aligned}
\sigma_{rr}(r,\theta)&=b_{1}(r)\cos{\theta},&\quad\sigma_{r\theta}(r,\theta)&=b_{2}(r)\sin{\theta},\\\
R_{rr}(r,\theta)&=\beta_{1}(r)\cos{\theta},&\quad
R_{r\theta}(r,\theta)&=\beta_{2}(r)\sin{\theta},\\\
m_{rrr}(r,\theta)&=\kappa_{1}(r)\cos{\theta},&\quad
m_{rr\theta}(r,\theta)&=\kappa_{2}(r)\sin{\theta}.\end{aligned}\right\\}$ (77)

In (75)–(77), $a_{1}(r)$, $a_{2}(r)$, $\alpha_{1}(r)$, $\alpha_{2}(r)$,
$c(r)$, $d(r)$, $\delta(r)$, $b_{1}(r)$, $b_{2}(r)$, $\beta_{1}(r)$,
$\beta_{2}(r)$, $\kappa_{1}(r)$ and $\kappa_{2}(r)$ are thirteen unknowns in
total, which need to be determined from thirteen equations (67)–(72).

Substitution of constitutive relations (73)–(74) along with solution ansatz
(75)–(77) in (67)–(72) results into a set of eight first-order ordinary
differential equations (emanating from (67)–(69)) and five second-order
ordinary differential equations (from (70)–(72)). The latter include second
derivatives of $\delta(r)$, $\beta_{1}(r)$, $\beta_{2}(r)$, $\kappa_{1}(r)$
and $\kappa_{2}(r)$ due to regularised constitutive relations (73)–(74). With
this, one obtains an analytic solution for the unknowns $a_{1}(r)$,
$a_{2}(r)$, $\alpha_{1}(r)$, $\alpha_{2}(r)$, $c(r)$, $d(r)$, $\delta(r)$,
$b_{1}(r)$, $b_{2}(r)$, $\beta_{1}(r)$, $\beta_{2}(r)$, $\kappa_{1}(r)$ and
$\kappa_{2}(r)$, which has been relegated to appendix E for better
readability. Expressions (E)–(E) contain eight integration constants
$C_{1,2,3}$ and $K_{1,2,3,4,5}$ that are determined from boundary conditions
(39)–(45).

Let us investigate a few qualitative features of these solutions. These
solutions are the sums of the bulk contributions (which are polynomials in
$1/r$) and the Knudsen layer contributions (terms containing the exponential
functions $\Theta_{i}(r)$); see appendix E. The bulk solutions contain three
constants of integration $C_{1,2,3}$ while Knudsen layers introduce five
additional constants—$K_{1,2}$ for vector variables and $K_{3,4,5}$ for scalar
variables. The bulk part can further be divided into the classical NSF
contribution and second-order contribution (in the Knudsen number) that is
neglected in the classical NSF theory. In the hydrodynamic limit (i.e. Stokes
flow without velocity slip and temperature jump), $C_{1}=-3$ and $C_{2}=1/2$
while $C_{3}=0$ (Lamb, 1932). The last relation ($C_{3}=0$) implies an
isothermal flow condition. Hence, the non-zero bulk contributions to the
temperature and heat flux stem from the second-order terms, i.e. the terms
with $C_{3}$ in equations (Ea) and (Ec,d). The NSF equations with velocity-
slip boundary condition predict $C_{1}$ and $C_{2}$ as a function of the
Knudsen number $\mathrm{Kn}$ and the accommodation coefficient $\chi$
(Torrilhon, 2010).

For the Knudsen layer contributions in (E)–(E), the LR13 theory (Torrilhon,
2010) approximates Knudsen layers with two exponential functions whilst the
LR26 theory give superposition of five exponential functions to describe the
Knudsen layers. On the contrary, the NSF theory does not predict Knudsen
layers at all. In the following, we shall show that the LR26 theory gives a
remarkably improved flow description with added details to the Knudsen layers.

###  6.2 Gas flow past a rigid sphere

First, let us consider the classical problem of a gas flow past a rigid
sphere. The same problem was investigated by Torrilhon (2010) using the LR13
and linear NSF equations with velocity slip.

Figure 8 illustrates the profiles of the radial and polar components of the
velocity, pressure and temperature as functions of the radial distance $r$
starting from the surface of the sphere $r=1$.

Figure 8: Profiles of the (a) radial velocity, (b) polar velocity (c) pressure
and (d) temperature of the gas past a rigid sphere as functions of the radial
distance from the surface of the sphere for various Knudsen numbers. Results
obtained from the LR26 theory (solid lines) are compared with those obtained
from the LR13 (dashed lines) and NSF (dotted lines) theories (Torrilhon, 2010)
and with the numerical solution of the LBE (symbols) from Takata et al.
(1993).

The exact solutions for the problem under consideration obtained with the LR26
equations in the present work are compared with those obtained with the LR13
and linear NSF solutions in Torrilhon (2010) and with the numerical solutions
of the LBE for three values of the Knudsen number $\mathsf{Kn}=0.1,0.4,1$
given in Takata et al. (1993). Note that the Knudsen number $\mathsf{Kn}$ used
in Takata et al. (1993) relates to the Knudsen number $\mathrm{Kn}$ used in
the present work via $\mathrm{Kn}=1.27\,\mathsf{Kn}/\sqrt{2}$. Therefore, the
Knudsen number values $\mathsf{Kn}=0.1,0.4,1$ in Takata et al. (1993)
correspond to the Knudsen number values $\mathrm{Kn}\approx 0.09,0.36,0.9$,
respectively, used in the present work. As expected, for small values of the
Knudsen number [red curves and symbols ($\mathrm{Kn}=0.09$)], all theories,
including the NSF theory, give a good match with the LBE data for the
velocities. However, for moderate values of the Knudsen number [blue curves
and symbols ($\mathrm{Kn}=0.36$)], the results from the NSF theory show
deviations from those from kinetic theory (LBE results) while the LR13 theory
still gives reasonably good results. The mismatch between the solutions from
the NSF equations and LBE is attributed to the lack of Knudsen layer
description by the former. For even larger values of the Knudsen number [green
curves and symbols ($\mathrm{Kn}=0.9$)], even the results from the LR13 theory
start deviating from the LBE results. Notably, the LR13 theory in this case
predicts $u_{\theta}/(u_{\infty}\sin{\theta})>1$ near the surface of the
sphere, which means that the slip velocity on the top of the sphere
($\theta=\pi/2$) is larger than the far field velocity—a result which is
clearly unphysical and contradicted by more accurate theories (LBE and LR26).
Hence, one may conclude that for large values of the Knudsen number, the
results from LR13 theory cannot be trusted. On the other hand, the LR26 theory
offers an excellent match upto $\mathrm{Kn}\sim 1$. Interestingly, the LR26
theory also predicts $u_{\theta}/(u_{\infty}\sin{\theta})>1$ for
$\mathrm{Kn}\gtrsim 5$, which deduces that the LR26 theory can predict
reasonably good results for the Knudsen numbers $0<\mathrm{Kn}\lesssim 5$.

Another remarkable feature in figure 8(d) is that the NSF theory predicts a
non-negative temperature profile whereas all other theories, including the
LBE, give a non-positive temperature profile.

Figure 9:  Heat flux lines drawn on top of the temperature contours for
$\mathrm{Kn}=0.09$ (top row) and $\mathrm{Kn}=0.36$ (bottom row) from the
LR26, LR13 and NSF theories.

To demonstrate the consequences of this, the heat flux lines and temperature
contours obtained from LR26, LR13 and NSF theories are illustrated in figure
9. It is evident from the figure that all the theories, including the NSF
theory, do predict the thermal polarisation—an effect that, for the problem
under consideration, shows different temperatures in the front and back sides
of the sphere in the absence of nonlinear dissipations and external
temperature differences (Aoki & Sone, 1987; Takata et al., 1993; Torrilhon,
2010). Nevertheless, the LR13 and LR26 theories show a hot region in the front
side of the sphere and a cold region in the back side of the sphere while the
temperature prediction by the NSF theory is the other way round. The
temperature predictions by the LR13 and LR26 theories are in agreement with
the findings of Takata et al. (1993) based on the LBE and of Torrilhon (2010)
based on the LR13 equations, and also make sense because the temperature is
likely to be more on the front side of the sphere due to compression of the
gas at the front side of the sphere and is likely to be comparatively lesser
on the back side of the sphere due to expansion of the gas at the back side of
the sphere. Furthermore, the LR13 and LR26 theories predict that the heat
flows from the cold region (back side of the sphere) to the hot region (front
side of the sphere), which is a non-Fourier effect. The NSF theory cannot
capture this effect and shows the heat flowing from the hot region to the cold
region. Such a non-Fourier effect is well established in lid-driven cavity
problems, (see e.g. Rana et al., 2013, and references therein). By comparing
the top and bottom rows of figure 9, it is also clear that the non-uniformity
in temperature (i.e., the temperature difference between the temperatures at
the front and back sides of the sphere) is more for a larger Knudsen number.
This is also in agreement with the findings of Takata et al. (1993).

Figure 10:  Normalised drag force from the NSF (dot-dashed red line), R13
(dashed green line) and R26 (solid blue line) theories as a function of the
Knudsen number. The circles denote the experimental data of Goldberg (1954)
(fitted by Bailey et al. (2005)). The results are normalised with the Stokes
drag $F_{\mathrm{Stokes}}=6\pi\mathrm{Kn}u_{\infty}$.

Figure 10 illustrates the drag force normalised with the Stokes drag,
$F_{\mathrm{Stokes}}=6\pi\mathrm{Kn}u_{\infty}$, plotted over the Knudsen
number. The experimental data from Goldberg (1954) (fitted by Bailey et al.
(2005)) is included for comparison. As expected, all theories agree with the
experimental data for small values of the Knudsen number. As the Knudsen
number increases, the drag decreases due to slip, which in not embedded in the
Stokes formula and hence the Stokes formula gives vanishing drag force for all
values of the Knudsen number. Notably, the drag force from the NSF theory with
velocity-slip boundary conditions attains a nonzero positive value as
$\mathrm{Kn}\to\infty$ while the other theories (LR13 and LR26) predict a
vanishing drag in this limit—an observation tantamount to experimental
results. Furthermore, the LR26 theory provides an excellent prediction for the
drag force for values of the Knudsen number beyond unity (upto
$\mathrm{Kn}\approx 10$). On the contrary, the LR13 theory underpredicts the
drag force significantly for $\mathrm{Kn}\gtrsim 0.2$.

###  6.3 Gas flow past an evaporating spherical droplet

Finally, we consider the case of an evaporating droplet placed in a saturated
vapour ($p_{\mathrm{sat}}=p_{\infty}=0$) flowing past the droplet with
velocity $u_{\infty}$, i.e. phase change is driven by the flow field, not by
the far field. The assumptions made in the previous section remain true here,
i.e. the shape and size of the droplet remain fixed and the internal motion of
the droplet is ignored. In addition, the evaporation coefficient $\vartheta$
is assumed to be unity.

Figures 11 again illustrates the profiles of the velocity components, pressure
and temperature as functions of the radial direction $r$ starting from the
surface of the sphere $r=1$. The exact solutions for the problem under
consideration obtained with the LR26 equations (solid lines) and LR13
equations (dashed lines) are compared with those obtained with the numerical
solutions of the LBE for three values of the Knudsen number
$\mathsf{Kn}=0.1,0.4,1$ given in Sone et al. (1994). Notice that the Knudsen
number $\mathsf{Kn}$ is the same as that used in Takata et al. (1993).
Therefore, similarly to the above, these Knudsen number values correspond to
$\mathrm{Kn}=0.09,0.26,0.9$, respectively.

Figure 11: Profiles of the (a) radial velocity, (b) polar velocity (c)
pressure and (d) temperature of the gas past an evaporating spherical droplet
as functions of the radial distance from the surface of the droplet for
various Knudsen numbers. Results obtained from the LR26 theory (solid lines)
are compared with the LR13 theory (dashed lines) and with those obtained from
the numerical solution of the LBE (symbols) from Sone et al. (1994).

As expected, when the Knudsen number is small [red lines and symbols
($\mathrm{Kn}=0.09$)], both the LR13 and LR26 theories (dashed and solid
lines) agree well with the results from the LBE (symbols). However, as the
Knudsen number increases, the differences between the results from the two
theories become noticeable—in particular, for the radial velocity and
temperature. At higher Knudsen numbers ($\mathrm{Kn}\gtrsim 0.5$), the Knudsen
layer description by the LR13 theory is inaccurate whilst the LR26 theory
remains valid even for $\mathrm{Kn}=1$.

Figure 12:  Normalised drag force from the LR13 (dashed green line) and LR26
(solid blue line) theories and from the numerical solution of the LBE
(circles) from Sone et al. (1994) as a function of the Knudsen number. The
results are normalised with the Stokes drag
$F_{\mathrm{Stokes}}=6\pi\mathrm{Kn}u_{\infty}$.

The normalised drag force in the case of a gas flow past an evaporating
droplet is illustrated in figure 12. The solutions obtained from the NSF, LR13
and LR26 theories are compared against the numerical solutions of the LBE
given in Sone et al. (1994). In contrast to the non-evaporating droplet case,
all theories in this case predict vanishing drag as $\mathrm{Kn}\to\infty$.
For large values of the Knudsen number, the NSF theory overpredicts the drag
while the LR13 theory underpredicts it. Again, the LR26 theory gives accurate
predictions for the drag force even for large values of the Knudsen number
($\mathrm{Kn}\lesssim 10$).

##  7 Conclusion

In this article, we have addressed the question of thermodynamically
admissible boundary conditions for the LR26 equations. It has been shown that
these equations are accompanied with a proper entropy inequality with non-
negative entropy generation (the $H$-theorem).

The entropy production rate at the interface has been written as a sum of the
products of the fluxes with odd degree in $n$ and the requirement of the non-
negative entropy generation at the wall has been imposed, which have
eventually yielded the desired boundary conditions for the LR26 equations.
Furthermore, considering the coupling between the cross-terms of the same
tensorial structure (the Curie principle) a positive entropy production has
been achieved by setting the unknown boundary values for the odd fluxes
proportional to their driving forces. The phenomenological coefficients
appearing in these boundary conditions have been obtained by comparing them
with those obtained via the Maxwell accommodation model (with Grad’s
approximation) to the leading order, thus allowing us to formulate the PBC
aligned with the Maxwell kinetic-wall model. We have demonstrated that the
solutions of the LR26 equations with the MBC either do not satisfy the Onsager
reciprocity relations or do not respect the second law of thermodynamics or do
not comply with both. On the other hand, the solutions of the LR26 equations
with the PBC satisfy the Onsager reciprocity relations and always respect the
second law of thermodynamics.

As an application of the PBC, we have computed the analytic solution of the
LR26 equations for a steady vapour passing over an evaporating spherical
droplet. To the best of authors’ knowledge, this is the first time this
problem has been considered within the extended fluid dynamics framework.
These analytic solution have been found to be useful in understanding some
intriguing rarefaction effects, such as the Knudsen layers and thermal
polarisation. A special case of this problem, when there is no evaporation,
has also been considered, which is a classic problem in fluid dynamics. It has
been found that the LR26 equations with the PBC are able to provide a more
detailed and accurate description of the flow features emanating from
rarefaction effects up to $\mathrm{Kn}\lesssim 1$. Furthermore, the analytic
results for the drag force obtained from the LR26 theory with the PBC have
been found to be in excellent agreement with the theory/experiments even for
very large values of the Knudsen number—in the case of both an evaporating
droplet and a non-evaporating droplet.

In § 4, it was noticed that matching the coefficients of the lower-order
moments (i.e. moments involved in the G26 equations) in the PBC and MBC yields
the matrices $\mathcal{L}$ to be symmetric and semi-positive definite. The
reason behind it is presently unclear and will be explored in the future work.
Furthermore, it is worthwhile noting that the results in the present work are
valid only for slow moving droplets, i.e. when the Reynolds number
$\mathrm{Re}$ ($=\mathrm{Ma}\sqrt{\gamma}/\mathrm{Kn})$ is much less than
unity. Here $\gamma$ (=5/3 for monatomic gases) is the ratio of specific heat
capacities and $\mathrm{Ma}$ is the Mach number. For a relatively large value
of the Reynolds number, the nonlinear convective terms play an important role
in evaporation dynamics, which have been neglected in the present work. In
classical fluid dynamics texts (see e.g. Lamb, 1932), Oseen correction terms
are utilised to obtain solutions for the flow past various bodies at moderate
Reynolds numbers [$\mathcal{O}(\mathrm{Re})\sim 1$]. Oseen correction terms
are linear in nature as oppose to the non-linear Navier–Stokes operator,
making these equations analytically tractable. The present work can be
extended to flows with a small but finite Reynolds number by including the
Oseen correction terms (with the partial inclusion of convective terms) in the
future to study the combined effects of inertia and rarefaction on evaporation
dynamics. Motivated by the present findings, the present work will also be
extended to investigate unsteady flow problems and to incorporate coupling
with liquid dynamics (see e.g. Rana et al., 2019) elsewhere in the future.
Thermophoresis on a sphere is a classic problem in rarefied gas flows, which
has been studied by many researchers; see e.g. Yamamoto & Ishihara (1988);
Takata et al. (1994); Takata & Sone (1995); Beresnev & Chernyak (1995); Takata
(2009); Padrino et al. (2019); Kalempa & Sharipov (2020). The model developed
in the present work can also be used to investigate thermophoresis on an
evaporating droplet in the future.

##### Acknowledgments.

ASR acknowledges the funding from the BITS, Pilani (India) through the
Research Initiation Grant. VKG gratefully acknowledges financial support
through the ‘MATRICS’ project MTR/2017/000693 funded by the SERB, India. This
work has also been financially supported in the UK by EPSRC Grants No.
EP/N016602/1, EP/P031684/1, EP/S029966/1.

##### Declaration of Interests.

The authors report no conflict of interest.

##  Appendix A The entropy density for the LR26 equations from the Boltzmann
entropy density

The Boltzmann entropy density $\hat{\eta}$ and entropy flux $\hat{\Gamma}_{i}$
(in dimensional form) are given by (24) and (25). We replace $\hat{f}$ in (24)
and (25) with the G26 distribution function (Struchtrup, 2005; Gupta &
Torrilhon, 2012; Gupta, 2015)

$\displaystyle\hat{f}_{\mathrm{G26}}=\hat{f}_{\mathrm{M}}\left(1+\Phi\right),$

where

$\displaystyle\hat{f}_{\mathrm{M}}=\frac{\hat{\rho}}{(2\pi\hat{\theta})^{3/2}}\exp\left(-\frac{\hat{C}^{2}}{2\hat{\theta}}\right),$
(78)

is the Maxwellian distribution function and

$\displaystyle\Phi$
$\displaystyle=\frac{\hat{\sigma}_{ij}}{2}\frac{\hat{C}_{i}\hat{C}_{j}}{\hat{\rho}\hat{\theta}^{2}}-\frac{\hat{q}_{i}\hat{C}_{i}}{\hat{\rho}\hat{\theta}^{2}}\left(1-\frac{\hat{C}^{2}}{5\hat{\theta}}\right)+\frac{\hat{m}_{ijk}}{6\hat{\rho}\hat{\theta}^{3}}\hat{C}_{i}\hat{C}_{j}\hat{C}_{k}$
$\displaystyle\quad+\frac{\hat{\Delta}}{8\hat{\rho}\hat{\theta}^{2}}\left(1-\frac{2\hat{C}^{2}}{3\hat{\theta}}+\frac{\hat{C}^{4}}{15\hat{\theta}^{2}}\right)-\frac{\hat{R}_{ij}}{4\hat{\rho}\hat{\theta}^{3}}\left(1-\frac{\hat{C}^{2}}{7\hat{\theta}}\right)\hat{C}_{i}\hat{C}_{j}$

is the dimensionless deviation of the distribution function from the
equilibrium state. Here, $\hat{C}_{i}=\hat{c}_{i}-\hat{v}_{i}$ is the
(dimensional) peculiar velocity with $\hat{c}_{i}$ and $\hat{v}_{i}$ being the
(dimensional) microscopic and macroscopic velocities, respectively.

For the entropy density and entropy flux of the LR26 equations, the
integrations in (24) and (25) need to be evaluated by expending
$\hat{f}_{\mathrm{G26}}$ up to second order in deviations. For this purpose,
we introduce a smallness parameter $\varepsilon$, whose power in a term tells
the order of the deviation. The parameter is introduced just for the
bookkeeping purpose and will be replaced with unity at the end of
computations. Each moment is replaced with its value in the equilibrium plus
its deviation from its equilibrium value. We write

$\displaystyle\left.\begin{gathered}\hat{\rho}:=\hat{\rho}_{0}(1+\varepsilon\rho),\quad\hat{v}_{i}:=\varepsilon\sqrt{\hat{\theta}_{0}}\,v_{i},\quad\hat{\theta}:=\hat{\theta}_{0}(1+\varepsilon
T),\quad\hat{\sigma}_{ij}:=\varepsilon\hat{\rho}_{0}\hat{\theta}_{0}\sigma_{ij}\\\
\hat{q}_{i}:=\varepsilon\hat{\rho}_{0}\hat{\theta}_{0}^{3/2}q_{i},\quad\hat{m}_{ijk}:=\varepsilon\hat{\rho}_{0}\hat{\theta}_{0}^{3/2}m_{ijk},\quad\hat{\Delta}:=\varepsilon\hat{\rho}_{0}\hat{\theta}_{0}^{2}\Delta,\quad\hat{R}_{ij}:=\varepsilon\hat{\rho}_{0}\hat{\theta}_{0}^{2}R_{ij},\end{gathered}\right\\}$
(81)

where $\rho$, $v_{i}$, $T$, $\sigma_{ij}$, $q_{i}$, $m_{ijk}$, $\Delta$ and
$R_{ij}$ are the dimensionless deviations of the respective quantities from
their equilibrium values. The microscopic velocity $\hat{c}_{i}$ is also
scaled to give the dimensionless microscopic velocity
$\hat{c}_{i}=\sqrt{\hat{\theta}_{0}}\,c_{i}$.

A Taylor series expansion of $\hat{f}_{\mathrm{M}}$, to second order in
$\varepsilon$, yields

$\displaystyle\hat{f}_{\mathrm{M}}\approx\hat{f}_{M}^{(0)}\left[1+\varepsilon
f_{M}^{(1)}+\varepsilon^{2}f_{M}^{(2)}\right],$ (82)

where

$\displaystyle\hat{f}_{M}^{(0)}=\frac{\hat{\rho}_{0}}{(2\pi\hat{\theta}_{0})^{3/2}}\exp\left(-\frac{c^{2}}{2}\right)$
(83)

is the Maxwellian distribution function in the equilibrium, and

$\displaystyle f_{M}^{(1)}$
$\displaystyle=\rho-\frac{3T}{2}+\frac{Tc^{2}}{2}+c_{j}v_{j},$ (84)
$\displaystyle f_{M}^{(2)}$ $\displaystyle=-\frac{3\rho
T}{2}+\frac{15T^{2}}{8}-\frac{v^{2}}{2}+\left(\rho
T-\frac{5T^{2}}{2}+v^{2}+\frac{T^{2}c^{2}}{4}\right)\frac{c^{2}}{2}+\left(\rho-\frac{5T}{2}+\frac{Tc^{2}}{2}\right)c_{j}v_{j}.$
(85)

Similarly, a Taylor series expansion of $\Phi$, to second order in
$\varepsilon$, yields

$\displaystyle\Phi\approx\varepsilon\Phi^{(1)}+\varepsilon^{2}\Phi^{(2)},$
(86)

where

$\displaystyle\Phi^{(1)}$
$\displaystyle=\frac{\sigma_{kl}}{2}c_{k}c_{l}-q_{k}c_{k}\left(1-\frac{c^{2}}{5}\right)+\frac{m_{klr}}{6}c_{k}c_{l}c_{r}$
$\displaystyle\quad\;+\frac{\Delta}{8}\left(1-\frac{2c^{2}}{3}+\frac{c^{4}}{15}\right)-\frac{R_{kl}}{4}\left(1-\frac{c^{2}}{7}\right)c_{k}c_{l},$
(87) $\displaystyle\Phi^{(2)}$
$\displaystyle=-\frac{\sigma_{kl}}{2}\Big{[}(\rho+2T)c_{k}c_{l}+c_{k}v_{l}+c_{l}v_{k}\Big{]}$
$\displaystyle\quad\;+q_{k}\left[(\rho+2T)c_{k}-\frac{1}{5}(\rho+3T)c^{2}c_{k}+\left(1-\frac{c^{2}}{5}\right)v_{k}-\frac{2}{5}c_{k}c_{l}v_{l}\right]$
$\displaystyle\quad\;-\frac{m_{klr}}{6}\Big{[}(\rho+3T)c_{k}c_{l}c_{r}+c_{k}c_{l}v_{r}+c_{k}c_{r}v_{l}+c_{l}c_{r}v_{k}\Big{]}$
$\displaystyle\quad\;-\frac{\Delta}{8}\left[\rho+2T-\frac{2}{3}(\rho+3T)c^{2}+\frac{1}{15}(\rho+4T)c^{4}-\frac{4}{3}\left(1-\frac{c^{2}}{5}\right)c_{k}v_{k}\right]$
$\displaystyle\quad\;\,+\frac{R_{kl}}{4}\left[\left(1-\frac{c^{2}}{7}\right)\Big{\\{}(\rho+3T)c_{k}c_{l}+c_{k}v_{l}+c_{l}v_{k}\Big{\\}}-\frac{T}{7}c^{2}c_{k}c_{l}-\frac{2}{7}c_{k}c_{l}c_{r}v_{r}\right].$
(88)

Therefore, to second order in $\varepsilon$, $\hat{f}_{\mathrm{G26}}$ can be
approximated as

$\displaystyle\hat{f}_{\mathrm{G26}}$
$\displaystyle\approx\hat{f}_{M}^{(0)}\left[1+\varepsilon\phi^{(1)}+\varepsilon^{2}\phi^{(2)}\right],$
(89)

where

$\displaystyle\phi^{(1)}=\Phi^{(1)}+f_{M}^{(1)}\quad\text{and}\quad\phi^{(2)}=\Phi^{(2)}+\Phi^{(1)}f_{M}^{(1)}+f_{M}^{(2)}.$
(90)

Let $\hat{y}=\hat{f}_{M}^{(0)}$. Hence, to second order in $\varepsilon$,

$\displaystyle\ln{\frac{\hat{f}_{\mathrm{G26}}}{\hat{f}_{M}^{(0)}}}\approx\ln{\left(1+\varepsilon\phi^{(1)}+\varepsilon^{2}\phi^{(2)}\right)}\approx\varepsilon\phi^{(1)}+\varepsilon^{2}\left(\phi^{(2)}-\frac{\big{\\{}\phi^{(1)}\big{\\}}^{2}}{2}\right).$
(91)

From (24), the dimensionless entropy
$\eta=\hat{\eta}/(\hat{k}_{B}\hat{\rho}_{0})$, to second order in
$\varepsilon$, is given by

$\displaystyle\eta$
$\displaystyle=-\frac{1}{\hat{\rho}_{0}}\int\hat{f}_{\mathrm{G26}}\left(\ln{\frac{\hat{f}_{\mathrm{G26}}}{\hat{f}_{M}^{(0)}}}-1\right)\mathrm{d}\hat{\bm{c}}$
$\displaystyle\approx-\frac{1}{\hat{\rho}_{0}}\int\hat{f}_{M}^{(0)}\left(1+\varepsilon\phi^{(1)}+\varepsilon^{2}\phi^{(2)}\right)\left[\varepsilon\phi^{(1)}+\varepsilon^{2}\left(\phi^{(2)}-\frac{\big{\\{}\phi^{(1)}\big{\\}}^{2}}{2}\right)-1\right](\hat{\theta}_{0}^{3/2}\mathrm{d}\bm{c})$
$\displaystyle=\frac{1}{(2\pi)^{3/2}}\int\exp{\left(-\frac{c^{2}}{2}\right)}\left(1-\varepsilon^{2}\frac{\big{\\{}\phi^{(1)}\big{\\}}^{2}}{2}\right)\mathrm{d}\bm{c}.$
(92)

The integration over the velocity space is performed in a straightforward way,
which yields

$\displaystyle\eta=1+\varepsilon^{2}\left(-\frac{1}{2}\rho^{2}-\frac{1}{2}v^{2}-\frac{3}{4}T^{2}-\frac{1}{4}\sigma^{2}-\frac{1}{5}q^{2}-\frac{1}{12}m^{2}-\frac{1}{56}R^{2}-\frac{1}{240}\Delta^{2}\right).$
(93)

Substituting $\varepsilon=1$, one obtains (26) with $a_{0}=1$.

From (25), the dimensionless entropy flux
$\Gamma_{i}=\hat{\Gamma}_{i}/(\hat{k}_{B}\hat{\rho}_{0}\sqrt{\hat{\theta}_{0}})$,
to second order in $\varepsilon$, is given by

$\displaystyle\Gamma_{i}$
$\displaystyle=-\frac{1}{\hat{\rho}_{0}\sqrt{\hat{\theta}_{0}}}\int\hat{f}_{\mathrm{G26}}\left(\ln{\frac{\hat{f}_{\mathrm{G26}}}{\hat{f}_{M}^{(0)}}}-1\right)\hat{c}_{i}\,\mathrm{d}\hat{\bm{c}}$
$\displaystyle\approx-\frac{\hat{\theta}_{0}^{3/2}}{\hat{\rho}_{0}}\int\hat{f}_{M}^{(0)}\left(1+\varepsilon\phi^{(1)}+\varepsilon^{2}\phi^{(2)}\right)\left[\varepsilon\phi^{(1)}+\varepsilon^{2}\left(\phi^{(2)}-\frac{\big{\\{}\phi^{(1)}\big{\\}}^{2}}{2}\right)-1\right]c_{i}\mathrm{d}\bm{c}$
$\displaystyle=\frac{1}{(2\pi)^{3/2}}\int\exp{\left(-\frac{c^{2}}{2}\right)}\left(1-\varepsilon^{2}\frac{\big{\\{}\phi^{(1)}\big{\\}}^{2}}{2}\right)c_{i}\,\mathrm{d}\bm{c}.$
(94)

Again, the integration over the velocity space is performed in a
straightforward way, which yields

$\displaystyle\Gamma_{i}$
$\displaystyle=-\varepsilon^{2}\left[(\rho+T)v_{i}+Tq_{i}+v_{j}\sigma_{ij}+\frac{2}{5}\sigma_{ij}q_{j}+\frac{1}{2}m_{ijk}\sigma_{jk}\right.$
$\displaystyle\left.\quad\;+\frac{1}{5}R_{ij}q_{j}+\frac{1}{15}q_{i}\Delta+\frac{1}{14}m_{ijk}R_{jk}\right].$
(95)

Substituting $\varepsilon=1$ and $\rho+T=p$, one obtains

$\displaystyle\Gamma_{i}$
$\displaystyle=-\left(pv_{i}+Tq_{i}+v_{j}\sigma_{ij}+\frac{2}{5}\sigma_{ij}q_{j}+\frac{1}{2}m_{ijk}\sigma_{jk}+\frac{1}{5}R_{ij}q_{j}+\frac{1}{15}q_{i}\Delta+\frac{1}{14}m_{ijk}R_{jk}\right).$
(96)

Comparing (96) and (3), it is evident that except for the terms containing the
higher-order moments, namely $\Phi_{ijkl}$, $\Psi_{ijk}$, $\Omega_{j}$,
$\Omega_{i}$, all others terms in (96) and (3) are exactly the same.

##  Appendix B The entropy density for the linearised R26 equations

Let us assume that the entropy density for the LR26 equations is given by
(27). Taking the partial derivative of (27) with respect to $t$, one obtains

$\displaystyle\frac{\partial\eta}{\partial t}=$ $\displaystyle
a_{1}\rho\frac{\partial\rho}{\partial t}+a_{2}v_{i}\frac{\partial
v_{i}}{\partial t}+a_{3}T\frac{\partial T}{\partial
t}+a_{4}\sigma_{ij}\frac{\partial\sigma_{ij}}{\partial
t}+a_{5}q_{i}\frac{\partial q_{i}}{\partial t}$
$\displaystyle+a_{6}m_{ijk}\frac{\partial m_{ijk}}{\partial
t}+a_{7}R_{ij}\frac{\partial R_{ij}}{\partial
t}+a_{8}\Delta\frac{\partial\Delta}{\partial t}.$ (97)

Replacing the time derivatives of field variables in the above equation using
the G26 equations (7)–(11) (with $F_{i}=0$ and $p\approx\rho+T$) and
(15)–(17), the time derivative of $\eta$ can be simplified as follows:

$\displaystyle\frac{\partial\eta}{\partial t}=$ $\displaystyle-
a_{1}\rho\frac{\partial v_{i}}{\partial
x_{i}}-a_{2}v_{i}\left(\frac{\partial\sigma_{ij}}{\partial
x_{j}}+\frac{\partial\rho}{\partial x_{i}}+\frac{\partial T}{\partial
x_{i}}\right)-\frac{2}{3}a_{3}T\left(\frac{\partial q_{i}}{\partial
x_{i}}+\frac{\partial v_{i}}{\partial x_{i}}\right)$ $\displaystyle-
a_{4}\sigma_{ij}\left(\frac{\partial m_{ijk}}{\partial
x_{k}}+\frac{4}{5}\frac{\partial q_{\langle i}}{\partial
x_{j\rangle}}+\frac{1}{\mathrm{Kn}}\sigma_{ij}+2\frac{\partial v_{\langle
i}}{\partial x_{j\rangle}}\right)$ $\displaystyle-
a_{5}q_{i}\left(\frac{1}{2}\frac{\partial R_{ij}}{\partial
x_{j}}+\frac{1}{6}\frac{\partial\Delta}{\partial
x_{i}}+\frac{\partial\sigma_{ij}}{\partial
x_{j}}+\frac{\mathrm{Pr}}{\mathrm{Kn}}q_{i}+\frac{5}{2}\frac{\partial
T}{\partial x_{i}}\right)$ $\displaystyle-
a_{6}m_{ijk}\left(\frac{\partial\Phi_{ijkl}}{\partial
x_{l}}+\frac{3}{7}\frac{\partial R_{\langle ij}}{\partial
x_{k\rangle}}+\frac{\mathrm{Pr}_{m}}{\mathrm{Kn}}m_{ijk}+3\frac{\partial\sigma_{\langle
ij}}{\partial x_{k\rangle}}\right)$ $\displaystyle-
a_{7}R_{ij}\left(\frac{\partial\Psi_{ijk}}{\partial x_{k}}+2\frac{\partial
m_{ijk}}{\partial x_{k}}+\frac{2}{5}\frac{\partial\Omega_{\langle i}}{\partial
x_{j\rangle}}+\frac{\mathrm{Pr}_{R}}{\mathrm{Kn}}R_{ij}+\frac{28}{5}\frac{\partial
q_{\langle i}}{\partial x_{j\rangle}}\right)$ $\displaystyle-
a_{8}\Delta\left(\frac{\partial\Omega_{i}}{\partial
x_{i}}+\frac{\mathrm{Pr}_{\Delta}}{\mathrm{Kn}}\Delta+8\frac{\partial
q_{i}}{\partial x_{i}}\right)$ (98)

or

$\displaystyle\frac{\partial\eta}{\partial t}=$
$\displaystyle-\left(a_{1}\rho\frac{\partial v_{i}}{\partial
x_{i}}+a_{2}v_{i}\frac{\partial\rho}{\partial
x_{i}}\right)-\left(a_{2}v_{i}\frac{\partial T}{\partial
x_{i}}+\frac{2}{3}a_{3}T\frac{\partial v_{i}}{\partial
x_{i}}\right)-\left(a_{2}v_{i}\frac{\partial\sigma_{ij}}{\partial
x_{j}}+2a_{4}\sigma_{ij}\frac{\partial v_{i}}{\partial x_{j}}\right)$
$\displaystyle-\left(\frac{2}{3}a_{3}T\frac{\partial q_{i}}{\partial
x_{i}}+\frac{5}{2}a_{5}q_{i}\frac{\partial T}{\partial
x_{i}}\right)-\left(a_{4}\sigma_{ij}\frac{\partial m_{ijk}}{\partial
x_{k}}+3a_{6}m_{ijk}\frac{\partial\sigma_{ij}}{\partial x_{k}}\right)$
$\displaystyle-\left(\frac{4}{5}a_{4}\sigma_{ij}\frac{\partial q_{i}}{\partial
x_{j}}+a_{5}q_{i}\frac{\partial\sigma_{ij}}{\partial
x_{j}}\right)-\left(\frac{1}{2}a_{5}q_{i}\frac{\partial R_{ij}}{\partial
x_{j}}+\frac{28}{5}a_{7}R_{ij}\frac{\partial q_{i}}{\partial x_{j}}\right)$
$\displaystyle+\left(\frac{1}{6}a_{5}q_{i}\frac{\partial\Delta}{\partial
x_{i}}+8a_{8}\Delta\frac{\partial q_{i}}{\partial
x_{i}}\right)-\left(\frac{3}{7}a_{6}m_{ijk}\frac{\partial R_{ij}}{\partial
x_{k}}+2a_{7}R_{ij}\frac{\partial m_{ijk}}{\partial x_{k}}\right)$
$\displaystyle-a_{6}m_{ijk}\frac{\partial\Phi_{ijkl}}{\partial
x_{l}}-a_{7}R_{ij}\frac{\partial\Psi_{ijk}}{\partial
x_{k}}-\frac{2}{5}a_{7}R_{ij}\frac{\partial\Omega_{\langle i}}{\partial
x_{j\rangle}}-a_{8}\Delta\frac{\partial\Omega_{i}}{\partial x_{i}}$
$\displaystyle-\frac{a_{4}}{\mathrm{Kn}}\sigma^{2}-a_{5}\frac{\mathrm{Pr}}{\mathrm{Kn}}q^{2}-a_{6}\frac{\mathrm{Pr}_{m}}{\mathrm{Kn}}m^{2}-a_{7}\frac{\mathrm{Pr}_{R}}{\mathrm{Kn}}R^{2}-a_{8}\frac{\mathrm{Pr}_{\Delta}}{\mathrm{Kn}}\Delta^{2}.$
(99)

The goal is to write the above equation in the form of (23), which can be
accomplished by taking

$\displaystyle\left.\begin{gathered}a_{1}=a_{2},\quad
a_{2}=\frac{2}{3}a_{3},\quad
a_{2}=2a_{4},\quad\frac{2}{3}a_{3}=\frac{5}{2}a_{5},\quad a_{4}=3a_{6},\\\
\frac{4}{5}a_{4}=a_{5},\quad\frac{1}{2}a_{5}=\frac{28}{5}a_{7},\quad\frac{1}{6}a_{5}=8a_{8},\quad\frac{3}{7}a_{6}=2a_{7}.\end{gathered}\right\\}$
(102)

These are nine equations in eight variables. Nevertheless, they are
consistent, and only seven of them are linearly independent. Thus, on taking
$a_{1}=-1$, the other coefficients are

$\displaystyle a_{2}=-1,\quad a_{3}=-\frac{3}{2},\quad
a_{4}=-\frac{1}{2},\quad a_{5}=-\frac{2}{5},\quad a_{6}=-\frac{1}{6},\quad
a_{7}=-\frac{1}{28},\quad a_{8}=-\frac{1}{120}.$ (103)

Note that the constant $a_{1}$ has been taken as negative in order to get the
non-negative entropy production rate. Consequently, the entropy (27) becomes
(26).

##  Appendix C Entropy flux at the interface

The entropy flux at the interface is given by (32). Substituting
$\Gamma_{i}^{\mathrm{(gas)}}=\Gamma_{i}$ from (3) and
$\Gamma_{i}^{\mathrm{(liquid)}}$ from (33) into (32), the entropy flux at the
interface reads

$\displaystyle\Sigma^{I}$
$\displaystyle=-\underline{\left(pv_{i}-p^{\ell}v_{i}^{\ell}+Tq_{i}-T^{\ell}q_{i}^{\ell}+v_{j}\sigma_{ij}-v_{j}^{\ell}\sigma_{ij}^{\ell}\right)n_{i}}$
$\displaystyle\quad\,-\left(\frac{2}{5}\sigma_{ij}q_{j}+\frac{1}{2}m_{ijk}\sigma_{jk}+\frac{1}{5}R_{ij}q_{j}+\frac{1}{15}q_{i}\Delta+\frac{1}{14}m_{ijk}R_{jk}\right.$
$\displaystyle\left.\quad\,+\frac{1}{6}\Phi_{ijkl}m_{jkl}+\frac{1}{28}\Psi_{ijk}R_{jk}+\frac{1}{70}R_{ij}\Omega_{j}+\frac{1}{120}\Omega_{i}\Delta\right)n_{i}.$
(104)

Let us refer the underlined term in (C) by $\underline{\Sigma}^{I}$, i.e.

$\displaystyle\underline{\Sigma}^{I}$
$\displaystyle=\left(pv_{i}-p^{\ell}v_{i}^{\ell}+Tq_{i}-T^{\ell}q_{i}^{\ell}+v_{j}\sigma_{ij}-v_{j}^{\ell}\sigma_{ij}^{\ell}\right)n_{i}$
$\displaystyle=pv_{i}n_{i}-p^{\ell}v_{i}^{\ell}n_{i}+Tq_{i}n_{i}-T^{\ell}q_{i}^{\ell}n_{i}+v_{i}\sigma_{ij}n_{j}-v_{i}^{\ell}\sigma_{ij}^{\ell}n_{j}.$
(105)

Equation (105) is simplified using the conservation laws at the interface. For
a simple interface (assumed here), the conservation laws for the mass,
momentum and energy apprise the continuity of the normal fluxes (Struchtrup et
al., 2017). With the gas (or vapour) properties on the left-hand sides and the
liquid properties on the right-hand sides, the conservation laws at the
interface in linear-dimensionless form read (Müller, 1985; Struchtrup et al.,
2017)

$\displaystyle\left(v_{i}-v_{i}^{I}\right)n_{i}$
$\displaystyle=\rho^{\ell}\left(v_{i}^{\ell}-v_{i}^{I}\right)n_{i},$ (106)
$\displaystyle\left(p\delta_{ij}+\sigma_{ij}\right)n_{j}$
$\displaystyle=\left(p^{\ell}\delta_{ij}+\sigma_{ij}^{\ell}\right)n_{j},$
(107) $\displaystyle\left[h_{0}\left(v_{i}-v_{i}^{I}\right)+q_{i}\right]n_{i}$
$\displaystyle=\left[\rho^{\ell}h_{0}^{\ell}\left(v_{i}^{\ell}-v_{i}^{I}\right)+q_{i}^{\ell}\right]n_{i},$
(108)

where $v_{i}^{I}$ is the velocity of the interface, and $h_{0}$ and
$h_{0}^{\ell}$ are the dimensionless enthalpies of the vapour and liquid,
respectively, at temperature $\hat{T}_{0}$. The variables in (106)–(108) have
been non-dimensionalised using $\hat{\rho}_{0}$, $\hat{\theta}_{0}$ and their
appropriate combinations. Substituting $\sigma_{ij}^{\ell}n_{j}$ and
$q_{i}^{\ell}n_{i}$ from (107) and (108) in (105), the underlined term in (C)
becomes

$\displaystyle\underline{\Sigma}^{I}$
$\displaystyle=p\left(v_{i}-v_{i}^{\ell}\right)n_{i}+\left(T-T^{\ell}\right)q_{i}n_{i}-T^{\ell}\left(h_{0}-h_{0}^{\ell}\right)\left(v_{i}-v_{i}^{I}\right)n_{i}+\left(v_{j}-v_{j}^{\ell}\right)\sigma_{ij}n_{i}$
$\displaystyle=p\mathscr{V}_{i}n_{i}+\mathcal{T}q_{i}n_{i}-T^{\ell}\left(h_{0}-h_{0}^{\ell}\right)\mathcal{V}_{i}n_{i}+\mathscr{V}_{j}\sigma_{ij}n_{i}.$
(109)

where $\mathscr{V}_{i}=v_{i}-v_{i}^{\ell}$, $\mathcal{T}=T-T^{\ell}$,
$\mathcal{V}_{i}=v_{i}-v_{i}^{I}$. On using (106), it turns out that

$\displaystyle\mathscr{V}_{i}\,n_{i}=\left(v_{i}-v_{i}^{\ell}\right)n_{i}=\left[\mathcal{V}_{i}-\left(v_{i}^{\ell}-v_{i}^{I}\right)\right]n_{i}=\left(1-\frac{1}{\rho^{\ell}}\right)\mathcal{V}_{i}n_{i}\approx\mathcal{V}_{i}\,n_{i}$
(110)

because the density of a gas/vapour is usually much smaller than the density
of a liquid/solid, and hence $\rho^{\ell}\gg 1$ or $1/\rho^{\ell}\ll 1$. The
enthalpy term in (109) is eliminated using the definition of the
(dimensionless) saturation pressure $p_{\mathrm{sat}}$ at temperature
$T^{\ell}$ (Rana et al., 2019),

$\displaystyle p_{\mathrm{sat}}\equiv
p_{\mathrm{sat}}(T^{\ell})=T^{\ell}\left(h_{0}-h_{0}^{\ell}\right).$ (111)

Using (110) and (111), equation (109) reduces to

$\displaystyle\underline{\Sigma}^{I}$
$\displaystyle=p\mathcal{V}_{i}\,n_{i}+\mathcal{T}q_{i}n_{i}-p_{\mathrm{sat}}\mathcal{V}_{i}\,n_{i}+\mathscr{V}_{j}\sigma_{ij}n_{i}=\left(\mathcal{P}\mathcal{V}_{i}+\mathcal{T}q_{i}+\mathscr{V}_{j}\sigma_{ij}\right)n_{i},$
(112)

where $\mathcal{P}=p-p_{\mathrm{sat}}$. Replacing the underlined terms in (C)
with the value of $\underline{\Sigma}^{I}$ obtained in (112), the entropy flux
at the interface is given by (4).

##  Appendix D Decomposition of symmetric-tracefree tensors

In the following, we decompose vectors and symmetric-tracefree tensors of rank
two and three into their components in the direction of a given normal
$\bm{n}$ and in the directions perpendicular to it (mutually perpendicular
tangential directions, let us say $\bm{t}_{1}$ and $\bm{t}_{2}$). Obviously,
$n_{i}{t_{1}}_{i}=n_{i}{t_{2}}_{i}={t_{1}}_{i}{t_{2}}_{i}=0$. When it is not
required to distinguish between $\bm{t}_{1}$ and $\bm{t}_{2}$, both are
referred to with a common notation $\bm{t}$ for simplicity, and $\bm{t}$ could
either be $\bm{t}_{1}$ and $\bm{t}_{2}$.

  * •

A vector $a_{i}$ is decomposed as

$\displaystyle a_{i}=a_{n}n_{i}+\bar{a}_{i},$ (113)

where $a_{n}n_{i}$ with $a_{n}:=a_{k}n_{k}$ is the normal component of $a_{i}$
while $\bar{a}_{i}:=a_{i}-a_{n}n_{i}$ is the tangential component of $a_{i}$.
By definition, $\bar{a}_{i}$ is such that $\bar{a}_{i}n_{i}=0$. Furthermore,
multiplication of (113) with $t_{i}\in\\{{t_{1}}_{i},{t_{2}}_{i}\\}$ gives

$\displaystyle a_{t}=\bar{a}_{t},$ (114)

where $a_{t}:=a_{i}t_{i}$, $\bar{a}_{t}:=\bar{a}_{i}t_{i}$.

  * •

A symmetric-tracefree rank-2 tensor $A_{ij}$ is decomposed as (Rana &
Struchtrup, 2016; Rana et al., 2018a)

$\displaystyle A_{ij}=\frac{3}{2}A_{nn}\,n_{\langle
i}n_{j\rangle}+\bar{A}_{ni}n_{j}+\bar{A}_{nj}n_{i}+\tilde{A}_{ij},$ (115)

where Einstein summation never applies to the indices ‘$n$’, and

$\displaystyle\left.\begin{aligned} A_{nn}&:=A_{ij}n_{i}n_{j},\\\
\bar{A}_{ni}&:=A_{ij}n_{j}-A_{nn}n_{i},\\\
\tilde{A}_{ij}&:=A_{ij}-\frac{3}{2}A_{nn}\,n_{\langle
i}n_{j\rangle}-\bar{A}_{ni}n_{j}-\bar{A}_{nj}n_{i}.\end{aligned}\right\\}$
(116)

It is easy to verify that $\bar{A}_{ni}n_{i}=0$ and
$\tilde{A}_{kk}=\tilde{A}_{ij}n_{i}=\tilde{A}_{ij}n_{j}=0$. Furthermore,
multiplication of (116)2 with $t_{i}\in\\{{t_{1}}_{i},{t_{2}}_{i}\\}$ gives

$\displaystyle\bar{A}_{nt}=A_{nt},$ (117)

where $\bar{A}_{nt}:=\bar{A}_{ni}t_{i}$, $A_{nt}=A_{tn}:=A_{ij}t_{i}n_{j}$;
multiplication of (116)3 with $t_{i}t_{j}$ gives

$\displaystyle\tilde{A}_{tt}$ $\displaystyle=A_{tt}+\frac{1}{2}A_{nn},$ (118)

where $\tilde{A}_{tt}:=\tilde{A}_{ij}t_{i}t_{j}$ and
$A_{tt}:=A_{ij}t_{i}t_{j}$; multiplication of (116)3 with
${t_{1}}_{i}{t_{2}}_{j}$ gives

$\displaystyle\tilde{A}_{{t_{1}}{t_{2}}}$ $\displaystyle=A_{{t_{1}}{t_{2}}},$
(119)

where $\tilde{A}_{{t_{1}}{t_{2}}}:=\tilde{A}_{ij}{t_{1}}_{i}{t_{2}}_{j}$ and
$A_{{t_{1}}{t_{2}}}:=A_{ij}{t_{1}}_{i}{t_{2}}_{j}$.

  * •

A symmetric-tracefree rank-3 tensor $A_{ijk}$ is decomposed as

$\displaystyle A_{ijk}$ $\displaystyle=\frac{5}{2}A_{nnn}n_{\langle
i}n_{j}n_{k\rangle}+\frac{15}{4}\bar{A}_{nn\langle
i}n_{j}n_{k\rangle}+3\tilde{A}_{n\langle ij}n_{k\rangle}+\check{A}_{ijk}.$
(120)

where

$\displaystyle\left.\begin{aligned} A_{nnn}&:=A_{ijk}n_{i}n_{j}n_{k},\\\
\bar{A}_{nni}&:=A_{ijk}n_{j}n_{k}-A_{nnn}n_{i},\\\
\tilde{A}_{nij}&:=A_{ijk}n_{k}-\frac{3}{2}A_{nnn}n_{\langle
i}n_{j\rangle}-\bar{A}_{nni}n_{j}-\bar{A}_{nnj}n_{i},\\\
\check{A}_{ijk}&:=A_{ijk}-\frac{5}{2}A_{nnn}n_{\langle
i}n_{j}n_{k\rangle}-\frac{15}{4}\bar{A}_{nn\langle
i}n_{j}n_{k\rangle}-3\tilde{A}_{n\langle
ij}n_{k\rangle}.\end{aligned}\right\\}$ (121)

It is again easy to verify that
$\bar{A}_{nni}n_{i}=\tilde{A}_{nij}n_{j}=\tilde{A}_{nii}=\check{A}_{ijk}n_{k}=\check{A}_{ikk}=0$.
Furthermore, multiplication of (121)2 with
$t_{i}\in\\{{t_{1}}_{i},{t_{2}}_{i}\\}$ gives

$\displaystyle\bar{A}_{nnt}=A_{nnt},$ (122)

where $\bar{A}_{nnt}:=\bar{A}_{nni}t_{i}$,
$A_{nnt}=A_{tnn}:=A_{ijk}t_{i}n_{j}n_{k}$; multiplication of (121)3 with
$t_{i}t_{j}$ gives

$\displaystyle\tilde{A}_{ntt}$ $\displaystyle=A_{ntt}+\frac{1}{2}A_{nnn},$
(123)

where $\tilde{A}_{ntt}:=\tilde{A}_{nij}t_{i}t_{j}$ and
$A_{ntt}=A_{ttn}:=A_{ijk}t_{i}t_{j}n_{k}$; multiplication of (121)3 with
${t_{1}}_{i}{t_{2}}_{j}$ gives

$\displaystyle\tilde{A}_{n{t_{1}}{t_{2}}}$
$\displaystyle=A_{n{t_{1}}{t_{2}}},$ (124)

where $\tilde{A}_{n{t_{1}}{t_{2}}}:=\tilde{A}_{nij}{t_{1}}_{i}{t_{2}}_{j}$ and
$A_{n{t_{1}}{t_{2}}}=A_{{t_{1}}{t_{2}}n}:=A_{ijk}{t_{1}}_{i}{t_{2}}_{j}n_{k}$;
multiplication of (121)4 with $t_{i}t_{j}t_{k}$ gives

$\displaystyle\check{A}_{ttt}=A_{ttt}+\frac{3}{4}A_{nnt},$ (125)

where $\check{A}_{ttt}:=\check{A}_{ijk}t_{i}t_{j}t_{k}$ and
$A_{ttt}:=A_{ijk}t_{i}t_{j}t_{k}$; multiplication of (121)4 with
${t_{1}}_{i}{t_{1}}_{j}{t_{2}}_{k}$ gives

$\displaystyle\check{A}_{{t_{1}}{t_{1}}{t_{2}}}$
$\displaystyle=A_{{t_{1}}{t_{1}}{t_{2}}}+\frac{1}{4}A_{nn{t_{2}}},$ (126)

where
$\check{A}_{{t_{1}}{t_{1}}{t_{2}}}:=\check{A}_{ijk}{t_{1}}_{i}{t_{1}}_{j}{t_{2}}_{k}$
and $A_{{t_{1}}{t_{1}}{t_{2}}}:=A_{ijk}{t_{1}}_{i}{t_{1}}_{j}{t_{2}}_{k}$.

##  Appendix E Analytic solution for the unknowns in ansatz (75)–(77)

The analytic solution for the unknowns appearing in (75)–(77) are as follows:

$\displaystyle a_{1}(r)$
$\displaystyle=\frac{C_{1}}{2r}+\frac{C_{2}}{3r^{3}}-K_{1}\Theta_{1}(r)\left(\frac{0.975701\mathrm{Kn}^{3}}{r^{3}}+\frac{0.497885\mathrm{Kn}^{2}}{r^{2}}\right)$
$\displaystyle\quad-
K_{2}\Theta_{2}(r)\left(\frac{1.4893\mathrm{Kn}^{3}}{r^{3}}+\frac{1.88528\mathrm{Kn}^{2}}{r^{2}}\right),$
(127) $\displaystyle a_{2}(r)$
$\displaystyle=\frac{C_{1}}{4r}-\frac{C_{2}}{6r^{3}}+K_{1}\Theta_{1}(r)\left(\frac{0.487851\mathrm{Kn}^{3}}{r^{3}}+\frac{0.248943\mathrm{Kn}^{2}}{r^{2}}+\frac{0.127032\mathrm{Kn}}{r}\right)$
$\displaystyle\quad+K_{2}\Theta_{2}(r)\left(\frac{0.74465\mathrm{Kn}^{3}}{r^{3}}+\frac{0.942641\mathrm{Kn}^{2}}{r^{2}}+\frac{1.19327\mathrm{Kn}}{r}\right),$
(128) $\displaystyle\alpha_{1}(r)$
$\displaystyle=\frac{\mathrm{Kn}^{2}}{r^{3}}\left(\frac{C_{3}}{6\mathrm{Kn}}-\frac{3C_{1}}{2}\right)+K_{1}\Theta_{1}(r)\left(\frac{1.50825\mathrm{Kn}^{3}}{r^{3}}+\frac{0.769634\mathrm{Kn}^{2}}{r^{2}}\right)$
$\displaystyle\quad-
K_{2}\Theta_{2}(r)\left(\frac{0.80299\mathrm{Kn}^{3}}{r^{3}}+\frac{1.01649\mathrm{Kn}^{2}}{r^{2}}\right),$
(129) $\displaystyle\alpha_{2}(r)$
$\displaystyle=\frac{\mathrm{Kn}^{2}}{r^{3}}\left(\frac{3C_{1}}{4}-\frac{C_{3}}{12\mathrm{Kn}}\right)$
$\displaystyle\quad-
K_{1}\Theta_{1}(r)\left(\frac{0.754123\mathrm{Kn}^{3}}{r^{3}}+\frac{0.384817\mathrm{Kn}^{2}}{r^{2}}+\frac{0.196366\mathrm{Kn}}{r}\right)$
$\displaystyle\quad+K_{2}\Theta_{2}(r)\left(\frac{0.401495\mathrm{Kn}^{3}}{r^{3}}+\frac{0.508246\mathrm{Kn}^{2}}{r^{2}}+\frac{0.643381\mathrm{Kn}}{r}\right),$
(130) $\displaystyle c(r)$
$\displaystyle=\frac{C_{3}}{45r^{2}}-K_{3}\Theta_{3}(r)\left(\frac{0.0172036\mathrm{Kn}^{2}}{r^{2}}+\frac{0.0200113\mathrm{Kn}}{r}\right)$
$\displaystyle\quad-
K_{4}\Theta_{4}(r)\left(\frac{0.0629047\mathrm{Kn}^{2}}{r^{2}}+\frac{0.0426083\mathrm{Kn}}{r}\right)$
$\displaystyle\quad-
K_{5}\Theta_{5}(r)\left(\frac{0.443094\mathrm{Kn}^{2}}{r^{2}}+\frac{0.200538\mathrm{Kn}}{r}\right),$
(131) $\displaystyle d(r)$
$\displaystyle=\frac{C_{1}\mathrm{Kn}}{2r^{2}}-K_{3}\Theta_{3}(r)\left(\frac{0.931108\mathrm{Kn}^{2}}{r^{2}}+\frac{1.08307\mathrm{Kn}}{r}\right)$
$\displaystyle\quad-
K_{4}\Theta_{4}(r)\left(\frac{0.254173\mathrm{Kn}^{2}}{r^{2}}+\frac{0.172163\mathrm{Kn}}{r}\right)$
$\displaystyle\quad-
K_{5}\Theta_{5}(r)\left(\frac{0.0501203\mathrm{Kn}^{2}}{r^{2}}+\frac{0.0226838\mathrm{Kn}}{r}\right),$
(132) $\displaystyle\delta(r)$
$\displaystyle=K_{3}\Theta_{3}(r)\left(\frac{14.0055\mathrm{Kn}^{2}}{r^{2}}+\frac{16.2913\mathrm{Kn}}{r}\right)-K_{4}\Theta_{4}(r)\left(\frac{1.71119\mathrm{Kn}^{2}}{r^{2}}+\frac{1.15907\mathrm{Kn}}{r}\right)$
$\displaystyle\quad+K_{5}\Theta_{5}(r)\left(\frac{3.75265\mathrm{Kn}^{2}}{r^{2}}+\frac{1.6984\mathrm{Kn}}{r}\right).$
(133) $\displaystyle b_{1}(r)$
$\displaystyle=\frac{C_{1}\mathrm{Kn}}{r^{2}}+\frac{\mathrm{Kn}^{4}}{r^{4}}\left(-\frac{10C_{1}}{\mathrm{Kn}}+\frac{2C_{2}}{\mathrm{Kn}^{3}}+\frac{2C_{3}}{5\mathrm{Kn}^{2}}\right)$
$\displaystyle\quad+K_{3}\Theta_{3}(r)\left(\frac{6.19341\mathrm{Kn}^{4}}{r^{4}}+\frac{7.20421\mathrm{Kn}^{3}}{r^{3}}+\frac{3.72443\mathrm{Kn}^{2}}{r^{2}}+\frac{1.08307\mathrm{Kn}}{r}\right)$
$\displaystyle\quad+K_{4}\Theta_{4}(r)\left(\frac{4.98597\mathrm{Kn}^{4}}{r^{4}}+\frac{3.37723\mathrm{Kn}^{3}}{r^{3}}+\frac{1.01669\mathrm{Kn}^{2}}{r^{2}}+\frac{0.172163\mathrm{Kn}}{r}\right)$
$\displaystyle\quad+K_{5}\Theta_{5}(r)\left(\frac{2.20218\mathrm{Kn}^{4}}{r^{4}}+\frac{0.996676\mathrm{Kn}^{3}}{r^{3}}+\frac{0.200481\mathrm{Kn}^{2}}{r^{2}}+\frac{0.0226838\mathrm{Kn}}{r}\right),$
(134) $\displaystyle b_{2}(r)$
$\displaystyle=\frac{\mathrm{Kn}^{4}}{r^{4}}\left(-\frac{5C_{1}}{\mathrm{Kn}}+\frac{C_{2}}{\mathrm{Kn}^{3}}+\frac{\text{C3}}{5\mathrm{Kn}^{2}}\right)$
$\displaystyle\quad+K_{3}\Theta_{3}(r)\left(\frac{3.09671\mathrm{Kn}^{4}}{r^{4}}+\frac{3.6021\mathrm{Kn}^{3}}{r^{3}}+\frac{1.39666\mathrm{Kn}^{2}}{r^{2}}\right)$
$\displaystyle\quad+K_{4}\Theta_{4}(r)\left(\frac{2.49299\mathrm{Kn}^{4}}{r^{4}}+\frac{1.68862\mathrm{Kn}^{3}}{r^{3}}+\frac{0.381259\mathrm{Kn}^{2}}{r^{2}}\right)$
$\displaystyle\quad+K_{5}\Theta_{5}(r)\left(\frac{1.10109\mathrm{Kn}^{4}}{r^{4}}+\frac{0.498338\mathrm{Kn}^{3}}{r^{3}}+\frac{0.0751805\mathrm{Kn}^{2}}{r^{2}}\right),$
(135) $\displaystyle\beta_{1}(r)$
$\displaystyle=\frac{\mathrm{Kn}^{4}}{r^{4}}\left(\frac{12C_{3}}{5\mathrm{Kn}^{2}}-\frac{228C_{1}}{7\mathrm{Kn}}\right)$
$\displaystyle\quad+K_{1}\Theta_{1}(r)\left(\frac{46.338\mathrm{Kn}^{4}}{r^{4}}+\frac{23.6455\mathrm{Kn}^{3}}{r^{3}}+\frac{4.02199\mathrm{Kn}^{2}}{r^{2}}\right)$
$\displaystyle\quad-
K_{2}\Theta_{2}(r)\left(\frac{4.00877\mathrm{Kn}^{4}}{r^{4}}+\frac{5.07465\mathrm{Kn}^{3}}{r^{3}}+\frac{2.14131\mathrm{Kn}^{2}}{r^{2}}\right)$
$\displaystyle\quad-
K_{3}\Theta_{3}(r)\left(\frac{42.868\mathrm{Kn}^{4}}{r^{4}}+\frac{49.8643\mathrm{Kn}^{3}}{r^{3}}+\frac{25.7788\mathrm{Kn}^{2}}{r^{2}}+\frac{7.49652\mathrm{Kn}}{r}\right)$
$\displaystyle\quad+K_{4}\Theta_{4}(r)\left(\frac{7.38704\mathrm{Kn}^{4}}{r^{4}}+\frac{5.00359\mathrm{Kn}^{3}}{r^{3}}+\frac{1.5063\mathrm{Kn}^{2}}{r^{2}}+\frac{0.255071\mathrm{Kn}}{r}\right)$
$\displaystyle\quad+K_{5}\Theta_{5}(r)\left(\frac{37.9774\mathrm{Kn}^{4}}{r^{4}}+\frac{17.1881\mathrm{Kn}^{3}}{r^{3}}+\frac{3.45738\mathrm{Kn}^{2}}{r^{2}}+\frac{0.391191\mathrm{Kn}}{r}\right),$
(136) $\displaystyle\beta_{2}(r)$
$\displaystyle=\frac{\mathrm{Kn}^{4}}{r^{4}}\left(\frac{6C_{3}}{5\mathrm{Kn}^{2}}-\frac{114C_{1}}{7\mathrm{Kn}}\right)$
$\displaystyle\quad+K_{1}\Theta_{1}(r)\left(\frac{23.169\mathrm{Kn}^{4}}{r^{4}}+\frac{11.8228\mathrm{Kn}^{3}}{r^{3}}+\frac{3.01649\mathrm{Kn}^{2}}{r^{2}}+\frac{0.51309\mathrm{Kn}}{r}\right)$
$\displaystyle\quad-
K_{2}\Theta_{2}(r)\left(\frac{2.00439\mathrm{Kn}^{4}}{r^{4}}+\frac{2.53732\mathrm{Kn}^{3}}{r^{3}}+\frac{1.60598\mathrm{Kn}^{2}}{r^{2}}+\frac{0.677662\mathrm{Kn}}{r}\right)$
$\displaystyle\quad-
K_{3}\Theta_{3}(r)\left(\frac{21.434\mathrm{Kn}^{4}}{r^{4}}+\frac{24.9321\mathrm{Kn}^{3}}{r^{3}}+\frac{9.66706\mathrm{Kn}^{2}}{r^{2}}\right)$
$\displaystyle\quad+K_{4}\Theta_{4}(r)\left(\frac{3.69352\mathrm{Kn}^{4}}{r^{4}}+\frac{2.50179\mathrm{Kn}^{3}}{r^{3}}+\frac{0.564861\mathrm{Kn}^{2}}{r^{2}}\right)$
$\displaystyle\quad+K_{5}\Theta_{5}(r)\left(\frac{18.9887\mathrm{Kn}^{4}}{r^{4}}+\frac{8.59403\mathrm{Kn}^{3}}{r^{3}}+\frac{1.29652\mathrm{Kn}^{2}}{r^{2}}\right),$
(137) $\displaystyle\kappa_{1}(r)$
$\displaystyle=\frac{24C_{1}\mathrm{Kn}^{2}}{5r^{3}}+\frac{\mathrm{Kn}^{5}}{r^{5}}\left(-\frac{7344C_{1}}{49\mathrm{Kn}}+\frac{16C_{2}}{\mathrm{Kn}^{3}}+\frac{208C_{3}}{35\mathrm{Kn}^{2}}\right)$
$\displaystyle\quad+K_{1}\Theta_{1}(r)\left(\frac{64.3578\mathrm{Kn}^{5}}{r^{5}}+\frac{32.8408\mathrm{Kn}^{4}}{r^{4}}+\frac{6.70326\mathrm{Kn}^{3}}{r^{3}}+\frac{0.570095\mathrm{Kn}^{2}}{r^{2}}\right)$
$\displaystyle\quad+K_{2}\Theta_{2}(r)\left(\frac{50.8419\mathrm{Kn}^{5}}{r^{5}}+\frac{64.3599\mathrm{Kn}^{4}}{r^{4}}+\frac{32.5889\mathrm{Kn}^{3}}{r^{3}}+\frac{6.87563\mathrm{Kn}^{2}}{r^{2}}\right)$
$\displaystyle\quad+K_{3}\Theta_{3}(r)\left(\frac{30.5159\mathrm{Kn}^{5}}{r^{5}}+\frac{35.4962\mathrm{Kn}^{4}}{r^{4}}+\frac{18.5802\mathrm{Kn}^{3}}{r^{3}}+\frac{5.60327\mathrm{Kn}^{2}}{r^{2}}+\frac{0.931108\mathrm{Kn}}{r}\right)$
$\displaystyle\quad+K_{4}\Theta_{4}(r)\left(\frac{72.4497\mathrm{Kn}^{5}}{r^{5}}+\frac{49.0736\mathrm{Kn}^{4}}{r^{4}}+\frac{14.9579\mathrm{Kn}^{3}}{r^{3}}+\frac{2.62674\mathrm{Kn}^{2}}{r^{2}}+\frac{0.254173\mathrm{Kn}}{r}\right)$
$\displaystyle\quad+K_{5}\Theta_{5}(r)\left(\frac{71.6731\mathrm{Kn}^{5}}{r^{5}}+\frac{32.4383\mathrm{Kn}^{4}}{r^{4}}+\frac{6.60653\mathrm{Kn}^{3}}{r^{3}}+\frac{0.775193\mathrm{Kn}^{2}}{r^{2}}+\frac{0.0501203\mathrm{Kn}}{r}\right),$
(138) $\displaystyle\kappa_{2}(r)$
$\displaystyle=\frac{4C_{1}\mathrm{Kn}^{2}}{5r^{3}}+\frac{\mathrm{Kn}^{5}}{r^{5}}\left(-\frac{3672C_{1}}{49\mathrm{Kn}}+\frac{8C_{2}}{\mathrm{Kn}^{3}}+\frac{104C_{3}}{35\mathrm{Kn}^{2}}\right)$
$\displaystyle\quad+K_{1}\Theta_{1}(r)\left(\frac{32.1789\mathrm{Kn}^{5}}{r^{5}}+\frac{16.4204\mathrm{Kn}^{4}}{r^{4}}+\frac{3.91023\mathrm{Kn}^{3}}{r^{3}}+\frac{0.570095\mathrm{Kn}^{2}}{r^{2}}+\frac{0.0484851\mathrm{Kn}}{r}\right)$
$\displaystyle\quad+K_{2}\Theta_{2}(r)\left(\frac{25.4209\mathrm{Kn}^{5}}{r^{5}}+\frac{32.18\mathrm{Kn}^{4}}{r^{4}}+\frac{19.0102\mathrm{Kn}^{3}}{r^{3}}+\frac{6.87563\mathrm{Kn}^{2}}{r^{2}}+\frac{1.45063\mathrm{Kn}}{r}\right)$
$\displaystyle\quad+K_{3}\Theta_{3}(r)\left(\frac{15.2579\mathrm{Kn}^{5}}{r^{5}}+\frac{17.7481\mathrm{Kn}^{4}}{r^{4}}+\frac{8.25788\mathrm{Kn}^{3}}{r^{3}}+\frac{1.60094\mathrm{Kn}^{2}}{r^{2}}\right)$
$\displaystyle\quad+K_{4}\Theta_{4}(r)\left(\frac{36.2249\mathrm{Kn}^{5}}{r^{5}}+\frac{24.5368\mathrm{Kn}^{4}}{r^{4}}+\frac{6.64796\mathrm{Kn}^{3}}{r^{3}}+\frac{0.750496\mathrm{Kn}^{2}}{r^{2}}\right)$
$\displaystyle\quad+K_{5}\Theta_{5}(r)\left(\frac{35.8366\mathrm{Kn}^{5}}{r^{5}}+\frac{16.2192\mathrm{Kn}^{4}}{r^{4}}+\frac{2.93623\mathrm{Kn}^{3}}{r^{3}}+\frac{0.221484\mathrm{Kn}^{2}}{r^{2}}\right).$
(139)

In (E)–(E), the coefficients

$\displaystyle\left.\begin{gathered}\Theta_{1}(r)=\exp\left[-\dfrac{0.510285}{\mathrm{Kn}}(r-1)\right],\quad\Theta_{2}(r)=\exp\left[-\dfrac{1.26588}{\mathrm{Kn}}(r-1)\right],\\\
\Theta_{3}(r)=\exp\left[-\frac{1.16321}{\mathrm{Kn}}(r-1)\right],\quad\Theta_{4}(r)=\exp\left[-\frac{0.677347}{\mathrm{Kn}}(r-1)\right]\\\
\textrm{and}\quad\Theta_{5}(r)=\exp\left[-\frac{0.452587}{\mathrm{Kn}}(r-1)\right]\end{gathered}\right\\}$
(143)

describe the Knudsen layer functions that vanish as $r\to\infty$.

## References

  * Aoki & Sone (1987) Aoki, K. & Sone, Y. 1987 Temperature field induced around a sphere in a uniform flow of a rarefied gas.  Phys. Fluids 30, 2286–2288.
  * Bailey et al. (2005) Bailey, C. L., Barber, R. W., Emerson, D. R., Lockerby, D. A. & Reese, J. M. 2005 A critical review of the drag force on a sphere in the transition flow regime.  AIP Conference Proceedings 762, 743–748.
  * Beckmann et al. (2018) Beckmann, A. F., Rana, A. S., Torrilhon, M. & Struchtrup, H. 2018 Evaporation boundary conditions for the linear R13 equations based on the Onsager theory.  Entropy 20, 680.
  * Beresnev & Chernyak (1995) Beresnev, S. & Chernyak, V. 1995 Thermophoresis of a spherical particle in a rarefied gas: Numerical analysis based on the model kinetic equations.  Phys. Fluids 7, 1743–1756.
  * Bobylev (1982) Bobylev, A. V. 1982 The Chapman–Enskog and Grad methods for solving the Boltzmann equation. Sov. Phys. Dokl. 27, 29–31.
  * Bobylev (2018) Bobylev, A. V. 2018 Boltzmann equation and hydrodynamics beyond Navier–Stokes.  Philos. Trans. R. Soc. A 376, 20170227.
  * Cai et al. (2014) Cai, Z., Fan, Y. & Li, R. 2014 Globally hyperbolic regularization of Grad’s moment system.  Comm. Pure Appl. Math. 67, 464–518.
  * Cai & Wang (2020) Cai, Z. & Wang, Y. 2020 Regularized 13-moment equations for inverse power law models.  J. Fluid Mech. 894, A12.
  * Cercignani (1975) Cercignani, C. 1975 Theory and Application of the Boltzmann Equation. Edinburgh: Scottish Academic.
  * Cercignani (2000) Cercignani, C. 2000 Rarefied Gas Dynamics: From Basic Concepts to Actual Calculations. Cambridge: Cambridge University Press.
  * Chapman & Cowling (1970) Chapman, S. & Cowling, T. G. 1970 The Mathematical Theory of Non-Uniform Gases. Cambridge: Cambridge University Press.
  * Chernyak & Margilevskiy (1989) Chernyak, V. G. & Margilevskiy, A. Y. 1989 The kinetic theory of heat and mass transfer from a spherical particle in a rarefied gas.  Int. J. Heat Mass Transfer 32, 2127–2134.
  * Claydon et al. (2017) Claydon, R., Shrestha, A., Rana, A. S., Sprittles, J. E. & Lockerby, D. A. 2017 Fundamental solutions to the regularised 13-moment equations: efficient computation of three-dimensional kinetic effects.  J. Fluid Mech. 833, R4.
  * Goldberg (1954) Goldberg, R. 1954 Slow flow of a rarefied gas past a spherical obstacle. PhD thesis, New York University.
  * Grad (1949) Grad, H. 1949 On the kinetic theory of rarefied gases.  Comm. Pure Appl. Math. 2, 331–407.
  * de Groot & Mazur (1962) de Groot, S. R. & Mazur, P. 1962 Non-Equilibrium Thermodynamics. Amsterdam: North-Holland.
  * Gu et al. (2019) Gu, X.-J., Barber, R. W., John, B. & Emerson, D. R. 2019 Non-equilibrium effects on flow past a circular cylinder in the slip and early transition regime.  J. Fluid Mech. 860, 654–681.
  * Gu & Emerson (2007) Gu, X. J. & Emerson, D. R. 2007 A computational strategy for the regularized 13 moment equations with enhanced wall-boundary conditions.  J. Comput. Phys. 225, 263–283.
  * Gu & Emerson (2009) Gu, X.-J. & Emerson, D. R. 2009 A high-order moment approach for capturing non equilibrium phenomena in the transition regime.  J. Fluid Mech. 636, 177–216.
  * Gu & Emerson (2014) Gu, X.-J. & Emerson, D. R. 2014 Linearized-moment analysis of the temperature jump and temperature defect in the Knudsen layer of a rarefied gas.  Phys. Rev. E 89, 063020.
  * Gu et al. (2010) Gu, X.-J., Emerson, D. R. & Tang, G.-H. 2010 Analysis of the slip coefficient and defect velocity in the Knudsen layer of a rarefied gas using the linearized moment equations.  Phys. Rev. E 81, 016313.
  * Gupta (2015) Gupta, V. K. 2015 Mathematical modeling of rarefied gas mixtures. PhD thesis, RWTH Aachen University, Germany.
  * Gupta (2020) Gupta, V. K. 2020 Moment theories for a $d$-dimensional dilute granular gas of Maxwell molecules.  J. Fluid Mech. 888, A12.
  * Gupta & Torrilhon (2012) Gupta, V. K. & Torrilhon, M. 2012 Automated Boltzmann collision integrals for moment equations.  AIP Conf. Proc. 1501, 67–74.
  * Kalempa & Sharipov (2020) Kalempa, D. & Sharipov, F. 2020 Drag and thermophoresis on a sphere in a rarefied gas based on the Cercignani–Lampis model of gas–surface interaction.  J. Fluid Mech. 900, A37.
  * Karlin (2018) Karlin, I. 2018 Derivation of regularized Grad’s moment system from kinetic equations: modes, ghosts and non-Markov fluxes.  Philos. Trans. R. Soc. A 376, 20170230.