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

Abstract

Experimental and theoretical studies on millimetre-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., J. Atmos. Sci., vol. 29, 1972, pp. 728–740). While the drag force on a spherical liquid droplet at high viscosity ratios (of the liquid to the gas), is approximately the same as that on a rigid sphere of the same size, the other quantities of interest (e.g. the temperature) in the case of a rarefied gas flow over a liquid droplet differ from the same quantities in the case of a rarefied gas flow over a rigid sphere. The goal of this article is to study the effects of internal motion within a spherical microdroplet/nanodroplet – 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 a 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 $\varLambda _\mu$ – is a non-zero 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 ($\varLambda _\mu \to \infty$)) and (ii) liquid flow over a gas bubble (case of zero inside-to-outside viscosity ratio ($\varLambda _\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, Grace & Weber 1978), which present 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 the case of a 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 microsized/nanosized liquid droplet that accounts for the internal flow dynamics, although the problem of rarefied gas flow past a microsized/nanosized evaporating/non-evaporating droplet without internal circulation has been a subject of some recent works (see, e.g. Rana, Lockerby & Sprittles 2018b, 2019; Rana et al. 2021a,b; Tiwari, Klar & Russo 2021; De Fraja et al. 2022). Some open-source software, like OpenFOAM, in combination with methods, such as the volume-of-fluid, level-set and direct simulation Monte Carlo (DSMC) methods, have also been utilised to investigate gas–liquid multiphase flow problems (Malekzadeh & Roohi 2015; Chakraborty 2019; Chakraborty et al. 2019). Nevertheless, the problems studied with this software are typically in a somewhat different direction than the one considered in this paper; for instance, the above references focus on bubble/droplet formation and its dynamics.

The existence of internal circulation in liquid droplets falling in air (for which $\varLambda _\mu \approx 56$) was already speculated 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 stream function 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 mm. However, for large drops (of diameters bigger than 1 mm), 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 small values of the Reynolds number $Re$, 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, Ryskin & Fishbein (1976) also solved the vorticity stream function 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 \varLambda _\mu < \infty$ and for external Reynolds numbers $0.5 \leq \textit {Re} \leq 100$ that cover flow over a solid sphere, over a liquid drop and over a small gas bubble. For $\textit {Re} \ll 1$, they found that 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 a fairly accurate drag coefficient for $\textit {Re} \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 viscosity ratio and external Reynolds number. However, a comparison of the drag coefficients obtained from the formulae of Abdel-Alim & Hamielec (1975) and Rivkind & Ryskin (1976) reveals that there could be differences up to 20 % (for $\textit {Re} \leq 20$) in the values of the drag coefficients obtained from them. Moreover, the drag coefficients computed from neither of the formulae of Abdel-Alim & Hamielec (1975) nor Rivkind & Ryskin (1976) could approach 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 formulae of Abdel-Alim & Hamielec (1975) and Rivkind & Ryskin (1976), Oliver & 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 semianalytical 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 semianalytical 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 \textit {Re} \leq 20$. Since the formulae 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 < \textit {Re} < 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 microsized/nanosized 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. The liquid phase inside the droplet can be modelled with the Navier–Stokes equations. However, it is important to note that the Navier–Stokes–Fourier (NSF) equations are not adequate for describing rarefied gas flow (Sone 2002; Struchtrup 2005) outside the droplet. Any fluid flow – including a rarefied gas flow – can, in principle, be described by the Boltzmann equation; nevertheless, its numerical solutions are computationally very expensive in general and particularly for flows in the so-called transition regime (Struchtrup 2005). The main source of problems in dealing with the Boltzmann equation is the Boltzmann collision operator appearing on the right-hand side of the Boltzmann equation. Thus there has been a significant amount of research in developing ways alternative to directly solving the Boltzmann equation for investigating rarefied gas flows. One of the most commonly used numerical techniques for investigating rarefied gas flows is the DSMC method, which is a probabilistic particle-based method developed by Bird (1994) to solve the Boltzmann equation numerically. Since its development, the DSMC method has been ameliorated and employed to investigate several canonical rarefied gas flow problems; see, e.g. Rana, Mohammadzadeh & Struchtrup (2015), Stefanov, Roohi & Shoja-Sani (2022), Taheri, Roohi & Stefanov (2022), Sadr & Hadjiconstantinou (2023) and references therein. Nevertheless, the DSMC method also demands a very high computational cost for processes in the transition regime. Aiming to substitute for the involved Boltzmann collision operator in the Boltzmann equation, some simplified models – generically referred to as kinetic models – have also been proposed. Some widely used kinetic models are the Bhatnagar–Gross–Krook model (Bhatnagar, Gross & Krook 1954), the ellipsoidal statistical Bhatnagar–Gross–Krook model (Holway 1966) and the Shakhov model (Shakhov 1968). These kinetic models have also been utilised with the DSMC method. However, each of these kinetic models has its own shortcomings/difficulties; the reader is referred to Struchtrup (2005) for details of these kinetic models. The widely accepted models for describing transition-regime flows are the extended macroscopic equations derived from the Boltzmann equation predominantly through two asymptotic-expansion based approaches, namely the Chapman–Enskog expansion method (Chapman & Cowling 1970) and the Grad moment method (Grad 1949). The models resulting from both methods again have their own merits and demerits. Let us skip the details of them for the sake of succinctness; the interested reader may refer to Struchtrup (2005) and Torrilhon (2016) for details. To circumvent the demerits associated with the above two methods, Struchtrup & Torrilhon regularised the equations resulting from the Grad moment method (referred to as the Grad moment equations) by performing a Chapman–Enskog-like expansion on the Grad moment equations and derived the so-called regularised 13-moment (R13) equations (Struchtrup & Torrilhon 2003; Struchtrup 2004). The R13 equations, since their derivation, have been remarkably successful in describing rarefied gas flows in the transition regime; see Torrilhon (2016) and reference therein. Since the R13 equations have been derived via an asymptotic expansion in the powers of a dimensionless parameter the Knudsen number, which is defined as the ratio of the mean free path of the gas to a characteristic length scale in the problem, it is not surprising that the R13 equations yield meaningful results mostly for small Knudsen numbers (i.e. for flows in the early transition regime). Aiming to cover more of the transition regime, Gu & Emerson (2009) derived the regularised 26-moment (R26) equations by extending the method proposed by Struchtrup & Torrilhon (2003). The R26 equations, in general, describe transition-regime flows better than the R13 equations, especially for relatively large Knudsen numbers. Notwithstanding, the R26 equations also have limitations due to their derivation also through an asymptotic expansion in powers of the Knudsen number. It can be stated empirically that the R26 equations yield very good results for transition-regime flows up to the Knudsen number close to unity (Gu & Emerson 2009; Rana et al. 2018b, 2021a) but may yield quantitatively different results for certain processes beyond the Knudsen number unity; see, e.g. Rana et al. (2018b, 2021a). Despite this, the system comprised of the R26 equations (Gu & Emerson 2009) is the best known macroscopic model up to date for describing transition-regime rarefied gas flows. Therefore, we shall model the gas phase (outside the droplet) with the system of the R26 equations (Gu & Emerson 2009) and the liquid phase inside the droplet with the Navier–Stokes equations. For comparison purpose, we shall also include the analytic solution obtained by solving the NSF equations for the gas phase. It is worthwhile noting that although the surface tension force is an important force that ought to be accounted for while investigating gas–liquid multiphase flows, considering the effect of the surface tension forces is beyond the scope of this paper and will be considered elsewhere in the future. Here, we shall assume that the surface tension forces on the droplet are strong enough to maintain its spherical shape. This assumption is justified at least for droplets made of some commonly used liquids, as discussed in § 4.4.

To find appropriate boundary conditions concomitant to the R13 and R26 equations is another challenging task; nevertheless, 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; Beckmann et al. 2018) for some boundary value problems, and proposed a new set of phenomenological boundary conditions (PBC) for the LR13 equations in Rana & Struchtrup (2016) and for the LR26 equations in Rana et al. (2021a). As a next step, the PBC valid 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 in Rana et al. (2021a) for the external flow. In summary, we solve the LR26 equations – and also the NSF equations for comparison purposes – for the gas phase (outside the droplet) along with the PBC and the linearised Navier–Stokes equations for the liquid phase (inside the droplet) along with the coupled boundary conditions to obtain the analytic solution for the flow fields. To validate the analytic solution obtained in the present work, we also compare the drag force on the liquid droplet computed analytically in the present work with that obtained from Millikan's famous oil-drop experiment. A comparison of the drag force obtained from the present theory with the results obtained from Millikan's oil-drop experiment reveals that the drag force on a spherical liquid droplet at high viscosity ratios is nearly the same as that on a rigid sphere of the same size. Hence, in many practical applications (wherein the viscosity ratio is usually large), for example a water droplet moving through air, the droplet can be treated as a rigid sphere for the drag force computations. However, we show that the internal motion of the liquid in the droplet does have effects on the other quantities of interest (e.g. the temperature).

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 $\hat {z}$-direction with a uniform velocity $\hat {u}_\infty$ over a spherical droplet made of an incompressible liquid and centred at origin, as depicted in figure 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 § 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 $(\hat {x},\hat {y},\hat {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,{\rm \pi} ]$ and $\phi \in [0,2{\rm \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 $\hat {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 restrict the present study to Stokes flows, i.e. to small Reynolds number flows ($\textit {Re} \ll 1$), and to slow flows, i.e. to small Mach number flows ($Ma \ll 1$), so that the linearised equations and linearised boundary conditions can be utilised in order to obtain an analytic solution of the problem. Such an analytic solution is valid for slow flows ($Ma \ll 1$) and in for low Reynolds numbers ($\textit {Re} \ll 1$). Given that rarefied gas flows encountered in microdevices and nanodevices 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. Modelling of the gas phase

The gas phase in the problem is modelled 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$, (trace-free) third velocity moment $\hat {m}_{ijk}$, partially contracted (trace-free) fourth velocity moment $\hat {R}_{ij}$ and fully contracted fourth velocity moment (also referred to as the scalar fourth moment) $\hat {\varDelta }$. The field variables with hats 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,

(2.1)\begin{equation} \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 \varDelta := \frac{\hat{\varDelta}}{\hat{p}_0 \hat{R}\hat{T}_0}, \end{gathered}\right\} \end{equation}

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, pressure, stress and heat flux of the gas, respectively; similarly, $m_{ijk}$, $R_{ij}$ and $\varDelta$ 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 (2.1) 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)

(2.2)\begin{gather} \frac{\partial v_r}{\partial r} +\frac{2 v_r}{r} +\frac{\mathscr{D} v_\theta}{r} = 0, \end{gather}

(2.3a)$$\begin{gather} \frac{\partial p}{\partial r} +\frac{\partial \sigma_{rr}}{\partial r} +\frac{3 \sigma_{rr}}{r} +\frac{\mathscr{D} \sigma_{r\theta}}{r} = 0, \end{gather}$$

(2.3b)$$\begin{gather}\frac{\partial \sigma_{r\theta}}{\partial r} +\frac{3 \sigma_{r\theta}}{r} -\frac{1}{2 r} \frac{\partial \sigma_{rr}}{\partial \theta} +\frac{1}{r} \frac{\partial p}{\partial \theta}= 0, \end{gather}$$

(2.4)\begin{gather} \frac{\partial q_r}{\partial r} +\frac{2 q_r}{r} +\frac{\mathscr{D} q_\theta}{r} = 0, \end{gather}

(2.5a)\begin{gather} \frac{\partial m_{rrr}}{\partial r} +\frac{4 m_{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}{Kn} \sigma_{rr}, \end{gather}

(2.5b)\begin{gather} \frac{\partial m_{rr\theta}}{\partial r} +\frac{4 m_{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}\nonumber\\ \quad\qquad -\,\frac{1}{2 r} \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}{Kn} \sigma_{r\theta}, \end{gather}

(2.6a)$$\begin{gather} \frac{1}{2} \left(\frac{\partial R_{rr}}{\partial r} +\frac{3 R_{rr}}{r} \right) +\frac{1}{2} \frac{\mathscr{D} R_{r\theta}}{r} +\frac{1}{6} \frac{\partial \varDelta}{\partial r} -\frac{\partial p}{\partial r} +\frac{5}{2} \frac{\partial T}{\partial r} ={-} \frac{Pr}{Kn} q_r, \end{gather}$$

(2.6b)$$\begin{gather}\frac{1}{2} \left(\frac{\partial R_{r\theta}}{\partial r} +\frac{3 R_{r\theta}}{r}\right) +\frac{1}{6r} \frac{\partial \varDelta}{\partial \theta} -\frac{1}{4 r} \frac{\partial R_{rr}}{\partial \theta} -\frac{1}{r}\frac{\partial p}{\partial \theta} +\frac{5}{2r} \frac{\partial T}{\partial \theta} ={-} \frac{Pr}{Kn} q_\theta, \end{gather}$$

(2.7a)\begin{gather} \frac{1}{r}\mathscr{D} \left(-\frac{6}{5}\sigma_{r\theta}+\varPhi_{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) \nonumber\\ \quad +\, \frac{\partial \varPhi_{rrrr}}{\partial r}+\frac{5 \varPhi_{rrrr}}{r} +\frac{9}{35}\left( \frac{\partial R_{rr}}{\partial r}-\frac{2 R_{rr}}{r}\right) ={-}\frac{{Pr}_{m}}{Kn}m_{rrr}, \end{gather}

(2.7b)\begin{gather} \frac{6}{5r}\frac{\partial \sigma_{rr}}{\partial \theta} +\frac{6}{35r} \frac{\partial R_{rr}}{\partial \theta }-\frac{1}{2r}\frac{\partial \varPhi_{rrrr}}{\partial \theta} +\frac{8}{5}\left( \frac{\partial \sigma_{r\theta}}{\partial r}-\frac{2 \sigma_{r\theta}}{r}\right)\nonumber\\ \quad+\,\frac{8}{35} \left( \frac{\partial R_{r\theta}}{\partial r}-\frac{2 R_{r\theta}}{r}\right) +\frac{\partial \varPhi_{rrr\theta}}{\partial r}+\frac{5 \varPhi_{rrr\theta}}{r} ={-}\frac{{Pr}_{m}}{Kn}m_{rr\theta}, \end{gather}

(2.8a)\begin{align} &\frac{1}{r}\mathscr{D} \left( 2m_{rr\theta} -\frac{2}{15}\varOmega_{\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)\nonumber\\ &\qquad+2\left(\frac{\partial m_{rrr}}{\partial r}+\frac{4 m_{rrr}}{r}\right) +\frac{\partial \psi_{rrr}}{\partial r}+\frac{4\psi_{rrr}}{r} +\frac{4}{15}\left( \frac{\partial \varOmega_{r}}{\partial r}-\frac{\varOmega_{r}}{r}\right)={-}\frac{{Pr}_{R}}{Kn}R_{rr}, \end{align}

(2.8b)\begin{align} &2\left(\frac{\partial m_{rr\theta}}{\partial r}+\frac{4 m_{rr\theta}}{r}\right) + \frac{\partial \psi_{rr\theta}}{\partial r}+\frac{4 \psi_{rr\theta}}{r} +\frac{1}{5}\left( \frac{\partial \varOmega _{\theta}}{\partial r}-\frac{\varOmega _{\theta}}{r} \right)\nonumber\\ &\qquad + \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 \varOmega_{r}}{\partial \theta }={-}\frac{{Pr}_{R}}{Kn}R_{r\theta}, \end{align}

(2.9)\begin{gather} \frac{1}{r}\mathscr{D} \left( 8 q_{\theta}+\varOmega _{\theta }\right) + 8 \left( \frac{\partial q_r}{\partial r}+\frac{2 q_r}{r}\right) + \frac{\partial \varOmega_{r}}{\partial r}+\frac{2\varOmega_{r}}{r} ={-}\frac{{Pr}_{\varDelta }}{Kn}\varDelta \end{gather}

with the additional unknowns in (2.7)–(2.9) being

(2.10a)$$\begin{gather} \varPhi_{rrrr} ={-} 4 \frac{Kn}{{Pr}_{\varPhi}} \left[\frac{4}{7}\left( \frac{\partial m_{rrr}}{\partial r}-\frac{3 m_{rrr}}{r}\right) -\frac{3}{7}\frac{\mathscr{D} m_{rr\theta}}{r}\right], \end{gather}$$

(2.10b)$$\begin{gather}\varPhi_{rrr\theta} ={-} 4 \frac{Kn}{{Pr}_{\varPhi}} \left[\frac{15}{28}\left( \frac{\partial m_{rr\theta}}{\partial r}-\frac{3 m_{rr\theta}}{r}\right) +\frac{5}{14 r} \frac{\partial m_{rrr}}{\partial \theta}\right], \end{gather}$$

(2.11a)$$\begin{gather} \psi_{rrr} ={-}\frac{27}{7} \frac{Kn}{{Pr}_{\psi}} \left[\frac{3}{5}\left( \frac{\partial R_{rr}}{\partial r}-\frac{2 R_{rr}}{r}\right) -\frac{2}{5} \frac{\mathscr{D} R_{r\theta}}{r}\right], \end{gather}$$

(2.11b)$$\begin{gather}\psi_{rr\theta} ={-}\frac{27}{7} \frac{Kn}{{Pr}_{\psi}} \left[\frac{8}{15}\left( \frac{\partial R_{r\theta}}{\partial r}-\frac{2 R_{r\theta}}{r}\right) + \frac{2}{5 r} \frac{\partial R_{rr}}{\partial \theta}\right], \end{gather}$$

(2.12a)$$\begin{gather} \varOmega_{r} ={-}\frac{7}{3} \frac{Kn}{{Pr}_{\varOmega}} \left[\frac{\partial \varDelta }{\partial r}+\frac{12}{7}\frac{\mathscr{D} R_{r\theta}}{r} + \frac{12}{7} \left(\frac{\partial R_{rr}}{\partial r}+\frac{3 R_{rr}}{r}\right) \right], \end{gather}$$

(2.12b)$$\begin{gather}\varOmega_{\theta} ={-}\frac{7}{3} \frac{Kn}{{Pr}_{\varOmega}} \left[\frac{1}{r}\frac{\partial \varDelta}{\partial \theta} + \frac{12}{7}\left( \frac{\partial R_{r\theta}}{\partial r}+\frac{3 R_{r\theta}}{r}\right) -\frac{6}{7 r}\frac{\partial R_{rr}}{\partial \theta}\right]. \end{gather}$$

Here $\mathscr {D}\equiv \cot \theta +{\partial }/{\partial \theta }$. Equations (2.2)–(2.12), henceforth, will be referred to as the linearised R26 (LR26) equations. The coefficients ${Kn}$, ${Pr}$, ${Pr}_m$, ${Pr}_R$, ${Pr}_\varDelta$, ${Pr}_{\varPhi }$, ${Pr}_{\psi }$, ${Pr}_{\varOmega }$ in the LR26 equations are the dimensionless numbers arising from the non-dimensionalisation of the equations. In particular, the numbers

(2.13a,b)\begin{equation} {Kn} = \frac{\hat{\mu}}{\hat{\rho}_0 \sqrt{\hat{R}\hat{T}_0} \hat{R}_0} \quad\text{and}\quad {Pr} = \frac{5}{2} \frac{\hat{\mu}}{\hat{\kappa}} \hat{R} \end{equation}

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 the linearisation, the viscosity $\hat {\mu }$ and thermal conductivity $\hat {\kappa }$ in (2.13a,b) are the viscosity and thermal conductivity of the gas at the reference state temperature $\hat {T}_0$; and hence both are constant. Let us denote the viscosity and thermal conductivity of the gas at the reference state temperature $\hat {T}_0$ by $\hat {\mu }_0$ and $\hat {\kappa }_0$, respectively. Thus, owing to the linearisation, $\hat {\mu } = \hat {\mu }_0$ and $\hat {\kappa } = \hat {\kappa }_0$ throughout this work. The values of the numbers ${Pr}$, ${Pr}_m$, ${Pr}_R$, ${Pr}_\varDelta$, ${Pr}_{\varPhi }$, ${Pr}_{\psi }$, ${Pr}_{\varOmega }$ 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 ${Pr} = 2/3$, ${Pr}_m = 3/2$, ${Pr}_R = 7/6$, ${Pr}_\varDelta = 2/3$, ${Pr}_{\varPhi } = 2.097$, ${Pr}_{\psi } = 1.698$, ${Pr}_{\varOmega } = 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 (2.2) can be identified as the equation of continuity for the gas, (2.3a) and (2.3b) as the momentum balance equations in the $r$- and $\theta$-directions, respectively, and (2.4) as the energy balance equation, (2.5a) and (2.5b) as the balance equations for the $rr$- and $r\theta$-components of the stress, (2.6a) and (2.6b) as the heat flux balance equations in the $r$- and $\theta$-directions, and so on.

For comparison purposes, we shall also model the gas phase with the linearised NSF equations. In this case, the linear, dimensionless, steady-state NSF equations are (2.2)–(2.4) along with the closure

(2.14a,b)$$\begin{gather} \sigma_{rr} ={-} 2 {Kn}\frac{\partial v_r}{\partial r}, \quad \sigma_{r\theta} ={-} {Kn}\left(\frac{\partial v_\theta}{\partial r} - \frac{v_\theta}{r} + \frac{1}{r} \frac{\partial v_r}{\partial\theta}\right), \end{gather}$$

(2.15a,b)$$\begin{gather}q_r ={-} \frac{5}{2} \frac{Kn}{Pr} \frac{\partial T}{\partial r} \quad\text{and}\quad q_\theta ={-} \frac{5}{2} \frac{Kn}{Pr} \frac{1}{r} \frac{\partial T}{\partial \theta}. \end{gather}$$

2.2. Modelling of the liquid phase inside the droplet

The liquid phase inside the spherical droplet is modelled 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, for example, the textbook by 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

(2.16)\begin{gather} \frac{\partial v_r^{(\ell)}}{\partial r}+\frac{2 v_r^{(\ell)}}{r}+\frac{\mathscr{D} v_\theta^{(\ell)}}{r}=0, \end{gather}

(2.17a)\begin{gather} \frac{\partial p^{(\ell)}}{\partial r}+\frac{\partial\sigma_{rr}^{(\ell)}}{\partial r}+\frac{3\sigma_{rr}^{(\ell)}}{r}+\frac{\mathscr{D} \sigma_{r\theta}^{(\ell)}}{r}=0, \end{gather}

(2.17b)\begin{gather}\frac{\partial\sigma_{r\theta}^{(\ell)}}{\partial r}+\frac{3\sigma_{r\theta}^{(\ell)}}{r}-\frac{1}{2r}\frac{\partial\sigma_{rr}^{(\ell)}}{\partial \theta}+\frac{1}{r}\frac{\partial p^{(\ell)}}{\partial \theta}=0, \end{gather}

(2.18)\begin{gather} \frac{\partial q_r^{(\ell)}}{\partial r}+\frac{2q_r^{(\ell)}}{r}+\frac{\mathscr{D}q_{\theta}^{(\ell)}}{r}=0{,} \end{gather}

with

(2.19a,b)$$\begin{gather} \sigma_{rr}^{(\ell)} ={-} 2 \varLambda_\mu {Kn} \frac{\partial v_r^{(\ell)}}{\partial r}, \quad \sigma_{r\theta}^{(\ell)} ={-} \varLambda_\mu {Kn} \left(\frac{\partial v_\theta^{(\ell)}}{\partial r} - \frac{v_\theta^{(\ell)}}{r} + \frac{1}{r} \frac{\partial v_r^{(\ell)}}{\partial \theta}\right), \end{gather}$$

(2.20a,b)$$\begin{gather}q_r^{(\ell)} ={-} \frac{5}{2} \varLambda_\kappa \frac{Kn}{Pr} \frac{\partial T^{(\ell)}}{\partial r} \quad\text{and}\quad q_\theta^{(\ell)} ={-} \frac{5}{2} \varLambda_\kappa \frac{Kn}{Pr} \frac{1}{r} \frac{\partial T^{(\ell)}}{\partial \theta}. \end{gather}$$

The superscript ‘$(\ell )$’ in (2.16)–(2.19a,b) 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 (2.16)–(2.19a,b) are as follows. Firstly,

(2.21a,b)\begin{equation} 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}} \end{equation}

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

(2.22)\begin{equation} p^{(\ell)} = \frac{\hat{p}^{(\ell)}-\hat{p}^{(\ell)}_0}{\hat{p}_0} \end{equation}

is the dimensionless deviation in the pressure of the liquid droplet from its pressure in the equilibrium state $\hat {p}_0^{(\ell )}$, which is given by the pressure inside a stationary droplet in a quiescent environment, i.e. $\hat {p}_0^{(\ell )} = \hat {p}_0 + 2\hat {\gamma }/\hat {R}_0$. Note that the reference pressure inside a stationary droplet is higher than the reference pressure outside the droplet because of the surface tension $\hat {\gamma }$. Furthermore, in (2.16)–(2.19a,b),

(2.23a,b)\begin{equation} \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} \end{equation}

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

(2.24a,b)\begin{equation} 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}} \end{equation}

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

(2.25)\begin{equation} T^{(\ell)} = \frac{\hat{T}^{(\ell)}-\hat{T}_0}{\hat{T}_0} \end{equation}

is the dimensionless deviation in the temperature of the liquid droplet from the temperature in the reference state $\hat {T}_0$; $\varLambda _\mu = \hat {\mu }^{(\ell )} / \hat {\mu }$ is the ratio of the viscosity of the liquid to the viscosity of the gas and $\varLambda _\kappa = \hat {\kappa }^{(\ell )} / \hat {\kappa }$ is the ratio of the thermal conductivity of the liquid to the thermal conductivity of the gas. Similarly to the above, owing to the linearisation, $\hat {\mu }^{(\ell )}$ and $\hat {\kappa }^{(\ell )}$ are the viscosity and thermal conductivity of the liquid at the reference temperature $\hat {T}_0$; consequently, $\hat {\mu }$, $\hat {\mu }^{(\ell )}$, $\hat {\kappa }$, $\hat {\kappa }^{(\ell )}$, $\varLambda _\mu$ and $\varLambda _\kappa$ are all constant in the present work. The field variables with hats and superscript ‘${(\ell )}$’ can be identified as the field variables of the original NSF equations. Furthermore, (2.16) can be identified as the equation of continuity, (2.17a) and (2.17b) as the momentum balance equations in the $r$- and $\theta$-directions, respectively, and (2.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, Gupta & Struchtrup (2018a) and 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,b) 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 (2.2)–(2.12) – for the problem under consideration – at the interface (i.e. at $r=1$) read (Rana et al. 2021a)

(2.26a)$$\begin{gather} v_r = v_r^{I} = 0, \end{gather}$$

(2.26b)$$\begin{gather}q_r ={-} \frac{\chi}{2 - \chi} \sqrt{\frac{2}{\rm \pi}} \left(2 \mathcal{T} + \frac{1}{2} \sigma_{rr} + \frac{5}{28} R_{rr} + \frac{1}{15} \varDelta -\frac{1}{6} \varPhi_{rrrr} \right), \end{gather}$$

(2.26c)$$\begin{gather}m_{rrr} = \frac{\chi}{2 - \chi}\sqrt{\frac{2}{\rm \pi}} \left(\frac{2}{5} \mathcal{T} - \frac{7}{5} \sigma_{rr} - \frac{1}{14} R_{rr} + \frac{1}{75} \varDelta - \frac{13}{15} \varPhi_{rrrr}\right), \end{gather}$$

(2.26d)$$\begin{gather}\varPsi_{rrr} = \frac{\chi}{2 - \chi}\sqrt{\frac{2}{\rm \pi}} \left(\frac{6}{5} \mathcal{T} + \frac{9}{5} \sigma_{rr} - \frac{93}{70} R_{rr} + \frac{1}{5} \varDelta + \frac{11}{15} \varPhi_{rrrr}\right), \end{gather}$$

(2.26e)$$\begin{gather}\varOmega_r = \frac{\chi}{2 - \chi}\sqrt{\frac{2}{\rm \pi}} \left(8 \mathcal{T} + 2 \sigma_{rr} - R_{rr} - \frac{4}{3} \varDelta - \frac{2}{3} \varPhi_{rrrr}\right), \end{gather}$$

(2.26f)$$\begin{gather}\sigma_{r\theta} ={-} \frac{\chi}{2 - \chi}\sqrt{\frac{2}{\rm \pi}} \left( \mathscr{V}_\theta +\frac{1}{5} q_\theta +\frac{1}{2}m_{rr\theta} - \frac{1}{14} \varPsi_{rr\theta} - \frac{1}{70} \varOmega_r\right), \end{gather}$$

(2.26g)$$\begin{gather}R_{r\theta} ={-} \frac{\chi}{2 - \chi}\sqrt{\frac{2}{\rm \pi}} \left( -\mathscr{V}_\theta + \frac{11}{5} q_\theta + \frac{1}{2} m_{rr\theta} + \frac{13}{14} \varPsi_{rr\theta} + \frac{13}{70} \varOmega_r\right), \end{gather}$$

(2.26h)$$\begin{gather}\varPhi_{rrr\theta} ={-} \frac{\chi}{2 - \chi}\sqrt{\frac{2}{\rm \pi}} \left( -\frac{4}{7} \mathscr{V}_\theta -\frac{12}{35} q_\theta +\frac{9}{7} m_{rr\theta} - \frac{2}{49} \varPsi_{rr\theta} - \frac{2}{245} \varOmega_r\right), \end{gather}$$

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$.

While solving the gas phase with the linearised NSF equations (for comparison purposes), the appropriate boundary conditions are (2.26a) and the ones obtained by ignoring the higher-order moments in boundary conditions (2.26b) and (2.26f). Thus, for the linearised NSF equations, the appropriate boundary conditions are (2.26a) and

(2.27)$$\begin{gather} \frac{\chi}{2 - \chi} \sqrt{\frac{2}{\rm \pi}} \left(2 \mathcal{T} + \frac{1}{2} \sigma_{rr}\right) + q_r = 0, \end{gather}$$

(2.28)$$\begin{gather}\frac{\chi}{2 - \chi}\sqrt{\frac{2}{\rm \pi}} \left( \mathscr{V}_\theta +\frac{1}{5} q_\theta\right) + \sigma_{r\theta} = 0. \end{gather}$$

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.

(2.29)\begin{equation} v_r^{(\ell)} = 0 \quad\text{at}\ r=1. \end{equation}

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

(2.30)\begin{equation} q_r = q_r^{(\ell)} \quad\text{and}\quad \sigma_{r\theta} = \sigma_{r\theta}^{(\ell)} \quad\text{at}\ r=1. \end{equation}

It is worthwhile noting that the density ratio (of the liquid to the gas) – despite being an important parameter in gas–liquid interfacial flows – appear neither in the governing equations (2.2)–(2.9) and (2.16)–(2.18), nor in the constitutive relations (2.10)–(2.12), (2.19a,b) and (2.19a,b) nor in boundary conditions (2.26), (2.29) and (2.30) due to our assumptions of no phase change and of the surface tension force being strong enough to keep the spherical shape and size of the droplet unchanged. Hence, the density ratio does not play any role in the present work. Nevertheless, when taking the effect of the surface tension force into account (i.e. when the droplet is allowed to change its shape) and/or taking the phase-change into account (i.e. when the droplet is allowed to change its size), boundary conditions (2.26a) and (2.29) need to be modified appropriately. The modified boundary conditions will consist of the densities of both liquid and gas, and hence the density ratio. In addition, when accounting for the effects of the surface tension forces, the stress boundary condition (the second of (2.30)) also needs to be changed. A general force balance condition at an interface between two fluids (labelled ‘1’ and ‘2’) is given by (Landau & Lifshitz 1959)

(2.31)\begin{equation} (\hat{p}_1-\hat{p}_2+\hat{\gamma}_c \hat{\varkappa}) n_i ={-}\left[\hat{\sigma}_{ij}^{(1)} - \hat{\sigma}_{ij}^{(2)}\right] n_j + \frac{\partial \hat{\gamma}_c}{\partial \hat{x}_i}, \end{equation}

where $\hat {p}_1$ and $\hat {p}_2$ are the pressures exerted on the interface by the fluids ‘1’ and ‘2’, respectively; $\hat {\gamma }_c$ is the surface tension coefficient; $\hat {\varkappa }$ is the local surface curvature; $n_i$ is the unit normal at the interface; and $\hat {\sigma }_{ij}^{(1)}$ and $\hat {\sigma }_{ij}^{(2)}$ are the stress tensors of the fluids ‘1’ and ‘2’, respectively. In the present work, we assume that the droplet remains spherical without any growth or shrinkage and that it does not deform either; in other words, the radius of the droplet $\hat {R}_0$ is assumed to remain constant. Therefore, we neglect the non-equilibrium force balance in the normal direction as described in (2.31). It is justified since the equilibrium aspect has already been considered when defining $\hat {p}_0^{(\ell )}$ above. Additionally, for the force balance in the tangential direction, we multiply (2.31) by the unit tangent vector $t_i$ at the interface. This multiplication makes the left-hand side of the resulting equation vanish. By disregarding the effect of surface tension forces, specifically the term $\partial \hat {\gamma }_c / \partial \hat {x}_i$ (commonly referred to as the Marangoni surface tension gradient) in (2.31), the boundary condition simplifies to $\hat {\sigma }_{nt}^{(1)} = \hat {\sigma }_{nt}^{(2)}$, which in the dimensionless form is the same as the second boundary condition in (2.30) considered in the present work.

3. Analytic solution methodology

We solve the equations derived in § 2 analytically by following a method proposed in Torrilhon (2010) and Rana et al. (2018b). 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

(3.1)\begin{align} \left. \begin{aligned} v_r(r,\theta) & = \mathbb{v}_1(r) \cos{\theta}, & \quad v_\theta(r,\theta) & ={-} \mathbb{v}_2(r) \sin{\theta}\\ p(r,\theta) & = \mathbb{p}(r) \cos{\theta}, & \quad T(r,\theta) & = \mathbb{T}(r) \cos{\theta},\\ \sigma_{rr}(r,\theta) & = \mathbb{s}_1(r) \cos{\theta}, & \quad \sigma_{r\theta}(r,\theta) & = \mathbb{s}_2(r) \sin{\theta},\\ q_r(r,\theta) & = \mathbb{q}_1(r) \cos{\theta}, & \quad q_\theta(r,\theta) & ={-}\mathbb{q}_2(r) \sin{\theta},\\ m_{rrr}(r,\theta) & = \mathbb{m}_1(r) \cos{\theta}, & \quad m_{rr\theta}(r,\theta) & = \mathbb{m}_2(r) \sin{\theta},\\ R_{rr}(r,\theta) & = \mathbb{R}_1(r) \cos{\theta}, & \quad R_{r\theta}(r,\theta) & = \mathbb{R}_2(r) \sin{\theta},\\ \varDelta(r,\theta) & = \mathbb{d}(r) \cos{\theta}, \end{aligned} \right\} \end{align}

(3.2)\begin{align} \left.\begin{aligned} v_r^{(\ell)}(r,\theta) & = \mathbb{v}_1^{(\ell)}(r) \cos{\theta}, & \quad v_\theta^{(\ell)}(r,\theta) & ={-} \mathbb{v}_2^{(\ell)}(r) \sin{\theta}, \\ p^{(\ell)}(r,\theta) & = \mathbb{p}^{(\ell)}(r) \cos{\theta}, & \quad T^{(\ell)}(r,\theta) & = \mathbb{T}^{(\ell)}(r) \cos{\theta}, \end{aligned} \right\} \end{align}

where the functions $\mathbb {v}_1, \mathbb {v}_2, \mathbb {p}, \mathbb {T}, \mathbb {s}_1, \mathbb {s}_2, \mathbb {q}_1, \mathbb {q}_2, \mathbb {m}_1, \mathbb {m}_2, \mathbb {R}_1, \mathbb {R}_2, \mathbb {d}, \mathbb {v}_1^{(\ell )}, \mathbb {v}_2^{(\ell )}, \mathbb {p}^{(\ell )}$ and $\mathbb {T}^{(\ell )}$ are the functions of $r$ alone.

The above ansatzes for the field variables are inserted in (2.2)–(2.12) and in (2.16)–(2.19a,b). After simplification ($\cos {\theta }$ and $\sin \theta$ in each equation get cancelled), 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

(3.3a–d)\begin{align} \mathbb{v}_1^{(\ell)}(r) = b_1 + \frac{b_2 r^2}{2}, \quad \mathbb{v}_2^{(\ell)}(r) = b_1 + b_2 r^2, \quad \mathbb{p}^{(\ell)}(r) = 5 b_2 \varLambda_\mu {Kn} \, r, \quad \mathbb{T}^{(\ell)}(r) = b_3 r, \end{align}

where $b_1$, $b_2$ and $b_3$ are the integration constants, which are computed using the interface conditions (2.29) and (2.30). While obtaining solution (3.3a–d), 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 the 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 (2.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

To access the validity of the findings of this work, we first present the results on the drag force and compare them with those obtained from an experiment. After validating the drag force, we shall present the results on the physical field variables, which are often difficult to measure through experiments.

4.1. Drag force

Before computing the drag force with the analytic solution obtained in the present work, let us first comment on the experimental data, which will be used to validate our analytical findings on the drag force. The experimental data is actually from the famous oil-drop experiment conducted in 1909 by R. A. Millikan, the Nobel laureate in Physics 1923, to measure the electric charge carried by a single electron. Through this experiment, he also computed the drag force on oil drops of different radii falling in the air at terminal velocity (Millikan 1923). Millikan's experimental data have been fitted by several researchers to obtain an empirical formula for the drag force in the Knudsen–Weber form (Knudsen & Weber 1911)

(4.1)\begin{equation} F= F_{Stokes}\left[\frac{1}{1+{Kn}\left(\mathfrak{a} + \mathfrak{b} \, \mathrm{e}^{-\mathfrak{c}/{Kn}}\right)}\right], \end{equation}

where $F_{Stokes}=6{\rm \pi} \,{Kn} \, u_{\infty }$ is the Stokes drag with $u_{\infty } = \hat {u}_\infty / \sqrt {\hat {R}\hat {T}_0}$ being the far-field dimensionless velocity of gas approaching the droplet (for simplicity, we have taken ${u_{\infty } = 1}$ in this article), and $\mathfrak {a}$, $\mathfrak {b}$ and $\mathfrak {c}$ are the experimentally determined constants. It turns out that for an oil drop of size comparable with the mean free path of the air (i.e. in the transition regime), the experimental data for oil drops in air from Millikan's oil-drop experiment requires these constants to be (Kennard 1938)

(4.2a–c)\begin{equation} \mathfrak{a} = 1.23, \quad \mathfrak{b} = 0.41 \quad\text{and}\quad \mathfrak{c} = 0.88. \end{equation}

The most accurate raw data from Millikan's experiment were reviewed later by Allen & Raabe (1982) using very precise values of the physical constants known at that time and the nonlinear least-squares fitting technique. With these, the new values of the constants $\mathfrak {a}$, $\mathfrak {b}$ and $\mathfrak {c}$ obtained by Allen & Raabe (1982) are

(4.3a–c)\begin{equation} \mathfrak{a} = 1.155 \pm 0.008, \quad \mathfrak{b} = 0.471 \pm 0.011 \quad\text{and}\quad \mathfrak{c} = 0.596 \pm 0.050. \end{equation}

An improved version of Millikan's experimental apparatus was designed and built by Allen & Raabe (1985) aiming to replace the oil drops by solid spheres in Millikan's experiment and to determine new values of the constants $\mathfrak {a}$, $\mathfrak {b}$ and $\mathfrak {c}$ corresponding to the drag force on a solid sphere suspended in air. By taking three different types of solid spherical particles – namely, polystyrene latex-divinylbenzene particles, polyvinyltoluene particles and polystyrene latex particles – with the Knudsen number ranging from $0.03$ to $7.2$ in this experiment, they found the values of the constants $\mathfrak {a}$, $\mathfrak {b}$ and $\mathfrak {c}$ to be

(4.4a–c)\begin{equation} \mathfrak{a} = 1.142 \pm 0.0024, \quad \mathfrak{b} = 0.558 \pm 0.0024 \quad\text{and}\quad \mathfrak{c} = 0.999 \pm 0.0212. \end{equation}

Using the modulated dynamic light scattering technique – a technique fundamentally different from the ones used by Millikan and by Allen and Raabe, Hutchins, Harper & Felder (1995) also measured the drag force on spherical polystyrene latex (solid) particles suspended in dry air for the Knudsen number values ranging from $0.06$ to $500$, and found the values of the constants $\mathfrak {a}$, $\mathfrak {b}$ and $\mathfrak {c}$ to be

(4.5a–c)\begin{equation} \mathfrak{a} = 1.2310 \pm 0.0022, \quad \mathfrak{b} = 0.4695 \pm 0.0037 \quad\text{and}\quad \mathfrak{c} = 1.1783 \pm 0.0091. \end{equation}

To compare the drag force computed analytically in the present work, we use the data resulting from formula (4.1) with the constants $\mathfrak {a}$, $\mathfrak {b}$ and $\mathfrak {c}$ given in (4.2a–c) (taken from the textbook Kennard (1938)) and in (4.2a–c) (determined by Allen & Raabe (1982)), which give accurate fit to the experimental data from Millikan's oil-drop experiment. In addition, we shall also include the data resulting from formula (4.1) on taking the constants $\mathfrak {a}$, $\mathfrak {b}$ and $\mathfrak {c}$ given in (4.2a–c) (determined by Allen & Raabe (1985)) and in (4.2a–c) (determined by Hutchins et al. (1995)), which give accurate fit to the data obtained from experiments performed with the (aforementioned) solid spheres suspended in air. Although one may remonstrate that comparing the results on the drag force obtained in the present work, which is on a monatomic gas flow over a liquid droplet, with those from Millikan's oil-drop experiment, which was on oil drops falling in air, would not be fair due to the fact that air is not a monatomic gas, the viscosity of the monatomic gas in the present work is taken to be the same as that of air used in Millikan's oil-drop experiment – making the comparison sensible. Furthermore, we are not aware of any theoretical/experimental work, which presents the data for monatomic gas flow over a liquid droplet. Therefore, it is reasonable to compare the results on the drag force obtained in the present work with the data from Millikan's oil-drop experiment.

In the present work, the drag force is computed analytically as follows. The drag force in the integral form can be written as (Torrilhon 2010)

(4.6)\begin{equation} F ={-}2{\rm \pi}\int_{0}^{\rm \pi} {\boldsymbol{z}}_{x}(\theta)\boldsymbol{\cdot}(\boldsymbol{P}(r,\theta)\boldsymbol{\cdot} {\boldsymbol{k}})\sin{\theta}\,\mathrm{d}\theta, \end{equation}

where ${\boldsymbol {z}}_{x}(\theta )=(\cos {\theta },-\sin {\theta },0)$, ${\boldsymbol {k}}=(1,0,0)$ and $\boldsymbol {P}=p\boldsymbol {I}+\boldsymbol {\sigma }$ is the pressure tensor. The force on the surface of the liquid droplet is represented by $\boldsymbol {P}\boldsymbol {\cdot }{\boldsymbol {k}}$ and its scalar product with ${\boldsymbol {z}}_{x}$ gives the component of the force in the direction of the flow. Substituting ansatz (3.1) in (4.6), we obtain

(4.7)\begin{equation} F=\frac{4{\rm \pi}}{3}\left[-\mathbb{p}(1)-\mathbb{s}_{1}(1)+2\mathbb{s}_{2}(1)\right]. \end{equation}

This, on simplifying further, yields

(4.8)\begin{equation} F={-}2{\rm \pi} C_{1} {Kn}, \end{equation}

where the constant $C_{1}$ is determined using the boundary conditions given in (2.3).

Figure 2 exhibits the drag force normalised with the Stokes drag $F_{Stokes}$ 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 $\varLambda _\kappa = 100$ and for different values of the viscosity ratio $\varLambda _\mu$. In addition, the explicit values of the normalised drag force for the same value of the thermal conductivity ratio $\varLambda _\kappa = 100$ and for different values of the Knudsen number and viscosity ratio are also given in table 1.

Figure 2. Drag force (normalised with the Stokes drag $F_{Stokes} =6{\rm \pi} \,{Kn}\,u_\infty$) on a liquid droplet from the linearised NSF and R26 theories as a function of the Knudsen number for a fixed thermal conductivity ratio $\varLambda _\kappa =100$. The square and star symbols denote the experimental data from Millikan's oil-drop experiment fitted by the empirical formulae of Kennard (1938) ((4.1) with (4.2a–c)) and Allen & Raabe (1982) ((4.1) with (4.2a–c)), respectively. The disk and diamond symbols denote the data from experiments performed with solid spherical particles in air fitted by the empirical formulae of Allen & Raabe (1985) ((4.1) with (4.2a–c)) and Hutchins et al. (1995) ((4.1) with (4.2a–c)), respectively, and are included just for comparison.

Table 1. Drag force normalised with the Stokes drag, $F_{Stokes}=6{\rm \pi} \,{Kn}\,u_\infty$, for different values of the Knudsen number and viscosity ratios at a fixed thermal conductivity ratio $\varLambda _\kappa =100$.

It is clear from figure 2 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 non-zero positive value) as the Knudsen number approaches infinity. On the other hand, it can be seen from figure 2 and table 1 that with the increase in the viscosity ratio $\varLambda _\mu$, the drag force approaches unity for smaller values of the Knudsen number. Since the viscosity ratio of an oil to air is somewhere in between $10^2$ to $10^3$ (depending on the oil), the figure clearly 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 data from Millikan's oil-drop experiment for the Knudsen number values as big as $1$. It is also worthwhile noting that the oil used in Millikan's experiment not only has a very high viscosity but also has a high thermal conductivity; therefore the results on the drag force from Millikan's experiment (and also those in the present work for high viscosity ratios) are very close to the drag force on a solid sphere (of the same size as the size of the liquid droplet) – as also manifested by the propinquity of symbols in figure 2.

Figure 3 depicts the drag force computed with the LR26 equations and normalised with the drag force given by (4.9) again for a fixed value of the thermal conductivity ratio $\varLambda _\kappa = 100$ and for different values of the viscosity ratio $\varLambda _\mu$. In addition, the explicit values of the normalised drag force for the same value of the thermal conductivity ratio $\varLambda _\kappa = 100$ and for different values of the Knudsen number and viscosity ratio are also given in table 2. 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

(4.9)\begin{equation} F_{HB}=F_{Stokes} \times \frac{1 + \dfrac{2}{3 \varLambda_\mu}}{1 + \dfrac{1}{\varLambda_\mu}}. \end{equation}

From figure 3 and table 2, the drag force normalised with the drag force $F_{HB}$ (given by (4.9)) for very small Knudsen numbers is very close to unity for all values of the viscosity ratio $\varLambda _\mu$, which is not the case when the drag force normalised with the Stokes drag (cf. figures 2 and 3 for small Knudsen numbers). However, for large Knudsen numbers, the drag force scaled with $F_{HB}$ is more or less same as the drag force scaled with $F_{Stokes}$. Indeed, from (4.9), as $\varLambda _\mu \to \infty$, $F_{HB} \to F_{Stokes}$; therefore for large viscosity ratios the scaling with $F_{HB}$ or $F_{Stokes}$ does not matter. Hence the dotted green lines and solid red lines in figures 2 and 3 coincide with each other.

Figure 3. Drag force on the liquid droplet normalised with the drag force given by (4.9). The thermal conductivity ratio is $\varLambda _\kappa =100$.

Table 2. Drag force normalised with the drag force given by (4.9) for different values of the Knudsen number and viscosity ratio with the thermal conductivity ratio being fixed at 100.

4.2. Flow fields: temperature and heat flux

Figures 4, 5 and 6 illustrate the heat flux lines and temperature contours in the $\hat {y}=0$ plane for a fixed viscosity ratio $\varLambda _\mu = 100$ and for the thermal conductivity ratios $\varLambda _\kappa =1$, $10$ and $100$, respectively. Figures 4(a,d), 5(a,d) and 6(a,d) display the results for the Knudsen number ${Kn}=0.09$. Figures 4(b,e), 5(b,e) and 6(b,e) display the results for the Knudsen number ${Kn}=0.36$. Figures 4(c, f), 5(c, f) and 6(c, f) display the results for the Knudsen number ${Kn}=0.9$. Figures 4(a–c), 5(a–c) and 6(a–c) illustrate the results obtained with the LR26 equations for the rarefied gas flow outside the droplet and with the linear NSF equations for the liquid inside the droplet. Figures 4(d–f), 5(d–f) and 6(d–f) are obtained with the NSF equations for both the liquid inside the droplet and gas flow outside the droplet. Comparing the respective panels in figures 4(a–c), 5(a–c) and 6(a–c) with those in figures 4(d–f), 5(d–f) and 6(d–f), the R26 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 figures 4(d–f), 5(d–f) and 6(d–f) show that the NSF equations cannot capture the non-Fourier heat transfer in the external (rarefied gas) flow. In fact, the R26 equations for the external flow give high (low) temperature on the front (back) side of the liquid droplet while the NSF equations for the external flow predict the opposite of this (cf. the corresponding panels in figures 4a–c, 5a–c and 6a–c, with figures 4d–f, 5d–f and 6d–f). In addition, comparing figures 4, 5 and 6 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 4. Heat flux lines plotted on top of the temperature contours for the viscosity ratio $\varLambda _\mu =100$ and for different values of the Knudsen number: ${Kn} = 0.09$ (a,d), ${Kn} = 0.36$ (b,e) and ${Kn} = 0.9$ (c, f). The results for the liquid phase (internal flow) have been computed with the linear NSF equations in all the cases while those for the gas phase (external flow) have been computed with the LR26 equations (a–c) and with the linear NSF equations (d–f). The thermal conductivity ratio is $\varLambda _\kappa =1$.

Figure 5. Same as figure 4 but for the thermal conductivity ratio $\varLambda _\kappa =10$.

Figure 6. Same as figure 4 but for the thermal conductivity ratio $\varLambda _\kappa =100$.

In order to get an insight of the temperature and heat flux at a given position, we plot the temperature and radial component of the heat flux (both divided by $\cos {\theta }$ to understand the results for any angle $\theta$) with respect to the position $r$ in figures 7 and 8, respectively, again for a fixed viscosity ratio $\varLambda _\mu = 100$ and for the thermal conductivity ratios $\varLambda _\kappa =1$, $10$ and $100$. Figures 7(a,b) and 8(a,b) display the results for the Knudsen number ${Kn}=0.09$. Figures 7(c,d) and 8(c,d) display the results for the Knudsen number ${Kn}=0.36$. Figures 7(e,f) and 8(e,f) display the results for the Knudsen number ${Kn}=0.9$. Figures 7(a,c,e) and 8(a,c,e) illustrate the results obtained with the LR26 equations for the rarefied gas flow outside the droplet and with the linear NSF equations for the liquid inside the droplet while figures 7(b,d,f) and 8(b,d,f) display the results obtained with the linear NSF equations for both liquid inside the droplet and gas outside the droplet. The droplet interface has been demarcated by a vertical black line at $r=1$ in all panels of figures 7 and 8. For the liquid phase, the analytic solution for the temperature and radial component of the heat flux can actually be written quite easily from the first equation in (2.19a,b), the fourth equation in (3.2) and the third equation in (3.3a–d), and they read

(4.10a,b)\begin{equation} T^{(\ell)} = b_3 r \cos{\theta} \quad\text{and}\quad q_r^{(\ell)} ={-} \frac{5}{2} \varLambda_\kappa \frac{Kn}{Pr} b_3 \cos{\theta}. \end{equation}

Therefore, $T^{(\ell )}/\cos {\theta }$ for the liquid phase is a linear function of $r$ while $q_r^{(\ell )}/\cos {\theta }$ for the liquid phase at a given Knudsen number and at a given thermal conductivity ratio is just a constant. Consequently, the curves on the left of the vertical black lines in figure 7 are straight lines and the curves on the left of the vertical black lines in each panel of figure 8 are horizontal lines.

Figure 7. Dimensionless deviations in the temperature (scaled with $\cos {\theta }$) as a function of position $r$ for a fixed viscosity ratio $\varLambda _\mu =100$ and for different values of the Knudsen number: ${Kn} = 0.09$ (a,b), ${Kn} = 0.36$(c,d) and ${Kn = 0.9}$ (e, f). The vertical black line at $r=1$ demarcates the interface between the liquid and gas. 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 (a,c,e) and with the NSF equations (b,d, f).

Figure 8. Dimensionless radial heat flux (scaled with $\cos {\theta }$) as a function of position $r$ for a fixed viscosity ratio $\varLambda _\mu =100$ and for different values of the Knudsen number: ${Kn} = 0.09$ (a,b), ${Kn} = 0.36$ (c,d) and ${Kn = 0.9}$ (e, f). The vertical black line at $r=1$ demarcates the interface between the liquid and gas. 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 (a,c,e) and with the NSF equations (b,d, f).

The temperature contours depicted in figures 4, 5 and 6 actually correspond to temperature profiles exhibited in figure 7 by red, black and green lines, respectively. To see this, let us, for example, fix $\theta = 0$ and let us focus on figure 4 and the red lines in figure 7 (figures 5 and 6 and the lines of other colours in figure 7 can be understood analogously). From (the red lines in) figure 7, the temperature of the liquid droplet decreases when moving from its centre towards its interface. However, on comparing figure 7(a,c,e) and figure 7(b,d, f), one finds that the temperature of the external gas predicted by the LR26 equations (figure 7a,c,e) increases on moving away from the interface whereas that predicted by the linear NSF equations (figure 7a,c,e) decreases on moving away from the interface. These are in agreement with the temperature contours portrayed in figure 4 (notice the contours for $\theta =0$ when moving away from the centre of the droplet). Furthermore, figure 7 clearly shows a non-zero temperature jump $\mathcal {T} = T - T^{(\ell )}$ at the interface. From figure 7(b,d, f), the temperature of the liquid (gas) predicted by the linear NSF equations is minimum (maximum and positive) at the interface, and hence the temperature jump $\mathcal {T} = T - T^{(\ell )}$ predicted by the NSF equations is always positive. This is not the case when the LR26 equations are used for the external flow (see figure 7a,c,e). In fact, when the LR26 equations are used for the external flow, the temperature jump $\mathcal {T} = T - T^{(\ell )}$ is negative in most of the cases but positive in some cases.

Figure 8 exhibits the radial component of the heat flux (scaled with $\cos {\theta }$) with respect to the position $r$. As described above, the radial component of the heat flux for the liquid phase is constant (notice the horizontal lines on the left of the vertical black line in figure 8). On the contrary, the radial component of the heat flux is non-constant for the gas phase (see the curves on the right of the vertical black line in figure 8), and is decreasing monotonically for $r>1$.

4.3. Flow fields: pressure and velocity

Figures 9, 10 and 11 illustrate the velocity streamlines plotted over the pressure contours in the $\hat {y}=0$ plane for a fixed thermal conductivity ratio $\varLambda _\kappa = 100$ and for the viscosity ratios $\varLambda _\mu = 1$, $10$ and $100$, respectively. Figures 9(a,d), 10(a,d) and 11(a,d) denote the results for the Knudsen number ${Kn}=0.09$. Figures 9(b,e), 10(b,e) and 11(b,e) denote the results for the Knudsen number ${Kn}=0.36$. Figures 9(c,f), 10(c,f) and 11(c,f) denote the results for the Knudsen number ${Kn=0.9}$. Figures 9(a–c), 10(a–c) and 11(a–c) illustrate the results obtained with the LR26 equations for the rarefied gas flow outside the droplet and with the linear NSF equations for the liquid inside the droplet while figures 9(d–f), 10(d–f) and 11(d–f) are obtained with the linear NSF equations for both liquid inside the droplet and gas flow outside the droplet. The streamlines in figures 9, 10 and 11 exhibit that, owing to the motion of the gas in the $\hat {z}$-direction, the liquid inside the droplet near the top (close to $(r,\theta,\phi ) = (1, {\rm \pi}/2, 0)$) and bottom (close to $(r,\theta,\phi ) = (1, {\rm \pi}/2, {\rm \pi})$) surfaces of the droplet starts moving in direction of the gas flow but since the liquid cannot move outside of the droplet, it flows in the opposite directions near the centre of the droplet, hence forming two counter-rotating vortices inside the droplet. The pressure contours in figures 9, 10 and 11 show that the magnitude of the pressure at a point (inside or outside of the liquid droplet) increases with an increase in the Knudsen number.

Figure 9. Velocity streamlines plotted over the pressure contours for the thermal conductivity ratio $\varLambda _\kappa =100$ and for different values of the Knudsen number: ${Kn} = 0.09$ (a,d), ${Kn} = 0.36$ (b,e) and ${Kn = 0.9}$ (c, f). 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 (a–c) and with the NSF equations (d–f). The viscosity ratio is $\varLambda _\mu = 1$.

Figure 10. Same as figure 9 but for the viscosity ratio $\varLambda _\mu =10$.

Figure 11. Same as figure 9 but for the viscosity ratio $\varLambda _\mu =100$.

Furthermore, a comparison of figures 9(a–c), 10(a–c) and 11(a–c) with figures 9(d–f), 10(d–f) and 11(d–f) reveals that there are practically no differences in the magnitudes of the pressure contours inside the liquid droplet but there are noticeable differences in the magnitudes of the pressure contours for the gas phase in the two cases (i.e. when using the LR26 equations ( figures 9a–c, 10a–c and 11a–c) and the linear NSF equations ( figures 9d–f, 10d–f and 11d–f) for the gas phase). Nonetheless, the differences in the magnitudes of the pressure contours for the gas phase decrease with an increase in the Knudsen number.

Similarly to the above, in order to get an insight of the pressure and velocity at a given position, we plot the pressure (divided by $\cos {\theta }$ to understand the results for any angle $\theta$) with respect to the position $r$ in figure 12 and the $z$-component of the velocity for a fixed $\theta ={\rm \pi} /2$ with respect to the position $r$ in figure 13. For both figures 12 and 13, the thermal conductivity ratio is again fixed to $\varLambda _\kappa = 100$ and the viscosity ratios are taken as $\varLambda _\mu =1$, $10$ and $100$. Figures 12(a,b) and 13(a,b) exhibit the results for the Knudsen number ${Kn}=0.09$. Figures 12(c,d) and 13(c,d) exhibit the results for the Knudsen number ${Kn}=0.36$. Figures 12(e,f) and 13(e,f) exhibit the results for the Knudsen number ${Kn=0.9}$. Figures 12(a,c,e) and 13 (a,c,e) depict the results obtained with the LR26 equations for the rarefied gas flow outside the droplet and with the linear NSF equations for the liquid inside the droplet while figures 12(b,d, f) and 13(b,d, f) depict the results obtained with the linear NSF equations for both liquid inside the droplet and gas outside the droplet. The droplet interface has again been demarcated by a vertical black line at $r=1$ in figures 12 and 13. For the liquid phase, the analytic solution for the pressure and $z$-component of the velocity can also be written quite easily from equations one to three (3.2) and equations one to three (3.3a–d), and they read

(4.11)\begin{equation} p^{(\ell)} = 5 b_2 \varLambda_\mu {Kn} \, r \cos{\theta} \end{equation}

and

(4.12)\begin{equation} v_z^{(\ell)} = v_r^{(\ell)} \cos{\theta} - v_\theta^{(\ell)} \sin{\theta} = b_1 + b_2 r^2 - \frac{b_2 r^2}{2}\cos^2{\theta}. \end{equation}

Therefore, inside the droplet, $p^{(\ell )}/\cos {\theta }$ is a linear function of $r$ while $v_z^{(\ell )}$ is parabolic in $r$ for $\theta = {\rm \pi}/2$. Consequently, the curves on the left of the vertical black lines in figure 12 are straight lines and the curves on the left of the vertical black lines in figure 13 are parabolas. Figure 12 delineates that the dependence of the pressure in the gas on the viscosity ratio is negligible and that the pressure in the liquid increases with increasing the viscosity ratio, although the increase in pressure in the liquid for large viscosity ratios is also insignificant. It is important to note that figure 12 illustrates the deviations in the pressures in the liquid droplet and in the gas from the pressures in their respective equilibrium states, which are different for the liquid and gas due to the Laplace pressure. The dimensional pressure in the liquid from (2.22) is given by

(4.13)\begin{equation} \hat{p}^{(\ell)} = \hat{p}_0 [p^{(\ell)} + 1 + 2\gamma \,{Kn}] \end{equation}

where

(4.14)\begin{equation} \gamma = \frac{\hat{\gamma}}{\hat{\mu}_0 \sqrt{\hat{R}\hat{T}_0}}, \end{equation}

is the dimensionless surface tension. The dimensionless deviation in the liquid pressure $p^{(\ell )}$ can be made as small as desired in comparison with $1+2\gamma \,{Kn}$ by taking $\hat {u}_\infty$ to be sufficiently small. This keeps our assumptions of the shape of the droplet being spherical and its size being fixed justifiable.

Figure 12. Dimensionless deviations in the pressure (scaled with $\cos {\theta }$) as a function of position $r$ for a fixed thermal conductivity ratio $\varLambda _\kappa =100$ and for different values of the Knudsen number: ${Kn} = 0.09$(a,b), ${Kn} = 0.36$ (c,d) and ${Kn} = 0.9$ (e, f). The vertical black line at $r=1$ demarcates the interface between the liquid and gas. 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 (a,c,e) and with the NSF equations (b,d, f).

Figure 13. The $z$-component of the (dimensionless) velocity as a function of position $r$ in the plane $\theta ={\rm \pi} /2$ for a fixed thermal conductivity ratio $\varLambda _\kappa =100$ and for different values of the Knudsen number: ${Kn} = 0.09$ (a,b), ${Kn} = 0.36$ (c,d) and ${Kn} = 0.9$ (e, f). The vertical black line at $r=1$ demarcates the interface between the liquid and gas. 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 (a,c,e) and with the NSF equations (b,d, f).

Figure 13 is actually a one-dimensional interpretation of the velocity streamlines shown in figures 9, 10 and 11. To see the interpretation, notice that figures 9, 10 and 11 show the results in the $\hat {y}=0$ plane while figure 13 displays the results in the $\hat {z}=0$ (or $\theta ={\rm \pi} /2$) plane, and hence we can see the results in both figures along a common line, which is the $x$-axis. Figure 13 reveals that as we start from the centre of the droplet along the positive $x$-axis, the $z$-component of the velocity of the liquid inside the droplet is initially negative, then at some point it becomes zero and then it turns positive as we approach towards the surface of the droplet, rendering a vortex in the upper half of the droplet. Further, on moving outside of the droplet (along the positive $x$-axis), the $z$-component of the velocity of the gas remains always positive and approaches $u_\infty =1$ as $r$ approaches $\infty$. This is exactly what figures 9, 10 and 11 show if we focus on the streamlines falling on the positive $x$-axis in figures 9, 10 and 11. It is evident from figure 13 that there is a discontinuity in the $z$-components of the velocities of the liquid and gas at the interface of the liquid droplet – referred to as the velocity slip – and that the magnitude of the discontinuity increases with increase in the Knudsen number. Furthermore, figure 13 also reveals that the velocity of the liquid droplet approaches zero as the viscosity ratio $\varLambda _\mu$ becomes larger and larger, which is consistent with the case of a rarefied gas flow over a solid sphere.

4.4. Discussion on results for some commonly used fluids

In the above subsections, we have analysed the drag force on the droplet and flow profiles (streamlines, heat flux lines, velocity, temperature, pressure and heat flux) of the liquid and gas over a range of dimensionless parameters, namely, the Knudsen number ${Kn}$, the viscosity ratio $\varLambda _\mu$ and the thermal conductivity ratio $\varLambda _\kappa$. In this section, we pick some commonly used fluids to gauge the applicability/limitations of the present work.

Table 3 exhibits the parameters, namely the viscosity ratio $\varLambda _\mu$, the thermal conductivity ratio $\varLambda _\kappa$, the dimensionless surface tension $\gamma$ and the reference pressure $\hat {p}_0$, for some commonly used liquids at the reference temperature $\hat {T}_0 = 300$ K. The viscosity and thermal conductivity ratios are with respect to the argon gas, which has the viscosity $\hat {\mu }_0 = 2.2725\times 10^{-5}$ Pa s at the reference temperature $\hat {T}_0 = 300$ K. The gas constant for argon is $\hat {R} = 208.1\ \textrm {J} (\textrm {kg}\ \textrm {K})^{-1}$. The parameters listed in table 3 have been determined using the data provided by the National Institute of Standards and Technology (2023). The reference pressure $\hat {p}_0$ is taken to be slightly above the saturation pressure of the fluid at the reference temperature $\hat {T}_{0}=300$ K.

Table 3. Parameters for some common engineering fluids at 300 K. The viscosity and thermal conductivity ratios are with respect to the argon gas.

Recall from § 4.1 that, for large viscosity ratios, the drag force on the liquid droplet approaches the drag force on a rigid sphere (the case of $\varLambda _\mu \to \infty$). As a matter of fact, Happel & Brenner's formula (4.9) also deduces that the percentage error between the drag force on the liquid droplet with respect to the drag force in the rigid sphere case is less than $10\,\%$ for all $\varLambda _\mu \geq 2.5$. For most common engineering fluids, the viscosity ratio $\varLambda _\mu$ is apparently bigger than $2.5$ (see table 3). Hence, for many practical applications (such as, a water droplet moving through air), the viscosity ratio is high enough to treat the liquid droplet to be a rigid sphere – as far as one is only concerned about the drag force. Notwithstanding, the coupled dynamics of the liquid droplet and gas leads to some interesting flow features in the temperature and heat profiles, as shown in §§ 4.2 and 4.3, that could potentially play important roles in processes involving phase change (Rana et al. 2018b). To emphasise the point, we plot in figure 14 the temperature profile in the case of argon gas flow past a rigid spherical ball made of glass, which has the thermal conductivity approximately twice of that of water, at ${Kn} = 0.3$ by solid red lines. The corresponding temperature profiles in the cases of argon gas flow past spherical liquid droplets of water (dashed blue lines), propyl alcohol (purple dot–dashed lines) and methanol (green solid lines) are also shown in figure 14. The vertical black line at $r=1$ demarcates the interface between the solid/liquid and gas. Evidently, there are significant differences in temperature profiles in all above cases. The differences are primarily due to their different thermal conductivity ratios and are more pronounced within the solid sphere or liquid droplet. Smaller heat conductivity ratio apparently leads to larger temperature gradients within and outside the liquid droplet or solid sphere.

Figure 14. Dimensionless deviations in the temperature (scaled with $\cos {\theta }$) as a function of position $r$ for ${Kn} = 0.3$ in the case of argon gas flow past a rigid spherical glass ball (solid red lines), in the cases of argon gas flow past spherical liquid droplets of water (dashed blue lines), propyl alcohol (purple dot–dashed lines) and methanol (green solid lines). The vertical black line at $r=1$ demarcates the interface between the solid/liquid and gas.

In order to highlight the effects of the internal motion on the temperature profiles, we plot in figure 15 the ratio of the temperature jump at the interface of a liquid droplet to the temperature jump at the interface of a corresponding solid sphere ($\varLambda _\mu \to \infty$) having the same thermal conductivity as that of the liquid against the Knudsen number. We choose the liquid droplets to be made of water (dashed blue lines in figure 15), propyl alcohol (dot–dashed purple lines in figure 15) and methanol (solid green lines in figure 15) so that our assumption of the liquid droplet being spherical holds good. The differences in the temperature jump ratios for different liquid droplets with respect to its solid counterpart (having the same thermal conductivity) are evident in figure 15, especially for relatively smaller Knudsen numbers, although the differences decrease with the increasing Knudsen number. Methanol having the viscosity ratio approximately $23$ shows the largest deviation (of approximately $30\,\%$).

Figure 15. The ratio of the temperature jump at the interface of a liquid droplet (made of water, propyl alcohol and methanol) to the temperature jump at the interface of a corresponding solid sphere ($\varLambda _\mu \to \infty$) having the same thermal conductivity as that of the liquid against plotted over the Knudsen number.

An important assumption made in this work is that the liquid droplet remains spherical throughout. This assumption, in other words, means that the surface tension force on the droplet is assumed to be larger than the pressure difference between inside and outside of the droplet at the interface, i.e.

(4.15)\begin{equation} \hat{p}^{(\ell)}-\hat{p} < \frac{2\hat{\gamma}}{\hat{R}_0}, \end{equation}

where $\hat {\gamma }$ is the surface tension and $\hat {R}_0$ is the radius of the droplet. Condition (4.15) in the dimensionless form reads

(4.16)\begin{equation} \frac{p^{(\ell)}-p}{Kn} < 2\gamma. \end{equation}

From figure 12, $(p^{(\ell )}-p)/{Kn}$ at the interface is approximately $4.44$ (for figure 12a,b), approximately $3.33$ (for figure 12c,d) and approximately $2$ (for figure 12e, f). Comparing these numbers with the twice of the dimensionless surface tension given in table 3, one can see that condition (4.16) holds true for water, propyl alcohol, methanol, ammonia and pentane but does not hold for isobutane, R134a and propane. Therefore, there are some liquids (e.g. water, propyl alcohol, methanol, etc.) for which our assumptions of the droplet being spherical is valid. Nevertheless, there are also many liquids for which the surface tension forces are not strong enough to maintain the liquid droplet spherical. For such liquids, it will be necessary to account for the effect of surface tension forces and to extend the present work.

5. Conclusion

In this article, we have studied the effects of internal motion within a spherical liquid droplet with its diameter being comparable to the mean free path of the surrounding (monatomic) gas under the assumptions of no phase change involved and the surface tension force being strong enough to maintain the liquid droplet spherical throughout. The rarefied gas phase in this article has been modelled with the LR26 equations while the flow in the liquid phase has been modelled with the incompressible NSF equations. Owing to our assumptions, the problem has become somewhat simplified – allowing for the analytic solution, which was still not so easy to obtain. The analytic solution plays a vital role in understanding the kinetic effects in rarefied gases and also helps in developing a deeper understanding on the effects, which the internal motion in a liquid droplet has on the flow pattern of the rarefied gas flowing past it. With the analytic solution, the effects of the liquid to gas viscosity ratio and the liquid to gas thermal conductivity ratio on the drag force and the overall flow dynamics have been investigated.

To validate the analytic solution obtained in the present work, the analytic results for the drag force obtained from the LR26 theory in the present work have been compared with the limited experimental measurements available in the literature. It turns out that the analytic results for the drag force obtained from the LR26 theory agree closely with experimental data even for significantly large values of the Knudsen number. After validating the analytic results on the drag force, the physical field variables, which are usually difficult to measure in experiments, have been presented. The analytic results obtained in the present work show a clear discontinuity in the velocities of the liquid and gas at the interface of the liquid droplet (figure 13). It is worthwhile noting that the viscosity ratio of the liquid to the gas in practical applications is usually large. Thus, from the findings of figure 2, it is acceptable in such applications to treat the liquid droplet as a rigid sphere of the same size if one is concerned only about the drag force. Nevertheless, § 4.4 shows that there could be significant differences in the other quantities (e.g. the temperature) when treating the liquid droplet as a rigid sphere. In addition, for some liquids and gases used in some other applications, the viscosity ratio could be close to unity even at standard conditions (at $25\,^\circ \textrm {C}$ and $1$ atmospheric pressure); for instance, the viscosity ratio between water (liquid) and ammonia (gas) at standard conditions is approximately 3.6. For such applications, the presented theory would yield meaningful results.

It is important to note that, owing to our assumptions, the density ratio does not appear in our analysis and hence does not affect the results in the present work. Nevertheless, the density ratio as well as the surface tension are the two most important parameters, whose effects are significant in problems of rarefied gas flow past liquid droplets and cannot be ignored in practical applications. The present work may be treated as a first step to achieving a good mathematical understanding of the (simplified) problem. More realistic problems, involving surface tension and/or phase change, will be subjects of future research.