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.


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 Λ µ -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 (Λ µ → ∞)] and (ii) liquid flow over a gas bubble [case of zero inside-to-outside viscosity ratio (Λ µ ≈ 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 R. Bhattacharjee, S. Saini, V. K. Gupta and A. S. Rana 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 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 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. 2019Rana et al. , 2021a,b;,b;Tiwari et al. 2021;De Fraja et al. 2022;Rana et al. 2018b).Some open-source software, like OpenFOAM, in combination with methods, such as the volume-of-fluid, level-set and the direct simulation Monte Carlo (DSMC) methods, have also been utilised to investigate gas-liquid multiphase flow problems (Malekzadeh & Roohi 2015;Chakraborty et al. 2019;Chakraborty 2019).Nevertheless, the problems studied with these 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 Λ µ ≈ 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 millimetre.However, for large drops (of diameters bigger than 1 millimetre), 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 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 ⩽ Λ µ < ∞ and for external Reynolds numbers 0.5 ⩽ Re ⩽ 100 that cover flow over a solid sphere, over a liquid drop and over a small gas bubble.For Re ≪ 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 Re ≪ 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 Re ⩽ 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) and 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 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 ⩽ Re ⩽ 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 < 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 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.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 et al. (2015); Stefanov et al. (2022); Taheri et al. (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 (BGK) model (Bhatnagar et al. 1954), the ellipsoidal statistical BGK (ES-BGK) 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 asymptoticexpansion 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 and 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 socalled 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. 2018bRana et al. , 2021a) ) but may yield quantitatively different results for certain processes beyond the Knudsen number unity; see, e.g.Rana et al. (2018bRana et al. ( , 2021a)).Despite this, the system comprised of the R26 equations (Gu & Emerson 2009) is the best known macroscopic model till 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), e.g., 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.

Problem formulation
We consider a slow steady uniform flow of a monatomic rarefied gas approaching from the negative ẑ-direction with a uniform velocity û∞ 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 (r, θ, ϕ), which is related to the Cartesian coordinate system (x, ŷ, ẑ) via (x, ŷ, ẑ) ≡ (r sin θ cos ϕ, r sin θ sin ϕ, r cos θ).Here, r ∈ [0, ∞), θ ∈ [0, π] and ϕ ∈ [0, 2π).The spherical symmetry of the droplet implies that the flow parameters are independent of the direction ϕ.Consequently, all the field variables pertaining to the problem are functions of r and θ only.
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 (Re ≪ 1), and to slow flows, i.e. to small Mach number flows (Ma ≪ 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 ≪ 1) and in for low Reynolds numbers (Re ≪ 1).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 ρ0 , a constant temperature T0 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.

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 ρ, velocity vi , temperature T , stress tensor σij , heat flux qi , (tracefree) third velocity moment mijk , partially contracted (tracefree) fourth velocity moment Rij and fully contracted fourth velocity moment (also referred to as the scalar fourth moment) ∆.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: where R is the gas constant; p0 = ρ0 R T0 is the pressure in the reference state; ρ, v i , T , p, σ 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 ∆ are the dimensionless deviations in the corresponding quantities.In addition, the droplet radius R0 is taken as the length scale for making the space variable r dimensionless, i.e. r := r/ R0 .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) with the additional unknowns in (2.7)-(2.9)being ) ∂R rr ∂θ . (2.12b) 2)-(2.12), henceforth, will be referred to as the linearised R26 (LR26) equations.The coefficients Kn, Pr, Pr m , Pr R , Pr ∆ , Pr Φ , Pr ψ , Pr Ω in the LR26 equations are the dimensionless numbers arising from the nondimensionalisation of the equations.In particular, the numbers Kn = μ ρ0 R T0 R0 and Pr = 5 2 μ κ R (2.13) are referred to as the Knudsen number and the Prandtl number, respectively, with μ being the viscosity of the gas and κ being the thermal conductivity of the gas.It should be noted that, owing to the linearisation, the viscosity μ and thermal conductivity κ in (2.13) are the viscosity and thermal conductivity of the gas at the reference state temperature T0 ; and hence both are constant.Let us denote the viscosity and thermal conductivity of the gas at the reference state temperature T0 by μ0 and κ0 , respectively.Thus, owing to the linearisation, μ = μ0 and κ = κ0 throughout this work.The values of the numbers Pr, Pr m , Pr R , Pr ∆ , Pr Φ , Pr ψ , Pr Ω 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 ∆ = 2/3, Pr Φ = 2.097, Pr ψ = 1.698,Pr Ω = 1 (Gu & Emerson 2009).The subscripts r and θ 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 σ rθ is the rθ-component of the deviation in the stress tensor.Equation (2.2) can be identified as the equation of continuity for the gas, equations (2.3a) and (2.3b) as the momentum balance equations in the r-and θ-directions, respectively, and equation (2.4) as the energy balance equation, equations (2.5a) and (2.5b) as the balance equations for the rrand rθ-components of the stress, equations (2.6a) and (2.6b) as the heat flux balance equations in the r-and θ-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.15)

Modelling of the liquid phase inside the droplet
The liquid phase inside the spherical droplet is modelled with the linearised, steadystate, 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 ρ0 and temperature T0 in the reference state and the droplet radius R0 as the length scale.After simplification, the linear, R. Bhattacharjee, S. Saini, V. K. Gupta and A. S. Rana dimensionless steady-state NSF equations for modelling the liquid phase of the problem under consideration read (2.20) The superscript '(ℓ)' in (2.16)-(2.20)has been used to indicate that the variables with the superscript '(ℓ)' belong to the liquid phase (i.e. to the liquid droplet).The variables in (2.16)-(2.20)are as follows.The variables are the dimensionless deviations in the r-and θ-components of the velocity of the liquid, respectively; is the dimensionless deviation in the pressure of the liquid droplet from its pressure in the equilibrium state p(ℓ) 0 , which is given by the pressure inside a stationary droplet in a quiescent environment, i.e. p(ℓ) 0 = p0 + 2γ/ R0 .Note that the reference pressure inside a stationary droplet is higher than the reference pressure outside the droplet because of the surface tension γ.Furthermore, in (2.16)-(2.20), are the dimensionless deviations in the rr-and rθ-components of the stress tensor for the liquid, respectively; are the dimensionless deviations in the r-and θ-components of the heat flux of the liquid, respectively; is the dimensionless deviation in the temperature of the liquid droplet from the temperature in the reference state T0 ; Λ µ = μ(ℓ) /μ is the ratio of the viscosity of the liquid to the viscosity of the gas and Λ κ = κ(ℓ) /κ 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, μ(ℓ) and κ(ℓ) are the viscosity and thermal conductivity of the liquid at the reference temperature T0 ; consequently, μ, μ(ℓ) , κ, κ(ℓ) , Λ µ and Λ κ are all constant in the present work.The field variables with hats and superscript '(ℓ)' can be identified as the field variables of the original NSF equations.Furthermore, equation (2.16) can be identified as the equation of continuity, equations (2.17a) and (2.17b) as the momentum balance equations in the r-and θ-directions, respectively, and equation (2.18) as the energy balance equation.

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. (2021aRana et al. ( , 2018a) ) 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 ϑ = 0 and the normal component of the flow velocity relative to the interface velocity 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 θ-and ϕ-directions are the two tangential directions; nevertheless, the present problem is independent of ϕ due to spherical symmetry.Therefore, we shall replace the subscripts n with r and i with θ 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) where v I r is the velocity of the interface and is zero for the present problem, χ is the accommodation coefficient is the velocity slip and T = T − T (ℓ) 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 ϕ.
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.28) 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. v (ℓ)  r = 0 at r = 1. (2.29) (ii) The heat flux and shear stress are continuous at the interface; this implies that 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.19) and (2.20) 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 (2.30) 2 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) where p1 and p2 are the pressures exerted on the interface by the fluids '1' and '2', respectively; γc is the surface tension coefficient; κ is the local surface curvature; n i is the unit normal at the interface; and σ(1) ij and σ(2) ij 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 too; in other words, the radius of the droplet R0 is assumed to remain constant.Therefore, we neglect the nonequilibrium force balance in the normal direction as described in equation (2.31).It is justified since the equilibrium aspect has already been considered when defining p(ℓ) above.Additionally, for the force balance in the tangential direction, we multiply equation (2.31) by the unit tangent vector t i at the interface.This multiplication makes the lefthand side of the resulting equation vanish.By disregarding the effect of surface tension forces, specifically the term ∂γ c /∂ xi (commonly referred to as the Marangoni surface tension gradient) in equation (2.31), the boundary condition simplifies to σ(1) nt = σ(2) nt , which in the dimensionless form is the same as boundary condition (2.30) 2 considered in the present work.

Analytic solution methodology
We solve the equations derived in § 2 analytically by following a method proposed in (Torrilhon 2010;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 θ 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 θ indices are taken to be proportional to cos θ, the vectorial/tensorial components having an odd number of θ indices are taken to be proportional to sin θ 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 where the functions ℓ) and T (ℓ)   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.20).After simplification (cos θ and sin θ 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 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.3), 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.
R. Bhattacharjee, S. Saini, V. K. Gupta and A. S. Rana 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 , . . ., 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.

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.

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) where F Stokes = 6π Kn u ∞ is the Stokes drag with u ∞ = û∞ / R T0 being the farfield dimensionless velocity of gas approaching the droplet (for simplicity, we have taken u ∞ = 1 in this article), and a, b and 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) To compare the drag force computed analytically in the present work, we use the data resulting from formula (4.1) with the constants a, b and c given in (4.2) [taken from the textbook Kennard (1938)] and in (4.3) [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 a, b and c given in (4.4) [determined by Allen & Raabe (1985)] and in (4.5) [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 oildrop 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) where ⃗ z x (θ) = (cos θ, − sin θ, 0), ⃗ k = (1, 0, 0) and P = pI + σ is the pressure tensor.The force on the surface of the liquid droplet is represented by P • ⃗ k and its scalar product with ⃗ z x gives the component of the force in the direction of the flow.Substituting ansatz (3.1) in (4.6), we obtain This, on simplifying further, yields 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 Λ κ = 100 and for different values of the viscosity ratio Λ µ .In addition, the explicit values of the normalised drag force for the same value of the thermal conductivity ratio Λ κ = 100 and for different values of the Knudsen number and viscosity ratio are also given in table 1.
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 nonzero 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 Λ µ , 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 Λ κ = 100 and for different values of the viscosity ratio Λ µ .In addition, the explicit values of the normalised drag force for the same value of the thermal conductivity ratio Λ κ = 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 R. Bhattacharjee, S. Saini, V. K. Gupta and A. S. Rana equations in the present work with the drag force computed with an explicit formula given by Happel & Brenner (1965), which-in our notations-reads 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 Λ µ , 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 Λ µ → ∞, F HB → 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.

Flow fields: temperature and heat flux
Figures 4, 5 and 6 illustrate the heat flux lines and temperature contours in the ŷ = 0 plane for a fixed viscosity ratio Λ µ = 100 and for the thermal conductivity ratios Λ κ = 1, 10 and 100, respectively.The panels in top, middle and bottom rows of each figure display the results for the Knudsen numbers Kn = 0.09, 0.36 and 0.9, respectively.The panels in the left columns of figures 4, 5 and 6 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 those on the right columns of figures 4, 5 and 6 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 4, 5 and 6, 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 the right column of figures 4, 5 and 6 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 the left and right columns of figures 4, 5 and 6).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.
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 θ to understand the results for any angle θ) with respect to the position r in figures 7 and 8, respectively, again for a fixed viscosity ratio Λ µ = 100 and for the thermal conductivity ratios Λ κ = 1, 10 and 100.The panels in top, middle and bottom rows in both figures display the results for the Knudsen numbers Kn = 0.09, 0.36 and 0.9, respectively.The panels in the left column of each of figures 7 and 8 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 those in the right column of each of these figures 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   To see this, let us, for example, fix θ = 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 the left and right panels in figure 7, one finds that the temperature of the external gas predicted by the LR26 equations (left panels) increases on moving away from the interface whereas that predicted by the linear NSF equations (left panels) decreases on moving away from the interface.These are in agreement with the temperature contours portrayed in figure 4 (notice the contours for θ = 0 when moving away from the centre of the droplet).Furthermore, each panel of figure 7 clearly shows a nonzero temperature , Kn = 0.36 (middle row) and Kn = 0.9 (bottom row).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 for the panels in the left column and with the NSF equations for the panels in the right column.
jump T = T − T (ℓ) at the interface.From the panels on the right column of figure 7, 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 T = T − T (ℓ)  predicted by the NSF equations is always positive.This is not the case when the LR26 equations are used for the external flow (see panels on the left column of figure 7).In fact, when the LR26 equations are used for the external flow, the temperature jump T = T − T (ℓ) 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 θ) 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 each panel of 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 each panel of figure 8), and is decreasing monotonically for r > 1. , Kn = 0.36 (middle row) and Kn = 0.9 (bottom row).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 for the panels in the left column and with the NSF equations for the panels in the right column.

Flow fields: pressure and velocity
Figures 9, 10 and 11 illustrate the velocity streamlines plotted over the pressure contours in the ŷ = 0 plane for a fixed thermal conductivity ratio Λ κ = 100 and for the viscosity ratios Λ µ = 1, 10 and 100, respectively.The panels in top, middle and bottom rows of each figure again denote the results for the Knudsen numbers Kn = 0.09, 0.36 and 0.9, respectively.The panels in the left columns of figures 9, 10 and 11 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 those on the right columns of figures 9, 10 and 11 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 ẑ-direction, the liquid inside the droplet near the top [close to (r, θ, ϕ) = (1, π/2, 0)] and bottom [close to (r, θ, ϕ) = (1, π/2, π)] 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.Furthermore, a comparison of the corresponding panels on the left and right columns of each of figures 9, 10 and 11 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 (panels in the left columns) and the linear NSF equations (panels in the right columns) 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 θ to understand the results for any angle θ) with respect to the position r in figure 12 and the z-component of the velocity for a fixed θ = π/2 with respect to the position r in figure 13.For both figures 12 and 13, the thermal conductivity ratio is again fixed to Λ κ = 100 and the viscosity ratios are taken as Λ µ = 1, 10 and 100.The panels in top, middle and bottom rows in both figures again exhibit the results for the Knudsen numbers Kn = 0.09, 0.36 and 0.9, respectively.The panels in the left column of each of figures 12 and 13 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 those in the right column of each of these ) and Kn = 0.9 (bottom row).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 for the panels in the left column and with the NSF equations for the panels in the right column.
figures 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 all panels of 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 (3.2) 1,2,3 and (3.3) 1,2,3 , and they read Therefore, inside the droplet, p (ℓ) / cos θ is a linear function of r while v z is parabolic in r for θ = π/2.Consequently, the curves on the left of the vertical black lines in each panel 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 for the panels in the left column and with the NSF equations for the panels in the right column.
of figure 12 are straight lines and the curves on the left of the vertical black lines in each panel of 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. is the dimensionless surface tension.The dimensionless deviation in the liquid pressure p (ℓ) can be made as small as desired in comparison with 1 + 2γKn by taking û∞ to be sufficiently small.This keeps our assumptions of the shape of the droplet being spherical and its size being fixed justifiable.
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 ŷ = 0 plane while figure 13 displays the results in the ẑ = 0 (or θ = π/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 ∞ = 1 as r approaches ∞.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 Λ µ becomes larger and larger, which is consistent with the case of a rarefied gas flow over a solid sphere.

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 Λ µ and the thermal conductivity ratio Λ κ .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 Λ µ , the thermal conductivity ratio Λ κ , the dimensionless surface tension γ and the reference pressure p0 , for some commonly used liquids at the reference temperature T0 = 300 K.The viscosity and thermal conductivity ratios are with respect to the argon gas, which has the viscosity μ0 = 2.2725 × 10 −5 Pa s at the reference temperature T0 = 300 K.The gas constant for argon is R = 208.1 J/(Kg K).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 p0 is taken to be slightly above the saturation pressure of the fluid at the reference temperature T0 = 300 K.
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 Λ µ → ∞).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 Λ µ ⩾ 2.5.For most common engineering fluids, the viscosity ratio Λ µ 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 about 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.
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 (Λ µ → ∞) 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 about 23 shows the largest deviation (of about 30%).
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.From figure 12, (p (ℓ) − p)/Kn at the interface is about 4.44 (for the top row), about 3.33 (for the middle row) and about 2 (for the bottom row).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.

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°C and 1 atmospheric pressure); for instance, the viscosity ratio between water (liquid) and ammonia (gas) at standard conditions is about 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.
Declaration of interests.The authors report no conflict of interest.
Appendix A. Analytic expressions for the unknowns in ansatzes (3.1) and (3.2) The analytic solution for the unknowns appearing in (3.1) and (3.2) are as follows.

Figure 1 .
Figure 1.Schematic of a rarefied gas flow past a spherical liquid droplet.
.23, b = 0.41 and c = 0.88.(4.2)The most accurate raw data from Millikan's experiment were reviewed later byAllen & 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 a, b and c obtained byAllen & Raabe (1982) are a = 1.155 ± 0.008, b = 0.471 ± 0.011 and c = 0.596 ± 0.050.(4.3)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 a, b and 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 a, b and c to be a = 1.142 ± 0.0024, b = 0.558 ± 0.0024 and c = 0.999 ± 0.0212.(4.4)Using the modulated dynamic light scattering technique-a technique fundamentally different from the ones used by Millikan and by Allen and Raabe, Hutchins et al. (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 a, b and c to be a = 1.2310 ± 0.0022, b = 0.4695 ± 0.0037 and c = 1.1783 ± 0.0091.(4.5)

Figure 2 .
Figure 2. Drag force (normalised with the Stokes drag F Stokes = 6π Kn u∞) on a liquid droplet from the linearised NSF and R26 theories as a function of the Knudsen number for a fixed thermal conductivity ratio Λκ = 100.The square and star symbols denote the experimental data from Millikan's oil-drop experiment fitted by the empirical formulae ofKennard (1938) [formula (4.1) with (4.2)] andAllen & Raabe (1982) [formula (4.1) with (4.3)], respectively.The disk and diamond symbols denote the data from experiments performed with solid spherical particles in air fitted by the empirical formulae ofAllen & Raabe (1985) [formula (4.1) with (4.4)] andHutchins et al. (1995) [formula (4.1) with (4.5)], respectively, and are included just for comparison.

Figure 3 .
Figure3.Drag force on the liquid droplet normalised with the drag force given by (4.9).The thermal conductivity ratio is Λκ = 100.

Figure 4 .
Figure 4. Heat flux lines plotted on top of the temperature contours for the viscosity ratio Λµ = 100 and for different values of the Knudsen number: Kn = 0.09 (top row), Kn = 0.36 (middle row) and Kn = 0.9 (bottom row).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 for the panels in the left column and with the linear NSF equations for the panels in the right column.The thermal conductivity ratio is Λκ = 1.

Figure 5 .
Figure 5. Same as figure 4 but for the thermal conductivity ratio Λκ = 10.

Figure 6 .
Figure 6.Same as figure 4 but for the thermal conductivity ratio Λκ = 100.

Figure 7 .
Figure7.Dimensionless deviations in the temperature (scaled with cos θ) as a function of position r for a fixed viscosity ratio Λµ = 100 and for different values of the Knudsen number: Kn = 0.09 (top row), Kn = 0.36 (middle row) and Kn = 0.9 (bottom row).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 for the panels in the left column and with the NSF equations for the panels in the right column.

Figure 8 .
Figure8.Dimensionless radial heat flux (scaled with cos θ) as a function of position r for a fixed viscosity ratio Λµ = 100 and for different values of the Knudsen number: Kn = 0.09 (top row), Kn = 0.36 (middle row) and Kn = 0.9 (bottom row).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 for the panels in the left column and with the NSF equations for the panels in the right column.

Figure 9 .
Figure9.Velocity streamlines plotted over the pressure contours for the thermal conductivity ratio Λκ = 100 and for different values of the Knudsen number: Kn = 0.09 (top row), Kn = 0.36 (middle row) and 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 Λµ = 1.

Figure 12 .
Figure12.Dimensionless deviations in the pressure (scaled with cos θ) as a function of position r for a fixed thermal conductivity ratio Λκ = 100 and for different values of the Knudsen number: Kn = 0.09 (top row), Kn = 0.36 (middle row) and Kn = 0.9 (bottom row).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 for the panels in the left column and with the NSF equations for the panels in the right column.

Figure 13 .
Figure13.The z-component of the (dimensionless) velocity as a function of position r in the plane θ = π/2 for a fixed thermal conductivity ratio Λκ = 100 and for different values of the Knudsen number: Kn = 0.09 (top row), Kn = 0.36 (middle row) and Kn = 0.9 (bottom row).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 for the panels in the left column and with the NSF equations for the panels in the right column.

Figure 14 .Figure 15 .
Figure14.Dimensionless deviations in the temperature (scaled with cos θ) 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.
the surface tension and R0 is the radius of the droplet.Condition (4.15) in the dimensionless form reads p (ℓ) − p Kn < 2γ.(4.16)

Table 1 .
Drag force normalised with the Stokes drag, F Stokes = 6π Kn u∞, for different values of the Knudsen number and viscosity ratios at a fixed thermal conductivity ratio Λκ = 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.

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