Evaporation boundary conditions for the R 13 equations of rarefied gas dynamics

The regularized 13 moment (R13) equations are a macroscopic model for the description of rarefied gas flows in the transition regime. The equations have been shown to give meaningful results for Knudsen numbers up to about 0.5. Here, their range of applicability is extended by deriving and testing boundary conditions for evaporating and condensing interfaces. The macroscopic interface conditions are derived from the microscopic interface conditions of kinetic theory. Tests include evaporation into a half-space and evaporation/condensation of a vapor between two liquid surfaces of different temperatures. Comparison indicates that overall the R13 equations agree better with microscopic solutions than classical hydrodynamics.


I. INTRODUCTION
2][3][4][5] Gas rarefaction leads to the occurrence of phenomena such as velocity slip and temperature jump at boundaries, Knudsen layers in front of boundaries, transpiration flow, thermal stresses, or heat transfer without temperature gradients, all of which are reproduced by solutions of the R13 equations.Proper modeling of boundary conditions (BCs) is essential to obtain a meaningful description of rarefied flows, and below we develop and test conditions for liquid-vapor boundaries with condensation and evaporation.With this, the range of application of the R13 equations is extended in particular towards micro-devices with phase change; possible applications include, e.g., gas supply channels in micro-fuel cells.
Indeed, it is well known that rarefied gas flows cannot be accurately described by the Navier-Stokes-Fourier (NSF) equations of classical hydrodynamics. 3,6,7While adding jump and slip boundary conditions to NSF introduces some of these effects into classical hydrodynamics, all bulk rarefaction effects remain elusive. 8A fully accurate description of rarefied flow behavior requires solution of the Boltzmann equation, e.g., by the direct simulation Monte Carlo (DSMC) method 9 or by direct numerical simulation, 10 but typical computational times are often significantly above those for hydrodynamics, in particular, when the degree of rarefaction is mild.Macroscopic models for rarefied flows, such as the R13 equations, extend the equations of hydrodynamics by accounting for a small number of additional flow variables, yet describe bulk a) struchtr@uvic.ca.URL: http://www.engr.uvic.ca/˜struchtr/rarefaction phenomena in good approximation, and thus offer an alternative tool for computation. 4,5he R13 equations are derived as approximations of the Boltzmann by means of the order of magnitude method, 11,12 which combines elements of the Chapman-Enskog 6 and Grad 13,14 methods.The resulting equations avoid the problems exhibited by the individual methods.Higher order Chapman-Enskog expansions yield the Burnett 15 and super-Burnett 16 equations, which are well known to be unstable in time dependent processes 17,18 and to give unphysical solutions in steady state. 19Grad's 13 moment equations, or systems with more moments, are stable but exhibit unphysical sub-shocks for high speed flows. 20Moreover, the Grad equations are not linked to the Knudsen number; hence, it is difficult to know a priori which set of moments should be considered for a given process.In contrast, the R13 equations are stable, have no sub-shocks, and are derived to be of third order in the Knudsen number (super-Burnett order). 12The method can be extended to higher moment numbers. 21,22etails of the R13 equations depend on the molecular model, and in the following we consider only the wellestablished equations for Maxwell molecules. 12Equations for hard sphere molecules, 23 polyatomic gases, 24 and gas mixtures 25 are discussed elsewhere in the literature.
Just as the transport equations are derived from the Boltzmann equation, the corresponding macroscopic boundary and interface conditions are derived from the microscopic boundary and interface conditions for the Boltzmann equation. 26or this we use an extended Maxwell boundary model 27 with a condensation/evaporation coefficient (probability that a vapor particle condenses at impact on the liquid) and an accommodation coefficient (probability for diffuse vs. specular reflection for a non-condensing particle). 28The derivation follows the same line as that of the wall boundary conditions for non-condensing interfaces.Naturally, for a vanishing condensation coefficient, the boundary conditions for solid walls are recovered.
][33][34][35] For both cases, the R13 and NSF equations can be analytically solved, and their solutions are compared with DSMC simulations.The ability of the macroscopic models to model rarefied gas flows with evaporating or condensing interfaces is discussed.
While we discuss only one-dimensional problems as applications, the boundary conditions are developed for full three-dimensional geometry, and multi-dimensional problems will be considered in the future.
We note that in the present paper, we consider problems that focus on transport in the vapor, with given liquid data at the interface.The full simulation of problems involving liquid-vapor interface requires a matching model for the transport in the liquid, which for completeness is outlined in the Appendix A.
The remainder of the paper proceeds as follows: Sec.II will briefly review the R13 equations and their relation to kinetic theory, as well as to classical hydrodynamics.Based on this, the evaporation and condensation boundary conditions for the R13 equations will be derived in Sec.III from the corresponding boundary conditions for the Boltzmann equation.The relation of the extended R13 evaporation/condensation conditions to those for hydrodynamics, such as Hertz-Knudsen-Schrage, will be discussed as well.In Sec.IV, we solve simple one-dimensional evaporation problems-a half-space evaporation problem and heat and mass transfer between two liquid layers-with the R13 equations and compare results to kinetic theory simulations and results from classical hydrodynamics.The paper ends with our conclusions.Some preliminary results were published before; 36 the present contribution gives all required details for the derivation of the boundary conditions, corrects several errors, and expands on the results presented and discussed.

A. Moments and moment equations
We briefly recall the origin of the R13 equations from kinetic theory of monatomic gases.The aim of kinetic theory is to find the velocity distribution function f (x i , t, c i ) which is defined such that f (x i , t, c i ) dcdx gives the number of particles in the phase space cell dcdx at time t; x i is the location in space and c i is the microscopic velocity of a particle.In a microscopic approach, the distribution function is the solution of the Boltzmann equation, 6,7 while in a macroscopic approach, one derives transport equations as a set of suitable moments of the Boltzmann equation, where the resulting moment equations are closed by means of an ansatz for the distribution.This so-called moment method goes back to the work of Grad, and his 13 moment equations are the best known set of macroscopic equations. 13,14The basic 13 variables are mass density ρ, macroscopic velocity i , temperature T, anisotropic stress tensor σ ij (with σ ii = 0), and heat flux vector q i , which are moments of the distribution function as ρ = m fdc, ( Here, u = 3 2 RT = 3 2 θ is the specific internal energy, θ = RT is temperature in specific energy units with the gas constant R, and C i = c i i is the peculiar velocity.Indices in angular brackets denote the symmetric and trace-free part of a tensor. 3he corresponding moment equations are the conservation laws for mass, momentum, and energy, which can be written as and the balance equations for stress and heat flux The collision terms of the above equations were determined for Maxwell molecules, µ is the shear viscosity and κ = 15 4 µ is the heat conductivity. 3or small Knudsen numbers, the Chapman-Enskog method 3,6 can be used to reduce the equations for stress and heat flux to the Navier-Stokes and Fourier laws,

B. Constitutive equations
In addition to the 13 variables introduced above, the 13 moment equations ( 2)-( 4) contain the higher moments, Phys.Fluids 29, 092004 (2017) In order to close the system of moment equations, constitutive equations for these must be provided.The classical Grad closure 3,13 simply leads to The regularized 13 moment equations arise from an alternative closure, 1 which accounts for parts of, but not the complete, transport equations for m ijk , ∆, R ij .The equations are best derived by means of the order of magnitude method, 11,12 which mixes elements of the Grad and Chapman-Enskog methods to find the relevant terms of the infinite moment system that are required to have accuracy at a certain order in the Knudsen number, which is defined as Kn = µ √ θ pL , where L is the relevant macroscopic length. 3The regularized 13 moment equations arise as the appropriate set of equations at 3rd order in the Knudsen number (super-Burnett order) and consist of Eqs. ( 2)-( 4) and the constitutive equations (pressure obeys the ideal gas law, p = ρθ), Some comments regarding the closure are in order: First, the closure depends on the molecular model employed in the Boltzmann equation and is performed easiest for Maxwell molecules, which results in the coefficients as given above.For other molecular models, the transport equations, and the coefficients in the closure, change.This was discussed in Ref. 23, where the linear transport equations for hard sphere molecules were derived.Non-linear second order (i.e., Burnett order) equations for hard sphere molecules were derived before; 3,11 these could be combined with the linear regularization at third order to give a non-linear system of equations for non-Maxwellian molecules; this will be discussed elsewhere.
Second, the order of magnitude method allows for some ambiguity of the equations, since it aims only for the proper terms at a given Knudsen order.Since the R13 equations are of third order in Kn, one can modify them in a way that the third order contributions remain unchanged, while higher order contributions change.Over the years, this ambiguity has led to several versions of the R13 equations presented in the literature.The first derivation 1 of the R13 equations relied on the artificial assumption that higher moments would change on a faster time scale and a Chapman-Enskog expansion for the fast scale; the resulting equations include some terms that contribute at the fourth order in Kn, which, one can argue, should not be present in a third order theory.The re-derivation by means of the order of magnitude method removed these fourth order terms. 12Even then, the number of boundary conditions required for the solution of the equations differed between the linear and the non-linear equations. 26For the closure presented (8), the equations were modified such that linear and non-linear equations require the same number of boundary conditions, 37 and this is the preferred version of the R13 equations for slow flows.For non-linear processes such as shock waves, various versions show marked differences in the results. 38ince the full Boltzmann collision term is difficult to handle, many authors base numerical simulation of rarefied gas flows on the Bhatnagar-Gross-Krook (BGK) model, 39 which approximates the Boltzmann collision term with a rate ansatz.Qualitatively, the BGK model gives an accurate description of rarefied gas flows, but due to its simplicity, it does not lead to quantitative agreement; most prominently, it gives a wrong value for the Prandtl number. 3,7In order to be able to compare our description of evaporation/condensation processes in rarefied flows with BGK solutions, we occasionally consider the R13 equations for the BGK model, where the heat conductivity is κ BGK = 5 2 µ, and the closure relations are similar to the above but have different coefficients in almost all terms, 12

C. Distribution function in the bulk
At a physical boundary, such as a wall or a phase boundary, the distribution function of incoming particles changes due to reflection, condensation, and evaporation. 26,28To derive boundary conditions for moments from the boundary conditions for the phase density, a detailed expression for the distribution function in terms of the moments is required.For this we use the same distribution that is used for the closure of the moment system, namely, the Grad distribution for the 26 moments (1, 6).The Grad-26 moment distribution is the local equilibrium distribution f M times a suitable polynomial, 3 where the equilibrium function is the local Maxwellian Insertion of this distribution into the definitions (1, 6) of moments gives identities.

III. MICROSCOPIC AND MACROSCOPIC BOUNDARY CONDITIONS A. Distribution function at evaporating interface
At an evaporating liquid interface, some vapor particles that hit the interface condense, while those that do not condense are reflected; moreover, some particles are injected into the vapor by evaporation from the liquid.We write the distribution function directly in front of the interface as where f is the distribution of incident particles (negative velocity C I n normal, and relative, to the interface) and f + is the distribution of emitted particles (positive velocity C I n normal to the interface).The prime at the velocity of incident particles simplifies to distinguish between incident and emitted particles.For finding boundary conditions for moments, we shall describe the vapor by the R13 distribution function (10), i.e., we set f = f |G26 .
The distribution of particles leaving the interface is due to evaporating particles and the reflection of non-condensing particles back into the vapor.For the latter we follow the classical Maxwell model, which assumes that particles are either specularly reflected or thermalized and leave in a Maxwellian. 27,28he distribution of evaporating particles is obtained from the assumption that the liquid side of the interface is in local equilibrium, 33,40,41 and thus evaporating particles leave in a Maxwellian distribution f M p sat (θ L ) , θ L , C I , where p sat (θ L ) is the saturation pressure corresponding to the temperature θ L of the liquid at the interface (see Appendix A).For thermal equilibrium, the outgoing distribution reduces to The emitted distribution function consists of three terms describing evaporating particles, specularly reflected particles, and thermalized reflecting particles, respectively, Here, ϑ is the evaporation/condensation probability and χ is the accommodation coefficient, defined as the probability that a reflected particle is thermalized.Note that in thermal equilibrium, where all distributions are equal to f M p sat (θ L ) , θ L , C I , the distribution f + must also be equal to the Maxwellian.Hence, the factors in Eq. ( 13) must add to unity.The velocity C I i is the velocity of a vapor particle as seen from an observer resting with the liquid at the interface; C I n = C I j n j is the velocity normal to the interface; C I i − 2C I n n i is the specular reflection velocity (flipping of normal part) and the notation is such that f |G26 C I i − 2C I n n i denotes the distribution of specularly reflected particles.The pressure p in the Maxwellian for thermalized particles is determined from the condition that non-condensing particles must return to the vapor, which gives Since the Maxwellian is linear in pressure, and both Maxwellians in Eq. ( 13) depend on the liquid temperature, we can introduce the abbreviation and rewrite the outgoing distribution in a more compact notation, Continuity conditions for fluxes of moments are used to find the boundary conditions for the moments. 4,5,14,26For a general description, we write the moments as where Ψ A (c k ) denotes the moment generating tensor polynomials in c k with a suitable multi-index A, and the flux of the moment u A then is Only some of the fluxes must be continuous at the interface, as will be discussed further below.

B. Frames of reference
For the evaluation of interface conditions, we found it best to consider the fluxes from the viewpoint of an observer who slips with the gas along the liquid-vapor interface.For proper bookkeeping, we have to deal with three different particle velocities, which differ by the chosen frame of reference: (a) the velocity relative to the liquid-vapor interface, which moves with velocity v I i (index n refers to the interface normal and indices t α with α = 1, 2 indicate the two tangential directions), (b) the usual peculiar velocity of kinetic theory, 3 which is the velocity relative to the motion of the gas, which moves with velocity v i , and (c) the velocity for the slipping observer Evaporation velocity V n and slip velocity V t α are defined as the relative macroscopic velocities between interface and gas, With the interface normal n i pointing into the gas, positive V n = V i n i implies evaporation and negative V n implies condensation.Here, we chose the tangential velocity of the interface as the velocity of the liquid.Note, however, that in evaporation the normal velocity of the liquid differs from the normal velocity of the interface (see Appendix A).
The various particle velocities are related as

C. Macroscopic boundary conditions
For an observer slipping with the gas along the liquidvapor interface, who observes the particle velocity Ĉi , the normal flux F Ak n k (18) of a moment computed with the distribution function directly at the interface, i.e., the distribution f int of Eq. ( 12), must be equal to the normal flux computed with the distribution function f |R13 of the gas just in front of the interface, 3 that is, or, since in Eq. ( 12) we substitute Here, following the work of Grad, 14 we have to use continuity only of those fluxes that are odd in Ĉn , which implies functions ΨA which are even in Ĉn .The 13 variables of the R13 equations are moments based on the weights, and of these the following tensor components are even in the normal velocity c n : For the observer slipping at the liquid interface, this translates into which must be inserted into Eq.( 25).We also insert f + from Eq. ( 13) and substitute ( Ĉn → − Ĉn ) in the second integral to find after some straightforward manipulations the interface conditions read (recall that ΨA is even in Ĉn ), where P * is given by Eqs. ( 14) and (15).Since all distributions are Gaussians time polynomials, the actual evaluation of the boundary conditions is straightforward, but somewhat cumbersome.Some care has to be taken to use the proper velocities for evaluation.For instance, the Maxwellian of evaporating and thermalizing particles is centered in the restframe of the liquid-vapor interface.For the evaluation, we have to consider that Ĉn = C I n and Ĉt α = To proceed, we collect results for the relevant integrals.For the Maxwellian half-space integrals, we find For the other integrals in Eq. ( 29), we have to re-write the distribution function f |R13 (10) for the slipping observer, i.e., with the velocity Ĉi .To make the procedure manageable, we consider only linear effects in non-equilibrium.For this we insert into Eq.( 10) and expand to linear terms in all non-equilibrium quantities V n , ∆, q k , σ ij , R ij , m ijk .In particular, this restricts the results to evaporation velocities V n well below the speed of sound.With this restriction, we find From this, the full fluxes of f |G26 are calculated as and the corresponding half fluxes as Here, we have introduced the effective pressure Π as

D. Evaporation boundary conditions for R13
The integrals of Sec.III C are all that is required to find the boundary conditions for the R13 equations at evaporating and condensing interfaces by the evaluation of Eq. ( 29).Before we list the results, we point out the deficiency of the R13 equations to fully resolve Knudsen layers, which was discussed, e.g., in Ref. 42.Indeed, full resolution of the Knudsen layer requires a larger number of moments, which is not feasible in practice. 22,43However, solutions of the R13 equations between solid walls, i.e., no evaporation or condensation, have shown that their ability not only to describe jump and slip at the boundaries but also to approximate Knudsen layers leads to a marked improvement of the overall simulation.
The usual way to improve the boundary conditions is to introduce correction coefficients for jump and slip into the interface conditions, which is usually done in jump and slip hydrodynamics. 7,26A more careful approach to correct the boundary conditions for R13 equations uses the ideas of thermodynamics of irreversible processes 44,45 to derive interface conditions with tunable Onsager coefficients from the entropy generation at the interface, based on the R13 approximation of the entropy flux. 42A similar approach could be useful for the evaporation case, but this is deferred to future work.
The expression for the evaporation mass flux results for ΨA = 1 as Hence, the evaporation flux is determined through the difference between the saturation pressure p sat (θ L ) of the liquid at the interface (weighted with the factor θ θ L ) and the effective pressure Π.This expression is a generalization of the classical Hertz-Knudsen-Schrage law [46][47][48][49] to the R13 equations.For non-evaporating surfaces, where ϑ = 0, this reduces to the statement that there is no flow through the surface, V n = 0.The classical Hertz-Knudsen law replaces Π with the pressure p of the vapor and replaces the factor √ 2/π by 1/ √ 2π.The correction of the latter factor due to Schrage accounts for the fact that the vapor is in a non-equilibrium state, as it is in the R13 equations.
The other conditions are generalizations of the established wall boundary conditions for the R13 equations, 26 to which they reduce for non-evaporating interfaces (ϑ = 0).Due to the evaporation and condensation processes, there are additional contributions with the evaporation flux ρV n , and the evaporation coefficient appears in the overall slip and jump coefficients ϑ+χ(1−ϑ) 2−ϑ−χ(1−ϑ) (which for ϑ = 0 reduces to the familiar form χ 2−χ ).The conditions for m t α t β n and R t α n are extensions to two-tangential directions, which were not stated in the literature before so that the interface conditions below allow us to formulate three dimensional problems.
In summary, we find the following interface relations: Generalized velocity slip condition, generalized temperature jump condition, generalized interface conditions for higher moments, In Sec.IV, we shall proceed with the evaluation of the above interface conditions and study applications to some simple evaporation problems.

E. Interface conditions for hydrodynamics
In the hydrodynamic limit, only first order contributions in the Knudsen number are retained in the equations.Analysis of the bulk transport equations by means of the order of magnitude method 12 reveals that stress σ ij and heat flux q n are of first order in Kn and are given by the Navier-Stokes-Fourier equations (5), while m ijk , ∆, and R ij are of second order.The evaporation mass flux ρV n is not linked to the Knudsen number, but we have already considered it to be small of first order in all equations.Thus, for the hydrodynamic limit, we consider only terms linear in V n , V t α , σ ij , and q i in the above, to obtain the interface conditions for hydrodynamics which are as follows: (a) evaporation/condensation mass flow, (b) slip condition, and (c) evaporation/condensation heat transfer, The subscript NSF indicates that these are the appropriate boundary conditions for classical hydrodynamics.It is implied, of course, that in the above, stress tensor and heat flux are given by the laws of Navier-Stokes and Fourier (5).7][48][49] The slip condition (43) is the typical slip condition for a gas, 7 only that the overall slip coefficient now depends on the evaporation coefficient ϑ as well.Typically, the contribution with σ nn in Eq. ( 42) does not appear in the literature, since it vanishes for NSF in planar geometry, as, e.g., in the problems below.It should be included, however, since it will contribute for curved phase interfaces, such as in droplets and bubbles.
Quite often in the discussion of evaporation mass flows, a condition for energy flux (44) is not given.It is clear, however, that the simulation of an evaporation experiment requires the full set of transport equations with the appropriate boundary conditions; hence, an interface condition for energy flux (and possibly higher moments) is necessary.We point to our earlier work 49 where an evaporation experiment was evaluated with a variant of the Schrage model, and it became clear that the evaporation flow was mainly determined by the amount of heat transferred to the interface through the liquid, while the mass flow condition determined the liquid temperature, and the heat flux condition determined the temperature jump between liquid and vapor at the interface.
For one-dimensional processes in planar geometry, with transport only normal to the interface, and for only small deviations from equilibrium, Eqs. ( 42) and ( 44) reduce to The above can be rewritten as which is the natural form in the theory of linear irreversible thermodynamics 28,50 with a symmetric matrix of Onsager coefficients.
Further reduction by assuming full accommodation ( χ = 1) and inversion of the matrix gives with the symmetric matrix of resistivities 28 rαβ Exact calculations based on kinetic theory (with full accommodation, χ = 1) yield explicit corrections to account for Knudsen layer effects 33,50 with the resistivities We shall examine both sets of resistivities in the examples below.

A. One-dimensional equations
To put the R13 equations with evaporation/condensation to test, we now consider flows in simple one-dimensional geometry and steady state.For simplicity, we ignore all details of mass and heat transfer through the liquid and consider the temperature of the liquid at the interface as given.The normal of the interface points into the x 1 = x direction, and all flow properties are functions only of this coordinate.
For these first tests, we consider small deviations from a rest state where vapor and liquid are in equilibrium at temperature θ 0 , hence the reference pressure is p 0 = p sat (θ 0 ) = ρ 0 θ 0 .We use the rest state data and the length scale L to make the variables dimensionless.With all flows only in the x-direction, the variable space is reduced to We simplify the R13 equations ( 2)-( 4) and ( 8) for geometry, keeping only first order terms in the velocity.The linearized dimensionless conservation laws (2) become (for a simpler notation, we omit the circumflexes that indicate dimensionless quantities) which are easily integrated to give constant mass flow J 0 , constant normal stress P 0 , and constant overall energy flux Q 0 , ρv = J 0 = const., p + σ = P 0 = const., With this, the (linearized, dimensionless, one-dimensional) R13 constitutive equations for Maxwell molecules (8) or the BGK model ( 9) reduce to where the coefficients β, γ distinguish between Maxwell molecules ( β MM = γ MM = 1) and the BGK model ( β BGK = 3 2 , γ BGK = 7 6 ).Note that, due to linearization, products such as J 0 ∂θ ∂x and J 0 ∂σ ∂x can be ignored.The linearized balance equations for the xx-component of stress (3) and the x-component of heat flux (4) reduce to (with These can be integrated easily to give where K, A, B are constants of integration.We note that σ is of the Knudsen layer type, i.e., it decays exponentially away from the interface on the scale of the mean free path.We also note that classical hydrodynamics gives σ = 0, and θ = K − 4q 0 x 15Kn , which is the case for A = B = 0.
For the full solution, we have to find the 6 constants of integration {J 0 , P 0 , q 0 , K, A, B}, which requires 6 boundary conditions in total, that is, 3 boundary conditions on either side of the domain.
When a boundary is a liquid-vapor interface, we will use the appropriate boundary conditions from ( 36)-( 41) which in linearized and dimensionless form reduce to (no tangential components) Here it is assumed that all dimensionless pressures p sat (θ L ), P 0 , and all dimensionless temperatures θ L , θ, are close to unity.J 0,n , q 0,n are the products of the flows J 0 , q 0 with the normal at the respective boundary.
For NSF, we consider the interface conditions (48), either with the simple resistivities (49) or with the corrected resistivities ( 50), which we write for easy comparison with the above as Since we only consider the temperature of the liquid at the interface, but not the detailed transport processes through the liquid layers, the relations for these given in Appendix A are not required.

B. Half-space problem
We first consider the classical problem of the steady evaporation into a half-space, consisting of an evaporating interface at temperature θ L with evaporation pressure p sat (θ L ) where evaporation is forced by prescribing the pressure p ∞ = P 0 at a large distance of the interface.Dimensionless temperature θ ∞ and velocity v ∞ = J 0 are controlled such that the vapor at a large distance is in a drifting equilibrium state.
A comprehensive account of the kinetic theory treatment of the problem is given in Ref. 30.Here, we compare the R13 solution with Pao's results from the linearized BGK model 31 and with the linearized form of non-linear jump formulas derived by Ytrehus by a moment method solution of the Boltzmann equation for Maxwell molecules. 29Both formulations assume perfect evaporation (ϑ = 1).
The conditions for this problem require q 0 = σ ∞ = 0 and hence the solution of the problem follows from Eq. ( 56) as With the pressure p ∞ prescribed, the unknowns of the problem are temperature θ ∞ , evaporation speed v ∞ , and the constant A for the normal stress.These follow from the boundary conditions at the interface (57), which for this problem can be brought into the form For hydrodynamics, the corresponding equations based on (58) must be used, with the resistivities ( 49) or (50); after linearization, For comparison of the various models, we follow the work of Ytrehus and consider the dimensionless ratios which are just the resistivities r11 and r12 of Eq. ( 48) multiplied by 2 √ π.We compare results from R13, Ytrehus (Y), Pao (P), and hydrodynamics [NSF and corr.NSF, with the coefficients (49) or (50), respectively], for which we find and The corrected NSF values (corr.NSF) agree with the exact data, since they were fitted to similar results.Pao's results agree with our NSF results (49), that is, without any Knudsen layer correction.Phys.Fluids 29, 092004 (2017)   For this problem, the R13 equations for Maxwell molecules ( β = 1) and the BGK model ( β = 3/2) give almost identical results.The coefficient α p agrees well with Ytrehus' result (relative error 0.7%), which in turn agrees well with our own numerical solution for the BGK model.For the temperature coefficient α θ , however, the R13 equations differ from the accurate value (also confirmed by numerical solution of BGK) by about 10%.This deviation is similar to the deviation of temperature jump and slip for the R13 equation in half-space problems which was pointed out in Ref. 51 and discussed in detail in Ref. 42.
The results presented next will show that for larger Knudsen numbers, all macroscopic models cannot accurately model heat flux and temperature profile simultaneously.The coefficient α θ is extracted from the temperature profile, that is, its deviation is related to this conflict.
To close this section, we point out that in this particular half-space problem, the non-convective heat flux, q, vanishes, which implies that the boundary conditions cannot be fully explored.This is evident from the observation that the resistivity r22 of Eq. ( 48) does not occur in the solution for the NSF problem.

C. Heat and mass transfer between two reservoirs
Next we study heat and mass transfer between two liquid reservoirs (see Fig. 1), comparing predictions of macroscopic equations (NSF, R13) with Direct Simulation Monte Carlo (DSMC) results.This problem is richer than the previous one, in that it has a non-vanishing heat flux.Also, we will look at detailed profiles, including resolved Knudsen layers, of temperature and normal stress.
To be specific, we consider one-dimensional transport between two liquid layers identified with sub-or super-scripts 0, 1 located at (dimensionless) locations x 0,1 = ∓ 1 2 .For the solution, we have to consider boundary conditions on both sides of the domain.The interface normal points into the vapor, i.e., into the positive direction at x 0 = − 1 2 and into the negative direction at x 1 = 1 2 ; accordingly V n (x 1 ) = −V n (x 0 ) = −J 0 , etc.We need to consider the boundary conditions (57) at both boundaries.The three pairs (57) of boundary conditions for V n , q n , m nnn are best applied by taking their pairwise sums and differences, respectively.After some calculation, the solution of the linear problems ( 53) and ( 56) assumes the form FIG. Setup for heat and mass transfer between two reservoirs.
Two of the pairwise sums of the BC give θ 0 L − θ 0 + θ 1 L − θ 1 = 0, which was used to simplify the above.The evaporation fluxes J 0 , the heat flux q 0 , and the amplitude of the Knudsen layer A are obtained from a linear system which results from subtracting the interface conditions at both sides, Inversion of this system yields the integration constants {J 0 , q 0 , A} and insertion of these into (65) gives detailed profiles of temperature and stress.
The NSF solution is obtained by setting A = 0 in Eq. ( 65), which implies σ = 0 so that p = P 0 is the homogenous pressure in the vapor.Subtraction of the interface conditions (48) on both sides yields two equations for mass and heat flux that, for comparison with the corresponding R13 expressions from above, are best written as The final solution of the problem is obtained by solving (67) for J 0 and q 0 and inserting the results (with A = 0) into Eq.( 65).R13 and NSF results will now be compared to DSMC simulations (for β = 1), for which we set the evaporation coefficient ϑ to unity so that ϑ For NSF, Phys.Fluids 29, 092004 (2017) we will use the classical result (50), and also the simple result (49), which does not account for Knudsen layer corrections.Since we are solving only the linearized equations, the DSMC simulations are performed for small deviations from equilibrium.Specifically, we prescribe the dimensionless temperatures of the liquid layers as θ 0,1 = 1 ± ∆θ.The solutions depend strongly on the corresponding saturation pressures and we will show results for psat θ 0,1 L = 1±∆p with ∆θ = 0.05, ∆p = 0.05 and ∆θ = 0.01, ∆p = 0.075, respectively, where the second case produces the often discussed inverted temperature profile. 31All results will be shown for a variety of Knudsen numbers in the range Kn ∈ {0.039, 1}.
With these data for temperature and pressure, the solutions of the linearized macroscopic equations yield profiles of temperature deviation and stress that are anti-symmetric relative to the midpoint x = 0.Although the deviations from equilibrium are small, the DSMC results show non-symmetric profiles due to small non-linearities.In order to have a proper comparison, we use not the original DSMC data but symmetrized temperature and stress profiles, which are obtained by appropriate averaging of the left and right parts of the profiles.

Standard temperature profile
When the dimensionless saturation pressures of the liquid layers are psat θ 0,1 L = 1 ± 0.05, the temperature profile shows the expected behavior, with a jump at both boundaries FIG. 2. Kn = 0.078: Temperature and stress profiles for ∆θ = 0.05 and ∆p = 0.05, for DSMC simulation (symmetrized; green, dashed), corrected NSF (blue, dashed), uncorrected NSF (black, dotted-dashed), and R13 (red).and, since the left boundary is hotter, a decreasing temperature in the bulk.Figures 2 and 3 show the temperature curve and the corresponding normal stress for Kn = 0.078 and Kn = 0.235, respectively, as functions of location x.Evidently, all macroscopic models have some difficulty in matching the full details of the temperature curve more so at larger Kn.Compared to DSMC, the NSF models give larger jumps at the boundaries (recall that the dimensionless liquid temperatures are θ 0 = 1.05 and θ 1 = 0.95, respectively) and a somewhat flatter temperature profile; they cannot describe normal stress (σ = 0).The R13 equations give jumps and temperature slopes close those of DSMC and give the normal stress in good approximation.
Maybe more important than the finer details of the profiles is the description of overall mass and heat transfer, that is, the dependence of the integration constants J 0 and q 0 on the Knudsen number.Figure 4 shows J 0 , q 0 , and the boundary stress σ (x 1 ) for Knudsen numbers in {0, 1}, for DSMC (green dots), corrected NSF (blue dashed), uncorrected NSF (black dashed-dotted), and R13 (red continuous).
From the figure, we see that the corrected NSF equations (blue, dashed) give, for this problem, a good description of mass and heat transfer for all Knudsen numbers considered.As seen before, they do not give any non-equilibrium stress (σ = 0) and the temperature profiles do not match.That is, the accuracy in mass and heat flux is offset by the inaccurate temperature and stress profiles; it is not possible to adjust the resistivities rαβ to fit both fluxes and temperature profile.The uncorrected NSF equations (black, dashed-dotted), which do not correct Knudsen layer effects, do not provide a good match for fluxes or profiles.
The R13 equations (red, continuous) give a reasonable description of fluxes and stresses for relatively small Knudsen numbers (Kn ≤ 0.25), with relative errors for mass flux of less than 2%.
The visual fit for heat flux q 0 is good for all macroscopic systems, but it should be noted that for small Knudsen numbers, the energy flux q 0 is rather small, and there are fluctuations in the DSMC simulations, which leads to larger relative errors.According to (53), the total energy flux is 2 J 0 θ 0 + q 0 , that is, energy transport-in particular at smaller Knudsen numbers-is dominated by convective transport.For the total energy flux Q 0 , the corrected R13 equations differ from DSMC by not more than 1% for all Knudsen numbers considered, which is close to the accuracy of corrected NSF, where the relative error is below 0.9%.

Inverted temperature profile
Increase of the saturation pressure difference ∆p for unchanged temperature difference ∆θ implies an increase in the enthalpy of vaporization.For the transport between two liquid layers, this implies increased convective transport of energy or smaller conductive contribution to transport.For a certain range of values, and for the left liquid layer hotter than the right, it is possible that the temperature gradient in the bulk is inverted so that the vapor in front of the hotter liquid on the left is colder than the vapor in front of the colder liquid on the right. 31,34,35,52It is worth mentioning that although a clear experimental evidence of temperature inversion during evaporation/condensation has not yet been given, inverted temperature profiles have been found in molecular dynamics simulations, 34 thus showing that their occurrence is not an artifact of kinetic boundary conditions.Temperature inversion occurs for DSMC simulations at ∆p = 0.075 and ∆θ = 0.01, for which we show the same curves as before.Figure 5 shows temperature and stress for Kn = 0.078.Note that the left and right temperatures are θ 0 = 1.01 and θ 1 = 0.95, respectively, which implies that the temperature jumps at the liquid-vapor interfaces are smaller than in the previous case (∆p = 0.05) and that the resolution of temperature is finer.Indeed, the R13 solutions show temperature Knudsen layer contributions due to the last term in Eq. (65) 2 , which are present in Fig. 2 as well but, due to a different scale, hardly visible.As can be seen from Eq. ( 65), the Knudsen layer amplitude in temperature and the normal stress σ are both determined through the integration constant A. At a larger Knudsen number, e.g., at Kn = 0.235 as shown in Fig. 6, the amplitude A becomes so large that R13 fails to predict the inversion of the temperature profile, although the conductive heat flux, q 0 , is inverted indeed, as can be seen in Fig. 7. NSF, on the other hand, describes the temperature profile fairly accurately but cannot describe the normal stress contribution.
The variation of mass flow J 0 and convective heat flux q 0 with the Knudsen number is shown in Fig. 7. Compared to Sec.IV C 1, we now have a larger mass flux and negative heat flux, that is, mass flux J 0 and total energy flux Q 0 point from the warm towards the cold liquid layer, but the conductive heat flux q 0 points in the opposite direction.As the figure shows, the boundary conditions for R13 (red) yield some deviation.Looking at the relative errors, we find that for Kn ≤ 0.4, corrected NSF predicts the mass flow with less than 0.5% deviation from the DSMC solutions and the heat flux q 0 with not FIG.6. Inverted temperature profile at Kn = 0.235: Temperature and stress profiles for ∆θ = 0.01 and ∆p = 0.075, for DSMC simulation (symmetrized; green, dashed), corrected NSF (blue, dashed), uncorrected NSF (black, dotteddashed), and R13 (red, continuous).Note that the liquid temperatures at the left and right are θ 0 = 1.01 and θ 1 = 0.99, respectively.FIG. 7. Mass and heat transfer for inverted profile: Mass flow J 0 and heat flux q 0 as a function of the Knudsen number Kn, for ∆θ = 0.01 and ∆p = 0.075: DSMC simulation (green dots), corrected NSF (blue, dashed), uncorrected NSF (black, dotted-dashed), and R13 (red).more than 10% deviation.The errors are larger for R13, with up to 2.1% for mass flux and up to 30% for heat flux; the total energy flux Q 0 = 5 2 J 0 + q 0 , however, is predicted within 1.5%.Uncorrected NSF exhibits the largest overall errors, up to 7.5% for the mass flux and total energy flux.

V. DISCUSSION AND CONCLUSIONS
In summary, our results indicate that the R13 equations with evaporation boundary conditions yield a qualitative description of the simple evaporation and condensation problems discussed here.In detail, however, they do not give perfect agreement.
From the shown data, one might conclude that the NSF equations with corrected boundary conditions give a more accurate description even for larger Knudsen numbers, but this is misleading.It must be kept in mind that flow problems in rarefied gases show a rich array of features, such as Knudsen layers, jump and slip, thermal stresses, thermal transpiration flow, heat flux without temperature gradient, etc., which are not accessible to the NSF theory, see Refs. 4 and 5 and references therein.The rich interplay between the different rarefaction effects comes into play in particular for multi-dimensional problems 8,37 but is mainly lost in the simple one-dimensional problems discussed above.Since R13 includes the rarefaction effects, and NSF does not, their predictions differ greatly in multidimensional problems.The good predictions of NSF for the problems discussed above are to some extent accidental, due to simple geometry.A case in point is the normal stress that R13 can reproduce in good accuracy but which is always zero for NSF.
By derivation, the NSF equations are valid only for rather small Knudsen numbers, and their good predictions for mass and energy fluxes for the problems considered here must be considered as circumstantial.Moreover, we note that the NSF interface conditions include fitted resistivities r αβ , while the R13 boundary conditions are taken directly from the derivation of kinetic theory, without any ad hoc fitting coefficients.The NSF equations with uncorrected resistivities (49) give significantly larger errors than R13, which must be considered as a fitting-free improvement of the NSF equations.
Just as adjustment of interface resistivities gives a marked improvement for NSF, it is expected that the R13 predictions can be improved by adjusting the interface conditions.In Ref. 42, we introduced an Onsager matrix in the kinetic boundary conditions for R13 at non-absorbing walls, which lead to reasonable improvement, see also Ref. 53.A similar procedure is possible here and will be explored in the future.

APPENDIX A: LIQUID PHASE AND CONSERVATION LAWS AT THE INTERFACE
In the bulk of the paper, we have considered only the gas phase (vapor), while the liquid phase only appeared in the interface conditions through its interface temperature θ L and the corresponding saturation pressure p sat (θ L ).For the problems that we have in mind for the future, such as evaporation of micro-droplets in micro-channels, the liquid phase can be assumed to be an incompressible liquid of constant mass density ρ L , with constant specific heat c L .The transport equations for the liquid phase then are simply the incompressible Navier-Stokes-Fourier equations, which read (sub-/super-script L denotes the liquid) with the respective laws for stress and heat flux, Here, µ L and κ L are the liquid's viscosity and heat conductivity, respectively.
The conservation laws for mass, momentum, and energy for the gas (2) and the vapor (A1) hold in the respective bulk phases.For a complete description, we also require the conservation laws at the interface.The easiest model of the interface is a massless interface without surface tension and surface energy.As before, we consider an interface that moves with the velocity v I i and is characterized by the normal vector n i that points from the liquid towards the gas.For such an interface, the conservation laws of mass, momentum, and energy declare the continuity of the normal fluxes for an observer resting on the interface.With the liquid properties on the left, and the gas properties on the right, these can be written as 54 ρ L v L j − v I j n j = ρ v j − v I j n j , (A3) where h = u + p ρ denotes the enthalpy.Surface tension effects must be added for curved surfaces, such as droplets or bubbles. 54t sufficiently low pressures, liquid enthalpy can be approximated as h L = c L θ − h 0 (A6) with an enthalpy constant h 0 that must be carefully chosen in order to incorporate the latent heat into the model.Indeed, in kinetic theory, the enthalpy of the (monatomic) vapor is Note that we have absorbed the gas constant R into temperature as θ = RT , hence the specific heats c L and c p = 5 2 are dimensionless; p sat (θ 0 ) is the saturation pressure at reference temperature θ 0 .

APPENDIX B: DETAILS OF DSMC CALCULATIONS
Reference flow field profiles, both for evaporation into a half space and for the evaporation/condensation between two parallel plates, have been obtained by DSMC simulations of a monatomic vapor.The code is based on a quite standard implementation of Koura's null-collision estimator of collision rate. 55The presented profiles are based on hard spheres simulations.A smaller number of computationally more expensive Maxwell molecule simulations showed very little variation with respect to hard spheres profiles, once the two model cross sections have been tuned on the same viscosity at the reference temperature, for the problem at hand.It should be observed that in both problems, temperature variations are extremely small; hence, the different viscosity temperature exponents of the two collision models are not seen in the simulations.
Boundary conditions at evaporating/condensing liquid surfaces have been assigned by assuming the unit evaporation/condensation coefficient and Maxwellian distribution function of evaporating molecules, characterized by the liquid surface temperature θ L and the saturated vapor density n sat (θ L ).In the case of the evaporation into a half space, the boundary condition at infinity requires some care since the asymptotic Maxwellian state is not known but it is a result of the simulation.Taking into account that the density and temperature at infinity are determined by the downstream Mach number (or flow speed), the correct equilibrium state can be determined by the boundary condition formulated in Ref. 56.
In order to match the linearized R13 regime, DSMC has been used to study small departures from equilibrium.The one-dimensional computational domain has been divided into a number of cells of the same size, not larger that λ 0 /10, λ 0 being the equilibrium mean free path.Not less than 3000 particles per cell have been used, in order to capture small deviations from equilibrium.Typically, each simulation has been started from a given spatially uniform equilibrium state.Sampling of macroscopic quantities has been started after the onset of steady flow conditions, the sampling time duration being determined by the estimated relative statistical error.The latter does not exceed 0.01.The highest computing time has been equal to about 26 h on a single processor in a case in which the Knudsen number has been set equal to 1/100, the number of particles is equal to 3 × 10 6 , and the number of time steps is equal to 3.2 × 10 5 .Smaller or much smaller computational effort is required for larger Knudsen numbers because of the smaller particle number and shorter transient duration.
H.S. and A.B. are supported by the Natural Sciences and Engineering Research Council (NSERC), A.S.R. is supported by the EPSRC Programme Grant No. EP/N016602/1, and the collaboration with A.F. is supported by the National Group of Mathematical Physics (INDAM-GNFM).

2 θ
= c p θ (A7) so that the heat of vaporization ish LV (θ) = h − h L = h 0 − c L − c p θ. (A8)Due to the simplicity of the model, only one value of the heat of vaporization is required to fit the constant h 0 .With a reference value h 0 LV = h LV (θ 0 ), we have h 0 = h 0 LV + c L − c p θ 0 , hence the enthalpy of the liquid and the enthalpy of vaporization areh L = c L (θ − θ 0 ) + c p θ 0 − h 0 LV , (A9) h LV (θ) = h 0 LV − c L − c p (θ − θ 0 ) .(A10)In Ref.49, it is shown that this simple model-vapor as ideal gas and liquid as incompressible ideal liquid, both with constant specific heats-leads to an explicit expression for saturation pressure that in the present notation reads as p sat (θ) = p sat (θ 0 ) exp c L − c p 1 −