Three-dimensional aspects of fluid flows in channels. I. Meniscus and Thin Film regimes
aa r X i v : . [ phy s i c s . f l u - dyn ] S e p Three-dimensional aspects of fluid flows in channels. I. Meniscus and Thin Filmregimes
R. Ledesma-Aguilar, ∗ A. Hern´andez-Machado, and I. Pagonabarraga Departament d’Estructura i Constituents de la Mat`eria. Universitat de Barcelona,Avinguda Diagonal 647, E-08028 Barcelona, Spain Departament de F´ısica Fonamental. Universitat de Barcelona,Avinguda Diagonal 647, E-08028 Barcelona, Spain (Dated: November 15, 2018)We study the forced displacement of a fluid-fluid interface in a three-dimensional channel formedby two parallel solid plates. Using a Lattice-Boltzmann method, we study situations in which aslip velocity arises from diffusion effects near the contact line. The difference between the slip andchannel velocities determines whether the interface advances as a meniscus or a thin film of fluid isleft adhered to the plates. We find that this effect is controlled by the capillary and P´eclet numbers.We estimate the crossover from a meniscus to a thin film and find good agreement with numericalresults. The penetration regime is examined in the steady state. We find that the occupation fractionof the advancing finger relative to the channel thickness is controlled by the capillary number andthe viscosity contrast between the fluids. For high viscosity contrast, Lattice-Boltzmann resultsagree with previous results. For zero viscosity contrast, we observe remarkably narrow fingers. Theshape of the finger is found to be universal.
I. INTRODUCTION
Advancing fronts in fluid systems involve the motionof a fluid-fluid interface, a surface that lives in a three-dimensional world, and which is often constrained by asolid boundary. A typical example is that of an interfacemoving in a channel[1, 2, 3].Examples of advancing fronts in channels are imbibi-tion, a process in which a wetting fluid invades the chan-nel due to an uncompensated capillary pressure, and theviscous fingering process[2, 3, 4], where a low-viscosity(orhigh-density) fluid penetrates a high-viscosity(or low-density) one.The problem of viscous fingering in a channel has beenwidely studied in the framework of Hele-Shaw theory. AHele-Shaw cell is the two-dimensional limiting case of avery thin channel, where the equations of motion are av-eraged over the channel thickness. This reduces the inter-face to a line, the leading interface , that lives in the planeof the cell. As a consequence, the approximation discardsany effects arising from the full three-dimensional struc-ture of the interface.Nonetheless, penetration in the gap of a Hele-Shawcell is a fundamental three-dimensional effect that hasimportant repercussions in the viscous fingering prob-lem. As theoretical studies have pointed out[5], a thinfilm of viscous fluid left adhered to the cell plates asthe front advances modifies the capillary pressure at theleading interface, thus altering the front morphology.This has been confirmed in experiments of steady vis-cous fingers[6], where the presence of a thin film led tofingers not predicted by the two-dimensional theory. In ∗ Electronic address: [email protected] the following paper, we will address the role of the thinfilm in viscous fingers.A thin wetting film is not the only consequence ofa three-dimensional interfacial structure. In the con-text of liquid films spreading over dry substrates[7, 8],where a two-dimensional approximation is typically ap-plied, three-dimensional effects are also important. Forinstance, the stability of a spreading front depends onthe wetting properties of the fluid. Experimentally, ithas been observed[9] that a crossover from unstable tostable fronts occurs when the dynamic contact angle ex-ceeds π/
2, a situation which renders the velocity fieldwithin the film three-dimensional.The problem of a moving interface in a three-dimensional channel must take into account a dynamiccontact line, the intersection point between the fluid-fluid interface and the channel walls. In classic fluidmechanics, a moving contact line violates the usual no-stick boundary condition, leading to a divergent viscousdissipation[10]. Hence, contact line dynamics must con-sider a regularizing mechanism of the viscous dissipa-tion singularity. A slip velocity in the vicinity of adriven contact line arises naturally in diffuse interfacemodels of binary fluids[11], regularizing the singularity.These models consist of the usual Navier-Stokes equa-tions coupled to a convection-diffusion equation of anorder parameter[12]. Diffuse interface effects enter theforce balance equations in the shape of order parametergradients that play the role of a Young force. As a re-sult, the contact line slips over the solid boundary[11, 13].Away from the contact line, order parameter gradientsvanish and the stick boundary condition is recovered.The size of the diffusion region, l D , is then a measureof how strong slip is for a given system and is clearly animportant parameter. This size was estimated by Briantand Yeomans[14], who characterize l D for the case of aninterface subjected to shearing walls. They focused onthe dependence of l D ( L in their notation) on the modelparameters, finding a scaling relation that was verifiednumerically.Important implications arising from a relatively largeor small slip velocity compared to the leading interfacevelocity in forced fronts can be foreseen. Whenever bothvelocities are comparable, the interface should maintaina meniscus shape. Conversely, as the slip velocity be-comes small compared to the channel velocity the inter-face shape should develop as a finger, leaving a thin filmof fluid adhered to the walls of the channel.In this paper we study the penetration process acrossthe channel thickness. We study the motion of the fullthree-dimensional interface between two viscous fluidswhen it is subjected to a gravitational body force. Wetreat the case of a strictly flat leading interface, focusingonly in the details that pertain to the channel thickness.We work with symmetric fluids as well as with fluids ofdifferent densities or viscosities.We focus on two principal matters. We first describehow the contact line and leading interface velocities arerelated, and propose the mechanisms that determine thevelocity ratio. We find that the velocity ratio is controlledby the force balance at the interface and by diffusioneffects localized at the contact line.Secondly, we study the thin film that forms inevitablyin the case of small slip. In that case the front decouplesfrom the contact line leading to the growth of a finger,even when the interface is linearly stable to the Rayleigh-Taylor instability. We find that the fraction of occupa-tion of the thin film relative to the channel thicknessis a function of the capillary number and the viscositycontrast between the fluids. The high viscosity contrastcase is validated by comparing our results to the numer-ical work of Halpern and Gaver[15] which is consistentwith previous results of Taylor[16]. For fluids with zeroviscosity contrast, it turns out that the finger width hasmuch lower values than for the high viscosity contrastcase at fixed capillary number.The morphology of our fingers is very much alike to theSaffman-Taylor finger shape, a prediction of the Hele-Shaw theory. Nevertheless, it is important to stressthat although the case of a flat leading interface is two-dimensional, the equations of motion are not equivalentto those of the Saffman-Taylor problem. Therefore, pen-etration in the channel thickness cannot be attributedto the Saffman-Taylor instability. Likewise, the selectionrule of the steady state, i.e. , the actual dependence ofthe thin film thickness with the front velocity, cannotbe mapped to the theoretical predictions of the viscousfingering theory.We will address these matters by means of numericalsimulations of the mesoscopic equations of the system.To do so, we take advantage of a powerful integration al-gorithm in fluid dynamics, the Lattice-Boltzmann schemefor binary fluids.The paper is organized as follows. In Sec. II wepresent the equations that govern the system in the meso- scopic regime. In Sec. III we briefly present the Lattice-Boltzmann algorithm for binary fluid flows. Sec. IV Ais dedicated to simulation results of the forced interface,from which two steady state regimes are found; a non-penetrating regime, in which the interface advances as ameniscus, and a penetrating one, in which a single fin-ger emerges and achieves steady state. In Sec. IV B wepresent a scaling argument of the equations of motionthat leads to an estimate of the ratio between the slip andfront velocities. Such argument explains the crossoverfrom one regime to the other. In Sec. IV C we extendour results to fluids of different densities or viscosities.Sec. IV D is devoted to the steady state finger. Finally,in Sec. V we present the conclusions of this work. II. GOVERNING EQUATIONS
We consider a channel formed by two solid plates par-allel to the xy plane, each of length L and infinite width,located at positions z = 0 and z = b . Initially, two flu-ids fill the channel and are separated by a flat interfaceperpendicular to the solid walls, as shown in Fig. 1. Theequilibrium contact angle is hence θ E = π/
2. Contactlines are located at z = 0 and z = b , while the leadinginterface is located at z = b/ φ ( ~r );and order parameter which is constant in the bulk ofeach fluid and varies smoothly across a diffuse interfacialregion. Within this approximation, the equilibrium stateof the system is described by a Helmholtz free energy[12] F{ ρ, φ } = Z d ~r (cid:16) V ( φ, ρ ) + κ ~ ∇ φ ) (cid:17) . The first term in the integrand is a volume contribution,given by V ( φ, ρ ) = Aφ / Bφ / ρ/ ρ . The ρ dependent term corresponds to an ideal gas contribution,while the φ dependent terms allow for the coexistence oftwo phases. The presence of an interface is accountedfor by the last term in the integrand, which penalizesspatial variations of the order parameter by a factor κ .Minimization of F leads to the chemical potential, µ = ∂ φ V − κ ∇ φ, (cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1) Lx b fluid 1 φ xz y fluid 2 ρgξ
FIG. 1: Schematic representation of the system. The leadinginterface and contact line positions are indicated. and total pressure tensor[17],¯¯ P T = (cid:16) ρ + φ∂ φ V − V − κ (cid:16) φ ∇ φ + | ~ ∇ φ | (cid:17)(cid:17) ¯¯ δ + κ~ ∇ φ~ ∇ φ, where ¯¯ δ is the diagonal matrix. The pressure tensor hasan ideal contribution given by ¯¯ P = ρ ¯¯ δ , and an orderparameter contribution. In equilibrium, the order pa-rameter profile for the flat interface (sketched in Fig. 1)is φ ∗ ( x, z ) = − φ eq tanh( x/ξ ), where φ eq = ( − A/B ) / isthe bulk equilibrium value of the order parameter and ξ = ( − κ/ A ) / is the length scale of the interfacial re-gion; this profile leads to the difference between equilib-rium values ∆ φ = 2 φ eq and the energy per unit area ofthe interface, σ = ( − κA / B ) / . Since the interfaceis diffuse, a choice for the nominal interface position hasto be made. We choose the level surface φ = 0.The divergence of the pressure tensor yields the forceper unit volume that acts on the fluid: − ~ ∇ P − φ~ ∇ µ .The first term is the pressure gradient, while the sec-ond arises from order parameter inhomogeneities. Con-sequently, the Navier-Stokes equations are[12], ρ (cid:16) ∂ t ~v + ~v · ~ ∇ ~v (cid:17) = − ~ ∇ P − φ~ ∇ µ + η ∇ ~v + ρ~g, (1)where ~v is the fluid velocity, η is the fluid viscosity and ~g is the acceleration of gravity.The dynamics of the order parameter are described bya convection-diffusion equation, ∂ t φ + ~v · ~ ∇ φ = M ∇ µ, (2)where M is a mobility. For small deviations from theequilibrium configuration, an expansion of the chemicalpotential in powers of φ − φ ∗ yields a first order diffusioncoefficient D = M ( A + Bφ eq ), so the relative importanceof the advective and diffusive terms can be estimatedthrough a P´eclet number, P e = | ~v · ~ ∇ φ | / | D ∇ φ | .The system can be represented as a sheet of fluid inthe xz plane with periodic boundary conditions appliedin the y direction. This is equivalent to a channel ofinfinite width in the y direction with a flat leading in-terface. Stick boundary conditions are imposed at thewalls, ~v ( x, z = 0) = ~v ( x, z = b ) = ~
0, while no flowboundary conditions are imposed for the order param-eter, φ~v ( x, z = 0) = φ~v ( x, z = b ) = ~
0. At both endsof the channel the flow is homogeneous. This is ensuredby setting ∂ x ρ~v ( x = 0 , z ) = ∂ x ρ~v ( x = L, z ) = ~ ∂ x φ~v ( x = 0 , z ) = ∂ x φ~v ( x = L, z ) = ~ . Contact line dynamics arise from the diffuse natureof the interface, which allows for slip in the interfacialregion by a diffusive mechanism. The size over whichslip takes place, l D , is a function of the fluid propertiesand has been estimated by Briant and Yeomans[14], whohave given a scaling relation, l D ∼ ( ηξ M/ ∆ φ ) / . III. LATTICE BOLTZMANN METHOD
We solve numerically Eqs. (1) and (2) by means ofthe Lattice-Boltzmann algorithm presented in Ref. [17].The dynamics are introduced by discretized Boltzmannequations of two distribution functions, f i ( ~r + ~c i , t + 1) − f i ( ~r, t ) = − τ f ( f i − f eqi ) + F fi , (3)and g i ( ~r + ~c i , t + 1) − g i ( ~r, t ) = − τ g ( g i − g eqi ) . (4)In these equations, f i and g i are distribution functions,where the index i counts over the model velocity set.Space is discretized as a cubic lattice where nodes arejoined by velocity vectors, ~c i . Space and time units inEqs. (3) and (4) are set to unity. Likewise, the densityof the fluids is set to one. We use the D3Q15 velocityset, which consists of fifteen velocity vectors: six of mag-nitude 1 that correspond to nearest neighbors, eight ofmagnitude √ c s = 1 / √
3. In Eqs. (3) and (4), distribution functions arefirst relaxed to equilibrium values, represented by f eqi and g eqi , with relaxation timescales τ f and τ g . The term F fi is related to the external forcing. Following the collisionstage, distribution functions are propagated to neighbor-ing sites.Hydrodynamic variables are defined through momentsof the f i and g i . The local density and order parame-ter are given by P i f i = ρ and P i g i = φ. The fluidmomentum and order parameter current, are defined as P i f i ~c i = ρ~v and P i g i ~c i = φ~v . Local conservationof mass and momentum is enforced through the con-ditions P i f eqi = ρ , P i g eqi = φ , P i f eqi ~c i = ρ~v and P i g eqi ~c i = φ~v . In equilibrium, the pressure tensor andchemical potential are defined as P i f eqi ~c i ~c i = ρ~v~v + ¯¯ P T and P i g eqi ~c i ~c i = ˆ M µ ¯¯ δ + φ~v~v .The equilibrium distribution functions and the forcingterm are written as expansions in powers of ~v [18], i.e. , f eqi = ρω ν (cid:18) A fν + 3 ~v · ~c i + 92 ~v~v : ~c i ~c i − v + ¯¯ G f : ~c i ~c i (cid:19) ,g eqi = ρω ν (cid:18) A gν + 3 ~v · ~c i + 92 ~v~v : ~c i ~c i − v + ¯¯ G g : ~c i ~c i (cid:19) and F fi = 4 ω ν (cid:18) − τ f (cid:19) h ~f · ~c i (1 + ~v · ~c i ) − ~v · ~f i . Here, ν stands for the three possible magnitudes of the ~c i set. Coefficient values are ω = 2 / , ω = 1 / ω √ = 1 / A f = 9 / − / P, A f = A f √ = 1 /ρ Tr ¯¯ P and ¯¯ G f = 9 / (2 ρ ) ¯¯ P − δ Tr ¯¯ P ; A g = 9 / − / Mµ, A g = A g √ = 3 ˆ M µ/ρ and ¯¯ G g = 9 / (2 ρ ) ˆ Mµ (¯¯1 − ¯¯ δ ) , where ¯¯1 isthe unit matrix.Eqs. (1) and (2) can be recovered as a Chapman-Enskog expansion of Eqs. (3) and (4)[18]. The Lattice-Boltzmann scheme maps to the hydrodynamic modelthrough the relaxation timescales, i.e. , η = (2 τ f − / M = ( τ g − /
2) ˆ M, and through the body force ~f = ρ~g .Solid boundaries in the Lattice-Boltzmann method areimplemented by means of the well known bounce-backrules[18, 19]. In the lattice nodes that touch the solid,the propagation scheme is modified so the distributionfunctions are reflected to the fluid rather than absorbedby the solid. As a consequence, a stick condition forthe velocity is recovered approximately halfway from thefluid node to the solid node. IV. RESULTS
We study the process of penetration across the channelthickness in the presence of a dynamic contact line. Aswe have explained above, fingering is expected wheneverthe slip velocity is small compared to the leading inter-face velocity. In our model, slip is controlled by diffusionin the vicinity of the contact line. To measure the im-portance of diffusivity we use a typical definition of theP´eclet number,
P e = U b/D , where U is the velocity ofthe leading interface. The other relevant control param-eter is the capillary number, which follows from the ratiobetween viscous and capillary forces, Ca = ηU/σ . Wefocus on flows governed by viscous and capillary forces.To enforce this situation we neglect the convective termin Eq. (1). To assure that we work on the low Mach num-ber regime, the fluid velocity is restricted to U ≤ . A. Effect of diffusivity, surface tension andviscosity
We first consider two fluids with equal viscosities anddensities. The size of the interface is set to ξ = 0 . φ and its spatial deriva-tives across the interface[17]. Starting from a flat in-terface configuration, we perform a set of five runs atfixed forcing, viscosity and surface tension. For each runwe choose a different diffusion coefficient, which we fixthrough the mobility. In terms of dimensionless numbersthis corresponds to fix Ca and vary P e . Parameter val-
FIG. 2: Interface evolution in the channel thickness directionfor varying diffusion strength. Time interval between inter-faces is δt ≃ .
17 in b/U units. The thick profile in eachfigure corresponds to the latest time. Meniscus regime: (a) D = 0 .
073 and (b) D = 0 . D = 0 . D = 0 .
012 and (e) D = 0 . ues are U = 5 × − , η = 10 − and σ = 4 . × − ,where U is the expected leading interface velocity, calcu-lated as U = b ρg/ (8 η ). Channel dimensions are b = 23and L = 500.In Fig. 2 we show a time sequence of the interface po-sition for each run. In our simulations, v y = 0 so a flatleading interface is located at z = b/
2. Sequences (a)and (b) correspond to runs with the highest diffusioncoefficients. In both cases a steady meniscus is clearlyobserved. It is also appreciable that the meniscus in se-quence (a), corresponding to the highest diffusivity, isless curved than the meniscus in sequence (b). The nextthree sequences, (c), (d) and (e), show an abrupt changein the interface configuration. Instead of a meniscus, weobserve a penetrating structure that emerges from thecenter of the channel leaving a thin film of fluid adheredto the solid plates. The finger width in runs (c)-(e) isapproximately 17 lattice spacings. For the size of theinterface used, the order parameter saturates to its equi-librium value at the solid surface. Nonetheless, to ruleout any effects associated to the size of the interface, wehave verified that the finger width (relative to the channelthickness) does not depend on b , as we will see bellow.All runs achieve a steady state in which the velocityof the leading interface is constant. This velocity is thesame for runs (a) and (b) and due to mass conservationis slightly larger (a few percent) for runs (c), (d) and (e).The capillary number is not affected much by this effect, TABLE I: Parameter values for η and σ varying runs, U ≃ × − , D = 7 . × − . σ = 4 . × − η = 10 − η shape σ shape0.1 meniscus 0.0044 meniscus0.2 finger 0.0037 meniscus0.4 finger 0.0032 meniscus0.6 finger 0.0027 finger and we will take it as constant. The relevant effect isassociated to the variation of diffusivity.The velocity of the contact line increases with increas-ing diffusivity, as can be deduced from the contact lineposition in sequences (c), (d) and (e). Nevertheless, thevelocity of the leading interface and the width of the pen-etrating finger are the same for all three runs. This is adirect confirmation of the fact that contact line dynam-ics are decoupled from leading interface dynamics in thepresence of a thin film, as proposed by Park and Homsyin Ref.[5].It is clear from these runs that the crossover for pen-etration is set by the difference between the leading in-terface velocity, U , and the slip velocity at the contactline, v s . For a meniscus, v s = U , while penetration oc-curs whenever v s < U . As v s depends on the strength ofdiffusivity, we can draw as a conclusion that penetrationcan be achieved by increasing P e .We now explore the effect of capillarity on the dynam-ics of the interface. To do so we force the interface at fixedvelocity, diffusivity and surface tension(resp. viscosity)while we vary the viscosity (resp. surface tension). As aconsequence,
P e is fixed while Ca is varied.Results are summarized in Table I. The first columnshows parameters for runs in which the viscosity is var-ied. We observe that penetration occurs as η increases.The second column in Table I shows results for varyingsurface tension. We observe that penetration occurs as σ is decreased. We can conclude that capillarity plays asimilar role as diffusivity, as penetration occurs for lowvalues of Ca . B. The Onset of Penetration
Our results suggest that the crossover from the menis-cus regime to the thin film regime is controlled at leastby two mechanisms. On the one hand, viscous stressesdeform the interface. As a result surface tension tendsto restore the interface shape to its equilibrium value.On the other hand, advection causes order parametergradients. As a consequence, diffusivity generates a slipvelocity at the contact line. In this section we will seethat the balance between these mechanisms is controlledby
P e and Ca .Let us write the force balance per unit volume of fluidin the frame of reference of the interface. We introduce orthogonal curvilinear coordinates, s , the arclength alongthe curve φ = 0, and u , the normal distance to a pointon this curve. In terms of these coordinates the normalcomponent of Eq. (1) (in absence of inertial terms) is ∂ u P = − φ∂ u µ + η ∇ v n + ρg n , (5)where the subscript n stands for the normal componentand the subscript u denotes differentiation with respectto u .The force per unit area acting on the interfacial regionis obtained by integrating (5) across the interface:∆ P = − σ ( κ Dσ − κ Eσ ) + (cid:0) η ∇ v n + ρg n (cid:1) ξ, (6)where the term σ ( κ Dσ − κ Eσ ) arises from the integration ofthe chemical potential term[12], with κ Dσ and κ Eσ beingthe dynamic and equilibrium curvatures, which are pos-itive for a bump protruding in the x direction. We haveassumed that neither of the last two terms in the righthand side vary appreciably across the interface. Eq. (6)should be interpreted as the usual Gibbs-Thomson con-dition plus a dynamic term proportional to ξ , which van-ishes either in equilibrium or in the sharp interface limit.We will now examine Eq. (6) in the vicinity of the con-tact line. The mesoscopic nature of the interface givesrise to a finite size region where diffusion is important.This results in a slip velocity, v s , for the contact line. Wenow reproduce the scaling argument presented in Ref.[14]to obtain the diffusion size, l D , and consequently v s . Wewill subsequently compare the slip velocity to the lead-ing interface velocity in terms of Ca and P e , which areparameters that can be linked to experiments.The slip velocity and the size of the diffusion region fixthe magnitude of viscous dissipation in Eq. (6),∆ P ∼ − σ ( κ Dσ − κ Eσ ) + (cid:18) ηv s l D + ρg n (cid:19) ξ. (7)Since in the contact line region the time variation of theorder parameter is ∂φ/∂t ≃ v s ∆ φ/ξ , the order parametervariation obeys v s ∆ φξ ∼ D ∆ φl D . (8)Using Eq. (8) to eliminate l D from Eq. (7) we get∆ P ∼ − σ ( κ Dσ − κ Eσ ) + ηv s D + ρg n ξ. The last term in this expression is order ξ , while therest of terms are of order ξ . The term in the left handside is the excess pressure drop caused by the curvaturedifference, which is small in our simulations. Neglectingboth the pressure gradient and the body force we extractthe following scaling law for the slip velocity: v s ∼ σ ( κ Dσ − κ Eσ ) Dη . aCa − Thin Film Meniscus Ca − P e
FIG. 3: P´eclet and capillary numbers for meniscus and thinfilm regimes. (cid:3)
Meniscus regime. (cid:4)
Thin film regime. ◦ Meniscus regime (different viscosities). • Thin film regime(different viscosities). △ Meniscus regime (different densi-ties). N Thin film regime (different densities). The solid linecorresponds to Eq. (9), with a ≃ . . The interface curvature is a consequence of the under-lying velocity profile, which is set by the thickness ofthe channel. Therefore, the curvature difference scales as( κ Dσ − κ Eσ ) ∼ a/b , with a being a typical amplitude. Usingthis expression and measuring v s in units of the leadinginterface velocity we arrive at the following expression: (cid:16) v s U (cid:17) ∼ aCa − P e − . This indicates that both
P e and Ca control how the slipvelocity compares to the leading interface velocity. For ameniscus v s = U , so we arrive at the following condition: P e = aCa − . (9)To check the validity of the prediction we perform sev-eral runs of forced interfaces varying simulation param-eters in a wide range(see Tables III and IV). We coverup to four decades in the P e and Ca until numerical sta-bility issues of the code begin to show up. In Fig. 3 weshow a plot of our results in the Ca − P e plane. Datashown in the figure sketches two regions; at high
P e and Ca values the thin film regime is observed, whereas themeniscus corresponds to low values of both parameters.We also show our prediction, for which we find a fittingvalue for a , namely a ≃ . . Typical experiments with molecular liquids correspondto high, O (10 ) − O (10 ), values of the P´eclet number[6,20], thus falling in the thin film region sketched in thediagram. However, menisci are expected in systems witha diffuse interface, such as colloid-polymer mixtures[21].In terms of experimental parameters, the diffusion lengthscales as l D ∼ ( ηξ D/σ ( κ Dσ − κ Eσ )) / . In colloid-polymermixtures the ratio ( ξ /σ ) / is about 10 times largerthan for molecular liquids. As a consequence, in suchsystems menisci should be observable at relatively highcapillary numbers. For molecular liquids, this effect canbe achieved by decreasing the system size, for instance, in microfluidic devices, where the typical size of the chan-nel is a few microns[22]. For such small sizes, the P´ecletnumber is O (1), about 10 times smaller than for tradi-tional channels. Hence, the transition from menisci tothin films would be observable in the regime of relativelyhigh capillary number, say Ca O (10 − ) − O (10 ). C. Asymmetric Fluids
We have shown that the ratio between the leading in-terface and contact line velocities is controlled by the in-terplay between the P´eclet and capillary numbers. Weexpect that this fact holds for fluids of either differ-ent densities or viscosities. A forced interface betweenasymmetric fluids can be destabilized by virtue of theRayleigh-Taylor instability, when the more dense fluiddisplaces the less dense one. We explore situations forwhich the instability is absent. To ensure this, we keepthe channel thickness below the first unstable wave-length, given by l c = 2 π ( σ/ ∆ ρg ) / [23].The forcing is set according to ~f =8 η ( φ ) U exp /b A ( φ )ˆ x . where the local viscos-ity is set according to the mixing rule η ( φ ) =( η + η )(1 − cφ/φ eq ) /
2, characterized by the vis-cosity contrast c = ( η − η ) / ( η + η ), and U exp isthe maximum expected velocity for a Poiseuille profile.The φ dependent part is set as A ( φ ) = 1 if c = 0 and A = ( φ + φ eq ) / ∆ φ otherwise. In experiments the typicalsituation is that an effectively inviscid fluid displaces aviscous one, which corresponds to c → . We approachthis situation by setting c = 0 . . Following the generalconvention in the literature, here we define the capillarynumber as Ca = η U/σ . For all cases, we consider thatboth fluids have the same diffusion coefficient.We have performed a set of runs in which we vary
P e at fixed Ca for fluids with finite density or viscositycontrast. The details of the runs are summarized in Ta-ble II. For each case, both menisci and fingers can beobtained depending on the value of the P´eclet number.As expected, penetration occurs for sufficiently high P e .We can conclude that the appearance of a thin film isindependent of the Rayleigh-Taylor instability. In Fig. 3we plot results of this section in the
P eCa − plane. Forfluids of different densities this value is consistent withthe symmetric estimate of a ≃ .
3. For fluids of differentviscosities the crossover occurs at a ≃ .
5. Anyhow, thequalitative behavior remains independent of the degreeof asymmetry between the fluids.
D. Steady state finger in the channel thicknessdirection
We now turn our attention to the steady state fingerthat appears for high values of the product
CaP e , whichis the usual situation in most experiments. As we haveshown in Sec. IV A, diffusion only affects the motion of
TABLE II: Parameter values for runs of asymmetric fluids.Varying densities Varying viscosities σ = 9 . × − , U ≃ × − σ = 4 . × − , U ≃ × − D shape D shape0.146 meniscus 0.0488 meniscus0.098 meniscus 0.0244 meniscus0.049 finger 0.0122 meniscus0.024 finger 0.006 meniscus0.018 finger 0.001 finger the contact line, and has a negligible effect in the steadystate finger. The relevant control parameters, as pro-posed in the literature, are the capillary number and theviscosity contrast between the fluids. The steady stateif often characterized by measuring the finger width, λ b ,which is the fraction of occupation of fluid 1 relative tothe channel thickness.We explore λ b as a function of Ca at a given c value.We consider three situations, c = 0 and zero densitycontrast, c = 0 with finite density contrast and c = 0with zero density contrast. For the last case we choose c = 0 .
90 and c = 0 .
95. The low Ca runs have been per-formed varying b to rule out lattice artifacts. We havefound that results do not depend on the channel thick-ness chosen, the smallest thickness considered here being b = 23. Tabs. V, VI and VII display the parameter valuesused in each run.In Fig. 4 we plot λ b as a function of Ca . We findthat λ b depends on the viscosity contrast, as the c = 0points fall in a clearly different curve than the c = 0 . c = 0 .
95 points.On the contrary, the density contrast does not play arelevant role. Points belonging to the c = 0 curve wereobtained using five different gap sizes; b = 147, b = 23, b = 35, b = 51 and b = 95, of which the b = 35 setwas done with fluids of different densities. Results showno difference between zero or finite density difference. Infact, the gap size for the latter was large enough for the Halpern and Gaver b = 147 c = 0 . b = 95 c = 0 . b = 51 c = 0 . b = 35 c = 0 . b = 23 c = 0 . b = 147 c = 0 . b = 51 c = 0 . b = 23 c = 0 . Caλ b FIG. 4: Finger width as a function of the Ca . The error in themeasured finger width is calculated from one lattice spacingand corresponds to approximately δλ b ≃ .
04 in the figure.
Rayleigh-Taylor instability to be present. This meansthat the finger can develop as a consequence of low diffu-sion or by virtue of the Rayleigh-Taylor instability. Still,the steady state remains insensitive to the mechanismthat leads to penetration and is selected by Ca and c .Previous analytic predictions correspond to the low Ca regime at c = 1 and were carried out first byBretherton[24], who found that the finger width decaysas λ b → − . Ca / , as Ca →
0. An extensionwas done by Taylor[16], up to
Ca < .
09, for whichhe reported a decaying exponent of one-half. Numeri-cally, Reinelt and Saffman[25] solved the Stokes equa-tions using a finite difference algorithm, and consideredvalues up to
Ca <
2, which match the one-half expo-nent of Taylor at small Ca . Halpern and Gaver[15] ex-tended the prediction beyond Ca = 2 by means of aboundary element analysis of the Stokes equations. Theirresults can be fitted to an exponential law λ b = 1 − . (cid:0) − exp( − . Ca . ) (cid:1) (shown in Fig. 4), whichreproduces Reinelt and Saffman results and matches thepower law prediction of Taylor for low Ca . For large Ca ,this law saturates to a limiting value of λ b = 0 . et al [26] studied the range 0 . ≤ Ca < . ≤ Ca ≤ . Ca <
2. For
Ca > − ≤ Ca ≤ [28]. They match Halpernand Gaver prediction as c →
1. Already at c = 0 .
90, wereproduce accurately the finger width saturation value,for which we find λ b = 0 . ± . Ca , theerror increases for the c = 0 . Ca ≃ .
09 the error for c = 0 . c = 0 .
95 it reduces to 2%. For Ca = 0 . c = 0 .
95. This agreement shows thatthe Lattice-Boltzmann approach gives accurate resultsfor a wide range of Ca , improving previous results[27].As can be seen from Fig. 4, fingers with zero viscositycontrast, a case that has not been studied previously,are much narrow than fingers with c = 0 . c = 0 . Ca has a powerlaw behavior for 0 . ≤ Ca ≤
1, with an exponent m =0 . ± .
02. For
Ca O (10), the finger width saturates toa notably small value, λ b ≃ . ± . xy plane of a Hele-Shawcell as a consequence of viscosity or density asymmetriesbetween the fluids. Our problem is fundamentally dif-ferent because penetration in the channel thickness canoccur for linearly stable interfaces in the context of hy-drodynamic stability, e.g., by virtue of low diffusivity atthe contact line. Moreover, even in the case where the in-terface is linearly unstable, it is due to a Rayleigh-Taylorinstability and not through the Saffman-Taylor one.Still, we compare our finger profiles with the Saffman-Taylor ones. To do so, we recall the semi-empiric shapefound by Pitts[29], which for our geometry readscos (cid:18) πz ′ λ b b (cid:19) = exp (cid:18) πx ′ λ b b (cid:19) , (10)where x ′ and z ′ measure the distance from the finger tip.For both c = 0 and c = 0 .
9, Eq. (10) is a good approx-imation to our profiles at low Ca . This agreement islost gradually as Ca increases. In Fig. 5 we show inter-face profiles for c = 0 . c = 0 . Ca considered. Profiles for all other Ca values liebetween the shown profiles and are omitted from the fig-ure. For c = 0 .
0, Eq. (10) describes better our profiles forlarge Ca , while for c = 0 . Ca increases. V. CONCLUSIONS
We have studied the forced motion of a fluid-fluid inter-face in a three-dimensional channel by means of a meso-scopic model that takes into account contact line dynam-ics.Our results describe two possible scenarios regardinginterface dynamics. A meniscus regime is found when-ever the contact line velocity is comparable to the lead-ing interface velocity. Conversely, when the contact linevelocity is smaller than the leading interface velocity themeniscus configuration is lost, leading to penetration ofone fluid into the other in a fingering fashion. A thin πxλ b b π z λ b b FIG. 5: Rescaled interface profiles for the lowest and highest Ca values of the c = 0 . c = 0 . c = 0 . ◦ Ca = 0 .
01 and • Ca = 74 .
68. For c = 0 . (cid:3) Ca = 0 .
17 and (cid:4) : Ca = 52 .
82. Solid Line: Pitts semi-empirical finger shape.The error bar is shown in the bottom-right and correspondsto one lattice spacing. The length of the diffuse interfacecorresponds to approximately one half of a unit in the figure. film of displaced fluid is hence left adhered to the platesof the channel.The crossover from the meniscus to the thin film regimeis controlled by the competition between surface and vis-cous stresses, as well as by the competition of diffusiveand advective timescales, on top of the usual hydrody-namic instabilities. These mechanisms can be accountedfor through simple scaling arguments. We find a pre-diction for the crossover in terms of the capillary andP´eclet numbers which describes accurately our numeri-cal results. Menisci are found for low capillary and P´ecletnumbers, when surface tension and diffusion dominateover viscous stresses and advection respectively. An ex-ample of a system where diffusion is important is thatof colloid-polymer mixtures. For such systems, the rel-atively large size of the interface together with low sur-face tensions leads to large diffusion regions near de con-tact line. For instance, in Ref.[21] the size of the in-terface is typically ξ ≃ µ m while the surface tensionis σ ≃ µ N. For molecular liquids the size of the inter-face is of the order of nanometers, while σ ≃
10 mN /m .As explained in Sec.IV B, the size of the diffusion re-gion scales as l D ∼ ( ξ /σ ) / , all other parameters keptconstant. With these values the ratio of the diffusionlength between colloid-polymer mixtures and molecularliquids is at least two orders of magnitude. Thin filmsare obtained for high values of both parameters, whenadvection and viscous stresses are dominant. Our pre-diction works for both symmetric and asymmetric flu-ids. From the crossover prediction, we propose that thinfilms can be assured as long as CaP e > . CaP e > . b = 7 . × − m, witha silicone oil with η = 9 . σ = 20 . D = 1 . × − cm / s[30]. Typical velocities in theexperiments ranged from U = 0 .
01 m/s to U = 0 . CaP e ≃ . × , which isconsistent with our prediction for the thin film regime.We have examined the steady state of the thin filmregime. We have found, in agreement with Ref.[5], thatcontact line dynamics does not affect the steady statefinger shape. The capillary number and the viscositycontrast between the fluids determine the shape of thefinger.For a low-viscosity fluid pushing a high-viscosity one,the finger narrows with increasing capillary number downto a limiting value of 0 .
57 in units of the channel thick-ness. These results agree with the numerical results ofRef.[15] for the wide range of capillary numbers consid-ered. Due to computational limitations we do not inves-tigate fingers with very low capillary numbers. Nonethe-less, as we recover results from Ref.[15], we expect thatthe low capillary number limit can be recovered by ourmethod as well.For fluids with equal viscosities we have found a curveof the finger width as a function of the capillary numberthat does not follow any previous results. The width ofthe finger decreases with increasing capillary number, anexpected observation, but to a remarkably limiting widthof 0 .
38 in units of the channel thickness. This contrastswith the saturation value of the asymmetric case.The steady state is independent of whether or not theinterface is linearly unstable to the Rayleigh-Taylor insta-bility. This reinforces the conjecture of the steady statebeing independent of the mechanism that first leads topenetration of one fluid to the other. For low capillarynumbers the shape of our fingers is universal and is con-sistent with the finger shape of Pitts, which suggests thatthe steady state can be described on simple mechanicalequilibrium grounds.In a future work we will address the problem of vis-cous fingering allowing for a non-flat leading interface.As we have shown, the shape of the interface across thechannel thickness can be controlled by tuning the diffu-sion strength in the contact line. Hence, it is possible todescribe situations in which the leading interface under-goes a fingering process both in the presence and absenceof a thin film in the channel thickness direction. In thepresence of a thin film, additional control over the in-terface shape can be gained by choosing between fluidsof equal or different viscosities. These features are veryconvenient to study in detail the three-dimensional effectsthat arise in the viscous fingering problem.
VI. ACKNOWLEDGMENTS
We acknowledge financial support from Direcci´on Gen-eral de Investigaci´on (Spain) under projects FIS 2006-12253-C06-05 and FIS 2005-01299. R.L.-A. wishes to ac-knowledge support from CONACyT (M´exico) and Fun-daci´on Carolina(Spain). Part of the computational workherein was carried on in the MareNostrum Supercom-puter at Barcelona Supercomputing Center.
APPENDIX A: PARAMETER VALUES
TABLE III: Parameter values which develop a meniscus. Forall runs b = 23. σ D η U Ca P e TABLE IV: Parameter values which develop a finger. For allruns b = 23. σ D η U Ca P e σ η U Ca D P e λ b Uniform Forcing, c = 0 . b = 230.0044 0.10 0.0084 0.192 0.00185 104.2 0.6610.0037 0.10 0.0085 0.228 0.00222 87.87 0.6500.0032 0.10 0.0085 0.270 0.00265 73.95 0.6410.0027 0.20 0.0089 0.662 0.00064 321.1 0.5610.0023 0.20 0.0090 0.790 0.00076 272.34 0.5090.0019 0.20 0.0091 0.936 0.00091 229.33 0.490 b = 510.0044 0.10 0.0085 0.194 0.00442 97.68 0.6780.0037 0.10 0.0085 0.230 0.00449 96.84 0.6530.0032 0.10 0.0086 0.272 0.00459 95.59 0.6400.0027 0.20 0.0087 0.646 0.00490 90.31 0.5320.0023 0.20 0.0090 0.789 0.00489 93.75 0.5240.0019 0.20 0.0091 0.941 0.00490 95.05 0.497 b = 950.0044 0.10 0.0084 0.192 0.01608 49.64 0.7060.0037 0.10 0.0085 0.229 0.01621 49.75 0.6820.0032 0.10 0.0086 0.272 0.01634 49.85 0.6610.0027 0.10 0.0087 0.323 0.01648 49.93 0.6580.0023 0.10 0.0087 0.383 0.01661 50 0.6220.0019 0.10 0.0088 0.455 0.01675 50.04 0.598 b = 1470.0046 0.005 0.0095 0.010 0.1464 9.54 0.8350.0046 0.01 0.0094 0.020 0.1464 9.54 0.828 TABLE VI: Parameter values for stationary fingers presentedin Sec. IV D. For all runs b = 35. σ η U Ca D P e λ b Non-Uniform Forcing, c = 0 . b = 350.0046 0.50 0.0049 0.535 0.00071 241.17 0.5520.0023 0.50 0.0051 1.109 0.00071 250.09 0.4930.0011 0.50 0.0052 2.270 0.00071 255.98 0.4500.0006 0.50 0.0053 4.604 0.00071 259.57 0.4230.0003 0.50 0.0053 9.272 0.00071 261.38 0.4040.00014 0.50 0.0054 18.617 0.00071 262.39 0.3930.00007 0.50 0.0054 37.306 0.00071 262.91 0.3880.00002 0.50 0.0016 43.688 0.00037 154.85 0.3830.00004 0.50 0.0054 74.685 0.00071 263.16 0.386TABLE VII: Parameter values for stationary fingers presentedin Sec. IV D. σ η U Ca D P e λ b Uniform Forcing, c = 0 . b = 230.0046 0.38 0.0092 0.756 0.00976 21.57 0.6700.0011 0.38 0.0097 3.218 0.00244 91.76 0.6030.0023 0.38 0.0095 1.569 0.00488 44.73 0.6440.0006 0.38 0.0099 6.519 0.00122 185.91 0.5850.0003 0.38 0.0099 13.124 0.00061 374.28 0.5780.0001 0.38 0.0100 26.373 0.00031 752.09 0.5740.0001 0.38 0.0100 52.817 0.00015 1506.22 0.572 b = 510.0046 0.095 0.0040 0.082 0.0061 33.3 0.8050.0046 0.095 0.0108 0.223 0.0061 90.3 0.743Uniform Forcing, c = 0 . b = 1470.0046 0.095 0.0041 0.084 0.14 4.07 0.8220.0046 0.049 0.0007 0.008 0.14 0.71 0.9310.0046 0.488 0.0008 0.087 0.01 9.89 0.810 [1] P.-G. de Gennes, F. Brochard-Wyart, and D. Qu´er´e, Capillarity and Wetting Phenomena: Drops, Bubbles,Pearls and Waves , Springer-Verlag, New York, 2003.[2] Y. Couder,
Perspectives in Fluid Dynamics , CambridgeUniversity Press, Cambridge, 2000.[3] P. Pelc´e, editor,
Dynamics of Curved Fronts , AcademicPress, San Diego, 1988.[4] P. Saffman and G. Taylor, The penetration of a fluid intoa porous medium or Hele-Shaw cell containing a moreviscous liquid, Proc. R. Soc. London, Ser. A , 312(1958).[5] C. Park and G. Homsy, 2-Phase displacement in Hele-Shaw cells: Theory, J. Fluid Mech. , 291 (1984).[6] P. Tabeling and A. Libchaber, Film draining and theSaffman-Taylor problem, Phys. Rev. E , 794 (1986).[7] S. M. Troian, X. L. Wu, and S. A. Safran, Fingeringinstability in thin wetting films, Phys. Rev. Lett. ,1496 (1989).[8] M. P. Brenner, Instability mechanism at driven contactlines, Phys. Rev. E , 4597 (1993).[9] I. Veretennikov, A. Indeikina, and H.-C. Chuang, Frontdynamics and fingering of a driven contact line, J. FluidMech. , 81 (1998).[10] P. de Gennes, Wetting: statics and dynamics, Rev. Mod.Phys. , 827 (1985).[11] D. Jacqmin, Contact-line dynamics of a diffuse fluid in-terface, J. Fluid Mech. , 57 (2000).[12] A. Bray, Theory of phase ordering-kinetics, Adv. Phys. , 357 (1994).[13] H.-Y. Chen, D. Jasnow, and J. Vi˜nals, Interface and con-tact line motion in a two phase fluid under shear flow,Phys. Rev. Lett. , 1686 (2000).[14] A. Briant and J. Yeomans, Lattice Boltzmann simula-tions of contact line motion. II. Binary fluids, Phys. Rev.E , 031603 (2004).[15] D. Halpern and D. Gaver, Boundary element analysis ofthe time-dependent motion of a semi-infinite bubble in achannel, J. Comp. Phys. , 366 (1994).[16] G. Taylor, Deposition of a viscous fluid on the wall of atube, J. Fluid Mech. , 161 (1961).[17] V. M. Kendon, M. E. Cates, I. Pagonabarraga, J. C. De-splat, and P. Bladon, Inertial effects in three-dimensionalspinodal decomposition of a symmetric binary mixture: a Lattice-Boltzmann study, J. Fluid Mech. , 147(2001).[18] A. Ladd and R. Verberg, Lattice-Boltzmann simulationsof particle-fluid suspensions, J. Stat. Phys. , 1191(2001).[19] J. C. Desplat, I. Pagonabarraga, and P. Bladon, Lud-wig: A parallel Lattice-Boltzmann code for complex flu-ids, Comp. Phys. Comm. , 273 (2001).[20] P. Petitjeans and T. Maxworthy, Miscible displacementsin capillary tubes, J. Fluid Mech. , 37 (1996).[21] D. Aarts, M. Schmidt, and H. Lekkerkerker, Direct visualobservation of thermal capillary waves, Science , 847(2004).[22] P. Tabeling, Introduction to Microfluidics , Oxford Uni-versity Press, Oxford, 2005.[23] S. Chandrasekhar,