Terminal states of thermocapillary migration of a planar droplet at moderate and large Marangoni numbers

TL;DR

This study analyzes the terminal states of thermocapillary droplet migration at various Marangoni numbers using analytical and numerical methods, revealing the dominant physical mechanisms and surface tension effects.

Contribution

It provides a detailed analysis of the terminal states and physical mechanisms of thermocapillary droplet migration at different Marangoni numbers, combining analytical and numerical approaches.

Findings

Steady/unsteady terminal states are identified at different Marangoni numbers.

Surface tension effects vary between the top and bottom surfaces of the droplet.

Thermal convection dominates over conduction at large Marangoni numbers.

Abstract

In this paper, thermocapillary migration of a planar droplet at moderate and large Marangoni numbers is investigated analytically and numerically. By using the dimension-analysis method, the thermal diffusion time scale is determined as the controlling one of the thermocapillary droplet migration system. During this time, the whole thermocapillary migration process is fully developed. By using the front-tracking method, the steady/unsteady states as the terminal ones at moderate/large Marangoni numbers are captured in a longer time scale than the thermal diffusion time scale. In the terminal states, the instantaneous velocity fields in the unsteady migration process at large Marangoni numbers have the forms of the steady ones at moderate Marangoni numbers. However, in view of the former instantaneous temperature fields, the surface tension of the top surface of the droplet gradually becomes the main component of the driving force on the droplet after the inflection point appears. It is different from that the surface tension of the bottom surface of the droplet is the main component of the driving force on the droplet for the latter ones. The physical mechanism of thermocapillary droplet migration can be described as the significance of the thermal convection around the droplet is higher than/just as the thermal conduction across the droplet at large/moderate Marangoni numbers.

Equations

$$v_o=-{\sigma }_{\Theta }G{R}_0/{\mu }_1,$$

$$\begin{array}[]{l}\frac{\partial{\rho_{i}}}{\partial t}+\nabla\cdot(\rho_{i}{\bf v}_{i})=0,\\ \frac{\partial\rho_{i}{\bf v}_{i}}{\partial t}+\nabla\cdot(\rho_{i}{\bf v}_{i}{\bf v}_{i})=-\nabla p_{i}+\frac{1}{Re}\nabla\cdot\mu_{i}(\nabla{\bf v}_{i}+\nabla{\bf v}_{i}^{T})+{\bf f}_{\sigma},\\ \frac{\partial\Theta_{i}}{\partial t}+\nabla\cdot({\bf v}_{i}\Theta_{i})=\frac{1}{Ma}\frac{\kappa_{i}}{k_{i}}\nabla\cdot({k_{i}\nabla\Theta_{i}}),\end{array}$$

$$Re=\frac{{\rho }_1{v}_0{R}_0}{{\mu }_1},Ma=\frac{v_0{R}_0}{{\kappa }_1}.$$

$${\mathbf{v}}_i=0,{\Theta }_i=z$$

$$\begin{array}[]{l}{\bf v}_{1}(x,z_{0})={\bf v}_{1}(x,z_{1})=0,\Theta_{1}(x,z_{0})=z_{0},\Theta_{1}(x,z_{1})=z_{1},\\ u_{i}(0,z)=0,\frac{\partial v_{i}}{\partial x}(0,z)=0,\frac{\partial{\Theta_{i}}}{\partial x}(0,z)=0,\\ {\bf v}_{1}(x_{1},z)=0,\Theta_{1}(x_{1},z)=z.\end{array}$$

$$w_{ij}({\mathbf{r}}_p)=d(x_p-i\Delta x)d(z_p-j\Delta z),$$

$$d(r)=\left\{\begin{array}[]{ll}(1/4\Delta r)[1+\cos(\pi r/2\Delta r)],&|r|<2\Delta r,\\ 0,&|r|\geq 2\Delta r,\end{array}\right.$$

$$\begin{array}[]{ll}\delta{\bf f}_{\sigma}&=\int_{\Delta s}\frac{\partial}{\partial s}(\sigma{\bf\tau})ds(R_{0}/\rho_{1}v^{2}_{0})/(R^{2}_{0}\delta x\delta z)\\ &=[(\sigma{\bf\tau})_{2}-(\sigma{\bf\tau})_{1}]/\rho_{1}v^{2}_{0}R_{0}\delta x\delta z\\ &=\{[(1/Ca-\Theta){\bf\tau}]_{2}-[(1/Ca-\Theta){\bf\tau}]_{1}\}/Re\delta x\delta z,\end{array}$$

$$\begin{array}[]{ll}F_{z}&=\int_{S}{\bf n}\cdot{\bf\Pi}_{1}\cdot{\bf i}_{z}dS\\ &=\int_{0}^{\pi}(\Pi_{1nn}\cos\theta-\Pi_{1n\tau}\sin\theta)d\theta\\ &=\int_{0}^{\pi}[\Pi_{1nn}\cos\theta-(\Pi_{2n\tau}-\frac{1}{ReCa}\frac{\partial\sigma}{\partial\theta})\sin\theta]d\theta\\ &=\int_{0}^{\pi}(\Pi_{1nn}\cos\theta-\Pi_{2n\tau}\sin\theta-\frac{1}{Re}\frac{\partial\Theta}{\partial\theta}\sin\theta)d\theta\\ &\approx-\frac{2\pi}{Re}(1+\mu_{2})V_{\infty}-\frac{1}{Re}\int_{0}^{\pi}\frac{\partial\Theta}{\partial\theta}\sin\theta d\theta\\ &\approx-\frac{2\pi}{Re}(1+\mu_{2})V_{\infty}+\frac{1}{Re}\Sigma_{i=1}^{\frac{N}{2}}|\frac{\partial\Theta}{\partial\theta}|_{t}\sin\theta_{i}\Delta\theta+\frac{1}{Re}\Sigma_{i=\frac{N}{2}+1}^{N}|\frac{\partial\Theta}{\partial\theta}|_{b}\sin\theta_{i}\Delta\theta\\ &=-\frac{2\pi}{Re}(1+\mu_{2})V_{\infty}+\frac{1}{Re}(|\frac{\Delta\Theta}{\Delta\theta}|_{t}+|\frac{\Delta\Theta}{\Delta\theta}|_{b})\Delta\theta\Sigma_{i=1}^{\frac{N}{2}}\sin\theta_{i},\end{array}$$

$$\begin{array}[]{l}\nabla\cdot{\bf v}_{i}=0,\\ \nabla p_{i}=\frac{\mu_{i}}{Re}\Delta{\bf v}_{i},\\ \Delta\Theta_{i}=0,\end{array}$$

$$(v_{1r},v_{1\theta })\to (-V_{\infty }\cos \theta ,V_{\infty }\sin \theta ),p_1\to 0,{\Theta }_1\to r\cos \theta$$

$$\begin{array}[]{l}v_{1r}(1,\theta)=v_{2r}(1,\theta)=0,\\ v_{1\theta}(1,\theta)=v_{2\theta}(1,\theta),\\ {\bf n}\cdot{\bf\Pi}_{1}-{\bf n}\cdot{\bf\Pi}_{2}=\frac{1}{ReCa}[\sigma{\bf n}-\nabla\sigma],\\ \Theta_{1}(1,\theta)=\Theta_{2}(1,\theta),\\ \frac{\partial\Theta_{1}}{\partial r}(1,\theta)=k_{2}\frac{\partial\Theta_{2}}{\partial r}(1,\theta),\end{array}$$

$$\begin{array}[]{l}\Psi_{1}(r,\theta)=-V_{\infty}(r-\frac{1}{r})\sin\theta,\\ \Psi_{2}(r,\theta)=-V_{\infty}r(r^{2}-1)\sin\theta,\\ p_{1}(r,\theta)=0,\\ p_{2}(r,\theta)=p^{\prime}_{0}-8\mu_{2}V_{\infty}r\cos\theta/Re,\\ \Theta_{1}(r,\theta)=(r+\frac{1-k_{2}}{1+k_{2}}\frac{1}{r})\cos\theta,\\ \Theta_{2}(r,\theta)=\frac{2}{1+k_{2}}r\cos\theta,\end{array}$$

$${\int }_{S_1}({\Pi }_{1n\tau }\sin \theta -{\Pi }_{1nn}\cos \theta )dS={\int }_0^{\pi }({\Pi }_{1n\tau }\sin \theta -{\Pi }_{1nn}\cos \theta ){|}_{r=1}d\theta =0.$$

$$V_{\infty }=-\frac{1}{2\pi (1+{\mu }_2)}{\int }_0^{\pi }\frac{\partial {\Theta }_1}{\partial \theta }(1,\theta )\sin \theta d\theta =\frac{1}{2(1+{\mu }_2)(1+k_2)}.$$

Full text

Keywords  Interfacial tension; Thermocapillary migration of a droplet; Large Marangoni numbers; Microgravity

1 Introduction

An immiscible droplet or bubble is placed in an ambient fluid with temperature gradient in microgravity environment, it will move in the direction of increasing temperature due to the surface tension. This phenomenon is called as thermocapillary migration in fundamental hydrodynamics and has many practical applications[1]. The pioneering work in [2] used a linear model to predict the migration speed of a droplet in the limits of zero Reynolds (Re) and Marangoni (Ma) numbers. To include inertial effects, the above analysis was extended to treat the weak nonlinear model in the range of small Ma numbers[3]. Under the quasi-steady state assumption, the analytical results were confirmed with the experimental results[4]. Since then, thermocapillary migration of a droplet in a large range of Ma numbers has been studied extensively by a series of the theoretical analyses[5], numerical simulations[6] and experimental investigations[7]. In particular, for large Ma numbers, it was reported that the migration speed of a droplet increases with increasing Ma number[8]. The theoretical results are in a qualitative agreement with the correspondent numerical simulations[9], but are qualitatively inconsistent with the experimental investigations[10,11]. In the above works for large Ma numbers, both the theoretical analyses and numerical simulations are based on the assumption of quasi-steady state. However, in the experiment investigations the droplet migrations are in an accelerating process and do not reach steady states. Therefore, although the thermocapillary droplet migration at small Ma numbers is well understood, it is unclear what kinds of states might happen in the thermocapillary droplet migration at large Ma numbers. Extension of the quasi-steady state assumption proposed at small Ma numbers to large Ma numbers is still a topic to be studied with emphasis laid on the physical mechanism. Moreover, in view of the mechanism of the varied surface tension with temperature, the thermocapillary migration of a droplet was extended to drive a film or a droplet on a substrate with a horizontal temperature gradient, such as the examples of films climbing planes by thermocapillary effect[12, 13, 14] and moving droplets on horizontal substrates[15, 16, 17, 18]. This subject called as thermocapillary actuation has the similar physical mechanism as thermocapillary droplet migration and potential applications on the chemical industry and biological engineering.

To investigate thermocapillary droplet migration in a large range of Ma numbers, steady-unsteady thermocapillary migrations of a planar droplet at moderate and large Ma numbers were observed in both the numerical and analytical studies[19], where the numerical results are in qualitative agreement with the experimental results. However, the controlling time scale of whole migration process is the thermal diffusion time scale and far greater than the convective one. It is expected that a steady state only appears in the relative time scale of the order $O(Ma)$ when the convective time scale is adopted to non-dimensionalize time. In this paper, using the dimension-analysis method, we first determine four different time scales in the thermocapillary migration system and connect them with the non-dimensional parameters of the system. Then, using the front-tracking method, we numerically study thermocapillary migration of a planar droplet at moderate and large Ma numbers and capture their terminal states in a longer time scale, which is beyond the relative thermal diffusion time scale $O(Ma)$. The physical mechanism of the steady/unsteady migration at the moderate/large Ma numbers is further analyzed by comparing the time variation of the temperature fields at the front and the rear of the droplet and determining the driving force on the droplet based on the terminal temperature distributions at the interface of the droplet.

The paper is organized as follows. In section 2, we present the governing equations of the problem of thermocapillary migration of a planar droplet and the numerical methods to solve them. In section 3, we make the dimension analysis of the above problem to obtain the controlling time scale of the whole thermocapillary migration process. In section 4, numerical results for two cases of temperature gradients $G=12$K/cm and 9 K/cm are analyzed. Finally, in section 5, some conclusions and discussions are given.

2 Governing equations and numerical methods

Consider the thermocapillary migration of a planar droplet with the radius $R_{0}$ in a continuous phase fluid of infinite extent (with the dynamical viscosity $\mu_{1}$) under a uniform temperature gradient $G$. The reference velocity is defined as

[TABLE]

where $\sigma_{\Theta}(=d\sigma/d\Theta)$ is the rate of change of surface tension with temperature. By taking $R_{0}$, $v_{o}$ and $GR_{0}$ as characteristic quantities to make coordinates, velocity and temperature dimensionless, the continuity, momentum and energy equations for the continuous phase fluid and the droplet in a laboratory coordinate system are written in the following non-dimensional form

[TABLE]

where ${\bf v}_{i}=(u_{i},v_{i})$, $p_{i}$ and $\Theta_{i}$ are velocity, pressure and temperature, respectively. The physical coefficients (density $\rho_{i}$, dynamic viscosity $\mu_{i}$, thermal conductivity $k_{i}$ and thermal diffusivity $\kappa_{i}$) are non-dimensionlized by the quantities of continuous fluid. Symbols with subscript 1 and 2 denote physical coefficients of the continuous fluid and the droplet, respectively. ${\bf f}_{\sigma}$ is the scaled surface tension acting on the interface. The Re and Ma numbers are respectively defined as

[TABLE]

As shown schematically in Fig. 1, only half of the velocity/temperature field is determined due to the mirror symmetry about $z$ axis involved in the system. The solutions of Eqs. (2) satisfy the following initial conditions of the half domain $x\in[0,x_{1}]$ and $z\in[z_{0},z_{1}]$

[TABLE]

and boundary conditions at the top and bottom walls ($z=z_{1}$ and $z=z_{0}$), on the central symmetric axis ($x=0$) and at the right wall ($x=x_{1}$)

[TABLE]

To simulate the thermocapillary droplet migration, a fixed regular staggered MAC grid in the computational domain is used. For discretizing Eqs. (2), a second-order central difference scheme for the spatial variables and an explicit predictor-corrector second-order scheme for time integration are adopted. By using the front-tracking method[20], the immiscible interface is considered to have a finite width, so that all physical coefficients across the interface are continuous. Here, a weighting function[21] is taken as

[TABLE]

where

[TABLE]

and $(x_{p},z_{p})$ is the interface node. With the updated physical coefficients, the velocity, pressure and temperature fields are computed by the Chorin’s projection method. Meanwhile, the non-dimensional surface tension[20] is written in the form of body force as

[TABLE]

where ${\bf\tau}$ is an unit tangent vector, $s$ is the arc length along the interface, $\sigma(=\sigma_{0}+\sigma_{\Theta}\Theta)$ is the surface tension coefficient and $Ca(=v_{0}\mu_{1}/\sigma_{0}$) is the Capillary number. The surface tension at the interface is distributed to the grid points by means of the weighting function (6). More details of the numerical methods were presented in [19].

3 Dimension analysis

The system of thermocapillary droplet migration, which has the five basic dimensions (mass $M$, length $L$, time $T$, temperature $K$ and quantity of heat $Q$), is governed by the seven non-dimensional parameters ($Re$, $Ma$, $Ca$, $\rho_{2}$, $\mu_{2}$, $k_{2}$ and $\kappa_{2}$). They are related to the eight independent characteristic quantities ($R_{0}$, $v_{0}$, $GR_{0}$, $\sigma_{0}$, $\rho_{1}$, $\mu_{1}$, $k_{1}$ and $\kappa_{1}$). Based on the above five basic dimensions, the eight characteristic quantities can be described as $R_{0}\sim L$, $v_{0}\sim L/T$, $GR_{0}\sim K$, $\sigma_{0}\sim MT^{-2}$, $\rho_{1}\sim ML^{-3}$, $\mu_{1}\sim ML^{-1}T^{-1}$, $k_{1}\sim QL^{-1}K^{-1}$ and $\kappa_{1}\sim L^{2}T^{-1}$. By using the dimension-analysis method, four time scales describing different physical processes are determined as the convective time scale ($T_{c}=R_{0}/v_{0}$), the momentum diffusion time scale ($T_{m}=R^{2}_{0}\rho_{1}/\mu_{1}$), the thermal diffusion time scale ($T_{t}=R^{2}_{0}/\kappa_{1}$) and the capillary action time scale ($T_{\sigma}=R_{0}\mu_{0}/\sigma_{0}$). When the convective time scale $T_{c}$ is taken as the scaled time in Sect. 2, the other three dimensionless time scales are written as $Re=T_{m}/T_{c}$, $Ma=T_{t}/T_{c}$ and $Ca=T_{\sigma}/T_{c}$, respectively. They are exactly the three non-dimensional parameters of the system.

In [19], the numerical studies on the thermocapillary droplet migration focused on the parameters of the system $0.66\leq Re\leq 53.4$, $44.7\leq Ma\leq 3622.8$ and $0.0044\leq Ca\leq 0.040$, which are regulated with the radius $R_{0}$. For each $R_{0}$, $Ma$ is the largest one in the three parameters of the system. It is clear that the thermal diffusion time scale $T_{t}$ is the largest one of the system and controls the whole migration process. So, a steady or unsteady state as the terminal one of the thermocapillary droplet migration will be expected to reach in a time scale of this order, or $t=O(Ma)$ in the dimensionless terms.

4 Results and analysis

To verify the accuracy of the numerical model in Sect. 2, we adopt very small Re and Ma numbers to determine the terminal migration velocity of the droplet and compare with the analytical solution in the limits of zero Re and Ma numbers given in the Appendix. The non-dimensional parameters are chosen as Re=0.01, Ma=0.01, Ca=0.01, $\rho_{2}=\mu_{2}=k_{2}=\kappa_{2}=0.5$. The computational domain is fixed to the size $6\times 12$. Based on $72\times 144$, $144\times 288$, $288\times 576$ and $432\times 864$ grid points, i.e.,12, 24, 48 and 72 grid points per droplet radius, migration velocities of the droplet are plotted in Fig. 2(a). For each grid resolution, the migration velocity of the droplet increased from zero reaches a steady value at last and approximates to the analytical result $V_{\infty}=0.22$. A convergent trend is found when increasing the grid resolution in the simulations. In Fig. 2(b), the migration velocities of the droplet at four Re(= 0.005,0.01,0.05 and 0.1), Ma=0.01, Ca=0.01, $\rho_{2}=\mu_{2}=k_{2}=\kappa_{2}=0.5$ with the grid resolution for 48 grid points per droplet radius exhibit a convergent approximation to the analytical result with an error (about 9%) when decreasing Re numbers. In the following calculations, we fix 48 grid points per droplet radius as the grid resolution. It is expected that thermal boundary layers around the droplet surface at large Ma numbers have a thickness of $O(Ma^{-1/2})$[22]. So the above grid resolution can be used to depict the thermocapillary droplet migration in a large range of Ma numbers.

The silicone oil of nominal viscosity 5cst and the FC-75 Fluorinert liquid, i.e., the working media in the space experiments[11], are adopted as the continuous phase fluid and the droplet, respectively. The physical parameters of the continuous fluid and the droplet at temperature $25^{o}$C are given in Table I. $\sigma_{T}$ is fixed as -0.044 dyn/cmK[11] and $\sigma_{0}\approx$ 6 dyn/cm[23] is adopted. From the values of the continuous fluid parameters, Prantdl number ($Pr=Ma/Re=\mu_{1}/\rho_{1}\kappa_{1}$) is determined as 67.8. The computational domain is chosen as $\{x,z\}\in\{[0,4],[-4,44]\}$. The initial droplet is placed at the position (0,0) and the integration time step is varied in $1\times 10^{-5}-2\times 10^{-4}$ depending on the migration velocity. To simulate the migration process in the longer time scale $O(Ma)$ at moderate and large Ma numbers, we regulate $R_{0}$ to approach the given Ma number.

4.1 Flow field with the temperature gradient G=12K/cm

To simulate the migration process with $G$= 12 K/cm, the droplets with $R_{0}=0.05$cm, 0.075cm and 0.10cm are taken to make the systematic parameters, which are presented in Table II and hereinafter referred to as Ma numbers. For $Ma\leq 178.9$, the thickness of thermal boundary layers is not smaller than 1/13.4. So the above grid resolution is sufficiently high to describe the thermal boundary layers for $Ma\leq 178.9$.

Figure 3(a) displays the time evolution of droplet migration velocities for Ma=44.7, 100.6 and 178.9. The maximal computational time is $T_{max}=220$. As given in Sect. 3, the controlling time scale of whole migration process is the relative thermal diffusion time scale $T_{t}\sim O(Ma)$. Since $T_{max}$ is beyond $T_{t}$, all physical processes including the thermal diffusion in the system are fully developed at $T_{max}$. The whole migration processes exhibit three stages based on the curve features . At the initial one($t\leq 3$), the migration velocities markedly increase from zero to local maximums. Then they have decreasing-increasing processes in the range $3<t\leq 70$. After the oscillating processes, the migration velocities trend to two different qualitative behaviors in the range $70<t\leq T_{max}$ depending on the Ma numbers. One is the approaching a steady value (for Ma=44.7), the other two are the monotonically increasing with time and do not reach any constants (for Ma=100.6 and 178.9). The increasing slope for Ma=178.9 is larger than that for Ma=100.6. To further depict the above migration processes, time evolution of temperature difference $\Delta\Theta[=\Theta(1,0)-\Theta(1,\pi)$] between the front ($\theta=0$) and the rear ($\theta=\pi$) of the droplet is shown in Fig. 3(b). $\Delta\Theta$ initially drops from 2 and then has an oscillating process. At last, it approaches a constant for Ma=44.7 and monotonically increases with time for Ma=100.6 and 178.9. The stationary temperature difference satisfies the requirement of steady thermocapillary migration of the droplet. However, the increasing temperature difference deviates from the requirement of steady thermocapillary migration of the droplet. Thus, it is concluded that the thermocapillary droplet migration at the moderate/large Ma numbers can/cannot reach a steady state and thus in a steady/unsteady process.

In Fig. 4(a), terminal streamlines in a reference frame moving with the droplets at $T_{max}$ for Ma=44.7, 100.6 and 178.9 are shown. The droplets are located at the positions $z_{c}=$24.7, 24.6 and 29.2, respectively. Two vortices symmetric about the $z$-axis appear within the droplet. For each Ma number, the external streamlines go around the droplet and the flow over the droplet does not separate. Although the terminal streamlines evolved in the range of Ma numbers correspond to both the steady and unsteady states, their patterns are similar and almost independent of Ma numbers. These properties reveal that the instantaneous velocity fields in the unsteady migration processes at large Ma numbers may have the forms of the steady ones at moderate Ma numbers.

Figure 4(b) displays terminal isotherms in the laboratory coordinate frame at $T_{max}$ for Ma=44.7, 100.6 and 178.9. For each Ma number, the whole temperature field is divided into external and internal domains of the droplet. On the one hand, the isotherms in the external domain exhibit three kinds of characteristics in the different areas. The first one is the unperturbed uniformly-spaced parallel lines above the droplet. The second one is the bending isotherms near the droplet to the contrary of the migration direction. More and more isotherms are concentrated to the interface of the droplet as Ma number increases, so that the thermal boundary layer is formed near the interface. The last one is the curved temperature field with a gradient below the droplet. In the thermal boundary layer, there appear two regions related to the thermal transfer across the interface. One is the most of interface with $\frac{\partial\Theta_{i}}{\partial r}(1,\theta)>0$, the other is a local part near the rear of the droplet with $\frac{\partial\Theta_{i}}{\partial r}(1,\theta)<0$. The thermal flux across the interface in the former/latter one goes from the outside/inside of the droplet to the inside/outside. The former/latter one is larger/smaller as Ma number increases. So, for large Ma numbers, the thermal flux across the interface at the terminal states almost goes from the outside of the droplet to the inside. On the other hand, the isotherms in the internal domain of the droplet exhibit different kinds of the patterns depending on Ma numbers. For Ma=44.7, the pattern of the isotherms has a small cap-type when the droplet migrates in the steady process. In this case, the temperature in the small cap-type isotherm ($\Theta_{2}=15.8$) is the lowest and about 2.8 lower than the front temperature. For Ma=100.6, the pea-type isotherms as the terminal state are formed in the unsteady migration process. At this time, the temperature in the pea-type isotherm ($\Theta_{2}=12.8$) is the lowest and about 4.4 lower than the front temperature. For Ma=178.9, the isotherms with two vortices are generated when the droplet is in the accelerating process. The lowest temperature within the droplet ($\Theta_{2}=11$) appears in the center of the vortex and is about 9.7 lower than the front temperature. In general, the thermal convection and conduction are two ways of heat transfer in the system, but the droplet can only obtain the thermal energy though the thermal conduction across the interface. For moderate Ma numbers, the thermal convection around the droplet and the thermal conduction across the droplet are equally important to heat transfer in the system, so that the external and internal temperature of the droplet may have a fixed difference. This leads to a steady migration process. For large Ma numbers, the thermal convection around the droplet is a stronger way of heat transfer in the system than the thermal conduction across the droplet. The internal temperature of the droplet is more difficult to increase than the external one of the droplet during the migration of droplet. In the whole migration process, although the internal temperature of the droplet increases, its increment is far lower than that of the external temperature. The increasing difference of temperature leads to an unsteady migration precess.

Figure 4(c) displays the relative temperature distributions $\Delta\Theta/Re=[\Theta(1,\theta)-\Theta(1,\pi)]/Re$ along the interface from the front ($\theta=0$) to the rear ($\theta=\pi$) of the droplet for the terminal states shown in Fig. 4(b). In the range of Ma numbers, $\Delta\Theta/Re$ monotonously decreases as $\theta$ increases, i.e., $\frac{\partial\Delta\Theta}{\partial\theta}/Re=\frac{\partial\Theta}{\partial\theta}/Re<0$. Although $-\frac{\partial\Theta}{\partial\theta}/Re>0$ holds on the whole surface, its averaged values $|\frac{\Delta\Theta}{\Delta\theta}|_{t}/Re\approx[\Theta(1,0)-\Theta(1,\pi/2)]/(\pi Re/2)$ on the top surface {$\theta\in[0,\frac{\pi}{2}]$} and $|\frac{\Delta\Theta}{\Delta\theta}|_{b}/Re\approx[\Theta(1,\pi/2)-\Theta(1,\pi)]/(\pi Re/2)$ on the bottom surface {$\theta\in(\frac{\pi}{2},\pi]$} are varied depending on Ma numbers. For Ma=44.7, $|\frac{\Delta\Theta}{\Delta\theta}|_{t}/Re<|\frac{\Delta\Theta}{\Delta\theta}|_{b}/Re$. For Ma=100.6, the above relation holds, but their difference decreases. For Ma=178.9, $|\frac{\Delta\Theta}{\Delta\theta}|_{t}/Re\approx|\frac{\Delta\Theta}{\Delta\theta}|_{b}/Re$, and an inflection point ($\frac{\partial^{2}\Theta}{\partial\theta^{2}}=0$) appears near $\theta=\pi/2$. The net force acting on the droplet exerted by the continuous phase fluid in the vertical direction[1] is written as

[TABLE]

where ${\bf n}$ is the unit normal to the surface, ${\bf\Pi}_{1}=-p_{1}{\bf I}+\frac{1}{Re}[\nabla{\bf v}_{1}+(\nabla{\bf v}_{1})^{T}]$ is the stress tensor of the continuous phase and $V_{\infty}$ is the instantaneous migration velocity of the droplet. In the above derivation, the tangential stress balance $\Pi_{1n\tau}-\Pi_{2n\tau}=-\frac{1}{ReCa}\frac{\partial\sigma}{\partial\theta}$ at the interface of the two-phase fluids is used. Due to the above investigations in the steady and unsteady migration processes, the instantaneous velocity fields of two-phase fluids are assumed as those in the limits of zero Re and zero Ma numbers given in the Appendix. This assumption is also supported with the same velocity fields at the zero-order approximation of thermocapillary migration of a spherical droplet at zero Ma(Re), small Ma(Re) and large Ma(Re) numbers[2, 3, 8].

In Eq. (9), the net force includes of the first term for the drag force and the second term for the driving force. When the driving force is larger than the drag force, i.e., $F_{z}>0$, the droplet migrates in an accelerating process. When the driving force is equal to the drag force, i.e., $F_{z}=0$, the droplet is in a steady state. The surface tensions of the top and bottom surfaces of the droplet simultaneously contribute to the driving force. However, for moderate Ma numbers, the surface tension of the bottom surface, which is larger than that of the top surface, is the main component of the driving force. For large Ma numbers, the surface tensions of the top and bottom surfaces of the droplet have the same order of magnitude in the driving force when the inflection point appears. The driving force on the droplet decreases as Ma number increases.

4.2 Flow field with the temperature gradient G = 9 K/cm

To simulate the migration process with $G$=9K/cm, the droplets with $R_{0}=0.075$cm, 0.10cm and 0.125cm are taken to make the systematic parameters, which are presented in Table III. For $Ma\leq 209.7$, the thickness of thermal boundary layer is not smaller than 1/14.5. So the above grid resolution is sufficiently high to describe the thermal boundary layers for $Ma\leq 209.7$.

Figure 5(a) displays the time evolution of droplet migration velocities for Ma=75.5, 134.2 and 209.7. The maximal computational time $T_{max}=220$ is beyond the relative thermal diffusion time scale $T_{t}\sim O(Ma)$. The initial migration processes exhibit increasing-decreasing-increasing processes in the range $0<t\leq 100$. Then the migration velocities trend to two different qualitative behaviors in the range $100<t\leq T_{max}$ depending on Ma numbers. One is the approaching a steady value (for Ma=75.5), the other two are the monotonically increasing with time and do not reach any constants (for Ma=134.2 and 209.7). The increasing slope for Ma=209.7 is larger than that for Ma=134.2. These migration processes can be understood in the time evolution of temperature difference $\Delta\Theta$ between the front and the rear of the droplet shown in Fig. 5(b). For Ma=75.5, the stationary temperature difference satisfies the requirement of steady thermocapillary migration of the droplet. However, for Ma=134.2 and 209.7, the monotonically increasing temperature difference deviates from the requirement of steady thermocapillary migration of the droplet. Thus, it can be concluded that the thermocapillary droplet migration at moderate/large Ma numbers can/cannot reach a steady state and thus in a steady/unsteady process.

In Fig. 6(a), terminal streamlines in a reference frame moving with the droplets at $T_{max}$ for Ma=75.5, 134.2 and 209.7 are shown. The droplets are located at the positions $z_{c}=$24.2, 26.3 and 34.3, respectively. The terminal streamlines evolved in the range of Ma numbers reveal that the instantaneous velocity fields are similar no matter they are the steady states or the unsteady ones. Figure 6(b) displays terminal isotherms in the laboratory coordinate frame at $T_{max}$ for Ma=75.5, 134.2 and 209.7. For each Ma number, the isotherms in the external domain of the droplet bend to the contrary of the migration direction and converge near the interface of the droplet to form a thermal boundary layer. In the thermal boundary layer, there appear two regions with $\frac{\partial\Theta_{i}}{\partial r}(1,\theta)>0$ for the most of interface and $\frac{\partial\Theta_{i}}{\partial r}(1,\theta)<0$ for a local part near the rear of the droplet. They exhibit two thermal fluxes with different directions across the interface and evolve as Ma number increases. For large Ma numbers, the thermal flux across the interface at the terminal states almost goes from the outside of the droplet to the inside. Meanwhile, the isotherms in the internal domain of the droplet exhibit different kinds of patterns depending on Ma numbers. For Ma=75.5, the pattern of the isotherms has a small cap-type when the droplet migrates in the steady process. The lowest temperature within the droplet is $\Theta_{2}=13.8$ in the small cap-type isotherm and about 3.5 lower than the front temperature. The constant temperature difference $\Delta\Theta=3.5$ between the outside and inside of the droplet implies a steady migration precess. It is based on the fact that the thermal convection and conduction are two equally important to heat transfer in the system at moderate Ma numbers. For Ma=134.2, the pea-type isotherms are formed in the unsteady migration process. The lowest temperature within the droplet is $\Theta_{2}=12$ in the pea-type isotherm and about 6.2 lower than the front temperature. For Ma=209.7, the isotherms with two vortices are generated when the droplet is in the accelerating process. The lowest temperature within the droplet is $\Theta_{2}=11.2$ in the center of the vortex and about 15.0 lower than the front temperature. The temperature difference $\Delta\Theta$ between the outside and inside of the droplet increases from 6.2 to 15 when Ma number increases from 134.2 to 209.7. This leads to an unsteady migration precess. The reason is that the thermal convection around the droplet is a stronger way of the heat transfer in the system than the thermal conduction across the droplet at large Ma numbers. As shown in Fig. 6(c), $\Delta\Theta/Re$ in the whole surface monotonously decreases as $\theta$ increases in the range of Ma numbers. For Ma=75.5, $|\frac{\Delta\Theta}{\Delta\theta}|_{t}/Re<|\frac{\Delta\Theta}{\Delta\theta}|_{b}/Re$. For Ma=134.2, the above relation holds, but their difference decreases. The surface tension of the bottom surface, which is larger than that of the top surface, is the main component of the driving force in Eq. (9). For Ma=209.7, $|\frac{\Delta\Theta}{\Delta\theta}|_{t}/Re>|\frac{\Delta\Theta}{\Delta\theta}|_{b}/Re$, and an inflection point ($\frac{\partial^{2}\Theta}{\partial\theta^{2}}=0$) appears near $\theta=\pi/2$. The surface tension of the top surface, which is larger than that of the bottom surface, is the main component of the driving force in Eq. (9) after the inflection points appears. The driving force on the droplet decreases first and then increases as Ma number increases. From Figs. 4(c) and 6(c), it can be concluded that the appearance of the inflection point implies that the driving force will increase as Ma number increases.

5 Conclusions and discussions

To sum up, thermocapillary migration of a planar droplet at moderate and large Ma numbers has been investigated analytically and numerically. Firstly, by using the dimension-analysis method, the thermal diffusion time scale has been determined as the controlling one of the thermocapillary droplet migration system. During this time, the whole thermocapillary migration process is fully developed. Then, by the aid of the front-tracking method, the numerical simulations in a longer time scale, which is beyond the thermal diffusion time scale, have exhibited that the terminal states of thermocapillary droplet migration at moderate and large Ma numbers are steady and unsteady, respectively. In the terminal states, the instantaneous velocity fields in the unsteady migration processes at large Ma numbers have the forms of the steady ones at moderate Ma numbers. However, in view of the former instantaneous temperature fields, the surface tension of the top surface of the droplet gradually becomes the main component of the driving force on the droplet after the inflection point appears. It is different from that the surface tension of the bottom surface of the droplet is the main component of the driving force on the droplet for the latter ones. The analysis based on the time evolution of velocity and temperature fields implies that the temperature difference between the front and rear of the droplet is a constant for moderate Ma numbers and has the increasing tread for large Ma numbers. It satisfies/deviates from the requirement of steady thermocapillary migration of the droplet at moderate/large Ma numbers. These phenomena originate from the evolution of the heat transfer across/around the droplet in the system depending on Ma numbers. The significance of the thermal convection around the droplet is higher than/just as the thermal conduction across the droplet at large/moderate Ma numbers.

Acknowledgments This work was partially supported by the National Science Foundation through the Grants No. 11172310 and No. 11472284. The author thanks the research computing facility of the Institute of Mechanics of the Chinese Academy of Sciences for assisting in the computation.

Appendix: A general solution of steady thermocapillary migration of a planar droplet in the limits of zero Re and zero Ma numbers

The steady incompressible continuous, momentum and energy equations for the continuous phase fluid and the droplet in the limits of zero Re and Ma numbers can be written in the following dimensionless form in a polar coordinate system ($r,\theta$) moving with the droplet velocity $V_{\infty}$

[TABLE]

where ${\bf v}_{i}=(v_{ir},v_{i\theta})$. The boundary conditions are written in the form of dimensionless as

[TABLE]

at infinity and

[TABLE]

at the interface of the two-phase fluids, ${\bf n}$ is the unit vector normal to the interface, ${\bf\Pi}_{i}=-p_{i}{\bf I}+\frac{\mu_{i}}{Re}[\nabla{\bf v}_{i}+(\nabla{\bf v}_{i})^{T}]$ is the stress tensors of the two-phase fluids. Following the methods for solving the linear model[2, 24], the solutions of the equations (10) with the boundary conditions (11)(12) can be determined as

[TABLE]

where the velocity fields $(v_{ir},v_{i\theta})=(\frac{\partial\Psi_{i}}{r\partial\theta},-\frac{\partial\Psi_{i}}{\partial r})$ are written in terms of the stream functions $\Psi_{i}(r,\theta)$ of the two-phase fluids. The pressure fields $p_{i}(r,\theta)$ with a constant $p^{\prime}_{0}$ are obtained by integrating the momentum equations in Eqs. (10).

In addition to the above boundary conditions (12) at the interface, the steady thermocapillary droplet migration requires that the total net force acting on the droplet is zero. In particular, the zero net force in the vertical direction[1] is expressed as

[TABLE]

When the shear stress $\Pi_{1r\theta}$ for the continuous fluid is replaced by using the stress boundary condition in Eqs. (12), the steady droplet migration speed is derived as

[TABLE]