1. Introduction
Droplet vaporization is important in many applications. For example, vaporizing droplets are present in sprays that are used for purposes such as extinguishing fire and liquid fuel combustion. A spray will have a range of droplet sizes, with some droplets small enough to be in the Stokes flow regime. The behaviour of vaporizing droplets under microgravity conditions is also of fundamental interest as it allows for phenomena such as capillary flows to be investigated for large droplets that are observable readily but where the masking effects of buoyancy are negligible.
This research is focused upon investigating analytically the influence of thermal Marangoni stresses and liquid and gas properties, including the gas-phase Lewis number of the vaporizing component, on translational motion of a slowly vaporizing droplet in the Stokes flow regime. The droplet is single component, and the gas phase is composed of only two species – the ambient gas and the molecules that vaporize from the droplet surface.
There have been numerous investigations, spanning decades, that focus on the behaviour of vaporizing (and burning) droplets. Some recent reviews are available in Raghavan (Reference Raghavan2019), Zang et al. (Reference Zang, Tarafdar, Tarasevich, Choudhury and Dutta2019), Sazhin (Reference Sazhin2017), Karbalaei, Kumar & Cho (Reference Karbalaei, Kumar and Cho2016), Erbil (Reference Erbil2012) and Sazhin (Reference Sazhin2006). A monograph by Sirignano (Reference Sirignano2010) is also available. To the best knowledge of this author, there have been no studies that focused on how thermocapillary stresses and the Lewis number of the vaporizing species interact to influence the translational motion of a vaporizing droplet in slow viscous flow without a far-field temperature gradient. There are prior analytical studies of droplet vaporization (including combustion) and motion in Stokes flow, e.g. Gogos & Ayyaswamy (Reference Gogos and Ayyaswamy1988), Jog, Ayyaswamy & Cohen (Reference Jog, Ayyaswamy and Cohen1996) and Ackerman & Williams (Reference Ackerman and Williams2005), but these studies did not consider the influences of thermocapillary stresses. Subramanian, Zhang & Balasubramaniam (Reference Subramanian, Zhang and Balasubramaniam1999) analysed mass transport from a droplet executing thermocapillary motion, but where the thermocapillary motion was driven by a temperature gradient that was imposed in the far field. Computational modelling results that account for thermocapillary effects have also been reported, e.g. Niazmand et al. (Reference Niazmand, Shaw, Dwyer and Aharon1994) and Shih & Megaridis (Reference Shih and Megaridis1996), but these studies focused on droplet Reynolds numbers that were too large to correspond to Stokes flow, and Lewis number effects were not discussed.
We next discuss the partial differential equations for conservation of mass, species, momentum and energy in the gas and liquid phases.
2. Governing equations
2.1. Conservation of energy and species in the gas phase
A schematic of the geometry and coordinate system is shown in figure 1. The centre of the droplet is at the origin of the coordinates shown.
Figure 1. Schematic of the geometry.
Equation (2.1) shows the energy conservation equation for the gas phase:
\begin{equation} \alpha_o\,\nabla^{2} T_o = v_{o,r}\,\frac{\partial T_o}{\partial r} + v_{o,\theta}\, \frac{1}{r}\,\frac{\partial T_o}{\partial \theta}, \end{equation}
where 
$T$ is temperature, 
$v$ is velocity, 
$\alpha$ is thermal diffusivity, the subscript 
$o$ denotes the gas phase, and the subscripts 
$r$ and 
$\theta$ denote components in the 
$r$ and 
$\theta$ directions. The species conservation equation is
\begin{equation} D_{io}\,\nabla^{2} Y = v_{o,r}\,\frac{\partial Y}{\partial r} + v_{o,\theta}\, \frac{1}{r}\,\frac{\partial Y}{\partial \theta},\end{equation}
where 
$D_{io}$ is the binary diffusion coefficient between the vaporizing component and the ambient gas, and 
$Y$ is the mass fraction of this component. The subscript 
$o$ is not included on 
$Y$ because there are only two species in the gas phase.
2.2. Conservation of energy and species in the liquid phase
The energy conservation in the liquid phase is
\begin{equation} \alpha_i\,\nabla^{2} T_i = v_{i,r}\,\frac{\partial T_i}{\partial r} + v_{i,\theta}\, \frac{1}{r}\,\frac{\partial T_i}{\partial \theta}, \end{equation}
where the subscript 
$i$ denotes the liquid. Because we are considering only single-component liquids, the species conservation equation in the liquid is simply that the mass fraction is unity there.
2.3. Conservation of mass and momentum in the gas and liquid phases
We assume that the fluids inside and outside the droplet are Newtonian and constant-density, though each fluid can have its own unique density. We also assume that gravity is negligible, that steady-state conditions apply, that the flow is axisymmetric, and that the Reynolds number based on the droplet diameter is small relative to unity. Our starting point is then the stream function for Stokes flow in spherical coordinates. This stream function automatically satisfies mass conservation as well as the 
$r$ and 
$\theta$ components of the momentum equations (neglecting nonlinear terms, of course).
The stream functions 
$\psi$ for the outer and inner fluids, i.e. the gas and the liquid, are
$$\begin{gather} \psi_o = \left(\frac{A_1}{r} + A_2 r + A_3 r^2 + A_4 r^4 \right) \sin^2(\theta), \end{gather}$$
$$\begin{gather}\psi_i = \left(\frac{B_1}{r} + B_2 r + B_3 r^2 + B_4 r^4 \right) \sin^2(\theta), \end{gather}$$
where 
$A_1$, 
$A_2$, 
$A_3$, 
$A_4$, 
$B_1$, 
$B_2$, 
$B_3$, and 
$B_4$ are constants. Given these stream functions, we may calculate the velocity components in the 
$r$ and 
$\theta$ directions by applying the expressions
$$\begin{gather} v_r = \frac{1}{r^2 \sin(\theta)}\,\frac{\partial\psi}{\partial\theta}, \end{gather}$$
$$\begin{gather}v_\theta = \frac{-1}{r \sin(\theta)}\,\frac{\partial\psi}{\partial r}, \end{gather}$$
to each phase separately. We require that the velocity components at the droplet centre are finite, leading to 
$B_1 = B_2 = 0$. We also require that the velocity field becomes uniform as 
$r \rightarrow \infty$, yielding 
$A_3 = v_\infty /2$ and 
$A_4$ = 0, where 
$v_\infty$ is the velocity component in the 
$+x$ direction at 
$r = \infty$. The velocity components then become
$$\begin{gather} v_{o,r} = 2 \left(\frac{A_1}{r^3} + \frac{A_2}{r} + \frac{v_\infty}{2} \right) \cos(\theta), \end{gather}$$
$$\begin{gather}v_{o,\theta} = \left(\frac{A_1}{r^3} - \frac{A_2}{r} - v_\infty\right) \sin(\theta), \end{gather}$$
$$\begin{gather}v_{i,r} = 2 \left(B_3 + B_4 r^2 \right) \cos(\theta), \end{gather}$$
$$\begin{gather}v_{i,\theta} ={-}2\left(B_3 + 2 B_4 r^2 \right) \sin(\theta). \end{gather}$$The shear stress components in the inner and outer fluids may now be evaluated as
$$\begin{gather} \tau_{i,r \theta}=\mu_i \left(\frac{1}{r}\,\frac{\partial v_{i,r}}{\partial \theta} + \frac{\partial v_{i,\theta}}{\partial r}\right)={-}\frac{2 B_3 + 10 B_4 r^2}{r}\,\mu_i \sin(\theta), \end{gather}$$
$$\begin{gather}\tau_{o,r \theta} = \mu_o \left(\frac{1}{r}\,\frac{\partial v_{o,r}}{\partial \theta} + \frac{\partial v_{o,\theta}}{\partial r}\right) ={-}\frac{5 A_1 + A_2 r^2 + v_\infty r^3}{r^4}\,\mu_o \sin(\theta), \end{gather}$$
where 
$\mu$ is viscosity and 
$\tau$ is shear stress. The gas pressure field in the outer fluid, i.e. 
$p_o$, which is evaluated by integrating the momentum equation, is given by
\begin{equation} p_o= p_\infty + \frac{2 A_2}{r^2}\,\mu_o \cos(\theta), \end{equation}
where 
$p_\infty$ is the pressure far from the droplet. By integrating the gas pressure and shear stress components at the gas–liquid interface, we may solve for the resultant drag force 
$F_D$ acting on the droplet:
\begin{align} F_D &={-} 2{\rm \pi} R^2 \int_{0}^{\rm \pi} \left( \tau_{o,r\theta} |_{r = R} \sin(\theta) + p_o |_{r = R} \cos(\theta)\right) \sin(\theta)\,{\rm d}\theta \nonumber\\ &= \frac{8 {\rm \pi}R^3 v_\infty + 40 {\rm \pi}A_1}{3 R^3}\,\mu_o. \end{align} Equation (2.15) shows that 
$F_D$ depends only on 
$A_1$, i.e. the constants 
$A_2$, 
$B_3$ and 
$B_4$ are not present. To proceed, we evaluate 
$A_1$ by considering the interface conditions.
2.4. Interface conditions
We enforce the following conditions at every location on the interface.
(i) The radial velocity component in the gas is negligible, which is consistent with the assumption of slow evaporation.
(ii) The radial velocity component in the liquid is negligible.
(iii) The gas and liquid tangential velocity components are equal.
(iv) There is a balance between the liquid and gas viscous shear stresses and the surface tension gradient, which in turn is caused by a surface temperature gradient.
(v) The liquid and gas temperatures are equal.
(vi) The difference between the liquid and gas heat fluxes is balanced by the energy per unit time per unit area required for vaporization.
(vii) The liquid and the gas are in thermodynamic equilibrium.
These interface conditions are expressed mathematically as
$$\begin{gather} v_{o,r} |_{r = R} = 0, \end{gather}$$
$$\begin{gather}v_{i,r} |_{r = R} = 0, \end{gather}$$
$$\begin{gather}v_{o,\theta} |_{r = R} - v_{i,\theta} |_{r = R} = 0, \end{gather}$$
$$\begin{gather}\tau_{i,r \theta} |_{r = R} - \tau_{o,r \theta} |_{r = R} - S \sin(\theta) = 0, \end{gather}$$
$$\begin{gather}T_{o} |_{r = R} - T_{i} |_{r = R} = 0, \end{gather}$$
$$\begin{gather}\lambda_o \left.\frac{\partial T_{o}}{\partial r}\right|_{r = R} - \lambda_i \left.\frac{\partial T_{i}}{\partial r}\right|_{r = R} + \rho_o D_{io} L \left.\frac{\partial Y}{\partial r}\right|_{r = R} = 0, \end{gather}$$
$$\begin{gather}Y |_{r = R} - f(T_i) |_{r = R} = 0, \end{gather}$$where
\begin{equation} S =\frac{\sigma_T}{R \sin(\theta)} \left.\frac{\partial T}{\partial \theta}\right|_{r = R} \end{equation}and
\begin{equation} \sigma_T = \frac{\partial \sigma}{\partial T}. \end{equation}
Here, 
$\lambda$ and 
$\sigma$ are thermal conductivity and surface tension, respectively. The phase equilibrium expression 
$f(T_i) |_{r = R}$ is to be determined from a phase equilibrium relation.
Substituting (2.8), (2.9), (2.10) and (2.11) into (2.16), (2.17), (2.18) and (2.19) allows for the constants 
$A_1$, 
$A_2$, 
$B_3$ and 
$B_4$ to be evaluated:
$$\begin{gather} A_{1} = \frac{R^3 v_\infty (4 \mu_i - \mu_o) - 2 R^4 S}{8(2 \mu_i + \mu_o)}, \end{gather}$$
$$\begin{gather}A_{2} ={-} \frac{3 R v_\infty (4 \mu_i + \mu_o) - 2 R^2 S}{8(2 \mu_i + \mu_o)}, \end{gather}$$
$$\begin{gather}B_{3} ={-}\frac{3 v_\infty\mu_{o} + 2 R S}{8 {\left(2 \mu_{i} + \mu_{o}\right)}}, \end{gather}$$
$$\begin{gather}B_{4} = \frac{3 v_\infty\mu_{o} + 2 R S}{8 R^{2} {\left(2 \mu_{i} + \mu_{o}\right)}}. \end{gather}$$Substituting (2.25) into (2.15) yields the following expression for the drag force:
\begin{equation} F_{D} = F_\mu + F_\sigma, \end{equation}where
\begin{equation} F_{\mu} = \frac{{\rm \pi} R v_\infty \mu_o (12 \tilde{\mu} + 1)} {2 \tilde{\mu} + 1} \end{equation}is the droplet drag that would exist in the absence of thermocapillary effects, and
\begin{equation} F_{\sigma} = \frac{ - 10 {\rm \pi}R^2 S}{3 (2 \tilde{\mu} + 1)} \end{equation}
accounts for the presence of surface tension gradients. The symbol 
$\tilde {\mu }$ is the ratio of the dynamic viscosities of the inner and outer fluids, i.e. 
$\tilde {\mu }$ = 
$\mu _i/\mu _o$. Note that the surface tension gradient influences 
$F_\sigma$ via the quantity 
$S$. To proceed, we will now consider how 
$S$ is influenced by evaporation.
3. Solving the energy and species equations
3.1. The unperturbed solution
Solving the gas-phase conservation equations for a stagnant environment, i.e. 
$v_\infty = 0$, leads to the expressions
$$\begin{gather} T_o = T_\infty + (T_s - T_\infty)\,\frac{R}{r}, \end{gather}$$
$$\begin{gather}Y = Y_\infty + (Y_s - Y_\infty)\,\frac{R}{r}, \end{gather}$$
$$\begin{gather}\dot{m}_0 = 4 {\rm \pi}R\,\frac{\lambda_o}{L}\,(T_\infty - T_s) = 4 {\rm \pi}R \rho_{o} D_{io}\,\frac{Y_s - Y_\infty}{1 - Y_s}, \end{gather}$$
where the subscript 
$s$ denotes the droplet surface, and the subscript ‘0’ on 
$\dot {m}_0$ indicates a spherically symmetrical solution for the mass flow rate off the droplet surface (
$\dot {m}$). The enthalpy of vaporization is denoted by the symbol 
$L$. To proceed, we introduce the gas-phase Lewis number, 
$Le = \alpha _o/D_{io}$. It then follows that
\begin{equation} \frac{c_p (T_\infty - T_s) }{L} = \frac{1}{Le}\,\frac{Y_s - Y_\infty}{1-Y_s}, \end{equation}
and for 
$Y_s \ll 1$, which is consistent with the assumption of slow vaporization rates, we obtain
\begin{equation} \frac{Y_s - Y_\infty}{T_\infty-T_s}\,\frac{L}{c_p} = Le. \end{equation}
The variables 
$T_s$ and 
$Y_s$ are related via a Clausius–Clapeyron phase equilibrium relation
\begin{equation} Y_s = \varGamma \exp \left[{-\frac{L}{R_{g}} \left(\frac{1}{T_s} - \frac{1}{T_\infty}\right)}\right], \end{equation}
where 
$R_g$ is the gas constant for the vaporizing species, and 
$\varGamma$ is the equilibrium mass fraction of the vaporizing species at temperature 
$T_\infty$. Linearizing the phase equilibrium expression by considering small variations of 
$T_s$ about 
$T_\infty$ yields
\begin{equation} Y_s = \varGamma + \frac{\varGamma L} {R_g T_{\infty}^{2}}\, (T_s - T_\infty ) + \cdots.\end{equation}3.2. Influence of convection
We define the dimensionless liquid temperature 
$t = (T_i - T_\infty )/(T_s - T_\infty )$, the dimensionless gas temperature 
$h = (T_o - T_\infty )/(T_s - T_\infty )$, and the rescaled gas mass fraction 
$m = (Y - Y_\infty ) / (Y_s - Y_\infty )$, leading to the conservation equations
$$\begin{gather} \nabla^{2} t = \phi \left(u_{i,z}\,\frac{\partial t}{\partial z} + u_{i,\theta}\,\frac{1}{z}\,\frac{\partial t}{\partial \theta}\right), \end{gather}$$
$$\begin{gather}\nabla^{2} h = \epsilon \left(u_{o,z}\,\frac{\partial h}{\partial z} + u_{o,\theta}\,\frac{1}{z}\,\frac{\partial h}{\partial \theta}\right), \end{gather}$$
$$\begin{gather}\nabla^{2} m = \beta \left(u_{o,z}\,\frac{\partial m}{\partial z} + u_{o,\theta}\,\frac{1}{z}\,\frac{\partial m}{\partial \theta}\right), \end{gather}$$
where 
$z = r/R$, 
$u_{i,z} = v_{i,r}/v_\infty$, 
$u_{i,\theta } = v_{i,\theta }/v_\infty$, 
$u_{o,z} = v_{o,r}/v_\infty$ and 
$u_{o,\theta } = v_{o,\theta }/v_\infty$. In addition, 
$\phi = v_\infty R \mu _o /(\alpha _i (2 \mu _i + \mu _o))$ is a thermal transport Péclet number for the liquid, 
$\epsilon =v_\infty R /\alpha _o$ is a gas-phase thermal transport Péclet number, and 
$\beta = v_\infty R/D_{io}$ is a gas-phase species transport Péclet number. It is noted that 
$\beta = \epsilon \,Le$ and that a rescaled velocity within the liquid has been used to define 
$\phi$, where this velocity accounts for the difference in viscosity between the liquid and the gas.
For analysis, we assume 
$\phi \ll 1$, 
$\epsilon \ll 1$ and 
$\beta \ll 1$ along with the expansions
$$\begin{gather} t = t_0 + \phi t_1 + \cdots, \end{gather}$$
$$\begin{gather}h = h_0 + \epsilon h_1 + \cdots, \end{gather}$$
$$\begin{gather}m= m_0 + \beta m_1 + \cdots, \end{gather}$$
where 
$t_0 = 1$, 
$h_0 = 1/z$ and 
$m_0 = 1/z$ are the unperturbed (spherically symmetrical) solutions, and 
$t_1$, 
$h_1$ and 
$m_1$ correct the spherically symmetrical solutions when convection is present. Substituting these expansions into the transport equations for 
$t$, 
$h$ and 
$m$ yields the following first-order problem:
$$\begin{gather} \nabla^{2} t_1 = 0, \end{gather}$$
$$\begin{gather}\nabla^{2} h_1 ={-} \frac{1}{z^2} {\left(1 - \frac{3 - \zeta}{2z} + \frac{1 - \zeta}{2z^{3}}\right)} \cos(\theta), \end{gather}$$
$$\begin{gather}\nabla^{2} m_1 ={-} \frac{1}{z^2} {\left(1 - \frac{3 - \zeta}{2z} + \frac{1- \zeta}{2z^{3}}\right)} \cos(\theta), \end{gather}$$
where 
$u_{o,z}$ has been evaluated by setting 
$S = 0$. It is also noted that 
$\zeta = 3 \tilde {\mu }/(4 + \tilde {\mu })$.
The general solutions of (3.14), (3.15) and (3.16), i.e. the sums of the homogeneous and inhomogeneous solutions, are
$$\begin{gather} t_1=\sum_{n=0}^{\infty} \left(a_{n} z^n + \frac{b_{n}}{z^{n+1}}\right)P_{n}, \end{gather}$$
$$\begin{gather}h_{1} = {\sum_{n=0}^{\infty} \left(c_{n} z^n + \frac{d_{n}}{z^{n+1}}\right)P_{n}} + \left(\frac{1}{2} + \frac{\zeta-3}{4z} - \frac{3}{8}\,\frac{\zeta-1}{z^2} + \frac{\zeta-1}{8 z^3}\right) \cos(\theta), \end{gather}$$
$$\begin{gather}m_{1} = \sum_{n=0}^{\infty} {\left(f_{n} z^n + \frac{g_{n}}{z^{n+1}}\right)P_{n}} + \left(\frac{1}{2} + \frac{\zeta-3}{4z} - \frac{3}{8}\,\frac{\zeta-1}{z^2} + \frac{\zeta-1}{8 z^3}\right) \cos(\theta), \end{gather}$$
where 
$a_n, b_n, c_n, d_n, f_n,g_n$ are arbitrary constants. We restrict our attention to situations where the temperature along the droplet surface varies as a constant added to another constant multiplied by 
$\cos (\theta )$, leading to
$$\begin{gather} t_1= a_0 + a_1 z \cos(\theta), \end{gather}$$
$$\begin{gather}h_{1}= c_0 + \frac{d_0}{z} + \frac{d_1}{z^2} \cos(\theta) +\left(\frac{1}{2} + \frac{\zeta-3}{4z} - \frac{3}{8}\, \frac{\zeta-1}{z^2} + \frac{\zeta-1}{8 z^3}\right) \cos(\theta), \end{gather}$$
$$\begin{gather}m_{1}= f_0 + \frac{g_0}{z} + \frac{g_1}{z^2} \cos(\theta) +\left(\frac{1}{2} + \frac{\zeta-3}{4z} - \frac{3}{8}\, \frac{\zeta-1}{z^2} + \frac{\zeta-1}{8 z^3}\right) \cos(\theta). \end{gather}$$
There are now a total of eight undetermined coefficients, i.e. 
$a_0, a_1, c_0, d_0, d_1, f_0, g_0,g_1$. The remaining interface conditions, i.e. (2.20), (2.21) and (2.22), allow six of these coefficients to be determined. There are no values of 
$c_0$ and 
$f_0$, however, that enable 
$h_1$ and 
$m_1$ to vanish as 
$z \rightarrow \infty$. As such, the solutions for 
$h_1$ and 
$m_1$ must be regarded as solutions that apply in an inner zone. The coefficients 
$c_0$ and 
$f_0$ are determined by matching to outer-zone solutions.
For analysis of the outer-zone energy equation, we define the variables 
$Z = \epsilon z$ and 
$H = h/\epsilon$, leading to
\begin{equation} \nabla^{2} H_0 = \frac{\partial H_0}{\partial Z} \cos(\theta) - \frac{1}{Z}\,\frac{\partial H_o}{\partial \theta} \sin(\theta), \end{equation}
where we have employed the approximations 
$v_{or} \approx v_\infty \cos (\theta )$ and 
$v_{o\theta } \approx - v_\infty \sin (\theta )$ for large 
$r$, and the subscript ‘0’ indicates that 
$H_0$ is the leading-order solution in the outer zone. Equation (3.23) has the solution, as discussed by Acrivos & Taylor (Reference Acrivos and Taylor1962),
\begin{align} H_0 &= C_0\,\frac{\rm \pi}{Z} \exp\left({\frac{Z \left[\cos(\theta)-1\right]}{2}}\right) P_0 +C_1\,\frac{\rm \pi}{Z} \left(1 + \frac{2}{Z}\right) \exp\left({\frac{Z \left[\cos(\theta)-1\right]}{2}}\right) P_1 \nonumber\\ &\quad +C_2\,\frac{\rm \pi}{Z} \left(1 + \frac{6}{Z} + \frac{12}{Z^2}\right) \exp\left({\frac{Z \left[\cos(\theta)-1\right]}{2}}\right) P_2+\cdots, \end{align}
where 
$C_0,C_1,C_2,\ldots$ are arbitrary constants that are evaluated by matching, which yields 
$c_0 = -1/2$, 
$C_o = 1/{\rm \pi}$ and 
$C_n = 0$ for 
$n \geq 1$. The leading-order solution in the outer zone is then
\begin{equation} H_0=\frac{1}{\epsilon z} \exp\left({\frac{\epsilon z \left[\cos(\theta)-1\right]}{2}}\right). \end{equation}
Equation (3.25) satisfies the boundary condition 
$h(z = \infty,\theta ) = 0$.
Analysis of the species equation follows a similar path. The leading-order solution in the outer zone is found to be given by
\begin{equation} M_0=\frac{1}{\beta z} \exp\left({\frac{ \beta z \left[\cos(\theta)-1\right]}{2}}\right), \end{equation}
where 
$M = m/\beta$. We also note that matching yields 
$f_0 = -1/2$.
Applying the interface conditions enables the remaining coefficients to be evaluated:
$$\begin{gather} a_0=\frac{1}{2}\,\frac{\epsilon}{\phi}\,\frac{Le(Le-1)}{Le + \eta}, \end{gather}$$
$$\begin{gather}a_1 = \frac{1}{16}\,\frac{\phi}{\epsilon}\,\frac{(\zeta+3)(1-Le)}{1 + \tilde{\lambda}/2 + {\varGamma \eta}/{Le}}, \end{gather}$$
$$\begin{gather}d_0 = \frac{1}{2}\,\frac{Le^2 + \eta}{Le + \eta}, \end{gather}$$
$$\begin{gather}d_1 = \frac{1}{16}\,\frac{(\zeta+3)(1-Le)}{1 + \tilde{\lambda}/2 + {\varGamma \eta}/{Le}}, \end{gather}$$
$$\begin{gather}g_0 = \frac{1}{2\,Le}\,\frac{Le^2 + \eta}{Le + \eta}, \end{gather}$$
$$\begin{gather}g_1 ={-} \frac{1}{16}\,\frac{\eta}{Le^2}\,\frac{(\zeta+3)(1-Le)}{1 + \tilde{\lambda}/2 + {\varGamma \eta}/{Le}}, \end{gather}$$where
\begin{equation} \eta=\frac{ L^2}{c_p R_g T_{\infty}^2}. \end{equation}
The quantity 
$S$ is evaluated by noting that the surface temperature gradient is given by 
$\epsilon (T_\infty - T_s)~a_1\sin (\theta )$. Substituting this gradient into (2.23) eventually yields
\begin{equation} S=\frac{(\zeta + 3) (1 - Le)}{16(1 + \tilde{\lambda}/2 + \varGamma \eta /Le)}\, \frac{\sigma_T (T_\infty - T_s) v_\infty}{\alpha_o}, \end{equation}
where 
$\tilde {\lambda } = \lambda _i/\lambda _o$, and the temperature difference 
$T_\infty - T_s$ is evaluated using
\begin{equation} T_\infty - T_s=\frac{L}{c_p}\,\frac{\varGamma - Y_\infty}{Le + \varGamma \eta}. \end{equation}The tangential velocity at the interface is given by
\begin{equation} v_{i,\theta} |_{r = R} ={-} \frac{3 v_\infty}{4 (1 + 2 \tilde{\mu})} \left(1 - \frac{R}{R_1}\right) \sin(\theta), \end{equation}where
\begin{equation} R_1 ={-} \frac{4 (4 + \tilde{\mu}) (1 + \tilde{\lambda}/2 + \varGamma \eta /Le) \mu_o \alpha_o}{(2 + \tilde{\mu}) \sigma_T (T_\infty - T_s)(1-Le)} \end{equation}
is the droplet radius at which the tangential interface velocity equals zero, but where 
$v_\infty \neq 0$. Because 
$\sigma _T < 0$ for common liquids, 
$R_1 > 0$ if 
$Le < 1$, and 
$R_1 < 0$ if 
$Le > 1$. The latter situation corresponds to cases where the tangential velocity is never zero (again, where 
$v_\infty \neq 0$).
We also have the following expression for 
$F_\sigma$:
\begin{equation} F_{\sigma}=\frac{5 {\rm \pi}}{4}\,\frac{2 + \tilde{\mu}}{(1 + 2 \tilde{\mu})(4 + \tilde{\mu})}\, \frac{(Le-1) \sigma_T (T_\infty - T_s)}{1 + \tilde{\lambda}/2 + \varGamma \eta / Le}\, \frac{v_\infty R^2}{\alpha_o}. \end{equation}
It is evident that if 
$\sigma _T < 0$, then 
$F_\sigma > 0$ (thermocapillary effects enhance drag) if 
$Le > 1$, and 
$F_\sigma < 0$ (thermocapillary effects decrease drag) if 
$Le < 1$. This behaviour is related to the temperature gradient along the droplet surface. Two quantities that determine the interface temperature profile are 
$d_0$ and 
$d_1$, where 
$d_0$ is related to the average interface temperature, and 
$d_1$ is related to the interface temperature gradient. The constant 
$d_0$ is always positive, so convection increases the average droplet surface temperature. However, 
$d_1$ is negative if 
$Le > 1$ and positive if 
$Le < 1$, such that the sign of the surface temperature gradient depends upon the Lewis number. If 
$Le < 1$, then the resulting thermocapillary stresses reduce the interfacial velocity, thus increasing the droplet drag. Conversely, if 
$Le > 1$, then the resulting thermocapillary stresses increase the interfacial velocity, which decreases the droplet drag. The reversal of the surface temperature gradient is discussed in Appendix A, where it is shown that this reversal is related to the effects of convection on the surface energy and mass balances.
We also define 
$\tilde {F} = F_D / F_\mu = 1 + F_\sigma / F_\mu$, which is independent of 
$v_\infty$, as shown:
\begin{equation} \tilde{F}=1 + \frac{5(2 + \tilde{\mu})}{4(4 + \tilde{\mu})(1 + 12 \tilde{\mu})}\, \frac{(Le-1) \sigma_T (T_\infty - T_s)}{1 + \tilde{\lambda}/2 + \varGamma \eta / Le}\, \frac{R}{\mu_o \alpha_o}. \end{equation} A critical situation arises when 
$\tilde {F} = 0$. Equation (3.40) shows the relationship for the radius 
$R_2$ that corresponds to 
$\tilde {F} = 0$:
\begin{equation} R_2=\frac{4 (4 + \tilde{\mu})(1 + 12 \tilde{\mu})}{5(2 + \tilde{\mu})}\, \frac{1 + \tilde{\lambda}/2 + \varGamma \eta / Le}{(1-Le) \sigma_T (T_\infty - T_s)}\, \mu_o \alpha_o = \frac{1 + 12 \tilde{\mu}}{5}\,R_1. \end{equation}Inserting (3.40) into (3.39) yields
\begin{equation} \tilde{F}=1 - \frac{R}{R_2}. \end{equation}
We interpret these results within the context of a droplet in a quiescent environment where gravity and other body forces are absent. Suppose that at time 
$t = 0$, a droplet with radius 
$R$ is slowly moving with speed 
$v_{\infty }$. If 
$R_2 < 0$, then 
$\tilde {F} > 0$ and the velocity of the droplet will decay continually with time. If 
$R_2 > 0$ and 
$R < R_2$, then 
$\tilde {F} > 0$ and the velocity of the droplet will also decay continually with time. If, however, 
$R = R_2$, then 
$\tilde {F} = 0$ and the droplet will move with a steady translational motion.
An interesting situation occurs when 
$R_2 > 0$ and 
$R > R_2$, in which case 
$\tilde {F} < 0$. Thermocapillary effects are then sufficiently strong to cause the droplet velocity to increase with time, where the droplet will accelerate in a direction that is opposite to the oncoming free stream flow. The creeping-flow assumption used in this analysis does not predict the existence of a velocity where 
$\tilde {F} = 0$ when 
$R > R_2$, so it appears to be possible for the droplet velocity to increase enough such that the assumption of Stokes flow no longer applies, at which point a velocity would likely be attained where 
$\tilde {F} = 0$. Further analysis is needed to investigate this situation. Eventually, however, the droplet velocity would decrease because of the reduction in 
$R$ as the droplet vaporizes. These analyses show that in the Stokes flow regime, the radius 
$R_2$ is a bifurcation point for a droplet in a microgravity environment. That is, if 
$R_2 > 0$, then the droplet velocity: (i) decreases when 
$R < R_2$; (ii) remains constant when 
$R = R_2$; and (iii) increases when 
$R > R_2$.
Table 1 shows representative values for 
$R_1$ and 
$R_2$ for different liquids and ambient gases with 
$p_\infty = 1$ atm (abs), 
$Y_\infty = 0$ and 
$T_\infty = 300$ K. The ambient gases are selected to allow for different Lewis numbers to be attained. If the Lewis number is less than unity, then 
$R_1$ and 
$R_2$ are both negative, while if the Lewis number is greater than unity, then 
$R_1$ and 
$R_2$ are both positive. For example, a CH
$_3$OH droplet in an N
$_2$ or He environment has 
$R_2 > 0$, but in an Xe environment, the droplet has 
$R_2 < 0$. The gas and liquid properties were evaluated using data from Linstrom & Mallard (Reference Linstrom and Mallard2021). Lewis numbers were obtained from Dandy (Reference Dandy2021).
Table 1. Values of 
$R_1$ and 
$R_2$ for 
$p = 1$ atm (abs), 
$Y_\infty =0$ and 
$T_\infty = 300$ K.
4. Convective modification of the vaporization rate
The mass flow rate off the droplet is calculated by integrating the mass flux over the droplet surface (4.1):
\begin{equation} \dot{m}={-} 2 {\rm \pi}R^2 \rho_o D_{io} \int_{0}^{\rm \pi} \left.\frac{\partial Y}{\partial r}\right|_{r = R} \sin(\theta)\,{\rm d}\theta. \end{equation}Evaluating the integral yields
\begin{equation} \frac{\dot{m}}{\dot{m}_0}=1 + \beta\,\frac{1}{2}\, \frac{Le^2 + \eta}{Le + \eta} + \cdots,\end{equation}where
\begin{equation} \dot{m}_0 = 4 {\rm \pi}R\,\frac{\lambda_o}{L}\,(T_\infty - T_s). \end{equation} The evaporation constant 
$K$, which is the negative of the time-rate-of-change of the square of the droplet diameter, is defined as
\begin{equation} K={-} 4\,\frac{{\rm d} R^2}{{\rm d} t} = \frac{2 \dot{m}}{{\rm \pi} R \rho_i}, \end{equation}
where 
$\rho _i$ is the droplet density. Inserting (4.4) and (4.6) into (4.2) yields
\begin{equation} \frac{K}{K_0}=1 + \beta\,\frac{1}{2}\,\frac{Le^2 + \eta}{Le + \eta} + \cdots,\end{equation}where
\begin{equation} K_0 = \frac{2 \dot{m}_0}{{\rm \pi} R \rho_i}. \end{equation}
Equations (4.2) and (4.5) illustrate the influence of the Lewis number on vaporization rates, i.e. increasing 
$Le$ increases the vaporization rate if 
$\beta$ is held constant. Similarly, if 
$Le$ is held constant, then increasing 
$\beta$, e.g. by increasing 
$R$ or 
$v_\infty$, increases the vaporization rate.
5. Droplet motion considerations
We now consider the rectilinear motion of a droplet in a quiescent environment and in the absence of any body forces, e.g. gravity. Equations (5.1) and (5.2) show the differential equations for the position 
$x$ and the velocity 
$v_\infty$ of a droplet, respectively, where 
$m$ is the instantaneous droplet mass, and 
$F_D$ is evaluated with (2.29):
$$\begin{gather} \frac{{\rm d}\kern0.7pt x}{{\rm d} t}= v_\infty, \end{gather}$$
$$\begin{gather}\frac{{\rm d} v_\infty}{{\rm d} t} ={-} \frac{F_D}{m}. \end{gather}$$
Note that here, 
$v_\infty$ is the droplet velocity in the reference frame of the gas.
Because the droplet is evaporating slowly, we neglect momentum transfer that may exist from non-uniform vaporization along the droplet surface. The droplet mass history is evaluated using
\begin{equation} \frac{{\rm d} m}{{\rm d} t}={-} \dot{m}, \end{equation}
where 
$\dot {m}$, which is the mass flow rate off the droplet surface, is evaluated using (4.2) and (4.3). Equations (5.1), (5.2) and (5.3) are integrated numerically using an available ordinary differential equation solver in the language R (Soetaert, Petzoldt & Setzer Reference Soetaert, Petzoldt and Setzer2010).
Figure 2 shows results for a C
$_8$H
$_{18}$ (n-octane) droplet in a quiescent N
$_2$ environment at 1 atm (abs) and 300 K. The initial droplet radius is 
$R = 0.53$ mm, which is larger than 
$R_2 = 0.48$ mm, and the initial droplet velocity is 0.1 mm s
$^{-1}$. The droplet radius decreases with time because of vaporization. The droplet velocity initially increases with time until the droplet radius decreases to 0.48 mm, at which time the droplet velocity is about 0.4 mm s
$^{-1}$, i.e. after about 50 s, after which the droplet velocity decreases with time. The total distance that the droplet moves is nearly 30 mm over a time scale of about 100 s. Also shown in figure 2 are histories for 
$Re$, 
$\epsilon$, 
$\beta$ and 
$\phi$. These quantities are all much smaller than unity over the calculation period. In addition, they all display peaks at a time slightly larger than 50 s, which is when 
$R$ becomes less than 
$R_2$.
Figure 2. Results for an n-octane droplet in a nitrogen environment with 
$R > R_2$.
For comparison, figure 3 shows results for another C
$_8$H
$_{18}$ (n-octane) droplet in a quiescent N
$_2$ environment at 1 atm (abs) and 300 K. Now the initial droplet radius is 
$R = 0.43$ mm, however, which is smaller than 
$R_2 = 0.48$ mm. The initial droplet velocity is 0.1 mm s
$^{-1}$. As before, the droplet radius decreases with time because of vaporization. Unlike figure 2, however, the droplet velocity decreases monotonically with time. The total distance that the droplet moves is slightly over 1 mm over a time scale of about 50 s. The histories for 
$Re$, 
$\epsilon$, 
$\beta$ and 
$\phi$ are also shown. These quantities are all much smaller than unity over the calculation period, and they all decrease monotonically with time. The velocity profile in figure 3 is similar to the velocity profile in figure 2 once the velocity has begun to decay.
Figure 3. Results for an n-octane droplet in a nitrogen environment with 
$R < R_2$.
We now consider the terminal velocity of a droplet as it falls vertically through a quiescent gas under the influence of gravity. We assume quasi-steady conditions, i.e. the weight of the droplet is balanced by buoyancy and drag forces. Under these conditions, we may derive
\begin{equation} \tilde{v}=\frac{v_\infty}{v_{\infty_{ref}}}=\frac{1}{1 + R/R_2}, \end{equation}where
\begin{equation} v_{\infty_{ref}}=\frac{4 R^2 (\rho_i - \rho_o) g}{3 \mu_o}\,\frac{1 + 2 \tilde{\mu}}{1 + 12 \tilde{\mu}} \end{equation}
is the terminal velocity in the absence of thermocapillary effects, and 
$g$ is the acceleration due to gravity. Equation (5.4) indicates that when 
$R_2=0$ (
$Le = 1$), a droplet will fall with speed 
$v_{\infty _{ref}}$. However, if 
$R_2<0$ (
$Le < 1$), then 
$v_{\infty } > v_{\infty _{ref}}$, while if 
$R_2 > 0$ (
$Le > 1$), then 
$v_{\infty } < v_{\infty _{ref}}$.
As an example, consider a water droplet of radius 0.01 mm falling under the influence of normal gravity (
$g = 9.81$ m s
$^{-1}$) in an N
$_2$ environment with the environmental conditions of table 1. Calculations yield 
$\tilde {v} = 1.000$, where the Reynolds number based on the droplet diameter is 
$Re = 0.014$. If the environment is changed to He, then it is found that 
$\tilde {v} = 0.998$ and 
$Re = 0.0018$, while for an Xe environment, 
$\tilde {v} = 1.019$ with 
$Re = 0.044$. These Reynolds numbers are all small relative to unity, so the assumption of Stokes flow is valid.
6. Conclusion
The influence of thermal Marangoni stresses and gas-phase properties on the motion of a slowly evaporating single-component droplet in Stokes flow has been investigated. The conservation equations were solved in the liquid and gas phases, and coupled at the gas–liquid interface by applying interface conditions for conservation of mass, species, momentum and energy. It was found that thermal Marangoni stresses influence droplet drag by changing the flow field in the vicinity of a droplet, and that the gas-phase Lewis number of the evaporating component determines whether Marangoni effects increase or decrease droplet drag. If the Lewis number is less than unity, then thermal Marangoni effects increase droplet drag, while if the Lewis number is greater than unity, then thermal Marangoni effects decrease droplet drag. This behaviour is related to the sign of the temperature gradient along the droplet surface that is induced by convection.
Conditions exist where a vaporizing droplet in a microgravity environment may exhibit constant or accelerating translational motion that is driven by thermocapillary effects. This behaviour will occur only if the droplet radius 
$R$ is larger than or equal to a critical radius 
$R_2$. If the droplet radius is smaller than 
$R_2$, then the droplet speed will decay with time.
The quasi-steady speed that a falling droplet would attain was also considered. It was found that relative to 
$Le = 1$, if 
$Le < 1$, then a droplet falls faster, while if 
$Le > 1$, then a droplet falls more slowly. This behaviour is caused by thermocapillary flow induced by an interfacial temperature gradient that develops from the droplet–gas relative motion.
Funding
This research received no specific grant from any funding agency, commercial or not-for-profit sectors.
Declaration of interests
The author reports no conflict of interest.
Appendix A
To expand further on the influence of the Lewis number on the surface temperature gradient, we express the outer-zone solutions in terms of 
$h$ and 
$m$ as
$$\begin{gather} h_{outer zone}\approx \frac{1}{z} \exp\left({\frac{\epsilon z \left[\cos(\theta)-1\right]}{2}}\right), \end{gather}$$
$$\begin{gather}m_{outer zone} \approx \frac{1}{z} \exp\left({\frac{\epsilon Le z \left[\cos(\theta)-1\right]}{2}}\right). \end{gather}$$
Equations (A1) and (A2) are plotted in figure 4 for 
$\epsilon = 0.1$ and Lewis numbers 0.8 and 1.2, where the solid lines correspond to (A1), and the dashed lines correspond to (A2). When 
$Le = 0.8$, the contours for 
$h$ are inside the contours of 
$m$, while for 
$Le = 1.2$, the opposite situation holds. When 
$Le = 1$, the contours overlap completely such that energy and mass are transported in exactly the same way (this situation is not shown in the figure). The solid and dashed lines in figure 4 overlap for 
$\theta = 0$, but they do not overlap for other values of 
$\theta$. The largest gap between the solid and dashed lines occurs when 
$\theta = {\rm \pi}$ such that Lewis number effects are strongest upstream of the droplet (for 
$Le \ne 1$). Finally, it is noted that convection increases gradients on the upstream side of the droplet, while gradients are decreased on the downstream side.
Figure 4 suggests an approximate approach to increase understanding of the reversal of the surface temperature gradient. To proceed, we divide the droplet into upstream and downstream halves, denoted by ‘
$+$’ and ‘
$-$’, respectively. Energy conservation is enforced on each half separately, and for simplicity we assume that energy transport through the droplet (between the two halves) is negligible.
An energy balance at the upstream surface of the droplet, i.e. 
$r = R$, yields
\begin{equation} \lambda_+\,\frac{{\rm d} T_+}{{\rm d} r} ={-} \rho_o D_+ L\,\frac{{\rm d} Y_+}{{\rm d} r}, \end{equation}
where 
$\lambda _+$ and 
$D_+$ are effective transport properties that account for both convection and diffusion. We characterize the derivatives as 
${\rm d} T_+/{\rm d} r = (T_\infty - T_{s+})/R$ and 
${\rm d} Y_+/{\rm d} r = (Y_\infty - Y_{s+})/R$, enabling the surface energy balance to be expressed as
\begin{equation} \frac{\lambda_+}{\rho_o D_+}\,(T_\infty - T_{s+}) = Y_{s+} - Y_\infty. \end{equation}Use of (3.7) allows the interface energy balance to be written as
\begin{equation} T_{s+} = T_\infty + \frac{\varGamma - Y_\infty}{\dfrac{\lambda_+}{\rho_o D_+} + \dfrac{\varGamma L}{R_g T_\infty^2}}. \end{equation}Applying this same procedure to the downstream half of the droplet yields
\begin{equation} T_{s-} = T_\infty + \frac{\varGamma - Y_\infty}{\dfrac{\lambda_-}{\rho_o D_-} + \dfrac{\varGamma L}{R_g T_\infty^2}}, \end{equation}
where 
$\lambda _- > \lambda _o$ and 
$D_- > D_o$ for 
$\epsilon >0$.
To evaluate 
$\lambda _+$, we characterize the energy flux to the upstream droplet surface as the sum 
$\lambda _o (T_\infty - T_{s+})/R + v_\infty c_p (T_\infty - T_{s+})$, i.e. diffusion plus convection, which can be written as 
$\lambda _o (1 + \epsilon ) (T_\infty - T_{s+})/R$. We thus define 
$\lambda _+ = \lambda _o (1 + \epsilon )$.
The mass flux from the upstream droplet surface to the environment is characterized as a sum of diffusive and convective components, i.e. 
$\rho _o D_{io} (Y_{s+} - Y_\infty )/R + \rho _o v_\infty (Y_{s+} - Y_\infty )$, which can be expressed as 
$\rho _o D_{io} (1 + \epsilon \,Le) (Y_{s+} - Y_\infty )/R$. This yields 
$D_+ = D_{io} (1 + \epsilon \,Le)$. Similar analyses for the downstream surface yield 
$\lambda _- = \lambda _o (1 - \epsilon )$, and 
$D_- = D_{io} (1 - \epsilon \,Le)$, where it is noted that the convective flow is away from the droplet in this case.
Substituting these variables into the interface energy balances and then subtracting 
$T_{s-}$ from 
$T_{s+}$ yields
\begin{equation} T_{s+} - T_{s-} = \frac{\varGamma - Y_\infty}{\dfrac{\lambda_o}{\rho_o D_o}\, \dfrac{1 + \epsilon}{1 + \epsilon Le} + \dfrac{\varGamma L}{R_g T_\infty^2}} - \frac{\varGamma - Y_\infty}{\dfrac{\lambda_o}{\rho_o D_{io}}\, \dfrac{1-\epsilon}{1-\epsilon Le} + \dfrac{\varGamma L}{R_g T_\infty^2}}. \end{equation}
We now expand (A7) in a Taylor series for 
$\epsilon \ll 1$, yielding
\begin{equation} T_{s+} - T_{s-} = C \epsilon (Le-1)+\cdots,\end{equation}
where 
$C$ is a constant. Equation (A8) yields the result that to leading order in 
$\epsilon$, the surface temperature gradient, which is characterized as 
$(T_{s+} - T_{s-})/{\rm \pi}$, changes its sign if 
$Le$ becomes smaller or larger than unity (for 
$\epsilon \neq 0$). If 
$Le = 1$, then to leading order the surface temperature gradient is zero.
Open access
