Modeling and simulation of thermocapillary flows using lattice Boltzmann method

Journal of Computational Physics 231 (2012) 4433–4453

Contents lists available at SciVerse ScienceDirect

Journal of Computational Physics journal homepage: www.elsevier.com/locate/jcp

Modeling and simulation of thermocapillary flows using lattice Boltzmann method Haihu Liu a, Yonghao Zhang b,⇑, Albert J. Valocchi a a b

Department of Civil & Environmental Engineering, University of Illinois at Urbana-Champagin, IL 61801, USA Department of Mechanical & Aerospace Engineering, University of Strathclyde, Glasgow G1 1XJ, UK

a r t i c l e

i n f o

Article history: Received 4 December 2011 Accepted 15 February 2012 Available online 25 February 2012 Keywords: Microfluidics Lattice Boltzmann model Thermocapillary Marangoni stresses Droplet manipulation

a b s t r a c t To understand how thermocapillary forces manipulate droplet motion in microfluidic channels, we develop a lattice Boltzmann (LB) multiphase model to simulate thermocapillary flows. The complex hydrodynamic interactions are described by an improved color-fluid LB model, in which the interfacial tension forces and the Marangoni stresses are modeled in a consistent manner using the concept of a continuum surface force. An additional convection–diffusion equation is solved in the LB framework to obtain the temperature field, which is coupled to the interfacial tension through an equation of state. A stress-free boundary condition is also introduced to treat outflow boundary, which can conserve the total mass of an incompressible system, thus improving the numerical stability for creeping flows. The model is firstly validated against the analytical solutions for the thermocapillary driven convection in two superimposed fluids at negligibly small Reynolds and Marangoni numbers. It is then used to simulate thermocapillary migration of three-dimensional deformable droplet at various Marangoni numbers, and its accuracy is once again verified against the theoretical prediction in the limit of zero Marangoni number. Finally, we numerically investigate how the localized heating from a laser can block the microfluidic droplet motion through the induced thermocapillary forces. The droplet motion can be completely blocked provided that the intensity of laser exceeds a threshold value, below which the droplet motion successively undergoes four stages: constant velocity, deceleration, acceleration, and constant velocity. When the droplet motion is completely blocked, four steady vortices are clearly visible, and the droplet is fully filled by two internal vortices. The external vortices diminish when the intensity of laser increases. Ó 2012 Elsevier Inc. All rights reserved.

1. Introduction Microfluidic devices consume small amount of materials and energy, and provide rapid analysis and diagnostics. Microfluidics therefore has potential to develop highly integrated and automated chemical and biological processes. Recently, droplet-based microfluidics has emerged as a promising flexible platform for microfluidic functions. As samples/reagents are confined in the droplets, it can avoid sample dilution caused by Taylor dispersion, and eliminate surface adsorption and cross sample contamination [1]. Numerous droplet-based applications including polymerase chain reaction [2,3], enzyme kinetic assays [4–6], protein crystallization [7,8], and synthesis of organic molecules or nanoparticles [9,10] have been demonstrated. Furthermore, simple Boolean logic functions have been performed in droplet microfluidic systems [11,12], a critical step towards the realization of a microfluidic computer. ⇑ Corresponding author. E-mail addresses: [email protected] (H. Liu), [email protected] (Y. Zhang), [email protected] (A.J. Valocchi). 0021-9991/$ - see front matter Ó 2012 Elsevier Inc. All rights reserved. doi:10.1016/j.jcp.2012.02.015

4434

H. Liu et al. / Journal of Computational Physics 231 (2012) 4433–4453

The manipulation of droplets with high precision and flexibility is essential in droplet-based microfluidic applications. The most common droplet manipulations are droplet formation, fission, fusion, mixing and sorting [13]. Various approaches have been employed to achieve these manipulations based on the sources of applied force, e.g., hydrodynamic stress [14–16], electrohydrodynamic force [17–19], thermocapillary force [20,21], surface acoustic wave [22], magnetic force [23,24], and optical forces [25,26]. Among these approaches, the optical forces are very attractive and may provide a promising route towards microactuation because light fields can be easily focused to micrometric spots, and they are contactless and dynamically reconfigurable, without any additional requirement on micro-chip fabrication [27]. However, optical forces are within the pN range, which are usually too weak to effectively manipulate fast moving droplets. In the recent years, an optical alternative based on the production of localized thermocapillary stresses (also known as the Marangoni effect [28]) at the droplet interface has been demonstrated. In the presence of surfactant, the magnitude of thermocapillary forces induced by the localized heating from a laser can reach lN range [29], so they are more effective than optical forces to counteract hydrodynamic forces acting on the moving droplets. Thermocapillary convection is a phenomenon of fluid movement induced by the interfacial tension difference due to the imposed temperature gradient at the interface. The thermocapillary convection was first investigated by Young et al. [30] who determined the terminal velocity of a spherical droplet in an infinite ambient fluid in the creeping flow limit when a linear temperature profile was imposed. Since then, many studies have been reported, but most of them have neglected deformation of the droplet [31–36]. The main motivation for these studies stems from microgravity applications where sedimentation and gravity-driven convection are largely eliminated. Thermocapillary convection is particularly suitable for driving flows at small scales where the interface effects dominate. Optically induced thermocapillary forces have now been used for manipulating droplets in microfluidic devices. It was observed that a floating droplet [37] and a non-wetting droplet deposited on a substrate [38] can move under laser-driven thermocapillary stresses. Baroud et al. [29] showed that a water droplet, carried by an immiscible oil flow in a microchannel, may be blocked when subjected to a focused laser spot. This blocking indicates that the drag force is balanced by the thermocapillary force, which is generated by the laser heating. This thermocapillary force was also combined with the microfluidic geometries to realize various droplet manipulations including mixing, sorting, fission, fusion, sampling and switching [39,40]. Quantitative understanding of thermocapillary flows requires detailed information of flowfield such as temperature and velocity distribution, as well as interfacial tension gradient along the interface. As experimental measurements of flowfield of droplet in microchannels are still formidable, we aim to develop numerical model to tackle this problem. However, numerical simulation of thermocapillary flows in microfluidic devices is a daunting task, where the capillary effect usually plays a dominant role. Discretization errors in calculation of interfacial forces may generate unphysical spurious velocities which can cripple the velocity field in the whole computational domain. Minimizing the spurious velocities at the interface still remains a major challenge for numerical models and algorithms. Due to the strong dependence of interfacial tension on temperature, the temperature fluctuations result in non-uniform interfacial tension forces and Marangoni stresses that affect the flowfield at the interface, which in turn alter the interfacial temperature distribution through the induced interfacial flows. While front tracking methods are not suitable for droplet breakup and coalescence, interface capturing methods such as volume-of-fluid and level set methods will suffer from numerical instability at the interface region when the interfacial tension becomes a dominant factor in microdroplet behavior [41]. Recently developed lattice Boltzmann (LB) method has shown great potential in modeling interfacial interactions while incorporating fluid flow as a system feature [42]. It is a pseudo-molecular method tracking evolution of the distribution function of an assembly of molecules and built upon microscopic models and mesoscopic kinetic equations [43]. Its mesoscopic nature can provide many advantages of molecular dynamics, making the LB method particularly useful for simulation of complex interfacial flows [44–48]. In this work, a multiphase lattice Boltzmann model is proposed to simulate thermocapillary flows. An improved colorfluid model is employed to simulate the flowfield and capture the phase interface, in which the interfacial forces, including the interfacial tension force and Marangoni stress, are modeled using the concept of a continuum surface force [49], and a recoloring algorithm proposed by Latva-Kokko and Rothman [50] is applied for phase segregation. An additional convection– diffusion equation is also solved in the LB framework to obtain temperature, which is coupled to the interfacial tension by an equation of state. We also introduce a stress-free boundary condition to treat the outflow boundary, which can conserve the total mass of an incompressible system, thus improving the numerical stability for creeping flows. After the capability and accuracy of this model are examined, it is used to simulate the thermocapillary migration of 3D deformable droplet at various Marangoni numbers. Finally, the droplet dynamical behavior is simulated for a droplet transported by a carrier fluid in a microchannel when a laser heating is applied. We investigate how the laser light can block the droplet motion and the effect of its intensity on the droplet motion. To the best of our knowledge, the present model is the first LB thermocapillary model, which is able to numerically simulate the dynamical blocking process of a microfluidic droplet subjected to the localized heating source from a laser. This study can provide useful information for understanding the thermocapillary flows in confined microfluidic devices and facilitate design to precisely control droplet behavior.

2. Numerical method In the LB method, a fluid is modeled as pseudo particles, whose distribution function fi is governed by the following LB equation with the Bhatnagar–Gross–Krook (BGK) collision operator [42]:

H. Liu et al. / Journal of Computational Physics 231 (2012) 4433–4453

fi ðx þ ei dt ; t þ dt Þ ¼ fi ðx; tÞ 

1

s

 fi ðx; tÞ  fieq ðx; tÞ ;

4435

ð1Þ

where fi(x, t) is the particle distribution function in the ith velocity direction at position x and time t, ei is the lattice velocity in the ith direction, s is the dimensionless relaxation time, and fieq is the equilibrium distribution function of fi given as a function of local density q and velocity u:

"

fieq

# ei  u ðei  uÞ2 u  u ; ¼ wi q 1 þ 2 þ  cs 2c2s 2c4s

ð2Þ

pffiffiffi where cs denotes the speed of sound which is given by c= 3 with c = dx/dt being the lattice speed and dx being the lattice length. In this work, we use D2Q9 and D3Q19 models for the 2D and 3D simulations, respectively. The lattice velocity ei and the weight coefficient wi for these two models are given as follows: D2Q9:

8 i ¼ 0; > < ð0; 0Þ; ei ¼ ð1; 0Þc; ð0; 1Þc; i ¼ 1; 2; 3; 4; > : ð1; 1Þc; i ¼ 5; 6; 7; 8;

ð3Þ

8 i ¼ 0; > < 4=9; wi ¼ 1=9; i ¼ 1; 2; 3; 4; > : 1=36; i ¼ 5; 6; 7; 8:

ð4Þ

D3Q19:

8 i ¼ 0; > < ð0; 0; 0Þ; ei ¼ ð1; 0; 0Þc; ð0; 1; 0Þc; ð0; 0; 1Þc; i ¼ 1; 2; . . . ; 6; > : ð1; ; 0Þc; ð1; 0; 1Þc; ð0; 1; 1Þc; i ¼ 7; 8; . . . ; 18;

ð5Þ

8 i ¼ 0; > < 1=3; wi ¼ 1=18; i ¼ 1; 2; . . . ; 6; > : 1=36; i ¼ 7; 8; . . . ; 18:

ð6Þ

The macroscopic local density and the momentum are related to the particle distribution function fi by

q¼

X i

fi ¼

X

fieq ;

i

qu ¼

X i

f i ei ¼

X

fieq ei :

ð7Þ

i

Using the Chapman–Enskog expansion for D2Q9 or D3Q19 model, Eq. (1) can lead to the Navier–Stokes equations in the long-wavelength and low-frequency limit [51]:

@ t q þ r  ðquÞ ¼ 0;

ð8Þ

@ t ðquÞ þ r  ðquuÞ  rp þ r  ½qmðru þ ruT Þ;

ð9Þ

where the pressure and the kinematic viscosity are given as p ¼ qc2s and m ¼ ðs  1=2Þc2s dt . Currently, the most applied LB multiphase methods are the color-fluid model [52], the pseudo-potential model [53], and the free-energy model [54,55]. In the color-fluid model, the procedure of redistribution of the colored density at each node to separate different phases requires time-consuming calculations of local maxima, and the perturbation step with redistribution of the colored distribution functions causes an anisotropic interfacial tension that induces high spurious velocities near interface [51]. The pseudo-potential model introduces the nearest-neighboring interaction between fluid particles to describe the intermolecular potential, and phase separation occurs with a properly chosen potential function. Although significant advances have recently made [56–58], further improvements are necessary for the pseudo-potential model to minimize spurious velocities at interface and control numerical instability for the flows with low capillary number. The free-energy model suffers from the lack of Galilean invariance [42], although the local conservation of mass and momentum is satisfied. In addition, small droplets are prone to dissolve since the multiphase system is always evolving towards the direction of minimal free energy [59]. Based on the original color-fluid model of Gunstensen et al. [52], the recent improvements have been made by Lishchuk et al. [60] as well as Latva-Kokko and Rothman [50], which enable simulational access to the flow regimes of low Reynolds and capillary numbers. However, the current LB color-fluid models only consider the effect of a constant interfacial tension, so they cannot simulate the thermocapillary flows around moving and stationary droplets, which are important for many microgravity applications and opto-thermal manipulation of droplets in microfluidic devices.

4436

H. Liu et al. / Journal of Computational Physics 231 (2012) 4433–4453

In this study, we will present a thermo-color-fluid model in the LB framework, which can accurately describe thermocapillary flows in an interfacial tension dominated multiphase system. The original color-fluid model of Gunstensen et al. [52] introduced red and blue particle distribution functions fiR and fiB to represent two different fluids. The total particle distribution function is then defined as: fi ¼ fiR þ fiB . The mixed fluid undergoes the collision and streaming operations as

fi ðx þ ei dt ; t þ dt Þ ¼ fi ðx; tÞ 

fi ðx; tÞ  fieq ðx; tÞ

s

þ Ui ;

ð10Þ

where Ui is an additional collision operator (also known as the perturbation step), which contributes to the mixed interfacial regions and generates an interfacial tension:

Ui ¼ AjGj cos 2ðhi  hf Þ;

ð11Þ

where A is a free parameter controlling the interfacial tension, hi is the polar angle of the lattice vector ei. And hf is the polar angle of the local color gradient G, which is calculated as

X

Gðx; tÞ ¼

wi ei ½qR ðx þ ei ; tÞ  qB ðx þ ei ; tÞ;

ð12Þ

i

where qR and qB are the nodal densities of red and blue fluids and defined by

qR ðx; tÞ ¼

X i

fiR ðx; tÞ;

qB ðx; tÞ ¼

X

fiB ðx; tÞ:

ð13Þ

i

To account for unequal viscosities of both fluids, the relaxation time s in Eq. (10) can be defined as



m¼ s

 1 2 qR qB c dt ¼ m þ m; 2 s qR þ qB R qR þ qB B

ð14Þ

where mR and mB are viscosities of the red and blue fluids, respectively. To prevent numerical diffusion of the two fluids across the interface, the so-called recoloring step is applied, which can keep the interface sharp and prevent the two fluids from mixing with each other. The colors are demixed by maximizing the work done against the color flux q(x, t), which is defined by

qðx; tÞ ¼

X

  wi ei fiR ðx; tÞ  fiB ðx; tÞ :

ð15Þ

i

As pointed out above, the perturbation step can introduce anisotropy and high spurious velocities near the interface. Additionally, for creeping flows, the recoloring step can produce lattice pinning [50], a phenomenon that the interface is pinned or attached to the simulation lattices, rendering an effective loss of Galilean invariance. It was also demonstrated that there is an increasing tendency for lattice pinning as the capillary and Reynolds numbers decrease [61]. Lishchuk et al. [60] used the concept of a continuum surface force (CSF) [49] to model a constant interfacial tension. In their algorithm, the perturbation step in the original Gunstensen model was replaced by a direct forcing term in the mixed region. In order to satisfy the stress boundary condition and the continuity of velocity, a local pressure gradient is forced throughout the interface as an additional body force, which is incorporated into the LB model as the forcing term Ui(x, t) (which is the additional collision operator in original Gunstensen model) at the collision step. It has been reported that this algorithm can greatly reduce the spurious currents and improve the isotropy of the interface. Extending the algorithm of Lishchuk et al. [60], we present a mathematical treatment responsible for the force contribution related to the interfacial tension gradient along the interface, i.e., Marangoni stress. It is well-known that the expression for the stress jump across the interface is given by

TR  n  TB  n ¼ r  ½rðI  n  nÞ;

ð16Þ T

where I is the second-order identity tensor, T = pI + qm(ru + ru ) is the stress tensor with the superscript ‘‘R’’ (‘‘B’’) denoting the red (blue) side of the interface, r is the local interfacial tension parameter, and n is the interfacial unit vector in the normal direction. In order to induce the local stress jump across the interface, a volume-distributed interfacial force F(x, t) given by Eq. (16), should be added in the momentum equation as an additional body force. The interfacial force is

Fðx; tÞ ¼ r  ½rðI  n  nÞdR  ¼ rjdR n þ rS rdR ;

ð17Þ

where rS = (I  n  n)  r is the surface gradient operator, dR is the Dirac delta function used to localize the force explicitly at the interface, and j is the local interface curvature, which is calculated by

j ¼ rS  n:

ð18Þ

The first term on the right-hand side of Eq. (17) is the interfacial tension force and the second term is the Marangoni stress. The unit normal vector and the Dirac delta function are determined by the gradient of the phase field rqN,

H. Liu et al. / Journal of Computational Physics 231 (2012) 4433–4453

n¼

rqN jrqN j

dR ¼

;

1 jrqN j; 2

4437

ð19Þ

where the phase field qN is defined as

qN ðx; tÞ ¼

qR ðx; tÞ  qB ðx; tÞ ; 1 6 qN 6 1: qR ðx; tÞ þ qB ðx; tÞ

ð20Þ

In a thermocapillary flow, an equation of state is required to relate the interfacial tension to the temperature, which may be linear or nonlinear. For the sake of simplicity, we only consider a linear relation between the interfacial tension and the temperature in this study, i.e.,

rðTÞ ¼ rref þ rT ðT  T ref Þ;

ð21Þ

where Tref is the reference temperature, rref is the interfacial tension at Tref, rT is the rate of change of interfacial tension with temperature, defined as rT = or/@T. Substituting Eqs. (19) and (21) into Eq. (17), we obtain the interfacial force as

Fðx; tÞ ¼

1 1 rjrqN þ rT jrqN jðI  n  nÞ  rT: 2 2

ð22Þ

The forcing term Ui in Eq. (10) can be calculated incorporating the spatially varying interfacial force as [62]



Ui ¼ 1 

   1 ei  u ei  u wi  Fðx; tÞdt ; þ e i 2s c2s c4s

ð23Þ

where the fluid momentum is given by:

quðx; tÞ ¼

X i

1 fi ðx; tÞei þ Fdt : 2

ð24Þ

Using the Chapman–Enskog multiscale analysis, Eq. (10) can recover the following momentum equation in the continuum limit with the forcing term and fluid velocity given by Eqs. (23) and (24) [62]:



@ t ðquÞ þ r  ðquuÞ ¼ r qc2s þ r  ½qmðru þ ruT Þ þ FðxÞ:

ð25Þ

The calculation of partial derivatives is required to evaluate the interface curvature and the normal vector. To minimize the discretization error, these derivatives are calculated using compact finite-difference stencils as follows:

@wðxÞ 1 X ¼ 2 wi wðx þ ei Þeia : @xa cs i

ð26Þ

In addition, the original recoloring step is modified by an anti-diffusion scheme proposed by Latva-Kokko and Rothman [50], which can solve the lattice pinning problem and create a symmetric distribution of particles around the interface to further reduce the spurious velocities. So the post-collision, post-segregation (recolored) particle distribution functions for the red and blue fluids can be calculated by

qR qq f þ þ b R B wi cos ujei j; qR þ qB i qR þ qB qB qq fiB ¼ f þ  b R B wi cos ujei j; qR þ qB i qR þ qB

fiR ¼

ð27Þ ð28Þ

where fiþ denotes the post-collision, pre-segregation value of the total particle distribution function along the ith lattice direction; b is the segregation parameter, fixed at 0.7, to maintain a narrow interface thickness and keep low spurious velocities [63] (our recent numerical study also indicates that b = 0.7 is necessary to reproduce correct droplet behavior [47]); u is the angle between the color gradient rqN and the lattice vector ei, which is defined as

cos u ¼

ei  rqN : jei krqN j

ð29Þ

3. Thermal lattice Boltzmann model using the passive-scalar approach Although LB method has been successful in isothermal flow simulations, significant effort is still required for thermal flows. The current thermal lattice Boltzmann (TLB) models can be classified into two major categories: the multi-speed approach and the passive-scalar approach. The multi-speed approach can be regarded as an extension of the isothermal LB models. It uses a large set of discrete velocities with higher order velocity terms in the equilibrium distribution function. Therefore, the macroscopic energy conservation equation can be obtained correctly. However, the multi-speed approach suffers from severe numerical instability and an unphysical fixed Prandtl number due to the single relaxation time. However, if

4438

H. Liu et al. / Journal of Computational Physics 231 (2012) 4433–4453

the viscous dissipation and compression work done by the pressure are negligible, the temperature evolution can be determined by a much simpler passive-scalar equation, which can be solved through an additional LB equation [64]. In order to simulate the evolution of temperature field (i.e., the passive-scalar field), another particle distribution function gi is used:

g i ðx þ ei dt ; t þ dt Þ ¼ g i ðx; tÞ 

g i ðx; tÞ  g eq i ðx; tÞ

sg

þ si ðx; tÞ;

ð30Þ



where sg is the dimensionless relaxation time related to the thermal diffusivity as k ¼ sg  12 c2s dt ; g eq i ðx; tÞ is the equilibrium distribution function in the ith direction at the location x and the time t, which is chosen to recover the macroscopic convection–diffusion equation:

" g eq i ¼ wi T 1 þ

# ei  u ðei  uÞ2 u  u ; þ  c2s 2c2s 2c4s

ð31Þ

where T is the temperature calculated by

T¼

X

gi :

ð32Þ

i

The heating source term ST is incorporated into the macroscopic convection–diffusion equation via a specific forcing term si of the form

si ¼ wi ST :

ð33Þ

With the macroscopic velocity given by Eq. (24), Eq. (30) can recover the following convection–diffusion equation (i.e. passive-scalar equation):

@ t T þ u  rT ¼ r  ðkrTÞ þ ST :

ð34Þ

Following the definition of fluid viscosity given by Eq. (14), we can define the thermal diffusivity of mixture as

k¼

qR qB k þ k ; qR þ qB R qR þ qB B

ð35Þ

where kR (kB) is the thermal diffusivity of red (blue) fluid. There are several important dimensionless parameters controlling the thermocapillary motion, i.e., Reynolds number (Re), Marangoni number (Ma), capillary number (Ca), and the ratios of viscosity and thermal conductivity of fluids. Without losing generality, the blue fluid is chosen as the continuous phase in this study, so the dimensionless parameters are defined as follows:

Re ¼

LU

mB

;

Ma ¼

LU ¼ Re  Pr; kB

Ca ¼

qU mB ; rref

ð36Þ

and

k¼

mR kR ; v¼ ; mB kB

ð37Þ

where U and L are the characteristic velocity and length of the system, Pr is the Prandtl number, and k and v are the viscosity ratio and the thermal diffusivity ratio. 4. Boundary conditions For LB simulation of fluid flow, no-slip boundary condition is usually imposed at the solid walls using the bounce-back scheme. However, the bounce-back scheme only gives first-order numerical accuracy. To improve it, many boundary conditions have been proposed [65–69]. In this work, we employ a boundary condition based on the original Zou and He boundary condition [68], which not only gives second-order numerical accuracy, but also prevents mass leakage across the solid boundary. For the sake of simplicity, we explain the boundary conditions in the context of a D2Q9 model. For example, if the solid wall is assumed to be perpendicular to the y-direction with the lattice velocity e4, e7 and e8 pointing into the computation domain, and f4, f7 and f8 are unknown after propagation. These unknown variables can be determined by satisfying no-slip boundary condition on walls and bounce-back of the non-equilibrium distribution rule, i.e.,

qu ¼ f1 þ f5 þ f8  ðf3 þ f6 þ f7 Þ ¼ 0; qv ¼ f2 þ f5 þ f6  ðf4 þ f7 þ f8 Þ ¼ 0;

ð39Þ

f4  f4eq ¼ f2  f2eq :

ð40Þ

ð38Þ

H. Liu et al. / Journal of Computational Physics 231 (2012) 4433–4453

4439

From the above equations, we can obtain the unknown distribution functions as

1 f 7 ¼ f5 þ ðf1  f3 Þ; 2

f4 ¼ f2 ;

1 f 8 ¼ f6  ðf1  f3 Þ: 2

ð41Þ

To prevent the mass leakage across the solid wall, f0 should be modified as

f0 ¼ f0 þ f2 þ f5 þ f6  ðf4 þ f7 þ f8 Þ:

ð42Þ

where the superscript ‘⁄’ denotes the post-segregation, and pre-propagation values. At the inlet, a parabolic velocity profile is specified directly following Zou and He [68]. At the outlet, we impose a stressfree boundary condition, which can conserve the total mass of the flow system and improve the numerical stability for the flows with low Reynolds numbers. To illustrate it, we assume that the inlet boundary is set at the left-hand end of the system, x = xmin; the outlet boundary is set at the right-hand end, x = xmax; and the inlet velocity of fluid is directed to the right. With the boundary condition of Zou and He [68], the net mass at the inlet x = (xmin, y) follows:

  min ðyÞ ¼ f1 ðxÞ þ f5 ðxÞ þ f8 ðxÞ  f3 ðxÞ þ f6 ðxÞ þ f7 ðxÞ :

ð43Þ

The stress-free boundary condition at the outlet x = (xmax, y) can be implemented as

fi ðxÞ ¼



fi ðx0 Þ

i ¼ 3; 6; 7

fi ðxÞ þ df ðxÞ i ¼ 0:

ð44Þ

where x0 = (xmax  dx, y), and df is introduced to conserve the total mass of the system, which is defined as

df ðxÞ ¼ f1 ðxÞ þ f5 ðxÞ þ f8 ðxÞ  ½f3 ðxÞ þ f6 ðxÞ þ f7 ðxÞ  min ðyÞ:

ð45Þ

To simulate temperature field evolution, two different thermal boundary conditions are applied: the specified temperature at wall surfaces and the adiabatic boundary condition, which can be found in Ref. [70]. Note that it is straightforward to extend these boundary conditions for 3D simulations. 5. Results and discussion 5.1. Model validation To verify our LB model, we first investigate the thermocapillary driven flow in a heated microchannel with two superimposed planar fluids [71]. The geometric setup is shown in Fig. 1. The heights of the upper red fluid and the lower blue fluid are a and b, respectively, while the fluids are of infinite extension in the horizontal direction. We impose a uniform temperature to the upper wall and a sinusoidal temperature (which is higher than that of the upper wall) to the lower wall as

T B ðx; bÞ ¼ T h þ T 0 cosðxxÞ;

ð46Þ

T R ðx; aÞ ¼ T c ;

ð47Þ

and

respectively, where Th > Tc > T0 > 0, and x ¼ 2lp is a wave number with l being a length scale. The above temperature boundary conditions establish a periodic temperature field in the horizontal direction with a period length scale of l. Therefore, it is sufficient to only consider the solution in one period, i.e., l/2 < x < l/2.

Fig. 1. The schematic diagram of two immiscible fluids in a microchannel. The temperatures of the lower and upper walls are TB(x, b) = Th + T0 cos (xx) and TR(x, a) = Tc, respectively, where Th > Tc > T0 and x ¼ 2lp is a wave number.

4440

H. Liu et al. / Journal of Computational Physics 231 (2012) 4433–4453

With the characteristic length and velocity defined by L = b and U ¼ jrTljT 0

b

qmB

respectively, Re, Ma and Ca can be easily

determined by Eq. (36). Assuming Re  1, Ma  1, and Ca  1, it is possible to ignore the convective transport of momentum and energy, and the interface is to remain flat. By solving the simplified governing equations with the stress boundary condition equation (16) at the interface, analytical solutions can be obtained as [71]

T R ðx; yÞ ¼

ðT c  T h Þy þ vT c b þ T h a ~ vÞ sinhða ~; b; ~  xyÞ cosðxxÞ; þ T 0 f ða a þ vb

uRx ðx; yÞ ¼ U max

nh

i

o C R1 þ xðC R2 þ C R3 yÞ coshðxyÞ þ C R3 þ xC R1 y sinhðxyÞ sinðxxÞ;

h

i uRy ðx; yÞ ¼ vU max C R1 y coshðxyÞ þ C R2 þ C R3 y sinhðxyÞ cosðxxÞ;

ð48Þ

ð49Þ ð50Þ

and

T B ðx; yÞ ¼

vðT c  T h Þy þ vT c b þ T h a ~ vÞ½sinhða ~; b; ~Þ coshðxy ~Þ cosðxxÞ; ~Þ  v sinhðxyÞ coshða þ T 0 f ða a þ vb

uBx ðx; yÞ ¼ U max

nh



i

o C B1 þ x C B2 þ C B3 y coshðxyÞ þ C B3 þ xC B1 y sinhðxyÞ sinðxxÞ;

h

i uBy ðx; yÞ ¼ vU max C B1 y coshðxyÞ þ C B2 þ C B3 y sinhðxyÞ cosðxxÞ;

ð51Þ

ð52Þ ð53Þ

in the upper red fluid and the lower blue fluid, respectively. In the above equations, the unknown constants are determined by

~ ¼ bx; ~ ¼ ax; b a ~ coshða ~ 1 ; ~ vÞ ¼ ½v sinhðbÞ ~Þ þ sinhða ~Þ coshðbÞ ~; b; f ða 2

C R1 ¼ C B1 ¼

~Þ sinh ða 2

~Þ  a ~2 sinh ða 2 ~ sinh ðbÞ 2 ~ ~2 b sinh ðbÞ

;

C R2 ¼

;

C B2 ¼

~ aa 2

~Þ  a ~2 sinh ða ~ bb 2 ~ ~2 b sinh ðbÞ

;

C R3 ¼

;

C B3 ¼

ð54Þ ð55Þ ~  sinhð2a ~Þ 2a 2

~Þ  a ~2  2½sinh ða ~ ~ sinhð2bÞ  2b 2 ~ ~2  b 2½sinh ðbÞ

; ;

ð56Þ

and

  T 0 rT ~ vÞhða ~ kÞ; ~; b; ~; b; gða U max ¼ 

ð57Þ

~ vÞ ¼ sinh f ða ~ vÞ; ~; b; ~; b; gða

ð58Þ

qmB

where

and

~ kÞ ¼ ~; b; hða

2 ~ b ~2  ~Þ  a ~2 ½sinh2 ðbÞ ½sinh ða : 2 ~ ~ ~  2b ~ 2 ~ Þ  2a ~ þ ½sinh ða ~Þ  a ~2 ½sinhð2bÞ k½sinh ðbÞ  b ½sinhð2a 2

ð59Þ

The simulations are carried out in a 2D lattice domain with the channel length l = 160 lattice unit (lu), and the initial heights of fluid layer a = b = 40 lu. The periodic boundary conditions are applied on the left and right sides of the domain. On the upper and lower walls, the no-slip boundary conditions are imposed, and the wall temperatures are specified through Eqs. (46) and (47) with Th = 20, Tc = 10 and T0 = 4. The fluid properties are chosen as rT = 5 104, rref = 2.5 102 (at Tref = Tc), mR = mB = 0.2, and kB = 0.2. Consequently, Re, Ma and Ca that are typically Oð0:01Þ, up to Oð0:1Þ. To show the influence of the thermal diffusivity ratio on the velocity and temperature fields, both v = 1 and v = 1/5 are simulated. Fig. 2 shows the contours of the temperature field for the thermal diffusivity ratios: (a) v = 1 and (b) v = 1/5. It can be clearly seen that our simulation results (the solid-line contours) are in good agreement with the analytical solutions (the dashed-line contours), given by Eqs. (48) and (51). We notice that our simulation results slightly deviate from the analytical solutions at v = 1/5, which is due to the finite interface thickness of our LB model and the change of thermal diffusivity across the interface. At small thermal diffusivity ratio, i.e., v = 1/5, the temperature contours become denser in the upper fluid. Also, the slopes of the contour lines from the lower fluid tend to be normal to the interface, reflecting the fact that heat transfer between the lower fluid and interface in the y-direction is close to zero. Fig. 3 shows the comparison of velocity vectors between our simulation results and the analytical solutions of Eqs. (49), (50), (52), and (53) for v = 1 and v = 1/5. The simulated velocity vectors agree well with the analytical solutions except those close to the interface at v = 1/5, which is consistent

H. Liu et al. / Journal of Computational Physics 231 (2012) 4433–4453

4441

Fig. 2. The temperature contours for the fluid systems with thermal diffusivity ratios of (a) v = 1 and (b) v = 1/5. The simulation results and the analytical solutions are represented by the solid and dashed lines, labeled with temperature values.

with the deviation in the corresponding temperature field. Finally, it can be found by comparing Fig. 3(a) and (b) that the thermocapillary-driven convection will strengthen with the decrease of the thermal diffusivity ratio. 5.2. Thermocapillary migration of deformable droplets Thermocapillary migration of droplets was first analyzed by Young et al. [30] for the flows with infinitesimal Reynolds and Marangoni numbers, in which the convective transport of momentum and energy can be neglected compared to molecular transport of these quantities. The governing equations can therefore be linearized. The terminal migration velocity (also known as YGB velocity) of a spherical droplet (red fluid) with radius R in a constant temperature gradient rT1, within an unbounded blue fluid medium can be obtained as

U YGB ¼

2U ; ð2 þ vÞð2 þ 3kÞ

ð60Þ

where the nominal thermocapillary velocity U = rTrT1R/(qmB), which is also used to define Reynolds and Marangoni numbers, i.e., Re = UR/mB and Ma = UR/kB. A droplet of radius R = 16dx is placed inside a 3D box of size of 7.5R 7.5R 15R. The center of droplet is initially located at the center of the 3D computational box. No-slip boundary conditions are imposed on the top and bottom walls, and periodic boundary conditions are used in the x and y directions (the side faces). A linear temperature field is imposed in

4442

H. Liu et al. / Journal of Computational Physics 231 (2012) 4433–4453

Fig. 3. The velocity vectors for the fluid systems with thermal diffusivity ratios of (a) v = 1 and (b) v = 1/5. The velocity vectors are shown at every third grid point. The simulation results and the analytical solutions are represented by the red and blue lines with arrow, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

the z-direction, with T = 0 on the bottom wall and T = 24 on the top wall, resulting in rT1 = 0.1. In order to estimate the accuracy of the 3D LB color-fluid model, we first carry out the numerical simulation with the fluid properties of mR = mB = 0.2, kR = kB = 0.2, Tref = 12, rref = 1.875 103, and rT = 1.5625 104. Using these values, the theoretical migration velocity  104 , and the Reynolds and Marangoni numbers are 0.1. In the simulations, the droplet of a spherical droplet UYGB is 1:666 velocity ud is calculated as

R

N >0Þ

Vðq ud ðtÞ ¼ R

uz dV

dV VðqN >0Þ

P ¼

N ðx; tÞÞ x uz ðx; tÞNð P ; N ðx; tÞÞ Nð x

q

q

ð61Þ

H. Liu et al. / Journal of Computational Physics 231 (2012) 4433–4453

4443

Fig. 4. The time evolution of normalized migration velocity of a 3D spherical droplet at Ma = Re = 0.1. The dashed line represents the theoretical prediction in the limit of vanishing Reynolds and Marangoni numbers, while the solid line is our simulation result.

where the function N(qN) is defined as

NðqN Þ ¼



1; ðqN > 0Þ; 0;

ðqN 6 0Þ:

ð62Þ

Fig. 4 shows the evolution of the numerically calculated droplet velocity normalized by UYGB for the test case of Ma = Re = 0.1. The dimensionless time is defined as t⁄ = Ut/R. Obviously, our numerical result is in good agreement with the theoretical prediction since the effects of convective transport of momentum and energy are negligible in this test case. For thermocapillary migration of deformable droplet at large Marangoni numbers, no theoretical model is available. However, our model is able to directly simulate thermocapillary flows. Fig. 5 presents the evolutions of droplet migration velocity for three different Marangoni numbers, i.e. Ma = 1, 10, and 100, at the fixed Reynolds number of 1. The Marangoni number is achieved by varying kR and kB with v = 1. We also choose r0 = 2.5 102, rT = 1.5625 103, and keep all the other simulation parameters the same as those in the above test case. The early transient of droplet motion exhibits the same characteristics as those reported by Oliver and De Witt [33] and Welch [72], which is caused by the initial conditions used in our numerical simulations, i.e., ujt=0 = 0 and Tjt=0 = zrT1. The thermocapillary migration can reach steady state for all Marangoni numbers. It is evidenced that the terminal migration velocity is a monotonically decreasing function of the Marangoni number, which is consistent with the previous theoretical and numerical studies for the case of nondeformable droplets/bubbles

Fig. 5. The time evolutions of droplet migration velocity at different Marangoni numbers for Re = 1.

4444

H. Liu et al. / Journal of Computational Physics 231 (2012) 4433–4453

Fig. 6. The temperature field around the droplet (the left figures) and the corresponding velocity vectors (the right figures) in a reference frame moving with the droplet in the x  z meridian plane at: (a) Ma = 1, (b) Ma = 10, and (c) Ma = 100.

[31,33,73,36]. The dependence of the terminal migration velocity on the Marangoni number can easily be explained by the isotherms surrounding the droplet, which are shown in Fig. 6 where the temperature value is labeled on each contour. Obviously, the enhanced convective transport of momentum and energy with the increase of the Marangoni number results in the wrapping of the isotherms around the front of the droplet, leading to a substantial reduction of the temperature gradient at the droplet interface. Small average temperature gradient at interface will reduce the driving force for the droplet

H. Liu et al. / Journal of Computational Physics 231 (2012) 4433–4453

4445

Fig. 7. The time evolution of droplet shape, position and streamlines surrounding the droplet for Qs = 0.06 with the times taken as: (a) t⁄ = 0.3, (b) t⁄ = 0.48, (c) t⁄ = 0.63, (d) t⁄ = 0.72, (e) t⁄ = 0.78, (f) t⁄ = 0.96, (g) t⁄ = 1.11, (h) t⁄ = 1.2, (i) t⁄ = 1.32, (j) t⁄ = 1.44, (k) t⁄ = 1.47, (l) t⁄ = 1.5, (m) t⁄ = 1.59, (n) t⁄ = 1.68, and (o) t⁄ = 1.86. The solid lines without arrows are zero contours of qN, and the lines with arrows are the streamlines.

migration. The velocity vectors are also plotted in Fig. 6 for each case in a frame moving with the droplet centroid. Relative to the migrating droplet, the flow pattern within the droplet exhibits recirculation flow that is similar to the Hill’s spherical vortex [74]. It is evidenced that the vortex intensity weakens as the Marangoni number increases. 5.3. Thermocapillary flow for droplet manipulation in a microchannel Historically, the research on thermocapillary flow of moving droplet has mainly been focused in the microgravity conditions where the interfacial effects are dominant. Microfluidics has opened up a new research direction where bulk phenomena are negligible in comparison with interfacial effects due to large surface-to-volume ratio and low Reynolds number. A recent experimental study [29] showed that the localized heating from a laser could block the motion of a water droplet in a microchannel, acting as a microfluidic valve for microdroplet manipulation. The authors further demonstrated that this droplet behavior could be explained as a result of the forces acting on the droplet due to thermocapillary flow. However, their explanation was based on the numerical solution of the depth-averaged Stokes equations for a static circular droplet in an infinite domain with a pre-specified temperature distribution [75]. It is well-known that no blocking will occur when

4446

H. Liu et al. / Journal of Computational Physics 231 (2012) 4433–4453

the laser source is not applied. One may ask whether a threshold value of laser intensity exists, beyond which the droplet is forced to stop. In addition, in order to better manipulate droplet dynamical behavior, it is important to understand the evolution of the velocity and temperature fields around droplet. To tackle these challenges, we apply our model to study thermocapillary flows of water droplet in an oil carrier fluid in a microchannel, in which a laser is used to locally heat the fluid. For simplicity, here we approximate the heating from the laser as a local heat source of a Gaussian form:

( ST ¼

Q se 0;



ðxxs Þ2 þðyys Þ2 w2

; if

2

½ðx  xs Þ2 þ ðy  ys Þ2  6 d ;

ð63Þ

otherwise;

where Qs is the maximum heat flux of the hot spot, xs and ys are the position of hot spot where the heat flux is taken to be of the maximum value, d is the size of the diffused hot spot, and w is a parameter controlling the profile of heat flux. We carry out the simulations in a 360 100 lattice domain with a velocity inlet at the left boundary and a stress-free outflow boundary condition at the right boundary, respectively. At the upper and lower walls, the no-slip boundary conditions are imposed. We specify a fixed temperature T = 0 at all the boundaries except an adiabatic boundary condition, i.e., @T/ @n = 0, at the outlet boundary. A circular droplet of radius R = 25 lattices is initially located at (100, 50). The fluid is locally heated by a laser lighting source with its distribution given by Eq. (63), in which xs = 200, ys = 50, w = 4.0, d = 6.0, and Qs varies from 0.06 to 0.16. To reproduce the thermocapillary blocking, the interfacial tension to temperature gradient rT should be positive [29]. Here we take rT = 5 104, and rref = 0.02 at the reference temperature Tref = 0. The dimensionless numbers for describing the droplet dynamical behavior are chosen as Ca = 2.5 103, Re = Ma = 2.0, and k = v = 1.0, where the characteristic length and velocity are the channel height H and the mean velocity of inlet, Uin. When a laser heating source is applied, we observe that the droplet motion can be blocked to some extent, which depends on the magnitude of the laser heating source. As Qs increases from 0.06 to 0.16, the droplet motion can undergo two states. At a small Qs, the droplet motion is partially blocked. When the Qs increases beyond a threshold value i.e. Qs = 0.085, the droplet motion is completely blocked, so the droplet finally rests at or before the center (xs, ys) of the heating source. Fig. 7 shows the time evolution of the droplet shape, position, and the streamlines around the moving droplet for Qs = 0.06, which corresponds to the state that the droplet motion is partially blocked. The dimensionless time is chosen as t⁄ = ct, where c = Uin/ H is the shear rate. At an early time, i.e., t⁄ = 0.3, the streamlines are almost straight as the droplet is far away from the laser heating source. Note, although we can observe the streamlines passing through the droplet interface in Fig. 7(a), there is no mass exchange between the droplet and the carrier fluid. This can be clearly seen from Fig. 8, in which the streamlines of t⁄ = 0.3 are plotted in the reference frame moving with the droplet centroid. Therefore, it is not surprising to see a mass-conserving droplet at the subsequent times in Fig. 7. As the droplet moves towards the heating source, the streamlines become warped, and a small vortex pair with two counter-rotating vortices appears at the front region inside the droplet (see Fig. 7(b)). At t⁄ = 0.63, the droplet moves further towards the heating source. The size of vortex pair inside the droplet increases, and two small vortices start to grow at the upper and lower walls close to the heating source. As the time evolves, the internal vortex pair and the wall vortices both grow (see Fig. 7(d)), and their sizes reach a maximum value when the front end of droplet just passes the center of the

Fig. 8. The streamlines (black lines) surrounding the droplet in a frame moving with the droplet for Qs = 0.06 at t⁄ = 0.3. The solid circular line is zero contours of qN, and the lines with arrows are the velocity vectors.

H. Liu et al. / Journal of Computational Physics 231 (2012) 4433–4453

4447

heating source (see Fig. 7(e)). At t⁄ = 0.96, the sizes of wall vortices reduce significantly, whereas the internal vortex pair has negligible change. At t⁄ = 1.11, we can notice that the wall vortices has vanished and the internal vortex pair also start to decrease. The internal vortex pair continues to decrease and completely vanishes at t⁄ = 1.32 (see Fig. 7(h) and (i)). Notice that the streamlines once again return to the straight lines at t⁄ = 1.32, which are similar to those in Fig. 7(a). This is because the centroid of droplet roughly coincides with the center of heating source, and a large part of droplet interface is well situated at the isotherm, which can be clearly seen in Fig. 9(i), so that the Marangoni convection induced by the temperature gradient is negligible. As the droplet leaves the center of heating source and moves to the right, the symmetry of temperature field is broken, and the two small wall vortices emerge at the upper and lower sides of droplet (see Fig. 7(j)). Quickly, both of the wall vortices grow up, and a small vortex pair arises at the tail of the droplet (see Fig. 7(k)). However, this tail vortex pair rapidly disappears, and the size of wall vortices reach a peak value at t⁄ = 1.47. After that, the wall vortices begin to decrease until they vanish (see Fig. 7(m) and (n)). As the droplet continues to move away from the heating source, the streamlines gradually become less warped (see Fig. 7(o)). Compared to the velocity field shown in Fig. 7, the temperature field has negligible change, as can be seen in Fig. 9. The main reason is that thermal diffusion is the dominant energy transport mechanism at small Marangoni numbers. From the far field to the Center of Heating Source (CHS), the temperature increases, and the isotherms become denser, leading to a larger temperature gradient near the CHS. However, one should notice that the temperature gradient is positive on the left

Fig. 9. The time evolution of the temperature field surrounding the moving droplet for Qs = 0.06 with the times taken to be the same as those in Fig. 7. At each time step, the zero contour of qN is indicated by the solid circle, and the temperature contours are plotted with every 0.5 spacing and the initial value of 0.5 starting from the far field towards the center of the heating source, which is fixed at (x = 200, y = 50).

4448

H. Liu et al. / Journal of Computational Physics 231 (2012) 4433–4453

Fig. 10. The x-position of droplet centroid as a function of dimensionless time (ct) for Qs = 0, 0.06 and 0.08.

Fig. 11. The time evolution of droplet shape, position and streamlines surrounding the droplet for Qs = 0.086 with the times taken as: (a) t⁄ = 0.3, (b) t⁄ = 0.48, (c) t⁄ = 0.63, (d) t⁄ = 0.72, (e) t⁄ = 2.01, and (f) t⁄ = 4.98. The solid lines without arrows are zero contours of qN, and the lines with arrows are the streamlines.

side of CHS but it is negative on the right side. Considering the Marangoni force that causes the fluid to move towards the regions of high interfacial tension (or high temperature) along the interface, it is expected that the Marangoni force acting on

H. Liu et al. / Journal of Computational Physics 231 (2012) 4433–4453

4449

Fig. 12. The position of droplet centroid as a function of dimensionless time ct for Qs = 0.086, 0.1, 0.12, 0.14, and 0.16.

Fig. 13. The flowfield around the resting droplet in the final steady state for (a) Qs = 0.086 and (b) Qs = 0.16. The droplet is indicated by the circular line without arrows located in the middle of the two plates.

the droplet serves to block its motion when the droplet is located at the left side of the CHS. Once the droplet moves across the CHS, the droplet interface located at the right side of CHS will produce an opposite Marangoni force that accelerates the droplet motion. When the droplet is far from the CHS, the temperature gradient is so small that the droplet moves at an approximately constant velocity. Therefore, one can observe that the droplet motion successively undergoes four stages:

4450

H. Liu et al. / Journal of Computational Physics 231 (2012) 4433–4453

constant velocity, deceleration, acceleration, and constant velocity. This can be clearly seen in Fig. 10, which shows the xposition of droplet centroid as a function of the dimensionless time (ct), calculated by

R

Vðq xd ðtÞ ¼ R

N >0Þ

xdV

VðqN >0Þ

dV

P ¼

N ðx; tÞÞ x xðx; tÞNð P : N ðx; tÞÞ Nð x

q

q

ð64Þ

While the droplet motion experiences the transition of constant velocity, deceleration, acceleration, and constant velocity for all nonzero Qs below 0.085, the droplet blocking becomes more significant as Qs increases, e.g., see Qs = 0.08 in Fig. 10. When the (partial) blocking becomes significant, an important feature is that the droplet motion lasts a longer time with a lower velocity, whereas it cannot be completely blocked. By differentiating xd with respect to time, we can obtain that, the droplet initially moves at a constant velocity ud = 1.38Uin, which is independent of Qs. At the late stage, we can observe that the slopes of lines for Qs = 0.06 and 0.08 are a little steeper than that for Qs = 0, which means that the droplet moves at a velocity slightly larger than 1.38Uin. This is because a temperature gradient remains at the grid points close to the outlet

Fig. 14. Time evolution of the highest temperature in the whole microchannel for various Qs.

Fig. 15. The highest temperature at the steady state as a function of Qs. The discrete points are our simulation results, and the solid line is the linear fitting.

H. Liu et al. / Journal of Computational Physics 231 (2012) 4433–4453

4451

boundary although its value is small, and the Marangoni force is quasi-balanced by the viscous force due to the relative motion between the droplet and the surrounding fluid. It should be noted, strictly speaking, the late-stage of constant velocity is a process of extremely slow deceleration. Fig. 11 shows the time evolution of droplet shape, position, and streamline patterns around the moving droplet for Qs = 0.086, where the droplet motion is completely blocked due to the thermo-induced Marangoni stresses. At t⁄ = 0.3, the streamlines are almost horizontal and straight, whereas small deformation occurs, indicating that the heat has diffused to this area. As the droplet moves forward, a vortex pair with two counter-rotating vortices appears at the front region inside the droplet at t⁄ = 0.48, but their sizes are bigger than those in Fig. 7(b) due to a stronger laser heating source. The internal vortex pair grows continuously, and two small vortices emerge at the upper and lower walls with the x-coordinate of vortex core close to xs ((see Fig. 11(c))). The internal vortex pair and two wall vortices grow (see Fig. 11(d)), until reaching the maximum value at t⁄ = 2.01. Since then, the positions of the droplet and vortices, as well as their sizes and shapes remain the same (Fig. 11(f)), indicating that the droplet has been motionless and reached the steady state. It can be clearly seen that the internal vortex pair completely fills the droplet and the droplet behaves as a solid sphere. The formation of these stable vortex/recirculation regions is a result of force balance between the viscous stresses and the thermocapillary force. However, the positions of two vortices outside the droplet are different from the experimental and numerical observations of Baround et al. [29,75], in which the four vortices including the vortex pair inside the droplet closely distribute around the laser spot. The difference attributes to the fact that there is no carrier fluid flowing through the droplet in their simulation and experiment (if the thin film flow between the droplet and the solid wall in the experiment is negligible), whereas in the present simulations, considering the continuous injection of the incompressible carrier fluid at the inlet, the carrier fluid needs to flow through the motionless droplet and then stream out of the outlet boundary. Therefore, the two vortices outside the droplet are ‘‘driven away from’’ the laser spot. We also note that, the confinement of solid wall was not considered in the simulation of Baround et al. [75], and the droplet was assumed to be spherical. When we continue to increase Qs, the droplet motion will be completely blocked before the droplet arrives at the laser spot. The final rest position xd(t ? 1) depends on the magnitude of Qs. Fig. 12 plots the position of droplet centroid as a function of dimensionless time ct for Qs = 0.086, 0.1, 0.12, 0.14, and 0.16. For all the Qs, the droplet first moves at a constant velocity, and then decelerates until it comes to rest. The final rest position xd decreases as Qs increases. Fig. 13 shows the comparison of flowfield when the droplet is finally at rest for Qs = 0.086 and Qs = 0.16. Typically, a pair of vortices completely fill the droplet for both Qs, whereas the external vortices become smaller as Qs increases. The different droplet rest positions are attributed to the variation of temperature field. This can be seen in Fig. 14, which gives the time evolution of the maximum temperature, Tmax, for various Qs in the whole computational domain. We note that the maximum temperature occurs at the center of laser heating source. At the initial stage, the maximum temperature increases dramatically, and then decreases slightly until it reaches a constant value. In Fig. 15, we also plot the final steady Tmax as a function of Qs. It can be seen that the final steady Tmax is proportional to Qs. This characteristics of heat transfer can be easily understood because the heat diffusion is dominant in a typical microchannel and the governing equation for temperature will become

r  ½krT þ Q s e

ðxxs Þ2 þðyys Þ2 w2

¼ 0;

ð65Þ

when the system reaches the steady state. From Eq. (65), it can be readily obtained that the temperature is proportional to Qs. Obviously, the proportional relationship between the final steady Tmax and Qs is consistent with this analysis.

6. Conclusion A lattice Boltzmann multiphase model has been presented to simulate thermocapillary flows. The interfacial forces, including the interfacial tension forces and the Marangoni stresses due to the gradient of interfacial tension, are modeled using the concept of continuum surface force; and a recoloring algorithm proposed by Latva-Kokko and Rothman is applied for separating two fluids, which can overcome ‘‘lattice pinning’’, a problem that may prevent the interface from moving. An additional convection–diffusion equation is also solved in the LB framework to obtain temperature field, which is related to the interfacial tension by the equation of state. A stress-free boundary condition is introduced to treat the outflow boundary, which can conserve the total mass of an incompressible system, thus improving the numerical stability for creeping flows. Excellent agreement with theoretical predictions in the limit of zero Marangoni number are obtained. It has been found that the migration velocity of droplet diminishes when the Marangoni number increases. We numerically demonstrate the laser heating can generate the thermocapillary forces to block the droplet motion. When the intensity of laser heating source Qs increases from 0.06 to 0.16, the droplet motion will undergo transition from partial blocking to complete blocking with the transition occurring at Qs = 0.085, and the blocking becomes increasingly pronounced. In the case of partial blocking, the droplet motion successively experiences four stages: constant velocity, deceleration, acceleration, and constant velocity. In the case of complete blocking, the droplet firstly moves at a constant velocity, and then decelerates until it is stopped. In addition, as the droplet moves towards the laser heating source, both of the internal vortex pair and the external wall vortices grow until a maximum size has been reached. When the droplet finally becomes motionless, the external wall vortices decrease as Qs increases, but the internal vortex pair always fills the droplet. We also observe the highest temperature in the steady state is proportional to Qs, which is consistent with the theoretical analysis.

4452

H. Liu et al. / Journal of Computational Physics 231 (2012) 4433–4453

Acknowledgments This work was supported in part by the Center for Energy and Sustainability and by the NASA University Research Center at the California State University, Los Angeles. References [1] H. Song, D.L. Chen, R.F. Ismagilov, Reactions in droplets in microfluidic channels, Angew. Chem. Int. Ed. 45 (44) (2006) 7336–7356. [2] S. Mohr, Y. Zhang, A. Macaskill, P.J.R. Day, R.W. Barber, N.J. Goddard, D.R. Emerson, P.R. Fielden, Numerical and experimental study of a droplet-based PCR chip, Microfluid. Nanofluid. 3 (2007) 611–621. [3] Y. Zhang, P. Ozdemir, Microfluidic DNA amplification-a review, Analytica Chimica Acta 638 (2009) 115–125. [4] H. Song, R.F. Ismagilov, Millisecond kinetics on a microfluidic chip using nanoliters of reagents, J. Am. Chem. Soc. 125 (47) (2003) 14613–14619. [5] L.S. Roach, H. Song, R.F. Ismagilov, Controlling nonspecific protein adsorption in a plug-based microfluidic system by controlling interfacial chemistry using fluorous-phase surfactants, Anal. Chem. 77 (2005) 785–796. [6] A. Huebner, M. Srisa-Art, D. Holt, C. Abell, F. Hollfelder, A.J. deMello, J.B. Edel, Quantitative detection of protein expression in single cells using droplet microfluidics, Chem. Commun. (2007) 1218–1220. [7] B.T.C. Lau, C.A. Baitz, X.P. Dong, C.L. Hansen, A complete microfluidic screening platform for rational protein crystallization, J. Am. Chem. Soc. 129 (3) (2007) 454–455. [8] L. Li, R.F. Ismagilov, Protein crystallization using microfluidic technologies based on valves, droplets, and slipchip, Ann. Rev. Biophys. 39 (2010) 139– 158. [9] L.-H. Hung, K.M. Choi, W.-Y. Tseng, Y.-C. Tan, K.J. Shea, A.P. Lee, Alternating droplet generation and controlled dynamic droplet fusion in microfluidic device for CdS nanoparticle synthesis, Lab Chip 6 (2006) 174–178. [10] T. Hatakeyama, D.L. Chen, R.F. Ismagilov, Microgram-scale testing of reaction conditions in solution using nanoliter plugs in microfluidics with detection by MALDI-MS, J. Am. Chem. Soc. 128 (8) (2006) 2518–2519. [11] L.F. Cheow, L. Yobas, D.-L. Kwong, Digital microfluidics: droplet based logic gates, Appl. Phys. Lett. 90 (2007) 054107. [12] M. Prakash, N. Gershenfeld, Microfluidic bubble logic, Science 315 (5813) (2007) 832–835. [13] C.-G. Yang, Z.-R. Xu, J.-H. Wang, Manipulation of droplets in microfluidic systems, Trend. Anal. Chem. 29 (2) (2010) 141–157. [14] T. Thorsen, R.W. Roberts, F.H. Arnold, S.R. Quake, Dynamic pattern formation in a vesicle-generating microfluidic device, Phys. Rev. Lett. 86 (18) (2001) 4163–4166. [15] P. Garstecki, M.J. Fuerstman, H.A. Stone, G.M. Whitesides, Formation of droplets and bubbles in a microfluidic T-junction-scaling and mechanism of break-up, Lab Chip 6 (2006) 437–446. [16] S.L. Anna, H.C. Mayer, Microscale tipstreaming in a microfluidic flow focusing device, Phys. Fluids 18 (2006) 121512. [17] M.G. Pollack, R.B. Fair, A.D. Shenderov, Electrowetting-based actuation of liquid droplets for microfluidic applications, Appl. Phys. Lett. 77 (2000) 1725. [18] J.A. Schwartz, J.V. Vykoukal, P.R.C. Gascoyne, Droplet-based chemistry on a programmable micro-chip, Lab Chip 4 (2004) 11–17. [19] D. Chatterjee, H. Shepherd, R.L. Garrell, Electromechanical model for actuating liquids in a two-plate droplet microfluidic device, Lab Chip 9 (2009) 1219–1229. [20] A.A. Darhuber, J.P. Valentino, J.M. Davis, S.M. Troian, S. Wagner, Microfluidic actuation by modulation of surface stresses, Appl. Phys. Lett. 82 (2003) 657. [21] J.Z. Chen, S.M. Troian, A.A. Darhuber, S. Wagner, Effect of contact angle hysteresis on thermocapillary droplet actuation, J. Appl. Phys. 97 (2005) 014906. [22] T. Franke, A.R. Abate, D.A. Weitz, A. Wixforth, Surface acoustic wave (SAW) directed droplet flow in microfluidics for PDMS devices, Lab Chip 9 (2009) 2625–2627. [23] M. Okochi, H. Tsuchiya, F. Kumazawa, M. Shikida, H. Honda, Droplet-based gene expression analysis using a device with magnetic force-based-droplethandling system, J. Biosci. Bioeng. 109 (2) (2010) 193–197. [24] Y. Zhang, S. Park, K. Liu, J. Tsuan, S. Yang, T.-H. Wang, A surface topography assisted droplet manipulation platform for biomarker detection and pathogen identification, Lab Chip 11 (2011) 398–406. [25] R. Di Leonardo, G. Ruocco, J. Leach, M.J. Padgett, A.J. Wright, J.M. Girkin, D.R. Burnham, D. McGloin, Parametric resonance of optically trapped aerosols, Phys. Rev. Lett. 99 (1) (2007) 010601. [26] D. McGloin, D.R. Burnham, M.D. Summers, D. Rudd, N. Dewar, S. Anand, Optical manipulation of airborne particles: techniques and applications, Faraday Discuss. 137 (2008) 335–350. [27] J.-P. Delville, M.R. de Saint Vincent, R.D. Schroll, H. Chraïbi, B. Issenmann, R. Wunenburger, D. Lasseux, W.W. Zhang, E. Brasselet, Laser microfluidics: fluid actuation by light, J. Opt. A-Pure Appl. Opt. 11 (3) (2009) 034015. [28] L.E. Sciven, C.V. Sternling, The Marangoni effects, Nature 187 (1960) 186–188. [29] C.N. Baroud, J.-P. Delville, F. Gallaire, R. Wunenburger, Thermocapillary valve for droplet production and sorting, Phys. Rev. E 75 (4) (2007) 046302. [30] N.O. Young, J.S. Goldstein, M.J. Block, The motion of bubbles in a vertical temperature gradient, J. Fluid Mech. 6 (03) (1959) 350–356. [31] N. Shankar, R. Shankar Subramanian, The Stokes motion of a gas bubble due to interfacial tension gradients at low to moderate Marangoni numbers, J. Colloid Interface Sci. 123 (2) (1988) 512–522. [32] J. Chen, Y. Lee, Effect of surface deformation on thermocapillary bubble migration, AIAA J. 30 (1992) 993–998. [33] D. Oliver, K.D. Witt, Transient motion of a gas bubble in a thermal gradient in low gravity, J. Colloid Interface Sci. 164 (2) (1994) 263–268. [34] R. Subramanian, R. Balasubramaniam, The Motion of Bubbles and Drops in Reduced Gravity, Cambridge University Press, London, 2001. [35] R. Subramanian, R. Balasubramaniam, G. Wozniak, Physics of Fluids in Microgravity, Taylor & Prancis, London, 2002. Ch. Fluid mechanics of bubbles and drops, pp. 149–177. [36] Z. Yin, P. Gao, W. Hu, L. Chang, Thermocapillary migration of nondeformable drops, Phys. Fluids 20 (2008) 082101. [37] S. Rybalko, N. Magome, K. Yoshikawa, Forward and backward laser-guided motion of an oil droplet, Phys. Rev. E 70 (2004) 046301. [38] K.T. Kotz, K.A. Noble, G.W. Faris, Optical microfluidics, Appl. Phys. Lett. 85 (2004) 2658. [39] C.N. Baroud, M.R. de Saint Vincent, J.-P. Delville, An optical toolbox for total control of droplet microfluidics, Lab Chip 7 (2007) 1029–1033. [40] M.R. de Saint Vincent, R. Wunenburger, J.-P. Delville, Laser switching and sorting for high speed digital microfluidics, Appl. Phys. Lett. 92 (2008) 154105. [41] W. Shyy, R.W. Smith, H.S. Udaykumar, M.M. Rao, Computational Fluid Dynamics with Moving Boundaries, Taylor & Francis, London, 1996. [42] S. Succi, The Lattice Boltzmann Equation For Fluid Dynamics and Beyond, Oxford University Press, Oxford, 2001. [43] X. He, L.-S. Luo, A priori derivation of the lattice Boltzmann equation, Phys. Rev. E 55 (1997) R6333–R6336. [44] M.M. Dupin, I. Halliday, C.M. Care, Simulation of a microfluidic flow-focusing device, Phys. Rev. E 73 (2006) 055701. [45] H. Kusumaatmaja, J. Léopoldès, A. Dupuis, J.M. Yeomans, Drop dynamics on chemically patterned surfaces, Europhys. Lett. 73 (5) (2006) 740. [46] H. Liu, Y. Zhang, Phase-field modeling droplet dynamics with soluble surfactants, J. Comput. Phys. 229 (24) (2010) 9166–9187. [47] H. Liu, Y. Zhang, Droplet formation in microfluidic cross-junctions, Phys. Fluids 23 (2011) 082101. [48] W. Wang, Z. Liu, Y. Jin, Y. Cheng, LBM simulation of droplet formation in micro-channels, Chem. Eng. J. 173 (3) (2011) 828–836. [49] J.U. Brackbill, D.B. Kothe, C. Zemach, A continuum method for previous modeling surface tension, J. Comput. Phys. 100 (1992) 335–354. [50] M. Latva-Kokko, D.H. Rothman, Diffusion properties of gradient-based lattice Boltzmann models of immiscible fluids, Phys. Rev. E 71 (5) (2005) 056702.

H. Liu et al. / Journal of Computational Physics 231 (2012) 4433–4453 [51] [52] [53] [54] [55] [56] [57] [58] [59] [60] [61] [62] [63] [64] [65] [66] [67] [68] [69] [70] [71] [72] [73] [74] [75]

4453

S. Chen, G.D. Doolen, Lattice Boltzmann method for fluid flows, Ann. Rev. Fluid Mech. 30 (1) (1998) 329–364. A.K. Gunstensen, D.H. Rothman, S. Zaleski, G. Zanetti, Lattice Boltzmann model of immiscible fluids, Phys. Rev. A 43 (8) (1991) 4320–4327. X. Shan, H. Chen, Lattice Boltzmann model for simulating flows with multiple phases and components, Phys. Rev. E 47 (3) (1993) 1815. M.R. Swift, W.R. Osborn, J.M. Yeomans, Lattice Boltzmann simulation of nonideal fluids, Phys. Rev. Lett. 75 (5) (1995) 830–833. M.R. Swift, E. Orlandini, W.R. Osborn, J.M. Yeomans, Lattice Boltzmann simulations of liquid-gas and binary fluid systems, Phys. Rev. E 54 (5) (1996) 5041–5052. P. Yuan, L. Schaefer, Equations of state in a lattice Boltzmann model, Phys. Fluids 18 (2006) 042101. M. Sbragaglia, R. Benzi, L. Biferale, S. Succi, K. Sugiyama, F. Toschi, Generalized lattice Boltzmann method with multirange pseudopotential, Phys. Rev. E 75 (2) (2007) 026702. A. Kupershtokh, D. Medvedev, D. Karpov, On equations of state in a lattice Boltzmann method, Comput. Math. Appl. 58 (5) (2009) 965–974. R. van der Sman, S. van der Graaf, Emulsion droplet deformation and breakup with lattice Boltzmann model, Comput. Phys. Commun. 178 (2008) 492– 504. S.V. Lishchuk, C.M. Care, I. Halliday, Lattice Boltzmann algorithm for surface tension with greatly reduced microcurrents, Phys. Rev. E 67 (2003) 036701. I. Halliday, R. Law, C.M. Care, A. Hollis, Improved simulation of drop dynamics in a shear flow at low Reynolds and capillary number, Phys. Rev. E 73 (5) (2006) 056708. Z. Guo, C. Zheng, B. Shi, Discrete lattice effects on the forcing term in the lattice Boltzmann method, Phys. Rev. E 65 (2002) 046308. I. Halliday, A.P. Hollis, C.M. Care, Lattice Boltzmann algorithm for continuum multicomponent flow, Phys. Rev. E 76 (2007) 026708. Y. Peng, C. Shu, Y.T. Chew, Simplified thermal lattice Boltzmann model for incompressible thermal flows, Phys. Rev. E 68 (2003) 026701. A.J.C. LADD, Numerical simulations of particulate suspensions via a discretized Boltzmann equation. Part 1: theoretical foundation, J. Fluid Mech. 271 (1994) 309–311. T. Inamuro, M. Yoshino, F. Ogino, A non-slip boundary condition for lattice Boltzmann simulations, Phys. Fluids 7 (1995) 2928. S. Chen, D. Martnez, R. Mei, On boundary conditions in lattice Boltzmann methods, Phys. Fluids 8 (1996) 2527. Q. Zou, X. He, On pressure and velocity boundary conditions for the lattice Boltzmann BGK model, Phys. Fluids 9 (1997) 1591. Z.-L. Guo, C.-G. Zheng, B.-C. Shi, Non-equilibrium extrapolation method for velocity and pressure boundary conditions in the lattice Boltzmann method, Chinese Phys. 11 (4) (2002) 366. C.-H. Liu, K.-H. Lin, H.-C. Mai, C.-A. Lin, Thermal boundary conditions for thermal lattice Boltzmann simulations, Comput. Math. Appl. 59 (7) (2010) 2178–2193. B. Pendse, A. Esmaeeli, An analytical solution for thermocapillary-driven convection of superimposed fluids at zero Reynolds and Marangoni numbers, Int. J. Therm. Sci. 49 (7) (2010) 1147–1155. S.W. Welch, Transient thermocapillary migration of deformable bubbles, J. Colloid Interface Sci. 208 (2) (1998) 500–508. Y. Wang, X. Lu, L. Zhuang, Z. Tang, W. Hu, Numerical simulation of drop Marangoni migration under microgravity, Acta Astronaut. 54 (5) (2004) 325– 335. M.J.M. Hill, On a spherical vortex, Proc. R. Soc. Lond. 55 (1894) 219–224. F. Gallaire, C.N. Baroud, J.-P. Delville, Thermocapillary manipulation of microfluidic droplets: theory and applications, Int. J. Heat Technol. 26 (2008) 161–166.