3.1. Asymptotic expansion solution

Consistent with Wray & Moore (Reference Wray and Moore2023), we select the amplitude of the deformation $ \varepsilon $ as the perturbation parameter, enabling us to express the vapour concentration as an asymptotic series

(3.1) \begin{equation} c\left(r,\theta ,\varepsilon \right)\sim c^{0}\left(r,\theta \right)+\varepsilon\, c^{1}\left(r,\theta \right), \end{equation}

and give the evaporation flux by

(3.2) \begin{equation} J\left(\theta ,\varepsilon \right)\sim J^{0}\left(\theta \right)+\varepsilon\, J^{1}\left(\theta \right)=\frac{\partial c^{0}\left(r,\theta \right)}{\partial \boldsymbol{n}_{{i}}}\Bigg|_{r={r_{s}}}+\varepsilon\, \frac{\partial c^{1}\left(r,\theta \right)}{\partial \boldsymbol{n}_{{i}}}\Bigg|_{r={r_{s}}} .\end{equation}

Substituting (3.1) into (2.2), we get

(3.3) \begin{equation} {\Delta} c^{0}+\varepsilon\, {\Delta} c^{1}=0 .\end{equation}

Since (3.3) holds for arbitrary values of $ \varepsilon $ , we have ${\Delta} c^{0}=0$ and ${\Delta} c^{1}=0$ . In other words, both c ⁰ and c ¹ satisfy Laplace’s equation separately. Substituting (3.1) into (2.5) produces

(3.4a,b) \begin{equation} c^{0}(r\rightarrow \infty ,\theta )=0,\quad c^{1}(r\rightarrow \infty ,\theta )=0, \end{equation}

(3.5) \begin{equation} c^{0}\left(r=r_{s},\theta \right)+\varepsilon\, c^{1}\left(r=r_{s},\theta \right)=1 .\end{equation}

Given that r _s is a function of $ \varepsilon $ as well, the explicit dependence on $ \varepsilon $ in (3.5) has not been presented yet. This is attributed to the fact that $c^{0}(r=r_{s},\theta )$ and $c^{1}(r=r_{s},\theta )$ rely on $ \varepsilon $ via the variable r. To address this issue, we expand each term in a Taylor series to first order in $ \varepsilon $ around r = 1. Consequently, we can express (3.5) as

(3.6) \begin{equation} c^{0}\left(1,\theta \right)+\varepsilon\, c_{r}^{0}\left(1,\theta \right)P_{n}\left(\cos \theta \right)+\varepsilon\, c^{1}\left(1,\theta \right)+\varepsilon ^{2}\,c_{r}^{1}\left(1,\theta \right)P_{n}\left(\cos \theta \right)=1 .\end{equation}

Subsequently, we can group the coefficients of each power of $ \varepsilon $ and set them to zero. This condition is then transferred back to the undeformed surface r = 1:

(3.7a,b) \begin{equation} c^{0}(1,\theta )=1,\quad c_{r}^{0}(1,\theta )\,P_{n}(\cos \theta )+c^{1}(1,\theta )=0. \end{equation}

Initially, we try to obtain c ⁰ from the Laplace equation ${\Delta} c^{0}=0$ . In our investigation, we consider fixed contact angle 90°, allowing us to transform this problem into the evaporation of hemispherical droplets. Given that the droplets exhibit symmetry with respect to the substrate, it is evident that the evaporation flux is uniform throughout; specifically,

(3.8) \begin{equation} J^{0}(\theta)=1\quad\text{at }r=1. \end{equation}

To derive the expression for $J^{0}(\theta )$ at the deformed surface r = r _s , we also expand $c^{0}(r=r_{s},\theta )$ in a Taylor series, and represent it by two orthogonal components in polar coordinates:

(3.9) \begin{align} &J^{0}(\theta)={}\nonumber\\ & \! -\frac{\partial }{\partial r}\cos \alpha \!\left[\!c^{0}\! \left(r,\theta \right)+c_{r}^{0}\! \left(r,\theta \right)\varepsilon\, P_{n}\left(\cos \theta \right)+\cdots +\frac{c_{r}^{\left(m\right)}\left(1,\theta \right)}{m!}\varepsilon ^{m}\,P_{n}\left(\cos \theta \right)^{m}+\cdots \!\right]\!\Bigg|_{r=1}\nonumber\\[6pt] & \! -\frac{\partial }{\partial \theta }\sin \alpha \!\left[\!c^{0}\! \left(r,\theta \right)+c_{r}^{0}\! \left(r,\theta \right)\varepsilon\, P_{n}\left(\cos \theta \right)+\cdots +\frac{c_{r}^{\left(m\right)}\left(1,\theta \right)}{m!}\varepsilon ^{m}P_{n}\left(\cos \theta \right)^{m}+\cdots \!\right]\!\Bigg|_{r=1}, \nonumber\\\end{align}

where the subscript r indicates the derivative with respect to r, superscript numbers in parentheses represent the order of the derivative, and $ \alpha $ is the angle starting from the vector in the r direction $\boldsymbol{n}_r$ to the normal vector $\boldsymbol{n}_i$ of the surface (see figure 2). In polar coordinates, the angle between the tangent line on a curve and the line between the tangent point and the pole (see figure 2) can be expressed as

Figure 2.

Definitions of the intersection angles $ \alpha $ and $ \psi $ . In polar coordinates, for a point P(r, $ \theta $ ) of the curve, the normal vector at that point is n _i , the tangent vector is $ \unicode{x1D6D5} $ _i , and the vector in the direction of the line connecting the pole and the point P is n _r . The intersection angle $ \alpha $ of the point P is from the vector n _i to the vector n _r , while the intersection angle $ \psi $ is from the vector n _r to the vector $ \unicode{x1D6D5} $ _i . As both $ \alpha $ and $ \psi $ start from n _r , and the direction of n _r is always between n _i and $ \unicode{x1D6D5} $ _i , $ \alpha $ and $ \psi $ are measured in opposite directions: $ \alpha $ is measured clockwise, while $ \psi $ is measured anticlockwise.

(3.10) \begin{equation} \tan \psi =\frac{r_{s}\left(\theta \right)}{r'_{s}\left(\theta \right)},\quad-\frac{{\unicode{x03C0}} }{2}\leq \psi \leq \frac{{\unicode{x03C0}} }{2}, \end{equation}

where $r'_{s}(\theta )$ represents the derivative of $r_{s}(\theta )$ with respect to $ \theta $ . We can determine the angle $ \alpha $ by

(3.11a,b) \begin{equation} \cos \alpha =\frac{r_{s}\left(\theta \right)}{\sqrt{r_{s}\left(\theta \right)^{2}+r'_{s}\left(\theta \right)^{2}}},\quad \sin \alpha =\frac{r'_{s}\left(\theta \right)}{\sqrt{r_{s}\left(\theta \right)^{2}+r'_{s}\left(\theta \right)^{2}}}. \end{equation}

From (3.7a ) and (3.8), we can deduce that

(3.12) \begin{equation} \frac{\partial c^{0}\left(r,\theta \right)}{\partial \theta }\Bigg|_{r=1}=\frac{\partial ^{2}c^{0}\left(r,\theta \right)}{\partial \theta ^{2}}\Bigg|_{r=1}=0 .\end{equation}

Substituting (3.12) into Laplace’s equation yields

(3.13) \begin{equation} \frac{\partial ^{2}c^{0}\left(r,\theta \right)}{\partial r^{2}}\Bigg|_{r=1}=-2\frac{\partial c^{0}\left(r,\theta \right)}{\partial r}\Bigg|_{r=1}=2 .\end{equation}

In order to obtain higher-order derivatives of c ⁰ with respect to r, we can differentiate both sides of Laplace’s equation with respect to r. All terms with partial derivatives with respect to $\theta$ are equal to 0. This gives $c_{r}^{(3)}(1,\theta )=-2c_{r}^{(2)}(1,\theta )$ . The remainder of the process can be carried out in a similar manner, leading us to conclude that

(3.14) \begin{equation} c_{r}^{\left(m\right)}\left(1,\theta \right)=\left(-1\right)^{m}2^{m-1} .\end{equation}

With these conclusions, (3.9) becomes

(3.15) \begin{align} & J^{0}\left(\theta \right) = \cos \alpha -2\varepsilon\, P_{n}\left(\cos \theta \right)\cos \alpha +\varepsilon\, P_{n}'\left(\cos \theta \right)\cos \alpha +\cdots\nonumber \\ &=\!\sum _{m=0}^{\infty } \!\left\{\!\frac{\left(-2\right)^{m}}{m!}\varepsilon ^{m}\left[P_{n}\left(\cos \theta \right)\right]^{m}\cos \alpha -\frac{\left(-2\right)^{m-1}}{m!}\varepsilon ^{m}m\left[P_{n}\left(\cos \theta \right)\right]^{m-1}P_{n}'\left(\cos \theta \right)\cos \alpha \!\right\}\!. \nonumber \\\end{align}

Then we proceed to consider c ¹. Given the solution for c ⁰, the boundary condition of c ¹ (3.7b ) is expressed as

(3.16) \begin{equation} c^{1}\left(1,\theta \right)=P_{n}\left(\cos \theta \right) .\end{equation}

This Laplace’s equation with a Dirichlet boundary value problem can be readily solved using the method of separation of variables. Let us express the right-hand side of (3.16) using a Fourier–Legendre series expansion, also known as a generalised Fourier series expansion:

(3.17) \begin{equation} P_{n}\left(\cos \theta \right)=\sum _{k=0}^{\infty }A_{k}\displaystyle\frac{1}{r^{k+1}}P_{k}\left(\cos \theta \right) \! .\end{equation}

Apparently, when k ≠ n, A _k = 0, and when k = n, A _n = 1. The solution forc ¹ is

(3.18) \begin{equation} c^{1}\left(r,\theta \right)=\displaystyle\frac{1}{r^{n+1}}P_{n}\left(\cos \theta \right) \!.\end{equation}

Similar to J ⁰, J ¹ is given by

(3.19) \begin{align} J^{1}\left(\theta \right)&=\left[-\frac{\partial }{\partial r}c^{1}\left(r,\theta \right)\cos \alpha -\frac{\partial }{\partial \theta }c^{1}\left(r,\theta \right)\sin \alpha \right]\Bigg|_{r=1}\nonumber\\[4pt] &=\frac{n+1}{r^{n+2}}P_{n}\left(\cos \theta \right)\cos \alpha -\frac{1}{r^{n+1}}P_{n}'\left(\cos \theta \right)\sin \alpha\nonumber \\ &=\frac{n+1}{\left[1+\varepsilon\, P_{n}\left(\cos \theta \right)\right]^{n+2}}P_{n}\left(\cos \theta \right)\cos \alpha -\frac{1}{\left[1+\varepsilon\, P_{n}\left(\cos \theta \right)\right]^{n+1}}P_{n}'\left(\cos \theta \right)\sin \alpha . \end{align}

Combining (3.15) with (3.19), the asymptotic expansion of the evaporation flux is

(3.20) \begin{align} & J (\theta )\sim{} \!\!\sum _{m=0}^{\infty }\frac{\left(-2\right)^{m}}{m!}\varepsilon ^{m}\left[P_{n}\left(\cos \theta \right)\right]^{m}\cos \alpha\nonumber\\ &\quad{} -\frac{\left(-2\right)^{m-1}}{m!}\varepsilon ^{m}m\left[P_{n}\left(\cos \theta \right)\right]^{m-1}P_{n}'\left(\cos \theta \right)\cos \alpha\nonumber\\ &\quad{} +\varepsilon \left[\frac{n+1}{r^{n+2}}P_{n}\left(\cos \theta \right)\cos \alpha -\frac{1}{r^{n+1}}P_{n}'\left(\cos \theta \right)\sin \alpha \right]. \end{align}

With this result, we can predict the evaporation flux of non-spherical-cap droplets that are described by Legendre polynomials in any even order within a certain range of deformation amplitude $ \varepsilon $ . As $ \varepsilon $ approaches zero, the shape of the droplet becomes a hemisphere, and (3.20) converges to a constant, which matches the result obtained by Popov (Reference Popov2005).

3.2. Results of asymptotic prediction

The evaporation fluxes of these non-spherical-cap droplets predicted by (3.20) are shown in figure 3. Results are displayed for some different deformation amplitudes across three lowest even-shaped modes: n = 2, n = 4 and n = 6. From figure 3, it is obvious that a larger deformation amplitude and higher shape modes both result in greater non-uniformity in the evaporation flux. For an oscillating droplet, a larger $ \varepsilon $ represents a higher voltage, which can cause the droplet to become unstable. According to the experiments of Oh et al. (Reference Oh, Ko and Kang2008), when the voltage applied to a droplet with the shape mode of n = 2 increased to 100 V, correspondingly $ \varepsilon $ = 3.42 × 10^–1. A relatively large value of $ \varepsilon $ means a large deformation of the droplet. Similarly, a droplet with a large deformation will be difficult to maintain in its original shape in our problem. It is unnecessary to discuss scenarios involving large values of $ \varepsilon $ , as they lack practical significance. Consequently, we limit the value of $ \varepsilon $ within a certain range, and the width of this range decreases as the mode order increases.

Figure 3.

Scaled plots of six droplets, each characterised by (3.20), with the same undeformed radius but varying deformation amplitudes $ \varepsilon $ and shape modes n. Specifically: (a) $ \varepsilon $ = 0.1, n = 2; (b) $ \varepsilon $ = 0.2, n = 2; (c) $ \varepsilon $ = 0.1, n = 4; (d) $ \varepsilon $ = 0.12, n = 4; (e) $ \varepsilon $ = 0.08, n = 6; (f) $ \varepsilon $ = 0.1, n = 6. The length of an arrow represents the value of the evaporation flux. In (a), the length of the reference arrow corresponds to the unit evaporation flux.

To confirm the validity of the prediction, we cross-reference with direct numerical simulations (DNS) of the complete system (2.4)–(2.6) conducted using COMSOL Multiphysics (5.6) (refer to Appendix A for additional information). Figure 4 demonstrates reasonable agreement in the evaporation flux between our asymptotic prediction and the results obtained from DNS.

Figure 4.

Comparisons of evaporation flux between results of the asymptotic prediction (3.20) and DNS results from COMSOL for droplets with different deformation amplitudes and shape modes. Specifically: (a) n = 2, $ \varepsilon $ = 0.08; (b) n = 4, $ \varepsilon $ = 0.05; (c) n = 6, $ \varepsilon $ = 0.03; (d) n = 8, $ \varepsilon $ = 0.02; (e) n = 2, $ \varepsilon $ = 0.05; (f) n = 2, $ \varepsilon $ = 0.01.

We use the root mean square (RMS) to assess the relative error as it offers a more transparent and comprehensive measure of the overall deviations, particularly for non-monotonic features. The average relative error of (3.20) can be calculated by

(3.21) \begin{equation} e\left(n,\varepsilon \right)=\frac{\sqrt{\sum _{i=1}^{N}\left[J_{theo}\left(\theta _{i},n,\varepsilon \right)-J_{num}\left(\theta _{i},n,\varepsilon \right)\right]^{2}}}{\sqrt{\sum _{i=1}^{N}J_{num}^{2}\left(\theta _{i},n,\varepsilon \right)}},\quad 0\leq \theta _{i}\leq \frac{{\unicode{x03C0}} }{2}, \end{equation}

where $ \theta $ _i represents the data points in DNS, N is the total number of data points, and J _num and J _theo are the evaporation fluxes of our theoretical result (3.20) and DNS result, respectively. The average relative errors of the six examples in figure 4 are given by

(3.22a–f) \begin{align} & e(2,0.08)=1.77\,\% ,\quad e(4,0.05)=2.48\,\% ,\quad e(6,0.03)=2.27\,\% ,\nonumber\\ & e(8,0.02)=2.13\,\% ,\quad e(2,0.05)=1.25\,\% ,\quad e(2,0.01)=0.31\,\% , \end{align}

respectively. Under the assumption of small deformation, the asymptotic solution (3.20) exhibits good accuracy, with the average relative error remaining within 5 %. As shown in figure 4(a,e,f) and the RMS errors in (3.22), the asymptotic predictions and numerical results increasingly align as $ \varepsilon $ decreases, both converging to spatially uniform flux.

4.1. Boundary conditions

Given that the contact angles of our droplets are fixed at 90° and the droplets are non-thin, it is evident that the lubrication approximation cannot be applied in our cases. When the evaporation rate of the droplet is relatively high but the viscosity is low, and the inertial forces induced by evaporation-driven internal flow significantly exceed the viscous forces (i.e. the Reynolds number is much greater than 1), an inviscid flow assumption becomes highly suitable for describing the internal flow field. Let us focus solely on potential flow. Masoud & Felske (Reference Masoud and Felske2009a ) made significant contribution to the inviscid flow problem for spherical-cap droplets. They obtained an analytical solution for inviscid flow inside evaporating sessile droplets and for distributions of evaporative flux, both with arbitrary contact angles. They examined both the pinned contact line (the CR mode) and the unpinned contact line (the CA mode). We assume that the droplet maintains its non-spherical-cap shape and shrinks steadily during the evaporation process. This is reasonable, as the external forces required to maintain the droplet’s non-spherical-cap shape typically do not diminish with evaporation. Therefore, the CA mode is adopted here.

The potential flow can be described by

(4.1) \begin{equation} {\Delta} \overline{\phi }\left(r,\theta \right)=0,\quad r\leq r_{s}\text{ and }0\leq \theta \leq {\unicode{x03C0}} /2, \end{equation}

where $\overline{\phi }$ is the potential of the velocity field. We will also solve this Laplace’s equation in spherical coordinates. The physical parameters in this section are non-dimensionalised as

(4.2a,b,c) \begin{equation} \phi =\frac{\overline{\phi }\rho ^{2}}{D\left(\overline{c}_{s}-\overline{c}_{\infty }\right)},\quad u=\frac{\overline{u}\rho ^{2}R}{D\left(\overline{c}_{s}-\overline{c}_{\infty }\right)},\quad t=\frac{\overline{t}\rho ^{2}}{D\left(\overline{c}_{s}-\overline{c}_{\infty }\right)}, \end{equation}

where $\overline{u}$ denotes the velocities, and $\overline{t}$ represents time.

It is notable that the velocity components normal to the axis of symmetry and normal to the solid surface vanish. Therefore, the boundary condition of potential function $\phi$ can be written as

(4.3) \begin{equation} \phi _{\theta }\left(r,0\right)=\phi _{\theta }\left(r,{\unicode{x03C0}} /2\right)=0, \end{equation}

where the subscript $ \theta $ indicates the derivative with respect to $ \theta $ . Another boundary condition arises from the conservation of fluid, which can be expressed as

(4.4) \begin{equation} u_{s}=u_{J}+u_{f}, \end{equation}

where u _s represents the velocity of the free surface, u _J denotes the evaporation-induced radial flow velocity, and u _f is the fluid radial flow velocity. The relationship between these three quantities is depicted more clearly in figure 5(a). As the liquid in the droplet evaporates along the normal direction of the free surface, the droplet surface takes on the shape depicted by the dotted line in figure 5(a), while the droplet surface must maintain the non-spherical-cap shape described by (2.1). As a result, the fluid moves radially to compensate for the radial difference between these two shapes. The expression of the fluid radial flow velocity u _f is the Neumann boundary condition, which we need to determine.

Figure 5.

(a) Relationship between the evaporation-induced radial flow velocity u _J (indicated by black arrows), the fluid radial flow velocity u _f (indicated by red arrows) and the free surface velocity u _s (the vector sum of the black and red arrows). The solid line represents the free surface of the droplet at the initial moment. The dashed line represents the free surface after a short time period dt, which can still be described by (2.1) and has a smaller undeformed radius. The dotted line represents the droplet’s surface shape after dt, considering only the normal evaporation-induced deformation without accounting for fluid flow. The fluid must flow radially to compensate for the radial difference between the dashed line and the dotted line, thereby maintaining the shape of the droplet as depicted by the dashed line. (b) Relationship between the evaporation-induced radial flow velocity u _J and the normal deformation velocity, whose value is J ₀. The angle $ \beta $ is defined as the intersection angle between these two velocities. When dt approaches zero, the intersection angle $ \beta $ approaches $ \alpha $ .

For the sake of result simplicity, we will begin by assuming the evaporation model to be the spatially uniform evaporation, where the evaporation flux J remains constant, denoted J ₀. The droplet experiences the mass loss due to evaporation

(4.5) \begin{equation} \int J\,{\mathrm{d}}S=J_{0}S=-\frac{\mathrm{d}V_{L}}{\mathrm{d}t}, \end{equation}

where S is the surface area of the droplet, and V _L is the droplet volume. These values can be calculated by

(4.6) \begin{equation} S=2{\unicode{x03C0}} \int _{0}^{{{\unicode{x03C0}} }/{2}}r_{s}\left(\theta ,t\right)\sin \theta\, \sqrt{r_{s}^{2}\left(\theta ,t\right)+r_{s}'^{2}\left(\theta ,t\right)}\,\mathrm{d}\theta, \end{equation}

(4.7) \begin{equation} V_{L}=\frac{2}{3}{\unicode{x03C0}} \int _{0}^{{{\unicode{x03C0}} }/{2}}r_{s}^{3}\left(\theta ,t\right)\sin \theta\, \mathrm{d}\theta =\frac{2}{3}{\unicode{x03C0}} \int _{0}^{{{\unicode{x03C0}} }/{2}}\left[1+\varepsilon\, P_{n}\left(\cos \theta \right)\right]^{3}r_{0}^{3}\left(t\right)\sin \theta\, \mathrm{d}\theta, \end{equation}

where r ₀(t) denotes the function of the radius of the undeformed droplet over time. Since we are only considering the initial moment of the evaporation process in this section, we set $r_{0}(t=0)=1$ . The rate of volume loss is expressed as

(4.8) \begin{equation} \frac{\mathrm{d}V_{L}}{\mathrm{d}t}=\frac{\mathrm{d}V_{L}}{\mathrm{d}r_{0}}\frac{\mathrm{d}r_{0}}{\mathrm{d}t}=-J_{0}S, \end{equation}

which yields

(4.9) \begin{equation} \frac{\mathrm{d}r_{0}}{\mathrm{d}t}=-\frac{J_{0}S}{\mathrm{d}V_{L}/\mathrm{d}r_{0}}, \end{equation}

where $\mathrm{d}V_{L}/\mathrm{d}r_{0}$ can be readily obtained by differentiating (4.7) with respect to r ₀. The free surface velocity u _s is given by

(4.10) \begin{equation} u_{s}=-\frac{\mathrm{d}r_{0}}{\mathrm{d}t}\left[1+\varepsilon\, P_{n}\left(\cos \theta \right)\right]=-\frac{J_{0}S}{\mathrm{d}V_{L}/\mathrm{d}r_{0}}\left[1+\varepsilon\, P_{n}\left(\cos \theta \right)\right] .\end{equation}

As shown in figure 5(b), when the time dt approaches zero, the intersection angle $ \beta $ approaches $ \alpha $ . Then we have

(4.11) \begin{equation} u_{J}=-\frac{J_{0}}{\cos \alpha } .\end{equation}

Substituting (4.10) and (4.11) into (4.4), we can write the Neumann boundary condition of (4.1) as

(4.12) \begin{equation} \phi _{r}\left(r=r_{s},\theta \right)=u_{f}=J_{0}\left\{\frac{1}{\cos \alpha }-\frac{S}{\mathrm{d}V_{L}/\mathrm{d}r_{0}}\left[1+\varepsilon\, P_{n}\left(\cos \theta \right)\right]\right\} .\end{equation}

4.2. Asymptotic analysis

Similarly, expand the potential function as an asymptotic series

(4.13) \begin{equation} \phi \left(r,\theta ,\varepsilon \right)\sim \phi ^{0}\left(r,\theta \right)+\varepsilon\, \phi ^{1}\left(r,\theta \right) .\end{equation}

Both $\phi ^{0}$ and $\phi ^{1}$ satisfy Laplace’s equation separately. Substituting the expansion series (4.13) into the Dirichlet boundary condition (4.3), we get

(4.14a,b) \begin{equation} \phi _{\theta }^{0}(r,0)=\phi _{\theta }^{1}(r,0)=0,\quad \phi _{\theta }^{0}(r,{\unicode{x03C0}} /2)=\phi _{\theta }^{1}(r,{\unicode{x03C0}} /2)=0, \end{equation}

since $ \varepsilon $ is independent. The Neumann boundary condition (4.12) becomes

(4.15) \begin{equation} \phi _{r}^{0}\left(r=r_{s},\theta \right)+\varepsilon\, \phi _{r}^{1}\left(r=r_{s},\theta \right)=u_{f} .\end{equation}

Also, expanding each term of the left-hand side in a Taylor series about r = 1 to first order in $ \varepsilon $ , we obtain

(4.16) \begin{equation} \phi _{r}^{0}\left(1,\theta \right)+\phi _{rr}^{0}\left(1,\theta \right)\varepsilon\, P_{n}\left(\cos \theta \right)+\varepsilon \left[\phi _{r}^{1}\left(1,\theta \right)+\phi _{rr}^{1}\left(1,\theta \right)\varepsilon\, P_{n}\left(\cos \theta \right)\right]=u_{f} .\end{equation}

Grouping the coefficients of each power of $ \varepsilon $ and setting them to zero, the boundary conditions for $\phi ^{0}$ and $\phi ^{1}$ are given by

(4.17a,b) \begin{equation} \phi _{r}^{0}(r=1,\theta )=u_{f},\quad \phi _{rr}^{0}(1,\theta )\,P_{n}(\cos \theta )+\phi _{r}^{1}(1,\theta )=0. \end{equation}

First, let us consider the mixed boundary condition of Laplace’s equation with respect to $\phi ^{0}$ . The solution for $\phi ^{0}$ can be expressed in the form

(4.18) \begin{equation} \phi ^{0}\left(r,\theta \right)=\sum _{k=0}^{\infty }A_{k}r^{k}\,P_{k}\left(\cos \theta \right) .\end{equation}

Thus

(4.19) \begin{equation} \phi _{r}^{0}\left(r,\theta \right)=\sum _{k=0}^{\infty }kA_{k}r^{k-1}\,P_{k}\left(\cos \theta \right) .\end{equation}

Now we need to determine the coefficient of each term A _k . Substituting (4.12) and (4.19) into (4.17a ) yields

(4.20) \begin{equation} \sum _{k=0}^{\infty }kA_{k}\,P_{k}\left(\cos \theta \right)=J_{0}\left\{\frac{1}{\cos \alpha }-\frac{S}{\mathrm{d}V_{L}/\mathrm{d}r_{0}}\left[1+\varepsilon\, P_{n}\left(\cos \theta \right)\right]\right\} .\end{equation}

We express the right-hand side of (4.20) using a Fourier–Legendre series expansion

(4.21) \begin{equation} J_{0}\left\{\frac{1}{\cos \alpha }-\frac{S}{\mathrm{d}V_{L}/\mathrm{d}r_{0}}\left[1+\varepsilon\, P_{n}\left(\cos \theta \right)\right]\right\}=\sum _{k=0}^{\infty }a_{k}\,P_{k}\left(\cos \theta \right), \end{equation}

where

(4.22) \begin{equation} a_{k}=\left(\frac{2k+1}{2}\right)\int _{0}^{{\unicode{x03C0}} }J_{0}\left\{\frac{1}{\cos \alpha }-\frac{S}{\mathrm{d}V_{L}/\mathrm{d}r_{0}}\left[1+\varepsilon\, P_{n}\left(\cos \theta \right)\right]\right\}P_{k}\left(\cos \theta \right)\sin \theta\, \mathrm{d}\theta . \end{equation}

Combining (4.20) with (4.21), we get

(4.23) \begin{equation} A_{k}=\frac{a_{k}}{k} .\end{equation}

Since (4.22) does not have a closed form, it can be handled only through numerical integration. However, the numerical results indicate that most values of A _k are nearly zero, except when k = 2m, where $m=1,2,\ldots,n$ . Therefore, we can express $\phi ^{0}$ as

(4.24) \begin{equation} \phi ^{0}\left(r,\theta \right)\approx \sum _{k=1}^{n}A_{2k}r^{2k}\,P_{2k}\left(\cos \theta \right), \end{equation}

where

(4.25) \begin{equation} A_{2k}=\left(1+\frac{1}{4k}\right)J_{0}\int _{0}^{{\unicode{x03C0}} }\left[\frac{1}{\cos \alpha }-\frac{S}{\mathrm{d}V_{L}/\mathrm{d}r_{0}}\varepsilon\, P_{n}\left(\cos \theta \right)\right]P_{2k}\left(\cos \theta \right)\sin \theta\, \mathrm{d}\theta . \end{equation}

This solution satisfies the Dirichlet boundary condition (4.14) naturally.

The solution for $\phi ^{1}$ has a form similar to (4.18):

(4.26) \begin{equation} \phi ^{1}\left(r,\theta \right)=\sum _{j=0}^{\infty }B_{j}r^{j}\,P_{j}\left(\cos \theta \right) .\end{equation}

We can determine the unknown coefficients through the boundary condition (4.17b ), which can be translated into

(4.27) \begin{equation} \sum _{j=0}^{\infty }B_{j}r^{j}\,P_{j}\left(\cos \theta \right)\approx \sum _{k=1}^{n}A_{2k}r^{2k}\,P_{2k}\left(\cos \theta \right)P_{n}\left(\cos \theta \right), \end{equation}

after substituting the result of $\phi ^{0}$ into (4.17b ). Utilizing a method similar to the one used to determine A _k , B _j can be expressed simply as

(4.28) \begin{equation} B_{j}=-\left(1+\frac{1}{2j}\right)\int _{0}^{{\unicode{x03C0}} }\sum _{k=1}^{n}A_{2k}r^{2k}\,P_{2k}\left(\cos \theta \right)P_{n}\left(\cos \theta \right)P_{j}\left(\cos \theta \right)\sin \theta\, \mathrm{d}\theta .\end{equation}

We define the integral

(4.29) \begin{equation} I=\int _{0}^{{\unicode{x03C0}} }P_{2k}\left(\cos \theta \right)P_{n}\left(\cos \theta \right)P_{j}\left(\cos \theta \right)\sin \theta\, \mathrm{d}\theta .\end{equation}

Due to the orthogonality of Legendre polynomials, the integral I is non-zero only when $j\in \left\{2i\,|\,i\in \mathrm{\mathbb{N}}_{0},\ ({2k-n}/{2})\leq i\leq ({2k+n}/{2})\right\}$ . Additionally, for (4.28), j ≠ 0. Therefore, the amount of computation is greatly reduced. We obtain the solution for $\phi ^{1}$ as

(4.30) \begin{equation} \phi ^{1}\left(r,\theta \right)\approx \sum _{j=1}^{3n/2}B_{2j}r^{2j}\,P_{2j}\left(\cos \theta \right), \end{equation}

where B _2j can be calculated numerically using (4.28), and the calculation process can be greatly simplified by the aforementioned orthogonality.

Combining the solutions for $\phi ^{0}$ and $\phi ^{1}$ by (4.13), we get a semi-analytical expression for potential function $\phi$ . This enables us to describe the internal inviscid flow of an evaporating droplet with the non-spherical-cap shape described by (2.1).

4.3. Results of asymptotic prediction

By utilizing the semi-analytical expression for the potential function $\phi$ , we are able to generate the streamline plots for the aforementioned non-spherical-cap droplets (see figures 6 a,c,e). From the streamline plots, we observe that the flow trend of the fluid is from the protrusions of the droplet to its depressions. This can be attributed to the fact that the droplet surface velocity u _s at its protrusions is higher than at the depressions, while the distribution of the evaporation-induced radial flow velocity u _J at the droplet surface does not change as dramatically. Therefore, the fluid needs to flow from the protrusions to the depressions to balance this velocity difference and maintain the non-spherical-cap shape of the droplet surface. In order to validate the asymptotic prediction, we can compare the radial derivative of potential function $\phi _{r}$ at the droplet shape boundary with the original boundary condition (4.12). This comparison for three droplets is shown in figures 6(b,d,f) indicating strong agreement with our asymptotic prediction, except for positions near the central apex of the droplets (specifically at $ \theta $ = 0), where the difference is also within 5 % of the original boundary condition value for n = 2, $ \varepsilon $ = 0.1 and n = 4, $ \varepsilon $ = 0.06, and within 10 % for n = 6, $ \varepsilon $ = 0.05.

Figure 6.

(a,c,e) Stream plots of non-spherical-cap droplets drawn according to the asymptotic prediction. (b,d,f) Comparisons between the radial derivative of potential function $\phi _{r}$ at the droplet shape boundary and the original boundary condition (4.12). The evaporation model is spatially uniform evaporation with J ₀ = 1. Parameters are: (a,b) n = 2, $ \varepsilon $ = 0.1; (c,d) n = 4, $ \varepsilon $ = 0.06; (e,f) n = 6, $ \varepsilon $ = 0.05.

If we assume the evaporation model to be the diffusion model rather than spatially uniform evaporation, then we can get the asymptotic prediction of potential function in the diffusion evaporation model by simply replacing the constant evaporation flux J ₀ with the asymptotic expansion J in (3.20). For $\phi ^{0}$ , the coefficients A _2k in the diffusion model become

(4.31) \begin{equation} A_{2k,\mathit{diff}}=\left(1+\frac{1}{4k}\right)\int _{0}^{{\unicode{x03C0}} }J\left(\theta \right)\left[\frac{1}{\cos \alpha }-\frac{S}{\mathrm{d}V_{L}/\mathrm{d}r_{0}}\varepsilon\, P_{n}\left(\cos \theta \right)\right]P_{2k}\left(\cos \theta \right)\sin \theta\, \mathrm{d}\theta , \end{equation}

and we still only consider the terms where k = 2m, with m = 1, 2, …, n. For $\phi ^{1}$ , the expression for B _j remains unchanged. In this way, we can also create streamline plots and compare radial derivatives of potential function $\phi _{r}$ under the diffusion evaporation model (see figure 7). From figure 7, We can observe that the pattern of fluid flow remains largely consistent when compared to that depicted in figure 6, with only slight variations in the maximum flow velocity. Moreover, the errors of our asymptotic prediction remain within the specified range overall.

Figure 7.

(a,c,e) Stream plots of non-spherical-cap droplets predicted by the asymptotic analysis under the diffusion evaporation model. (b,d,f) Comparisons of radial derivatives of potential function $\phi _{r}$ . The blue solid line represents the original boundary condition, and the line of black circles represents our asymptotic prediction, under the diffusion evaporation model. The line of red points represents the boundary condition under spatially uniform evaporation when J ₀ = 1. Parameters are: (a,b) n = 2, $ \varepsilon $ = 0.1; (c,d) n = 4, $ \varepsilon $ = 0.06; (e,f) n = 6, $ \varepsilon $ = 0.05.