Numerical Study of Droplet Dynamics on a Solid Surface with

May 23, 2019 - LEARN ABOUT THESE METRICS ... An increase in Re is able to increase not only droplet deformation but .... We for the first time explore...
0 downloads 0 Views 2MB Size
Download PDF
Article Cite This: Langmuir 2019, 35, 7858−7870

pubs.acs.org/Langmuir

Numerical Study of Droplet Dynamics on a Solid Surface with Insoluble Surfactants Jinggang Zhang,† Haihu Liu,*,† and Yan Ba‡ †

School of Energy and Power Engineering, Xi’an Jiaotong University, 28 West Xianning Road, Xi’an 710049, China School of Astronautics, Northwestern Polytechnical University, 127 West Youyi Road, Xi’an 710072, China

Downloaded via KEAN UNIV on July 21, 2019 at 05:34:10 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

‡

ABSTRACT: Surfactants are widely used in many industrial processes, where the presence of surfactants not only reduces the interfacial tension between fluids but also alters the wetting properties of solid surfaces. To understand how the surfactants influence the droplet motion on a solid surface, a hybrid method for interfacial flows with insoluble surfactants and contact-line dynamics is developed. This method solves immiscible two-phase flows through a lattice Boltzmann colorgradient model and simultaneously solves the convection− diffusion equation for surfactant concentration through a finite difference method. In addition, a dynamic contact angle formulation that describes the dependence of the local contact angle on the surfactant concentration is derived, and the resulting contact angle is enforced by a geometrical wetting condition. Our method is first used to simulate static contact angles for a droplet resting on a solid surface, and the results show that the presence of surfactants can significantly modify surface wettability, especially when the surface is more hydrophilic or more hydrophobic. This is then applied to simulate a surfactantladen droplet moving on a substrate subject to a linear shear flow for varying effective capillary number (Cae), Reynolds number (Re), and surface wettability, where the results are often compared with those of a clean droplet. For varying Cae, the simulations are conducted by considering a neutral surface. At low values of Cae, the droplet eventually reaches a steady deformation and moves at a constant velocity. In either a clean or surfactant-laden case, the moving velocity of the droplet linearly increases with the moving wall velocity, but the slope is always higher (i.e., the droplet moves faster) in the surfactantladen case where the droplet exhibits a bigger deformation. When Cae is increased beyond a critical value (Cae,c), the droplet breakup would happen. The presence of surfactants is found to decrease the value of Cae,c, but it shows a non-monotonic effect on the droplet breakup. An increase in Re is able to increase not only droplet deformation but also surfactant dilution. The role of surfactants in the droplet behavior is found to greatly depend upon the surface wettability. For a hydrophilic surface, the presence of surfactants can decrease the wetting length and enables the droplet to reach a steady state faster; while for a hydrophobic surface, it increases the wetting length and delays the departure of the droplet from the solid surface.

■

INTRODUCTION The motion of droplets adhering to a solid surface is of both fundamental and practical importance because of their wide range of applications, such as remediation of nonaqueousphase liquids,1−3 enhanced oil recovery (EOR),4−6 and droplet-based microfluidics.7−9 In these applications, surfactants are often present as impurities or are added deliberately to the bulk fluid, which preferentially absorb at the interfaces and can critically alter the dynamics of multiphase flow systems. Take the EOR by chemical flooding as an example; the addition of surfactants can significantly reduce the interfacial tension between oil and water and alter wetting properties of reservoir rocks,4,5 thereby influencing the flooding pattern and performance. As a complement to theoretical and experimental studies, numerical simulations allow access to complete information about flowfield, interface, and surfactant concentration, providing important insight into the droplet motion on a solid surface with surfactants, and thus, helping design and optimize these multiphase flow systems. © 2019 American Chemical Society

Numerical simulation of contact-line dynamics with surfactants in a multiphase flow system is a challenging task. The surfactant concentration at the interface alters the interfacial tension locally.10 The non-uniform surfactant concentration creates non-uniform interfacial tension forces and Marangoni stresses, which in turn affect the flowfield in a complicated manner. Meanwhile, the flowfield will influence the surfactant distribution. The interaction between surfactants and the flowfield is highly nonlinear. Numerically, the surfactant evolution equation must be solved together with the hydrodynamic equations at moving and deforming interfaces that may undergo topological changes such as breakup and coalescence. In addition, when the droplet touches a solid substrate, the contact angle should dynamically vary with the local surfactant concentration and the moving contact lines (MCLs) need to be properly dealt with, which are Received: February 19, 2019 Revised: April 22, 2019 Published: May 23, 2019 7858

DOI: 10.1021/acs.langmuir.9b00495 Langmuir 2019, 35, 7858−7870

Article

Langmuir

length and time scales. Thus, the LB method is particularly useful for simulation of multicomponent multiphase flows. Many multicomponent multiphase LB models have been proposed, and they can be classified into four major types: color-gradient model,19−21 phase-field-based model,22−25 interparticle-potential model, 26−29 and mean-field theory model.30−32 Among these models, the color-gradient model was extensively used to simulate immiscible multiphase flow problems33−35 because of its strengths such as low spurious velocities, high numerical accuracy, strict mass conservation for each fluid, and good numerical stability for a broad range of fluid properties. In particular, it was extended to model the thermocapillary flows36,37 and contact-line dynamics in the absence of surfactants.34,38 With the aid of the color-gradient model, Liu et al.33 recently developed a hybrid method for interfacial flows with insoluble surfactants. In this method, the immiscible two-phase flows are solved using an improved LB color-gradient model that incorporates dynamic interfacial tension and Marangoni stress, while the convection−diffusion equation which describes the evolution of surfactant concentration in the entire fluid domain is solved using a finite difference (FD) method. The LB and FD simulations are coupled through an equation of state (EOS), which describes how interfacial tension varies with surfactant concentration. However, such a hybrid method can only simulate surfactant dynamics with droplets suspended in a carrier fluid, away from the solid wall. In this work, we will present a hybrid LB−FD method for the simulation of interfacial flows with both insoluble surfactants and contact-line dynamics. The hybrid method is essentially the same as the one by Liu et al.,33 except that the original Bhatnagar−Gross−Krook (BGK) collision operator is replaced by its multiple-relaxation-time (MRT) counterpart. In addition, a dynamic contact angle formulation that describes the dependence of the local contact angle on the surfactant concentration is derived, and the resulting contact angle is enforced through the geometrical formulation of Ding and Spelt39 with implementation following our recent work.34 The developed hybrid method is first used to investigate the effect of surfactants on the static contact angle for a droplet resting on a solid surface. This is then used to simulate a surfactantladen droplet moving on a substrate subject to a linear shear flow for varying effective capillary number (Cae), Reynolds number, and surface wettability, where the results are often compared with those of a clean droplet. We for the first time explore how the presence of surfactants influences the droplet motion, deformation, and breakup for varying Cae.

still a pending physical problem to be resolved, posing an additional numerical challenge. Because of the challenges mentioned above, there have been only a few studies devoted to interfacial flows with both surfactants and contact-line dynamics. Lai et al.11 presented an immersed boundary method for simulating the motion of a droplet on a solid substrate under the influence of insoluble surfactants. The unbalanced Young force, which results from the deviation of the contact angle from its static value, was incorporated into the momentum equations as a driving force for the contact-line motion. Xu and Ren12 developed a level-set method for studying the droplet deformation with insoluble surfactants on a horizontal surface and the detachment of a pendant droplet from a wall under gravity. A contact angle condition, which relates the unbalanced Young force to the slip velocity at the contact line, was introduced to account for the effect of surfactants on the surface wettability. Based on the domain decomposition for representing the interface, af Klinteberg et al.13 developed an explicit Eulerian method to simulate two-phase flows with insoluble surfactants and contact line dynamics in 2D. To drive the movement of the contact line, they defined the tangential velocity in the vicinity of the contact line and used the Navier slip condition on the walls away from the contact line. Based on thermodynamic principles, Zhang et al.14 derived a continuous model for the dynamics of two immiscible fluids with MCLs and insoluble surfactants. The condition for the dynamic contact angle was obtained from the consideration of energy dissipations. Yon and Pozrikidis15 studied the effect of insoluble surfactants on a droplet attached to a plane wall and subjected to an overpassing shear flow in the limit of a vanishing Reynolds number. They assumed that the contact line retains a circular shape during the droplet deformation. To simplify the problem, they pinned the contact line on the wall and, thus, did not account for the sliding motion. Their numerical results revealed that the overall behavior of the droplet is determined by the magnitudes of the capillary number and the viscosity ratio, and also by the sensitivity of the interfacial tension to variations in surfactant concentration. In addition, an arbitrary Lagrangian−Eulerian finite element scheme for computations of soluble surfactant droplet impingement on a horizontal surface was presented by Ganesan.16 In his work, the dynamic contact angle was expressed as a function of the static contact angle and the local surfactant concentration. Clearly, these studies are all based on traditional computational fluid dynamics (CFD) methods, which solve the conservation equations of macroscopic properties (e.g., mass and momentum) numerically. However, these methods either are not suitable for dealing with droplets undergoing topological changes such as breakup and coalescence or require an unphysical re-initialization process to represent the interface, which may violate mass conservation for each fluid. Moreover, an empirical slip model with a slip length at the molecular scale often has to be introduced in these methods to avoid stress singularities at MCLs.17,18 By contrast, the lattice Boltzmann (LB) method, as a mesoscopic numerical method, has gained great success in simulating complex fluid flows in the past few decades. This possesses several attractive advantages over traditional CFD methods, such as intrinsic parallelism, simplicity of coding, and capability of handling complex geometries. Besides, its kinetic nature allows the simple incorporation of microscopic physics without the constraints of molecular dynamics simulations of

■

NUMERICAL METHOD In this section, we present a numerical approach for the simulation of contact line dynamics and insoluble surfactants. This approach is built upon a recently developed hybrid method, which describes the immiscible two-phase flows through an improved LB color-gradient model and solves the surfactant concentration transport through a FD method. In addition, a wetting boundary condition is developed to model the fluid−surface interactions, in which a dynamic contact angle formulation that describes the dependence of the local contact angle on the surfactant concentration is derived. Improved LB Color-Gradient Model for Immiscible Two-Phase Flows. We adopt the MRT color-gradient model for the immiscible two-phase flows with variable interfacial tension, which is developed on the basis of the work by Liu et 7859

DOI: 10.1021/acs.langmuir.9b00495 Langmuir 2019, 35, 7858−7870

Article

Langmuir al.36 In this model, two immiscible fluids, namely red and blue fluids, are considered and their distribution functions are f Ri and f Bi , where the subscript i denotes the ith direction of the discrete velocity. The total distribution function is defined by f i = f Ri + f Bi . This model consists of three steps in each time step, that is, collision, recoloring, and streaming. The recoloring and streaming steps are identical to those described by Liu et al.,36 so they are not shown here. For the collision step, the total distribution function evolves as fi† (x , t ) = fi (x , t ) + Ωi(x , t ) + Fi ̅

where νR and νB are the viscosities of red and blue fluids, and ρN is the phase field function responsible for differentiating different fluids, which is defined by ρ N (x , t ) =

1 y i F̅ = M−1jjjI − SzzzMF̃ 2 { k

(1)

where x and t are the position and time, respectively, is the post-collision total distribution function, Ωi is the collision operator, and F̅ i is the forcing term. Instead of the BGK approximation used in previous studies,33,38,40 the MRT scheme is adopted for the collision operator. The MRT scheme has a number of advantages over the BGK approximation in simulating multiphase flow problems, including enhanced numerical instability,41,42 minimization of spurious velocities,43,44 and reproducing the correct contact angle. With the MRT scheme, the collision operator Ωi is given by

In eq 9, the fluid velocity is defined by ρu(x,t) = ∑i f i(x,t)ei + FS(x,t)δt/2, and the interfacial force FS includes not only the interfacial tension force but also the Marangoni stress caused by the non-uniform distribution of surfactants along the interface. Using the concept of continuum surface force, the interfacial force FS can be expressed as33

(2)

1 FS = − σ0[1 + E0 ln(1 − ψ /ψ∞)]κ∇ρ N 2 σE 1 − |∇ρ N | 0 0 [∇ψ − (n ·∇ψ)n] 2 ψ∞ − ψ

σ(ψ ) = σ0[1 + E0 ln(1 − ψ /ψ∞)]

(4)

(5)

∂φ(x) 1 = 2 ∑ ωiφ(x + eiδt )eiα ∂xα cs i

where s1, s2, and s4 are three independently adjustable relaxation rates, and they can be freely adjusted to ensure numerical stability and low spurious velocities. In this study, these free parameters are chosen as s1 = 1.63, s2 = 1.14, and s4 = 1.92, which were demonstrated to produce stable and accurate results in the simulation of contact line dynamics.46 For the sake of simplicity, red and blue fluids are assumed to have similar densities. To account for different viscosities of the two fluids, if needed, the kinematic viscosity ν is determined by the harmonic mean as47 1 + ρN 1 − ρN 1 = + R ν 2ν 2ν B

(11)

has been introduced to describe the relation between the interfacial tension σ and surfactant concentration ψ. In the above two equations, n is the interfacial unit normal vector defined by n = −∇ρN/|∇ρN|, κ is the local interface curvature related to n by κ = −∇s·n, ψ∞ is the surfactant concentration at the maximum packing, σ0 is the interfacial tension of a clean interface (i.e., ψ = 0), and E0 is the elasticity number. As previously performed in the literature,33,48 the initial average surfactant concentration ψin is given by a dimensionless surfactant coverage, that is, xin = ψin/ψ∞. To minimize the discretization errors, the compact FD is used to evaluate the partial derivatives involved in the normal vector and curvature calculations. Taking variable φ as an example, its partial derivatives can be evaluated using

where s7 and s8 are related to the kinematic viscosity ν of the fluid mixture by 1 1y i , and ν = jjjτ − zzzcs 2δt 2{ τ k

(10)

in which the Langmuir EOS, that is,

where ρ = ρR + ρB is the total density, with the superscripts “R” and “B” referring to the red and blue fluids, respectively; ei is the discrete velocity, ωi is the weighting factor, and cs is the speed of sound. In the two-dimensional 9-velocity lattice model, the speed of sound cs = δx /( 3 δt ) = 1/ 3 where the lattice spacing δx and the time step δt are both taken as 1, and the explicit expressions for ei, ωi, and M can be found in Lallemand and Luo.45 The diagonal relaxation matrix S reads as34

s7 = s8 =

(8)

where I is a 9 × 9 unit matrix, F̅ = [F̅ 0, F̅ 1, F̅ 2, ..., F̅ 8]T, and F̃ = [F̃ 0, F̃ 1, F̃ 2, ...,F̃ 8]T is given by ÄÅ ÉÑ ÅÅÅ ei − u (ei ·u)ei ÑÑÑ ÑÑ·FSδt Fi ̃ = ωiÅÅÅ + ÅÅÅ cs 2 cs 4 ÑÑÑÑÖ (9) Ç

where M is the transformation matrix and S is the diagonal relaxation matrix. f eq i is the equilibrium distribution function of f i and can be obtained by a second-order Taylor expansion of the Maxwell−Boltzmann distribution with respect to the local fluid velocity u: ÅÄÅ ÑÉ Å e ·u (e ·u)2 u 2 ÑÑÑÑ f ieq = ρωiÅÅÅÅ1 + i 2 + i 4 − Ñ ÅÅÅ 2cs cs 2cs 2 ÑÑÑÑÖ (3) Ç

S = diag[0, s1 , s2 , 0, s4 , 0, s4 , s7 , s8]

(7)

The forcing term in eq 1 contributes to the mixed interfacial regions and generates an interfacial force FS. In the MRT framework, the forcing term reads as44

f †i

Ωi(x , t ) = −(M−1SM)ij (f j (x , t ) − f jeq (x , t ))

ρ R (x , t ) − ρ B (x , t ) ρ R (x , t ) + ρ B (x , t )

(12)

FD Method for Surfactant Transport. The surfactant evolution equation of the diffuse-interface form, originally proposed by Teigen et al.,49 is used to describe the surfactant transport, as recently done by Liu et al.33 This not only allows for the solution of surfactant concentration in the entire fluid domain without additional extension or initialization procedures, but also offers great simplicity in dealing with topological changes such as interface rupturing and merging. By taking the Dirac function δΓ as |∇ρN|/2, the surfactant evolution equation of the diffuse-interface form can be written as33

(6) 7860

DOI: 10.1021/acs.langmuir.9b00495 Langmuir 2019, 35, 7858−7870

Article

Langmuir

contact angle. According to Young’s equation, the contact angle θ, in the presence of surfactants, satisfies

∂t(|∇ρ N |ψ ) + ∇·(ψ |∇ρ N |u) = Ds[(∇|∇ρ N |) ·∇ψ + |∇ρ N |∇2 ψ ]

(13)

σ0[1 + E0 ln(1 − ψ /ψ∞)]cos θ = σSB − σSR

where Ds is the surfactant diffusivity. It can be shown using the asymptotic analysis that eq 13 converges to the most used surfactant evolution equation of the sharp-interface form as the interface thickness tends to zero.50 Following Liu et al.,33 eq 13 is solved by the FD method. Specifically, a modified Crank−Nicolson scheme is used for the time discretization, and the resulting spatial derivatives are all discretized using the standard central difference schemes except the convection term u·∇(ψ|∇ρN|), which is discretized using a third-order weighted essentially non-oscillatory (WENO) scheme.33,51,52 In order to apply the third-order WENO scheme, at a fluid node (I,J), we need to know the values of u, ψ, and |∇ρN| at the nodes (I − 2:I + 2,J) in the xdirection and (I,J − 2:J + 2) in the y-direction. This means that two layers of solid nodes (also termed ghost nodes) neighboring to the solid wall need to be considered. Without loss of generally, we take the bottom wall, which is located halfway between the first layer of solid nodes (y = 0) and the first layer of fluid nodes (y = 1), as an example to explain how to obtain the values of u, ψ, and |∇ρN| at the ghost nodes, that is, y = −1 and y = 0. Assuming that the velocity continuously varies across the no-slip wall, one easily obtains ux,−1 = −ux,2 and ux,0 = −ux,1. For the surfactant concentration, nopenetration condition is imposed on the solid wall, which implies ψx,−1 = ψx,2 and ψx,0 = ψx,1. In addition, we simply set |∇ρN|x,−1 = |∇ρN|x,2 and |∇ρN|x,0 = |∇ρN|x,1 to ensure the continuity of ψ|∇ρN| across the wall, thus minimizing the discretization errors. Finally, the resulting linear system for ψt+δt from time and spatial discretization is solved by the successive over relaxation method with a relaxation factor of 1.2. Wetting Boundary Condition. The contact angle is a quantitative measure of the surface wettability, and represented by the angle made by an interface with a solid surface. As shown in Figure 1, a droplet of the red fluid, surrounded by the

(15)

where we have replaced the local interfacial tension σ by σ0[1 + E0 ln(1 − ψ/ψ∞)]. By combining eqs 14 and 15, we can obtain a dynamic contact angle formulation that describes the dependence of the contact angle on the local surfactant concentration cos θ(ψ ) =

cos θ0 1 + E0 ln(1 − ψ /ψ∞)

(16)

Once the local contact angle θ is obtained on the solid surface, it can be enforced through the geometrical formulation proposed by Ding and Spelt39 iπ y n w ·∇ρ N = −tanjjj − θ zzz|τw ·∇ρ N | k2 {

(17)

where nw is the unit normal vector to the wall pointing toward the fluids and τw is the unit vector tangential to the wall. Following Liu et al.,34 the enforcement of eq 17 is realized through one layer of solid nodes, which is a half lattice spacing away from the wall. Again, we consider the case of the bottom wall, where the no-slip boundary condition is imposed by the halfway bounce-back scheme.53 Upon discretization of eq 17, one has iπ y ρxN,0 = ρxN,1 + tanjjj − θ zzz|τw ·∇ρ N |δx 2 k {

(18)

where the tangential component of ∇ρN on the wall is determined by the following extrapolation scheme54 τw ·∇ρ N = 1.5∂xρ N |x ,1 − 0.5∂xρ N |x ,2

(19)

In eq 19, the derivatives on the right-hand side can be easily evaluated using the second-order central difference approximation. Using the values of ρNx,0, we are then able to compute the partial derivatives of ρN at all the fluid nodes through eq 12. In such a way, the dynamic contact angle is implicitly imposed on the solid surface.

■

RESULTS AND DISCUSSION As a preliminary study, we first simulate a droplet with and without surfactants sitting on a solid surface in a quiescent flow, and the obtained results are compared to examine whether the presence of surfactants can significantly modify the contact angle. We then simulate a surfactant-laden droplet on a substrate subject to a simple shear flow for varying capillary number, Reynolds number, and surface wettability, where the simulation results are often compared with those in the case of a clean droplet to show the role of surfactants. A Droplet Resting on a Solid Surface with Different Contact Angles. Consider a semi-circular droplet (red fluid) with a radius R = 100 initially deposited on the bottom wall. The size of the computational domain is 6R × 3R. The periodic boundary condition is used in the horizontal direction, while the proposed wetting boundary condition is imposed at the top and bottom walls in the vertical direction. Both fluids are assumed to have equal densities of 1 and equal viscosities of 0.25. The Langmuir EOS, eq 11, is employed, in which E0 = 0.5 and σ0 = 4 × 10−3, and the surfactant diffusivity is taken as Ds = 0.03. Two different values of surfactant

Figure 1. Schematic diagram of the contact angle for a droplet of the red fluid on an ideal surface. θ is the equilibrium contact angle measured from the red fluid side; σ is the interfacial tension between red and blue fluids, while σSR (σSB) is the interfacial tension between the red (blue) fluid and the solid surface.

second (blue) fluid, is placed on an ideal surface. In the absence of surfactants, the static contact angle θ0 is related to the interfacial tensions by Young’s equation σ0 cos θ0 = σSB − σSR

(14)

where σSR or σSB is the interfacial tension between the red or blue fluid and the solid surface. In the surfactant-laden case, for the sake of simplicity, we assume that the surfactants are only present at the fluid−fluid interface, which suggests that the presence of surfactants does not affect the values of σSR and σSB, but decreases the interfacial tension σ, thus, modifying the 7861

DOI: 10.1021/acs.langmuir.9b00495 Langmuir 2019, 35, 7858−7870

Article

Langmuir coverage are considered, that is, xin = 0 and 0.3, which correspond to the clean droplet and the surfactant-laden droplet, respectively. For either the clean droplet or the surfactant-laden droplet, the static contact angle θ0 is varied from 45° to 135°. Each simulation is run with a zero initial velocity throughout the domain until the droplet reaches the steady state. Figure 2 shows the equilibrium droplet shapes for both clean and surfactant-laden droplets at several typical values of θ0. In

droplets increases, with the biggest difference of around 11°. These results suggest that the presence of surfactants can significantly alter the surface wettability, even for a static problem. A Droplet on a Solid Surface Subject to a Shear Flow. In this section, we consider a droplet attached on a solid surface subjected to a linear shear flow. As shown in Figure 3,

Figure 3. Schematic diagram of a droplet meniscus on a solid surface subject to a linear shear flow.

two parallel walls are separated by a distance 2H. The bottom wall is stationary, and the shear flow is introduced to the system by moving the top wall towards the right with a constant velocity uw. The periodic boundary conditions are used in the horizontal direction, while the proposed wetting boundary condition is imposed at the top and bottom walls in the vertical direction. Initially, a droplet of the circular segment (red fluid) with an area Ad = πH2/2 and a static contact angle θ0 is deposited on the bottom wall. The droplet is assumed to have equal density and viscosity to those of the matrix (blue) fluid, and the velocity field is initialized by the imposed shear flow. In the presence of surfactants, the droplet behavior can be characterized by several important dimensionless numbers, including the capillary number (Ca), Reynolds number (Re), surface Péclet number (Pe), which are defined as34,55

Figure 2. Equilibrium droplet shapes for xin = 0 (represented by black lines) and xin = 0.3 (represented by red lines) at several typical values of θ0. The horizontal coordinate relative to the center of the computational domain is normalized by the initial droplet radius R.

this figure, the clean and surfactant-laden droplets are represented by black and red lines, respectively. In comparison with the clean droplets, the surfactant-laden droplets spread more, thus, having smaller equilibrium contact angles with the hydrophilic surfaces (θ0 < 90°), whereas they recoil more, thus, having larger equilibrium contact angles with the hydrophobic surfaces (θ0 > 90°). The difference in equilibrium droplet shapes for clean and surfactant-laden droplets can be explained as follows. The addition of surfactants leads to 0 < 1 + E0 ln(1 − ψ/ψ∞) < 1, so one can obtain from eq 16 that cos θ > cos θ0 > 0 and further θ < θ0 < 90° for the hydrophilic surfaces, and that cos θ < cos θ0 < 0 and further θ > θ0 > 90° for the hydrophilic surfaces. Using the measured droplet height and base diameter, we are able to compute the equilibrium contact angles, which are shown in Table 1. As expected, the simulated

Ca =

θ0 xin = 0 xin = 0.3

45

60

90

120

135

44.52 35.55

59.80 53.52

89.96 89.96

120.44 127.11

136.21 147.18

Re =

γeH ̇ , νB

Pe =

γeH ̇ Ds (20)

Here, e is the initial height of the droplet and γ̇ is the shear rate calculated using γ̇ = uw/(2H). It is known that a reduction of the interfacial tension and, thus, an increase of the effective capillary number, because the presence of surfactants can increase the droplet deformation. In order to rule out the effect of average surfactant concentration on reducing the interfacial tension, the effective capillary number

Table 1. Equilibrium Contact Angles of Clean and Surfactant-Laden Droplets at Various θ0a θ

ρ B ν Bγė , σ0

Cae =

a

All angles are shown in degrees.

ρ B ν Bγė Ca = σ0[1 + E0 ln(1 − xin)] 1 + E0 ln(1 − xin) (21)

is used instead of Ca in both clean and surfactant-laden systems. We compare the simulation results of clean and surfactant-laden droplets at the same value of Cae, and the difference in the simulation results is attributed to two factors: (1) non-uniform effects from non-uniform capillary pressure at the droplet interface and Marangoni stresses along the interface, and (2) surfactant dilution due to the increase of the interfacial length. The size of the computational domain is set as L × 2H = 1500 × 200, where the domain length L is large enough to avoid the occurrence that the droplet touches itself via the periodic boundary conditions before breakup. In

equilibrium contact angles are all in good agreement with their analytical solutions for the clean droplets, with the maximum difference as low as 1.21°. Consistent with the observations in Figure 2, the presence of surfactants is found to decrease the equilibrium contact angle for θ0 < 90°, but increase the equilibrium contact angle for θ0 > 90°. Note that such findings were also reported in recent studies, for example, by Xu and Ren12 and Lai et al.11 In addition, it is noticed that as θ0 is increased or decreased from 90°, the difference in the equilibrium contact angle between clean and surfactant-laden 7862

DOI: 10.1021/acs.langmuir.9b00495 Langmuir 2019, 35, 7858−7870

Article

Langmuir the Langmuir EOS, the elasticity number E0 = 0.5 and σ0 = 4.0 × 10−3. Unless otherwise stated, the simulations are run for Re = 1, Pe = 5, and θ0 = 90° and two different values of the surfactant coverage, that is, xin = 0 and 0.3. Note that in the clean case (xin = 0), the constant interfacial tension σ should be taken as σ0[1 + E0 ln(1 − xin)] = 4 × 10−3[1 + 0.5 ln(1 − 0.3)] = 3.29 × 10−3 to achieve the same effective interfacial tension as in the surfactant-laden case (xin = 0.3). Because the grid size may influence the simulation results of the present method, it is important to minimize this numerical error. As an example, we examine the effect of the grid size on the numerical results for Cae = 0.15 and xin = 0.3. Table 2

deformation and moves faster than the clean one. The increased deformation in the surfactant-laden case can be explained by the surfactant concentration distributions shown in Figure 5. In the presence of surfactants, the shear flow

Table 2. Steady State Results of Different Grid Sizes for Cae = 0.15 and xin = 0.3 grid size

S/H

ud/uw

ψ*min

ψ*max

L × 2H = 750 × 100 L × 2H = 1500 × 200 L × 2H = 2250 × 300

3.3044 3.3391 3.3487

0.1333 0.1280 0.1289

0.1216 0.1160 0.1140

1.8947 1.9134 1.9542

Figure 5. The distribution of normalized surfactant concentration (ψ*) at γ̇t = 30 for (a) Cae = 0.1, (b) Cae = 0.15, and (c) Cae = 0.2. The x coordinate relative to the center of the computational domain x0 is normalized by H. Note that the droplet has reached a stable deformation at γ̇t = 30 in these cases.

shows the simulation results in the steady state for three different grid sizes, that is, L × 2H = 750 × 100, L × 2H = 1500 × 200, and L × 2H = 2250 × 300. In this table, S and ud are the arc length and moving velocity of the equilibrium droplet, and ψmin * and ψmax * are the minimum and maximum values of ψ*, where ψ* = ψ/ψin is the normalized surfactant concentration at the steady droplet interface (represented by ρN = 0). Clearly, the grid sizes of L × 2H = 1500 × 200 and L × 2H = 2250 × 300 produce nearly identical results (with a relative difference below 2.5%), which confirms that the grid size of L × 2H = 1500 × 200 is able to provide gridindependent numerical results. The droplet deformation is first studied by increasing Cae from 0 with an increment of 0.05. In the deformation mode, the droplet can eventually reach a steady shape and slides on the solid surface at a constant velocity. It is found that the droplet deformation occurs for Cae ≤ 0.2 in the surfactantladen case, but occurs for Cae ≤ 0.25 in the clean case. Figure 4 shows the snapshots of the droplet sliding on the solid surface for Cae = 0.1, 0.15, and 0.2, where the clean and surfactantladen droplets are represented by black and red lines, respectively. It is seen that the droplet deformation increases with Cae in both clean and surfactant-laden cases. For each value of Cae, the surfactant-laden droplet exhibits a bigger

sweeps the surfactants from the rear of the droplet toward the front, resulting in a high surfactant concentration and low interfacial tension at the front. Because the interfacial tension acts to resist against the droplet deformation, the low interfacial tension at the front will contribute to the increased droplet deformation in the surfactant-laden case. A similar result was also obtained for a droplet subject to a simple shear flow and away from the solid wall.33 Figure 6a plots the normalized surfactant concentration ψ* at the steady droplet interface as a function of the arc length s for various values of Cae. Note that the arc length s, measured clockwise along the droplet interface from the receding contact point (i.e. the intersection point between the rear interface of the droplet and the solid surface) to the advancing contact point (i.e. the intersection point between the front interface of the droplet and the solid surface), is normalized by the initial droplet radius R. As Cae increases, ψ*min decreases in the entire range of Cae considered, and ψ*max increases for Cae ≤ 0.175. When Cae increases from 0.175 to 0.2, ψmax * does not increase but decreases, which is caused by the pronounced surfactant dilution because of interfacial stretching (see Figure 6b). In addition, we can notice that the minimum and maximum values of surfactant concentration are located close to the receding contact point and the advancing contact point, respectively. On the other hand, at high values of Cae, for example, Cae = 0.175 and 0.2, the position of the maximum surfactant concentration gradually moves away from the advancing contact point with an increase in Cae. This is because the increased droplet deformation leads to a big curvature of the droplet tip where the surfactants prefer to accumulate. When a droplet is suspended in a matrix fluid and subject to a shear flow, the droplet eventually deforms into more or less an ellipsoidal shape, which is usually characterized by the deformation parameter D, that is, D = (a − b)/(a + b), where a and b are the half-lengths of the major and minor axes of the droplet. However, the droplet does not exhibit an ellipsoidal shape because of the wall effects in the present study. Thus, the droplet deformation is characterized by the relative arc length Sr instead of the deformation parameter. The relative arc length

Figure 4. The snapshots of the droplet sliding on a solid surface at γ̇t = {20,30,40} for (a) Cae = 0.1, (b) Cae = 0.15, and (c) Cae = 0.2. For each Cae, the clean and surfactant-laden droplets are represented by the black and red lines, respectively. The x coordinate relative to the center of the computational domain x0 is normalized by H. 7863

DOI: 10.1021/acs.langmuir.9b00495 Langmuir 2019, 35, 7858−7870

Article

Langmuir

Figure 6. (a) The distribution of normalized surfactant concentration ψ* along the arc length s for the surfactant-laden droplet at various Cae, and (b) the relative arc length Sr vs the effective capillary number Cae for clean and surfactant-laden droplets. In (a), the inset shows the plot of the maximum surfactant concentration ψ*max and minimum surfactant concentration ψ*min as a function of Cae.

is defined as Sr = (S − S0)/S0, where S0 and S are the arc lengths of the initial droplet and the equilibrium droplet, respectively. It is clear that the greater the Sr is, the more the droplet will deform. Figure 6b shows the relative arc length as a function of Cae for both clean and surfactant-laden droplets. Consistent with the previous observations in Figure 4, the droplet deformation increases with Cae for both clean and surfactant-laden droplets, and the surfactant-laden droplet undergoes a bigger deformation than the clean one at the same value of Cae. In addition, it is seen that the difference in Sr between the clean droplet and the surfactant-laden droplet dramatically increases with Cae, which can be explained as follows. As mentioned before, the increased droplet deformation in the surfactant-laden case is attributed to the non-uniformity of surfactant distribution, and thus, the more non-uniform the surfactant distribution is, the bigger the difference in Sr between the clean droplet and the surfactantladen droplet will be. As shown in the inset of Figure 6a, the difference between ψ*max and ψ*min increases with Cae for Cae ≤ 0.175, suggesting an increased non-uniformity of the surfactant distribution. On the other hand, when Cae is increased from * and ψmin * seems to 0.175 to 0.2, the difference between ψmax decrease, but the difference in Sr between the clean droplet and the surfactant-laden droplet still increases. This is because the * and ψmin * does not account for the nondifference between ψmax uniformity of surfactants in the vicinity of the advancing contact point, that is, the interface from the advancing contact point to the location where ψ*max occurs; such a region is able to significantly increase the droplet deformation due to the relatively low interfacial tension and keeps growing with Cae (see Figure 6a). In the deformation mode, the droplet eventually slides on the solid surface at a constant velocity ud, and the droplet motion is quantified by the contact-line capillary number, which is defined as Cacl =

ρ B ν Bud . σ0[1 + E0 ln(1 − xin)]

Figure 7. Contact-line capillary number Cacl as a function of the effective capillary number Cae. The simulation data of the clean droplet and the surfactant-laden droplet are fitted separately using a linear relationship, and the resulting fitting lines are represented by the black and red lines, respectively.

of the droplet in the present study. However, previous results showed that the presence of surfactants retards the motion of the bubble rising under buoyancy56−58 and the droplet settling under gravity.59 So why do the surfactants play different roles in the present case and the bubble rising case? We shall answer this question in the following paragraph. We first consider the bubble rising under buoyancy, which is sketched in Figure 8a. When a bubble rises upward in a surfactant solution, the surfactants will be swept toward the

It should be noted

that ud equals to the contact-line velocity in the equilibrium state. Figure 7 shows the contact-line capillary number Cacl as a function of Cae. It is clearly seen that in either the clean or the surfactant-laden case, the contact-line capillary number linearly increases with Cae especially at low values of Cae, but the slope is always higher in the surfactant-laden case. The slope K refers to the ratio of Cacl to Cae, which can be further expressed as Ca u 2u K = Cacl = γdė = u d . By the linear fitting, we obtain the e

w

terminal velocities ud = 0.1125uw and ud = 0.1312uw for the clean and surfactant-laden droplets, respectively. This suggests that the presence of surfactants increases the moving velocity

Figure 8. Schematic diagrams of (a) a surfactant-laden bubble rising under buoyancy and (b) a surfactant-laden droplet sliding on a solid surface under shear. 7864

DOI: 10.1021/acs.langmuir.9b00495 Langmuir 2019, 35, 7858−7870

Article

Langmuir bottom of the bubble, where the interfacial tension is low. The non-uniform interfacial tension results in a Marangoni stress that is directed along the interface from a low interfacial tension region to a high interfacial tension region (see the dashed line with an arrow in Figure 8a). To balance the Marangoni stress, a viscous drag force is induced and exerted on the surface of the moving bubble, which is directed opposite to the Marangoni stress. Clearly, the induced net drag force is pointed downward, opposite to the rising direction of the bubble, so the bubble rising will be retarded in the presence of surfactants. We then turn to the present case, in which the droplet subjected to a shear flow moves to the right on the solid surface. As shown in Figure 8b, the surfactants are swept to the neighborhood of the advancing contact point (also see Figure 5), and the resulting Marangoni stress is directed anticlockwise along the droplet interface (see the dashed line with an arrow in Figure 8b). The viscous drag force, which is induced to balance the Marangoni stress, is pointed to the right, that is, in the same direction as the moving direction of the droplet, thus promoting the droplet motion. The wetting area or length is an important parameter influencing mass and energy transfer in many industrial applications, and it is defined as the contact area or length between the droplet and the surface. In this work, we introduce the relative wetting length lr, which is defined by the initial wetting length l0 and the steady-state wetting length l as lr = (l − l0)/l0, and quantify the effect of surfactants on the relative wetting length. In Figure 9, we plot the relative wetting length

Figure 10. Droplet shapes and streamlines for (a) clean and (b) surfactant-laden droplets at γ̇t = 30 and Cae = 0.1. The clean droplet and the surfactant-laden droplet are represented by the black and red lines, respectively.

visible inside the surfactant-laden droplet. The difference can be explained as follows: in the surfactant case, the Marangoni flow would induce an additional viscous flow inside the moving droplet, which is in the same direction as the droplet motion and, thus, suppresses the formation of the vortex. Consistent with our findings, a recent numerical study also showed that the recirculation present in the low-surfactant concentration droplet gradually weakens and even vanishes with increasing surfactant concentration during the droplet generation in a microfluidic T-junction.60 In addition, we notice more streamlines passing through the surfactant-laden droplet than the clean one, which suggests that the surfactant-laden droplet moves faster, in agreement with the results shown in Figure 7. As the effective capillary number increases above a critical value, Cae,c, the droplet continuously deforms without reaching a steady state, and finally breaks up into daughter droplets. Here, we investigate how the presence of surfactants influences the droplet breakup. First, we increase the effective capillary number to 0.25 and keep all the other parameters the same as in the last simulations. Figure 11 shows the snapshots of the clean and surfactant-laden droplets for Cae = 0.25. It is seen that the interface is stretched and eventually reaches a steady shape for the clean droplet; whereas for the surfactant-laden droplet, after the initial stretching, a neck forms and then gradually elongates and thins until the droplet finally breaks up into two unequal-sized parts: the big part migrates away from the wall while the small part keeps adhering to the wall. The difference indicates that the presence of surfactants promotes the droplet breakup, thus decreasing the critical capillary number for the droplet breakup. Specifically, the critical capillary number is changed from Cae,c > 0.25 for the clean droplet to 0.2 < Cae,c ≤ 0.25 for the surfactant-laden droplet (Cae,c should be greater than 0.2 because the surfactant-laden droplet can reach a steady deformation at Cae = 0.2, as previously shown in Figure 5). To understand the effect of surfactants on the process of droplet breakup, Figure 12 plots the time evolution of normalized surfactant concentration at the droplet interface for Cae = 0.25. In the early stage, the surfactants are swept by the shear flow to the downstream of the droplet (Figure 12b) where the interfacial tension is extremely low. The extremely low interfacial tension remarkably increases the local interface deformation, and thus, a neck is formed near the advancing contact point (see Figure 12c). In the necking stage, ψ*max decreases and ψmin * almost remains a constant (see Figure 12d,e), exhibiting a surfactant dilution. Although the surfactant dilution favors the prevention of the droplet deformation to some extent, the necking or elongation cannot be completely dismissed once the neck has reached a critical size. The above analysis shows that the decrease of Cae,c in the presence of surfactants is attributed to the increased droplet deformation, caused by the non-uniform effects associated with non-uniform capillary pressure, in the early stage.

Figure 9. Relative wetting length lr as a function of the effective capillary number Cae for clean and surfactant-laden droplets.

lr as a function of Cae. It is seen that in the clean case, the relative wetting length lr continuously increases with Cae. This is easy to understand because as Cae increases, the shear flow strengthens; in order to balance such a shear flow, an increased wall friction and, thus, an increased wetting length is needed. However, in the surfactant-laden case, the relative wetting length keeps nearly a constant except when Cae is close to the critical value for droplet breakup, for example, Cae = 0.2. This is because as Cae increases, the interfacial tension in the neighborhood of the advancing contact point decreases, leading to a decrease in the curvature radius of the local interface (see Figure 4) and, thus, offsetting the increase of the wetting length. Figure 10 plots the instantaneous streamlines and droplet shapes in both clean and surfactant-laden cases. Clearly, a single vortex is formed inside the clean droplet, but no vortex is 7865

DOI: 10.1021/acs.langmuir.9b00495 Langmuir 2019, 35, 7858−7870

Article

Langmuir

Figure 12. Distribution of the normalized surfactant concentration at the droplet interface for Cae = 0.25 at (a) γ̇t = 0, (b) γ̇t = 10, (c) γ̇t = 20, (d) γ̇t = 30, (e) γ̇t = 36.7, (f) γ̇t = 37.3, and (g) γ̇t = 40.

Figure 11. The snapshots of the droplet moving on a solid wall when subjected to a shear flow for Cae = 0.25 at (a) γ̇t = 0, (b) γ̇t = 10, (c) γ̇t = 20, (d) γ̇t = 30, (e) γ̇t = 36.7, (f) γ̇t = 37.3, and (g) γ̇t = 40. The clean and surfactant-laden droplets are represented by the black and red lines, respectively. The x coordinate relative to the center of the computational domain x0 is normalized by H.

this stage, the non-uniform effects dominate the droplet deformation due to big surfactant concentration gradients. As the droplet is quickly stretched, the surfactants exhibit a severe dilution at the interface (see Figure 14c−f), which weakens the non-uniform effects and becomes dominating the droplet deformation. Since the surfactant dilution increases the interfacial tension, the droplet deforms less than the clean one and finally the breakup is retarded in the presence of surfactants (see Figure 13c−g). In summary, the presence of surfactants is found to decrease the value of Cae,c, but it shows a non-monotonic effect on the droplet breakup: it may promotes or prevents the droplet breakup depending on the competition between non-uniform effects and surfactant dilution. When the non-uniform effects are dominant, the presence of surfactants promotes the droplet breakup; on the contrary, it prevents the droplet breakup. Another important parameter governing droplet behavior is the Reynolds number and its influence is investigated for a surfactant-laden droplet on a neutral surface at Cae = 0.15 and xin = 0.3. Three different values of Re are considered, that is, Re = 0.1, 1 and 10, and each simulation is run until a steady droplet deformation is reached (γ̇t = 20). Figure 16 shows the equilibrium droplet shapes and surfactant concentration distributions at the droplet interface. As Re increases, the droplet exhibits a larger deformation and the value of ψ*max decreases. We also quantify the droplet deformation and surfactant concentration distribution through Sr and ψmin * and ψ*max, and the corresponding results are listed in Table 3. It is seen that the droplet deformation increases but the values of ψmin * and ψmax * both decrease with increasing Re. The decrease in both ψmin * and ψmax * is because the surfactants become

Then, we further increase the effective capillary number to 0.3, and the corresponding snapshots for both clean and surfactant-laden droplets are shown in Figure 13. It can be seen in the early stage that a neck forms quickly for both droplets because of a strong shear flow, but the surfactant-laden droplet exhibits a bigger deformation than clean one (see Figure 13b). As the interface is stretched, the clean droplet later deforms more and finally breaks up earlier than the surfactant-laden droplet (see Figure 13c−e). In addition, we also notice that the clean droplet breaks up into three daughter droplets, known as ternary breakup, while the surfactant-laden droplet only breaks up into two daughter droplets, known as binary breakup (see Figure 13f,g). This suggests that the presence of surfactants can influence not only the breakup time but also the breakup mode. Finally, it should be noted that, as Cae continues to grow, for example, Cae = 0.35, the clean droplet always breaks up earlier than the surfactant-laden droplet (which is shown in Figure 15). As mentioned above, the clean droplet breaks up earlier than the surfactant-laden one, suggesting that the presence of surfactants prevents the droplet breakup at Cae = 0.3, which is contrary to the observation at Cae = 0.25. To clarify why the presence of surfactants prevents the droplet breakup, we plot the evolution of surfactant concentration for the surfactantladen droplet at Cae = 0.3 in Figure 14. Clearly, the surfactants first accumulates in the vicinity of the droplet tip where the interface curvature is relatively high (see Figure 14b), thus reducing the interfacial tension and leading to a more elongated droplet than the clean case (see Figure 13b). In 7866

DOI: 10.1021/acs.langmuir.9b00495 Langmuir 2019, 35, 7858−7870

Article

Langmuir

Figure 13. The snapshots of the droplet moving on a solid wall when subjected to a shear flow for Cae = 0.3 at (a) γ̇t = 0, (b) γ̇t = 10, (c) γ̇t = 20, (d) γ̇t = 25.1, (e) γ̇t = 25.7, (f) γ̇t = 28.3, and (g) γ̇t = 33. The clean and surfactant-laden droplets are represented by the black and red lines, respectively.

Figure 14. Distribution of the normalized surfactant concentration at the droplet interface for Cae = 0.3 at (a) γ̇t = 0, (b) γ̇t = 10, (c) γ̇t = 20, (d) γ̇t = 25.1, (e) γ̇t = 25.7, (f) γ̇t = 28.3, and (g) γ̇t = 33.

increasingly diluted as the droplet deformation significantly increases. The above studies on the droplet deformation and breakup all focus on the neutral surface where the contact angle equals to 90° no matter if the surfactants are present or not. Unlike the droplet on a neutral surface, the presence of surfactants changes not only the interfacial tension between two fluids but also the contact angle for a droplet on a hydrophilic or hydrophobic surface. Thus, it is also interesting to investigate how the presence of surfactants affects the droplet behavior on a hydrophilic surface and on a hydrophobic surface. In an attempt to do so, we consider a droplet with different contact angles (i.e., θ0 = 45° and 135°) attached on a solid wall subject to a shear flow, where the clean (xin = 0) and surfactant-laden (xin = 0.1) droplets are compared at Cae = 0.1. Figure 17 shows the snapshots of the clean and surfactantladen droplets for Cae = 0.1 and θ0 = 45°. Clearly, the surfactant-laden droplet recoils and quickly reaches a steady shape, while the clean droplet keeps spreading over the surface and cannot reach a steady shape even at γ̇t = 60. According to this result, it is expected in the oil production from oil-wet reservoirs by water flooding that the addition of surfactants is able to reduce the adhesion of oil droplets to rock surface, thus enhancing the rate of oil recovery and saving energy. The reduced wetting length (or adhesion) in the surfactant-laden case can be explained by the distributions of normalized surfactant concentration, which are shown in Figure 18. It is seen that the surfactants are swept from the rear of the droplet to the front by the shear flow, which results in a higher interfacial tension than the clean droplet at the rear interface

Figure 15. The snapshots of the droplet moving on a solid wall when subjected to a shear flow for Cae = 0.35 at (a) γ̇t = 0, (b) γ̇t = 10, (c) γ̇t = 15, (d) γ̇t = 21.4, (e) γ̇t = 26.0, and (f) γ̇t = 28.2. The clean and surfactant-laden droplets are represented by the black and red lines, respectively.

and a lower interfacial tension at the front interface. Since the interfacial tension acts to resist against the droplet deformation, the higher interfacial tension enables the surfactant-laden droplet to deform less at the rear where the interface is closer to its static shape at γ̇t = 0; and on the other 7867

DOI: 10.1021/acs.langmuir.9b00495 Langmuir 2019, 35, 7858−7870

Article

Langmuir

Figure 16. Droplet shapes and the distributions of normalized surfactant concentration at γ̇t = 20 for (a) Re = 0.1, (b) Re = 1, and (c) Re = 10. The x coordinate relative to the center of the computational domain x0 is normalized by H. Note that the droplet has reached a steady deformation at γ̇t = 20 in these cases.

Table 3. Steady State Results at Various Values of Re for Cae = 0.15 and xin = 0.3 Re

Sr

ψmin *

ψmax *

0.1 1 10

0.0549 0.0665 0.1448

0.1303 0.1143 0.0835

1.9592 1.9143 1.7468

Figure 19. The snapshots of the droplet moving on a solid wall when subjected to a shear flow for Cae = 0.1 and θ0 = 135° at (a) γ̇t = 0, (b) γ̇t = 2, (c) γ̇t = 4, (d) γ̇t = 6.5, (e) γ̇t = 10, and (f) γ̇t = 14. The clean and surfactant-laden (xin = 0.1) droplets are represented by the black and red lines, respectively.

Figure 20. Distributions of normalized surfactant concentration at several typical times for Cae = 0.1, xin = 0.1 and θ0 = 135°.

than its initial value, leading to a higher local interfacial tension than in the clean case, which can prevent the deformation of the local interface (see Figure 19b−d) and thus its departure from the solid surface (see Figure 19e,f).

■

Figure 17. The snapshots of the droplet moving on a solid wall when subjected to a shear flow for Cae = 0.1 and θ0 = 45° at (a) γ̇t = 0, (b) γ̇t = 20, (c) γ̇t = 40, and (d) γ̇t = 60. The clean and surfactant-laden (xin = 0.1) droplets are represented by the black and red lines, respectively.

CONCLUSIONS In this paper, a recently developed hybrid method, which couples the LB color-gradient model and FD method for interfacial flows with insoluble surfactants, is extended to the incorporation of contact-line dynamics. The color-gradient model is implemented in the framework of MRT, which is helpful to minimize spurious velocities, and potentially enhance numerical instability and reproduce correct contact angle. A dynamic contact angle formulation that describes the dependence of local contact angle on the surfactant concentration is derived, and the resulting contact angle is enforced by a geometrical wetting condition with implementation following our recent work.34 We first simulate the static contact angles for a clean droplet and a surfactant-laden droplet, and the results show that the presence of surfactants can significantly modify the contact angle, especially when the surface is more hydrophilic or more hydrophobic. We then simulate a surfactant-laden droplet moving on a solid surface subject to a linear shear flow for varying effective capillary number (Cae), Reynolds number (Re) and surface wettability, in which the obtained results are often compared with those in

Figure 18. Distributions of normalized surfactant concentration at several typical times for Cae = 0.1, xin = 0.1 and θ0 = 45°.

hand, the surfactant-laden droplet deforms more at the front where the interface curvature is higher (see Figure 17b−d). Figure 19 plots the snapshots of the clean and surfactantladen droplets for Cae = 0.1 and θ0 = 135°. In spite of similar shapes between the clean droplet and the surfactant-laden droplet, the wetting length decreases more slowly and the droplet departs from the solid surface later in the surfactantladen case. This is because, as shown in Figure 20, the surfactant concentration in the vicinity of contact points is less 7868

DOI: 10.1021/acs.langmuir.9b00495 Langmuir 2019, 35, 7858−7870

Article

Langmuir

Wettability Alteration. J. Dispersion Sci. Technol. 2016, 37, 1259− 1267. (5) Howe, A. M.; Clarke, A.; Mitchell, J.; Staniland, J.; Hawkes, L.; Whalan, C. Visualising Surfactant Enhanced Oil Recovery. Colloids Surf., A 2015, 480, 449−461. (6) Delshad, M.; Najafabadi, N. F.; Anderson, G.; Pope, G. A.; Sepehrnoori, K. Modeling Wettability Alteration by Surfactants in Naturally Fractured Reservoirs. SPE Reservoir Eval. Eng. 2009, 12, 361−370. (7) Yang, L.; Li, S.; Liu, J.; Cheng, J. Fluid Mixing in Droplet-Based Microfluidics with T Junction and Convergent-Divergent Sinusoidal Microchannels. Electrophoresis 2018, 39, 512−520. (8) Wang, J.; Wang, J.; Feng, L.; Lin, T. Fluid Mixing in DropletBased Microfluidics with a Serpentine Microchannel. RSC Adv. 2015, 5, 104138−104144. (9) Shembekar, N.; Chaipan, C.; Utharala, R.; Merten, C. A. Droplet-Based Microfluidics in Drug Discovery, Transcriptomics and High-Throughput Molecular Genetics. Lab Chip 2016, 16, 1314− 1331. (10) Dey, M.; Vivek, A. S.; Dixit, H. N.; Richhariya, A.; Feng, J. J. A Model of Tear-Film Breakup with Continuous Mucin Concentration and Viscosity Profiles. J. Fluid Mech. 2019, 858, 352−376. (11) Lai, M.-C.; Tseng, Y. H.; Huang, H. Numerical Simulation of Moving Contact Lines with Surfactant by Immersed Boundary Method. Commun. Comput. Phys. 2010, 8, 735−757. (12) Xu, J.-J.; Ren, W. A Level-Set Method for Two-Phase Flows with Moving Contact Line and Insoluble Surfactant. J. Comput. Phys. 2014, 263, 71−90. (13) af Klinteberg, L.; Lindbo, D.; Tornberg, A.-K. An Explicit Eulerian Method for Multiphase Flow with Contact Line Dynamics and Insoluble Surfactant. Comput. Fluids 2014, 101, 50−63. (14) Zhang, Z.; Xu, S.; Ren, W. Derivation of a Continuum Model and the Energy Law for Moving Contact Lines with Ensoluble Surfactants. Phys. Fluids 2014, 26, 062103. (15) Yon, S.; Pozrikidis, C. Deformation of a Liquid Drop Adhering to a Plane Wall: Significance of the Drop Viscosity and the Effect of an Insoluble Surfactant. Phys. Fluids 1999, 11, 1297−1308. (16) Ganesan, S. Simulations of Impinging Droplets with SurfactantDependent Dynamic Contact Angle. J. Comput. Phys. 2015, 301, 178−200. (17) Karapetsas, G.; Sahu, K. C.; Matar, O. K. Effect of Contact Line Dynamics on the Thermocapillary Motion of a Droplet on an Inclined Plate. Langmuir 2013, 29, 8892−8906. (18) Karapetsas, G.; Sahu, K. C.; Matar, O. K. Evaporation of Dessile Droplets Laden with Particles and Insoluble Surfactants. Langmuir 2016, 32, 6871−6881. (19) Gunstensen, A. K.; Rothman, D. H.; Zaleski, S.; Zanetti, G. Lattice Boltzmann Model of Immiscible Fluids. Phys. Rev. A: At., Mol., Opt. Phys. 1991, 43, 4320−4327. (20) Halliday, I.; Law, R.; Care, C. M.; Hollis, A. Improved Simulation of Drop Dynamics in a Shear Flow at Low Reynolds and Capillary Number. Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys. 2006, 73, 056708. (21) Liu, H.; Valocchi, A. J.; Kang, Q. Three-Dimensional Lattice Boltzmann Model for Immiscible Two-Phase Flow Simulations. Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys. 2012, 85, 046309. (22) Swift, M. R.; Osborn, W. R.; Yeomans, J. M. Lattice Boltzmann Simulation of Non-Ideal Fluids. Phys. Rev. Lett. 1995, 75, 830−833. (23) Pooley, C. M.; Kusumaatmaja, H.; Yeomans, J. M. Contact Line Dynamics in Binary Lattice Boltzmann Simulations. Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys. 2008, 78, 056709. (24) Wang, Y.; Shu, C.; Huang, H. B.; Teo, C. J. Multiphase Lattice Boltzmann Flux Solver for Incompressible Multiphase Flows with Large Density Ratio. J. Comput. Phys. 2015, 280, 404−423. (25) Kusumaatmaja, H.; Hemingway, E. J.; Fielding, S. M. Moving Contact Line Dynamics: From Diffuse to Sharp Interfaces. J. Fluid Mech. 2016, 788, 209−227.

the clean case to show the role of surfactants. For varying Cae, the simulations are conducted by considering a neutral surface. It is found at low values of Cae that the droplet eventually reaches a steady deformation and moves at a constant velocity. In either clean or surfactant-laden case, the contact-line capillary number (Cacl) of moving droplet linearly increases with Cae, but the slope is always higher in surfactant-laden case where the droplet exhibits a bigger deformation. This suggests that the presence of surfactants can not only increase the droplet deformation but also speed up the droplet motion, which are attributed to low interfacial tension at downstream of the droplet and Marangoni-induced drag force on the droplet surface, respectively. When Cae is increased beyond a critical value (Cae,c), the droplet cannot reach a steady deformation and the breakup would occur. The presence of surfactants is found to decrease the value of Cae,c, but it shows a non-monotonic effect on the droplet breakup due to the competition between the non-uniform effects associated with non-uniform capillary pressure and the surfactant dilution by interfacial stretching. Increasing Re is able to increase the droplet deformation remarkably, thus leading to an increased surfactant dilution. The effect of surfactants on the droplet behavior is shown to heavily depend upon the surface wettability. For a hydrophilic surface, the presence of surfactants can decrease the wetting length and enables the droplet to reach the steady state earlier; whereas for a hydrophobic surface, it increases the wetting length and thus delays the departure of the droplet from the solid surface.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Haihu Liu: 0000-0002-0295-1251 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work is supported by the National Natural Science Foundation of China (nos. 51876170, 51506168 and 51711530130), and the National Key Research and Development Project of China (no. 2016YFB0200902). H.L. gratefully acknowledges the financial supports from Thousand Youth Talents Program for Distinguished Young Scholars, and the Young Talent Support Plan of Xi’an Jiaotong University.

■

REFERENCES

(1) Asadollahfardi, G.; Darban, A. K.; Noorifar, N.; Rezaee, M. Mathematical Simulation of Surfactant Flushing Process to Remediate Diesel Contaminated Sand Column. Adv. Environ. Res. 2016, 5, 213− 224. (2) Besha, A. T.; Bekele, D. N.; Naidu, R.; Chadalavada, S. Recent Advances in Surfactant-Enhanced In-Situ Chemical Oxidation for the Remediation of Non-Aqueous Phase Liquid Contaminated Soils and Aquifers. Environ. Technol. Innovat. 2018, 9, 303−322. (3) Wu, B.; Li, H.; Du, X.; Zhong, L.; Yang, B.; Du, P.; Gu, Q.; Li, F. Correlation between DNAPL Distribution Area and Dissolved Concentration in Surfactant Enhanced Aquifer Remediation Effluent: A Two-Dimensional Flow Cell Study. Chemosphere 2016, 144, 2142− 2149. (4) Hou, B.; Wang, Y.; Cao, X.; Zhang, J.; Song, X.; Ding, M.; Chen, W. Mechanisms of Enhanced Oil Recovery by Surfactant-Induced 7869

DOI: 10.1021/acs.langmuir.9b00495 Langmuir 2019, 35, 7858−7870

Article

Langmuir (26) Carmeliet, J.; Chen, L.; Kang, Q.; Derome, D. Beyond-Cassie Mode of Wetting and Local Contact Angles of Droplets on Checkboard-Patterned Surfaces. Langmuir 2017, 33, 6192−6200. (27) Zhao, H.; Ning, Z.; Kang, Q.; Chen, L.; Zhao, T. Relative Permeability of Two Immiscible Fluids Flowing through Porous Media Determined by Lattice Boltzmann Method. Int. Commun. Heat Mass Transfer 2017, 85, 53−61. (28) Chen, L.; Kang, Q.; Robinson, B. A.; He, Y.-L.; Tao, W.-Q. Pore-Scale Modeling of Multiphase Reactive Transport with Phase Transitions and Dissolution-Precipitation Processes in Closed Systems. Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys. 2013, 87, 043306. (29) Randive, P.; Dalal, A.; Sahu, K.; Biswas, G.; Mukherjee, P. P. Wettability Effects on Contact Line Dynamics of Droplet Motion in an Inclined Channel. Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys. 2015, 91, 053006. (30) He, X.; Doolen, G. D. Thermodynamic Foundations of Kinetic Theory and Lattice Boltzmann Models for Multiphase Flows. J. Stat. Phys. 2002, 107, 309−328. (31) Yiotis, A. G.; Psihogios, J.; Kainourgiakis, M. E.; Papaioannou, A.; Stubos, A. K. A Lattice Boltzmann Study of Viscous Coupling Effects in Immiscible Two-Phase Flow in Porous Media. Colloids Surf., A 2007, 300, 35−49. (32) Zhang, J.; Li, B.; Kwok, D. Y. Mean-Field Free-Energy Approach to the Lattice Boltzmann Method for Liquid-Vapor and Solid-Fluid Interfaces. Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys. 2004, 69, 032602. (33) Liu, H.; Ba, Y.; Wu, L.; Li, Z.; Xi, G.; Zhang, Y. A Hybrid Lattice Boltzmann and Finite Difference Method for Droplet Dynamics with Insoluble Surfactants. J. Fluid Mech. 2018, 837, 381−412. (34) Liu, H.; Ju, Y.; Wang, N.; Xi, G.; Zhang, Y. Lattice Boltzmann Modeling of Contact Angle and its Hysteresis in Two-Phase Flow with Large Viscosity Difference. Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys. 2015, 92, 033306. (35) Xu, M.; Liu, H. Prediction of Immiscible Two-Phase Flow Properties in a Two-Dimensional Berea Sandstone Using the PoreScale Lattice Boltzmann Simulation. Eur. Phys. J. E 2018, 41, 124. (36) Liu, H.; Zhang, Y.; Valocchi, A. J. Modeling and Simulation of Thermocapillary Flows Using Lattice Boltzmann Method. J. Comput. Phys. 2012, 231, 4433−4453. (37) Liu, H.; Zhang, Y. Modelling Thermocapillary Migration of a Microfluidic Droplet on a Solid Surface. J. Comput. Phys. 2015, 280, 37−53. (38) Ba, Y.; Liu, H.; Sun, J.; Zheng, R. Color-Gradient Lattice Boltzmann Model for Simulating Droplet Motion with Contact-Angle Hysteresis. Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys. 2013, 88, 043306. (39) Ding, H.; Spelt, P. D. M. Wetting Condition in Diffuse Interface Simulations of Contact Line Motion. Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys. 2007, 75, 046708. (40) Ba, Y.; Liu, H.; Sun, J.; Zheng, R. Three Dimensional Simulations of Droplet Formation in Symmetric and Asymmetric TJunctions Using the Color-Gradient Lattice Boltzmann Model. Int. J. Heat Mass Transfer 2015, 90, 931−947. (41) d’Humières, D.; Ginzburg, I.; Krafczyk, M.; Lallemand, P.; Luo, L.-S. Multiple-Relaxation-Time Lattice Boltzmann Models in Three Dimensions. Philos. Trans. R. Soc., A 2002, 360, 437−451. (42) Porter, M. L.; Coon, E. T.; Kang, Q.; Moulton, J. D.; Carey, J. W. Multicomponent Interparticle-Potential Lattice Boltzmann Model for Fluids with Large Viscosity Ratios. Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys. 2012, 86, 036701. (43) Pooley, C. M.; Furtado, K. Eliminating Spurious Velocities in the Free-Energy Lattice Boltzmann Method. Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys. 2008, 77, 046702. (44) Yu, Z.; Fan, L.-S. Multirelaxation-Time Interaction-PotentialBased Lattice Boltzmann Model for Two-Phase Flow. Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys. 2010, 82, 046708.

(45) Lallemand, P.; Luo, L.-S. Theory of the Lattice Boltzmann Method: Dispersion, Dissipation, Isotropy, Galilean Invariance, and Stability. Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys. 2000, 61, 6546−6562. (46) Liu, H.; Valocchi, A. J.; Zhang, Y.; Kang, Q. Lattice Boltzmann Phase-Field Modeling of Thermocapillary Flows in a Confined Microchannel. J. Comput. Phys. 2014, 256, 334−356. (47) Xu, Z.; Liu, H.; Valocchi, A. J. Lattice Boltzmann Simulation of Immiscible Two-Phase Flow with Capillary Valve Effect in Porous Media. Water Resour. Res. 2017, 53, 3770−3790. (48) Xu, J.-J.; Li, Z.; Lowengrub, J.; Zhao, H. A Level-Set Method for Interfacial Flows with Surfactant. J. Comput. Phys. 2006, 212, 590− 616. (49) Teigen, K. E.; Song, P.; Lowengrub, J.; Voigt, A. A DiffuseInterface Method for Two-Phase Flows with Soluble Surfactants. J. Comput. Phys. 2011, 230, 375−393. (50) Teigen, K. E.; Li, X.; Lowengrub, J.; Wang, F.; Voigt, A. A Diffuse-Interface Approach for Modeling Transport, Diffusion and Adsorption/Desorption of Material Euantities on a Deformable Interface. Commun. Math. Sci. 2009, 4, 1009. (51) Jiang, G.-S.; Shu, C.-W. Efficient Implementation of Weighted ENO Schemes. J. Comput. Phys. 1996, 126, 202−228. (52) Xu, J.-J.; Zhao, H.-K. An Eulerian Formulation for Solving Partial Differential Equations Along a Moving Interface. J. Sci. Comput. 2003, 19, 573−594. (53) Ladd, A. J. C. Numerical Simulations of Fluid Particulate Suspensions Via a Discretized Boltzmann Equation. Part 2. Numerical Results. J. Fluid Mech. 1994, 271, 311−339. (54) Huang, J.-J.; Huang, H.; Wang, X. Numerical Study of Drop Motion on a Surface with Stepwise Wettability Gradient and Contact Angle Hysteresis. Phys. Fluids 2014, 26, 062101. (55) Kawasaki, A.; Onishi, J.; Chen, Y.; Ohashi, H. A lattice Boltzmann Model for Contact-Line Motions. Comput. Math. Appl. 2008, 55, 1492−1502. (56) Palaparthi, R.; Papageorgiou, D. T.; Maldarelli, C. Theory and Experiments on the Stagnant Cap Regime in the Motion of Spherical Surfactant-Laden Bubbles. J. Fluid Mech. 2006, 559, 1−44. (57) Bel Fdhila, R.; Duineveld, P. C. The Effect of Surfactant on the Rise of a Spherical Bubble at High Reynolds and Peclet Numbers. Phys. Fluids 1996, 8, 310−321. (58) Griffith, R. M. The Effect of Surfactants on the Terminal Velocity of Drops and Bubbles. Chem. Eng. Sci. 1962, 17, 1057−1070. (59) Chen, J.; Stebe, K. J. Marangoni Retardation of the Terminal Velocity of a Settling Droplet: The Role of Surfactant PhysicoChemistry. J. Colloid Interface Sci. 1996, 178, 144−155. (60) Riaud, A.; Zhang, H.; Wang, X.; Wang, K.; Luo, G. Numerical Study of Surfactant Dynamics During Emulsification in a T-Junction Microchannel. Langmuir 2018, 34, 4980−4990.

7870

DOI: 10.1021/acs.langmuir.9b00495 Langmuir 2019, 35, 7858−7870