Laplace pressure versus Marangoni convection in photothermal manipulation of micro droplet

We have numerically demonstrated micro droplet migrations by photothermal manipulation, in order to investigate the driving force and accompanying flows. We focus on an oil (oleic acid) droplet in water that provides a positive temperature dependence of interfacial tension. The present direct numerical simulation employs the volume-of-fluid method and the continuum-surface-force method with consideration of temperature dependency. The driving velocity and force magnitude, which we measured quantitatively, agree with experimental and theoretical data. We have found that the dominant driving force exerted on the microdroplet is the interface normal force, i.e., the Laplace pressure, on a curved interface rather than the tangential force. The tangential component directly triggers the Marangoni convection; however, the induced flows inside the droplet are opposite between the cases of “an oil droplet in water” and “a water droplet in oil”. This may be caused by the large difference between the viscosities of water and oil.


Introduction
Non-contact manipulation of micro-scale droplets in liquid, which is relevant to medical and biological engineering applications and chemical processes in microgravity, has increasingly attracted attention in the field of microfluidics. In recent decades, researchers have proposed several techniques of droplet manipulation by utilizing a local variation of interfacial tension, since the effects of interfacial phenomena are dominant in microfluidics relative to the inertia and buoyancy forces. For instance, microelectromechanical system devices were developed to drive liquid-metal droplets in an electrode channel [1][2][3], and successfully demonstrated electrowetting, which can change local interfacial tension with electric potential. In addition, laser-induced optical force may be used as a non-invasive and precise tool for droplet-based controls [4][5][6]. Its force magnitude is on the order of a pico newton. Compared to the optical force, an optically-induced thermal Marangoni convection may provide a larger resultant force, which would be in the nano-newton order in the micrometer scale [7][8][9][10][11]. Therefore, the photothermal Marangoni convection can be a powerful technique of on-demand bubble/droplet handling in a micro-channel liquid.
We briefly review studies on photothermal manipulation using the optically induced local temperature dependence of surface tension, ∂σ/∂T = 0. Baroud et al. [7] demonstrated that a laser-spot local heating on a liquid interface can induce a thermocapillary force sufficient to block the progress of a water droplet in a microchannel. Takeuchi et al. [8] applied the photothermal technique to move a bubble actively with a scanning laser. They estimated the thermocapillary force on the bubble by observing the behavior of the bubble following a laser spot, and determined the transport limit in terms of the supplied optical power, scanning speed, and bubble size. Muto and Motosuke [11] experimentally demonstrated the photothermal manipulation of a microscale picoliter oil droplet by laser-spot heating. They indirectly determined a driving force that acts on the target droplet, on the assumption that the force balances with the viscous drag, which can be estimated from the driving velocity u d given by a theoretical model. Their model was the Young-Goldstein-Block (YGB) theoretical equation [12,13], which was modified by neglecting the buoyancy, as follows: where λ, μ, and D are the thermal conductivity, viscosity, and droplet diameter, respectively. The subscripts ( w ) and ( o ) denote water and oil, respectively. Muto and Motosuke [11] could directly measure the driving velocity; however, it was actually difficult to precisely obtain the velocity because of limited flow geometry and nonideal, non-uniform, temperature gradient, i.e., ∂T /∂x was not constant. Moreover, detailed measurements of the pressure and force fields in such microfluids are practically impossible.
In this paper, we report a three-dimensional numerical study by means of DNS (direct numerical simulation) to investigate quantitatively the photothermal droplet manipulation. The interfacial tension in the liquid-liquid two-phase flow is simulated by the VoF (volume of fluid) method [14] and CSF (continuous surface force) model [15] with consideration of the temperature dependence. This DNS allows us to simulate systematically both/either the Laplace pressure near the curved interface and/or the thermal Marangoni convection. In this context, we discuss the different contributions of the Laplace pressure and Marangoni convection, both of which should be significant phenomena for microscale fluid dynamics. To the authors' knowledge, there are only a few studies on a flow driven by non-uniform Laplace pressure under a temperature gradient, although many researchers numerically investigated the force balance on a droplet subjected to a Marangoni flow, e.g., [9,[16][17][18][19][20]].
2 Physical model and numerical method 2.1 Configuration of an oil droplet in water Figure 1 shows the analysis object, which is an oil droplet in a water pool. The droplet form is initially spherical with a diameter D in a rectangular water reservoir. The dimensions of the computational domain are L x × L y × L z = 4D × 2.4D × 2.4D with uniform Cartesian grids of N x × N y × N z = 70 × 42 × 42. The grid dependency will be discussed in Section 2.3. This domain size was chosen by considering the required grid resolution and to prevent the droplet from touching the side walls.  The boundary condition is non-slip on the reservoir solid wall. The initial position of the droplet center is set at (x, y, z) = (0.7L x , 0.5L y , 0.5L z ). As for the thermal field, an initially uniform temperature gradient is given in the x direction. One of the control parameters for the thermal field is the constant temperature difference between the cold and hot walls. Similarly, the thermal boundary condition on the other reservoir solid walls is the Dirichlet condition, i.e., the temperature of the cold and hot walls are respectively fixed at T cold and T hot , while the temperature T wall (x) of the side walls changes linearly in the x direction In this study, we simply use the temperature difference between the x-wise ends of the droplet, ΔT ≡ D(∂T /∂x) t=0 = D(T hot − T cold )/L x . The fluid properties we used here are shown in Table 1. Those properties are set with reference to the experimental condition [11]; however, we assume constant properties except for the interfacial tension, i.e., σ = σ 0 + ∂σ/∂T (T − T 0 ). Here, the base value of σ 0 = 1.0 mN/m at T 0 = T cold was determined by reference to an actual measurement [11]. It is observed that the interface between oleic acid (oil) and water exhibits a positive ∂σ/∂T . Under such condition of positive dependence, the droplet will be driven toward the cold-wall side, according to the driving principle and experimental demonstration [11].

Numerical procedure
The governing equations of three-dimensional velocity field u = (u, v, w) are the equation of continuity for incompressible fluid the Navier-Stokes equation with an interfacial tension term f σ by the CSF model and the energy equation for passive-scalar temperature field In addition, the advection equation of the VoF function, F , is also calculated: The value of F indicates a volume ratio of two fluids in each computational cell. The density ρ and the viscosity μ at each cell depend on F , while σ depends on T . The last term on the right-hand side of equation (4) represents the force due to the interfacial tension σ. This term may be separated into two components, as follows where the force directions of f σ⊥ and f σ are normal and parallel to the interface, respectively. By calculating equation (6) using the VoF and donor accepter method [14], we track the liquid-liquid interface by SLIC (simple line interface calculation) method [21]. To compute the surface tension, we use the CSF model, which allows us to calculate the interface curvature κ as well as the unit vector n normal to interface: Here, ∇ s is the gradient tangential to the interface. Then, both the normal force f σ⊥ and the tangential force f σ can be calculated as The Laplace pressure should be induced by the normal force f σ⊥ , while the Marangoni convection may be triggered by the tangential force f σ . The numerical methodology for the DNS itself is briefly described as follows. An HSMAC (SOLA) method [22] was used for the coupling of equations (3) and (4). The second-order Adams-Bashforth method was employed for the time advancement. A finite difference method was adopted for the spatial discretization and we chose a second-order central scheme.

Validation of DNS code and grid resolution
The simulations presented in this paper were run using our developed code. As a test case to examine the validity of our code, a cubic drop in zero gravity and its transformation to spherical shape were chosen under the condition analogous to a two-dimensional square-drop simulation by Brackbill et al. [15]. The fluid properties were the same between the two simulations; however, our result was computed threedimensionally. Figure 2 shows the cross-sectional views of the droplet, which was deformed by unbalanced surface-tension forces (those force vectors are also shown) and oscillated about its equilibrium shape. The drop shape at t = 2 s is nearly circular in cross section, and the normal forces emerge uniformly at the drop surface. These normal forces induce a pressure increase corresponding to the Laplace pressure. The result in Figure 2 matches closely with that obtained by Brackbill et al. [15], although there was a deviation from their simulation at t = 0.2 s. The mismatches in shape and force magnitude at t = 0.2 s are reasonable because of the different degrees of the simulated dimension: the square drop in 2D simulation becomes square cylindrical again at t = 0.2 s; the cubic drop in our 3D simulation takes the form of a regular  octahedron, and as a result, its sides undergo the strong normal forces of surface tension.
As for the present main target, we carried out several DNSs with different numbers of grids, in order to determine the appropriate grid resolution. From Figure 3, in which the obtained |u d | (defined later in Sect. 3.2) attains a saturated value when N x ≥ 70, we conclude that the present grids of N x × N y × N z = 70 × 42 × 42 sufficiently provide the appropriate resolutions for the present study. As will be discussed later, the driving velocity u d fluctuates with time, maybe because of a numerical error that cannot be avoided in the standard VoF method with Cartesian grids. Adopting a high-order VoF method such as the PLIC (piecewise linear interface construction) method and/or body-fitted coordinate method would be a promising countermeasure against this issue. However, we decide to employ the present method of VoF/SLIC, which is rather simpler than the PLIC method and provides a performance that is sufficient to achieve the desired results of this study.  condition in reference [11]. The field in the x-z plane across the center of droplet is visualized. Two snapshots are taken immediately before the migration, as shown in Figure 4a, and after reaching an equilibrium state, in (b). At equilibrium, the driving force should be balanced with the viscous drag force, which increases with increasing velocity of the droplet. We successfully demonstrated the steady migration of a droplet from the hot to cold side, as expected. However, the flow fields seem to be disturbed by spurious currents. A spurious current is a numerical artifact, which occurs as a result of tracking a finite-thickness interface and simulating the Laplace pressure as a body force. We also found that even without temperature gradient, the DNS on a spherical droplet with a constant surface tension exhibited similar spurious currents; however, the shape and position of the droplet were completely static (figure not shown). As observed in Figure 4, although the droplet may be slightly deformed by the flows and spurious currents in this simulation, its spherical shape basically remains owing to the low Weber numbers, refer to Table 2. Moreover, in the thermal fields, heat conduction is dominant in the present configuration, although convective transport was actually considered in our DNS. Therefore, we may neglect the influence of spurious currents on the driving velocity of a thermally manipulated droplet.

Driving velocity and driving force
We systematically tested various conditions in terms of the droplet diameter (D = 30-100 μm) and temperature gradient (ΔT ≡ D∂T /∂x = 1.5-4.0 K). Most conditions shown in Table 2 were referred from previous experimental conditions [11]. The non-dimensional parameters in Table 2 are based on D, σ 0 , and u d predicted by equation (1), and the water property values given in Table 1. is defined as time-and volume-averaged speed within the droplet: where V is the volume of the droplet, and the time integration is performed after reaching the equilibrium state. Shown as well in Figure 5 is a predicted curve of the YGB theoretical equation [12]: u d ≈ −0.379ΔT [mm/s] under the present condition given in Table 1. Here, as mentioned above, the influence of spurious currents on the driving velocity may be neglected, because the currents take a point-symmetric form about the center of the droplet and its spatial integration should be zero. Figure 5 reveals that the driving velocity is indeed proportional only to ΔT , indicating consistency with the YGB theory in equation (1). Even with different values of D, the magnitude of u d is constant if ΔT (≡ D∂T /∂x) is unchanged, as shown in Figure 5a. For instance, a ΔT of 2 K always results in u d ≈ −0.6, yielding u d /ΔT ≈ −0.3. According to Figure 5b, by including the other ΔT , the present result of u d /ΔT is slightly lower compared with that predicted by the YGB theory. This is because we simulated a droplet in a finite closed system surrounded by non-slip solid walls. A larger viscous drag might be exerted on the present moving droplet, and thus, u d is rather damped relative to the YGB theory. For a fixed ΔT and much larger D, the driving velocity is maintained up to D = 1 mm, at which the Reynolds number reaches the value of unity. We confirmed that u d decreases monotonically as D increases beyond 1 mm.
The magnitude of the driving force exerted on the droplet is shown in Figure 6. We estimated the magnitude F d by assuming that the driving force and viscous drag are balanced. The viscous drag can be determined based on the Hadamard-Rybczynski (HR) equation [23], which uses the fluid properties and the driving velocity: Then, the following combination of the YGB and HR equations can be drawn as Here, the drag coefficient ζ is a function of the blockage ratio D/L y and the viscosity ratio μ w /μ o [24]. In the present DNS, both D/L y and μ w /μ o are unchanged among the Temperature difference, (   Fig. 6. Dependencies of droplet driving force on droplet diameter and temperature gradient. The experimental results [11] and DNS data are derived from equation (11) and measured |u d |, while the theoretical values are from equation (12). In the experiment, the drag coefficient ζ is not constant because of the changing D with a fixed channel width L y .
test cases, and thereby, we employ ζ = 3 (refer to Chen et al. [24]). In the counterpart experiment [11] shown in Figure 6, ζ varied slightly as the droplet size changed relative to the fixed microchannel gap. To compare with experimental data, Figure 6 includes limited DNS results, which we obtained under the same condition of ∂T /∂x with the experiment for each D. This means that ΔT is not fixed in the figure. Hence, the nonlinear profiles observed in Figure 6a must not be simply a function of D 2 . The driving force F d is observed to increase rather linearly with respect to ΔT , as shown in Figure 6b, agreeing well with the theoretical prediction of equation (12). We confirm the driving force in the order of several nano-Newton (nN) under a temperature gap of several K, which reasonably agrees with the experimental evaluation. This may validate an insight of the experimental approach to determine the magnitude of the acting force, although some disparities occur between the experiment and DNS, in ranges of relatively large and small D.

Discussion
We would discuss the main mechanism of thermal droplet manipulation, based on additional DNS including hypothetical separation calculations regarding the normal and tangential forces of the interfacial tension. By omitting either f σ⊥ or f σ in the simulation, we may ignore one of the Laplace pressure or the Marangoni convection and simply extract the effects of the other component. For instance, a simulation without f σ⊥ would provide no Laplace pressure field; however, a Marangoni convection may be induced by a temperature gradient via the term f σ . Three cases of DNS considering f σ⊥ only, f σ only, and both components are performed in order for us to compare the separate contributions of the Laplace pressure and Marangoni convection. Although thermofluid dynamics including the present system are nonlinear physical phenomena, this hypothetical calculation would not misguide our understanding of photothermal manipulation because of the dominance of thermal conduction and the very low Reynolds and Weber numbers (refer to Tab. 2). Furthermore, we performed DNSs of a water droplet in oil in addition to the previously discussed oil droplet in water. The physical model of the water droplet in oil has similar geometry as in Figure 1; however, the fluid positional relationship is changed inside out.
The driving velocity and flow fields we obtained from the hypothetical DNSs are shown in Figures 7 and 8  equation (10). Its positive value corresponds to a migration toward the hot side. Only one case (without f σ⊥ for oil droplet) results in a positive driving velocity, while the other cases lead to movements in a common direction (toward the cold side).
First, let us focus on the magnitude of driving velocity and its temporal variation shown in Figure 7. Obviously, the driving velocity, by considering only f σ⊥ , is greater than that with only f σ . In the former case, |u d | is roughly four times larger than in the latter case. From this result, it can be conjectured that the dominant contribution for thermally manipulated droplet should be the normal force f σ⊥ , which induces directly the Laplace pressure, not the Marangoni convection. It should be noted that the oscillation that occurs when we consider f σ⊥ , as observed in Figure 7, should be a unphysical phenomenon similarly to the spurious current. We have confirmed that the amplitude and wave length of the oscillation reduced by adopting fine grids; however, the mean magnitude of u d did not change significantly, as shown in Figure 3. In contrast, the DNS with only f σ gives a smooth time variation of u d and reaches an equilibrium state before, at least, t = 10 ms. As mentioned above, it is interesting to note that the sign of u d , induced by f σ only, is reversed between the oil droplet and the water droplet. In particular, the sign of u d for the oil droplet is opposite to the driving direction in the actual case (where both f σ⊥ and f σ are considered). Therefore, we must regard the tangential component f σ as the driving force of the Marangoni convection, but not necessarily as a force driving the oil droplet. Figure 8 visualizes flow fields in equilibrium state, in the x-z plane across the center of the droplet, under the condition of D = 40 μm and ΔT = 2 K. The white contour corresponds to the water, and blue the oil. Within the oil droplet shown in Figures 8a1 and 8c1, the flow with the same order of the driving velocity can be observed. However, Figure 8b1 exhibits an opposite flow, which is in the same direction as the Marangoni convection at the interface. On the other hand, the flow within the "water" droplet develops in the same direction with the droplet migration, while the Marangoni convection toward the positive x direction actually occurs at the interface: see Figure 8b2. If we plot the relative velocity against the droplet, we would clearly observe circulating motions inside the droplet (not shown here). Owing to the matching direction of the two flows induced by f σ⊥ and f σ , the resultant "actual" flow inside the water droplet is larger than each flow, the obtained u d for the water droplet is also consistent with the YGB theory, as shown in Figure 8c2. From these results, we conclude that the Marangoni convection would play a role in either inhibition or enhancement of a droplet migration.
In Figure 9, we depicted the mechanisms of three different driven flows inside a drop: the Laplace-pressure-induced flow with temperature gradient, large-volumetric flow by Marangoni convection, and small-volumetric one. These schematics of flow  fields are views from a reference frame/coordinate at rest, similar to Figure 8. Figure 9a explains that the normal force f σ⊥ , the strength of which depends on the temperature, would produce a non-uniform Laplace pressure field and pressureinduced flow inside the droplet. This Laplace-pressure-induced flow causes a droplet migration toward the cold side. Figures 9b and 9c show that the Marangoni convection by the tangential force f σ produces a "boundary layer" of velocity in each side of the water-oil interface. Because of the large difference in viscosity, the induced flow volume in the oil side is greater than that in the water side. This difference in "viscous boundary layer" may cause the opposite results between oil and water droplets. It should be noted that in the cases of Figures 9b and 9c, the droplet cannot maintain a spherical form owing to the ignored normal force, i.e., without the Laplace pressure.
In such a case without the normal force, the initially-spherical droplet becomes a pancake-like shape after a long period of simulation.

Summary
In the framework of three-dimensional DNS, we demonstrated the oil-droplet migration in water toward the cold side with a positive ∂σ/∂T . The obtained driving velocity and force of the droplet manipulation were in good agreement with the experimental data and theoretical prediction. With fixed ΔT , the driving velocity is constantly independent of the droplet size, at least, in the present range (D = 30-100 μm). We found that in photothermal manipulation of microdroplet, the non-uniform Laplace pressure field would induce a dominant flow for the droplet migration, while the Marangoni convection at the droplet interface would work on the migration positively or negatively depending on the fluid property combination.
In this study, we determined indirectly the integrated driving force based on its balancing with the viscous drag. To measure locally the force and avoid the spurious current, we may adopt a high-order VoF method [25,26] in a future study.