Efficient simulation of non-classical liquid-vapour phase-transition flows: a method of fundamental solutions

When a liquid evaporates it loses energy; that is why sweating cools us down. Evaporation is an effective cooling mechanism widely employed by nature; for example, as a body-temperature control for mammals (Sherwood, 2005), and in engineering applications, such as spray drying, spray cooling and microelectronics cooling (Dhavaleswarapu et al., 2012; Plawsky et al., 2014; Ling et al., 2014). The use of pot-in-pot coolers by mankind can be tracked back to Bronze Age civilisation (Evans, 2000); a device that relies on the principle of evaporative cooling.

The description of phase-transition processes in micro/nano scales is intriguing, mainly due to the inability of classical continuum theories to accurately capture the rarefied vapour-flow characteristic. Specifically, as the representative physical length scale ($\ell$) of the flow becomes comparable to the mean free path ($\lambda$) in the vapour, i.e., the Knudsen number $(\mathrm{Kn}=\lambda/\ell)\approx 1$ (Bird, 1994; Cercignani, 2000; Sone, 2007; Struchtrup, 2005; Torrilhon, 2016). The classical continuum description based on the Euler or Navier–Stokes–Fourier (NSF) equations is valid only in regimes for which $\mathrm{Kn}\lesssim 0.001$—referred to as the hydrodynamic regime. There are, of course, computational atomistic techniques, such as molecular dynamics simulations (MD), which are capable of modelling phase transition at nano-scales. However, due to computational time and memory limitations, it is impractical to use them for multi-scale processes, which span over a wide range of time and length scales.

The behaviour of a dilute vapour phase can readily be described by the Boltzmann equation, which is an evolution equation for the (velocity) distribution function of vapour molecules (Bird, 1994; Cercignani, 2000). The Boltzmann equation offers the complete mesoscopic description of the vapour for all range of $\mathrm{Kn}$; however, the exact analytical solution of the Boltzmann equation is intractable in a general case while its numerical solutions are computationally very expensive. Grad—in his seminal work (Grad, 1949a, b)—proposed an asymptotic solution procedure for the Boltzmann equation through the method of moment. A variety of other extended fluid dynamics equations, aiming to describe processes under rarefied conditions, have been developed, see e.g. Eu (1992); Müller & Ruggeri (1998); Kremer (2010); Struchtrup (2005); Chakraborty & Durst (2007); Dongari et al. (2010); Myong (2014); Torrilhon (2016). These models close the gap between classical fluid dynamics and kinetic theory; that is, they aim at a good description in the transition regime ($0.001\lesssim\mathrm{Kn}\lesssim 1$).

The classical phase transition models, found in the literature, assume that the temperature and the velocity tangential to the interface are continuous across it. However, it is now well established, both experimentally and theoretically (McGaughey & Ward, 2002; Bond & Struchtrup, 2004; Rana et al., 2019), that such assumptions are invalid at a nano-confined liquid-vapour interface. For such small length scales, the Knudsen number lies in the transition regime (or beyond) and a difference in velocity and temperature jump are observed across the interface—thus it is important to employ models capable of describing strong non-equilibrium effects that occur at high Knudsen numbers. In Rana et al. (2019), the lifetime of an evaporating nanodroplet was studied. The study revealed that, when the drop size is large ($>$micron), its diameter-squared decreases in proportion to time; known as the $d^{2}$-law. However, as the diameter approaches sub micron and nano-scales, a crossover to a new behavior is observed, with the diameter now reducing in proportion to time (following a ‘$d$–law’). By taking into account the temperature discontinuity at the liquid-vapour interface,the drop’s lifetime was correctly predicted.

Notably, a full theory on boundary conditions is still missing for these nonlinear extended fluid dynamics equations, for which the notion of solvability and stability is more intricate. The approach suggested by Grad (1958), and recently explored by researchers (Gu & Emerson, 2007; Torrilhon & Struchtrup, 2008; Gu & Emerson, 2009; Struchtrup et al., 2017; Rana et al., 2018b), violates Onsager symmetry conditions (Rana & Struchtrup, 2016; Beckmann et al., 2018; Sarna & Torrilhon, 2018), leading to instabilities which influence the convergence of numerical schemes, and prevent convergence of moments (Torrilhon & Sarna, 2017; Zinner & Öttinger, 2019).

The second law of thermodynamics (entropy inequality) plays an important role in finding constitutive relations in the bulk (de Groot & Mazur, 1962; Gyarmati, 1970; Harten, 1983; Lebon et al., 2008), in designing numerical schemes (Tadmor & Zhong, 2006; Kumar & Mishra, 2012; Chandrashekar, 2013), as well as, in developing physically admissible boundary conditions (Bond & Struchtrup, 2004; Struchtrup & Torrilhon, 2007; Kjelstrup & Bedeaux, 2008; Schweizer et al., 2016; Rana & Struchtrup, 2016; Rana et al., 2018a; Sarna & Torrilhon, 2018; Beckmann et al., 2018). For the latter, one determines the entropy generation at the boundary and finds the boundary conditions as phenomenological laws that guarantee positivity of the entropy generation at the boundary. For example, in Struchtrup & Torrilhon (2007), Rana & Struchtrup (2016), Sarna & Torrilhon (2018) and Beckmann et al. (2018) the linearized second law (taking a quadratic entropy) was employed to determine boundary conditions for the linearized version of the moment equations. The resulting phenomenological boundary conditions were thermodynamically consistent for all processes and complied with the Onsager reciprocity relations (Onsager, 1931).

Over the last few years, an impressive body of research work has been devoted to the development of constitutive theories whose closed conservation laws guarantee the second law (Öttinger, 2010; Struchtrup & Torrilhon, 2010; Torrilhon, 2012; Liu et al., 2013; Paolucci & Paolucci, 2018). The Burnett equations, which are obtained perturbatively from the kinetic theory, are unstable (Bobylev, 2006) and due to lack of any coherent second law these equations lead to physically inadmissible solutions (Struchtrup & Nadler, 2020). The Burnett equation are unstable due to lack of a proper entropy inequality. On the other hand, the moment equations (Grad 13-moment equations, regularized 13-moment equations, etc.) are accompanied by proper entropy inequalities, and are stable, but only in the linear cases (Struchtrup & Torrilhon, 2007; Rana & Struchtrup, 2016; Torrilhon & Sarna, 2017). In Rana et al. (2018a), the coupled constitutive relations (CCR) were obtained—as an enhancement to the NSF equations—by considering a correction term to the non-convective entropy flux (in addition to the classic entropy flux: heat-flux/thermodynamic temperature). As a result, the CCR adds several second-order terms to the NSF equations, in the bulk and in the boundary conditions, capturing important rarefaction effects in moderately rarefied conditions—such as the Knudsen minimum, and non-Fourier heat transfer—which cannot be described by classical hydrodynamics (Rana et al., 2018a).

Through this article, we are going to utilize the CCR theory, thereby some comments about these equations are in order. The CCR theory originates from the arguments of irreversible thermodynamics. In CCR theory, additional terms appear in the constitutive relations of NSF equations, where the coupling between constitutive relations for the heat-flux vector and stress tensor is introduced via a coupling coefficient ($\alpha_{0}$). For $\alpha_{0}=0$, the coupling vanishes and one obtains the classical NSF relations. Besides, for the flow conditions considered in the article (steady state and linearized equations), the CCR system (conservation laws plus CCR) reduces to the linearized Grad 13-moment equations for $\alpha_{0}=2/5$; a value obtained for Maxwell molecule (MM) gases (Struchtrup, 2005). Thus, under suitable assumptions, this article provides a unified framework for NSF, Grad-13 and the CCR system. More precisely, we explore the method of fundamental solutions (MFS) for the CCR system (hence, NSF, Grad-13) that can be useful for constructing solutions over a wide range of free-stream profiles and complex geometries.

The MFS (Kupradze & Aleksidze, 1964) is a mesh-free numerical technique for solving partial differential equations based on using the fundamental solutions (Green’s functions). The main advantage that the MFS has over the more classical numerical methods, such as finite difference method, finite element method, is that its solution procedure does not require discretisation of the interior of the computational domain, thus a significant amount of computational effort and time is saved. Furthermore, the MFS does not require the (very computationally expensive) 3D mesh generation, as a result, the MFS is ideally suited for the problems involving moving interfaces and phase-change processes. The first steps towards utilizing MFS for rarefied gas flows were taken by Lockerby & Collyer (2016), who derived Green’s functions for the Grad 13 system and utilized them to solve some canonical problems. The present study further develops a MFS-based framework (MFS–CCR) for phase-transition boundary-value problems, which, as we shall show later, will involve derivation of new Green’s functions generated by a singular mass source/sink term. In this article, a set of temperature-jump and velocity-slip boundary conditions for a liquid-vapour interface will be derived for the linearised CCR system and implemented within the framework of MFS. However, it should be noted that the developed MFS methodology in this article can also be applied to other type of liquid-vapour interface boundary conditions, such as the classical Schrage law (Schrage, 1953), statistical rate theory (Ward & Fang, 1999) and phenomenological approach based on Hertz-Knudsen-Schrage relation (Liang et al., 2017). We do not address these boundary conditions in this article, instead leaving this for a future investigation.

The remainder of this paper is organized as follows. In Section 2, we introduce the linearized and dimensionless CCR model. In Section 3, the derivation of Green’s functions for the CCR model is presented, followed by the derivation of thermodynamically consistent boundary conditions at the phase interface in Section 4. A brief summary of the method of fundamental solutions and implementation is given in Section 5. In Section 6.1 the numerical scheme is applied to evaporation from a spherical droplet and low-speed gas flow around a sphere (for which the analytical solutions exist) in order to validate our numerical scheme. To illustrate the utility of fundamental solutions for complex geometries, in Sections 6.2, 6.3 and 6.4 we solve the motion of two spherical non-evaporating droplets with different orientations, evaporation of two interacting droplets, and evaporation of a deformed droplet, respectively. Concluding remarks are made in Section 8.

In this article we shall only consider flow conditions where the deviations from a constant equilibrium state—given by a constant reference density $\hat{\rho}_{0}$, a constant reference temperature $\hat{T}_{0}$, and all other fields zero—are small. Thereby, the governing equations can be linearized with respect to the reference state. The governining equations are put into dimensionless form by introducing

$$\mathbf{r}=\frac{\mathbf{\hat{r}}}{\hat{\ell}}\text{, \quad}\rho=\frac{\hat{\rho}-\hat{\rho}_{0}}{\hat{\rho}_{0}}\text{,\quad and }T=\frac{\hat{T}-\hat{T}_{0}}{\hat{T}_{0}}\text{,}$$

(1)

where $\mathbf{r}$ is the position vector in Cartesian coordinates, $\hat{\ell}$ is the reference length scale, and $\rho$ and $T$ are the dimensionless perturbations in the density and temperature from their reference values, respectively; hats denote dimensional quantities. The dimensionless velocity vector $\mathbf{v}$, the stress tensor $\bm{\Pi}$ and the heat-flux vector $\mathbf{q}$ are

$$\mathbf{v}=\frac{\mathbf{\hat{v}}}{\sqrt{\mathrm{\hat{R}}\hat{T}_{0}}}\text{, }\quad\bm{\Pi}=\frac{\bm{\hat{\Pi}}}{\hat{\rho}_{0}\mathrm{\hat{R}}\hat{T}_{0}}\text{,\quad and }\mathbf{q}=\frac{\mathbf{\hat{q}}}{\hat{\rho}_{0}\left(\mathrm{\hat{R}}\hat{T}_{0}\right)^{3/2}}\text{,}$$

(2)

where $\mathrm{\hat{R}}$ is the gas constant. The perturbation in pressure $p$ is given by a linearized equation of state

$$p=\frac{\hat{\rho}\hat{R}\hat{T}}{\hat{\rho}_{0}\hat{R}\hat{T}_{0}}-1=\left(1+\rho\right)\left(1+\theta\right)-1\approx\rho+\theta$$

(3)

where the last equation assumes small perturbations $\rho$ and $\theta$, hence $\rho\theta$ is assumed to be negligible in comparison to $\rho+\theta$.

Accordingly, dimensionless and steady-state (linearized) conservation laws for mass, momentum, and energy read

$$\nabla\cdot\mathbf{v}=0\text{,\quad}\nabla p+\nabla\cdot\bm{\Pi}=0\text{,\quad and }\nabla\cdot\mathbf{q}=0\text{.}$$

(4)

It should be noted that conservation laws (4) contain the stress tensor and heat flux vector as unknowns, hence constitutive equations are required for these quantities. If the CCR closure is adopted, the (linearized) constitutive relations for $\bm{\Pi}$ and $\mathbf{q}$ read (Rana et al., 2018a)

$$\bm{\Pi}=-2\mathrm{Kn}\langle\nabla\mathbf{v}\rangle-2\alpha_{0}\mathrm{Kn}\langle\nabla\mathbf{q}\rangle\text{,\quad and }\mathbf{q}=-\frac{c_{p}\mathrm{Kn}}{\Pr}\left(\nabla T+\alpha_{0}\nabla\cdot\bm{\Pi}\right)$$

(5)

where, $c_{p}=5/2(=\hat{c}_{p}/\hat{R})$ is the specific heat for mono-atomic gases, $\Pr$ is the Prandtl number and the Knudsen number $\mathrm{Kn}=\hat{\mu}_{0}/\left(\hat{\rho}_{0}\sqrt{R\hat{T}_{0}}\hat{\ell}\right)$, where $\hat{\mu}_{0}$ is the viscosity of the gas at the reference state. Furthermore, the angular brackets in (5) denote the traceless and symmetrical component of the tensor, for instance

$$\langle\nabla\mathbf{v}\rangle=\frac{1}{2}\left(\nabla\mathbf{v+}\nabla\mathbf{v}^{T}\right)-\frac{1}{3}\left(\nabla\cdot\mathbf{v}\right)\mathbf{I}\text{,}$$

(6)

where $\mathbf{I}$ is the identity matrix.

Since the constitutive relations (5) for the fluxes $\bm{\Pi}$ and $\mathbf{q}$ are coupled (in contrast to NSF), these equations are referred as the coupled constitutive relations (CCR). The benefit of the CCR model compared to extended moment equations (Grad-13 equations, for instance) is that it retains full thermodynamic structure without any restriction (Rana et al., 2018a).

In the above equations the phenomenological coefficient $\alpha_{0}$ appears, which is chosen such that the CCR model agrees with the Burnett coefficients (in the sense of the Chapman–Enskog expansion of the Boltzmann equation); the NSF equations are obtained when $\alpha_{0}=0$. The value of $\alpha_{0}$ for the Maxwell molecule gases is 2/5; for other power potentials it can be computed from the relation $\alpha_{0}=\frac{\Pr\varpi_{3}}{5}$, where $\varpi_{3}$ is the Burnett coefficient (Struchtrup, 2005; Rana et al., 2018a).

It is important to note here, in steady state, the linearized CCR model reduces to the linearized Grad’s-13 moment equations for $\alpha_{0}=2/5$ (i.e., for Maxwell molecule gases). Hence, using the naming convention of Lockerby & Collyer (2016); Claydon et al. (2017), the solution of the CCR equations (4–5) will still be referred as the Generalised Gradlet.

In this section, we shall write the fundamental solutions (Green’s functions) for the CCR system (4–5), that will be used further in constructing the new solutions for complex geometries and phase-change problems.

In this section, we shall derive the phase interface boundary conditions for the CCR system within the framework of irreversible thermodynamics (Kjelstrup & Bedeaux, 2008). In order to establish (thermodynamically consistent) boundary conditions at the phase interface, one determines the entropy generation at the interface, and finds the boundary conditions as phenomenological laws that guarantee positivity of the entropy generation rate. For the linearized CCR equations, the entropy generation rate at the interface is (Beckmann et al., 2018; Rana et al., 2018a)

	$\displaystyle\Sigma_{\text{surface}}$	$\displaystyle=-\left(\mathbf{v}-\mathbf{v}^{I}\right)\cdot\mathbf{n}\left[p-p_{sat}+\left(\mathbf{n}\cdot\bm{\Pi}\cdot\mathbf{n}\right)\right]-\mathbf{q}\cdot\mathbf{n}\left[T-T^{l}+\alpha_{0}\left(\mathbf{n}\cdot\bm{\Pi}\cdot\mathbf{n}\right)\right]$
		$\displaystyle\quad-\sum_{i=1}^{2}\left(\mathbf{t}^{(i)}\cdot\bm{\Pi}\cdot\mathbf{n}\right)\left[\mathbf{t}^{(i)}\cdot\left(\mathbf{v-v}^{l}\right)+\alpha_{0}\mathbf{t}^{(i)}\cdot\mathbf{q}\right]\geq 0\text{,}$		(15)

where $\mathbf{n}$ is a unit normal pointing from a boundary point into the gas, and $\left\{\mathbf{t}^{(i)}\right\}_{i=1,2}$are orthonormal tangent vectors to the boundary. Here, $\mathbf{v}^{I}$ denotes the velocity of the interface, $T^{I}$ denotes the temperature of the liquid at the interface and $p_{sat}{\left(T^{I}\right)}$ is the saturation pressure.

The positivity of entropy generation rate $\Sigma_{surface}$ is ensured by adopting the following boundary equations

	$\displaystyle\left(\mathbf{v}-\mathbf{v}^{I}\right)\cdot\mathbf{n}$	$\displaystyle=$	$\displaystyle-\eta_{11}\left(p-p_{sat}+\mathbf{n}\cdot\bm{\Pi}\cdot\mathbf{n}\right)+\eta_{12}\left(T-T^{I}+\alpha_{0}\mathbf{n}\cdot\bm{\Pi}\cdot\mathbf{n}\right)$		(16a)
	$\displaystyle\mathbf{q}\cdot\mathbf{n}$	$\displaystyle=$	$\displaystyle\eta_{12}\left(p-p_{sat}+\mathbf{n}\cdot\bm{\Pi}\cdot\mathbf{n}\right)-(\eta_{22}+2\tau_{0})\left(T-T^{I}+\alpha_{0}\mathbf{n}\cdot\bm{\Pi}\cdot\mathbf{n}\right)$		(16b)

and

	$\displaystyle\mathbf{t}^{(1)}\cdot\bm{\Pi}\cdot\mathbf{n}$	$\displaystyle=$	$\displaystyle-\varsigma\left(\mathbf{v-v}^{I}+\alpha_{0}\mathbf{q}\right)\cdot\mathbf{t}^{(1)}\text{,}$		(17a)
	$\displaystyle\mathbf{t}^{(2)}\cdot\bm{\Pi}\cdot\mathbf{n}$	$\displaystyle=$	$\displaystyle-\varsigma\left(\mathbf{v-v}^{I}+\alpha_{0}\mathbf{q}\right)\cdot\mathbf{t}^{(2)}\text{.}$		(17b)

Through boundary conditions (16), the evaporative mass and heat flux are governed by the difference between pressure and saturation pressure, and temperature difference across the interface. Furthermore, the boundary conditions (17a-17b) relates the shear stress to the tangential velocity slip and thermal transpiration—a flow induced by a tangential heat flux—which is a second-order effect (Sone, 2007; Hadjiconstantinou, 2003; Cercignani & Lorenzani, 2010). The boundary conditions (16-17) are an extension to the classical Hertz-Knudsen-Schrage law for evaporation (Schrage, 1953), taking into account the temperature jump, velocity slip and some second-order effects.

In boundary conditions (16), $\eta_{ij}$’s are Onsager resistivitity coefficients, which are obtained from the asymptotic ($\mathrm{Kn}\rightarrow 0$) kinetic theory (Sone, 2007), as (see appendix B for details)

$$\eta_{11}=0.9134\sqrt{\frac{2}{\pi}}\frac{\vartheta}{2-\vartheta}\text{, }\eta_{12}=0.3915\sqrt{\frac{2}{\pi}}\frac{\vartheta}{2-\vartheta}\text{, and }\eta_{22}=0.1678\sqrt{\frac{2}{\pi}}\frac{\vartheta}{2-\vartheta}\text{.}$$

(18)

These coefficients were obtained assuming evaporation/condensation coefficient $\vartheta$ is independent of the impact energy of molecules and that all vapour molecules that are not condensing are thermalized, i.e., the accommodation coefficient, $\chi=1$. The temperature-jump coefficient $\tau_{0}=0.8503\sqrt{2/\pi}$ and velocity-slip coefficient $\varsigma=0.8798\sqrt{2/\pi}$.. Throughout this article, we shall take $\vartheta=0$ for the canonical boundaries and $\vartheta=1$ (a value largely accepted in literature) for the phase-change boundaries.

Although in the development of our MFS approach, we use equation (16-17) as the boundary conditions, this methodology is in principle capable of accommodating other type of phase-interface boundary conditions, such as statistical rate theory (Ward & Fang, 1999) and phenomenological approach based on Hertz-Knudsen-Schrage relation (Liang et al., 2017). We do not address these boundary conditions in this article but should be worthy of future investigation.

In next sections, we will introduce, and use, the method of fundamental (MFS) solutions for CCR (CCR-MFS), which is computationally economical, and show that it provides reliable solutions of the CCR and NSF equations in good agreement to the known closed-form analytical results.

Let us consider $N_{c}$ collocation points on the boundary at $\left\{\mathbf{r}^{c(j)}\right\}_{j=1}^{N_{c}}$, while $N_{s}$ singularity points are located outside the computational domain (i.e. inside the solid/liquid body) with position vectors $\left\{\mathbf{r}^{s(j)}\right\}_{j=1}^{N_{s}}$; see Figure 1.

The flow-field variables are given by a superposition of the Gradlets, thermal Gradlets and sourcelets, formulated in (13-14), as

	$\displaystyle\mathbf{v}^{(i)(Gr)}$	$\displaystyle=$	$\displaystyle\sum_{j=1}^{N^{s}}\frac{1}{8\pi}\mathbf{J}{\left(\mathbf{r}^{(ij)}\right)}\cdot\mathbf{f}^{(j)}+\frac{3\alpha_{0}^{2}c_{p}\mathrm{Kn}^{2}}{4\pi\Pr}\mathbf{K}{\left(\mathbf{r}^{(ij)}\right)}\cdot\mathbf{f}^{(j)}$		(19a)
	$\displaystyle\mathbf{v}^{(i)(S)}$	$\displaystyle=$	$\displaystyle\sum_{j=1}^{N^{s}}\frac{h^{(j)}}{4\pi}\frac{\mathbf{r}^{\left(ij\right)}}{|\mathbf{r}^{\left(ij\right)}|^{3}}\text{,}$		(19b)
	$\displaystyle\mathbf{v}^{(i)}$	$\displaystyle=$	$\displaystyle\mathbf{v}^{(i)(Gr)}+\mathbf{v}^{(i)(S)}\text{,}$		(19c)
	$\displaystyle p^{(i)}$	$\displaystyle=$	$\displaystyle\sum_{j=1}^{N^{s}}\frac{\mathrm{Kn}}{4\pi}\frac{\mathbf{f}^{(j)}\mathbf{\cdot r}^{\left(ij\right)}}{|\mathbf{r}^{\left(ij\right)}|^{3}}\text{,}$		(19d)
	$\displaystyle\bm{\Pi}^{(i)}$	$\displaystyle=$	$\displaystyle\sum_{j=1}^{N^{s}}\frac{3\mathrm{Kn}}{4\pi}\left(\mathbf{f}^{(j)}\mathbf{\cdot r}^{\left(ij\right)}+2\alpha_{0}\mathrm{Kn}g^{(j)}+2h^{(j)}\right)\mathbf{K}{\left(\mathbf{r}^{(ij)}\right)}\text{,}$		(19e)
and

	$\displaystyle T^{(i)}$	$\displaystyle=$	$\displaystyle\sum_{j=1}^{N^{s}}\frac{\Pr}{4\pi c_{p}}\frac{g^{(j)}}{|\mathbf{r}^{\left(ij\right)}|}\text{,}$		(20a)
	$\displaystyle\mathbf{q}^{(i)}$	$\displaystyle=$	$\displaystyle\sum_{j=1}^{N^{s}}\frac{\mathrm{Kn}}{4\pi}\frac{g^{(j)}\mathbf{r}^{\left(ij\right)}}{|\mathbf{r}^{\left(ij\right)}|^{3}}-\frac{3\alpha_{0}c_{p}\mathrm{Kn}^{2}}{4\pi\Pr}\mathbf{K}{\left(\mathbf{r}^{(ij)}\right)}\cdot\mathbf{f}^{(j)}\text{.}$		(20b)

are fundamental solutions for the CCR equations. Here, $\mathbf{r}^{\left(ij\right)}=\mathbf{r}^{c\left(i\right)}-\mathbf{r}^{s\left(j\right)}$ are the displacement vectors from the $j$th singularity site to the $i$th collocation point.

There are total of $5\times N^{s}$ unknowns in (19-20)—$5$ at each singularity site—which need to be specified via suitable boundary conditions. The $4$ boundary conditions in (16–17) need to be satisfied at every collocation point. Hence, it turns out, one additional boundary condition is required (at each collocation point), we use the following boundary conditions:

	$\displaystyle\mathbf{v}^{(i)(S)}\cdot\mathbf{n}^{(i)}$	$\displaystyle=-\eta_{11}\left(p^{(i)}-p_{sat}^{(i)}+\mathbf{n}^{(i)}\cdot\bm{\Pi}^{(i)}\cdot\mathbf{n}^{(i)}\right)$
		$\displaystyle\quad+\eta_{12}\left(T^{(i)}-T^{(i)I}+\alpha_{0}\mathbf{n^{(i)}}\cdot\bm{\Pi}^{(i)}\cdot\mathbf{n}^{(i)}\right)\text{,}$		(21a)
	$\displaystyle\left(\mathbf{v}^{(i)(Gr)}-\mathbf{v}^{\left(i\right)I}\right)\cdot\mathbf{n}^{\left(i\right)}$	$\displaystyle=0\text{,}$		(21b)

Notably, we have split the boundary condition (16a) into two components (21a) and (21b). Clearly, these two boundary conditions are sufficient for condition (16a) to hold. The boundary conditions for the normal heat-flux (16b) and shear-stress (17) read

	$\displaystyle\mathbf{q}^{(i)}\cdot\mathbf{n}^{(i)}$	$\displaystyle=\eta_{12}\left(p^{(i)}-p_{sat}^{(i)}+\mathbf{n}^{(i)}\cdot\bm{\Pi}^{(i)}\cdot\mathbf{n}^{(i)}\right)$
		$\displaystyle\quad-\left(\eta_{22}+2\tau_{0}\right)\left(T^{(i)}-T^{(i)I}+\alpha_{0}\mathbf{n^{(i)}}\cdot\bm{\Pi}^{(i)}\cdot\mathbf{n}^{(i)}\right)\text{,}$		(22a)
	$\displaystyle\mathbf{t}^{(1)\left(i\right)}\cdot\bm{\Pi}^{\left(i\right)}\cdot\mathbf{n}^{\left(i\right)}$	$\displaystyle=-\varsigma\left(\mathbf{v}^{\left(i\right)}\mathbf{-v}^{\left(i\right)I}+\alpha_{0}\mathbf{q}^{\left(i\right)}\right)\cdot\mathbf{t}^{(1)\left(i\right)}\text{,}$		(22b)
	$\displaystyle\mathbf{t}^{(2)\left(i\right)}\cdot\bm{\Pi}^{\left(i\right)}\cdot\mathbf{n}^{\left(i\right)}$	$\displaystyle=-\varsigma\left(\mathbf{v}^{\left(i\right)}\mathbf{-v}^{\left(i\right)I}+\alpha_{0}\mathbf{q}^{\left(i\right)}\right)\cdot\mathbf{t}^{(2)\left(i\right)}\text{.}$		(22c)

For convenience, we write the set of equations (21-22) in matrix form :

$$\sum_{j=1}^{N^{s}}\mathbf{\mathscr{L}}^{ij}\mathbf{u}^{\left(j\right)}=\mathbf{b}^{\left(i\right)}\text{,}$$

(23)

where $\mathbf{u}^{\left(j\right)}$ is the vector containing the $5$ degrees of freedom at each singularity point, i.e.,

$$\mathbf{u}^{\left(j\right)}=\left\{g^{\left(j\right)}\text{, }f_{1}^{\left(j\right)}\text{, }f_{2}^{\left(j\right)}\text{, }f_{3}^{\left(j\right)}\text{, }h^{\left(j\right)}\right\}\text{,}$$

Where $\mathscr{L}^{ij}$ is a $5\times 5$ coefficient matrix, and $\mathbf{b}^{\left(i\right)}$ is a $5\times 1$ vector, which can be readily extracted using a symbolic software; we used Mathematica® for this. The inhomogeneous part $\left\{\mathbf{b}^{\left(i\right)}\right\}_{i=1}^{N^{c}}$ contains the properties of the interface, namely, $T^{(i)I}$, $p^{(i)}_{sat}$ and $\mathbf{v}^{\left(i\right)I}$. The $5N^{c}\times 5N^{s}$ linear system (23) is solved for $\left\{\mathbf{u}^{\left(i\right)}\right\}_{i=1}^{N^{s}}$ using an iterative quasi-minimal residual method (Freund & Nachtigal, 1991). Moreover, the normal and the tangent vectors $\left\{\mathbf{n}^{(j)}\text{, }\mathbf{t}^{(1)\left(j\right)}\text{, }\mathbf{t}^{(1)\left(j\right)}\right\}_{j=1}^{N^{c}}$ are obtained using Householder reflection formula for vector orthogonalization (Lopes et al., 2013).

The CCR-MFS described above involves some numerical parameters specifically, how many collocation ($N^{c}$) and singularity ($N^{s}$) points to take, and the location of singularity points outside of the computational domain. These parameters govern the overall efficiency of the numerical scheme, by striking a balance between numerical error and computational time. Throughout this study, we consider an equal number of collocation and singularity points, i.e., $N^{c}=N^{s}$ and take $\gamma\in\left(0,1\right)$ as a geometry-dependent parameter, which governs how far the singularity points are from the boundary. For instance, for the spherical geometry shown in Figure 1, we choose $\gamma=R_{s}/R_{c}$.

In order to validate our code and to establish the numerical accuracy of the CCR-MFS scheme, we validate our results for a case of an evaporating spherical droplet, the analytic solutions to this problems are readily available in (Rana et al., 2018b) for the case of Grad-13 equations. Additionally, in order to establish the accuracy of our models, we consider a slow flow around a sphere (evaporative and non-evaporative) for various values of Knudsen number and compare our solutions with results from kinetic theory.