6.1. Validation and verification of CCR–MFS

6.1.1. Validation case: evaporation from a spherical droplet

First, we consider a liquid droplet of fixed radius with a given interface temperature $T^{I}$ and the corresponding saturation pressure $p_{sat}$, immersed in its own vapour, see figure 2. The far-field conditions are given by $T_{\infty }=0$, $p_{\infty }=0$, that is, we consider the far-field conditions and the droplet radius ($R_0$) for non-dimensionalisation. Throughout this paper we do not consider dynamics within the droplets, surface tension, etc., see Rana et al. (Reference Rana, Lockerby and Sprittles2019).

Figure 2.

Schematic representation of a liquid droplet surrounded by its own vapour. The equations are made dimensionless with respect to the far-field conditions ($T_{\infty }$, $p_{\infty }$). The temperature and the pressure deviations vanish as $r\rightarrow \infty$.

Assuming the spherical symmetry of this problem, the analytic solutions for the radial velocity $v_{r}$, heat flux $q_{r}$, normal stress $\varPi _{rr}$ and the temperature $T$ are obtained from (2.4a–c) and (2.5a,b) as (Rana et al. Reference Rana, Lockerby and Sprittles2018b)

(6.1a–e)\begin{equation} v_{r}=\frac{c_{1}}{r^{2}}{,\quad }q_{r}=\frac{c_{2}}{r^{2}},\quad \varPi _{rr}=\frac{4{Kn}\left( c_{1}+\alpha _{0}c_{2}\right) }{r^{3}}{,}\quad T=\frac{Pr }{c_{p}{Kn}}\frac{c_{2}}{r},\quad \text{and}\quad p=0{.} \end{equation}

Here, $r$ is the radial direction and $c_{1}$ and $c_{2}$ are constants of integration. Solutions (6.1a–e) can also be derived from the fundamental solutions (3.7a)–(3.7c), which are given in Appendix C.

The flow is driven by (i) a unit pressure difference while the temperature of the liquid is equal to the far-field temperature (i.e. $p_{sat}=1$ and $T^{I}=T_{\infty }=0$) or (ii) a unit temperature difference while the saturation pressure in the liquid is the same as the far-field pressure (i.e. $p_{sat}=p_{\infty }=0$ and $T^{I}=1$). In both cases, the constants of integration, namely, the mass flux ($c_{1}$) and heat flux ($c_{2}$), per unit area, are obtained from boundary conditions (4.2) at $r=1$. For Grad-13 (i.e. $\alpha _0=2/5$) these are obtained as

(6.2a)\begin{gather} c_{1}^{p}=0.220562+\frac{0.0613424+0.0960246{Kn}}{0.138548+{Kn}(0.716085+{Kn})}{,} \end{gather}

(6.2b)\begin{gather}c_{2}^{p}=c_{1}^{\tau}={-}\frac{{Kn}(0.108868+0.551404{Kn})}{ 0.138548+{Kn}(0.716085+{Kn})}{,} \end{gather}

(6.2c)\begin{gather}c_{2}^{\tau}=\frac{0.152568{Kn}+{Kn}^{2}(0.924355+1.37851 {Kn}))}{0.0406847+{Kn}(0.348827+{Kn}(1.00974+{Kn}))}{.} \end{gather}

Here, superscripts ‘$p$’ and ‘$\tau$’ denote the pressure-driven case ($p_{sat}=1$ and $T^{l}=T_{\infty }=0$) and temperature-driven case ($p_{sat}=p_{\infty }=0$ and $T^{I}=1$), respectively. The results of these two problems can be combined to evaluate the total evaporative mass and heat flux from the droplet (Rana et al. Reference Rana, Lockerby and Sprittles2019). Owing to the microscopic reversibility of the evaporation and condensation processes the Onsager reciprocity relations hold, which give $( c_{2}^{p}=c_{1}^{\tau })$ (Chernyak & Margilevskiy Reference Chernyak and Margilevskiy1989; Rana et al. Reference Rana, Lockerby and Sprittles2018b).

The mass and the heat flux predicted by the CCR–MFS method are obtained by integrating $\boldsymbol {v}^{(Gr)}$ (5.1a), (5.1b) and $\boldsymbol {q}^{(Gr)}$ (5.2b) over the droplet surface, i.e.

(6.3)\begin{gather} c_{1}^{(MFS)} =\frac{1}{4{\rm \pi} }\int_{0}^{{\rm \pi} }\int_{0}^{2{\rm \pi} }\boldsymbol{v} ^{(Gr)}\boldsymbol{\cdot} \boldsymbol{n}\,\textrm{d}\varphi \,\textrm{d}\theta =\frac{1}{4{\rm \pi} } \sum_{i=1}^{N^{s}}h^{(i)}{,} \end{gather}

(6.4)\begin{gather}c_{2}^{(MFS)} =\frac{1}{4{\rm \pi} }\int_{0}^{{\rm \pi} }\int_{0}^{2{\rm \pi} }\boldsymbol{q} ^{(Gr)}\boldsymbol{\cdot} \boldsymbol{n}\,\textrm{d}\varphi \,\textrm{d}\theta =\frac{{Kn}}{4{\rm \pi} } \sum_{i=1}^{N^{s}}g^{(i)}{.} \end{gather}

Let us first briefly consider the convergence characteristics of the CCR–MFS method. The configuration for the collocation and singularity points is shown in figure 1. We define errors $\epsilon$ between the numerical prediction of $c_1^{p}$, $c_1^{T}$ and $c_2^{T}$ to the analytical values (6.2) as

(6.5a–c)\begin{equation} \epsilon\left(c_1^{p}\right) = c_1^{p\left(MFS\right)}-c_1^{p},\quad \epsilon\left(c_2^{\tau}\right) = c_2^{\tau\left(MFS\right)}-c_2^{\tau} \quad \text{and} \quad \epsilon\left(c_1^{\tau}\right) = c_1^{\tau\left(MFS\right)}-c_1^{\tau}{.} \end{equation}

Figure 3 shows the numerical errors of the MFS (with $\alpha _0=2/5$) in $c_1^{p}$, $c_2^{\tau }$ (a) and $c_1^{\tau }$ (b) for three different singularity site locations ($\gamma =R_s/R_c=0.05, 0.1, 0.25$) and various values of the Knudsen number. The result indicates very high accuracy with a moderate number of points ($N^c=112$). The best results are obtained when the singularity sites are farther away from the collocation nodes, i.e. $\gamma$ is small. However, it leads to poorer conditioning of the equation system (5.5), which in turn results in larger least-square errors and computational time.

Figure 3.

Log–log plot of errors against Knudsen number with three different singularity locations $\gamma =R_s/R_c=0.05$, 0.1, 0.25 and with $N^c=N^s=112$: (a) error in $c_1^{p}$ and $c_2^{\tau }$ and (b) error in $c_1^{\tau }$.

As Lockerby & Collyer (Reference Lockerby and Collyer2016), numerical tests were also conducted by increasing the number of collocation points, as well as, by shifting the singularities sites (i.e.breaking the symmetry of the collocation and singularity points), the results were found to be satisfactory. Henceforth, the numerical simulations are performed with $\gamma =0.1$, unless otherwise stated.

Figure 4 shows the variations in the mass- and heat-flux coefficients ($c_1$ and $c_2$) for a range of Knudsen numbers, the numerical parameters are included in the caption. In order to put the CCR models into perspective, we compare our predictions with different theories, namely, NSF ($\alpha _0=0$), Grad-13 ($\alpha _0=2/5$) and ($\alpha _0=3/5$), against the results (symbols) from the linearised Boltzmann equations (LBEs) by Takata et al. (Reference Takata, Sone, Lhuillier and Wakabayashi1998). The solutions of the NSF equations with classical boundary conditions, i.e. assuming the (linearised) Hertz–Knudsen–Schrage law for evaporation (Schrage Reference Schrage1953) and the continuity of temperature across the interface, are also included (thin green lines) in figure 4 for the comparison. These boundary conditions can be obtained from (4.2a) by taking $\eta _{12}=0$ and $T=T^{I}$ instead of (4.2a), which are only valid for ${Kn}\rightarrow 0$ (Sone Reference Sone2007; Rana et al. Reference Rana, Lockerby and Sprittles2018b).

Figure 4.

Mass- and heat-flux coefficients as a function of the Knudsen number: (a) mass flux ($c_1^{p}$) in pressure-driven case and (b) normalised mass ($c_1^{\tau }/{Kn}$) and heat flux ($4c_2^{\tau }/15{Kn}$) coefficients in the temperature-driven case. The results obtained from MFS with different theories are compared. The symbols denoting the results from Takata et al. (Reference Takata, Sone, Lhuillier and Wakabayashi1998) and (green) thin line representing the analytic results from classical boundary conditions. The number of collocation nodes used in each case are $N^c=N^s=112$ and $\gamma =0.1$.

For the pressure-driven case (cf. figure 4a) mass flux goes from liquid to vapour (i.e. $c_1^{p}\geqslant 0$) and heat flows from vapour to liquid (i.e. $c_2^{p} = c_1^{\tau }\leqslant 0$) because the enthalpy of phase change must be supplied to the droplet to keep its temperature constant. All models agree in the limit of ${Kn}\rightarrow 0$ for the mass flux, the NSF equations with temperature-jump boundary conditions give better results over classical theories as the Knudsen number grows. Interestingly, the NSF equations with jump boundary conditions predicts a zero mass flux in the free-molecular regime, these results contradict kinetic-theory predictions. On the other hand, the Grad-13 theory gives a non-zero mass flux, but still predicts almost half of the mass flux of the LBE result. By performing an asymptotic expansion and matching the mass flux as ${Kn}\rightarrow \infty$, we find $\alpha _0=3/5$ (a value extensively used in this study), with this model consequently providing the best agreement with the LBEs.

For the temperature-driven case (cf. figure 4b), Fourier's law with a no-jump boundary condition gives heat flux $c_2^{\tau } = 15{Kn}/4$, while the cross effects, i.e. the heat flux owing to pressure difference or the mass flux owing to temperature difference, vanish. We have used this value to normalise the heat-flux computation in figure 4(b). All models, except the no-jump boundary conditions (thin green lines), match reasonably well with the LBE data, with $\alpha _0=3/5$ giving a slight improvement over other models. However, the cross effects ($c_1^{\tau }$) are not well captured by CCR theory, which requires resolution of the Knudsen layer, these are beyond the capabilities of the models considered in this paper, see Rana et al. (Reference Rana, Lockerby and Sprittles2018b) for a discussion on this.

6.1.2. Slow flow around a rigid sphere

Let us now consider the case of low-speed rarefied gas flow around a rigid ($\vartheta =0$) spherical stationary particle of radius $R_{0}$. The particle is assumed to be isothermal, that is, the solid-to-gas conductivity ratio is large, with $T^{I}=0$, i.e. the far-field temperature and the temperature of the particle are same. The net force exerted on the sphere by gas is defined as

(6.6)\begin{equation} \boldsymbol{F}^{(Gr)}=\int_{0}^{{\rm \pi} }\int_{0}^{2{\rm \pi} }\left( p^{(Gr)}\boldsymbol{I}+ \boldsymbol{\varPi }^{(Gr)}\right) \boldsymbol{\cdot} \boldsymbol{n}\,\textrm{d}\varphi \,\textrm{d}\theta ={Kn} \sum_{i=1}^{N^{s}}\boldsymbol{f}^{(i)}{.} \end{equation}

The drag force is the projection of the net force (6.6) on to the streamwise direction. The Stokes formula for drag force (in the direction of flow) exerted by a spherical particle on the fluid flow has the form Lamb (Reference Lamb1945)

(6.7)\begin{equation} F^{S}={-}6{\rm \pi} {Kn}u_{\infty }{,} \end{equation}

where $u_\infty$ is the far-field velocity.

The Stokes formula is only valid for ${Kn}\rightarrow 0$, and requires corrections at finite Knudsen number. Figure 5 shows the normalised (with Stokes’ drag (6.7)) drag coefficient versus the Knudsen number. As expected, all theories agree in the small Knudsen limit. As the value of Knudsen number increases, the normalised drag decreases owing to increasing slip (which is not embedded in the Stokes formula). Notably, the normalised drag force reaches a finite value from the NSF equations with velocity-slip boundary conditions as $Kn\rightarrow \infty$, whereas the extended theories (CCR) predict a vanishing drag in this limit, an observation tantamount to experimental results. Nevertheless, between Grad-13 ($\alpha _0=2/5$) and the CCR, with phenomenological coefficient $\alpha _0=3/5$, the latter gives the best quantitative agreement for values of the Knudsen number beyond unity (i.e. well beyond its apparent limits of applicability).

Figure 5.

The normalised drag force computed from the CCR–MFS, with $\alpha _0 = 0$ (NSF), $\alpha _0 = 2/5$ (Grad-13) and $\alpha _0 = 3/5$ against Knudsen number. The experimental data of Millikan (1923) (fitted by Allen & Raabe (1982)) is plotted (circles) for comparison. The results are normalized with Stokes’ drag (6.7).

6.1.3. Slow flow around a spherical liquid droplet

In this section, we consider uniform flow (say, along the $z$-direction) of a saturated vapour over a spherical liquid droplet with a uniform temperature. The far-field temperature of the surrounding vapour is same as that of the droplet, i.e. $T^{I}=T_{\infty }=0$, $p_{sat}=p_{\infty }=0$ and $\boldsymbol {v}_{\infty }=-\boldsymbol {v}^{I} = \{0,0,u_{\infty }\}$. We further assume that the droplet size remains fixed (quasi-steady-state assumption) and that deformation of the droplet is negligible, that is, the capillary number is very small.

Figure 6 shows the normalised drag force on a liquid droplet versus Knudsen number, while allowing phase-change at the interface. This problem has been considered by Sone et al. (Reference Sone, Takata and Wakabayashi1994) from the LBE, and is clearly important from an engineering point of view. Owing to the motion of the vapour phase, evaporation and condensation take place on the surface of the droplet, thereby slightly reducing the drag on the sphere, cf. figures 5 and 6. Curiously, all theories, except for that of Stokes, give reasonably accurate results. Note that, in this case, the net mass and heat flux from the droplet are zero with evaporation on the front side of the sphere and condensation on the back.

Figure 6.

The normalised drag force on a liquid droplet versus Knudsen number. The results computed from the CCR–MFS, with $\alpha _0 = 0$ (NSF), $\alpha _0 = 2/5$ (Grad-13) and $\alpha _0 = 3/5$ are compared with LBE data (symbols) of Sone, Takata & Wakabayashi (Reference Sone, Takata and Wakabayashi1994).

The problems considered in previous sections typically allow for analytic solutions, which are often not available for various practical problem of interest. In the following sections, we develop cases of increasing complexity, starting first with the motion over two particles at a finite distance from each other.

6.2. Motion over two solid spheres

Let us consider a flow over two solid spheres of equal radius with both fixed in shape and size. The flow is characterised by the following parameters: radii of the spheres (i.e. the Knudsen number), centre-to-centre distance $2l_C$ and the angle between the flow direction and the centreline; see figure 7.

Figure 7.

Schematic for the flow around two solid spheres. Here, $2l_C$ is the centre-to-centre distance and is the angle between flow direction and the line joining the centres.

For an axisymmetric case (), the theoretical result for the drag force (which is equal for both the spheres) calculated by Stimson, Jeffery & Filon (Reference Stimson, Jeffery and Filon1926) reads

(6.8)\begin{equation} F^{S}={-}6{\rm \pi} Kn u_{\infty}\beta, \end{equation}

where $\beta$ is a correction coefficient which may be written in the form

(6.9)\begin{equation} \beta =\frac{4}{3}\sinh{\alpha} \sum_{k=1}^{\infty }\frac{n\left( n+1\right) }{\left( 2n-1\right) \left( 2n+3\right) }\left( 1-\frac{4\sinh ^{2}{\left( n+ \dfrac{1}{2}\right) \alpha} -\left( 2n+1\right) ^{2}\sinh ^{2}\alpha }{2\sinh \left( 2n+1\right) \alpha +\left( 2n+1\right) \sinh 2\alpha }\right), \end{equation}

with $\alpha = \cosh l_{C}$.

In figure 8(a), the drag force (normalised with Stokes’ drag (6.7) for the single sphere) for an axisymmetric case () is plotted as a function of centre-to-centre distance for different values of the Knudsen number. The results from Stimson et al. (Reference Stimson, Jeffery and Filon1926), which are only valid in the limiting case $Kn \rightarrow 0$, are included for comparison ($\diamond$), and as expected, for small values of Knudsen number ($Kn = 10^{-3}$), our results match with Stimson et al. (Reference Stimson, Jeffery and Filon1926). The experimental results for $l_C = 1/2$ obtained by Cheng et al. (Reference Cheng, Allen, Gallegos, Yeh and Peterson1988) (and given by the empirical formula (6.10)) for $Kn = 10^{-3}$, $10^{-1}$ and $0.5$ are also included ($\boldsymbol {\times }$) for comparison. For small Knudsen numbers ($Kn \lesssim 10^{-1}$), all models (NSF, Grad-13, CCR) agree, whereas at large Knudsen numbers these are markedly different, with the NSF model over-predicting the drag compared with Grad-13 and CCR with $\alpha _0=3/5$. The effect of proximity is clearly visible on drag force, which decreases as $l_C$ is reduced (i.e. the spheres gets closer). Furthermore, in the limit $l_C\rightarrow \infty$ the drag approaches that of the single sphere case (cf. figure 5). In addition, as with the case of a single sphere (cf. figures 5 and 8a), the force monotonically decreases with increases in Knudsen number, for all models.

Figure 8.

The drag force on two equal spheres versus centre-to-centre distance: (a) flow direction parallel to their line of centres and (b) flow direction perpendicular to their line of centres. The results computed from the CCR–MFS, with $\alpha _0 = 0$ (NSF), $\alpha _0 = 2/5$ (Grad-13) and $\alpha _0 = 3/5$ are compared with the Stokes solution ($\diamond$) derived by Stimson et al. (Reference Stimson, Jeffery and Filon1926). The experimental results (denoted by $\boldsymbol {\times }$) for a doublet ($l_C = 1/2$) obtained by Cheng et al. (Reference Cheng, Allen, Gallegos, Yeh and Peterson1988) are included for comparison.

For the non-axisymmetric case (), the normalised drag force is plotted in figure 8(b) as a function of $l_C$, with $Kn = 10^{-3}$, $10^{-1}$ and $0.5$. The experimental results of Cheng et al. (Reference Cheng, Allen, Gallegos, Yeh and Peterson1988) for a doublet ($l_C = 1/2$) are denoted by symbols ($\boldsymbol {\times }$). The overall drag force in the non-axisymmetric case is higher compared with the axisymmetric case (cf. figure 8a,b) and as before, the NSF theory over-predicts the drag.

6.2.1. Drag on a doublet

An interesting case arises when one considers $l_C\rightarrow 1/2$, i.e. when both spheres are in contact with one another, known as a doublet. Cheng et al. (Reference Cheng, Allen, Gallegos, Yeh and Peterson1988) carried out experiments to study the effects of orientation (modulated by an electric field) on the measured drag force.

Figure 9 shows a comparison of the different models, computed using CCR–MFS, over a wide range of the Knudsen number. Here, the experimental data for the drag force acting on a doublet are taken from the empirical formula from Cheng et al. (Reference Cheng, Allen, Gallegos, Yeh and Peterson1988), as

(6.10)\begin{equation} F ={-}\frac{4{\rm \pi} Knu_{\infty}}{2^{1/3}}\frac{\phi_0}{1+\dfrac{Kn}{\phi_0}\left[ 1.142+0.558\exp\left({-}0.999\dfrac{\phi_0}{Kn}\right) \right]}{,} \end{equation}

where $\phi _0 = 1.21$ or $\phi _0 = 1.37$ for flow moving in parallel () or in perpendicular () direction to the centreline, respectively.

Figure 9.

The drag force on a doublet: (a) flow parallel to its line of centres (axisymmetric case) and (b) flow perpendicular to its line of centres (non-axisymmetric case). The experimental data are taken from the empirical formula (6.10) obtained in Cheng et al. (Reference Cheng, Allen, Gallegos, Yeh and Peterson1988).

As shown in figure 9, reassuringly for $Kn \rightarrow 0$, the results converges to the without-slip solution of the Stokes equation. For larger values of the Knudsen number, the force decreases owing to slip, and as $Kn \rightarrow \infty$, it reaches a finite value for the NSF theory; while from (6.10) it vanishes. Remarkably, from the results of the Grad-13 and the CCR models, we find that the force indeed vanishes for large Knudsen numbers. The results obtained from the Grad-13 and CCR models are both in good agreement with the experimental results, a sightly better match at larger Knudsen number is obtained by the CCR with $\alpha _0=3/5$.

Figures 10 and 11 show the streamlines and the speed contours for the case $Kn = 0.1$ and $Kn = 0.5$, respectively. The horizontal and vertical axes represent $x$- and $z$-axis, respectively, with flow in the $x$-direction. The left-hand side of these figures shows the Grad-13 ($\alpha _0=2/5$) results for the case , and the right-hand side shows the results computed from the NSF ($\alpha _0=0$) equations. Owing to the linearity of the equations, the speed contours are symmetric about the $x=0$ plane, resulting in equal drag forces on both spheres.

Figure 10.

Streamlines and speed contours from the (a) Grad-13 and (b) NSF equations for Knudsen number 0.1: flow over a doublet with flow direction along the centreline (). Vertical and horizontal axes represent $x$- and $y$-directions ($N_c = N_s = 650$ and $\gamma = R_s/R_c = 0.1$).

Figure 11.

Streamlines and speed contours from the (a) Grad-13 and (b) NSF equations for Knudsen number 0.5: flow over a doublet with flow direction along the centreline (). Vertical and horizontal axes represent $x$- and $y$-directions ($N_c = N_s = 650$ and $\gamma = R_s/R_c = 0.1$).

These contours show that, as to be expected, the slip velocity at the bottom and top surfaces increases with Knudsen number, resulting in reduced drag forces. For $Kn = 0.1$ (figure 10), the flow line pattern predicted by the Grad-13 (a) and NSF (b) are very similar. The slip velocity is maximum near the top and bottom surface, which is virtually the same for the NSF and Grad-13 theories, thus giving similar predictions for the drag reduction (cf. figure 8). However, for $Kn = 0.5$ (figure 11), the Grad-13 theory predicts a larger slip and thus greater drag reduction. The results from the CCR with $\alpha _0=3/5$ are similar to those obtained from the Grad-13 theory, hence these are not shown here.

In figures 12 and 13, we again show the streamlines and speed contours predicted by both theories at $Kn=0.1$ and $Kn=0.5$, respectively, for the case when flow (in the $z$-direction) is in a perpendicular direction to the centreline, i.e. . Unlike the previous case, this is a truly three-dimensional set-up, where no symmetry exists in the azimuthal direction, with the MFS it is remarkably simple to solve this fully three-dimensional flow. The slip velocity increases with increasing the Knudsen number thus reducing the drag coefficient.

Figure 12.

Streamlines and speed contours from the (a) Grad-13 and (b) NSF equations for Knudsen number 0.1: flow over a doublet with flow direction perpendicular to the centreline (). Vertical and horizontal axes represent $x$- and $y$-directions ($N_c = N_s = 650$ and $\gamma = R_s/R_c = 0.1$).

Figure 13.

Streamlines and Mach contours from the (a) Grad-13 and (b) NSF equations for Knudsen number 0.5: flow over a non-evaporative ($\vartheta = 0$) doublet with flow direction perpendicular to the centreline (). Vertical and horizontal axes represent $x$- and $y$-directions ($N_c = N_s = 650$ and $\gamma = R_s/R_c = 0.1$).

6.3. Interaction of two evaporating droplets

The evaporation dynamics of two interacting droplets suspended in vapour will be investigated here; a problem relevant to dense-spray applications. We consider two droplets placed is a saturated vapour with centre-to-centre distance $2l_C$. Again, we consider the shape, size and the surface temperature of the droplet to be fixed.

As before, we consider two cases, i.e. (i) evaporation is driven by temperature difference (i.e. $p_{sat}=p_{\infty }=0$ and $T^{I}=1$) and (ii) evaporation owing to pressure difference ($p_{sat}=1$ and $T^{I}=T_{\infty }=0$). The distance between the two droplets and the value of the Knudsen number are varied in order to investigate the proximity and rarefaction effects. The results from three models are compared with the NSF, Grad-13 and CCR with $\alpha _0=3/5$.

Temperature-driven case: In figure 14, we show the mass flux and the heat flux per unit area as a function of Knudsen number for the temperature-driven case. The mass flux reduces monotonically with the Knudsen number whereas the NSF theory predicts a higher mass flux than the CCR. Figures 14(a) and 14(b) illustrate the mass flux with centre-to-centre distance $2l_C=0.001$ (i.e. when the droplets are next to each other) and $2l_C= 10$ (i.e. at a distance where proximity effects are negligible), respectively. The corresponding heat flux is shown in figures 14(c) and 14(d).

Figure 14.

Two interacting droplets for the temperature-driven case ($T^{I}=1, p_{sat}=0$): profiles of the (a) mass flux with centre-to-centre distance $2l_C=0.001$, (b) mass flux with $2l_C=10$ (c) heat flux with $2l_C=0.001$ and (d) heat flux with $2l_C=10$ as functions of the Knudsen number for NSF, Grad-13 and CCR with $\alpha _0=3/5$. Symbols denote MFS and curves represent analytic solution from a single droplet.

When the droplets are next to each other a shielding effect is observed, and, as a result, the mass flux reduces (cf. figure 14a,b). At $Kn=0.05$, all three models predict approximately $29\,\%$ reduction in the mass flux as compared with a single droplet (i.e. $l_C\rightarrow \infty$). The reduction in the mass flux decreases with increase in Knudsen number. For $Kn=1$, NSF, Grad-13 and CCR with $\alpha _0=3/5$ predicts approximately $9\,\%$, $13\,\%$ and $12\,\%$ reduction in mass flux, respectively.

The mass fluxes computed from the MFS (symbols) for centre-to-centre distance $10$ are shown in figure 14(b). The analytic results from a single droplet are also plotted (continuous lines) for comparison. As expected, as the droplets are moved farther apart, the results converge to the single-droplet case; the mass flux is within $5\,\%$ (for $Kn=0.05$) and $1\,\%$ (for $Kn=1$) for the single-droplet case.

Similarly, the proximity effects on the heat flux are visible in figures 14(c) and 14(d), for $2l_C=0.001$ and $2l_C=10$, respectively. Again, the heat flux reduces as $l_C\rightarrow 1/2$, and as the value for the Knudsen number increases the proximity effects diminish. The percentage difference in the heat flux from the single-droplet cases varies from approximately $29\,\%$ to $12\,\%$ for $2l_C= 1.001$ and from approximately $4\,\%$ to $1\,\%$ for $2l_C= 10$, as the value of Knudsen number varies from $0.05$ to 1.

The vaporisation dynamics of a pair of droplets placed adjacent to each other is markedly different from the single-droplet counterpart. In figure 15, (a) the streamlines superimposed on the velocity magnitude contours and (b) the heat-flux lines and the magnitude contours are shown. The heat leaves the hot droplets owing to negative temperature gradients while condensation occurs in order to compensate for the heat loss. Figure 15 shows that the mass flux and the heat flux at the confined region between the droplets are sufficiently lower than at the outer regions. In this region, the vapour temperature is relatively high, thus a smaller temperature jump which leads to a reduced mass flux and heat flux.

Figure 15.

Two droplets placed adjacent to each other ($2l_c=1.001$) for the temperature-driven case ($T^{I}=1, p_{sat}=0$): (a) streamlines and speed contours and (b) heat-flux lines and magnitude, obtained from the Grad-13 MFS. Vertical and horizontal axes represent $x$- and $y$-directions ($N_c = N_s = 650$ and $\gamma = R_s/R_c = 0.1$).

Pressure-driven case: In figure 16, we show (a) mass flux with centre-to-centre distance $2l_C=0.001$, (b) mass flux with $2l_C=10$, (c) heat flux with $2l_C=0.001$ and (d) heat flux with $2l_C=10$ as functions of the Knudsen number for the pressure-driven case. Again, in order to gain insights into the proximity effects, the numerical solution obtained via the MFS (symbols) and the analytic results from a single-droplet case (continuous lines) are compared in figure 16. In this case, the mass flux is slightly reduced (${\lesssim }3\,\%$ for $2l_C=0.001$ and ${\lesssim } 0.06\,\%$ for $2l_C=10$) in the range ($0.05\leq Kn\leq 1$). On the other hand, the proximity effects on the heat flux are observed to be significant. For $2l_C=0.001$, the magnitude of the heat flux is reduced approximately $29.3\,\%$ for NSF ($28.7\,\%$ for Grad-13 and CCR) at $Kn=0.05$ and approximately $10\,\%$ for NSF (${\leq } 15\,\%$ for Grad-13 and CCR) at $Kn=1$. For $2l_C=0.001$, all models give less than a $4\,\%$ reduction in heat flux, for all ranges of Knudsen number considered.

Figure 16.

Two interacting droplets for the pressure-driven case ($T^{I}=0, p_{sat}=1$): profiles of (a) mass flux with centre-to-centre distance $2l_C=0.001$, (b) mass flux with $2l_C=10$, (c) heat flux with $2l_C=0.001$ and (d) heat flux with $2l_C=10$ as functions of the Knudsen number for NSF, Grad-13 and CCR with $\alpha _0=3/5$. Symbols denote MFS and curves represent analytic solution from a single droplet.

Another interesting point to be noted from figure 16(c,d) is that the heat flux predicted by the NSF equations is significantly lower in magnitude compared with the CCR and Grad-13 models, especially in the transition regime ($Kn \geqslant 0.1$). Within this regime, the heat flux is not only driven by the temperature gradient (i.e.Fourier's law) but also owing to the pressure gradient, this second-order effect can easily be understood from the constitutive relations (2.5a,b) and the momentum balance equation (2.4a–c), which give

(6.11)\begin{equation} \boldsymbol{q}={-}\frac{c_{p}Kn}{Pr}\left( \boldsymbol{\nabla} T-\alpha _{0}\boldsymbol{\nabla} p \right). \end{equation}

Here, the second term on the right-hand side is the non-Fourier contribution to the heat flux owing to pressure gradient (Rana et al. Reference Rana, Gupta and Struchtrup2018a). In figure 17, we show two different contributions to the scaled normal heat flux ($\boldsymbol {q}\boldsymbol {\cdot }\boldsymbol {n}/Kn$) along the surface of the droplet (in the $x$–$z$ plane) for the pressure-driven case with $Kn=0.1$. As predicted from the Grad-13 theory, the heat flux owing to the pressure gradient (non-Fourier) and heat flux owing to the temperature gradient (Fourier) have opposing effects, and consequently the net heat flux is reduced in the second-order theories.

Figure 17.

Two interacting droplets for the pressure-driven case ($T^{I}=0, p_{sat}=1$): Fourier and non-Fourier contributions to the normal heat flux ($\boldsymbol {q}\boldsymbol {\cdot }\boldsymbol {n}/Kn$) given by the Grad-13 theory plotted along the top surface of the droplet in $x$–$z$ plane.

In figure 18, we again show the streamlines and the velocity magnitude contours (a), and the heat-flux lines and the magnitude contours (b) for the pressure-driven case. In this case, evaporation occurs owing to a negative pressure gradient and the heat flows inwards. The mass flux and the heat flux at the confined region between the droplets are reduced owing to a high vapour pressure.

Figure 18.

Two droplets placed next to each other ($2l_c=1.001$) for the pressure-driven case ($T^{I}=0, p_{sat}=1$): (a) stream lines and speed contours (b) heat-flux lines and magnitude, obtained from the Grad-13 MFS. Vertical and horizontal axis represent $x$- and $y$-directions. ($N_c = N_s = 650$ and $\gamma = R_s/R_c = 0.1$).

The high-pressure region created between two droplets pushes them away from each other. The net drag force (in the $x$-direction) acting on each droplet is plotted in figure 19 as a function $Kn$. The results for various values of centre-to-centre distance shown in figure 19 are computed from the Grad-13 model, the drag force predicted by the NSF model is within $10\,\%$ and not shown. The force decreases exponentially as droplets become farther apart, i.e. as $l_C$ increases. Interestingly, the drag force attains a minimum value at a certain Knudsen number, which depends on how far the droplets are from each other. As $Kn\rightarrow \infty$, the drag force vanishes.

Figure 19.

The drag force acting on the droplets versus Knudsen numbers for various values of the separation distance as predicted by the Grad-13 theory for the pressure-driven case.

6.4. Evaporation of a deformed droplet

In this section, we study evaporation/condensation from a deformed droplet with a fixed shape and size, in order to gain some initial insight into evaporation occurring during oscillatory and/or deformed droplet events, leaving the step to unsteady flow as a focus for future work. The effects of heat conduction and the flow inside the droplet are not considered. The parametric equation for the droplet surface in the second harmonic mode is given by (see, e.g., Sprittles & Shikhmurzaev Reference Sprittles and Shikhmurzaev2012)

(6.12)\begin{equation} \boldsymbol{X}=\eta _{0}\left[ 1+\frac{\eta }{2}\left( 3\cos ^{2}\phi -1\right) \right] \left\{ \sin \phi \cos \theta,\sin \phi \sin \theta,\cos \phi \right\}, \end{equation}

where $\phi \in \lbrack 0,{\rm \pi} ]$, $\theta \in \lbrack 0,2{\rm \pi} ]$ and $\eta$ is the shape factor. The volume of the droplet is

(6.13)\begin{equation} V=\eta _{0}\frac{4{\rm \pi} }{3}\left(\frac{35+21\eta ^{2}+2\eta ^{3}}{35}\right)^{1/3}{.} \end{equation}

Therefore, in order to have the volume $V$ equal to a spherical droplet of radius $1$, the normalising factor $\eta _{0}$ is given by

(6.14)\begin{equation} \eta _{0}=\left(\frac{35}{35+21\eta ^{2}+2\eta ^{3}}\right)^{1/3}{.} \end{equation}

The droplet shapes with different values of shape parameter $\eta =0$, $0.5$ and $1$ are depicted in figure 20; the corresponding (non-dimensional) surface area of the droplets are $4{\rm \pi} , 4{\rm \pi} \times 1.078$ and $4{\rm \pi} \times 1.192$, respectively.

Figure 20.

Droplet shapes given by second spherical harmonics with shape parameter $\eta = 0$, $0.5$ and $1$ and the normalising factor $\eta _0 = 35^{1/3}/(35+21 \eta ^2+2\eta ^3)^{1/3}$ ensuring that the droplet has fixed volume.

We shall again consider two cases: (i) flow caused by a unit pressure difference while the temperature of the liquid is equal to the far-field temperature (i.e. $p_{sat}-p_{\infty }=1$ and $T^{l}-T_{\infty }=0$) and (ii) flow caused by a unit temperature difference while the saturation pressure in the liquid is same as the far-field pressure (i.e. $p_{sat}-p_{\infty }=0$ and $T^{l} -T_{\infty }=1$).

Pressure-driven case: In figure 21(a–c), we show the mass flux per unit area as a function of $Kn$ computed using the NSF, Grad-13 and $\alpha _0=3/5$ models, respectively for the shape parameter $\eta = 0, 0.5$ and $1$. Here, the numerical results from the MFS are denoted by symbols and the analytic results for the spherical droplet ($\eta = 0$) are represented by the solid line. For small deformity ($\eta = 0, 0.5$), the mass flux is nearly same, however when the droplet is more deformed ($\eta = 1$) the mass flux is reduced. The percentage reduction in the mass flux (compared with a spherical droplet) is approximately $11\,\%$ for $Kn=0.05$ and approximately $21\,\%$ for $Kn=1$. In a similar manner to the case of two droplets placed next to each other, considered in the previous section, the reduction in the mass flux is caused by the formation of a high pressure in the high curvature (for $\eta = 1$) cusp region near $z=0$, which is not present at smaller $\eta$.

Figure 21.

Pressure-driven case ($p_{sat}-p_{\infty }=1$, $T^{l}-T_{\infty }=0$); mass flux per unit area versus Knudsen number with $\eta = 0$, $0.5$ and $1$: (a) NSF, (b) Grad-13 and (c) $\alpha _0=3/5$. The symbols denote the MFS and the solid lines represent analytic results for a spherical droplet.

Similarly, in figure 22(a–c) we plot the heat flux per unit area as a function of $Kn$ for the same shape parameters. Note that the heat flux owing to pressure difference is a first-order quantity in terms of $Kn$ (that is, it vanishes as $Kn\rightarrow 0$), hence in these figures we have scaled the heat flux with inverse $Kn$. Again, the heat-flux magnitude reduces as $\eta \rightarrow 1$. The reduction in the heat flux is approximately 13 %. Furthermore, the Grad-13 and CCR theories predict a larger heat flux compared with the NSF results owing to non-Fourier heat-flux contribution (6.11).

Figure 22.

Pressure-driven case ($p_{sat}-p_{\infty }=1$, $T^{l}-T_{\infty }=0$); heat flux per unit area/${Kn}$ versus Knudsen number with $\eta = 0$, $0.5$ and $1$: (a) NSF, (b) Grad-13 and (c) $\alpha _0=3/5$. The symbols denote the MFS and the solid lines represent analytic results for a spherical droplet.

Temperature-driven case: In figure 23, we present the heat flux per unit area for the temperature-driven case for different shape factors (see legend in the figures) for NSF (a), Grad-13 (b) and CCR with $\alpha _0=3/5$ (c). Note that, owing to the Onsager reciprocity relations (Chernyak & Margilevskiy Reference Chernyak and Margilevskiy1989), the mass flux for the temperature-driven case is the same as the heat flux for the pressure-driven case (figure 22), hence it is not shown. Evidently, from figure 23, the magnitude of the heat flux reduces with the deformation, i.e. as $\eta \rightarrow 1$. For NSF (figure 23a), the reduction is approximately 3 % for $\eta =0.5$ and approximately 12 % for $\eta =1$.

Figure 23.

Temperature-driven case ($p_{sat}-p_{\infty }=0$, $T^{l}-T_{\infty }=1$); heat flux per unit area/Kn versus Knudsen number with $\eta = 0$, 0.5 and 1: (a) NSF, (b) Grad-13 and (c) $\alpha _0=3/5$. The symbols denote the MSF solutions and the solid lines represent analytic results for a spherical droplet.

Curvature dependence of the saturation pressure: Throughout this paper, we have taken $T^{I}$ and $p_{sat}$ as independent parameters. However, in general, the Clausius–Clapeyron–Kelvin formula provides a relation between the saturation pressure and the temperature at the interface (linearised and dimensionless) (McElroy Reference McElroy1979; Bond & Struchtrup Reference Bond and Struchtrup2004):

(6.15)\begin{equation} p_{sat}\left( T^{I}\right) =p_{sat}^{p}\left( T^{I}\right) \underline{\exp \left( 2\gamma \kappa \right)} {.} \end{equation}

Here, $\gamma$ is the dimensionless surface tension ($\gamma =\hat {\gamma }/ \hat {\rho }_{l}\hat {\ell }\hat {R}\hat {T}_{0}$), $\kappa$ is the dimensionless curvature ($\kappa =\hat {\kappa }\hat {\ell }$) and $p_{sat}^{p}( T^{I}) =H_{0}T^{I}$ is the saturation pressure for a planar surface with $H_{0}$ (${=}\hat {H}_{0}/\hat {R}\hat {T}_{0}$) being the dimensionless heat of evaporation.

Thus far, the modelling assumption has been that $H_{0}$ (and $\gamma$) can be independently chosen to give a desired value of the saturation pressure. Notably, for the noble gases (vapour/liquid) the Kelvin correction term (underlined) in (6.15) are often small, unless the radius of curvature $\hat {\kappa }$ is in nanometers. For example, for argon, if we choose $\hat {T}_0=103\ \mathrm {K}$, $\hat {p}_0=p^p_{sat} (\hat {T}_0) = 0.4\ \mathrm {MPa}$, $\hat {\rho }_l = 1294.6\ \mathrm {kg} \ \mathrm {m}^{-3}$ and $\hat {\gamma } = 8.75\ \mathrm {mN}\ \mathrm {m}^{-1}$ (a case considered in Rana et al. Reference Rana, Lockerby and Sprittles2019), the Kelvin correction term gives $\exp (0.6305/a)$, where $a$ (in nanometres) is the inverse of local curvature $\hat {\kappa }$. Therefore, for a drop of size greater than $10\ \mathrm {nm}$, the effect of the Kelvin correction term will be less then $6.1\,\%$.

In order to compare our results with existing literature, in § 7 we use the dimensionless heat of evaporation $H_{0}$ as an independent parameter.