Evaporation-driven vapour microflows: analytical solutions from moment methods

Macroscopic models based on moment equations are developed to describe the transport of mass and energy near the phase boundary between a liquid and its rarefied vapour due to evaporation and hence, in this study, condensation. For evaporation from a spherical droplet, analytic solutions are obtained to the linearised equations from the Navier–Stokes–Fourier, regularised 13-moment and regularised 26-moment frameworks. Results are shown to approach computational solutions to the Boltzmann equation as the number of moments are increased, with good agreement for Knudsen number ${\lesssim}1$ , whilst providing clear insight into non-equilibrium phenomena occurring adjacent to the interface.


Introduction
The study of two phase systems is both of practical significance in engineering and fundamental interest in the pure sciences. For example, continued miniaturisation of mechanical and electrical components has led to steep increases in power density, making thermal management a critical aspect of microelectronics and microelectromechanical systems (MEMS) design today. Evaporation and condensation are very efficient modes of heat transfer and as such are employed in component cooling and in energy conversion systems (Ling et al. 2014;Plawsky et al. 2014) leading to considerable research into the study of phase transition at the micro and nano scales. Understanding the evaporation of liquid droplets is critical to a broad range of natural processes, such as body temperature control of mammals (Sherwood 2005), and many industrial cooling processes, such as spray drying, spray cooling and microelectronics cooling (Dhavaleswarapu, Murthy & Garimella 2012;Plawsky et al. 2014;Hodes et al. 2015).
One of the main difficulties in the description of phase-change microflows is that the classical Navier-Stokes-Fourier (NSF) model can fail to accurately capture the 964 A. S. Rana, D. A. Lockerby and J. E. Sprittles (Ma 1.65). Gu & Emerson (2007) presented a set of wall boundary conditions for the R13 equations derived from the Maxwell accommodation model. However, in their work unphysical oscillations were reported near the boundary due to the over-prescription of the boundary conditions, as identified and rectified by . Thus far, the R13 equations have been considered for canonical boundary-value problems, such as planar and cylindrical Couette and Poiseuille flows ), transpiration flows and gas flow past a sphere (Torrilhon 2010), among others, in one-and two-dimensional numerical simulations.
To extend the applicability of the moment method further into the transition regime, Gu & Emerson (2009) used the regularisation process to obtain the regularised 26moment (R26) equations. The kinetic wall boundary conditions for the R26 equations were also derived in the same paper using the Maxwell accommodation model. The R26 theory describes Knudsen layers (regions of strongly non-equilibrium behaviour within a few mean free paths of the boundary) more accurately than the R13 equations (Gu & Emerson 2014); the NSF and G13 equations are not able to predict these at all. The R13/R26 systems have produced accurate results for small to moderate Kn and have the advantage that they offer fast, sometimes analytic, solutions to enable us to understand, and develop intuition for the general flow behaviour. In this article, we shall consider whether similar results can be derived for phase-change phenomena.

Modelling evaporation and condensation
Flows involving evaporation/condensation processes are usually modelled using a sharp interface approach, where the interface separating the vapour and liquid phases is assumed to be infinitely thin. The phase transition process itself is modelled by setting the appropriate boundary conditions on the interface. Different variants of these boundary conditions have been employed in the literature. The standard assumptions, found in numerous classical textbooks and incorporated into most computational fluid dynamics software, is that the temperature and the velocity tangential to the interface are continuous across it. However, it is now well established, both experimentally and theoretically, that such assumptions are invalid at small length scales, where a velocity slip and temperature jump are observed (McGaughey & Ward 2002;Rahimi & Ward 2005). The first mathematical model for the evaporation and condensation processes was developed by Maxwell (1867) for quasi-steady flow, where the flow at the interface can be assumed equal to the diffusion vapour flow far away from the interface. However, the classical diffusion flow of vapour starts beyond the Knudsen layer. The Hertz-Knudsen-Schrage (HKS) relation (Schrage 1953) can be applied to model the mass flux across the interface, along with the NSF equations in the bulk, provided Kn is sufficiently small (Persad & Ward 2016). Fuchs (1959) obtained a semi-empirical formula to account for the Knudsen layer by matching the free molecular solution, valid at the interface, with the diffusion solution valid at the edge of the Knudsen layer, the width of which is used as a fitting parameter to give a good match for all Kn. Young (1991) furthered this approach by deriving expressions for the mass and energy fluxes across a fluid droplet in a pure vapour assuming the G13 distribution inside the Knudsen layer. However, Young's approach also involves an empirical parameter (related to thickness of the Knudsen layer) and does not resolve the Knudsen layer properly. Bond & Struchtrup (2004) considered interface conditions using the kinetic theory of gases and irreversible thermodynamics (using NSF equations) to develop predictive Evaporation-driven vapour microflows 965 expressions for the mass and energy fluxes across the interface and compared the results between irreversible thermodynamics, statistical rate theory and experimental data. In particular, they considered a condensation coefficient which depends on the impact energy of the condensing particle and also pointed out the effect of interfacial curvature on temperature jump. More recently, steps towards using higher-order moment methods for mass and energy transport across the phase boundary were made by Struchtrup & Frezzotti (2016), with a full set of boundary conditions for the R13 equations derived and solved for a planar case. In this article, we extend this work to both curved interfaces, using spherical drops, and also derive phase-change boundary conditions for the R26 moment equations. This will allow us to compare the NSF, R13 and R26 models for evaporation and condensation for both a spherical droplet and a planar interface, in order to establish their accuracy against benchmark solutions from the Boltzmann equation. In this article, due to the linearity of the equations and boundary conditions considered, the results obtained for the evaporation process can simply be translated to condensation by reversing the signs of thermodynamics forces.
The article is organised as follows. A brief summary of basic elements of kinetic theory and the moment method is given in § 2, followed by the derivation of moment equations, namely the NSF, R13 and R26 equations for a spherical geometry, in § 3. In § 4, analytical solutions to these equations are found. In § 5, boundary conditions for the R26 equations are derived and applied for the evaporation problem. Throughout § §6-8, Boltzmann solutions provide the benchmark for our analytic results; beginning, in § 6, by considering Kn → 0. In § 9, we analyse the magnitude of competing higherorder contributions to the total normal stress. Finally, in § 10 we conclude and discuss future directions of research.

Moment methods in kinetic theory
Sufficiently far from the critical point, the vapour can be treated as an ideal gas and one can use the kinetic theory of dilute gases for this phase. A monoatomic ideal gas can be characterised by a one-particle distribution function f which depends on time t, spatial coordinates x i and molecular velocity c i . The distribution function is governed by the Boltzmann equation (Kogan 1969;Cercignani 1975) where G k is the external force per unit mass acting on the gas and is assumed to be independent of the microscopic velocity c k . The term S is the collision operator that describes the change of the distribution function due to interaction between particles. For most engineering processes, the main interest is not in detailed information about the distribution function f , but in the macroscopic quantities, such as mass density , macroscopic velocity v i , and temperature T. These are defined as moments of the distribution function: where C i = c i − v i is the peculiar velocity, R is the gas constant and dc = dc 1 dc 2 dc 3 . The temperature is related to the pressure ℘ and density by the ideal-gas law ℘ = RT. The pressure tensor p ij , and the heat-flux vector q i are the second and contracted third-order moments of the distribution function f , respectively, i.e.
Furthermore, the pressure tensor can be separated as p ij = ℘δ ij + σ ij , where δ ij is the Kronecker delta, σ ij = p ij is the deviatoric stress tensor, and ℘ = p kk /3 is the pressure. The angular brackets are used to denote the traceless part of a symmetric tensor. Instead of attempting a direct solution of the Boltzmann equation (2.1), the moment method provides a bridge between kinetic theory and classical hydrodynamics via evolution equations, known as moment equations, for the macroscopic moments. A set of moment equations is obtained by multiplication of the Boltzmann equation (2.1) with an arbitrary function Ψ A and subsequent integration over the microscopic velocity space. For example, by choosing Ψ A = 1, c i and C 2 /2, we get the conservation laws for mass, momentum and energy densities, respectively. However, in pronounced non-equilibrium situations, it is necessary to extend the set of macroscopic variables beyond the 5 conservative variables, so as to include higher-order moments. For instance, in the 13-moment approximation Ψ A = 1, c i , C i C j and C 2 C i /2, corresponding to the 13 moments (5 conservative variables, 5 components of σ ij and 3 components of q i ). A further extension of the moment equations contain even higher-order moments (see below) and here we follow the same notation as used in Gu & Emerson (2009) for the R26 equations, where ∆, Ω i , R ij , m ijk , ψ ijk and φ ijkl are scalar, first-, second-, third-and fourth-order symmetric trace-free tensors, respectively, defined as Here, the higher moments have been constructed in such a way that they vanish in the 13-moment theory. The governing equations for the moments can be readily obtained from the Boltzmann equation (2.1), see e.g. Struchtrup (2005); however, they cannot be solved as they stand, since they do not form a closed set of equations. The regularisation process presented by Struchtrup & Torrilhon (2003) and Struchtrup (2004) provides a framework to derive a closed set of moment equations, giving the R13 equations at third-order accuracy while at the fifth order we arrive at the Evaporation-driven vapour microflows 967 Vapour Far-field Droplet r FIGURE 1. Schematic representation of a liquid droplet surrounded by its own vapour. T l is the temperature of liquid at the interface which corresponds to the saturation pressure ℘ l . The temperature and pressure of the vapour at distance far from the surface of the droplet are given as T ∞ and ℘ ∞ , respectively. R26 equations, see Struchtrup (2005) and Gu & Emerson (2009). The derivation of the R13 and R26 moment equations for a monatomic gas of Maxwell molecules (i.e. molecules repelling each other with an intermolecular force ∝ r −5 ) and their linearised version can be found in , Gu et al. (2010) and Gu & Emerson (2014). The equations of these models will now be presented for the spherically symmetric case.

Problem statement
Consider a liquid droplet of fixed radius R 0 at a given temperature T l , with corresponding saturation pressure in the liquid ℘ l , immersed in its own vapour; see figure 1. This problem has been studied fairly extensively, e.g. Sone & Onishi (1978), Chernyak & Margilevskiy (1989), Takata et al. (1998), which allows us to test the accuracy of the derived evaporation/condensation boundary conditions for the extended moment equations. Let the temperature and pressure of the vapour at a distance far from the surface of the droplet be given by T ∞ and ℘ ∞ , respectively. In general, the saturation pressure ℘ l is a function of T l given by the Clausius-Clapeyron relation; however, when the droplet has fixed properties this is not invoked so that ℘ l and T l can be varied independently. A spherical coordinate system (r, θ , ϕ) is considered, and because of the spherical symmetry, the flow is independent of the azimuthal and polar directions. Two cases will be considered.
(i) Pressure-driven flow. The process in the vapour is driven by a dimensionless pressure difference p l = (℘ l − ℘ ∞ )/℘ ∞ while the temperature of the liquid is equal to the far-field temperature. (ii) Temperature-driven flow. The process is driven by a dimensionless temperature difference θ l = (T l − T ∞ )/T ∞ while the saturation pressure in the liquid is same as the far-field pressure.
In general, any combination of these two cases can be superimposed, owing to the linearity of the equations and boundary conditions considered. As a result of the flow configuration, the velocity vector v, the heat-flux vector q and stress tensor σ simplify to v ≡ {v(r), 0, 0} , q ≡ {q(r), 0, 0} and σ ≡ Furthermore, the only non-zero components of the higher-order tensors Ω i , R ij , m ijk , ψ ijk and φ ijkl are with the remaining components either zero or following from the symmetry and tracefree conditions.

Relation to the unsteady evaporation of a drop
The steady process considered in this article will form the basis of future work into the unsteady problem, where one is interested in how long it takes a drop to completely evaporate. In this case, the temperature distribution inside the liquid drop T l (t, r) must be calculated as part of the solution, as must the position of the drop's radius R(t), see e.g. the molecular dynamics study by Holyst & Litniewski (2008) and the experimental and theoretical study by Holyst et al. (2013). The pressure inside the liquid is also now part of the solution with the saturation pressure ℘ l a function of T l and R given by the Clausius-Clapeyron relation with corrections (Kelvin equation) that account for the curvature of the liquid-vapour interface (Young 1991). The problem can be significantly simplified owing to the high density ratio between the liquid and the vapour phase, which creates a quasi-steady process in the vapour phase.
To determine the radius of the drop and the temperature distribution in the liquid the mass flux j and heat flux q into the vapour must be calculated. For small deviations from equilibrium, i.e. p l 1 and θ l 1, the pressure-driven and temperature-driven cases can be combined to give the mass flux j and the heat flux q as j = j p p l + j τ θ l , (2.10) q = q p p l + q τ θ l , (2.11) where, j p (q p ) and j τ (q τ ) are the mass (heat) flux for the pressure-driven and the temperature-driven cases, respectively, and are derived in this article.

Extended moment equations in spherical geometry
Let R 0 , ℘ ∞ , T ∞ be the reference length, pressure and temperature, respectively. The equations are made dimensionless by introducinĝ whereθ andp are the dimensionless deviation of the temperature and pressure from their far-field values, respectively. The Knudsen number is defined as where µ ∞ is the gas viscosity coefficient at the reference state. The mean free path, λ, is related to the usual definition of the mean free path 0 (Sone 2007) by Here, A = 1.27 for the hard-sphere model and A = 1 for the BGK model, which gives the same fluid viscosity for different models and allows for a consistent comparison. For brevity, the hat will be removed hereafter, and unless otherwise stated, all variables will be given in dimensionless form.
For the linearised equations, only terms that are linear in deviations from the reference equilibrium state are considered. Accordingly, dimensionless and linearised conservation laws for mass, momentum and energy, in spherically symmetric coordinates, read which contain the normal component of the viscous stress, σ and the radial heat flux, q as unknowns. The various theories for transport provide equations for stress and heat flux which are presented in the following subsections.

Linearised NSF equations
In the dimensionless form and one-dimensional spherical geometry, the NSF constitutive laws for stress and heat flux are given by Here, the Prandtl number Pr = 2/3 for Maxwell molecules (MM) model and Pr = 1 for the BGK model. The linearised form of the mass and momentum conservation laws (3.4) along with the constitutive relation for stress tensor (3.5) give the Stokes equations. In this case, the flow problem decouples from the Fourier-based energy problem, so that non-equilibrium cross-effects, such as thermal stress and non-Fourier heat flux (Sone 2007;Rana et al. 2013), cannot be predicted by the NSF system. High-order macroscopic models based on moment equations are developed to describe these effects.

Higher-order moment equations
The balance equations for the heat-flux vector and stress tensor follow from the Boltzmann equation as 8r 15 Kn Pr ∂θ ∂r . (3.7b) The right-hand side terms of these balance equations contain NSF constitutive laws (3.5-3.6), which can also be obtained by a first-order Chapman-Enskog expansion (Chapman & Cowling 1970) of the stress and heat-flux balance equations (3.7).

R13 constitutive relations
The balance equations (3.7) do not form a closed set, since they contain additional higher moments m, R and ∆. In G13 theory, m, R and ∆ all are set to zero. A non-vanishing closure for these higher-order moments is given by the R13 constitutive relations (Struchtrup 2005). The linear R13 constitutive relations, for the considered geometry, take the following form where transport coefficients Pr M , Pr R and Pr ∆ depend upon the choice of intermolecular potential function appearing in Boltzmann collision operator. For the MM model they are given by Pr = 2/3, Pr M = 3/2, Pr R = 7/6, Pr ∆ = 2/3 and all are 1 for BGK approximations (Truesdell & Muncaster 1980;Gu & Emerson 2009). The conservation laws (3.4), the balance equations for heat-flux vector and stress tensor (3.7) along with the constitutive relations (3.8) form the R13 equations.

R26 constitutive relations
The 26-moment theory consists of the conservation laws, the balance equations for the heat flux and stress and the balance equations for the higher moments, m ijk , R ij and ∆. The evolution equations for m ijk , R ij and ∆ were obtained from the Boltzmann equation in Gu & Emerson (2009), which for the present problem read 9r 2 35 As for the NSF and the R13 cases, constitutive relations are needed to express the fluxes Φ, ψ and ω, to form a closed system. The R26 equations include the constitutive relations for these higher-order terms as (Gu & Emerson 2009 where more transport coefficients appear as Pr Φ = 2.097 and Pr ψ = 1.698 for MM and Pr Φ = Pr ψ = 1 for the BGK model. The constitutive relations (3.10) with balance equations (3.9) and (3.7) and conservation laws (3.4) form the system of the R26 equations.
Comparison of the R13 constitutive relations (3.8) with the balance equations (3.9) reveals that the R13 constitutive relations (3.8) stem from a Chapman-Enskog expansion of (3.9) in Kn, i.e. taking the right-hand side of (3.9) to be zero. Likewise, the R26 constitutive relations (3.10) follow from the Chapman-Enskog expansion of the evolution equations for the moments φ ijkl , ψ ijk and ω i from the 45-moment theory.
In the next section, we shall derive analytic solutions to the NSF, R13 and R26 equations for the spherical geometry and then in § 5 formulate and apply the boundary conditions required for the evaporation problem.

Analytical solutions of the moment equations
Due to the linearity of the equations and the simple geometry, we have found that all equations can be solved analytically. By integrating the mass and energy conservation laws (3.4a,c), we immediately find v = c 1 r 2 and q = where c 1 and c 2 are integration coefficients. These solutions are obtained from the conservation laws without any assumption on constitutive relations and therefore they are valid for all values of the Knudsen number. The integration coefficients c 1 and c 2 , however, depend on other field variables through boundary conditions, and consequently, these depend on the constitutive relations. Therefore, for the given flow configurations, the solutions of the NSF equations include two integration constants, c 1 and c 2 . These constants need to be fixed from two boundary conditions at the interface at r = 1.

Solution to the R13 equations
The integration of the heat-flux and stress balance equations (3.7) and use of the R13 constitutive relations (3.8) with the solutions (4.1), give where β 1 = √ 5Pr M /3. In (4.5) we have already used the far-field conditions (4.4). Similarly, the momentum equations (3.4b) and (4.5a) give the pressure The R13 equations therefore allow for a non-uniform pressure in the gas that decays exponentially with the scaled distance rβ 1 from the interface. It is a purely non-equilibrium effect caused by the Knudsen layer. The non-uniform pressure has been observed in molecular dynamics ( Solutions (4.5)-(4.6) of the R13 equations have three integration constants, c 1,2 and a 1 , which need to be determined from the three boundary conditions at r = 1.

Solution to the R26 equations
The integration of the R26 equations is more involved and can be found in the appendix; here we present only the final results. The pressure and normal stress are given as the superposition of three decaying exponentials expressed as where, γ i ∈ R + for the MM model and BGK model are given in table 1. Interestingly, the normal stress σ (4.7b) has three contributions: (i) Newtonian stress (with Kn c 1 ), (ii) thermal stress (with Kn c 2 ) and (iii) Knudsen layer contributed stress (with exponentials). The origin of these different contributions to the stress are best explained by the stress balance equation (3.7a) from which the Newtonian stress and thermal stress arise due to the velocity gradient and the heat-flux gradient terms, respectively, and derivative of the higher-order moment m produces the Knudsen layer.  The only non-zero contribution to the pressure (4.7a) is due to Knudsen layer, given by the radial momentum conservation (3.4b).
As the NSF equations are only accurate up to first order in Kn (by Chapman-Enskog expansion), they capture the Newtonian stress (first order) but not the thermal stress (second order). The effect of these different contributions on the vapour flow will be studied later in § 9.
The temperature from the R26 equations reads θ = 2Pr 5Kn (4.8) The temperature θ has a Fourier contribution, (2Pr/5Kn)(c 2 /r), which comes from the right-hand side terms in the heat-flux balance equation (3.7b), plus an additional contribution due to the Knudsen layer (superposition of exponential functions), which describes the influence of higher moments in (3.7b). The coefficient k 1 , k 2 and k 3 in (4.8) depend on the collision model, they are given in equation (A 6) of the appendix. The solutions of the R26 equations for p, σ and θ given in (4.7) and (4.8) and other variables (given in appendix A) contain 5 integration constants, c 1,2 and a 1,2,3 . These constants are to be evaluated from interface conditions at r = 1. The boundary conditions for NSF, R13 and R26 equations are derived in the next section.

Formulation of the boundary conditions
The boundary conditions for evaporating and condensing interfaces for the R13 equations have been derived by Struchtrup & Frezzotti (2016), and here we follow the same steps for the derivation of the evaporating/condensing interface conditions for the R26 equations.
A phase interface, where particles approaching the interface from the vapour phase (c I k n k < 0) are described by the Grad 45-moment (G45) distribution function f | G45 and molecules entering the gas (c I k n k 0) are distributed according to a known distribution f + , see figure 2. Note that f | G45 reconstructs the distribution function in Hermite polynomials such that it reproduces all 45 moments appearing in the R26 equations. Similarly, for the NSF and R13 equations, G13 and G26 distribution functions are used, respectively.
Here, n k denotes the normal vector to the wall pointing into the vapour phase, and c I k is the velocity of vapour molecules in the reference frame of the interface. The explicit form of f | G45 can be found in the literature, for example by Gu & Emerson (2009). The distribution of molecules entering the gas is written as (Struchtrup & Frezzotti 2016)   where, ϑ is defined as the evaporation/condensation coefficient, and χ is the accommodation coefficient, defined as the fraction of reflected molecules being thermalised. Here c * k represents the specularly reflected values of c I k with respect to n k . In (5.1), θ l is the temperature of the liquid phase at the interface and p l is saturation pressure at θ l , following the definitions and the notations in (3.1). Note that, for non-evaporating interfaces (ϑ = 0), the above model reduces to the Maxwell accommodation model Gu & Emerson 2009).
At the interface, the distribution functionf of the vapour molecules is therefore written asf = f | G45 , c I k n k < 0 f + , c I k n k 0. (5.2) Once the distribution function at the interface is defined, the process of deriving the interface boundary conditions is the same as that for the wall boundary conditions Gu & Emerson 2009). Again, assuming small deviations, i.e. |p l | 1 and |θ l | 1, linearisation of the distribution (5.2) can be performed. The distribution function (5.2) is further simplified when relations (2.8) and (2.9) specifying the flow configuration are introduced.

Boundary conditions for NSF, R13 and R26
The evaporation and condensation boundary conditions for the R26 equations read Evaporation-driven vapour microflows and T = θ − θ l denotes the temperature jump at the interface. As expected, for nonevaporating interfaces (ϑ = 0), the above interface boundary conditions are reduced to the wall boundary conditions for the R26 equations derived by Gu & Emerson (2009). The R13 equations are obtained from the first three equations (5.3a-c), with R, m and ∆ replaced by the R13 constitutive laws (3.8) and Φ = 0.
The NSF boundary conditions, permitting a temperature jump, require that both the stress and the heat flux in (5.3a) and (5.3b) are replaced by the NSF laws (3.1) alongside R = m = ∆ = Φ = 0.

Comparison with classical theories
At thermodynamic equilibrium, the chemical potential and temperature are assumed continuous across the liquid/vapour interface. These conditions lead to the Clausius-Clapeyron relations, i.e. a relation between the saturation pressure and temperature at equilibrium. However, when the interface is not under local equilibrium conditions, it is usually assumed that the mass flux j (in dimensional form) obeys the Hertz-Knudsen-Schrage (HKS) relation, derived from the kinetic theory of gases as (Schrage 1953 where σ e and σ c denote the evaporation and condensation coefficients, respectively, representing the fraction of molecules that strike the interface and change phases from their initial liquid or vapour states, respectively. When written in the dimensionless and linearised form, the HKS relation (5.5) for σ c = σ e = ϑ reduces to which is the linear Hertz-Knudsen-Schrage (LHKS) relation. In the standard approach, the HKS is combined with the equality of liquid and vapour temperatures at the interface, i.e. θ = θ l . (5.7) Comparing with our boundary conditions for the R26 equations, we see that (5.3a) is a generalisation of the LHKS relation which contains the higher-order moments (R, Φ) and an additional term (σ /2) in (5.3a) due to the geometric curvature. In what follows, the results obtained from macroscopic models are compared with the linearised Boltzmann solutions in order to determine the gains in accuracy from using the R13 or R26 theories over the classical one.

Summary of models
Below, we summarise the hierarchy of models which are to be investigated in the forthcoming sections. Each is composed of bulk equations for r > 1 supplemented with boundary conditions at the droplet surface (r = 1).
A Mathematica script can be found in the online supplementary material (available at https://doi.org/10.1017/jfm.2018.85), which computes and lists these integration constants for the NSF, R13 as well as the R26 equations.
6. Evaporation of a spherical drop as Kn → 0 As a starting point we consider Kn → 0 for the two cases of pressure-driven flow (p l = 1, θ l = 0) and temperature-driven flow (p l = 0, θ l = 1), which allows us to compare the results with previously obtained solutions of the linearised Boltzmann equation (LBE).
To establish the limits of applicability of our analytic solutions, we will compare them with numerical solution of the LBE based on both the S-model (Chernyak & Margilevskiy 1989), and hard-sphere molecules (Sone & Onishi 1978;Takata et al. 1998;Sone 2007), for the case of complete evaporation (ϑ = 1). Therefore, henceforth, ϑ = 1.
6.1. Pressure-driven flow (p l = 1, θ l = 0) The LBE with S-model predicts that the mass-flux coefficient c 1 (= vr 2 ) → 0.6654 as Kn → 0 (Chernyak & Margilevskiy 1989) whilst Sone & Onishi (1978) applied the asymptotic method to the LBE for hard-sphere molecules to obtain c 1 → 0.6633 as Kn → 0. The value 0.66-0.67 can be directly compared to those obtained from the macroscopic models, as tabulated in table 2, for a gas composed of Maxwell molecules and for the BGK model.
For both the R13 and R26 theories, c 1 is predicted to within two significant figures of the LBE result with very little dependence on the collision model used. In contrast, the NSF (+jump) model only agrees to one significant figure and the NSF fails to even match this. Interestingly, even in the limit Kn → 0 where, naively, one may expect all of the models to coincide, including higher moments enables more accurate prediction of pressure-driven evaporative flow. This poor behaviour of macroscopic boundary conditions for the NSF equations stems from the approximation f gas f | G13 for molecules impinging the liquid from within the Knudsen layer.  TABLE 3. Results of pressure-driven calculations (p l = 1, θ l = 0). Values of lim Kn→0 (c 2 Pr/Kn) for different macroscopic models with complete evaporation (ϑ = 1) for BGK and MM models. These should be compared to −0.5250 for the S-model (Chernyak & Margilevskiy 1989).
Notably, many studies use fitting parameters, such as effective condensation coefficients or effective temperature jump coefficients, which are fitted to the asymptotic solutions to the kinetic equation for a specific problem. However, the boundary conditions used in this study contain no fitting parameters; hence the results presented are fully predictive and the model is not problem specific.
As one can see from table 3, the NSF fails to predict a gradient in c 2 as Kn → 0 whilst all the NSF(+jump)/R13/R26 theories provide reasonable agreement with the LBE results. Notably, the R26 theory gives the best prediction out of these models.
6.2. Temperature-driven flow (p l = 0, θ l = 1) Our benchmark results from the LBE equation (Chernyak & Margilevskiy 1989) as The value for lim Kn→0 (c 1 Pr/Kn) predicted by the moment equations for the temperature-driven case is the same as lim Kn→0 (c 2 Pr/Kn) for the pressure-driven case. These values were discussed in § 6.1 and are tabulated in table 3. The equality of c 1 for temperature-driven case and c 2 for pressure-driven case follows from Onsager's reciprocal relations (De Groot & Mazur 1984;Takata et al. 1998).
All macroscopic models describe the correct asymptotic limit for lim Kn→0 (c 2 Pr/Kn) = (5/2), which follows from Fourier's law (3.6). : curves of the normalised pressure p (a) and the temperature defect Θ (b) are plotted against scaled distance from the interface η, for the NSF, R13 and R26 theories with complete evaporation (ϑ = 1). Note that NSF curves do not deviate from zero. Also plotted are the results for the LBE from Sone & Onishi (1978).
In summary, our results indicate that the R26 equations with evaporation boundary conditions yield an excellent quantitative description in all cases, for both MM and BGK collision models, which are not matched by the NSF or R13 theories.

Evaporation from a planar surface: Knudsen layer structure and Ytrehus' problem
Whilst no flow profiles within the Knudsen layer are available from kinetic theory for the case of an evaporating sphere, Sone and coauthors have examined the structure of Knudsen layers for evaporation/condensation from a planar interface (Sone & Onishi 1978;Takata et al. 1998). The planar interface can be derived as a particular case of the solutions in the spherical geometry obtained in § 4, recovered by defining and taking an asymptotic limit Kn → 0. Here, η is the dimensionless distance from the interface, which is made dimensionless with the mean free path, λ = µ ∞ √ 2RT ∞ /p ∞ . It is convenient to introduce a dimensionless temperature defect Θ as which has been defined so as to vanish for the NSF and NSF(+jump) models so that Θ is a purely non-Fourier contribution to the temperature that only appears inside Knudsen layer.
7.1. Knudsen layer: pressure-driven flow (p l = 1, θ l = 0) Figure 3 compares the pressure p and temperature defect Θ for the different macroscopic theories with the results by Sone & Onishi (1978). The results indicate that NSF and NSF(+jump) cannot describe non-zero pressure and temperature defects, whilst the Knudsen layer profile for pressure (figure 3a) is captured quite accurately by both the R13 and R26 equations. However, the temperature defect Θ near the wall (figure 3b) is only accurately predicted by the R26 equations, with the LBE result around three times smaller than that predicted by the R13 equations at η = 0. are plotted against scaled distance from the interface η, for the NSF, R13 and R26 theories with complete evaporation (ϑ = 1). Also plotted are the results for the LBE from Sone & Onishi (1978).

7.2.
Knudsen layer: temperature-driven flow (p l = 0, θ l = 1) Unlike the pressure-driven case, in temperature-driven flow ( figure 4a,b), the pressure and temperature defect both vanish as Kn → 0 with O(Kn). Therefore, to obtain the next-order terms in Kn, the results are divided with the Knudsen number before the limit is taken to recover the most relevant values in this limit. Figure 4(a) shows that the LBE benchmark for the pressure profile exhibits a maxima at a finite distance from the boundary, inside the Knudsen layer, which is only captured by the R26 model. According to our analytical solutions, the ultimate source of the fascinating behaviour can be traced to the combination of unique exponentials appearing in the solution of the R26 equations in (4.7a). The actual Knudsen layer is a linear superposition of many exponential functions e −rγ i /Kn (Struchtrup 2008) -with different coefficients γ i -the R13 equations give only one such contribution so that non-monotonic behaviour cannot be captured. In contrast, the R26 system provides three exponential functions to form the Knudsen layer, thus capturing the more complex behaviour. For the temperature defect (figure 4b) the R26 solution still provides far greater accuracy in the Knudsen layer.

Summary of the behaviour of the Knudsen layer
The Knudsen layer is the kinetic boundary layer, found in the region of gas flow adjacent to the boundary. Whilst the R13 and R26 theories pick up the existence of Knudsen layer, it is the R26 theory which has access to three exponential functions to approximate this region. Therefore, R26 can (i) capture intricate qualitative features of the LBE solution and (ii) provide a good quantitative accuracy. Consequently, from here on, we focus on comparing the R26 equations with the classical approach of the NSF and its extension NSF(+jump). Henceforth, we consider only MM collision model. 7.4. Comparison of predicted mass-and heat-flux coefficients for evaporation from a planar surface (Ytrehus' problem) As a special case, we can find relevant flow coefficients for the steady evaporation from a planar surface into a half-space by using the reduced distance (7.1) and taking 980 A. S. Rana, D. A. Lockerby and J. E. Sprittles NSF NSF(+jump) R26 Ytrehus & Ostmo (1996) Ytrehus & Ostmo (1996).
an asymptotic limit Kn → 0. The half-space evaporation problem, often referred to as Ytrehus' problem (Ytrehus & Ostmo 1996), concerns an evaporating liquid surface kept at a temperature θ l and evaporation pressure p l . Far from the interface, the flow is in a uniform equilibrium state with bulk velocity v = v ∞ > 0 and q = 0. With the v ∞ and q prescribed, the coefficients c 1 and c 2 in (4.1) are (7.3a,b) and the temperature θ l and evaporation pressure p l are determined from the boundary conditions. Defining where the coefficients α p and α θ are determined by solving different macroscopic models. In table 4, we compare our results obtained from the NSF, NSF(+jump), R26 models with Ytrehus' solution using half-range moment method (Ytrehus & Ostmo 1996). As expected, NSF without the temperature jump boundary condition gives α θ = 0, but also gives 16 % deviation for the pressure coefficient α p . The NSF(+jump) with temperature jump boundary condition provide excellent agreement for the temperature coefficient α θ but underpredicts the pressure coefficient by 6 %. The agreement of the R26 equation for both the pressure and temperature coefficients is excellent ( 1 %), highlighting the importance of kinetic effects in the Knudsen layer next to the evaporating interface and reinforcing our conclusions from § 7.3. Figures 5 and 6 show the variations in flow coefficients (c 1 and c 2 ) with Knudsen number, for the pressure-driven and temperature-driven flows, respectively. In all the figures, the solid (red) curves show the R26 solutions, the NSF(+jump) solutions are denoted by the dashed (blue) curves and the NSF with the classical boundary conditions without temperature jump are depicted as dotted (green) curves. The results are compared with the LBE solutions from Chernyak & Margilevskiy (1989) and Takata et al. (1998), which are denoted by the symbols.
because the enthalpy of phase change must be supplied to the droplet to keep its temperature constant. The LHKS relation (5.6) used for NSF gives a constant mass flow coefficient for all Knudsen numbers, as shown in figure 5(a). The LHKS relation is derived for a planar interface; therefore, it does not take the curvature of the interface into account. Furthermore, the heat-flux coefficient is zero (figure 5b) for the NSF model, which follows from the no temperature jump condition at the interface. Consequently, the NSF model is insensitive to variations in Kn and loses validity for microflows. On the other hand, the NSF(+jump) model includes the temperature jump and the curvature effects, and gives a non-constant mass flux and non-zero heat flux, with an offset of approximately 10 % for Kn 0.4 in mass-flux coefficient and Kn 0.1 in heat-flux coefficient.
As illustrated in figure 5, the R26 results for c 1 are within 10 % of kinetic data for Kn 3 (figure 5a), while for c 2 they are valid for Kn 0.5 (figure 5b).
8.2. Temperature-driven flow (p l = 0, θ l = 1) In this case, figure 6 shows that c 1 < 0 and c 2 > 0; that is, a mass flux from vapour to liquid (condensation) and a positive heat flux from liquid are induced. This flow behaviour is governed by the boundary conditions (5.3a) and (5.3b). For θ l > 0, the heat goes from liquid to vapour giving c 2 > 0 from the boundary condition (5.3b); therefore the mass-flux equation (5.3b) leads to condensation at the interface, i.e. c 1 < 0.
From figure 6(a,b), it can be seen that all models agree well with the LBE simulations for small Kn, but in both cases the R26 model shows considerable improvement over the NSF-based models at moderate Kn. The R26 results are within 10 % of kinetic data for Kn 0.8 for both c 1 and c 2 , while the NSF(+jump) model shows similar agreement for Kn 0.3.

Onsager reciprocity relations
Due to the microscopic reversibility of the evaporation and condensation processes, the Onsager reciprocity relations should hold (Xu et al. 2006), which in this case state heat-flux coefficient c 2 = qr 2 as a function of the Knudsen number for the temperature-driven case (p l = 0, θ l = 1). The symbols, circles and diamonds, denote the results obtained in Chernyak & Margilevskiy (1989) and Takata et al. (1998), respectively, from the linearised Boltzmann equation.
that the heat-flux coefficient (c 2 ) driven by the pressure difference should be equal to the mass flow coefficient (c 1 ) driven by the temperature difference. In figure 7, our theoretical results for the heat-flux coefficient in the pressure-driven case and the mass flow coefficient in the temperature-driven case are compared to kinetic data from Takata et al. (1998). As shown, in the kinetic simulations Onsager reciprocity relations are valid for the entire range of Knudsen numbers. On the other hand, macroscopic theories derived here give agreement extending up to Kn 0.5 for R26 and Kn 0.03 for NSF(+jump), with approximately 10 % deviation. Note that for the classical NSF both c 1 and c 2 vanish. Therefore, within the range which we expect R26 to be accurate, Onsager reciprocity is approximately satisfied. It has been shown by Rana & Struchtrup (2016) that the macroscopic wall boundary conditions, derived using Grad's closure at the boundary, violate the reciprocity relation for high Knudsen numbers. Violation of the reciprocity relation is a serious concern for macroscopic boundary conditions; therefore, a more careful study of the thermodynamically admissible boundary conditions based on a full second law analysis with proper Onsager coefficients is required, which is beyond the scope of this paper. However, as observed in Rana & Struchtrup (2016), the evaporation/condensation boundary conditions derived here should allow us to identify and estimate the values of the unknown coefficients appearing in such thermodynamically admissible boundary conditions.
To summarise our findings in § 8, it has been shown in this section that the classical NSF model fails to even qualitatively describe LBE results. The R26 model shows significant quantitative improvements over the popular NSF(+jump) approach and, impressively, is within 10 % of all LBE data for Kn 1 as compared to Kn 0.1 for NSF(+jump). With respect to practical computations of microflows, this means the R26 theory can be relied upon for an order of magnitude smaller characteristic scales than the NSF(+jump) theory.

Non-Newtonian total normal stress
Having an analytic solution allows us to analyse the magnitude of competing contributions to this non-equilibrium flow. Here we study the total normal stress Downloaded from https://www.cambridge.org/core. (p + σ ) experienced by the liquid droplet, whose sign will be seen to change depending on which theory is used. Our primary focus is to understand how this happens by studying the competing contributions from the Newtonian stress, thermal stress and stress due to the Knudsen layer, as defined in § 4.3. Figure 8 shows the total normal stress obtained from three macroscopic models (NSF, NSF(+jump) and R26) for the pressure-driven case (figure 8a) and the temperature-driven case (figure 8b). The results evaluated are for a relatively small Knudsen number Kn = 0.1.
In pressure-driven flow (figure 8a), all data, including the classical NSF theory, predict essentially the same qualitative behaviour: a maximum at the interface and decay to zero in the far field. Nonetheless, there is a substantial gap in case of the classical NSF theory, in particular, near the interface, whereas a somewhat better agreement is observed between the NSF(+jump) and R26 models.
Notably, there is a fundamental difference in case of temperature-driven flow in figure 8( the NSF(+jump) model and positive for the R26 model. The ultimate origin of this behaviour can be explained in a concise way by the non-Newtonian contributions to the total normal stress in the R26 theory. Figure 9 shows the different contributions to the total normal stress for the pressure-driven case (figure 9a) and the temperature-driven case (figure 9b) obtained from the R26 equations for Kn = 0.1. In both of these figures, due to the small Kn, the Knudsen layer (denoted by the dotted purple lines) is restricted close to the interface. Furthermore, in both the cases, the thermal stress (dashed green lines) is of second order in Knudsen number, as the term Kn c 2 ∼ O(Kn 2 ), see § 6.
For the pressure-driven case (figure 9a), the hydrodynamic stress (denoted by the dot-dashed grey lines) is of first order in Kn. Hence, for this case, the thermal stress contribution as is one order of magnitude smaller than the Newtonian stress; so that the total normal stress (solid red lines) away from the interface is dominated by the Newtonian stress.
More interestingly, for the temperature-driven case (figure 9b), the Newtonian stress is of second order, as Kn c 1 ∼ O(Kn 2 ), i.e. of the same order as the thermal stress, but, critically, of opposite sign. As a result the thermal stress forces the total stress to become positive whereas without its contribution, as in the NSF, this contribution would be negative.

Summary and conclusion
This article has developed the macroscopic modelling of phase-change processes and shown that the R26 equations provide a significant improvement over conventional models, accurately approximating benchmark Boltzmann equation results up to Kn ≈ 1. Notably, the analytic solutions provide unique insight, by allowing us to decompose our solution and determine the relative contributions of different physical effects, particularly in and around the Knudsen layers.
The R26 equations have been considered for planar wall boundary-value problems, for example in Gu & Emerson (2009), Gu et al. (2010), Gu & Emerson (2014). However, the analytical solution in a spherical geometry and more importantly for a liquid-vapour phase boundary was attempted for the first time here. Our analytic solution highlights why the R26 model is superior to the R13: it has access to three exponential functions to describe the Knudsen layer rather than just one. These findings are consistent with those of Gu & Emerson (2014), where superposition of three exponentials from the R26 theory produced an improved description of the thermal Knudsen layer.
The general interest of the present study is to explore boundary-value problems for moderately rarefied gas flows, with an emphasis on evaporation from nanostructures. For example, evaporative cooling from nano pores ∼O(100 nm) can lead to thermal management in electrical systems directly at the source of hot spots. Such a device operates around atmospheric pressure, with mean free path λ ∼ O(10 nm) (based upon the saturation temperature) in the vapour phase, giving a Knudsen number Kn O(1). For the purpose of design optimisation, which requires repetition of simulations with different parameters, Boltzmann solvers can become impractical due to the large computational time involved. Macroscopic models, despite limitations on the their accuracy, give a computationally fast and detailed access to quantities of interest. The R26 theory developed here, which we have benchmarked against LBE, gives us a framework to study evaporation/condensation for Kn O(1), i.e. down to nanoscales.
Very recently, fundamental solutions (Green's functions) were derived for the linearised R13 equations by Claydon et al. (2017). The numerical framework based upon these fundamental solutions, generalised to account for phase-change phenomena, should allow for three-dimensional multiphase microflows to be modelled at remarkably low computational cost. Extension of the method of fundamental solutions for the R26 equations and development of a simulation tool for moment equations to capture geometrically complex flows will be the subject of future work.
Another line of inquiry is to study thermodynamically admissible boundary conditions for higher-order moment equations. It has been shown by Rana & Struchtrup (2016) that the macroscopic boundary conditions, derived using Grad's closure at the boundary, violate the reciprocity relation for high Knudsen numbers, i.e. there arises the problem of boundary conditions for a finite system of equations that approximate the microscopic boundary conditions for the Boltzmann equation. The boundary conditions must be consistent with the moment equations and the resulting problem must be resolved. Looking forward, a more careful study of the boundary conditions based on a full second law analysis with proper Onsager coefficients should lead to improved agreement between macroscopic theories and exact solutions of the Boltzmann equation.

986
A. S. Rana, D. A. Lockerby and J. E. Sprittles of exponential functions given by (4.7a). In Step 2, by substituting this ansatz in the momentum balance equation (3.4b) and integrating, we readily obtain the expression for σ given in (4.7b).
Step 2 introduces a constant integration b 1 . In Step 3, we shall use the stress balance equation ( where, b 2 is a new constant of integration. By substituting m from the previous equation into (3.10a), in Step 4, we obtain Φ. Note that, Step 4 does not introduces any additional integration constants. In Step 5, we substitute σ , m and Φ, from previous steps into the balance equation for m (3.9a), to obtain R. This step brings one new integration constant b 3 , scaling with r 2 , as where, expressions inside the brackets (·) do not depend on r. For any physically meaningful solutions, R must vanish as r → ∞, therefore b 3 = 0. An expression for ψ is obtained from the equation (3.10b) in Step 6. In Step 7, ω is obtained from the balance equation for R (3.9b), which again gives an integration constant b 4 scaling with r, therefore b 4 = 0. At this point, in Step 8, we have to obtain a solution for ∆ which satisfies equation (3.9c) and equation (3.10c), at the same time. We get ∆ from (3.9c) and substitute ω, ∆ and R in (3.10c) to get b 1 = 4 5 (5c 1 + 2c 2 )Kn, (A 3a) b 2 = 2c 2 5Pr M Pr R + 5c 1 + 2c 2 5Pr M 36Kn 2 , (A 3b) and a polynomial in γ as 120374γ 6 − 242756γ 4 + 119400γ 2 − 1 5306.4 = 0, (A 4) for the MM model, and for the BGK model. The positive roots of these polynomials were listed in table 1. Finally, the expression for the temperature θ is obtained by substituting ∆, R, σ and q into (3.7b) and then solving for θ. The new constant c 3 following from this integration vanishes due to (4.4). The coefficients k 1 , k 2 and k 3 appearing in (4.8) depend on the transport coefficients in the collision model and are given by k 1 = 9(63Pr φ + 80) 140Pr ∆ Pr ψ Pr φ , (A 6a) k 2 = − 81Pr M Pr φ + Pr ψ (64Pr ∆ + 9Pr φ (4Pr ∆ + 7Pr R + 2) + 80Pr R ) 36Pr ∆ Pr ψ Pr φ , (A 6b)