# An Efficient Moment Method for Modelling Nanoporous Evaporation

Thomas C. De Fraja<EMAIL_ADDRESS>Mathematics Institute, University of
Warwick, UK  School of Engineering, University of Warwick, UK  Anirudh S. Rana
Department of Mathematics, Birla Institute of Technology, India  Ryan Enright
Nokia Bell Labs, New Jersey, USA  Laura J. Cooper  Mathematics Institute,
University of Warwick, UK  Duncan A. Lockerby<EMAIL_ADDRESS>School
of Engineering, University of Warwick, UK  James E. Sprittles
<EMAIL_ADDRESS>Mathematics Institute, University of Warwick, UK

###### Abstract

Thin-film-based nanoporous membrane technologies exploit evaporation to
efficiently cool microscale and nanoscale electronic devices. At these scales,
when domain sizes become comparable to the mean free path in the vapour,
traditional macroscopic approaches such as the Navier-Stokes-Fourier (NSF)
equations become less accurate, and the use of higher-order moment methods is
called for. Two higher-order moment equations are considered; the linearised
versions of the Grad 13 and Regularised 13 equations. These are applied to the
problem of nanoporous evaporation, and results are compared to the NSF method
and the method of direct simulation Monte Carlo (i.e. solutions to the
Boltzmann equations). Linear and non-linear versions of the boundary
conditions are examined, with the latter providing improved results, at little
additional computational expense, compared to the linear form. The outcome is
a simultaneously accurate and computationally efficient method, which can
provide simulation-for-design capabilities at the nanoscale.

Evaporative cooling, nanoscale flow, gas kinetic effects, higher order moment
methods

##  I Introduction

Recent developments in technologies, such as laser diodes, power amplifiers,
and the upcoming 6G, require more sophisticated thermal-management systems
than traditional cooling techniques, as well as more compact solutions. For
instance gallium nitride based power amplifiers generate hot spots on the
micro scale, with heat fluxes of over 1kW/cm2, temperatures reaching over
180∘C, and temperature fluctuations of around 40∘C over the space of a few
microns [1]. Such high temperatures and heat fluxes, if unmanaged, lead to a
decreased device efficiency and short operating lifetimes.

One promising emerging technology aims to improve heat dissipation
capabilities, while also being substantially more compact than traditional
cooling techniques. It uses a thin membrane (thickness $\sim$1$\mu$m) covered
in nanometre sized pores (pore diameter $\sim$100nm) to transfer heat away via
evaporation [2, 3]. An advantage of this cooling mechanism is that the
evaporated mass lost is naturally replaced by the capillary action of the
liquid (i.e. via wicking). This saves considerable power and space, as
substantial pumping is no longer needed.

Different factors can affect the evaporative capability of nanoporous
membranes, including the size of the nanopores [4], the thickness of the
membrane [4], and the working liquid used [2]. The liquid could also have an
effect on the shape of the meniscus formed inside the nanopores, which will
affect the evaporative surface area, potentially affecting evaporation rates
[3]. Optimising a system over such a large parameter space can be challenging
to explore experimentally, costly in both time and money, and difficult to
measure due to the small spatio-temporal scales of interest. Computational
models offer a cheaper alternative to experiments that enable access to ‘real
time’ flow characteristics. However, at the nanoscale, additional physics
renders the widely used Navier-Stokes-Fourier model inaccurate, motivating the
use of both molecular dynamics [5, 6] and methods based on the full Boltzmann
equation [7, 8, 9]. Unfortunately, whilst highly accurate, such approaches can
be computationally expensive, particularly when simulating device-scale
dynamics.

Besides the Boltzmann equation, rarefaction effects that are beyond the
resolution of the NSF system can be predicted by extended macroscopic moment
equations [10, 11, 12]. The moment equations have the favourable form of a set
of partial differential equations describing the evolution of macroscopic
quantities, such as mass density, temperature, velocity, heat flux, stress
tensor and so on, defined as moments of the distribution function. These
equations are obtained by an asymptotic reduction of the Boltzmann equation at
different levels of approximation. The moment method was introduced to gas
kinetic theory by Grad [13], who expressed the distribution function in terms
of Hermite polynomials. More recently, the regularisation of Grad’s 13-moment
(G13) equations have been obtained by Struchtrup and Torrilhon [14]. The
regularised 13-moment (R13) equations introduce additional terms to the G13
equations that overcome various deficiencies, such as the prediction of sub-
shocks at high Mach number (Ma $\geqslant$ 1.65). Notably, other approaches
for reducing the Boltzmann equation have been proposed [15, 16], but these
focus on non-linear/hyperbolic regimes beyond the remit of this article.
Importantly, to capture micro/nano-flows, where surface effects can be
dominant, Gu and Emerson [12] presented a set of wall boundary conditions for
the R13 equations derived from the Maxwell accommodation model, with
corrections to this work presented by Torrilhon and Struchtrup [11].
Subsequently, the R13 equations have been considered for canonical boundary-
value problems, such as planar and cylindrical Couette and Poiseuille flows
[17, 18], transpiration flows and gas flow past a sphere [19], among many
others, in one-, two- and three-dimensional numerical simulations [20, 21, 22,
12, 23, 24, 25, 20, 26]. Building on these successes, we are now able to focus
these equations on a specific technological application.

Section II introduces the equations, of which we consider three sets (or
‘levels’ of accuracy/complexity): the Navier-Stokes-Fourier equations (NSF);
the Grad 13 equations; and the regularised Grad 13 equations, known as the R13
equations. These are all moment-based approximations of the Boltzmann
equation, each more accurate than the previous [27]. We compare solutions of
these equations to computational results obtained by John _et al._ [28, 29],
who utilised the standard Direct Simulation Monte Carlo (DSMC) method of
solving the Boltzmann equation.

In Sections IV, V and VI, we present the results of nanoporous membrane
simulations based on the moment equations being solved using the finite
element method. The different models are compared to analytic results and DSMC
for both a fully linearised model in Section IV, and a model incorporating
non-linear boundary conditions in Section V; the latter providing the most
accurate representation of the DSMC simulation of nanoporous evaporation, and
a marked improvement over the fully linear model (see e.g. Figure 7). Section
VI analyses the effect of scale on the membrane simulations, extending the
work done in sections IV and V to higher Knudsen numbers.

Finally, in Section VII, conclusions and discussion are presented.

##  II The Moment Equations

The moment equations are derived as approximations of the Boltzmann equation
by Grad’s moment method [13]. Details of the derivation are well documented,
and are outlined in Appendix A. Index notation is used throughout, with
indices used for the spatial dimension. The Grad 13 and R13 equations are a
set of PDEs in terms of the following moments; the density $\varrho$, the
velocity $\mathcal{V}_{i}$, temperature given in specific energy units
$\Theta$, the trace-free and symmetric stress tensor $\Sigma_{ij}$, and the
heat flux $\mathcal{Q}_{i}$. For a formal introduction to these moments, see
the Appendix A and [10]. This forms a set of 13 moments, which we denote

$\Phi^{[13]}=\varrho\left\\{1,\mathcal{V}_{i},\Theta,\Sigma_{ij},\mathcal{Q}_{i}\right\\}.$

The governing equations for the two 13-moment systems, describing the
evolution of the set of moments $\Phi^{[13]}$, consist of the conservation
laws;

$\displaystyle D_{t}\varrho+\varrho\partial_{k}\mathcal{V}_{k}=0,$ (1)
$\displaystyle\varrho
D_{t}\mathcal{V}_{i}+\partial_{i}(\varrho\Theta)+\partial_{k}\Sigma_{ik}=G_{i},$
$\displaystyle\frac{3}{2}\varrho
D_{t}\Theta+\varrho\Theta\partial_{k}\mathcal{V}_{k}+\partial_{k}\mathcal{Q}_{k}+\Sigma_{kl}\partial_{k}\mathcal{V}_{l}=0,$

with the additional balance equations for stress and heat flux;

$\displaystyle
D_{t}\Sigma_{ij}+\Sigma_{ij}\partial_{k}\mathcal{V}_{k}+2\varrho\Theta\partial_{\langle
i}\mathcal{V}_{j\rangle}+2\Sigma_{k\langle
i}\partial_{k}\mathcal{V}_{j\rangle}$ (2)
$\displaystyle+\frac{4}{5}\partial_{\langle
i}\mathcal{Q}_{j\rangle}+\partial_{k}\mathcal{M}_{ijk}=-\varrho\Theta\frac{\Sigma_{ij}}{\mu},$
$\displaystyle
D_{t}\mathcal{Q}_{i}+\frac{5}{2}\Sigma_{ik}\partial_{k}\Theta-\Sigma_{ik}\Theta\partial_{k}\ln\varrho+\Theta\partial_{k}\Sigma_{ik}+\frac{7}{5}\mathcal{Q}_{i}\partial_{k}\mathcal{V}_{k}$
(3)
$\displaystyle+\frac{7}{5}\mathcal{Q}_{k}\partial_{k}\mathcal{V}_{i}+\frac{2}{5}\mathcal{Q}_{k}\partial_{i}\mathcal{V}_{k}+\frac{1}{2}\partial_{k}\mathcal{R}_{ik}+\frac{1}{6}\partial_{i}\mathfrak{D}+\mathcal{M}_{ikl}\partial_{k}\mathcal{V}_{l}$
$\displaystyle-\frac{\Sigma_{ik}}{\varrho}\partial_{l}\Sigma_{kl}+\frac{5}{2}\varrho\Theta\partial_{i}\Theta=-\frac{5}{2}\varrho\Theta\frac{\mathcal{Q}_{i}}{\kappa},$

where $\partial_{i}=\frac{\partial}{\partial x_{i}}$,
$D_{t}=\partial_{t}+v_{k}\partial_{k}$ is the material derivative, $G_{i}$ are
external forces, $\mu$ is the dynamic viscosity, and $\kappa=\frac{15}{4}\mu$
is the heat conductivity. Indices in angular brackets denote the trace-free-
symmetric part of a tensor [10]. We also assume the gas to be ideal, with the
pressure $\mathcal{P}=\varrho\Theta$.

The five-moment NSF system for the moments
$\Phi^{[5]}=\varrho\left\\{1,\mathcal{V}_{i},\Theta\right\\}$ applies closure
to equations (1), with well-known constitutive laws for stress and heat flux

$\displaystyle\Sigma_{ij}=-2\mu\partial_{\langle
i}\mathcal{V}_{j\rangle},\quad\mathcal{Q}_{i}=-\kappa\partial_{i}\Theta,$ (4)

which give the Navier-Stokes equations and the heat equation (based on
Fourier’s law). While the NSF system only requires the five moments
$\Phi^{[5]}$, it is worth noting that when solving this system numerically, we
consider the full 13 moments $\Phi^{[13]}$, in order to retain a first-order
system, replacing (2)-(3) by (4).

#### Closure of the 13 Moment System

The equations (1)-(3) give 13 equations for the 13 moments $\Phi^{[13]}$, but
there are the additional moments $\mathcal{M}_{ijk},\mathcal{R}_{ij}$ and
$\mathfrak{D}$. For the Grad 13 system, the higher moments
$\mathcal{M}_{ijk},\mathcal{R}_{ij}$ and $\mathfrak{D}$ are set to zero [13,
30, 10] in order to close (1)-(3).

The R13 system is derived via the order-of-magnitude method as an
approximation to the Boltzmann equation from its infinite set of corresponding
moment equations [31, 10]. In the order-of-magnitude method, first the leading
order of all moments is determined by means of the Chapman-Enskog expansion
(in Knudsen number) and then the moment equations are systematically reduced
by cancelling terms of higher order. The R13 equations are obtained via third
order closure, i.e., an order equivalent to Super-Burnett (a model which is
derived via a direct Chapman-Enskog expansion of the Boltzmann equation [10]).
Notably, over the years, a number of different variants of these equations has
been suggested with differences occurring in the non-linear terms of the
higher moments, which will not concern us here, as we will focus on the
linearised form.

The closure for the R13 system is given by

$\displaystyle 2\mu\Theta\partial_{\langle
i}\frac{\Sigma_{jk\rangle}}{\mathcal{P}}-\frac{4}{3}\frac{\Sigma_{\langle
ij}\mathcal{Q}_{k\rangle}}{\mathcal{P}}$ $\displaystyle=-\mathcal{M}_{ijk},$
(5) $\displaystyle\frac{24}{5}\mu\Theta\partial_{\langle
i}\frac{\mathcal{Q}_{j\rangle}}{\mathcal{P}}-\frac{64}{25}\frac{\mathcal{Q}_{\langle
i}\mathcal{Q}_{j\rangle}}{p}-\frac{20}{7}\frac{\Sigma_{k\langle
i}\Sigma_{j\rangle k}}{\varrho}$ $\displaystyle=-\mathcal{R}_{ij},$
$\displaystyle
12\mu\Theta\partial_{k}\frac{\mathcal{Q}_{k}}{\mathcal{P}}-\frac{56}{5}\frac{\mathcal{Q}_{k}\mathcal{Q}_{k}}{\mathcal{P}}-5\frac{\Sigma_{kl}\Sigma_{lk}}{\varrho}$
$\displaystyle=-\mathfrak{D}.$

In a similar vein to the NSF equations, though the R13 equations are a closed
set of PDEs for the moments $\Phi^{[13]}$, the full set of 26 moments

$\Phi^{[26]}=\varrho\left\\{1,\mathcal{V}_{i},\Theta,\Sigma_{ij},\mathcal{Q}_{i},\mathcal{M}_{ijk},\mathcal{R}_{ij},\mathfrak{D}\right\\}$

are considered for numerical simulation.

In line with the order-of-magnitude method, at the zeroth order closure the
Euler equations are recovered, while the second order closure leads to the NSF
equations. The Grad 13 moment equations (for the Maxwellian molecules) were
shown to be second order (i.e., Burnett order), and R13 third order (i.e.,
super-Burnett order). The R26 equations obtained by Gu and Emerson [32] are
fifth order accurate.

Apart from the accuracy of the macroscopic theories, the equations also differ
in their mathematical and physical nature. The Grad 13 system, due to their
hyperbolic character, produces unphysical sub-shocks for the flows with Mach
number, Ma$\geq 1.65$ [33]. On the other hand, the NSF and the R13 equations
give smooth shock structures for all Mach number [33, 34]. Moreover, the
results predicted by the R13 theory are closer to kinetic theory for Ma
$\lesssim 5$ [34]. The Burnett equations are known to show linear unstablites
for time-dependent problems [35] while the NSF, Grad 13 and R13 equations are
stable [31].

###  II.1 Boundary Conditions

The boundary conditions for the 13-moment systems can be derived by continuity
of fluxes of moments of the distribution function at the boundary [36]. The
liquid is assumed to be in equilibrium, with temperature $\Theta_{L}$ and
saturation pressure $\mathcal{P}_{s}$. The probability that a particle
evaporates/condenses at a liquid-vapour portion of the boundary is
$\vartheta$. The accommodation coefficient is $\chi$, defined as the
probability that a particle is thermalised. We consider only 2-dimensional
flows, leaving us with the 2-dimensional boundary coordinate system, with
normal direction denoted $n$, and tangential direction denoted $t$. The normal
coordinate representation of a tensor $T_{i_{1},...,i_{m}}$ is

$T_{n...n}=T_{i_{1}...i_{m}}n_{i_{1}}...n_{i_{m}},$

and similarly for tangential and mixed parts. The effective pressure is

$\Pi=\mathcal{P}+\frac{1}{2}\Sigma_{nn}-\frac{1}{120}\frac{\mathfrak{D}}{\Theta}-\frac{1}{28}\frac{\mathcal{R}_{nn}}{\Theta}.$

The resulting boundary conditions [36] are given in (6)-(11), with the
evaporation mass flow condition

$\varrho\mathcal{V}_{n}=\frac{\vartheta}{2-\vartheta}\sqrt{\frac{2}{\pi}}\left[\frac{\mathcal{P}_{s}\left(\Theta_{L}\right)}{\sqrt{\Theta_{L}}}-\frac{\Pi}{\sqrt{\Theta}}\right],$
(6)

tangential slip condition

$\displaystyle\Sigma_{tn}=-\frac{\vartheta+\chi(1-\vartheta)}{2-\vartheta-\chi(1-\vartheta)}\sqrt{\frac{2}{\pi\Theta}}\left[\Pi\mathcal{V}_{t}+\frac{1}{5}\mathcal{Q}_{t}+\frac{1}{2}\mathcal{M}_{tnn}\right]-\varrho\mathcal{V}_{n}\mathcal{V}_{t},$
(7)

and the normal heat transfer condition

$\displaystyle\mathcal{Q}_{n}=$
$\displaystyle-\frac{\vartheta+\chi(1-\vartheta)}{2-\vartheta-\chi(1-\vartheta)}\sqrt{\frac{2}{\pi\Theta}}\left[2\Pi\left(\Theta-\Theta_{L}\right)-\frac{\Pi}{2}\mathcal{V}_{t}^{2}+\frac{1}{2}\Theta\Sigma_{nn}+\frac{\mathfrak{D}}{15}+\frac{5}{28}\mathcal{R}_{nn}\right]$
(8)
$\displaystyle+\left[\frac{1}{2}\left(\mathcal{V}_{t}^{2}-\Theta_{L}\right)-\frac{5}{2}\left(\Theta-\Theta_{L}\right)\right]\varrho\mathcal{V}_{n}.$

The R13 equations require additional boundary conditions on top of (6)-(8),
obtained by considering higher moment fluxes [36], and are given by

$\displaystyle\mathcal{M}_{nnn}=$
$\displaystyle\frac{\vartheta+\chi(1-\vartheta)}{2-\vartheta-\chi(1-\vartheta)}\sqrt{\frac{2}{\pi\Theta}}\left[\frac{2}{5}\Pi\left(\Theta-\Theta_{L}\right)-\frac{3}{5}\Pi\mathcal{V}_{t}^{2}-\frac{7}{5}\Theta\Sigma_{nn}+\frac{\mathfrak{D}}{75}-\frac{1}{14}\mathcal{R}_{nn}\right]$
(9)
$\displaystyle-\frac{2}{5}\left[\Theta_{L}+\frac{3}{2}\mathcal{V}_{t}^{2}\right]\varrho\mathcal{V}_{n},$
$\displaystyle\mathcal{M}_{ttn}=$
$\displaystyle-\frac{\vartheta+\chi(1-\vartheta)}{2-\vartheta-\chi(1-\vartheta)}\sqrt{\frac{2}{\pi\Theta}}$
(10)
$\displaystyle\times\left[\Theta\Sigma_{tt}-\Pi\mathcal{V}_{t}^{2}+\frac{\mathcal{R}_{tt}}{14}+\frac{1}{5}\Pi\left(\Theta-\Theta_{L}\right)+\frac{1}{5}\Pi\mathcal{V}_{t}^{2}-\frac{1}{5}\Theta\Sigma_{nn}+\frac{\mathfrak{D}}{150}\right]$
$\displaystyle+\frac{1}{5}\left[4\mathcal{V}_{t}^{2}+\Theta_{L}\right]\varrho\mathcal{V}_{n},$
$\displaystyle\mathcal{R}_{tn}=$
$\displaystyle\frac{\vartheta+\chi(1-\vartheta)}{2-\vartheta-\chi(1-\vartheta)}\sqrt{\frac{2}{\pi\Theta}}$
(11)
$\displaystyle\times\left[\Pi\Theta\mathcal{V}_{t}-\frac{11}{5}\Theta\mathcal{Q}_{t}-\frac{1}{2}\Theta\mathcal{M}_{tnn}-\Pi\mathcal{V}_{t}^{3}+6\Pi\mathcal{V}_{t}\left(\Theta-\Theta_{L}\right)\right]$
$\displaystyle+\left[7\left(\Theta-\Theta_{L}\right)+\Theta_{L}-\mathcal{V}_{t}^{2}\right]\varrho\mathcal{V}_{t}\mathcal{V}_{n}.$

Notably, the mass-flux condition (6) is the generalisation of the classical
Hertz-Knudsen-Schrage law to the higher moment equations, with the vapour
pressure $\mathcal{P}$ being replaced by the effective vapour pressure $\Pi$,
and the pre-factor being twice as large [37]. The condition (7) can be seen as
a generalised slip condition, relating the tangential velocity
$\mathcal{V}_{t}$ with the shear stress $\Sigma_{tn}$, c.f. the Navier slip
model [38, 39]. Similarly, (8) relates the temperature jump
$\Theta-\Theta_{L}$ with the heat flux $\mathcal{Q}_{n}$, extending the
temperature continuity condition at the boundary. Equations (9)-(11) are
interface conditions for the higher moments. For non-evaporating surfaces
(e.g. at a solid boundary), where $\vartheta=0$, equation (6) reduces to the
impermeability condition $\mathcal{V}_{n}=0$, and reproduces the well known
jump coefficient $\frac{\chi}{2-\chi}$ for the remaining equations (7)-(11),
and removes additional contributions of mass flux, giving the vapour-solid R13
boundary conditions [40].

In dealing with fluid flow through micro-devices, one is faced with the
question of which model to use, which boundary conditions to apply and how to
proceed to obtain solutions to the problem at hand. Surface effects dominate
in small devices: the surface-to-volume ration for a device with character
length of 1 m is 1 m-1, while that for microelectromechanical systems (MEMS)
device of size 1 $\mu$m is $10^{6}$ m-1. The million fold increase in surface
area relative to the mass of the minute device substantially affects the
transport of mass, momentum and energy through the surface [41]. As such,
proper modelling of the boundary conditions is essential, and it is often the
case, particularly when approaching the transition regime, that much of the
extra behaviour captured by the higher moment equations is dominated by
modifications to the boundary conditions [42], as will be confirmed in Section
E. Still, gas rarefaction also leads to novel bulk phenomena, such as the
Knudsen layer, a region of few mean free path lengths thickness where non-
equilibrium effects dominates. A detailed examination of the moment models
shows that in higher order moment equations the Knudsen layer appears as
superpositions of several exponential layers with different coefficients in
their exponents [26, 17]. Notably, the NSF and Grad 13 theories cannot capture
the Knudsen layer, whilst the R13 theory does, with more moments (e.g. R26)
giving a better description of it [17, 43].

Over the years, a variety of methods have been devised for the simulation of
different moment-based systems [44, 45, 46] and, in particular, the R13
equations have been successfully applied to a range of boundary value problems
[17, 20]. Different numerical schemes have been developed for them, for
example finite difference [21], finite volume [22, 12], finite element method
[23, 24, 25], and mesh-free methods such as the method of fundamental
solutions [20, 26, 47].

###  II.2 Linearisation

We simplify the above equations and boundary conditions by linearising about a
homogeneous state described by the temperature of the liquid $\Theta_{L}$ and
saturation pressure $\mathcal{P}_{s}$. The saturation density $\varrho_{s}$ is
such that $\mathcal{P}_{s}=\varrho_{s}\Theta_{L}$. Therefore, we shall
consider small perturbations from the equilibrium and introduce characteristic
scales (i.e. non-dimensionalise) such that [14]

$\displaystyle\varrho$
$\displaystyle=\varrho_{s}(1+\rho),\quad\Theta=\Theta_{L}(1+\theta),\quad\mathcal{P}=\mathcal{P}_{s}(1+\rho+\theta),$
(12) $\displaystyle\mathcal{V}_{i}$
$\displaystyle=\sqrt{\Theta_{L}}v_{i},\quad\Sigma_{ij}=\varrho_{s}\Theta_{L}\sigma_{ij},\quad\mathcal{Q}_{i}=\varrho_{s}\sqrt{\Theta_{L}}q_{i},$
$\displaystyle\mathfrak{D}$
$\displaystyle=\varrho_{s}\Theta_{L}^{2}\Delta,\quad\mathcal{R}_{ij}=\varrho_{s}\Theta_{L}^{2}R_{ij},\,\,\,\mathcal{M}_{ijk}=\varrho_{s}\sqrt{\Theta_{L}}^{3}m_{ijk}.$

The linear (dimensionless) moments

$\varphi^{[26]}=\rho\left\\{1,v_{i},\theta,\sigma_{ij},q_{i},\Delta,R_{ij},m_{ijk}\right\\}$

are taken to be small, thus only linear terms in these will appear in the
equations. In this regime, the ideal gas law becomes $p=\rho+\theta.$

In Sections IV-VI, we assess parameters in comparison to DSMC results, for
which we require the dimensional moments $\Phi^{[26]}$. We note that the
equations are solved both analytically and numerically in terms of the linear
non-dimensional moments $\varphi^{[26]}$, after which, for comparison to DSMC,
they are plotted in terms of the dimensional moments $\Phi^{[26]}$.

We non-dimensionalise the dependent variables;

$\displaystyle x_{i}=L\hat{x}_{i},\quad t=\frac{L}{\sqrt{\Theta_{L}}}\hat{t},$

where $L$ is the characteristic length scale. The hat notation for the
dimensionless-dependent variables is henceforth dropped. The Knudsen number is
defined as

$\displaystyle\mathrm{Kn}=\frac{\mu\sqrt{\Theta_{L}}}{\mathcal{P}_{s}L}=\frac{\lambda}{L}$

where $\lambda$ is the mean free path of a particle. Therefore, in this setup,
variations in length (e.g. pore diameter) are obtained by varying the Knudsen
number.

The saturation pressure and liquid temperature can also be considered as
perturbations from their spatially homogeneous equilibrium, with

$\displaystyle\mathcal{P}_{s}=\mathcal{P}_{0}(1+p_{s}),$
$\displaystyle\Theta_{L}=\Theta_{0}(1+\theta_{L}).$

For the most part we will consider constant liquid temperature and saturation
pressure, thus $p_{s}=\theta_{L}=0$. However, these terms will be left in the
boundary conditions to more easily indicate pressure and temperature jumps. We
also consider stationary solutions, i.e. $\partial_{t}=0$.

Substitution into (1)-(11) yields the dimensionless and linearised
conservation laws

$\displaystyle\partial_{k}v_{k}$ $\displaystyle=0,$ (13)
$\displaystyle\partial_{i}\theta+\partial_{i}\rho+\partial_{k}\sigma_{ik}$
$\displaystyle=G_{i},$ $\displaystyle\partial_{k}q_{k}$ $\displaystyle=0,$

and the expressions for stress and heat flux

$\displaystyle-\frac{\sigma_{ij}}{\text{Kn}}=\frac{4}{5}\partial_{\langle
i}q_{j\rangle}+2\partial_{\langle i}v_{j\rangle}+\partial_{k}m_{ijk},$ (14)
$\displaystyle-\frac{2}{3}\frac{q_{i}}{\text{Kn}}=\frac{5}{2}\partial_{i}\theta+\partial_{k}\sigma_{ik}+\frac{1}{2}\partial_{k}R_{ik}+\frac{1}{6}\partial_{i}\Delta.$
(15)

The linearised NSF expressions for stress and heat flux are

$\displaystyle\sigma_{ij}=-2\text{Kn}\partial_{\langle i}v_{j\rangle},\quad
q_{i}=-\frac{15}{4}\text{Kn}\partial_{i}\theta.$ (4′)

The higher moments $m_{ijk},R_{ij},\Delta$ all vanish for the Grad 13 system,
while the closure of the R13 system is given by

$\displaystyle-\frac{m_{ijk}}{\text{Kn}}=2\partial_{\langle
i}\sigma_{jk\rangle},$
$\displaystyle-\frac{R_{ij}}{\text{Kn}}=\frac{24}{5}\partial_{\langle
i}q_{j\rangle},$ $\displaystyle-\Delta=0$ (16)

The linearised boundary conditions are given by

$v_{n}=\frac{\vartheta}{2-\vartheta}\sqrt{\frac{2}{\pi}}\left(p_{s}-\Pi-\frac{1}{2}\theta_{L}\right),$
(17) $\displaystyle\sigma_{tn}=$
$\displaystyle-\tilde{\chi}\left(v_{t}+\frac{1}{5}q_{t}+\frac{1}{2}m_{tnn}\right),$
(18) $\displaystyle q_{n}=-\tilde{\chi}$
$\displaystyle\left(2\left(\theta-\theta_{L}\right)+\frac{1}{2}\sigma_{nn}+\frac{5}{28}R_{nn}\right)-\frac{1}{2}v_{n},$
(19) $\displaystyle m_{nnn}=\tilde{\chi}$
$\displaystyle\left(\frac{2}{5}\left(\theta-\theta_{L}\right)-\frac{7}{5}\sigma_{nn}-\frac{1}{14}R_{nn}\right)-\frac{2}{5}v_{n},$
(20) $\displaystyle m_{ttn}=-\tilde{\chi}$
$\displaystyle\left(\sigma_{tt}+\frac{R_{tt}}{14}+\frac{1}{5}\left(\theta-\theta_{L}\right)-\frac{1}{5}\sigma_{nn}\right)+\frac{1}{5}v_{n},$
(21) $\displaystyle R_{tn}=\tilde{\chi}$
$\displaystyle\left(v_{t}-\frac{11}{5}q_{t}-\frac{1}{2}m_{tnn}\right),$ (22)

where

$\displaystyle\tilde{\chi}=\frac{\vartheta+\chi(1-\vartheta)}{2-\vartheta-\chi(1-\vartheta)}\sqrt{\frac{2}{\pi}},$

and the effective pressure $\Pi$ is

$\Pi=\rho+\frac{1}{2}\theta+\frac{1}{2}\sigma_{nn}-\frac{1}{28}R_{nn}.$

For all cases considered in this work, fully diffuse molecular re-emissions
are considered at the interface (i.e. $\chi$ = 1), as this is used for the
DSMC results [28, 29]. The evaporation coefficient is zero at a solid-vapour
interface, and we assume perfect evaporation at a liquid-vapour interface,
i.e. $\vartheta=1$.

For comparison, we solve the Navier-Stokes-Fourier equations with conventional
boundary conditions; that is, (18) is replaced by $v_{t}=0$ (no slip), and
(15) is replaced by $\theta=\theta_{L}=0$ (no temperature jump). The
evaporation condition (17) remains.

For the avoidance of doubt, we solve the following systems: NSF, with
conservation laws (13) and closure (4′), with the evaporation boundary
condition (17) accompanied by no slip and no temperature jump; Grad 13, with
conservation laws (13), balance equations (14)-(15) and closure
$m_{ijk}=0,R_{ij}=0,\Delta=0$, with boundary conditions (17)-(19); R13, with
conservation laws (13), balance equations (14)-(15) and closure (16), with
boundary conditions (17)-(22).

##  III Problem Formulation for Nanoporous Evaporation

The nanoporous membrane is assumed to be two-dimensional in order to allow
comparison with the DMSC results of John _et al._ [28, 29]. It is also
necessary to validate the moment equations for two-dimensional flows before
modelling a three-dimensional nanoporous membrane, where the computational
tractability of the model becomes more important.

The geometry of a single pore of the nanoporous membrane is shown in Figure 1,
with this work focussing on the unconventional modelling of the vapour flow,
where we can compare to DSMC, with future work coupling this to the liquid’s
(conventional) dynamics. The vapour part (green) is the domain of computation
where we solve the moment equations. It is a rectangle, with a circular
segment attached when considering a curved meniscus. The origin is marked
$\mathcal{O}$, and the base of the rectangle is the interface with the
nanoporous membrane at $y=0,$ and the top of the rectangle represents a far-
field at $y=H$. The height $H$ of the rectangle is set to be large enough so
that parameters decay sufficiently to an equilibrium before reaching the far-
field. Parameters at the far-field will be denoted

$\psi_{\infty}:=\psi(x,y=H).$

We assume the geometry to be periodic in $x$, with the length of periodicity –
the width of the rectangle – being the combined width of the meniscus, $L$,
and the width of the walls111The physical width of the wall would be $2W$ when
the wall from the adjacent pore is considered., $W$, either side. The porosity
$\phi$ is defined as the proportion of length consisting of the meniscus, so
$\phi=L/(2W+L)$. We take $L$ as the characteristic length scale, so that in
dimensionless parameters, we have $L=1$. A curved meniscus is modelled by
adding a segment of a circle whose corners are at $x=-\frac{1}{2}$ and
$x=\frac{1}{2}$. The meniscus shape is controlled by varying the radius of the
circle, $R$. The curvature is then $\eta=\frac{1}{2R}$, so that a flat
meniscus corresponds to $\eta=0$, and a semicircular meniscus corresponds to
$\eta=1$.

Figure 1: Geometry of a single nanopore. Two cases of meniscus curvature
$\eta$ are shown, with $\eta=0$ corresponding to the flat meniscus, and
$\eta=1$ corresponding to the semi-circular meniscus. The porosity shown here
is $\phi=0.5$, though porosity can take any value between 0 and 1.

In accordance with the DSMC approach [28, 29], we impose a constant
dimensionless saturation pressure $p_{s}$ and liquid temperature $\theta_{L}$
at the liquid-vapour interface. At the far-field, we have three unknown
parameters, $v_{i,\infty},\rho_{\infty},\theta_{\infty}$. We are free to set
one of these [48]; we set the velocity field
$v_{i,\infty}=\\{0,v_{\infty}\\}$, and solve for the other two. We also
consider the process as adiabatic, so that $q_{i,\infty}=0$.

The linearised 13 moment equations (13)-(16) are solved on the interior of the
domain. At the liquid-vapour interface, we require the boundary conditions
(17)-(19) for the Grad 13 equations, and additionally the boundary conditions
(20)-(22) for the R13 equations. At the solid-vapour interface, we have
$\vartheta=0$, so that (17) becomes $v_{n}=0$.

The equations are solved using COMSOL Multiphysics v5.6, a commercial finite-
element-method (FEM) software (COMSOL Inc, Sweden) [49]. Guided by standard
approaches to incompressible Navier-Stokes simulation [50], we use the Taylor-
Hood elements for momentum, with linear element discretizations for density
and quadratic elements for velocity (and stress). Building on this methodology
for FEM we use linear elements for the temperature, and heat flux and higher
order moments have quadratic elements. In earlier works [51, 31, 52], the
derivation of entropy consistent boundary conditions were presented for the
linearised moment equations. For the problems considered in this paper we
adopted the boundary conditions proposed in [36], and while work has been done
recently to obtain stabilised R13 FEM solvers for this system [23, 24, 25], no
instabilities were encountered here.

The height is set to be $H=5+W+5\text{Kn}$, as this was found to always leave
enough room for parameters to reach free-stream equilibrium conditions. As an
example, when $\phi=0.5$ and $\text{Kn}=0.05$, the mesh consists of 12,369
triangular domain elements and 580 boundary elements for the flat meniscus
case, and 14,115 triangular domain elements and 637 boundary elements for the
semicircular meniscus case. These meshes are shown in Figure 15 in Appendix B.

The DSMC results [28, 29] have been carried out by John et al. using SPARTA
[53], which is a highly scalable parallel open-source DSMC code [54, 55]. The
gas was assumed to be argon and the variable hard sphere (VHS) model was
employed. It was shown that the DSMC method can reproduce results of the
Boltzmann equation, and so it serves as a good benchmark for the moment method
[56, 57, 58].

##  IV Results - Linear Equations

Simulations allow us to study the dynamics of the evaporative process for
various membrane porosities and meniscus shapes. To do so, we evaluate the
free-stream parameters $\varrho_{\infty},\Theta_{\infty}$ and mass flux
$J_{\infty}=\mathcal{V}_{\infty}\varrho_{\infty}$ with respect to the free-
stream Mach number,
$\mathrm{Ma}_{\infty}=\mathcal{V}_{\infty}/\sqrt{\gamma\Theta_{\infty}}=v_{\infty}/\sqrt{\gamma\left(1+\theta_{\infty}\right)},$
where $\gamma=\frac{5}{3}$ is the specific heat ratio for a monatomic gas. Our
choice of linearisation restricts us to small Mach numbers.

#### Analytic Results for a Planar Evaporative Interface (1D Case)

The limiting case $W\to 0$ with a flat meniscus corresponds a planar
evaporative interface. Analytic (one-dimensional) results have previously been
obtained in this limit by various moment methods [36, 59], and are useful as a
benchmark for our computational results, particularly given the complexity of
the R13 PDE system, even once linearised. Here, we solve all three systems
analytically, first finding the general solution to the R13 equations, of
which the Grad 13 and NSF are seen as special cases that can be obtained upon
setting various terms to zero. The boundary condition system derived is then
solved to find far-field parameter expressions for each set of equations. The
R13 boundary condition system was previously obtained by Struchtrup _et al._
[36].

Here, we make the assumption that all parameters do not have $x$ or $z$
dependence, yielding a one-dimensional system which is straightforward to
solve, as considered by Struchtrup _et al._ [36]. The conservation equations
(13) are the same for all three models, and reduce to

$\displaystyle\partial_{y}v_{2}=0,\quad\partial_{y}\sigma_{12}=0,$
$\displaystyle\quad\partial_{y}q_{2}=0,$
$\displaystyle\partial_{y}(\rho+\theta+\sigma_{22})$ $\displaystyle=0.$

The R13 expressions for stress are

$\displaystyle\begin{aligned}
\sigma_{11}&=\underline{\frac{2}{3}\text{Kn}^{2}\partial_{yy}\sigma_{11}}-\underline{\frac{4}{15}\text{Kn}^{2}\partial_{yy}\sigma_{22}},\\\
\sigma_{12}&=\frac{2}{5}\partial_{y}q_{1}+\overline{\partial_{y}v_{1}},\\\
\sigma_{22}&=\underline{\frac{6}{5}\text{Kn}^{2}\partial_{yy}\sigma_{22}},\end{aligned}$
(23)

and the expressions for heat flux,

$\displaystyle\begin{aligned}
q_{1}&=\underline{\frac{9}{5}\text{Kn}^{2}\partial_{yy}q_{1}},\\\
q_{2}&=-\overline{\frac{15}{4}\text{Kn}\partial_{y}\theta}-\frac{3}{2}\text{Kn}\partial_{y}\sigma_{22}.\end{aligned}$
(24)

In (23) and (24), terms underlined indicate those that only appear in the R13
system, while terms overlined indicate the only (derivative) terms that appear
in the NSF system.

Integration of the above gives the general solution

$\displaystyle\begin{aligned} v_{2}&=C_{1},\qquad
v_{1}+\frac{2}{5}q_{1}=-\frac{C_{2}}{\text{Kn}}y+C_{5},\\\
\rho+\theta+\sigma_{22}&=C_{3},\qquad\frac{5}{2}\theta+\sigma_{22}=-\frac{2}{3}\frac{C_{4}}{\text{Kn}}y+C_{6}\\\
q_{2}&=C_{4},\qquad
q_{1}=C_{9}\exp\left(-\frac{\sqrt{5}}{3}\frac{y}{\text{Kn}}\right),\\\
\sigma_{11}=C_{7}\exp&\left(-\sqrt{\frac{3}{2}}\frac{y}{\text{Kn}}\right)-\frac{1}{2}C_{8}\exp\left(-\sqrt{\frac{5}{6}}\frac{y}{\text{Kn}}\right),\\\
\sigma_{12}&=C_{2},\qquad\sigma_{22}=C_{8}\exp\left(-\sqrt{\frac{5}{6}}\frac{y}{\text{Kn}}\right).\\\
\end{aligned}$ (25)

The positive exponential terms have been neglected since we may take
$y\to\infty$, and require that parameters are bounded in the far-field. This
can also be used to conclude that $C_{2}=C_{4}=0$. The Grad 13 and NSF systems
have $C_{7,8,9}=0$.

At the far field, $y\to\infty$, we have
$v_{i,\infty}=\left\\{0,v_{\infty}\right\\},q_{i,\infty}=\left\\{0,0\right\\}$.
Thus $C_{1}=v_{\infty}$ and $C_{5}=0$. Therefore, for NSF and Grad 13 we have
$v_{1}=0$, as expected.

With far-field notation, we find that
$\theta_{\infty}=\frac{2}{5}C_{6},\rho_{\infty}=C_{3}-\theta_{\infty},$ i.e.

$\displaystyle\theta$
$\displaystyle=\theta_{\infty}-\frac{2}{5}C_{8}\exp\left(-\sqrt{\frac{5}{6}}\frac{y}{\text{Kn}}\right),$
$\displaystyle\rho$
$\displaystyle=\rho_{\infty}-\frac{3}{5}C_{8}\exp\left(-\sqrt{\frac{5}{6}}\frac{y}{\text{Kn}}\right).$

At the meniscus, $y=0$, we have the boundary conditions (17)-(22). For R13
there are 5 unknowns, $\theta_{\infty},\rho_{\infty},C_{7},C_{8},C_{9},$
whereas for NSF and Grad 13 there are 2, $\theta_{\infty}$ and
$\rho_{\infty}.$ The NSF no temperature jump condition immediately gives
$\theta_{\infty}=\theta_{L}=0$ for the NSF system.

For the R13 system, equation (18) gives $C_{9}=0$, i.e. $v_{1}=q_{1}=0$, as
expected. The constant $C_{7}$ only appears in equation (21), and so can
easily be expressed in terms of $\theta_{\infty},\rho_{\infty}$ and $C_{8}$.

Boundary condition (17) is needed for the NSF system to obtain
$\rho_{\infty}$; (17) and (19) are needed for the Grad 13 system to obtain
$\theta_{\infty},\rho_{\infty}$; (17), (19) and (20) are required to solve for
$\theta_{\infty},\rho_{\infty},C_{8}$ in the R13 system. In all cases, the
boundary system may be written as

$\displaystyle\mathbf{b}=A\mathbf{U}$ ($*$)

with
$\mathbf{U}^{[\text{NSF}]}=\left\\{\rho_{\infty}\right\\}^{T},\mathbf{U}^{[\text{G}13]}=\left\\{\rho_{\infty},\theta_{\infty}\right\\}^{T},\mathbf{U}^{[\text{R}13]}=\left\\{\rho_{\infty},\theta_{\infty},C_{8}\right\\}^{T},$

$\displaystyle\mathbf{b}^{[\text{NSF}]}=\left[\begin{matrix}v_{\infty}\end{matrix}\right],$
$\displaystyle\mathbf{b}^{[\text{G}13]}=\left[\begin{matrix}v_{\infty}\\\
\frac{v_{\infty}}{2}\end{matrix}\right],$
$\displaystyle\mathbf{b}^{[\text{R}13]}=\left[\begin{matrix}v_{\infty}\\\
\frac{v_{\infty}}{2}\\\\[3.0pt] \frac{2v_{\infty}}{5}\end{matrix}\right],$
$\displaystyle
A^{[\mathrm{R}13]}=\left[\begin{matrix}-\sqrt{\frac{2}{\pi}}&-\frac{1}{\sqrt{2\pi}}&\frac{3}{5\sqrt{2\pi}}\\\
0&-\frac{4}{\sqrt{2\pi}}&\frac{3}{5\sqrt{2\pi}}\\\
0&\frac{4}{5\sqrt{2\pi}}&-\frac{39}{25}\sqrt{\frac{2}{\pi}}-\frac{\sqrt{30}}{5}\end{matrix}\right],$

$A^{[\text{G}13]}$ is the upper left $2\times 2$ sub-matrix of
$A^{[\text{R}13]}$, and $A^{[\text{NSF}]}$ is the upper left entry of
$A^{[\text{R}13]}$.

Struchtrup _et al._ [36] previously arrived at the above boundary-condition
system for the R13 equations, however explicit expressions for the solved
parameters were not given. Here, we solve ($*$ ‣ IV) to give explicit
expressions for discussion: the Grad 13 solution has

$\displaystyle\theta_{\infty}^{[\text{R}13]}=-\frac{1}{4}\sqrt{\frac{\pi}{2}}v_{\infty},$
$\displaystyle\rho_{\infty}^{[\text{R}13]}=-\frac{7}{8}\sqrt{\frac{\pi}{2}}v_{\infty},$

and R13 has

$\displaystyle\theta_{\infty}^{[\text{R}13]}$
$\displaystyle=-K\left(2\sqrt{30}\pi+18\sqrt{2\pi}\right)v_{\infty},$
$\displaystyle\rho_{\infty}^{[\text{R}13]}$
$\displaystyle=-K\left(7\sqrt{30}\pi+57\sqrt{2\pi}\right)v_{\infty},$
$\displaystyle C_{8}$ $\displaystyle=-10K\sqrt{2\pi}v_{\infty}$

where $K=\frac{1}{16\sqrt{15\pi}+120}$. The temperature expression for Grad
13, $\theta_{\infty}^{[\text{G}13]}$, is in-line with the analytic expression
for temperature obtained by Labuntsov and Kryukov [59], obtained via an
alternate approximation of the Boltzmann equation.

For classical NSF we have

$\displaystyle\theta_{\infty}^{[\text{NSF}]}=0,$
$\displaystyle\rho_{\infty}^{[\text{NSF}]}=-\sqrt{\frac{\pi}{2}}v_{\infty}.$

The expressions (25) give insight into some of the differences between the
three moment systems. All solutions for the NSF and Grad 13 systems are
spatially homogeneous, whereas for the R13 system there is an exponential
decay on the order of the Knudsen number. This represents the fact that the
R13 system is able to capture Knudsen-boundary-layer effects [14, 27], which
typically extend a few mean free paths, $\lambda$, from the boundary [60].
This phenomenon also manifests in the two-dimensional flows seen in the
numerical results, and is discussed in Section VI.

The far-field solutions for all systems, however, are independent of Knudsen
number, as the characteristic length $L$ is irrelevant in the evaluation of
far-field parameters when considering an infinitely long plane of liquid and a
one-dimensional flow. In Section VI, we see that far-field solutions are
dependent on Knudsen number when two-dimensionality is introduced.

Figure 2: Left: Far-field mass flux $J_{\infty}$, normalised with respect to
$J_{T}=\varrho_{s}\sqrt{\Theta_{L}/2\pi}$, against Mach number. Right: Far-
field density $\varrho_{\infty}$ and temperature $\theta_{\infty}$ against
Mach number. Comparison between DSMC, NSF, Grad 13 and R13, for the one-
dimensional case. Both numerical and analytic results for the moment equations
are shown, with numerical results matching analytic results well.

Numerical results for the planar case, computed in a 2D domain, match analytic
results, which gives us confidence in our computational approach. These are
shown in comparison to DSMC results in Figure 2, showing the far-field mass
flux $J_{\infty}$ normalised with respect to the total emitted (evaporative)
mass flux from a given surface, $J_{T}$, given by

$\displaystyle J_{T}=\varrho_{s}\sqrt{\frac{\Theta_{L}}{2\pi}},$

as well as the ratios $\varrho_{\infty}/\varrho_{s}$ and
$\Theta_{\infty}/\Theta_{L}$. The Grad 13 mass flux slightly closer
approximates the DSMC results for this case than that of the R13. While we
expect R13 to be a better approximation than Grad 13 in the mathematical limit
as $\text{Kn}\to 0$, as it is higher-order in Knudsen number, at finite
Knudsen numbers there is nothing to prevent Grad 13 fortuitously giving better
agreement with the DSMC than R13. Notably, from an engineering perspective,
what we will see throughout the article is that there is a general trend for
R13 to be the most accurate method, most notably in Section VI and Figure 10.

#### Numerical Results in Comparison to DSMC for the General (2D) Case

Figure 3: Flow-fields near the interface showing the normalized $y$-velocity,
$v/v_{\infty}$, and velocity streamlines for the flat meniscus cases,
$\eta=0$, with porosity $\phi=0.5$, at Kn=0.05. The far-field velocity is
$v_{\infty}=0.2172$, corresponding to $\mathcal{V}_{\infty}=60$ m/s. Three
models are shown; (a) DSMC from [29], (b) R13, and (c) NSF. Figure 4: Surface
plots near the interface showing the density ratio $\varrho/\varrho_{s}$, for
the flat meniscus case, $\eta=0$, with porosity $\phi=0.5$, at Kn=0.05. The
far-field velocity is $v_{\infty}=0.2172$, corresponding to
$\mathcal{V}_{\infty}=60$ m/s. Four models are shown; (a) DSMC from [28], (b)
R13, and (c) Grad 13, and (d) NSF.

A plot of the normalised y-component of velocity over a portion of the domain
is shown in Figure 3 for a representative case where both $v_{\infty}=0.2172$
and $\text{Kn}=0.05$ are relatively low, so we expect the models to agree; for
all cases the meniscus is flat ($\eta=0$) and the porosity is a half
($\phi=0.5$). Results for three different models are presented (DSMC from [28,
29], R13 and NSF from left to right), overlaid with velocity stream lines.
Grad 13 gives similar results to the R13 model, and is therefore not shown.
All moment models reproduce the DSMC velocity field well.

Figure 5: Normalised far-field parameters against far-field Mach number for
the full-porosity case ($\phi=1$), at Kn$=0.05$. Shown are
$J_{\infty}/J_{\infty}^{1\text{D}}$ for (a) the flat meniscus, (b) the semi-
circular meniscus; and $\varrho_{\infty}/\varrho_{s}$,
$\Theta_{\infty}/\Theta_{L}$ are shown for (c) the flat meniscus, (d) the
semi-circular meniscus. Comparison between DSMC, NSF, Grad 13 and R13. Figure
6: Normalised far-field parameters against far-field Mach number for the half-
porosity case ($\phi=0.5$), at Kn$=0.05$. Shown are
$J_{\infty}/J_{\infty}^{1\text{D}}$ for (a) the flat meniscus, (b) the semi-
circular meniscus; and $\varrho_{\infty}/\varrho_{s}$,
$\Theta_{\infty}/\Theta_{L}$ are shown for (c) the flat meniscus, (d) the
semi-circular meniscus. Comparison between DSMC, NSF, Grad 13 and R13.

The density flow field is shown in Figure 4, for all four models, for the same
case. Grad 13 and R13 show good agreement with DSMC for the variation of
density, while NSF densities are clearly inaccurate. Though the plots do not
show it, NSF sees negative values of density around the contact line between
the meniscus, wall and vapour. This is caused by a singularity on the wall
side of the wall-liquid-vapour contact line, and is discussed in Appendix D.
The singularity, however, is a local artifact, and does not affect far-field
parameters. Promisingly, the higher moment equations, Grad 13 and R13, do not
appear to exhibit singular behaviour; see Figure 18 in Appendix D for a
comparison of NSF and Grad 13.

For the assessment of far-field mass flux, it is convenient to normalise
computed mass flux values with respect to analytic results from the one-
dimensional case derived by Labuntsov and Kryukov [59]:

$J_{\infty}^{\mathrm{1D}}=v_{\infty}\left(2C_{s}^{2}+v_{\infty}^{2}-\sqrt{4C_{s}^{2}v_{\infty}^{2}+v_{\infty}^{4}}\right)\Big{/}2C_{s}^{2},$
(26)

where $C_{s}=0.6\sqrt{2\Theta_{L}}$. We again assess the ratios
$\varrho_{\infty}/\varrho_{s}$ and $\Theta_{\infty}/\Theta_{L}$.

Figure 5 show the results for the limiting case of $\phi\to 1$, corresponding
to a meniscus suspended by plate-like posts, for two cases of curvature,
$\eta=0$ and $\eta=1$. Computed far-field mass flux values for the two
13-moment systems (Grad 13 and R13), shown in Figures 5 (a) and (b), are in
relatively good agreement with DSMC compared to NSF, particularly at low Mach
numbers. For $\eta\neq 0$, where a two-dimensional flow is generated (and thus
parameters have $x$-dependence) the Knudsen number becomes relevant, and so
rarefied-gas dynamics modelling is required, even at relatively low Knudsen
numbers. For this reason, the 13-moment systems see excellent agreement for
the semicircular meniscus, while NSF deviates significantly (for example,
showing the wrong qualitative trend in Fig. 5 (b)).

Far-field density and temperature both show good agreement with DSMC for the
13-moment systems, and reasonable agreement with DSMC for the NSF system; see
Figures 5 (c) and (d). In particular, the initial slopes, corresponding to
$\mathrm{Ma}_{\infty}\approx 0$, of both densities and the temperatures for
both R13 and Grad 13 are in excellent agreement with DSMC. As a result, low-
Mach number results have near-perfect agreement in both density and
temperature. The NSF system suffers in this respect, especially for $\eta=1$,
since the initial slope in density doesn’t agree with the higher-order systems
(fig. 5 (d)). Thus, even at low-Mach numbers we have significant error for the
NSF system. This in turn means that higher Mach number results are even
further away from DSMC results.

Figure 6 shows the results for half porosity, $\phi=0.5$. We see large
discrepancies from DSMC values for computed far-field density, and thus mass
flux, for all three moment models. The introduction of walls into the
simulation produces a strongly two-dimensional flow, particularly near the
interface, as seen in Figures 3 and 4. This calls for more accuracy in both
bulk and boundary condition equations: for the one-dimensional case, only a
small amount of terms are removed upon linearisation, while in the two-
dimensional case, linearisation removes many more terms. This might explain
the discrepancies between the moment equations and DSMC seen for the linear
results in Figure 6, as DSMC inherently inlcudes the non-linearities.
Moreover, with linear boundary conditions and bulk equations, the parameters
$\rho_{\infty}$ and $\theta_{\infty}$ appear to decrease linearly with Mach
number. We saw for the one-dimensional case in the form of the analytic
solutions that these parameters are linear in $v_{\infty}$, which, for small
temperatures variations, gives an almost linear relation with Mach number. By
means of a Fourier expansion, it can be shown that for the two-dimensional
flat meniscus case far-field parameters are also linear in $v_{\infty}$, see
Appendix C. In contrast, the DSMC results for density begin to plateau as the
Mach number increases. This makes it impossible to retain good agreement for
high Mach numbers in the present set up. This is particularly prominent in
Figures 6 (c) and (d), where we see negative densities at moderate Mach
numbers, for all moment models.

The far-field temperature for the two 13-moment systems show good agreement
with DSMC values for all cases of porosity and meniscus shape, see (c) and (d)
in Figures 5 and 6. The lack of a temperature jump in the NSF system means no
temperature variation is seen at the far-field, and so accuracy is lost with
respect to the DSMC.

For all simulations in Figures 5 and 6, Grad 13 and R13 show very similar
results, as we would expect at the fairly low Knudsen number of 0.05. Much
larger differences are seen at higher Knudsen numbers, seen in Section VI. We
also see that an increase in either porosity or curvature gives an increase in
far-field mass flux, apparently due to the fact that there is increased
surface area for evaporation.

#### Discussion

The aim of this paper is to utilise the moment equations to simulate
evaporation from a nanoporous membrane. An important region for accuracy for
this is that of low Mach number flows at low-to-mid Knudsen numbers,
$0.05<\text{Kn}<1$, as DSMC is generally computationally expensive in this
regime, sometimes taking weeks to compute two-dimensional flows [61, 62, 63].
Macroscopic equations present an efficient alternative, with our simulations
from this section taking less than a minute for a given value of $v_{\infty}$
on a standard laptop. For the linear results presented above, we have accurate
results up to $\mathrm{Ma}_{\infty}\approx 0.1$, but beyond this the results
are too inaccurate to be used as a substitute for DSMC. We would therefore
like to increase the range in which we can use the moment equations and be
confident in the results.

Our 2D model is aimed at giving a simplistic, analyzable representation of a
nanoporous membrane. Since we have considered the process to be adiabatic,
energy is removed from the system solely by evaporative mass flux, and so the
accurate calculation of far-field densities is critical to understanding the
energy-dissipating capabilities of a nanoporous membrane for a given
configuration. The mass flux-driven energy dissipation would also be seen for
the working device [2], and so while our simplified 2D model will not provide
exact data corresponding to a cooling device, particularly as a confining
’ceiling’ is likely to be present in this case, it will provide general trends
which may aid the design of such a device.

To improve the accuracy of the macroscopic approach, we could consider even
more moments, since these would theoretically approximate the DSMC model to
higher accuracy [10, 64, 65]. While work has been done to understand the
number of moments required for accurate simulation of a given problem [66,
67], the computational expense and modelling complexity added for an increased
number of moments, particularly for two- and three-dimensional flows, makes
this strategy hard to justify [10, 68]. Notably, a promising new approach,
based on using a different number of moments in different regions of the
domain [69], has the potential to overcome these limitations and should be the
focus of future work in this field.

However, any linearised system of higher moment equations will still exhibit
linear relations between far-field density and $v_{\infty}$ for the flat
meniscus case: for the one-dimensional case with unitary porosity, we will
always be able to obtain a boundary condition system of the form
$\mathbf{b}=A\mathbf{U},$ where $A$ has no $v_{\infty}$-dependence, and the
linear, two-dimensional analysis done for arbitrary porosity in Appendix C
works for arbitrary moment systems. Thus high Mach number accuracy will be
sacrificed for low Mach number accuracy, or vice versa. Since the problem
seems to stem from our linearity assumption, a first-step improvement is to
introduce non-linearity in the boundary conditions, which we now discuss.

##  V Extension to non-linear boundary conditions

Here we argue the case for relaxing the linearity assumption of a component of
velocity at the boundary, and then give the modified boundary conditions.
Notably, the bulk equations remain linear. We then present the results using
these boundary conditions, including a comparison to the linear case.

At the far-field we have a constant velocity of $v_{\infty}$. When walls are
introduced at the interface ($\phi<1$), the continuity equation dictates that
there must be some value of $v_{2}$ such that $v_{2}>v_{\infty}$ at the
meniscus. Therefore, to retain accuracy at higher values of $v_{\infty}$ (and
hence higher Mach numbers) it is prudent to assume that $v_{2}$ is no longer
small along this boundary.

Insertion of (12) into (6)-(11) with the above assumption, and setting
$\chi=1$, gives the following non-linear boundary conditions:

$v_{n}(1+\rho)=\frac{\vartheta}{2-\vartheta}\sqrt{\frac{2}{\pi}}\left(p_{s}-\Pi-\frac{1}{2}\theta_{L}\right),$
(27) $\displaystyle\sigma_{tn}=$
$\displaystyle-\sqrt{\frac{2}{\pi}}\left(v_{t}+\frac{1}{5}q_{t}+\frac{1}{2}m_{tnn}\right)-v_{n}v_{t},$
(28) $\displaystyle q_{n}=-\sqrt{\frac{2}{\pi}}$
$\displaystyle\left(2\left(\theta-\theta_{L}\right)+\frac{1}{2}\sigma_{nn}+\frac{5}{28}R_{nn}\right)-\frac{1}{2}v_{n}\Big{[}\wp_{L}+5(\theta-\theta_{L})\Big{]},$
(29) $\displaystyle m_{nnn}=\sqrt{\frac{2}{\pi}}$
$\displaystyle\left(\frac{2}{5}\left(\theta-\theta_{L}\right)-\frac{7}{5}\sigma_{nn}-\frac{1}{14}R_{nn}\right)-\frac{2}{5}v_{n}\wp_{L},$
(30) $\displaystyle m_{ttn}=-\sqrt{\frac{2}{\pi}}$
$\displaystyle\left(\sigma_{tt}+\frac{R_{tt}}{14}+\frac{1}{5}\left(\theta-\theta_{L}\right)-\frac{1}{5}\sigma_{nn}\right)+\frac{1}{5}v_{n}\wp_{L},$
(31) $\displaystyle R_{tn}=\sqrt{\frac{2}{\pi}}$
$\displaystyle\left(v_{t}-\frac{11}{5}q_{t}-\frac{1}{2}m_{tnn}\right)+v_{n}v_{t},$
(32)

where $\wp_{L}=1+\rho+\theta_{L}$ and
$\Pi=\rho+\frac{1}{2}\theta+\frac{1}{2}\sigma_{nn}-\frac{1}{28}R_{nn}.$

The bulk equations remain the same for each of the three moment systems; only
the boundary conditions are changed, with (17)-(22) replaced by (27)-(32).

#### Analytic Result for a Planar Evaporative Interface (1D Case) with Non-
linear Boundary Conditions

We can again solve the system analytically for the case $\phi=1,\eta=0$. The
general solution (containing unknown constants of integration) is unchanged
from (25) from the linear case. The boundary system

$\displaystyle\mathbf{b}=A_{\text{n.l.}}\mathbf{U}$ ($**$)

changes from the linear case with the new boundary conditions. The vectors
$\mathbf{b}$ and $\mathbf{U}$ are unchanged, however the matrix $A$ becomes

$\displaystyle
A^{[\mathrm{R}13]}_{\text{n.l.}}=\left[\begin{matrix}-\sqrt{\frac{2}{\pi}}-v_{\infty}&-\frac{1}{\sqrt{2\pi}}&\frac{3}{5}v_{\infty}+\frac{3}{5\sqrt{2\pi}}\\\
-\frac{v_{\infty}}{2}&-\frac{5}{2}v_{\infty}-\frac{4}{\sqrt{2\pi}}&\frac{13}{10}v_{\infty}+\frac{3}{5\sqrt{2\pi}}\\\
-\frac{2v_{\infty}}{5}&\frac{4}{5\sqrt{2\pi}}&\frac{6}{25}v_{\infty}-\frac{39}{25}\sqrt{\frac{2}{\pi}}-\frac{\sqrt{30}}{5}\end{matrix}\right]$

and $A^{[\mathrm{G}13]}_{\text{n.l.}}$ is the upper left $2\times 2$ submatrix
of $A^{[\mathrm{R}13]}_{\text{n.l.}}$, and $A^{[\mathrm{NSF}]}_{\text{n.l.}}$
the upper left entry of $A^{[\mathrm{R}13]}_{\text{n.l.}}$. The non-linear
assumption does not introduce a non-constant system, but introduces
$v_{\infty}$-dependence on the right of ($**$ ‣ V), extending the boundary
condition system obtained by Struchtrup [36]: the matrix
$A^{[\mathrm{R}13]}_{\text{n.l.}}$ can be written as

$\displaystyle
A^{[\mathrm{R}13]}_{\text{n.l.}}=A^{[\mathrm{R}13]}_{\text{lin.}}+Cv_{\infty},\qquad
C=\left[\begin{matrix}-1&-0&\frac{3}{5}\\\\[3.0pt]
-\frac{1}{2}&-\frac{5}{2}&\frac{13}{10}\\\\[3.0pt]
-\frac{2}{5}&0&\frac{6}{25}\end{matrix}\right]$

where $A^{[\mathrm{R}13]}_{\text{lin.}}$ is the matrix from ($*$ ‣ IV) in
Section IV. This leads to a non-linear relation between the solutions and
$v_{\infty}$, given in Appendix E. The solutions exhibit asymptotic behaviour,
approaching a finite value as $v_{\infty}\to\infty$, which leads to greatly
improved results in the subsonic region we are concerned with,
$\text{Ma}_{\infty}<1$. Notably, for these solutions, densities remain
positive for all $\mathrm{Ma}_{\infty}>0$. In contrast, analytic results for
the linear boundary condition system from Section IV see negative far-field
densities for all Mach numbers above a certain value
($\mathrm{Ma}_{\infty}>0.618$ for NSF, $\mathrm{Ma}_{\infty}>0.8357$ for Grad
13, and $\mathrm{Ma}_{\infty}>0.8092$ for R13).

Numerical results once again match well with the analytic results, as now
described.

#### Numerical Results in Comparison to DSMC for the General (2D) Case with
Non-linear Boundary Conditions

We assess all parameters in the same way as for the linear case, with mass
flux normalised by $J_{\infty}^{1D}$, given by (26). For a direct comparison
with the linear case, we look at the density profiles for the case of half
porosity and a flat meniscus, $\phi=0.5,\eta=0$. Figure 7 shows this result.

Figure 7: Left: Far-field mass flux $J_{\infty}$, normalised with respect to
$J_{\infty}^{1\text{D}}$, against Mach number. Right: Far-field density
$\varrho_{\infty}$ as a ratio of saturation density, against Mach number.
Comparison between DSMC, and the linear and non-linear boundary conditions
models of the two 13-moment system, for the half-porosity case ($\phi=0.5$),
at Kn$=0.05$.

We can see clearly the improvement of the new model and thus the importance of
the boundary conditions to the system’s behaviour, as we may expect at Kn=0.05
where bulk non-equilibrium effects are less pronounced. We see the effect non-
linearity has on the relationship between Mach number and density, with
asymptotic behaviour similar to that discussed in the one-dimensional case
seen as Mach number increases. Moreover, the initial slope of the density plot
for $\mathrm{Ma}_{\infty}\approx 0$ is in very good agreement with that of
DMSC, meaning, along with better behaviour for larger Mach numbers, we see
even better agreement at lower Mach numbers than with the linear case.

Figure 8: Normalised far-field parameters against far-field Mach number for
the full-porosity case ($\phi=1$) and Kn$=0.05$. Shown are
$J_{\infty}/J_{\infty}^{1\text{D}}$ for (a) the flat meniscus, (b) the semi-
circular meniscus; and $\varrho_{\infty}/\varrho_{s}$,
$\Theta_{\infty}/\Theta_{L}$ for (c) the flat meniscus, (d) the semi-circular
meniscus. Comparison between DSMC, and non-linear boundary condition models
for NSF, Grad 13 and R13. Figure 9: Normalised far-field parameters against
far-field Mach number for the half-porosity case ($\phi=0.5$), at Kn$=0.05$.
Shown are $J_{\infty}/J_{\infty}^{1\text{D}}$ for (a) the flat meniscus, (b)
the semi-circular meniscus; and $\varrho_{\infty}/\varrho_{s}$,
$\Theta_{\infty}/\Theta_{L}$ for (c) the flat meniscus, (d) the semi-circular
meniscus. Comparison between DSMC, and non-linear boundary condition models
for NSF, Grad 13 and R13.

Further simulation data for combinations of $\phi=0.5,1$ and $\eta=0,1$, can
be seen in Figures 8 and 9. The improvement in accuracy allows us to look at
Mach numbers up to $\text{Ma}_{\infty}\approx 0.45$, and still retain
reasonable accuracy. In contrast, the fully linear case saw negative far-field
densities at $\text{Ma}_{\infty}\approx 0.25$. Figures 8 (a) and (c) show the
results from the one-dimensional case where analytic results were obtained.
Analytic results are not shown, but match well with the numerical model, which
was computed in a two-dimensional domain. The results give good agreement with
the DSMC results for all moment systems. The non-linearity of the solution
allows calculated values to closely follow the DSMC values for higher Mach
numbers, compared to the linear one-dimensional case, Figure 5 (c), where
density started with good agreement, but fell away from DSMC values as the
Mach number increased.

For the linear case, when walls were introduced, we saw far-field density
results diverge greatly from DSMC values at higher Mach numbers. The non-
linear boundary conditions give greatly improved results here over the linear
case, for both cases of curvature, seen in Figure 9 (compared to the linear
boundary condition case, Figure 6). We also see that for all results, when
$\mathrm{Ma}_{\infty}\approx 0$ the slopes of density are very accurate
compared to the DSMC plots, yielding excellent results for small Mach numbers.
As discussed in the linear case, low Mach number flows are an important region
of accuracy for the higher order moment systems.

The results presented give good results up to modest Mach numbers for all
cases of porosity and meniscus shape considered. This allows us to greatly
extend the range of Mach numbers in which we can have confidence in the moment
equations results compared to the linear case, where we could only be
confident up to $\mathrm{Ma}_{\infty}\approx 0.1$. For non-linear boundary
conditions, the 13-moment systems give accurate results up to
$\mathrm{Ma}_{\infty}\approx 0.3$. Therefore, at least for $\text{Kn}=0.05$,
the 13-moment systems make a good, efficient substitute for DSMC to reduce
computational expense on simulating nanoporous membrane dynamics.

##  VI Knudsen-Number Analysis

So far, we have focussed on $\text{Kn}=0.05$. The moment equations are
approximations of the Boltzmann equation with increasing accuracy in Knudsen
number, with NSF being first order accuracy, Grad 13 second order, and R13
third order [10, 31]. For the flat meniscus case, in the limiting case of
infinitely thin walls, $\phi\to 1$, the Knudsen number is irrelevant when
evaluating far-field parameters, and so we expect these models to be similar,
with differences coming from the different boundary condition models for each
set of equations. When either $\eta>0$ or $\phi<1$, the solutions become
dependent on the Knudsen number. This is seen in Figure 10, which gives the
results of the normalised far-field mass flux for increasing Knudsen number
for both linear and non-linear boundary condition models, for the
representative flat meniscus, half porosity case, and a far-field velocity of
$\mathcal{V}_{\infty}=30$m/s, corresponding222Assuming a molar mass of
$M=0.02896$kg/mol and temperature of $273$ K. to $v_{\infty}=0.1072$. All
macroscopic models retain accuracy up to around Kn=0.1. For Kn$>0.1$, the NSF
solutions diverge greatly from DSMC. Meanwhile, the Grad and R13 equations,
which are of higher-order accuracy in the Knudsen number, maintain reasonable
accuracy to the DSMC values. We do see the Grad 13 far-field mass flux
gradually decreasing, dropping from $J_{\infty}/J_{\infty}^{1D}=0.882$ at
$\text{Kn}=0.005$, to $J_{\infty}/J_{\infty}^{1D}=0.687$ at $\text{Kn}=1$ for
the non-linear case. On the other hand, R13 maintains a value of
$J_{\infty}/J_{\infty}^{1D}$ between 0.868 and 0.831 for all Knudsen numbers
for the non-linear boundary conditions case, hugging the DSMC values between
0.848 and 0.801. We also see that as $\text{Kn}\to 0$, all moment solutions
converge, as expected.

Figure 10: Far-field mass flux $J_{\infty}$ normalised with respect to
$J_{\infty}^{1D}$, against Knudsen number Kn. Comparison between DSMC, NSF,
Grad 13 and R13, for the half porosity ($\phi=0.5$), flat meniscus ($\eta=0$)
case, with $\mathcal{V}_{\infty}=30$m/s. Solid lines represent non-linear
boundary conditions models, dashed lines represent linear models.

All far-field parameters are assessed again with respect to Mach number for
various membrane configurations, this time for $\text{Kn}=1$, in the same way
as was done for $\text{Kn}=0.05$: mass flux $J_{\infty}$ is normalised with
respect to (26), and we assess the density and temperature ratios
$\varrho_{\infty}/\varrho_{s}$ and $\Theta_{\infty}/\Theta_{L}$. Simulations
are carried out using the non-linear boundary-condition model. Results are
shown in Figures 11 and 12.

Figure 11: Normalised far-field parameters against far-field Mach number for
the full-porosity case ($\phi=1$), at Kn$=1$. Shown are
$J_{\infty}/J_{\infty}^{1\text{D}}$ for (a) $\eta=0$, (b) $\eta=1$; and
$\varrho_{\infty}/\varrho_{s}$, $\Theta_{\infty}/\Theta_{L}$ are shown for (c)
$\eta=0$, (d) $\eta=1$. Comparison between DSMC, and non-linear boundary
condition models for NSF, Grad 13 and R13. Very large negative densities and
mass flux values are seen for NSF, and most of the computed results are not
shown. Figure 12: Normalised far-field parameters against far-field Mach
number for the half-porosity case ($\phi=0.5$), at Kn$=0.05$. From top to
bottom are $J_{\infty}/J_{\infty}^{1\text{D}}$, $\varrho_{\infty}/\varrho_{s}$
and $\Theta_{\infty}/\Theta_{L}$ for both cases of meniscus shape. Comparison
between DSMC, and non-linear boundary condition models for NSF, Grad 13 and
R13. Very large negative densities and mass flux values are seen for NSF, and
most of the computed results are not shown.

As expected, Figures 11 (a) and (c) show no change to the $\text{Kn}=0.05$
case, as Knudsen number is irrelevant for far-field parameter computation for
the case of a flat meniscus and no walls. For the semi-circular meniscus case
with no walls, shown in Figures 11 (b) and (d), the R13 system is still
showing good results at $\text{Kn}=1$, with all three parameters showing good
accuracy in comparison to DSMC. Grad 13 is also reasonably accurate in
comparison to DSMC, but less so than the R13. This is in-line with what we
would expect to see at higher Knudsen numbers, since R13 is of higher order
accuracy in Knudsen number than Grad 13. The NSF system, which is even lower
in accuracy in Knudsen number than Grad 13, is completely inaccurate in this
regime. The computed density and mass flux values are seen to deviate by large
factors from the DSMC results, even for small Mach numbers.

The case $\phi=0.5$ paints a similar picture; the higher Knudsen number has a
large impact on the free stream calculation of density and mass flux, with
very large negative values seen in Figures 12 (a), (b), (c) and (d) for the
NSF system. For free-stream velocities $v_{\infty}>0.4$, the FEM solver for
NSF failed to converge altogether. The reason for this is unclear, but a
failure to capture the flow near the singularity discussed in Appendix D
appears to be the most likely cause. The 13-moment systems again fair much
better at the higher Knudsen number than NSF, however Grad 13 sees negative
densities at $\text{Ma}_{\infty}\approx 0.25$ in Figure 12 (c), and even
sooner in Figure 12 (d). The R13 system once again is better than the other
moment systems, retaining good accuracy up to reasonable Mach numbers.

Changes in the porosity, $\phi$, seem to highlight the difference between the
Grad 13 and R13 systems: when $\phi=1$, far-field densities of the two
13-moment systems are comparable, as seen in Figure 11 (d); when $\phi=0.5$,
however, there are large differences in calculated densities between Grad 13
and R13, as seen in Figures 12 (c) and (d). This could be due to the fact that
Knudsen layers have a large effect when walls are introduced, as parameters
see large variations in values, particularly near the interface; see the flow-
field plots in Figure 3. An increased number of moments must be considered
when dealing with sharp gradients [65, 64, 69].

We can see some of the differences in the two 13-moment systems with the non-
linear boundary conditions by looking at the behaviour of temperature and heat
flux in the domain. The temperature surface plot, overlaid with heat flux
streamlines, is shown in Figure 13. With $\text{Kn}=1$, Figure 13 (a) shows
that the Grad 13 system has a large temperature drop adjecent to the wall,
while R13 exhibits larger temperature drops next to the meniscus, Figure 13
(b). The latter is in line with DSMC results, which show larger temperature
drops at the meniscus for all Knudsen numbers; see [28], page 14. Both moment
systems have flux lines emanating from the wall and absorbed at the meniscus.
Therefore, while streamlines of heat flux (in general) are directed from hot
to cold for the R13 system, we see an inverse Fourier-law prediction by the
Grad 13 system, with streamlines of heat flux directed from cold to hot.

Another difference in the two 13-moment systems is the presence of Knudsen
layers in the R13 system, which are absent in the Grad 13 system. The
existence of Knudsen layers is seen in Figure 14, where the temperature ratio
$\Theta/\Theta_{L}$ is plotted along the central vertical axis, $x=0$, for
$\text{Kn}=0.05$, $\eta=1$ and $\phi=0.5$. The Grad 13 profile for temperature
is mostly monotonically increasing, especially near the interface. Meanwhile,
both R13 and DSMC see a decrease in temperature near the interface, and then,
after reaching a minimum, we see the monotonically increasing behaviour, where
the calculated values plateau to equilibrium. The temperature minimum for R13
is attained at $y=-0.35$, corresponding to roughly three mean-free-path
lengths from the meniscus. Knudsen layers typically extend a few mean-free-
path lengths from the boundary, and we also saw that the one-dimensional case
had exponential decays on the order of the Knudsen number for R13, and so the
drop in temperature can be attributed to Knudsen boundary layer effects.

Figure 13: Temperature plots near the meniscus for (a) the Grad 13 systems,
and (b) the R13 system, overlaid with heat flux streamlines, for the flat
meniscus case, with $\phi=0.5$ and $\text{Kn}=1$.

Heat flux in the R13 system is also subject to Knudsen layer effects, similar
to those seen in the temperature just discussed. Heat flux profiles along the
lines $x=\pm\frac{1}{2}$ see an extremum at 0.5 mean-free-path lengths and 0.8
mean-free-path lengths from the interface for the tangential and normal
components of heat flux respectively. This results in a sign change inside the
domain, leading to the formation of heat flux vortices, seen in Figure 13 (b).
These vortices are not seen for Grad 13. Similar behaviour is discussed in the
problem of flow past a sphere, where Torrilhon shows the existence of heat
flux vortices for the R13 model, which are attributed to the Knudsen boundary
layers [19]. Moreover, Torrilhon showed that such vortices do not exist for a
hybrid Stokes-R13 flow, consisting of Stokes bulk equations with R13 boundary
conditions.

Figure 14: Temperature variation along the extent of of the domain, for the
semi-circular meniscus, at $\phi=0.5$ and $\text{Kn}=0.05$.

We have seen the R13 system able to produce accurate results up to Knudsen
numbers 1, while Grad 13 gives poor results in this regime, and NSF is
completely invalid. We saw evidence of Knudsen boundary layers from the R13
system, which are also in DSMC results, while Grad 13 saw no such behaviour,
which rationalises R13’s superiority, as rarefaction effects play a larger
role in the overall flow when the mean-free-path length of particles becomes
comparable to the flow’s characteristic length scale. We can have confidence
in the results from R13 with non-linear boundary conditions to be accurate for
all $0<\mathrm{Ma}_{\infty}<0.3$ and $0\leq\text{Kn}\leq 1$. This gives us a
large and useful range of membrane configurations which can be simulated
accurately and efficiently.

##  VII Conclusion

In this work, the linearised moment equations were used to describe
evaporation from a nanoporous membrane. Initially, a fully linear system was
considered, in particular with linear boundary conditions, however this failed
to give accurate results for Mach numbers above $\text{Ma}_{\infty}\approx
0.1$. When retaining non-linearity for velocity in the boundary conditions, we
derived an efficient method capable of reproducing DSMC simulation results,
effective up to reasonable Mach numbers and Knudsen numbers, which see huge
improvements over Navier-Stokes-Fourier results. The non-linear boundary
condition model for the R13 equations is able to produce good results for Mach
numbers up to 0.3 and up to $\text{Kn}=1$. The NSF is completely unviable in
this regime.

With the interfacial evaporation processes validated for the moment equations
versus DSMC simulations, future work can be devoted to extending the geometry
to a more representative version of a nanoporous membrane. This could include
the introduction of the liquid and solid elements of the membrane using
appropriate macroscopic equations, something not possible with DSMC
simulations due to the large variations in mean free path length of the
different materials. Introduction of new materials to the simulation requires
proper treatment of conservation laws at the interface, including mass,
momentum and energy balance, the basics of which can be found in [36] and
[37]. For a full description of the interface, particularly for curved
menisci, the surface tension must be introduced, relating the stresses to the
meniscus shape [70, 71].

In our model we assumed the existence of a far-field where parameters are
given space to decay. For a cooling device designed for micro-scale devices, a
more compact geometry is likely, in which a ceiling is brought in, and the gas
escapes from the sides due to a cross flow. Another direction this project
could be taken is to attempt to model this situation. The ceiling would
complicate the model since the evaporation would no longer be driven by the
far-field velocity, and may need to be driven by pressure gradients across the
interface. Introducing the full geometry with liquid and solid below the
interface discussed above could be introduced for the calculation of
saturation pressure, since the interplay between saturation pressure and
prescribed pressure gradients would play a large role in the evaporation from
the meniscus.

Another area of interest would be using the moment equations to simulate the
full three-dimensional geometry of a nanopore, something not done yet using
DSMC. The efficiency of the macroscopic approach would be very important for
this, since large parameter spaces in three dimensions is extremely
computationally expensive for DSMC. The two-dimensional results given here
provide justification and validation for the use of the 13-moment systems in
three-dimensions.

Future research may focus on analysis of the behaviour in pressure/density at
the wall-liquid-vapour contact line for the higher moment equations, as they
do not appear to exhibit singularities, while the NSF equations do (see Figure
18 in Appendix D). Singularity-free equations are desirable from a physical
viewpoint as well as for numerical methods, since mesh refinement can
drastically alter local parameter values when a singularity is present. This
study could involve a polar coordinate analysis similar to that performed by
Taylor [72], Moffat [73], and Nitsche and Parthasarthi [74].

More generally, our findings motivate the hybrid non-linear-boundary linear-
bulk equation system which reduces computational complexity from the fully-
non-linear approach and yet significantly extends the window of accuracy of
the moment equations to higher Mach numbers. Such a hybrid allows application
of the method of fundamental solutions, a highly-efficient numerical scheme
that employs fundamental solutions to the linear bulk equations (Lockerby and
Collyer [47], Claydon _et al._ [20]). This approach could be useful for a
variety of vapour/gas flows where moment equations are considered useful.

##  Appendix A Derivation of the Moment Equations

The equations are derived from basic considerations of particles in a gas in
the 6-dimensional phase space of positions and velocities representing the
system. Notation follows that of [10]. We consider the gas to be monatomic,
with each particle having mass $m$. The phase space consists of 3 positional
directions, $x_{i}$, and 3 velocity directions, $c_{i}$. The distribution
function $f(x_{i},c_{i},t)$ is defined such that
$f(x_{i},c_{i},t)\,d\mathbf{x}d\mathbf{c}$ gives the number of particles in
the cell $d\mathbf{x}d\mathbf{c}$ at time $t$. One considers the evolution of
the distribution function in a given subset $\Omega\subset\mathbb{R}^{3}$,
from which one obtains the Boltzmann equation,

$\frac{\partial f}{\partial t}+c_{k}\frac{\partial f}{\partial
x_{k}}+G_{k}\frac{\partial f}{\partial c_{k}}=\mathcal{S},$ (33)

where $G_{i}$ are external forces, and $\mathcal{S}=-\frac{\partial}{\partial
c_{k}}(W_{k}f)$ is the particle interaction term, $W_{k}$ denoting the
intermolecular forces.

Appropriate integration of $f$ over the velocity space gives the moments. For
example, multiplication of $f$ by $m$, the particle mass, and subsequent
integration gives the density

$\varrho=m\int f\,d\mathbf{c}.$

Multiplication of $f$ by $c_{i}$ and subsequent integration gives the momentum

$\varrho\mathcal{V}_{i}=m\int c_{i}f\,d\mathbf{c},$

The moments are generated by multiplication of $f$ by polynomials in $c_{i}$.
Following Struchtrup [10], we use the vector of polynomials

$\Psi^{[13]}=m\left\\{1,c_{i},\frac{1}{3}C^{2},C_{\langle
i}C_{j\rangle},\frac{1}{2}C^{2}C_{i}\right\\}$

to generate our 13 moments. Indices in angular brackets denote the trace-free-
symmetric part of a tensor [10], and $C_{i}=c_{i}-\mathcal{V}_{i}$ is the
peculiar velocity. This gives the 13 moments333The stress tensor is trace-free
and symmetric and so has 5 independent components. as

$\Phi^{[13]}=\varrho\left\\{1,\mathcal{V}_{i},\Theta,\Sigma_{ij},\mathcal{Q}_{i}\right\\}.$

In the above, $\Theta$ is the temperature in specific energy units, i.e.
$\Theta=\frac{k}{m}T$, where $T$ is the temperature, and $k$ is the Boltzmann
constant. Additionally, $\Sigma_{ij}$ is the stress tensor, and
$\mathcal{Q}_{i}$ is the heat flux tensor. We also study a five-moment system,
the NSF system, consisting of the moments
$\Phi^{[5]}=\varrho\left\\{1,\mathcal{V}_{i},\Theta\right\\}$.

Next, we derive the equations describing the moments $\Phi^{[13]}$. We perform
a similar process to the one carried out to obtain the moments $\Phi^{[13]}$,
but carried out on the Boltzmann equation (33), rather than the distribution
function $f$. That is, we multiply (33) by a polynomial in $c_{i}$, and then
integrate over the velocity space. For instance, multiplication of (33) by $m$
and subsequent integrating gives the conservation of mass equation

$D_{t}\varrho+\varrho\partial_{k}\mathcal{V}_{k}=0.$

Performing this process for all polynomials in the vector $\Psi^{[13]}$ gives
the equations (1)-(3).

##  Appendix B Domain finite element mesh

The domain mesh for a representative case of $\text{Kn}=0.05$ and with half
porosity ($\phi=0.5$) is shown in Figure 15, for both the flat meniscus and
semicircular meniscus case. The mesh consists of 12369 triangular domain
elements and 580 boundary elements for the flat meniscus case, and 14115
triangular domain elements and 637 boundary elements for the semicircular
meniscus case. The mesh is made finest near the interface where most two-
dimensional-type flow occurs. Near the far-field, where the flow is
essentially one-dimensional, only a coarse mesh is required

Figure 15: Domain mesh for nanoporous evaporation.

##  Appendix C Fourier decomposition for the two-dimensional case

Here we examine the structure of the analytic solutions for the flat meniscus
with arbitrary porosity. For porosity $\phi=1$, we saw a linear relationship
between density and far-field velocity. We will see that this is the case also
for arbitrary $\phi$. We assume parameters take the form

$\displaystyle\psi(x,y)=\sum_{k\in\mathbb{Z}}\psi_{k}(y)e^{-ik\pi\frac{x}{\ell}},$

where $\ell=\frac{1}{2}L+W$ is the (half) length of periodicity (the half
coming from the factor of 2 in front of $\pi$).

With this, the linearised moment equations (13)-(14) become an infinite set of
ordinary differential systems, one for each mode $k$. For each mode, the
(arbitrary $N$) moment system can be reduced to a matrix differential system
of the form

$\frac{d}{dy}\varphi^{[N]}_{k}=Q_{k}\varphi^{[N]}_{k},\qquad\text{(no sum over
}k)$ (34)

for the vector of $N$ moments

$\varphi^{[N]}_{k}=\left\\{\rho_{k}(y),v_{i,k}(y),\theta_{k}(y),\sigma_{ij,k}(y),q_{i,k}(y),...\right\\}^{t}.$

Here, $Q_{k}$ is a matrix dependent on $k$, $\ell$ and Kn.

The general solution to (34) has the form

$\displaystyle\varphi^{[N]}_{k}=\exp(Q_{k}y)\mathbf{C}^{\prime}_{k}=P_{k}\exp(J_{k}y)\mathbf{C}_{k}$

where $J_{k}$ is the Jordan matrix of $Q_{k}$, $P_{k}$ is the matrix of
generalised eigenvalues, and
$\mathbf{C}_{k}=P_{k}^{-1}\mathbf{C}^{\prime}_{k}$ is a column vector of
constants dependent on $k$.

The boundary condition (17) becomes

$v_{2}(x,0)=\sqrt{\frac{2}{\pi}}\left(p_{s}(x,0)-\Pi(x,0)-\frac{1}{2}\theta_{L}(x,0)\right)\chi_{\ell}(x),$
(35)

where

$\displaystyle\chi_{\ell}(x)=\begin{cases}1,&|x|<\frac{1}{2}\\\
0,&\frac{1}{2}\leq|x|\leq\frac{1}{2}+W,\end{cases}$

and other moment boundary conditions remain the same, as well as the moment
boundary conditions for the general $N$ moment system. The effective pressure
$\Pi$ is also the generalisation for the $N$ moment system.

The remaining boundary conditions do not involve the function $\chi_{\ell}$,
and so can be written as an infinite list of $N-1\times N-1$ matrix boundary
condition systems of the form

$\displaystyle\mathbf{b}_{k}=A_{k}\mathbf{C}_{k},$

one for each $k$, similar to ($*$ ‣ IV). These can be solved directly to write
all but one of the constants in $\mathbf{C}_{k}$ linearly in terms of the
remaining one, leaving a single constant to be solved for, say $z_{k}$, for
each $k$. The $z_{k}$ are found using (35).

The function $\chi_{\ell}$ can be written as a Fourier expansion, with

$\displaystyle\chi_{\ell}=\sum_{j\in\mathbb{Z}}a_{j}e^{-i\pi jx},$

where

$\displaystyle a_{j}=\frac{1}{2}\int_{-\frac{1}{2}}^{\frac{1}{2}}e^{i\pi
jx}dx=\begin{cases}\frac{1}{\pi j}\sin\left(\frac{\pi j}{2}\right)&j\neq 0\\\
\frac{1}{2}&j=0\end{cases}$

The $k$-th mode of (35) becomes

$\displaystyle
v_{2,k}(0)+\sqrt{\frac{2}{\pi}}\sum_{r\in\mathbb{Z}}\Pi_{r}(0)a_{k-r}=0$ (36)

The parameters $v_{2,k}$ and $\Pi_{k}$ are linear in $z_{k}$, thus we can
write (36) as

$\displaystyle
f(k,\ell)z_{k}+\sqrt{\frac{2}{\pi}}\sum_{r\in\mathbb{Z}}g(r,\ell)a_{k-r}z_{r}=0,$
(37)

for all $k$, for some known functions $f$ and $g$ linear in $k$.

Looking more closely at the differential equations for each mode, the zeroth
mode corresponds to the one-dimensional system, so we deduce that
$v_{2,0}=v_{\infty}$, since the continuity equation dictates the velocity is
constant in $y$ for this mode. This means that any of the moments
$\varphi^{[N]}_{k}$ will have at most linear dependence on $v_{\infty}$ and
$z_{n}$, with no products. Therefore, the zero-th moment system will take the
form

$\displaystyle\sqrt{\frac{2}{\pi}}\sum_{r\in\mathbb{Z}}g(r,\ell)a_{-r}z_{r}=-v_{\infty}.$
(38)

The entire system (37) can be brought into the form

$\displaystyle F_{rk}z_{r}=-v_{\infty}\delta_{0k},$ (39)

where $\delta_{ij}$ is the Kronecker delta symbol, and summation is performed
over $r$. The right hand side has linear dependence on $v_{\infty}$, and since
$F_{ij}$ has no dependence on $v_{\infty}$, we must have that $z_{k}$ is at
most linear in $v_{\infty}$ for every $k$.

To verify this solution, this process is done for isothermal NSF system (with
no slip or temperature jump). The matrix $Q_{k}$ in this case is

$\displaystyle Q_{k}=\left(\begin{matrix}0&i\frac{\pi
k}{\ell}&0&\frac{-1}{\text{Kn}}\\\ i\frac{\pi k}{\ell}&0&0&0\\\
0&-2\text{Kn}\left(\frac{\pi k}{\ell}\right)^{2}&0&-i\frac{\pi k}{\ell}\\\
-2\text{Kn}\left(\frac{\pi k}{\ell}\right)^{2}&0&i\pi k&0\end{matrix}\right),$
(40)

and then for $k\neq 0$ we have

$\displaystyle J_{k}=\left(\begin{array}[]{cccc}-\frac{\pi k}{\ell}&1&0&0\\\
0&-\frac{\pi k}{\ell}&0&0\\\ 0&0&\frac{\pi k}{\ell}&1\\\ 0&0&0&\frac{\pi
k}{\ell}\\\ \end{array}\right),$ $\displaystyle
P_{k}=\left(\begin{array}[]{cccc}\frac{\ell}{2\text{Kn}\pi
k}&0&-\frac{\ell}{2\text{Kn}\pi k}&0\\\ -\frac{i\ell}{2\text{Kn}\pi
k}&-\frac{i\ell^{2}}{2\text{Kn}\pi^{2}k^{2}}&-\frac{i\ell}{2\text{Kn}\pi
k}&\frac{i\ell^{2}}{2\text{Kn}\pi^{2}k^{2}}\\\ 0&-\frac{i\ell}{\pi
k}&0&-\frac{i\ell}{\pi k}\\\ 1&0&1&0\\\ \end{array}\right),$

and

$\displaystyle J_{0}=\left(\begin{array}[]{cccc}0&0&0&0\\\ 0&0&0&0\\\
0&0&0&1\\\ 0&0&0&0\\\ \end{array}\right),$ $\displaystyle
P_{0}=\left(\begin{array}[]{cccc}0&0&1&0\\\ 0&1&0&0\\\ 1&0&0&0\\\
0&0&0&-\text{Kn}\\\ \end{array}\right).$

The $F_{jk}$ are given by

$\displaystyle
F_{jk}=\begin{cases}\sqrt{\frac{2}{\pi}}\frac{\ell}{\pi(j-k)}\sin\left(\frac{\pi(j-k)}{2\ell}\right)&j\neq
k,\\\ \frac{1}{2}\sqrt{\frac{2}{\pi}}+\frac{\ell}{2\text{Kn}\pi|j|}&j=k\neq
0,\\\ \frac{1}{2}\sqrt{\frac{2}{\pi}}&j=k=0.\end{cases}$

The system (39) is solved by truncating the system, considering only $|i|\leq
M$ for some $M\in\mathbb{N}$. Normal velocity and density profiles across the
interface are shown in Figure 16, and are seen to converge to the numerical
solution as $M$ increases. We also see evidence of the singularity in density
present in the NSF equations along the wall-liquid-vapour contact line, which
we discuss in the next appendix.

Figure 16: Velocity and density profiles along the interface for isothermal
NSF, with $v_{\infty}=0.1$, $\phi=0.5$ (i.e. $\ell=2$), $\text{Kn}=0.05$.
Comparisons between the FEM approach and the Fourier approach, showing
convergence in $M$, for three values of $M$; $M=10,100$ and $1000$.

##  Appendix D A Brief Discussion on the Singularities at the Wall-Liquid-
Vapour Contact Line

As mentioned in the main text, the density profile seen in Figure 4 for NSF is
caused by a singularity in the interfacial density at the Wall-Liquid-Vapour
(WLV) contact point. This singularity can also be seen in Figure 16 in
Appendix C, where the density drops dramatically on the wall side of the
interface. The behaviour of pressure (and thus density) at such contact lines
for Stokes flow is discussed by Nitsche and Parthasarathi [74]. For the Hertz-
Knudsen-Schrage-like boundary condition

$\displaystyle v_{n}=-p,$

they found a singularity of order $r^{-\frac{1}{2}}$ in pressure, compared to
the usual $r^{-1}$ type singularity found in Stokes flow under the
prescription of uniform evaporative normal velocity. Validation of this
singularity class is seen in Figure 17 (recall that we assume an ideal gas,
and in this instance temperature is constant, so density is proportional to
pressure), where numerical results from our model are seen to be approximately
proportional to $r^{-\frac{1}{2}}$ as $r\to 0$ (there will be other exponents
of $r$ contributing to pressure, though $-\frac{1}{2}$ is the dominant
exponent for small $r$). While $r^{-\frac{1}{2}}$ is an integrable
singularity, the singularity still manifests in the numerical model for NSF,
and causes significant mesh-size-dependence for the pressure profile at the
WLV contact point, as well as large spikes in values. For the singularity-free
higher moment methods, finer meshes produce graphically indistinguishable
curves in all figures. For instance, the relationship with mesh size is
examined in Figure 18 for Grad 13 and NSF. NSF is clearly dependent on the
mesh size, whereas Grad 13 converges on a value as the mesh becomes finer. For
context, the mesh size at the WLV contact point for all simulations carried
out in the main sections of this paper was $5\times 10^{-3}$.

Figure 17: Numerically calculated density for the NSF model for increasing
radius from the WLV contact point, suggesting a singularity of type $r^{-1/2}$
Figure 18: Numerical values for density at the WLV contact point, for
increasing mesh size, and for the NSF model and Grad 13 model. Curves are
normalised with respect to the highest attained (absolute) value of density,
with $N_{\text{NSF}}=22.73$, and $N_{\text{G13}}=0.3093$

##  Appendix E Analytic Solutions for the Non-linear Boundary Condition Model
in the 1D Case

Inversion of ($**$ ‣ V) yields

$\displaystyle\rho_{\infty}^{[\text{G}13]}=-\frac{v_{\infty}\left(10\pi
v_{\infty}+7\sqrt{2\pi}\right)}{10\pi
v_{\infty}^{2}+17\sqrt{2\pi}v_{\infty}+16},\quad\theta_{\infty}^{[\text{G}13]}=-\frac{2\sqrt{2\pi}v_{\infty}}{10\pi
v_{\infty}^{2}+17\sqrt{2\pi}v_{\infty}+16},$

and

$\displaystyle\rho_{\infty}^{[\text{R}13]}$
$\displaystyle=-\frac{\sqrt{\pi}v_{\infty}\left[\left(5\pi\sqrt{30}+36\sqrt{2\pi}\right)v_{\infty}+7\sqrt{15\pi}+57\right]}{D},$
$\displaystyle\theta_{\infty}^{[\text{R}13]}$
$\displaystyle=-\frac{2\sqrt{\pi}v_{\infty}\left(2\sqrt{2\pi}v_{\infty}+\sqrt{15\pi}+9\right)}{D},$
$\displaystyle\mathbf{C}_{8}^{[\mathrm{R13}]}$
$\displaystyle=\frac{10\sqrt{\pi}v_{\infty}\left(\sqrt{2\pi}v_{\infty}+2\right)}{D},$

where

$\displaystyle D=$ $\displaystyle
5\pi\sqrt{2}\left(\sqrt{15\pi}+6\right)v_{\infty}^{2}+\left(17\pi\sqrt{15}+115\sqrt{\pi}\right)v_{\infty}+60\sqrt{2}+8\sqrt{30\pi}.$

Finally, the density for the NSF system is

$\displaystyle\rho_{\infty}^{[\mathrm{NSF,1D}]}=-\frac{v_{\infty}}{\sqrt{\frac{2}{\pi}}+v_{\infty}}.$

All expressions are seen to converge as $v_{\infty}\to\infty$.

###### Acknowledgements.

This work was supported by the EPSRC (Grants No. EP/N016602/1, No.
EP/P020887/1, and No. EP/P031684/1) and an EPSRC IAA Acceleration Account
(EP/R511808/1).

## References

  * [1] K. R. Bagnall, Device-level thermal analysis of GaN-based electronics.  PhD thesis, Massachusetts Institute of Technology, 2013.
  * [2] Z. Lu, T. R. Salamon, S. Narayanan, K. R. Bagnall, D. F. Hanks, D. S. Antao, B. Barabadi, J. Sircar, M. E. Simon, and E. N. Wang, “Design and modeling of membrane-based evaporative cooling devices for thermal management of high heat fluxes,” IEEE transactions on components, packaging and manufacturing technology, vol. 6, 2016.
  * [3] D. F. Hanks, Z. Lu, J. Sircar, T. R. Salamon, D. S. Antao, K. R. Bagnall, B. Barabadi, and E. N. Wang, “Nanoporous membrane device for ultra high heat flux thermal management,” Microsystems & Nanoengineering, vol. 4, Dec. 2018.
  * [4] Z. Lu, Design and modeling of a high flux cooling device based on thin film evaporation from thin nanoporous membranes.  PhD thesis, Massachusetts Institute of Technology, 2014.
  * [5] C. G. Gray, K. E. Gubbins, and C. G. Joslin, Theory of Molecular Fluids: Volume 2: Applications, vol. 10.  Oxford University Press, 2011.
  * [6] Y. Akkuş, “Modeling of evaporation from nanoporous membranes using molecular dynamics simulation,” Isı Bilimi ve Tekniği Dergisi, vol. 39, no. 1, pp. 91–99, 2019.
  * [7] G. A. Bird, Molecular gas dynamics and the direct simulation of gas flows.  1994\.
  * [8] G. Dimarco, R. Loubère, J. Narski, and T. Rey, “An efficient numerical method for solving the Boltzmann equation in multidimensions,” Journal of Computational Physics, vol. 353, pp. 46–81, 2018.
  * [9] W. Su, L. Zhu, and L. Wu, “Fast convergence and asymptotic preserving of the general synthetic iterative scheme,” SIAM Journal on Scientific Computing, vol. 42, no. 6, pp. B1517–B1540, 2020.
  * [10] H. Struchtrup, “Macroscopic transport equations for rarefied gas flows,” in Macroscopic Transport Equations for Rarefied Gas Flows, Springer, 2005.
  * [11] H. Struchtrup and M. Torrilhon, “Higher-order effects in rarefied channel flows.,” Physical review. E, Statistical, nonlinear, and soft matter physics, vol. 78 4 Pt 2, p. 046301, 2008.
  * [12] X. Gu and D. Emerson, “A computational strategy for the regularized 13 moment equations with enhanced wall-boundary conditions,” Journal of computational physics, vol. 225, no. 1, pp. 263–283, 2007.
  * [13] H. Grad, “On the kinetic theory of rarefied gases,” Communications on pure and applied mathematics, vol. 2, no. 4, pp. 331–407, 1949.
  * [14] H. Struchtrup and M. Torrilhon, “Regularization of Grad’s 13 moment equations: Derivation and linear analysis,” Physics of Fluids, vol. 15, no. 9, pp. 2668–2680, 2003.
  * [15] Y. Fan, J. Koellermeier, J. Li, R. Li, and M. Torrilhon, “Model reduction of kinetic equations by operator projection,” Journal of Statistical Physics, vol. 162, no. 2, pp. 457–486, 2016.
  * [16] J. Koellermeier, R. P. Schaerer, and M. Torrilhon, “A framework for hyperbolic approximation of kinetic equations using quadrature-based projection methods,” Kinetic & Related Models, vol. 7, no. 3, p. 531, 2014.
  * [17] P. Taheri, A. S. Rana, M. Torrilhon, and H. Struchtrup, “Macroscopic description of steady and unsteady rarefaction effects in boundary value problems of gas dynamics,” Continuum Mechanics and Thermodynamics, vol. 21, no. 6, pp. 423–443, 2009.
  * [18] P. Taheri and H. Struchtrup, “Effects of rarefaction in microflows between coaxial cylinders,” Physical Review E, vol. 80, no. 6, p. 066317, 2009\.
  * [19] M. Torrilhon, “Slow gas microflow past a sphere: analytical solution based on moment equations,” Physics of Fluids, vol. 22, no. 7, p. 072001, 2010.
  * [20] R. Claydon, A. Shrestha, A. S. Rana, J. E. Sprittles, and D. A. Lockerby, “Fundamental solutions to the regularised 13-moment equations: efficient computation of three-dimensional kinetic effects,” Journal of Fluid Mechanics, vol. 833, 2017.
  * [21] A. Rana, M. Torrilhon, and H. Struchtrup, “A robust numerical method for the r13 equations of rarefied gas dynamics: Application to lid driven cavity,” Journal of Computational Physics, vol. 236, pp. 169–186, 2013.
  * [22] M. Torrilhon, “Two-dimensional bulk microflow simulations based on regularized grad’s 13-moment equations,” Multiscale Modeling & Simulation, vol. 5, no. 3, pp. 695–728, 2006.
  * [23] L. Theisen and M. Torrilhon, “fenicsr13: A tensorial mixed finite element solver for the linear R13 equations using the fenics computing platform,” ACM Transactions on Mathematical Software (TOMS), vol. 47, no. 2, pp. 1–29, 2021.
  * [24] A. Westerkamp and M. Torrilhon, “Finite element methods for the linear regularized 13-moment equations describing slow rarefied gas flows,” Journal of Computational Physics, vol. 389, pp. 1–21, 2019.
  * [25] A. Westerkamp, A Continuous Interior Penalty Method for the Linear Regularized 13-Moment Equations Describing Rarefied Gas Flows.  PhD thesis, Universitätsbibliothek der RWTH Aachen, 2017.
  * [26] A. S. Rana, S. Saini, S. Chakraborty, D. A. Lockerby, and J. E. Sprittles, “Efficient simulation of non-classical liquid–vapour phase-transition flows: a method of fundamental solutions,” Journal of Fluid Mechanics, vol. 919, 2021.
  * [27] M. Torrilhon, “Modeling nonequilibrium gas flow based on moment equations,” Annual review of fluid mechanics, vol. 48, pp. 429–458, 2016.
  * [28] B. John, R. Enright, J. E. Sprittles, L. Gibelli, D. R. Emerson, and D. A. Lockerby, “Numerical investigation of nanoporous evaporation using direct simulation Monte Carlo,” Physical Review Fluids, vol. 4, Nov. 2019\.
  * [29] B. John, L. Gibelli, R. Enright, J. E. Sprittles, D. A. Lockerby, and D. R. Emerson, “Evaporation from arbitrary nanoporous membrane configurations: An effective evaporation coefficient approach,” Physics of Fluids, vol. 33, no. 3, p. 032022, 2021.
  * [30] H. Grad, “Principles of the kinetic theory of gases,” in Thermodynamik der Gase/Thermodynamics of Gases, pp. 205–294, Springer, 1958.
  * [31] H. Struchtrup, “Stable transport equations for rarefied gases at high orders in the Knudsen number,” Physics of Fluids, vol. 16, no. 11, pp. 3921–3934, 2004.
  * [32] X.-j. Gu and D. R. Emerson, “A high-order moment approach for capturing non-equilibrium phenomena in the transition regime,” Journal of fluid mechanics, vol. 636, pp. 177–216, 2009.
  * [33] M. Torrilhon and H. Struchtrup, “Regularized 13-moment equations: shock structure calculations and comparison to burnett models,” Journal of Fluid Mechanics, vol. 513, pp. 171–198, 2004.
  * [34] M. Y. Timokhin, H. Struchtrup, A. Kokhanchik, and Y. A. Bondar, “Different variants of r13 moment equations applied to the shock-wave structure,” Physics of Fluids, vol. 29, no. 3, p. 037105, 2017.
  * [35] A. Bobylev, “Instabilities in the chapman-enskog expansion and hyperbolic burnett equations,” Journal of statistical physics, vol. 124, no. 2, pp. 371–399, 2006.
  * [36] H. Struchtrup, A. Beckmann, A. Rana, and A. Frezzotti, “Evaporation boundary conditions for the R13 equations of rarefied gas dynamics,” Physics of Fluids, Sept. 2017.
  * [37] M. Bond and H. Struchtrup, “Mean evaporation and condensation coefficients based on energy dependent condensation probability,” Physical Review E, vol. 70, no. 6, p. 061605, 2004.
  * [38] C. Navier, “Mémoire sur les lois du mouvement des fluides,” Mémoires de l’Académie Royale des Sciences de l’Institut de France, vol. 6, no. 1823, pp. 389–440, 1823.
  * [39] K. Wang, Z. Chai, G. Hou, W. Chen, and S. Xu, “Slip boundary condition for lattice Boltzmann modeling of liquid flows,” Computers & Fluids, vol. 161, pp. 60–73, 2018.
  * [40] M. Torrilhon and H. Struchtrup, “Boundary conditions for regularized 13-moment-equations for micro-channel-flows,” Journal of Computational Physics, vol. 227, no. 3, pp. 1982–2011, 2008.
  * [41] M. Gad-el Hak, “The fluid mechanics of microdevices—the freeman scholar lecture,” 1999.
  * [42] A. Rana, D. Lockerby, and J. Sprittles, “Lifetime of a nanodroplet: kinetic effects and regime transitions,” Physical review letters, vol. 123, no. 15, p. 154501, 2019.
  * [43] X.-J. Gu, D. R. Emerson, and G.-H. Tang, “Kramers’ problem and the knudsen minimum: a theoretical analysis using a linearized 26-moment approach,” Continuum Mechanics and Thermodynamics, vol. 21, no. 5, pp. 345–360, 2009.
  * [44] Z. Cai and M. Torrilhon, “Numerical simulation of microflows using moment methods with linearized collision operator,” Journal of Scientific Computing, vol. 74, no. 1, pp. 336–374, 2018.
  * [45] J. Koellermeier and M. Torrilhon, “Numerical study of partially conservative moment equations in kinetic theory,” Communications in Computational Physics, vol. 21, no. 4, pp. 981–1011, 2017.
  * [46] J. Koellermeier and M. Torrilhon, “Two-dimensional simulation of rarefied gas flows using quadrature-based moment equations,” Multiscale Modeling & Simulation, vol. 16, no. 2, pp. 1059–1084, 2018.
  * [47] D. A. Lockerby and B. Collyer, “Fundamental solutions to moment equations for the simulation of microscale gas flows,” Journal of Fluid Mechanics, vol. 806, pp. 413–436, 2016.
  * [48] T. Ytrehus, “Molecular-flow effects in evaporation and condensation at interfaces,” Multiphase Science and Technology, vol. 9, no. 3, 1997.
  * [49] COMSOL Multiphysics, “Comsol multiphysics v5.6,” 2020.
  * [50] P. M. Gresho and R. L. Sani, “Incompressible flow and the finite element method. volume 1: Advection-diffusion and isothermal laminar flow,” 1998.
  * [51] A. S. Rana and H. Struchtrup, “Thermodynamically admissible boundary conditions for the regularized 13 moment equations,” Physics of Fluids, vol. 28, no. 2, p. 027105, 2016.
  * [52] N. Sarna and M. Torrilhon, “On stable wall boundary conditions for the hermite discretization of the linearised boltzmann equation,” Journal of Statistical Physics, vol. 170, no. 1, pp. 101–126, 2018.
  * [53] S. J. Plimpton, and M. A. Gallis, “Sparta direct simulation monte carlo (DSMC) simulator,” 2015.
  * [54] J. Fan, “Preface: 29th international symposium on rarefied gas dynamics,” in AIP Conference Proceedings, vol. 1628, pp. 1–2, American Institute of Physics, 2014.
  * [55] M. A. Gallis, T. Koehler, J. R. Torczynski, and S. J. Plimpton, “Direct simulation monte carlo investigation of the rayleigh-taylor instability,” Physical Review Fluids, vol. 1, no. 4, p. 043403, 2016.
  * [56] K. Nanbu, “Direct simulation scheme derived from the Boltzmann equation. ii. multicomponent gas mixtures,” Journal of the Physical Society of Japan, vol. 49, no. 5, pp. 2050–2054, 1980.
  * [57] H. Babovsky and R. Illner, “A convergence proof for nanbu’s simulation method for the full Boltzmann equation,” SIAM journal on numerical analysis, vol. 26, no. 1, pp. 45–65, 1989.
  * [58] W. Wagner, “A convergence proof for Bird’s direct simulation monte carlo method for the Boltzmann equation,” Journal of Statistical Physics, vol. 66, no. 3, pp. 1011–1044, 1992.
  * [59] D. Labuntsov and A. Kryukov, “Analysis of intensive evaporation and condensation,” International Journal of Heat and Mass Transfer, vol. 22, no. 7, pp. 989–1002, 1979.
  * [60] C. Cercignani, Rarefied gas dynamics: from basic concepts to actual calculations, vol. 21.  Cambridge University Press, 2000.
  * [61] N. Hadjiconstantinou and A. Garcia, “Statistical error in particle simulations of low mach number flows,” Elsevier, December 2000.
  * [62] F. Sharipov, “Modeling and calculations of rarefied gas flows: DSMC vs kinetic equation,” Invited talk, vol. 45, 2011.
  * [63] D. W. Mackowski, D. H. Papadopoulos, and D. E. Rosner, “Comparison of Burnett and DSMC predictions of pressure distributions and normal stress in one-dimensional, strongly nonisothermal gases,” Physics of Fluids, vol. 11, no. 8, pp. 2108–2116, 1999.
  * [64] I. Müller and T. Ruggeri, Rational extended thermodynamics, vol. 37.  Springer Science & Business Media, 2013.
  * [65] H. Struchtrup, “Heat transfer in the transition regime: Solution of boundary value problems for Grad’s moment equations via kinetic schemes,” Physical Review E, vol. 65, no. 4, p. 041204, 2002.
  * [66] J. Koellermeier, “Error estimators for adaptive simulation of rarefied gases using hyperbolic moment models,” in AIP Conference Proceedings, vol. 2132, p. 120004, AIP Publishing LLC, 2019.
  * [67] M. Abdelmalik and E. van Brummelen, “Error estimation and adaptive moment hierarchies for goal-oriented approximations of the boltzmann equation,” Computer Methods in Applied Mechanics and Engineering, vol. 325, pp. 219–239, 2017.
  * [68] M. Torrilhon, J. Au, and H. Struchtrup, “Explicit fluxes and productions for large systems of the moment method based on extended thermodynamics,” Continuum Mechanics and Thermodynamics, vol. 15, no. 1, pp. 97–111, 2003.
  * [69] M. Torrilhon and N. Sarna, “Hierarchical boltzmann simulations and model error estimation,” Journal of Computational Physics, vol. 342, pp. 66–84, 2017\.
  * [70] J. Young, “The condensation and evaporation of liquid droplets in a pure vapour at arbitrary Knudsen number,” International journal of heat and mass transfer, vol. 34, no. 7, pp. 1649–1661, 1991.
  * [71] I. Müller, Thermodynamics.  Pitman, Boston, 1985.
  * [72] G. Taylor, “On scraping viscous fluid from a plane surface,” The Scientific Papers of Sir Geoffrey Ingram Taylor, pp. 410–413, 1962.
  * [73] H. K. Moffatt, “Viscous and resistive eddies near a sharp corner,” Journal of Fluid Mechanics, vol. 18, no. 1, pp. 1–18, 1964.
  * [74] L. C. Nitsche and P. Parthasarathi, “Stokes flow singularity at the junction between impermeable and porous walls,” Journal of Fluid Mechanics, vol. 713, pp. 183–215, 2012.