FormalPara Greek symbols

σSurface tension N/m

ηDynamic viscosity kg ⋅ (m ⋅ s)−1

αThermal diffusivity m2/s

τRelaxation time for stream field

τTRelaxation time for temperature field

ωiWeight coefficients in i direction

φPhase change term

ρDensity kg ⋅ m−3

λHeat conductivity W ⋅ (m ⋅ K)−1

υKinetic viscosity m2 ⋅ s−1.

ψPseudo-potential of components.

FormalPara Super- and sub- scripts

iLattice direction

σLiquid/vapor component

sSolid phase

1 Introduction

Droplet evaporation has attracted increasing attention in recent decades, with its wide application in the industrial field. Commonly seen applications of droplet evaporation are Inkjet printing [1, 2], spray cooling [3], material processing [4], coffee-ring effect [5,6,7,8] and particle synthesis [9, 10]. In recent years, more medical applications have been developed for droplet evaporation, like DNA/RNA arrangement [11, 12] and medical diagnosis [13, 14].

Among the factors affecting the droplet evaporation process, Marangoni flow plays a vital role in influencing the inner dropwise flow pattern, hence affects the way heat is transferred and also the evaporation behaviour [15]. In recent years, many scholars have been working on studying the relationship between the Marangoni effect and the droplet evaporation process. Hu et al. [16, 17] studied the effect of Marangoni flow on the well-known “coffee-ring” effect and found that one of the key requirement for the formation of “coffee-ring” is the suppression of Marangoni flow. They also found that surfactant contamination has a good performance in suppressing the Marangoni effect. Anna-Lena Ljung et al. [18] adopted a numerical approach in investigating the role of boundary conditions on flow pattern both inside and outside the droplet. They found that the existence of the Marangoni effect increases the velocity within the droplet with up to three orders of magnitude. F. Girard et al. [19] studied the effect of Marangoni flow on the evaporation of volatile heated drops numerically when casting heat plate heating. They concluded that the impact brought by Marangoni flow on the evaporation rate could be neglected. The effect of Marangoni flow on the evaporation process is also investigated by R. Savino et al. [20] using silicone oils with different viscosities and hydrocarbons, and it’s believed that a larger, more uniform surface temperature can be caused by Marangoni flow thus the evaporation rate is increased eventually.

As from the formula of the Marangoni number \( {Ma}_T=-\frac{\mathrm{d}\sigma }{\mathrm{d}T}\frac{l\Delta T}{\eta \alpha} \), it can be easily concluded that the temperature gradient generates Marangoni flow on the fluid interface, hence when local heating is applied for a droplet, the more complicated surface temperature gradient will lead to a more complex flow pattern in the dropwise, changing the way heat is transferred and thus the evaporation rate. However, in most researches, the heating conditions are either heat plate or ambient heating, while local heating for a droplet is yet to be studied. Moreover, due to the small size of the droplet, there’s a significant limitation on studying this experimentally. As a result, the numerical approach is chosen to investigate the effect of local heating on the droplet evaporation process in this paper.

Lattice Boltzmann method (LBM) has been adopted by many researchers in their researches toward droplets for its easy parallelism, easy processing of boundary conditions, and easy implementation of programs. Yan et al. [21] did a simulation on liquid droplet behaviour on partial wetting surfaces while the liquid-gas density ratio is large. By combining Inamuro’s [22] and Briant’s [23] model, the study of a droplet falling onto a hydrophilic surface with hydrophobic strips was done. The results indicated that the current LBM approach could be used as a reliable way to study fluidic control on heterogeneous surfaces and other wetting related subjects. Mohammad Taghilou [24] et al. applied a model which uses the Cahn-Hilliard diffuse interface theory, to capture the liquid-gas interface. The passive scalar model is used to simulate the thermal effects. It’s found that by increasing the Prandtl number ratio between the droplet and its surrounding, thermal diffusion within the droplet will be delayed and this causes reduction in the droplet average temperature. Qing Guo [25] et al. studied the effects of ambient humidity of air and wall temperatures on sessile droplet evaporation by applying the multi-component multi-phase (MCMP) LBM model. The evaporation behaviour near the Leidenfrost temperature was studied, including the flow pattern and the heat flux. Salman et al. [26] gave a short review of the way density ratio and achievable temperature broke their limits during the past decades. The relationship between temperature, reference relaxation time, density ratios, a reduced parameter for the equation of state (EOS) and the interface thickness was given, to help lift the limitations when dealing with cases at a lower temperature than the Leidenfrost point.

On the other hand, surface wettability determines the droplet behaviour in various ways, from affecting the inner droplet flow patterns to affecting the evaporation rate. Y. Takata [27] et al. did some experiments on the evaporation of water drops on plasma-irradiated hydrophilic surfaces. The contact angle was measured, and it was found that the evaporation time decreases while the wetting limit and the Leidenfrost temperatures increase with the decrease of contact angle. Droplet morphology during the evaporation process on superhydrophobic surfaces is also researched by many researchers [28,29,30], and their unique behaviour is compared to that on a common surface. B. Sobac et al. [31] tested the evaporation of droplet with a wide range of contact angles and found that when the contact angle is small, the contact radius is the only factor to the evaporation rate.

In this paper, a pseudo-potential multi-component multi-phase (MCMP) lattice Boltzmann method (LBM) model is adopted, and a simulation of droplet evaporation with different local heating conditions is done. Substrate wettability is also considered to study the effect of contact angle on the evaporation process. Flow patterns and temperature distribution inside the droplet, as well as the droplet morphology and evaporation rate, are studied, and the attention is also given to the droplet slipping caused by asymmetric heating source.

2 Numerical methodology

The multi-component multiphase lattice Boltzmann method, which is adopted in this paper, is introduced in this section. Given that the droplet is an ideally perfect hemispherical, the flow pattern caused by Marangoni flow should be accordingly centrosymmetric when the heat source is evenly distributed on the substrate, hence on any section through the centre point the flow pattern should be constant. Based on that assumption, a 2-D model is adopted considering its simplicity and low demand for computing resource.

The Lattice Boltzmann Equation (LBE) for density and temperature are listed below, respectively, which contains collision and streaming steps:

$$ {f}_{\sigma, i}\left(x+{e}_i\delta t,t+\delta t\right)-{f}_{\sigma, i}\left(x,t\right)=-\frac{1}{\tau_{\sigma }}\left({f}_{\sigma, i}\left(x,t\right)-{f}_{\sigma, i}^{\mathrm{eq}}\left(x,t\right)\right) $$
(1)
$$ {g}_{\sigma, i}\left(x+{e}_i\delta t,t+\delta t\right)-{g}_{\sigma, i}\left(x,t\right)=-\frac{1}{\tau_{\sigma, \mathrm{T}}}\left({g}_{\sigma, i}\left(x,t\right)-{g}_{\sigma, i}^{\mathrm{eq}}\left(x,t\right)\right)+\delta t{\omega}_i{\varphi}_{\sigma } $$
(2)

where τσ is the dimensionless collision relation time of the σth component, ei is the lattice velocity vector and i stands for the lattice velocity direction. σ stands for different components, as 1 stands for base fluid while 2 stands for nanoparticle phase. fσ, i(x, t) and \( {f}_{\sigma, i}^{eq}\left(x,t\right) \) are the distribution function and equilibrium distribution function of the σth component with the velocity ei at lattice x and time t, respectively. The discrete velocity for each direction is:

$$ {e}_{\mathrm{i}}=\Big\{{\displaystyle \begin{array}{cc}\left(0,0\right)& i=0\\ {}c\left(\cos \frac{\left(i-1\right)\pi }{2},\sin \frac{\left(i-1\right)\pi }{2}\right)& i=1,\dots, 4\\ {}\sqrt{2}c\left(\cos \frac{\left(2i-9\right)\pi }{4},\sin \frac{\left(2i-9\right)\pi }{4}\right)& i=5,\dots, 8\end{array}} $$
(3)

where c is the reference lattice velocity.

Corresponding to the D2Q9 model, the velocity equilibrium distribution function and temperature equilibrium distribution function for each component is shown as follows:

$$ {f}_{\sigma, i}^{\mathrm{eq}}={\omega}_i{\rho}_{\sigma}\left[1+3\frac{e_i\cdot {\mathrm{u}}_{\sigma}^{eq}}{c^2}+\frac{9}{2}\frac{{\left({e}_i\cdot {\mathrm{u}}_{\sigma}^{eq}\right)}^2}{c^4}-\frac{3}{2}\frac{{\left({\mathrm{u}}_{\sigma}^{eq}\right)}^2}{c^2}\right] $$
(4)
$$ {g}_{\sigma, i}^{\mathrm{eq}}={\omega}_i\mathrm{T}\left[1+3\frac{e_i\cdot {\mathrm{u}}_{\sigma}^{eq}}{c^2}+\frac{9}{2}\frac{{\left({e}_i\cdot {\mathrm{u}}_{\sigma}^{eq}\right)}^2}{c^4}-\frac{3}{2}\frac{{\left({\mathrm{u}}_{\sigma}^{eq}\right)}^2}{c^2}\right] $$
(5)

where ωi is the weight coefficients.

The phase change source term is given as follows as obtained by Chaoyang et al. [32]:

$$ \varphi =\mathrm{T}\left[1-\frac{1}{\rho {c}_v}{\left(\frac{\mathrm{\partial p}}{\mathrm{\partial T}}\right)}_{\rho}\right]\nabla \cdot U+\left[\frac{1}{\rho {c}_v}\nabla \left(\lambda \nabla \mathrm{T}\right)-\nabla \cdot \left(\frac{\lambda }{\rho {c}_p}\nabla \mathrm{T}\right)\right] $$
(6)

where cv and cpare the thermal specific heat at constant volume and pressure, respectively.

For macroscopic density and velocity of each component, they can be calculated as

$$ {\rho}_{\sigma}\left(x,t\right)=\sum \limits_{i=0}^8{f}_{\sigma, i}\left(x,t\right) $$
(7)
$$ {\mathrm{u}}_{\sigma}\left(x,t\right)=\frac{\sum_i{e}_i{f}_{\sigma, i}\left(x,t\right)}{\rho_{\sigma}\left(x,t\right)}+\frac{\delta t{F}_{\sigma}\left(x,t\right)}{2{\rho}_{\sigma}\left(x,t\right)} $$
(8)

And the total density and velocity are written as

$$ \rho \left(x,t\right)=\sum \limits_{\sigma }{\rho}_{\sigma}\left(x,t\right) $$
(9)
$$ \mathrm{U}\left(x,t\right)=\frac{1}{\rho \left(x,t\right)}\sum \limits_{\sigma }{\rho}_{\sigma}\left(x,t\right){\mathrm{u}}_{\sigma}\left(x,t\right) $$
(10)

For the macroscale temperature, the temperature for each component is set as below:

$$ {\mathrm{T}}_{\sigma}\left(x,t\right)=\sum \limits_{i=0}^8{g}_{\sigma, i}\left(x,t\right) $$
(11)

and the temperature for the mixture is

$$ \mathrm{T}\left(x,t\right)=\frac{1}{\rho \left(x,t\right){c}_v}\sum \limits_{\sigma }{\rho}_{\sigma}\left(x,t\right){c}_{\sigma, v}{\mathrm{T}}_{\sigma}\left(x,t\right) $$
(12)

The kinetic viscosity and thermal diffusion coefficient are defined by

$$ {\nu}_{\sigma }=\frac{1}{3}{c}^2\left({\tau}_{\sigma, f}-\frac{1}{2}\right)\delta t $$
(13)
$$ {\alpha}_{\sigma }=\frac{1}{3}{c}^2\left({\tau}_{\sigma, \mathrm{T}}-\frac{1}{2}\right)\delta t $$
(14)

The Fσ mentioned above in Eq. (8) is the total force acting on the corresponding component σ . In this simulation, four different forces are considered: attractive force between the same components, the repulsive force between different forces, the interactive force between the solid substrate and the fluid phase, and gravity and buoyance force. The thermophoretic force was also considered but was too small compared to other forces hence neglected in this case. The interaction forces between fluid phases are calculated according to Gong and Cheng’s method [33]:

$$ {F}_{\sigma, \sigma }=-{\beta \psi}_{\sigma }(x)\sum \limits_{x\hbox{'}}{G}_{\sigma }{\psi}_{\sigma}\left(x\hbox{'}\right)+\frac{\left(1-\beta \right)}{2}\sum \limits_{x\hbox{'}}{G}_{\sigma }{\psi}_{\sigma}^2\left(x\hbox{'}\right) $$
(15)
$$ {F}_{\sigma, \sigma \prime }=-{\psi}_{\sigma }(x)\sum \limits_{x^{\prime }}{G}_{\sigma \prime }{\psi}_{\sigma \prime}\left(x^{\prime}\right) $$
(16)

whereσ' stands for the other component and x' stands for the neighbouring lattice. And the solid-fluid interaction force is shown below

$$ {F}_{\sigma, s}=-{\psi}_{\sigma }(x)\sum \limits_{x^{\prime }}{G}_ss\left(x^{\prime}\right) $$
(17)

while for s(x') the value varies from 0 to 1, depending on whether the lattice on x' is a fluid phase or solid. The gravity force is given as

$$ {F}_{\sigma, g}=\left({\rho}_{\sigma }-{\rho}_{avg}\right)g $$
(18)

Gσ, Gσ' and Gs are the coefficient of each interaction force. There are calculated by

$$ {G}_{\sigma, \sigma \prime, s}\left(x^{\prime}\right)=\Big\{{\displaystyle \begin{array}{cc}2{g}_{\sigma, \sigma \prime s}& \mid x^{\prime }-x\mid =1\\ {}{g}_{\sigma, \sigma \prime s}/2& \mid x^{\prime }-x\mid =\sqrt{2}\\ {}0& otherwise\end{array}} $$
(19)

The change in temperature affects the pseudo-potential of each component, which is determined by [34]

$$ {\psi}_{\sigma }(x)=\sqrt{\frac{2\left({p}_{\sigma }-{c}_s^2{\rho}_{\sigma}\right)}{c_0{g}_{\sigma }}} $$
(20)

The Peng-Robinson equation of state (P-R EOS) which is very popular is adopted in this paper and is written as

$$ p=\frac{\rho RT}{1- b\rho}-\frac{a{\rho}^2\varepsilon (T)}{1+2 b\rho -{b}^2{\rho}^2} $$
(21)

where \( a=\frac{0.457235{R}^2{T}_c^2}{p_c} \), \( b=\frac{0.077796{RT}_c}{p_c} \) and \( \varepsilon (T)=\left[1+\left(0.37464+1.54226\omega -0.26992{\omega}^2\right)\left(1-\sqrt{T/{T}_c}\right)\right] \), where the subscript c denotes critical state. With ω being the acentric factor, it’s chosen as ω = 0.344, and a, b, R are selected as 3/49, 2/21 and 1 respectively. In this case, the β in Eq. (15) equals to 1.16.

3 Results and discussion

3.1 Model validation

A simulation of a droplet evaporation process on a superhydrophobic surface is cast to provide validation for this numerical model. The contact angle of the droplet is set at around 150 degrees, and the heat source is evenly distributed on the substrate, as shown in Figs. 1, 2. The decrease of droplet volume, as well as the contact angle, is measured and compared to the results obtained by M.J. Gibbons et al. [35] in his paper Local heat transfer to an evaporating superhydrophobic droplet and theoretical result from the Young-Laplace solution.

Fig. 1
figure 1

Schematic of a droplet placed on a superhydrophobic surface. (CA = 150°)

Full size image
Fig. 2
figure 2

Model verification results. (a) Comparison of the volume change throughout the evaporation process. (b) Comparison of the contact angle change throughout the evaporation process

Full size image

It can be observed that the results are in a very close match, with the maximum error between this simulation and Gibbon’s experiment result being 2.53% when regarding the volume and 1.83% when regarding contact angle. Note that in the initial period of the evaporation process, a shape change occurs to the droplet as the droplet is heating up, and with the temperature conducting through the dropwise, the interactive force inside the droplet changes, leading to a decrease in the contact angle. While during the major evaporation process, the deviation of contact angle is less than 10%, so according to S. Dash et al. [36], this period could be considered as constant contact angle (CCA) state. Also, near the end of evaporation, the evaporation mode switches from the CCA model to constant contact radius (CCR) model, and the contact angle decreases rapidly while the curve of the volume is flatting. The morphology of droplet is in an unstable state in those two periods and cannot represent the common tendency of the evaporation process, so the initiate and the last period which covers up about 15% of the evaporation process is not taken into the comparison.

3.2 Simulation of the standard droplet evaporation process

We have also cast a simulation of the standard droplet evaporation process to offer a reference for the local heating study (Fig. 3). In this simulation, a water droplet is placed on substrates which provide three different contact angles, respectively 30, 60 and 90 degrees. The heat source is distributed evenly on the substrate, which means the substrate is at a constant higher temperature Th. The initial ambient temperature is set as the saturated temperature Ts with is lower than the substrate temperature. The temperature setting follows the work of Gong et al. [37], as both Th and Ts are calculated based on the critical temperature Tc, here we set Ts = 0.86Tc, and Th is set at three different temperatures, respectively 0.90Tc, 0.91Tc and 0.92Tc. The top and bottom boundaries are set as bounce-back boundaries, and the left and right boundaries are set as periodic boundaries.

Fig. 3
figure 3

Schematic of the simulation of the standard evaporation process

Full size image

In the initial state of the evaporation, due to the difference in the heat capacity and density of liquid and gas phase, it takes longer for the dropwise to reach the equilibrium temperature. Shown in Fig. 4 is the evolution of isotherms at the initial state of the evaporation process. To allow the system to become converge enough and reduce the error, gravity is added after 40,000 timesteps, and the heat source is added after 80,000 timesteps. The data is collected from 80,000 timesteps on in this research. The isotherm inside and outside the droplet are both nearly horizontal and becomes hierarchical at the droplet interface. With the evaporation proceeding, the isotherm both inside and outside of the dropwise rises from the substrate, representing the heat conduct from the substrate upwards. The height deviation of an isotherm tends to decrease as the evaporation takes place.

Fig. 4
figure 4

Evolution of isotherms at the initial state of evaporation. (a) Time step = 1000. (b) Time step = 5000. (c) Time step = 10,000. (d) Time step = 20,000

Full size image

Shown in Fig. 5 is the morphology of the droplet during the evaporation process. The droplet volume, contact angle and height/diameter ratio (h/d ratio) are studied. The data is recorded every 10,000 timesteps. The unit of the x-axis is ×10000timestep, and the y-axis of Fig.5(a) is the percentage of droplet volume to the initial volume when the heating starts. The unit of y-axis for Fig. 5(b) is degrees. Figure 5(a) shows the evolution of droplet volume to the time of evaporation. It can be observed that apart from the rate of decreasing in the droplet volume, the tendencies are all similar, as the temperature didn’t exceed the Leidenfrost point. For instance, we fitted the curve of dT = 0.5 ⋅ Tc, and the equation for variation in volume is v =  − 0.0003t3 + 0.0082t2 − 0.1379t + 0.9997. That shows the volume decreasing curve keeps flatting, as the contact area between the substrate and the droplet keeps decreasing. Shown in Fig.5(b) and (c), The contact angle and the h/d ratio share a very similar tendency, as when the contact angle decreases, given the droplet volume remains constant, the droplet grows flattered and thus the h/d ratio also decreases. At the initial state of evaporation, due to the increase of the droplet temperature, the interactive force within the droplet decreases as well as the surface tension, leading both the contact angle and the h/d ratio to decrease. While in the rare part of evaporation, with the evaporation mode transferring from the CCA model to the CCR model, the contact angle and the h/d ratio decrease rapidly to complete the evaporation.

Fig. 5
figure 5

Morphology of the droplet in the evaporation process under different heating temperatures. (contact angle = 90 degrees, substrate heating) (a) the volume evolution to time. (b) the contact angle to time. (c) the h/d ratio to time

Full size image

3.3 Simulation of the droplet evaporation under local heating

Marangoni flow has a decisive impact on the inner droplet flow pattern. The flow field can affect the way heat is transferred within the dropwise, thus influences the temperature distribution, droplet motion and morphology, and eventually evaporation rate. In this part, the effect of different local heating conditions on evaporation behaviour is investigated, and different surface wettability is adopted to give the droplet different initial contact angles.

Shown in the figure below is the schematic of this simulation. The droplet is placed on a substrate and surrounded by vapour gas. The initial temperature of the whole cavity simulated is the saturated temperature Ts = 0.86Tc, whereTc is the critical temperature calculated from the equation of state. The top and bottom boundaries are set as bounce-back boundaries, and the temperature of the top boundary is fixed at Ts, while two different temperature are set for three different areas on the bottom boundary. The bottom boundary is divided into three different parts, marked as A, B and C, the temperature of A part and C part are set at Ts while that of the B part is set at Th . Here, to maintain an acceptable evaporation rate to save the computational time, the heating temperatureTh is set at Th = 0.92Tc (Fig. 6).

Fig. 6
figure 6

Schematic of simulation of droplet evaporation under local heating

Full size image

Four different surface wettability are adopted in this study, with the major contact angle in the CCA period of evaporation being at 30°, 60°, 90° and 120° respectively. The size of the heating area B is fixed in this study, while different locations are chosen for the centre point of the heating area. With a fixed droplet volume, the heating area is set at 0.44 times the radius of the droplet when the contact angle is 90°. For simplicity, in this study, we mark the position of the centre of the heating area as Lc, and the length of the heating area is marked as a dimensionless 1. for example Lc = 0 means the heating source is located at the centre of the droplet contact section, and when Lc = 2, the centre of the heating source is located at twice the length of the heating area from the droplet centre to the triple-phase contact point.

Figure 7 shows the streamlines inside the droplet on substrates with different wettability when different values are picked for Lc. Figure 7(a)-(d) show the streamlines in droplets with a contact angle of 30 degrees, while the contact angle for Fig. 7(e)-(h), (i)-(l) and (m)-(o) are 60 degrees, 90 degrees and 120 degrees respectively. Lc is chosen as 0, 0.5, 1 and 1.5 for contact angle from 30° to 90°, while for contact angle 120°, due to the limitation of liquid-solid contact size, three different values 0, 0.25, 0.75 are picked for Lc. Data at time step 100,000 is chosen to ensure the droplet evaporation is in a steady CCA mode. It can be observed that despite different contact angles and different heating source locations, there are always two vortexes aside the heating source. These two are caused by the surface tension difference on the bottom surface, as they are close to the heat source, and the direction of flow near the bottom surface is from where the temperature is high to where the temperature is lower. According to the mechanism of the Marangoni effect, the surface tension decreases when the temperature increases [38] and the flow directions on the droplet interface are determined by surface tension (from where the surface tension is lower to where the surface tension is higher). The same phenomenon can be observed on the top surface of the droplet, as right above the heating source, the temperature is higher than surrounding areas, so the flow is driven from were above the heating source towards the triple-phase contact point on the two sides. Due to the existence of capillary flow the direction of which is from the liquid phase to the solid phase [39], the flow pattern near the triple-phase contact point is always from the bottom upwards, thus resulting in the flow caused by the Marangoni effect kept away from the droplet edge. By comparing the size of the upper vortexes and the lower ones, it can be found that for droplets with smaller contact angles, the flow caused by upper surface tension gradient is more substantial. Also, it can be observed from the lower flow field that for droplets with larger contact angles, when the heat source area is constant, the Marangoni effect is relatively prevailing as the vortexes caused by them are occupying larger area in the lower flow field. This is because for those droplets with larger contact angles, as the h/d ratio of them increases, the capillary flow is weakened compared to droplets with smaller h/d ratios. One more conclusion is that when applying asymmetric local heating conditions, the shape of droplet tend to be asymmetric. This is more obvious for droplets on super-hydrophilic substrates (figure c and d), and the hemisphere where the heat source is located at has a smaller average height than the opposite hemisphere.

Fig. 7
figure 7

Streamlines inside the droplet on substrates with different wettability when applied local heating

Full size image

As captured in the figure below is the temperature distribution of the droplet when local heating is applied. The data is also captured at 100000 timesteps after the heat distribution reaches a steady state. The contact angle and heat source location remain the same as the corresponding figures in Fig. 7. As predictable, the temperature distribution moves with the movement of the heat source, and unlike heat plate heating condition, the major part of the bottom of the droplet remains at a low temperature as the heat mainly transfers upwards, so overall the upper part of the droplet has a higher mean temperature than the bottom part. The highest temperature on the droplet surface is located correspondingly above the heating source, and there is a noticeable temperature jump on the droplet interface, as a difference of heat conductivity and heat capacity exists between the two fluids. More detailed temperature distribution on the droplet interface is presented in Fig. 8.

Fig. 8
figure 8

The temperature distribution in the flow field of droplets with different contact angles under different local heating condition

Full size image

Marangoni number for each instance is also calculated, as is an effective indicator of the strength of Marangoni effect. The results are shown in the figure below. It can be observed that droplet with a smaller contact angel has a bigger Marangoni number than those with bigger contact angles. Meanwhile, as the heating source moves towards the triple-phase contact point, the Marangoni number increases. This is because the Marangoni number is mostly determined by the maximum temperature difference of the droplet surface, and with a smaller contact angle, a bigger temperature difference can be obtained, as shown in Fig. 8. Also, when the heat source is closer to the droplet edge, a similar phenomenon happens.

As presented in Fig. 10, the temperature on the droplet interface changes with the change of the heat source location. The droplet is placed on a substrate with a contact angle = 90°, and four different heat source locations are applied, respectively Lc = 0, 0.5, 1and 1.5. Here on the x-axis, 0.5 indicates the centre of the droplet subface, while 1.0 means the right triple-phase contact point. Temperature data are taken every lattice length unit on the droplet interface, and the average surface temperature is calculated by dividing the sum of temperature data taken by the number of units counted. Results show that with the translation of the heat source position, the curve of temperature profile becomes asymmetric with the peak of the maximum value moving in the direction of the centre to the heat source. Also, due to the two sides of the droplet is thinner compared to the centre, the vertical distance from the heat source to the droplet interface is decreased, leading to a higher maximum interface temperature. This can be observed in Fig. 9. Thus, the average temperature of the droplet interface is also increased, as the average interface temperature for Lc = 0, 0.5, 1and 1.5 are respectively 0.0947907, 0.094803, 0.094823 and 0.094889 lattice temperature unit.

Fig. 9
figure 9

Marangoni number of systems with different contact angles under different heat source locations

Full size image
Fig. 10
figure 10

the droplet interface temperature profile for a droplet with contact angle = 90°

Full size image
Fig. 11
figure 11

evolution of droplet volume when applying local heating. (a) contact angle = 30°; (b) contact angle = 60°; (c) contact angle = 90°

Full size image
Fig. 12
figure 12

evolution of the h/d ratio and contact angle for droplets on substrates with different surface wettability under different local heating conditions. (a)-(c): evolution of h/d ratios of the droplet with basic contact angles at 30°, 60° and 90°; (d)-(f): evolution of contact angles of the droplet with basic contact angles at 30°, 60° and 90°.

Full size image
Fig. 13
figure 13

droplet position evolution under asymmetric local heating conditions when placed on smooth substrate. (a) Contact angle = 30°; (b) Contact angle = 60°; (c) Contact angle = 90°.

Full size image

As the evaporation process majorly takes place at the droplet interface, it’s fair to think that a higher droplet interface temperature will lead to a higher evaporation rate. This is further validated in Fig. 11, which gives the droplet volume evolution to time for droplets under different local heat sources on substrates with different wettability. Note that in this part, due to that the difference in contact angle has an effect on the contact area, and the contact area of the droplet with a contact angle of 120 degrees is too small to fit in Lc = 1.5, so the case of contact angle equals to 120 degrees is not taken into consideration. Also, for the reason that when the heat source position Lc = 1.5 is applied for the droplet with its contact angle being 30°, the shape of the droplet will be asymmetric, leading to errors when measuring the contact angle and calculating the droplet volume, Lc = 1.5 is neither taken into consideration for contact angle = 30°. As could be observed from the figure, despite different surface wettability, the droplet evaporation rate increases when the heat source moves closer to the droplet edge, and this tendency is more evident in those cases with larger contact angle, when the temperature difference on the droplet interface is larger. Moreover, droplets with a smaller contact angle still process a higher evaporation rate, in the same way when applied substrate heating.

The evolution of droplet morphology when applying local heating is also investigated by studying the evolution of h/d ratio and contact angle. The results are presented in Fig. 12. Because the evaporation time is much longer when applying local heating, the results shown in Fig. 12 includes just the front 360,000 timesteps since heating starts. Apart from the same regulation that the h/d ratio shares largely the same tendency as the contact angle, it can be concluded that, while both the h/d ratio and the contact angle go through a decrease in the heating-up state, the decrease is more severe when the heat source is closer to the droplet edge. For hydrophilic surfaces, with the increase of contact angle, the difference in droplet morphology between different heat source locations increases. However, since when asymmetric heating conditions are applied, the droplet shape also more or less goes asymmetric, and this is more obvious for droplets on super-hydrophilic surfaces, errors are likely to exist in this part, and more research could be done towards this in the future.

Considering that the substrate being a smooth surface, when the flow pattern is asymmetric (as presented in Fig. 7), the force acting on the droplet will lead the droplet to slip on the substrate. On the one hand, the flow on the droplet interface does interact with the ambient gas; on the other hand, the asymmetric temperature distribution can result in horizontal flows in the surrounding gas. The effect of asymmetric local heating on the droplet slip is studied in this paper, and the results are as follows:

Note that the case Lc = 1.5 is still excluded in this simulation for better accuracy. It can be concluded from Fig. 13 that the droplet slip phenomenon is more evident on super-hydrophilic substrates; in other words, more evident for smaller contact angles. Also, the slip distance increases when the heating source is closer to the triple-phase contact point. Another phenomenon is that when the contact angles are 30 degrees and 60 degrees, the droplet moves in the opposite direction to which the heat source is placed against the droplet centre, while when the contact angle is 90 degrees, the droplet slip direction is towards the way from the droplet centre to the heat source. Hence, there is a point between 60° and 90°, which when the droplet contact angle exceeds, the droplet slip direction reverses.

4 Conclusion

In this paper, we adopted a pseudo-potential multi-component multiphase lattice Boltzmann method model to investigate the effect of local heating on the droplet evaporation process. The model was successfully validated, and simulations consisting of different heating sources as well as different surface wettability are done. Streamlines, temperature distribution, Marangoni number and droplet morphology are studied. Moreover, an initial study towards the droplet slipping behaviour caused by asymmetric local heating conditions is proposed. Based on the present numerical modelling and analysis, the following can be concluded:

  1. 1.

    Asymmetric local heating conditions can cause asymmetric flow patterns and temperature distributions in the dropwise, and also causes the shape of the droplet being asymmetric. This is more obvious for droplets on super-hydrophilic substrates;

  2. 2.

    The gradient of surface intension is generated both on the upper surface and lower surface of the droplet when the droplet is heated locally, and Marangoni flow is generated on both surfaces. Flow caused by the surface tension gradient on the upper surface is stronger when the contact angle is smaller;

  3. 3.

    Marangoni number increases with the decrease of the droplet contact angle. For the same droplet, when the heating source is closer to the triple-phase contact point, a higher Marangoni number can be achieved.

  4. 4.

    The higher evaporation rate is achieved when the heat source is located closer to the triple-phase contact point. While for the same heat source location, a droplet with smaller contact angles has higher evaporation rates;

  5. 5.

    The droplet slips horizontally on the substrate when an asymmetric local heating condition is applied. The slipping behaviour is affected by both the droplet contact angle and heating source location, and there is a contact angle between 60° and 90°, which is the demarcation point of the slip direction.