# Rarefied gas flow past a liquid droplet: Interplay between internal and
external flows

Rahul Bhattacharjee1  Vinay Kumar Gupta2  Anirudh S. Rana1
<EMAIL_ADDRESS>1 Department of Mathematics, Birla
Institute of Technology & Science Pilani, Rajasthan 333031, India 2 Discipline
of Mathematics, Indian Institute of Technology Indore, Indore 453552, India

###### Abstract

Experimental and theoretical studies on millimeter-sized droplets suggest that
at low Reynolds number the difference between the drag force on a circulating
water droplet and that on a rigid sphere is very small (less than 1%) (LeClair
et al., 1972), which is not the case when the droplet is of
micrometer/nanometer size. The goal of this article is to study the effects of
internal motion within a spherical micro/nano droplet—such that its diameter
is comparable to the mean free path of the surrounding gas—on the drag force
and its overall dynamics. To this end, the problem of a slow rarefied gas
flowing over an incompressible liquid droplet is investigated analytically by
considering the internal motion of the liquid inside the droplet and also by
accounting for kinetic effects in the gas. Detailed results for different
values of the Knudsen number, the ratio of the thermal conductivities and the
ratio of viscosities are presented for the pressure and temperature profiles
inside and outside the liquid droplet. The results for the drag force obtained
in the present work are in good agreement with the theoretical and
experimental results existing in the literature.

##  1 Introduction

Liquid droplets moving in a gaseous medium are frequently encountered in
nature; for instance, in rainfalls, in coughing or sneezing, in irrigation
mist, etc.; they also have tremendous industrial applications, such as in
aerosol spray, spray cooling, thermal spray coating, agricultural spraying,
spray painting, food processing, fuel injection, fuel combustion, etc.
Therefore, a clear understanding of the motion of liquid droplets in a gaseous
medium (or, in other words, gas flow over a liquid droplet) is crucial in
designing devices involving droplet motion in gaseous medium and/or in
improving their performance.

For a gas flow over a liquid droplet, the ratio of the viscosities of the
liquid inside the droplet to that of the surrounding fluid—referred to as the
inside-to-outside viscosity ratio or, simply, the viscosity ratio and denoted
by $\Lambda_{\mu}$—is a nonzero finite number. The two limiting cases of the
problem are (i) fluid flow over a solid sphere [case of infinite inside-to-
outside viscosity ratio ($\Lambda_{\mu}\to\infty$)] and (ii) liquid flow over
a gas bubble [case of zero inside-to-outside viscosity ratio
($\Lambda_{\mu}\approx 0$)]. In the former case, there is no question of
internal motion and in the latter case, the internal fluid motion has
negligible effect on the shape of the gas bubble as well as on the dynamics of
the external flow (Oliver & Chung, 1985, 1987; Pozrikidis, 1989). Thus, it is
not surprising that these two limiting cases have been explored extensively in
the literature [see the references given in chapters 3 and 5 of the textbook
Clift et al. (1978), which presents a comprehensive review of these two
limiting cases], and they are now seemingly well understood. The internal
fluid motion in the case of liquid droplets, however, has a significant impact
on the dynamics of the external flow and, hence, should not be disregarded
(Oliver & Chung, 1987; Pozrikidis, 1989). From a mathematical standpoint, the
coupling of the external flow in case of liquid droplet with the internal flow
is through complicated boundary conditions on the droplet interface that makes
the theoretical and computational methods of analysis considerably involved.
From an experimental point of view, it is very challenging to measure the flow
inside the droplet without disturbing the shape of the droplet or without
changing the physical properties of the fluids. Consequently, only a handful
of studies have looked into the problem of gas flow over a liquid droplet
considering the internal fluid motion so far. Indeed, the authors could not
find any paper on a rarefied gas flow over a _micro-/nano-size_ liquid droplet
that accounts for the internal flow dynamics, although the problem of rarefied
gas flow past a micro-/nano-size evaporating/non-evaporating droplet without
internal circulation has been a subject of some recent works (see, e.g., Rana
et al., 2019, 2021a, 2021b; Tiwari et al., 2021; De Fraja et al., 2022).

The existence of internal circulation in liquid droplets falling in air (for
which $\Lambda_{\mu}\approx 56$) was speculated already in the beginning of
the last century by Lenard (1904) and was confirmed later through wind tunnel
experiments by Garner & Lane (1959) for large drops and by Pruppacher & Beard
(1970) for small drops. To the best of authors’ knowledge, the first attempt
to explain the effects of almost all the factors, including internal
circulation, on the shape and dynamics of large raindrops falling in air was
made by McDonald (1954). Through this study, McDonald (1954) concluded that,
for large drops, the internal circulation plays only a negligible role in
controlling the shape of the drop. To investigate the internal circulation in
water drops falling at terminal velocity in air, LeClair et al. (1972)
proposed four theoretical approaches based on (i) creeping flow assumption for
both internal and external flows, (ii) the assumptions of irrotational
external flow and inviscid internal flow, (iii) the boundary layer theory and
(iv) solving the vorticity-streamfunction formalism of the Navier–Stokes
equations numerically for both internal and external flows together. In the
same paper, they also presented a wind tunnel experimental study, similarly to
that of Pruppacher & Beard (1970), to gauge the validity of the results from
their theoretical approaches. By comparing the results obtained from the
theoretical approaches with those obtained from the wind tunnel experiment,
LeClair et al. (1972) found that the first approach markedly underestimated
the internal velocity while the second approach markedly overestimated it and
that the results from the third and fourth approaches were in reasonably good
agreement with the experimental data for drops of diameters smaller than 1
millimeter. However, for large drops (of diameters bigger than 1 millimeter),
they found that even their third approach overestimated the internal velocity
significantly showing a completely wrong trend and that their numerical
approach, although overestimating the internal velocity slightly, was able to
capture the trend of the internal velocity qualitatively. Furthermore, LeClair
et al. (1972) also concluded that, for low Reynolds number $\Rey$, the drag
force on the drop is practically the same as the drag force on a solid sphere
of the same Reynolds number. Following the numerical approach of LeClair et
al. (1972), Abdel-Alim & Hamielec (1975) investigated the effect of internal
circulation on the drag on a spherical droplet falling at terminal velocity
and presented an empirical formula—obtained by fitting their numerical
results—for the drag coefficient as a function of the viscosity ratio and
external Reynolds number. In another similar numerical study, Rivkind et al.
(1976) also solved the vorticity-streamfunction formalism of the Navier–Stokes
equations numerically via the method of finite differences to determine the
drag on a spherical fluid drop falling in another fluid for viscosity ratios
$0\leq\Lambda_{\mu}<\infty$ and for external Reynolds numbers $0.5\leq\Rey\leq
100$ that cover flow over a solid sphere, over a liquid drop and over a small
gas bubble. For $\Rey\ll 1$, they found that the the drag coefficient of the
drop can be expressed as a convex combination of the drag coefficients of the
solid sphere and that of the gas bubble, with the coefficients in the
combination being functions of the viscosity ratio. Their formula for the drag
coefficient of the drop turned out to yield fairly accurate drag coefficient
for $\Rey\ll 1$. Rivkind & Ryskin (1976) furthered the study to moderate
Reynolds numbers and also gave another (empirical) formula for the drag
coefficient of the drop as a function of the the viscosity ratio and external
Reynolds number. However, a comparison of the drag coefficients obtained from
the formulas of Abdel-Alim & Hamielec (1975) and Rivkind & Ryskin (1976)
reveals that there could be differences up to 20% (for $\Rey\leq 20$) in the
values of the drag coefficients obtained from them. Moreover, the drag
coefficients computed from neither of the formulas of Abdel-Alim & Hamielec
(1975) and Rivkind & Ryskin (1976) could approach to the drag coefficient
obtained from the Hadamard and Rybczynski relation (Clift et al., 1978) in the
vanishing Reynolds number limit. Aiming to decipher the discrepancies arising
from the formulas of Abdel-Alim & Hamielec (1975) and Rivkind & Ryskin (1976),
Oliver and Chung performed two studies on flows inside and outside of a fluid
sphere—first one for low Reynolds numbers and the second one for moderate
Reynolds numbers. In their first study, Oliver & Chung (1985) employed a
hybrid semi-analytical method comprising of the series-truncation technique
and the finite-difference method to study the effect of internal circulation
on bubble and droplet dynamics at low Reynolds numbers. They found that the
density difference has no significant effect on the drag coefficient at low
Reynolds numbers and that the drag coefficient increases with increasing
viscosity ratio. In their second study, Oliver & Chung (1987) employed another
hybrid semi-analytical method comprising of the series-truncation technique
and the finite-element method to predict the flows inside and outside a fluid
droplet at low to moderate Reynolds numbers. In this study, they found that
the formula of the drag coefficient from Abdel-Alim & Hamielec (1975) is
actually dubious while that from Rivkind & Ryskin (1976) is good in predicting
the drag coefficient for $2\leq\Rey\leq 20$. Since the formulas of both Abdel-
Alim & Hamielec (1975) and Rivkind & Ryskin (1976) for the drag coefficient
are inadequate in the zero Reynolds number limit, Oliver & Chung (1987) also
gave a predictive formula for the drag coefficient valid for $0<\Rey<2$. They
also concluded that the strength of the internal circulation increases with
increasing Reynolds number.

As aforementioned, we have found neither any theoretical work nor any
experimental work on rarefied gas flow around a micro-/nano-size liquid
droplet—especially when accounting for the internal circulation—in the
literature. Given that setting up an experiment at such a small scale is even
more challenging, the objective of this paper is to investigate the aforesaid
problem theoretically. It is important to note that the Navier–Stokes–Fourier
(NSF) equations are not adequate for rarefied gas flows, which demands to
employ a better model to the external flow (Struchtrup, 2005). Although the
Boltzmann equation is capable of solving rarefied gas flows very well, its
numerical solutions are in general computationally expensive. Therefore, we
shall model the gas phase with the system of the regularised 26-moment (R26)
equations (Gu & Emerson, 2009), which is the best known macroscopic model to
describe rarefied gas flows and is an extension of the system of the
regularised 13-moment (R13) equations derived from the Boltzmann equation with
two different approaches detailed in Struchtrup & Torrilhon (2003) and
Struchtrup (2004). It is, however, sufficient to model the liquid phase inside
the droplet with the Navier–Stokes equations, although we assume that the
surface tension forces on the droplet are strong enough to maintain it
spherical.

To find appropriate boundary conditions concomitant to the R13 and R26
equations is another challenging task; nevertheless, a remarkable progress has
been made in this direction since the pioneering work of Gu & Emerson (2007)
on deriving the boundary conditions for the R13 equations. Torrilhon &
Struchtrup (2008) noticed some inconsistencies in the boundary conditions
derived by Gu & Emerson (2007) and presented improved boundary conditions for
the R13 equations based on physical and mathematical requirements for the
problem under consideration. The boundary conditions of Torrilhon & Struchtrup
(2008) may generically be referred to as the macroscopic boundary conditions
(MBC). Following the approach of Torrilhon & Struchtrup (2008), Gu & Emerson
(2009) derived the MBC for the R26 equations. Recently, the MBC for the R13
equations have also been combined with the discrete velocity method in a
hybrid approach by Yang et al. (2020) to make the computations faster in the
near-wall region. Notwithstanding, Rana & Struchtrup (2016) and Rana et al.
(2021a) showed that the MBC for the linearised R13 (LR13) equations as well as
for the linearised R26 (LR26) equations are thermodynamically inconsistent and
violate the Onsager reciprocity relations (Onsager, 1931a, b) for some
boundary value problems, and proposed a new set of phenomenological boundary
conditions (PBC) for the LR13 equations in Rana & Struchtrup (2016) and that
for the LR26 equations in Rana et al. (2021a). As a next step, the PBC valid
for for processes involving phase change were also derived by Beckmann et al.
(2018) for the R13 equations and by Rana et al. (2021a) for the R26 equations.
In this paper, we shall employ the PBC derived Rana et al. (2021a) for the
external flow.

The remainder of this paper is organised as follows. The governing equations
in spherical coordinates along with the boundary conditions are presented in §
2. The methodology for solving the problem analytically is outlined in § 3.
The main results on the effect of the internal flow inside the liquid droplet
on the motion of the rarefied gas flow are illustrated in § 4. Finally,
concluding remarks are made in § 5.

##  2  Problem formulation

We consider a slow steady uniform flow of a monatomic rarefied gas approaching
from the negative $z$-direction over a spherical droplet made of an
incompressible liquid and centered at origin, as depicted in Fig. 1. Since
incompressible liquid flows can be described accurately by the most celebrated
equations of fluid dynamics—the Navier–Stokes equations, we model the flow
inside the liquid droplet with the Navier–Stokes equations. However, as stated
in Sec. 1, the Navier–Stokes equations are not adequate for describing
rarefied gas flows; therefore, we model the gas flow using the R26 equations,
which describe rarefied gas flows remarkably well. To exploit the spherical
symmetry of the droplet, we shall express all the equations in the spherical
coordinate system $(\hat{r},\theta,\phi)$, which is related to the Cartesian
coordinate system $(\hat{x},\hat{y},\hat{z})$ via
$(x,y,z)\equiv(\hat{r}\sin{\theta}\cos{\phi},\hat{r}\sin{\theta}\sin{\phi},\hat{r}\cos{\theta})$.
Here, $\hat{r}\in[0,\infty)$, $\theta\in[0,\pi]$ and $\phi\in[0,2\pi)$. The
spherical symmetry of the droplet implies that the flow parameters are
independent of the direction $\phi$. Consequently, all the field variables
pertaining to the problem are functions of $r$ and $\theta$ only.

Figure 1: Schematic of a rarefied gas flow past a spherical liquid droplet.

For the aforesaid problem, an analytic solution to the full Navier–Stokes
equations and the fully nonlinear R26 equations is seemingly impossible.
Therefore, we are aiming to obtain an analytic solution to this problem with
the linearised equations and linearised boundary conditions. Such an analytic
solution is valid for slow flows. Given that rarefied gas flows encountered in
micro and nano devices are usually slow flows, an analytic solution obtained
by solving the linearised equations along with the linearised boundary
conditions is very plausible for all practical purposes.

To obtain the linearised equations, the governing equations and boundary
conditions are linearised around a reference state, given by a constant
density $\hat{\rho}_{0}$, a constant temperature $\hat{T}_{0}$ and all other
field variables as zero. For simplicity, we shall work with the dimensionless
equations and boundary conditions, which are obtained by introducing the
dimensionless deviations from their reference state values. Here, the
deviations are assumed to be sufficiently small so that flow description with
the linearised equations remains valid.

###  2.1 Modeling of the gas phase

The gas phase in the problem is modeled with the linear, dimensionless,
steady-state R26 equations. The (fully nonlinear) R26 equations in the
Cartesian coordinate system have been propounded in Gu & Emerson (2009). The
field variables in the R26 equations are the density $\hat{\rho}$, velocity
$\hat{v}_{i}$, temperature $\hat{T}$, stress tensor $\hat{\sigma}_{ij}$, heat
flux $\hat{q}_{i}$, (tracefree) third velocity moment $\hat{m}_{ijk}$,
partially contracted (tracefree) fourth velocity moment $\hat{R}_{ij}$ and
fully contracted fourth velocity moment (also referred to as the scalar fourth
moment) $\hat{\Delta}$. The field variables with hat are the usual quantities
with dimensions, like the ones taken in Gu & Emerson (2009) but without hats.
The R26 equations are linearised around the reference state described above
and are made dimensionless by introducing the dimensionless deviations in the
field variables from their respective reference state values:

$\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}},\quad
p:=\frac{\hat{p}-\hat{p}_{0}}{\hat{p}_{0}},\quad\sigma_{ij}:=\frac{\hat{\sigma}_{ij}}{\hat{p}_{0}},\\\
q_{i}:=\frac{\hat{q}_{i}}{\hat{p}_{0}\sqrt{\hat{R}\hat{T}_{0}}},\quad
m_{ijk}:=\frac{\hat{m}_{ijk}}{\hat{p}_{0}\sqrt{\hat{R}\hat{T}_{0}}},\quad
R_{ij}:=\frac{\hat{R}_{ij}}{\hat{p}_{0}\hat{R}\hat{T}_{0}},\quad\Delta:=\frac{\hat{\Delta}}{\hat{p}_{0}\hat{R}\hat{T}_{0}},\end{gathered}\right\\}$
(3)

where $\hat{R}$ is the gas constant;
$\hat{p}_{0}=\hat{\rho}_{0}\hat{R}\hat{T}_{0}$ is the pressure in the
reference state; $\rho$, $v_{i}$, $T$, $p$, $\sigma_{ij}$ and $q_{i}$ are the
dimensionless deviations in the density, velocity, temperature, stress and
heat flux of the gas, respectively; similarly, $m_{ijk}$, $R_{ij}$ and
$\Delta$ are the dimensionless deviations in the corresponding quantities. In
addition, the droplet radius $\hat{R}_{0}$ is taken as the length scale for
making the space variable $\hat{r}$ dimensionless, i.e.
$r:=\hat{r}/\hat{R}_{0}$. We insert the dimensionless deviations (3) in the
original R26 equations and drop all the nonlinear terms in deviations along
with the time derivative terms. Finally, on transforming the resulting
equations from the Cartesian coordinate system to the spherical coordinate
system, we obtain the linear, dimensionless, steady-state R26 equations, which
read (Rana et al., 2021a)111There is a typo in the last term on the left-hand
side of equation (6.9$b$) of Rana et al. (2021a).

(4) () () (6) $\displaystyle\frac{\partial m_{rrr}}{\partial
r}+\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}$
$\displaystyle=-\frac{1}{\mathrm{Kn}}\sigma_{rr},$ (7a)
$\displaystyle\frac{\partial m_{rr\theta}}{\partial
r}+\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}$
$\displaystyle=-\frac{1}{\mathrm{Kn}}\sigma_{r\theta},$ (7b)
$\displaystyle\frac{1}{2}\left(\frac{\partial R_{rr}}{\partial
r}+\frac{3R_{rr}}{r}\right)+\frac{1}{2}\frac{\mathscr{D}R_{r\theta}}{r}+\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},$ (8a)
$\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},$ (8b)
$\displaystyle\frac{1}{r}\mathscr{D}\left(-\frac{6}{5}\sigma_{r\theta}+\Phi_{rrr\theta}-\frac{6}{35}R_{r\theta}\right)+\frac{9}{5}\left(\frac{\partial\sigma_{rr}}{\partial
r}-\frac{2\sigma_{rr}}{r}\right)$
$\displaystyle+\frac{\partial\Phi_{rrrr}}{\partial
r}+\frac{5\Phi_{rrrr}}{r}+\frac{9}{35}\left(\frac{\partial R_{rr}}{\partial
r}-\frac{2R_{rr}}{r}\right)$
$\displaystyle=-\frac{\mathrm{Pr}_{m}}{\mathrm{Kn}}m_{rrr},$ (9a)
$\displaystyle\frac{6}{5r}\frac{\partial\sigma_{rr}}{\partial\theta}+\frac{6}{35r}\frac{\partial
R_{rr}}{\partial\theta}-\frac{1}{2r}\frac{\partial\Phi_{rrrr}}{\partial\theta}+\frac{8}{5}\left(\frac{\partial\sigma_{r\theta}}{\partial
r}-\frac{2\sigma_{r\theta}}{r}\right)$
$\displaystyle+\frac{8}{35}\left(\frac{\partial R_{r\theta}}{\partial
r}-\frac{2R_{r\theta}}{r}\right)+\frac{\partial\Phi_{rrr\theta}}{\partial
r}+\frac{5\Phi_{rrr\theta}}{r}$
$\displaystyle=-\frac{\mathrm{Pr}_{m}}{\mathrm{Kn}}m_{rr\theta},$ (9b)
$\displaystyle\frac{1}{r}\mathscr{D}\left(2m_{rr\theta}-\frac{2}{15}\Omega_{\theta}+\psi_{rr\theta}-\frac{28}{15}q_{\theta}\right)+\frac{56}{15}\left(\frac{\partial
q_{r}}{\partial r}-\frac{q_{r}}{r}\right)\hskip 27.0301pt$
$\displaystyle+2\left(\frac{\partial m_{rrr}}{\partial
r}+\frac{4m_{rrr}}{r}\right)+\frac{\partial\psi_{rrr}}{\partial
r}+\frac{4\psi_{rrr}}{r}+\frac{4}{15}\left(\frac{\partial\Omega_{r}}{\partial
r}-\frac{\Omega_{r}}{r}\right)$
$\displaystyle=-\frac{\mathrm{Pr}_{R}}{\mathrm{Kn}}R_{rr},$ (10a)
$\displaystyle 2\left(\frac{\partial m_{rr\theta}}{\partial
r}+\frac{4m_{rr\theta}}{r}\right)+\frac{\partial\psi_{rr\theta}}{\partial
r}+\frac{4\psi_{rr\theta}}{r}+\frac{1}{5}\left(\frac{\partial\Omega_{\theta}}{\partial
r}-\frac{\Omega_{\theta}}{r}\right)\hskip 17.07164pt$
$\displaystyle+\frac{14}{5}\left(\frac{\partial q_{\theta}}{\partial
r}-\frac{q_{\theta}}{r}\right)-\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{1}{5r}\frac{\partial\Omega_{r}}{\partial\theta}$
$\displaystyle=-\frac{\mathrm{Pr}_{R}}{\mathrm{Kn}}R_{r\theta},$ (10b)
$\displaystyle\frac{1}{r}\mathscr{D}\left(8q_{\theta}+\Omega_{\theta}\right)+8\left(\frac{\partial
q_{r}}{\partial r}+\frac{2q_{r}}{r}\right)+\frac{\partial\Omega_{r}}{\partial
r}+\frac{2\Omega_{r}}{r}=-\frac{\mathrm{Pr}_{\Delta}}{\mathrm{Kn}}\Delta$ (11)

with the additional unknowns in (9)–(11) being

$\displaystyle\Phi_{rrrr}$
$\displaystyle=-4\frac{\mathrm{Kn}}{\mathrm{Pr}_{\Phi}}\left[\frac{4}{7}\left(\frac{\partial
m_{rrr}}{\partial
r}-\frac{3m_{rrr}}{r}\right)-\frac{3}{7}\frac{\mathscr{D}m_{rr\theta}}{r}\right],$
(12a) $\displaystyle\Phi_{rrr\theta}$
$\displaystyle=-4\frac{\mathrm{Kn}}{\mathrm{Pr}_{\Phi}}\left[\frac{15}{28}\left(\frac{\partial
m_{rr\theta}}{\partial
r}-\frac{3m_{rr\theta}}{r}\right)+\frac{5}{14r}\frac{\partial
m_{rrr}}{\partial\theta}\right],$ (12b) $\displaystyle\psi_{rrr}$
$\displaystyle=-\frac{27}{7}\frac{\mathrm{Kn}}{\mathrm{Pr}_{\psi}}\left[\frac{3}{5}\left(\frac{\partial
R_{rr}}{\partial
r}-\frac{2R_{rr}}{r}\right)-\frac{2}{5}\frac{\mathscr{D}R_{r\theta}}{r}\right],$
(13a) $\displaystyle\psi_{rr\theta}$
$\displaystyle=-\frac{27}{7}\frac{\mathrm{Kn}}{\mathrm{Pr}_{\psi}}\left[\frac{8}{15}\left(\frac{\partial
R_{r\theta}}{\partial
r}-\frac{2R_{r\theta}}{r}\right)+\frac{2}{5r}\frac{\partial
R_{rr}}{\partial\theta}\right],$ (13b) $\displaystyle\Omega_{r}$
$\displaystyle=-\frac{7}{3}\frac{\mathrm{Kn}}{\mathrm{Pr}_{\Omega}}\left[\frac{\partial\Delta}{\partial
r}+\frac{12}{7}\frac{\mathscr{D}R_{r\theta}}{r}+\frac{12}{7}\left(\frac{\partial
R_{rr}}{\partial r}+\frac{3R_{rr}}{r}\right)\right],$ (14a)
$\displaystyle\Omega_{\theta}$
$\displaystyle=-\frac{7}{3}\frac{\mathrm{Kn}}{\mathrm{Pr}_{\Omega}}\left[\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}\right].$ (14b)

Here $\mathscr{D}\equiv\cot{\theta}+\frac{\partial}{\partial\theta}$.
Equations (4)–(14), henceforth, will be referred to as the linearised R26
(LR26) equations. The coefficients $\mathrm{Kn}$, $\mathrm{Pr}$,
$\mathrm{Pr}_{m}$, $\mathrm{Pr}_{R}$, $\mathrm{Pr}_{\Delta}$,
$\mathrm{Pr}_{\Phi}$, $\mathrm{Pr}_{\psi}$, $\mathrm{Pr}_{\Omega}$ in the LR26
equations are the dimensionless numbers arising from the non-
dimensionalisation of the equations. In particular, the numbers

$\displaystyle\mathrm{Kn}=\frac{\hat{\mu}}{\hat{\rho}_{0}\sqrt{\hat{R}\hat{T}_{0}}\hat{R}_{0}}\quad\text{and}\quad\mathrm{Pr}=\frac{5}{2}\frac{\hat{\mu}}{\hat{\kappa}}\hat{R}$
(15)

are referred to as the Knudsen number and the Prandtl number, respectively,
with $\hat{\mu}$ being the viscosity of the gas and $\hat{\kappa}$ being the
thermal conductivity of the gas. It should be noted that, owing to
linearisation, the viscosity $\hat{\mu}$ and thermal conductivity
$\hat{\kappa}$ in (15) are the viscosity and thermal conductivity of the gas
at the reference state temperature $\hat{T}_{0}$; and hence both are constant.
The values of the numbers $\mathrm{Pr}$, $\mathrm{Pr}_{m}$, $\mathrm{Pr}_{R}$,
$\mathrm{Pr}_{\Delta}$, $\mathrm{Pr}_{\Phi}$, $\mathrm{Pr}_{\psi}$,
$\mathrm{Pr}_{\Omega}$ depend on the choice of the interaction potential
between two gas molecules. For the Maxwell interaction potential used in the
present work, the values of these numbers are $\mathrm{Pr}=2/3$,
$\mathrm{Pr}_{m}=3/2$, $\mathrm{Pr}_{R}=7/6$, $\mathrm{Pr}_{\Delta}=2/3$,
$\mathrm{Pr}_{\Phi}=2.097$, $\mathrm{Pr}_{\psi}=1.698$,
$\mathrm{Pr}_{\Omega}=1$ (Gu & Emerson, 2009). The subscripts $r$ and $\theta$
with the vectors/tensors in the LR26 equations denote their respective
components; for instance, $v_{r}$ is the $r$-component of the deviation in the
velocity vector and $\sigma_{r\theta}$ is the $r\theta$-component of the
deviation in the stress tensor. Equation (4) can be identified as the equation
of continuity for the gas, equations (() ‣ 2.1) and (() ‣ 2.1) as the momentum
balance equations in the $r$\- and $\theta$-directions, respectively, and
equation (6) as the energy balance equation, equations (7a) and (7) as the
balance equations for the $rr$\- and $r\theta$-components of the stress,
equations (8a) and (8b) as the heat flux balance equations in the $r$\- and
$\theta$-directions, and so on.

###  2.2 Modeling of the liquid phase inside the droplet

The liquid phase inside the spherical droplet is modeled with the linearised,
steady-state, incompressible NSF equations. The complete (fully nonlinear and
unsteady) NSF equations in the spherical coordinate system can be found in
some standard textbooks on fluid dynamics; see, e.g., the textbook Batchelor
(1967). The linear, dimensionless, steady-state NSF equations are obtained by
dropping the time derivative terms in the full NSF equations, linearising the
field variables around the reference state defined above and making them
dimensionless using the density $\hat{\rho}_{0}$ and temperature $\hat{T}_{0}$
in the reference state and the droplet radius $\hat{R}_{0}$ as the length
scale. After simplification, the linear, dimensionless steady-state NSF
equations for modelling the liquid phase of the problem under consideration
read

(16) (18)

with

$\displaystyle\sigma_{rr}^{(\ell)}=-2\mathrm{Kn}^{(\ell)}\frac{\partial
v_{r}^{(\ell)}}{\partial
r},\quad\sigma_{r\theta}^{(\ell)}=-\mathrm{Kn}^{(\ell)}\left(\frac{\partial
v_{\theta}^{(\ell)}}{\partial
r}-\frac{v_{\theta}^{(\ell)}}{r}+\frac{1}{r}\frac{\partial
v_{r}^{(\ell)}}{\partial\theta}\right),$ (19) $\displaystyle
q_{r}^{(\ell)}=-\frac{5}{2}\Lambda_{\kappa}\frac{\mathrm{Kn}}{\mathrm{Pr}}\frac{\partial
T^{(\ell)}}{\partial r}\quad\text{and}\quad
q_{\theta}^{(\ell)}=-\frac{5}{2}\Lambda_{\kappa}\frac{\mathrm{Kn}}{\mathrm{Pr}}\frac{1}{r}\frac{\partial
T^{(\ell)}}{\partial\theta}.$ (20)

The superscript ‘$(\ell)$’ in (16)–(20) has been used to indicate that the
variables with the superscript ‘$(\ell)$’ belong to the liquid phase (i.e. to
the liquid droplet). The variables in (16)–(20) are as follows.

$\displaystyle p^{(\ell)}=\frac{\hat{p}^{(\ell)}-\hat{p}_{0}}{\hat{p}_{0}}$
(21)

is the dimensionless deviation in the pressure of the liquid droplet from the
pressure in the reference state $\hat{p}_{0}$;

$\displaystyle
v_{r}^{(\ell)}=\frac{\hat{v}_{r}^{(\ell)}}{\sqrt{\hat{R}\hat{T}_{0}}}\quad\text{and}\quad
v_{\theta}^{(\ell)}=\frac{\hat{v}_{\theta}^{(\ell)}}{\sqrt{\hat{R}\hat{T}_{0}}}$
(22)

are the dimensionless deviations in the $r$\- and $\theta$-components of the
velocity of the liquid, respectively;

$\displaystyle\sigma_{rr}^{(\ell)}=\frac{\hat{\sigma}_{rr}^{(\ell)}}{\hat{p}_{0}},\quad\text{and}\quad\sigma_{r\theta}^{(\ell)}=\frac{\hat{\sigma}_{r\theta}^{(\ell)}}{\hat{p}_{0}}$
(23)

are the dimensionless deviations in the $rr$\- and $r\theta$-components of the
stress tensor for the liquid, respectively;

$\displaystyle
q_{r}^{(\ell)}=\frac{\hat{q}_{r}^{(\ell)}}{\hat{p}_{0}\sqrt{\hat{R}\hat{T}_{0}}}\quad\text{and}\quad
q_{\theta}^{(\ell)}=\frac{\hat{q}_{\theta}^{(\ell)}}{\hat{p}_{0}\sqrt{\hat{R}\hat{T}_{0}}}$
(24)

are the dimensionless deviations in the $r$\- and $\theta$-components of the
heat flux of the liquid, respectively;

$\displaystyle\mathrm{Kn}^{(\ell)}=\Lambda_{\mu}\mathrm{Kn}$ (25)

is the Knudsen number for the liquid phase;
$\Lambda_{\mu}=\hat{\mu}^{(\ell)}/\hat{\mu}$ is the ratio of the viscosity of
the liquid to the viscosity of the gas; and
$\Lambda_{\kappa}=\hat{\kappa}^{(\ell)}/\hat{\kappa}$ in (20) is the ratio of
the thermal conductivity of the liquid to the thermal conductivity of the gas.
It is worthwhile noting that, owing to linearisation,
$p^{(\ell)}\approx\rho^{(\ell)}+T^{(\ell)}$ with $\rho^{(\ell)}$ and
$T^{(\ell)}$ being the dimensionless deviations in the density and temperature
of the liquid. For the same reason, $\hat{\mu}$, $\hat{\mu}^{(\ell)}$,
$\hat{\kappa}$, $\hat{\kappa}^{(\ell)}$ are the viscosities and thermal
conductivities at the reference temperature $\hat{T}_{0}$; consequently,
$\hat{\mu}$, $\hat{\mu}^{(\ell)}$, $\hat{\kappa}$, $\hat{\kappa}^{(\ell)}$,
$\Lambda_{\mu}$ and $\Lambda_{\kappa}$ are all constant in the present work.
The field variables with hat and superscript ‘${(\ell)}$’ can be identified as
the the field variables of the original NSF equations. Furthermore, equation
(16) can be identified as the equation of continuity, equations (2.2) and
(2.2) as the momentum balance equations in the $r$\- and $\theta$-directions,
respectively, and equation (18) as the energy balance equation.

###  2.3 Boundary conditions

The physically admissible boundary conditions, which respect the second law of
thermodynamics and satisfy the Onsager reciprocity relations, for the LR26
equations have been derived in Rana et al. (2021a). It may be noted that the
boundary conditions derived in Rana et al. (2021a) are general—they consider
the motion of both gas and boundary in the normal direction along with
evaporation. In the problem under consideration, evaporation has been ignored
for simplicity, the interface between the liquid and gas has been assumed to
be fixed, and that neither the liquid nor the gas can penetrate the interface
is assumed. Therefore, we would take the evaporation/condensation coefficient
$\vartheta=0$ and the normal component of the flow velocity relative to the
interface velocity $\mathcal{V}_{n}=0$ in the boundary conditions given in
Rana et al. (2021a). We skip more details of the boundary conditions for the
sake of succinctness and present them here directly; interested readers are
referred to Rana et al. (2021a) for more details. For the problem under
consideration, the $r$-direction is the normal direction while the $\theta$\-
and $\phi$-directions are the two tangential directions; nevertheless, the
present problem is independent of $\phi$ due to spherical symmetry. Therefore,
we shall replace the subscripts $n$ with $r$ and $i$ with $\theta$ in the
boundary conditions derived in Rana et al. (2021a). Consequently, the boundary
conditions complementing the LR26 equations (4)–(14)—for the problem under
consideration—at the interface (i.e. at $r=1$) read (Rana et al., 2021a)

$\displaystyle v_{r}$ $\displaystyle=v_{r}^{I}=0,$ (26a) $\displaystyle q_{r}$
$\displaystyle=-\frac{\chi}{2-\chi}\sqrt{\frac{2}{\pi}}\left(2\mathcal{T}+\frac{1}{2}\sigma_{rr}+\frac{5}{28}R_{rr}+\frac{1}{15}\Delta-\frac{1}{6}\Phi_{rrrr}\right),$
(26b) $\displaystyle m_{rrr}$
$\displaystyle=\frac{\chi}{2-\chi}\sqrt{\frac{2}{\pi}}\left(\frac{2}{5}\mathcal{T}-\frac{7}{5}\sigma_{rr}-\frac{1}{14}R_{rr}+\frac{1}{75}\Delta-\frac{13}{15}\Phi_{rrrr}\right),$
(26c) $\displaystyle\Psi_{rrr}$
$\displaystyle=\frac{\chi}{2-\chi}\sqrt{\frac{2}{\pi}}\left(\frac{6}{5}\mathcal{T}+\frac{9}{5}\sigma_{rr}-\frac{93}{70}R_{rr}+\frac{1}{5}\Delta+\frac{11}{15}\Phi_{rrrr}\right),$
(26d) $\displaystyle\Omega_{r}$
$\displaystyle=\frac{\chi}{2-\chi}\sqrt{\frac{2}{\pi}}\left(8\mathcal{T}+2\sigma_{rr}-R_{rr}-\frac{4}{3}\Delta-\frac{2}{3}\Phi_{rrrr}\right),$
(26e) $\displaystyle\sigma_{r\theta}$
$\displaystyle=-\frac{\chi}{2-\chi}\sqrt{\frac{2}{\pi}}\left(\mathscr{V}_{\theta}+\frac{1}{5}q_{\theta}+\frac{1}{2}m_{rr\theta}-\frac{1}{14}\Psi_{rr\theta}-\frac{1}{70}\Omega_{r}\right),$
(26f) $\displaystyle R_{r\theta}$
$\displaystyle=-\frac{\chi}{2-\chi}\sqrt{\frac{2}{\pi}}\left(-\mathscr{V}_{\theta}+\frac{11}{5}q_{\theta}+\frac{1}{2}m_{rr\theta}+\frac{13}{14}\Psi_{rr\theta}+\frac{13}{70}\Omega_{r}\right),$
(26g) $\displaystyle\Phi_{rrr\theta}$
$\displaystyle=-\frac{\chi}{2-\chi}\sqrt{\frac{2}{\pi}}\left(-\frac{4}{7}\mathscr{V}_{\theta}-\frac{12}{35}q_{\theta}+\frac{9}{7}m_{rr\theta}-\frac{2}{49}\Psi_{rr\theta}-\frac{2}{245}\Omega_{r}\right),$
(26h)

where $v_{r}^{I}$ is the velocity of the interface and is zero for the present
problem, $\chi$ is the accommodation coefficient,
$\mathscr{V}_{\theta}=v_{\theta}-v_{\theta}^{(\ell)}$ is the velocity slip and
$\mathcal{T}=T-T^{(\ell)}$ is the temperature jump. Note that the remaining
boundary conditions for the rank-2 and rank-3 tensors given in Rana et al.
(2021a) are not needed due to the fact that the present problem is independent
of $\phi$.

In order to solve the equations corresponding to the liquid and gas phases
together, we need additional boundary conditions, which are as follows. $(i)$
Similarly to the gas, the liquid cannot penetrate the interface. This means
that the normal component of the velocity of the liquid should also vanish at
the interface, i.e.

$\displaystyle v_{r}^{(\ell)}=0\quad\text{at}\quad r=1.$ (27)

$(ii)$ The heat flux and shear stress are continuous at the interface; this
implies that

$\displaystyle
q_{r}=q_{r}^{(\ell)}\quad\text{and}\quad\sigma_{r\theta}=\sigma_{r\theta}^{(\ell)}\quad\text{at}\quad
r=1.$ (28)

##  3 Analytic solution methodology

We solve the equations derived in § 2 analytically by following a method
proposed in Torrilhon (2010). The key idea of this method is to convert the
system of partial differential equations to a system of ordinary differential
equations by presuming the dependence of the field variables on the azimuthal
angle $\theta$ through the sine and cosine functions alone. In this method,
the scalar variables and vectorial/tensorial components of a field variable
having an even number of $\theta$ indices are taken to be proportional to
$\cos{\theta}$, the vectorial/tensorial components having an odd number of
$\theta$ indices are taken to be proportional to $\sin{\theta}$ and the
proportionality constants are taken to be functions of $r$ alone (Torrilhon,
2010; Rana et al., 2021a). Using this idea, we assume that the field variables
are given by

$\displaystyle\left.\begin{aligned}
v_{r}(r,\theta)&=\mathbbm{v}_{1}(r)\cos{\theta},&\quad
v_{\theta}(r,\theta)&=-\mathbbm{v}_{2}(r)\sin{\theta}\\\
p(r,\theta)&=\mathbbm{p}(r)\cos{\theta},&\quad
T(r,\theta)&=\mathbbm{T}(r)\cos{\theta},\\\
\sigma_{rr}(r,\theta)&=\mathbbm{s}_{1}(r)\cos{\theta},&\quad\sigma_{r\theta}(r,\theta)&=\mathbbm{s}_{2}(r)\sin{\theta},\\\
q_{r}(r,\theta)&=\mathbbm{q}_{1}(r)\cos{\theta},&\quad
q_{\theta}(r,\theta)&=-\mathbbm{q}_{2}(r)\sin{\theta},\\\
m_{rrr}(r,\theta)&=\mathbbm{m}_{1}(r)\cos{\theta},&\quad
m_{rr\theta}(r,\theta)&=\mathbbm{m}_{2}(r)\sin{\theta},\\\
R_{rr}(r,\theta)&=\mathbbm{R}_{1}(r)\cos{\theta},&\quad
R_{r\theta}(r,\theta)&=\mathbbm{R}_{2}(r)\sin{\theta},\\\
\Delta(r,\theta)&=\mathbbm{d}(r)\cos{\theta},\end{aligned}\right\\}$ (29)
$\displaystyle\left.\begin{aligned}
v_{r}^{(\ell)}(r,\theta)&=\mathbbm{v}_{1}^{(\ell)}(r)\cos{\theta},&\quad
v_{\theta}^{(\ell)}(r,\theta)&=-\mathbbm{v}_{2}^{(\ell)}(r)\sin{\theta},\\\
p^{(\ell)}(r,\theta)&=\mathbbm{p}^{(\ell)}(r)\cos{\theta},&\quad
T^{(\ell)}(r,\theta)&=\mathbbm{T}^{(\ell)}(r)\cos{\theta},\end{aligned}\right\\}$
(30)

where the functions
$\mathbbm{v}_{1},\mathbbm{v}_{2},\mathbbm{p},\mathbbm{T},\mathbbm{s}_{1},\mathbbm{s}_{2},\mathbbm{q}_{1},\mathbbm{q}_{2},\mathbbm{m}_{1},\mathbbm{m}_{2},\mathbbm{R}_{1},\mathbbm{R}_{2},\mathbbm{d},\mathbbm{v}_{1}^{(\ell)},\mathbbm{v}_{2}^{(\ell)},\mathbbm{p}^{(\ell)}$
and $\mathbbm{T}^{(\ell)}$ are the functions of $r$ alone.

The above ansatz for the field variables are inserted in (4)–(14) and in
(16)–(20). After simplification ($\cos{\theta}$ and $\sin{\theta}$ in each
equation get canceled), one obtains two systems of ordinary differential
equations—one for the liquid and the other for the gas. These systems of
ordinary differential equations are solved independently and analytically. The
analytic solution for the system associated with the liquid phase is easy to
obtain. It turns out to be

$\displaystyle\mathbbm{v}_{1}^{(\ell)}(r)=b_{1}+\frac{b_{2}r^{2}}{2},\quad\mathbbm{v}_{2}^{(\ell)}(r)=b_{1}+b_{2}r^{2},\quad\mathbbm{p}^{(\ell)}(r)=5b_{2}\Lambda_{\mu}\mathrm{Kn}\,r,\quad\mathbbm{T}^{(\ell)}(r)=b_{3}r,$
(31)

where $b_{1}$, $b_{2}$ and $b_{3}$ are the integration constants, which are
computed using the interface conditions (27) and (28). While obtaining
solution (31), we have also used the fact that the solution should remain
bounded as $r$ approaches zero. The analytic solution for the system
associated with the gas phase is, however, not so straightforward to obtain.
We have used the computer algebra software Mathematica® to obtain the analytic
solution for the system associated with the gas phase. For better readability,
the solution has, however, been relegated to appendix A. It may be noted that
this solution contains eight integration constants, namely $C_{1}$, $C_{2}$,
$C_{3}$ and $K_{1},K_{2},\dots,K_{5}$, which are computed using the eight
boundary conditions (26). After applying the boundary conditions, the solution
for all field variables for both the liquid and gas phases become known.

##  4 Results and discussion

Figures 2, 3 and 4 illustrate the heat flux lines and temperature contours for
the viscosity ratio $\Lambda_{\mu}=100$ and for the thermal conductivity
ratios $\Lambda_{\kappa}=1,10$ and $100$, respectively. The panels in top,
middle and bottom rows of each figure denote the results for the Knudsen
numbers $\mathrm{Kn}=0.09,0.36$ and $0.9$, respectively. The panels in the
left columns of figures 2, 3 and 4 illustrate the results obtained with the
R26 equations for the rarefied gas flow outside the droplet and with the NSF
equations for the liquid inside the droplet while those on the right columns
of figures 2, 3 and 4 are obtained with the NSF equations for both the liquid
inside the droplet and gas flow outside the droplet. Comparing the respective
panels on the left and right columns of figures 2, 3 and 4, the LR26 equations
show that the heat in the external flow transfers from the cold region (back
side of the liquid droplet) to the hot region (front side of the liquid
droplet), which is a non-Fourier effect and for the liquid droplet inside as
predicted by the NSF equations it is just the other way round, i.e. the heat
flows from the hot region to cold region inside the liquid droplet. On the
other hand, panels in the right column of figures 2, 3 and 4 show that the NSF
equations cannot capture the non-Fourier heat transfer in the external
(rarefied gas) flow. Furthermore, comparing figures 2, 3 and 4 with those
obtained by Rana et al. (2021a), the difference in the values of temperature
profiles is evident, which is caused due to the effect of the internal
circulation inside the liquid droplet.

Figure 2: Heat flux lines plotted on top of the temperature contours for the
viscosity ratio $\Lambda_{\mu}=100$ and for different values of the Knudsen
number: $\mathrm{Kn}=0.09$ (top row), $\mathrm{Kn}=0.36$ (middle row) and
$\mathrm{Kn}=0.9$ (bottom row). The results for the liquid phase (internal
flow) have been computed with the NSF equations in all the cases while those
for the gas phase (external flow) have been computed with the LR26 equations
for the panels in the left column and with the NSF equations for the panels in
the right column. The thermal conductivity ratio is $\Lambda_{\kappa}=1$.

Figure 3: Same as figure 2 but for the thermal conductivity ratio
$\Lambda_{\kappa}=10$.

Figure 4: Same as figure 2 but for the thermal conductivity ratio
$\Lambda_{\kappa}=100$.

Figures 5, 6 and 7 illustrate the velocity streamlines plotted over the
pressure contours for the thermal conductivity ratio $\Lambda_{\kappa}=100$
and for the viscosity ratios $\Lambda_{\mu}=1,10$ and $100$, respectively. The
panels in top, middle and bottom rows of each figure again denote the results
for the Knudsen numbers $\mathrm{Kn}=0.09,0.36$ and $0.9$, respectively. The
panels in the left columns of figures 5, 6 and 7 illustrate the results
obtained with the R26 equations for the rarefied gas flow outside the droplet
and with the NSF equations for the liquid inside the droplet while those on
the right columns of figures 5, 6 and 7 are obtained with the NSF equations
for both the liquid inside the droplet and gas flow outside the droplet. In
each of figures 5, 6 and 7, we observe that with the increase in the Knudsen
number the magnitude of the maximum value of the pressure contour increases
for both the gas and liquid droplet. Further comparing the corresponding
panels on the left and right columns of figures 5, 6 and 7, we observe that
there is no significant difference in the magnitude of the pressure contours
on either side, in the case of the liquid droplet. On the other hand, we can
clearly see from figures 5, 6 and 7 that the difference in the magnitude of
the pressure contours for the gas computed with the NSF equations is about
half of the magnitude of those computed with the LR26 equations for
$\mathrm{Kn}=0.09$. Moreover, as the Knudsen number increases, the difference
in the magnitude of the pressure contours computed from both the LR26 and NSF
equations is clearly noticeable.

Figure 5:  Velocity streamlines plotted over the pressure contours for the
thermal conductivity ratio $\Lambda_{\kappa}=100$ and for different values of
the Knudsen number: $\mathrm{Kn}=0.09$ (top row), $\mathrm{Kn}=0.36$ (middle
row) and $\mathrm{Kn}=0.9$ (bottom row). The results for the liquid phase
(internal flow) have been computed with the NSF equations in all the cases
while those for the gas phase (external flow) have been computed with the LR26
equations for the panels in the left column and with the NSF equations for the
panels in the right column. The viscosity ratio is $\Lambda_{\mu}=1$.

Figure 6: Same as figure 5 but for the viscosity ratio $\Lambda_{\mu}=10$.

Figure 7: Same as figure 5 but for the viscosity ratio $\Lambda_{\mu}=100$.

Figure 8 illustrates the variation of the $z$-component, velocity at the
interface for different values of the Knudsen number. The interface (at $r=1$)
between the liquid and gas has been demarcated by a red line in the figure. It
is evident from the figure that there is a discontinuity at the interface of
the liquid droplet, and the magnitude of the discontinuity increases with
increasing value of the Knudsen number.

Figure 8:  Slip velocity measured along the radial direction for Knudsen
numbers $\mathrm{Kn}=0.1,0.4$ and $1$, respectively.

Figure 9 exhibits the drag force normalised with the Stokes drag,
$F_{\mathrm{Stokes}}=6\pi\mathrm{Kn}u_{\infty}$, plotted over the Knudsen
number. The figure illustrates the drag force estimated by the LR26 equations
for a fixed value of the thermal conductivity ratio $\Lambda_{\kappa}=100$ and
for different values of the viscosity ratio $\Lambda_{\mu}$. In addition, the
explicit values of the normalised drag force for the same value of the thermal
conductivity ratio $\Lambda_{\kappa}=100$ and for different values of the
Knudsen number and viscosity ratio are also given in table 1.

Figure 9:  Normalised drag force on a liquid droplet from the NSF and R26 theories as a function of the Knudsen number. The symbols denote the experimental data of Goldberg (1954). The results are normalised with the Stokes drag $F_{\mathrm{Stokes}}=6\pi\mathrm{Kn}u_{\infty}$. The thermal conductivity ratio is $\Lambda_{\kappa}=100$. Knudsen number |  $F/F_{\mathrm{Stokes}}$ ($\Lambda_{\kappa}=100$)
---|---
$\Lambda_{\mu}=1$ | $\Lambda_{\mu}=5$ | $\Lambda_{\mu}=10$ | $\Lambda_{\mu}=100$ | $\Lambda_{\mu}=1000$
0.01 | 0.826315 | 0.929546 | 0.952284 | 0.976313 | 0.978936
0.1 | 0.772654 | 0.843830 | 0.857370 | 0.870973 | 0.872416
0.5 | 0.544453 | 0.570468 | 0.574446 | 0.578199 | 0.578583
1.0 | 0.396739 | 0.408531 | 0.410211 | 0.411767 | 0.411926
5.0 | 0.154088 | 0.155314 | 0.155473 | 0.155617 | 0.155631
10.0 | 0.091075 | 0.091425 | 0.091470 | 0.091510 | 0.091514
Table 1: Drag force normalised with the Stokes drag,
$F_{\mathrm{Stokes}}=6\pi\mathrm{Kn}u_{\infty}$, for different values of the
Knudsen number and viscosity ratio with the thermal conductivity ratio being
fixed at 100.

It is clear from figure 9 and table 1 that as the Knudsen number increases,
the drag force decreases and approaches zero (in general, except for that from
the NSF equations with which the drag force seemingly attains a nonzero
positive value) as Knudsen number approaches infinity. On the other hand, it
can be seen from figure 9 and table 1 that with the increase in the viscosity
ratio $\Lambda_{\mu}$, the drag force approaches unity for smaller values of
the Knudsen number. The figure also shows that the drag force predicted by the
LR26 equations for higher values of the viscosity ratios is in an excellent
agreement with the experimental results of Goldberg (1954).

We also compare the drag force computed analytically with the LR26 equations
in the present work with the drag force computed with an explicit formula
given by Happel & Brenner (1965), which—in our notations—reads

$\displaystyle
F_{\text{HB}}=6\pi\mathrm{Kn}u_{0}\left(\frac{1+\frac{2}{3\Lambda_{\mu}}}{1+\frac{1}{\Lambda_{\mu}}}\right).$
(40)

Figure 10:  Drag force on the liquid droplet normalised with the drag force
given by (40). The thermal conductivity ratio is $\Lambda_{\kappa}=100$.

Figure 10 depicts the drag force computed with the LR26 equations and
normalised with the drag force given by (40) again for a fixed value of the
thermal conductivity ratio $\Lambda_{\kappa}=100$ and for different values of
the viscosity ratio $\Lambda_{\mu}$. In addition, the explicit values of the
normalised drag force for the same value of the thermal conductivity ratio
$\Lambda_{\kappa}=100$ and for different values of the Knudsen number and
viscosity ratio are also given in table 2.

Knudsen number |  $F/F_{\mathrm{HB}}$ ($\Lambda_{\kappa}=100$)
---|---
$\Lambda_{\mu}=1$ | $\Lambda_{\mu}=5$ | $\Lambda_{\mu}=10$ | $\Lambda_{\mu}=100$ | $\Lambda_{\mu}=1000$
0.01 | 0.991579 | 0.984226 | 0.982044 | 0.979546 | 0.979262
0.1 | 0.927186 | 0.893468 | 0.884164 | 0.873857 | 0.872707
0.5 | 0.653344 | 0.604026 | 0.592398 | 0.580114 | 0.578777
1.0 | 0.476088 | 0.432563 | 0.423030 | 0.413131 | 0.412063
5.0 | 0.184906 | 0.164450 | 0.160332 | 0.156133 | 0.155684
10.0 | 0.109291 | 0.096803 | 0.094328 | 0.091813 | 0.091545
Table 2: Drag force normalised normalised with the drag force given by (40)
for different values of the Knudsen number and viscosity ratio with the
thermal conductivity ratio being fixed at 100.

From figure 10 and table 2, the drag force normalised with the drag force
given by (40) for very small Knudsen numbers and corresponding values of the
viscosity ratio $\Lambda_{\mu}$ is very close to unity which is not the case
for the drag force normalised with the Stokes drag as can be observed from
Fig. 9. But gradually we can observe from figures 9 and 10 that with the
increase in the values of the Knudsen number both the results analytically are
in excellent agreement with each other.

##  5 Conclusion

In this article we analytically analyse the problem of a slow rarefied gas
flow over an incompressible liquid droplet with internal circulation. To the
best of authors’ knowledge this is the first time that we are working on a
multi-physics problem, i.e. a problem where we are studying the physics of
both the liquid and gas. The flow in the rarefied gas phase is studied by
solving the LR26 equations and the flow in the liquid phase is studied by
solving the incompressible Navier–Stokes equations. The analytic solutions
obtained play a vital role in understanding the kinetic effects in rarefied
gases and also helps in developing a deeper understanding on the effects, the
internal motion on a liquid droplet has on the flow pattern of the rarefied
gases flowing past it. Moreover, the analytic results obtained from the LR26
theory with the PBC for the drag coefficient agree closely with the limited
experimental measurements available in the literature for a wide range of
Knudsen numbers. Our future goal would be to consider the same multi-physics
set up and study the multi-physics both way coupling numerically using
different numerical schemes available.

##  Appendix A Analytic expressions for the unknowns in ansatz (29) and (30)

The analytic solution for the unknowns appearing in (29) and (30) are as
follows.

$\displaystyle\mathbbm{v}_{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-
K_{2}\Theta_{2}(r)\left(\frac{1.4893\mathrm{Kn}^{3}}{r^{3}}+\frac{1.88528\mathrm{Kn}^{2}}{r^{2}}\right),$
(49) $\displaystyle\mathbbm{v}_{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+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),$
(50) $\displaystyle\mathbbm{p}(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-
K_{4}\Theta_{4}(r)\left(\frac{0.254173\mathrm{Kn}^{2}}{r^{2}}+\frac{0.172163\mathrm{Kn}}{r}\right)$
$\displaystyle-
K_{5}\Theta_{5}(r)\left(\frac{0.0501203\mathrm{Kn}^{2}}{r^{2}}+\frac{0.0226838\mathrm{Kn}}{r}\right),$
(51) $\displaystyle\mathbbm{T}(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-
K_{4}\Theta_{4}(r)\left(\frac{0.0629047\mathrm{Kn}^{2}}{r^{2}}+\frac{0.0426083\mathrm{Kn}}{r}\right)$
$\displaystyle-
K_{5}\Theta_{5}(r)\left(\frac{0.443094\mathrm{Kn}^{2}}{r^{2}}+\frac{0.200538\mathrm{Kn}}{r}\right),$
(52) $\displaystyle\mathbbm{s}_{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+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+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+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),$
(53) $\displaystyle\mathbbm{s}_{2}(r)={}$
$\displaystyle\frac{\mathrm{Kn}^{4}}{r^{4}}\left(-\frac{5C_{1}}{\mathrm{Kn}}+\frac{C_{2}}{\mathrm{Kn}^{3}}+\frac{C_{3}}{5\mathrm{Kn}^{2}}\right)$
$\displaystyle+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+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+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),$
(54) $\displaystyle\mathbbm{q}_{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-
K_{2}\Theta_{2}(r)\left(\frac{0.80299\mathrm{Kn}^{3}}{r^{3}}+\frac{1.01649\mathrm{Kn}^{2}}{r^{2}}\right),$
(55) $\displaystyle\mathbbm{q}_{2}(r)={}$
$\displaystyle\frac{\mathrm{Kn}^{2}}{r^{3}}\left(\frac{3C_{1}}{4}-\frac{C_{3}}{12\mathrm{Kn}}\right)$
$\displaystyle-
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+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),$
(56) $\displaystyle\mathbbm{m}_{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+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+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+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+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+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),$
(57) $\displaystyle\mathbbm{m}_{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+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+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+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+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+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).$
(58) $\displaystyle\mathbbm{R}_{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+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-
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-
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+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+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),$
(59) $\displaystyle\mathbbm{R}_{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+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-
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-
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+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+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),$
(60) $\displaystyle\mathbbm{d}(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+K_{5}\Theta_{5}(r)\left(\frac{3.75265\mathrm{Kn}^{2}}{r^{2}}+\frac{1.6984\mathrm{Kn}}{r}\right).$
(61)

In (A)–(A), the coefficients

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

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

## References

  * Abdel-Alim & Hamielec (1975) Abdel-Alim, A. H. & Hamielec, A. E. 1975 A theoretical and experimental investigation of the effect of internal circulation on the drag of spherical droplets falling at terminal velocity in liquid media.  Ind. Eng. Chem. Fundam. 14, 308–312.
  * Bailey et al. (2004) Bailey, C., Barber, R. W., Emerson, D., Lockerby, D. A. & Reese, J. 2004 A critical review of the drag force on a sphere in the transition flow regime. In 24th International Symposium on Rarefied Gas Dynamics, Bari, Italy, pp. 743–748.
  * Batchelor (1967) Batchelor, G. K. 1967 An Introduction to Fluid Dynamics. Cambridge: Cambridge University Press.
  * 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.
  * Clift et al. (1978) Clift, R., Grace, J. R. & Weber, M. E. 1978 Bubbles, Drops, and Particles. New York: Dover.
  * De Fraja et al. (2022) De Fraja, T. C., Rana, A. S., Enright, R., Cooper, L. J., Lockerby, D. A. & Sprittles, J. E. 2022 Efficient moment method for modeling nanoporous evaporation.  Phys. Rev. Fluids 7, 024201.
  * Garner & Lane (1959) Garner, F. H. & Lane, J. J. 1959 Mass transfer to drops of liquid suspended in a gas stream. Part II: Experimental work and results. Trans. Inst. Chem. Eng. 37, 162–172.
  * Goldberg (1954) Goldberg, R. 1954 The slow flow of a rarefied gas past a spherical obstacle. PhD thesis.
  * 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.
  * Happel & Brenner (1965) Happel, J. & Brenner, H. 1965 Low Reynolds number hydrodynamics prentice-hall, , vol. 331.
  * LeClair et al. (1972) LeClair, B. P., Hamielec, A. E., Pruppacher, H. R. & Hall, W. D. 1972 A theoretical and experimental study of the internal circulation in water drops falling at terminal velocity in air.  J. Atmos. Sci. 29, 728–740.
  * Lenard (1904) Lenard, P. 1904 Über regen. Meteor. Z. 21, 249–260.
  * McDonald (1954) McDonald, J. E. 1954 The shape and aerodynamics of large raindrops.  J. Atmos. Sci. 11, 478–494.
  * Oliver & Chung (1985) Oliver, D. L. R. & Chung, J. N. 1985 Steady flows inside and around a fluid sphere at low Reynolds numbers.  J. Fluid Mech. 154, 215–230.
  * Oliver & Chung (1987) Oliver, D. L. R. & Chung, J. N. 1987 Flow about a fluid sphere at low to moderate Reynolds numbers.  J. Fluid Mech. 177, 1–18.
  * Onsager (1931a) Onsager, L. 1931a Reciprocal relations in irreversible processes. I.  Phys. Rev. 37, 405–426.
  * Onsager (1931b) Onsager, L. 1931b Reciprocal relations in irreversible processes. II.  Phys. Rev. 38, 2265–2279.
  * Pozrikidis (1989) Pozrikidis, C. 1989 Inviscid drops with internal circulation.  J. Fluid Mech. 209, 77–92.
  * Pruppacher & Beard (1970) Pruppacher, H. R. & Beard, K. V. 1970 A wind tunnel investigation of the internal circulation and shape of water drops falling at terminal velocity in air.  Q. J. R. Meteorol. Soc. 96, 247–256.
  * Rana et al. (2021a) Rana, A. S., Gupta, V. K., Sprittles, J. E. & Torrilhon, M. 2021a $H$-theorem and boundary conditions for the linear R26 equations: application to flow past an evaporating droplet.  J. Fluid Mech. 924, A16.
  * Rana et al. (2019) Rana, A. S., Lockerby, D. A. & Sprittles, J. E. 2019 Lifetime of a nanodroplet: Kinetic effects and regime transitions.  Phys. Rev. Lett. 123, 154501.
  * Rana et al. (2021b) Rana, A. S., Saini, S., Chakraborty, S., Lockerby, D. A. & Sprittles, J. E. 2021b Efficient simulation of non-classical liquid–vapour phase-transition flows: a method of fundamental solutions.  J. Fluid Mech. 919, A35.
  * Rana & Struchtrup (2016) Rana, A. S. & Struchtrup, H. 2016 Thermodynamically admissible boundary conditions for the regularized 13 moment equations.  Phys. Fluids 28, 027105.
  * Rivkind & Ryskin (1976) Rivkind, V. Ya. & Ryskin, G. M. 1976 Flow structure in motion of a spherical drop in a fluid medium at intermediate Reynolds numbers.  Fluid Dyn. 11, 5–12.
  * Rivkind et al. (1976) Rivkind, V. Ya., Ryskin, G. M. & Fishbein, G. A. 1976 Flow around a spherical drop at intermediate Reynolds numbers.  J. Appl. Math. Mech. 40, 687–691.
  * Struchtrup (2004) Struchtrup, H. 2004 Stable transport equations for rarefied gases at high orders in the Knudsen number.  Phys. Fluids 16, 3921–3934.
  * Struchtrup (2005) Struchtrup, H. 2005 Macroscopic Transport Equations for Rarefied Gas Flows. Berlin: Springer.
  * Struchtrup & Torrilhon (2003) Struchtrup, H. & Torrilhon, M. 2003 Regularization of Grad’s 13 moment equations: Derivation and linear analysis.  Phys. Fluids 15, 2668–2680.
  * Tiwari et al. (2021) Tiwari, S., Klar, A. & Russo, G. 2021 Modelling and simulations of moving droplet in a rarefied gas.  Int. J. Comput. Fluid Dyn. 35, 666–684.
  * Torrilhon (2010) Torrilhon, M. 2010 Slow gas microflow past a sphere: Analytical solution based on moment equations.  Phys. Fluids 22, 072001.
  * Torrilhon & Struchtrup (2008) Torrilhon, M. & Struchtrup, H. 2008 Boundary conditions for regularized 13-moment equations for micro-channel-flows.  J. Comput. Phys. 227, 1982–2011.
  * Yang et al. (2020) Yang, W., Gu, X.-J., Wu, L., Emerson, D. R., Zhang, Y. & Tang, S. 2020 A hybrid approach to couple the discrete velocity method and Method of Moments for rarefied gas flows.  J. Comput. Phys. 410, 109397.