Three-dimensional aspects of fluid flows in channels. II. Effects of Meniscus and Thin Film regimes on Viscous Fingers
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. II. Effects of Meniscus and ThinFilm regimes on Viscous Fingers
R. Ledesma-Aguilar, ∗ I. Pagonabarraga, and A. Hern´andez-Machado 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: October 31, 2018)We perform a three-dimensional study of steady state viscous fingers that develop in linear chan-nels. By means of a three-dimensional Lattice-Boltzmann scheme that mimics the full macroscopicequations of motion of the fluid momentum and order parameter, we study the effect of the thick-ness of the channel in two cases. First, for total displacement of the fluids in the channel thicknessdirection, we find that the steady state finger is effectively two-dimensional and that previous two-dimensional results can be recovered by taking into account the effect of a curved meniscus acrossthe channel thickness as a contribution to surface stresses. Secondly, when a thin film develops inthe channel thickness direction, the finger narrows with increasing channel aspect ratio in agreementwith experimental results. The effect of the thin film renders the problem three-dimensional andresults deviate from the two-dimensional prediction.
I. INTRODUCTION
Interfacial instabilities in three-dimensional channelsgive rise to a rich phenomenology in systems that rangefrom nano and microscales[1] to macrometric channels[2,3, 4], and from which a number of practical applicationscan be drawn.For instance, controlled drop breakup in micro-channels has proved useful in the fabrication of low poly-dispersity micro-emulsions[5] and in the enhancementof micro-reaction processes[6, 7]. In the latter, three-dimensional effects are crucial, as they are responsibleof a vortex flow structure within the droplet[8] that en-hances the mixing process of the reactants.A widely studied interfacial instability in channels isthat of fingering, which occurs whenever a low-viscosity(or high-density) fluid drives a high-viscosity(or low-density) one. The instability, first studied by Saffmanand Taylor[9], leads to interface dynamics where finger-like structures emerge and compete. The problem has asteady-state solution, composed by a single finger of con-stant velocity U and occupies a fraction λ of the widthof the channel.Experimentally, finger growth has been studied mainlyin Hele-Shaw cells. These consist of a pair of plates oflength L and width W separated by a thickness b . Forsuch systems, it has been pointed out[10] that the station-ary finger is determined by a single control parameter, amodified capillary number defined as 1 /B = 12 Ca/ǫ .For a fluid with viscosity η and surface tension σ , thecapillary number, Ca = ηU/σ , measures the competi-tion between driving forces, such as viscous stresses andgravity, and restoring forces, like surface tension. 1 /B ∗ Electronic address: [email protected] also includes the degree of asymmetry of the cell, givenby the aspect ratio ǫ = b/W . If 1 /B is the only con-trol parameter of the system, all experimental data, i.e. all finger widths, should be described by a single curvewhen plotted as a function of this parameter. Contraryto this view, experiments show that there exists a familyof curves λ vs. /B for different aspect ratios[11, 12].This fact suggests that a three-dimensional effect, givenby the interplay between the dynamics in the channel-thickness and in the channel-width, is determinant forthe steady-state solution.Theoretically, fluid-flow in a channel at small veloci-ties pertains to the lubrication regime, in which the flowoccurs mainly along the direction of L given that it ismuch larger than both W and b . Hele-Shaw flows area limiting case in lubrication theory, where b is muchsmaller than W . Owing to the smallness of b , the prob-lem is rendered effectively two-dimensional by averagingall fields over the thickness of the channel. Averaging theequations of motion also reduces the interface from a sur-face to a line, often called the leading interface . In viewsof the averaged model, three-dimensional effects enter asperturbative corrections to the boundary conditions thathold at the leading interface in terms of Ca and ǫ , par-ticularly to the Gibbs-Thomson condition, which relatesthe pressure drop across the interface to the interfacecurvature and surface tension.Progress towards a three-dimensional description ofthe problem has been made since the pioneering workof Saffman and Taylor[9], who solved the problem of astationary finger in the absence of surface tension in twodimensions. McLean and Saffman[13] included the ef-fect of surface tension and were the first to obtain a λ vs. /B prediction by solving numerically the two-dimensional model. According to their results, λ is amonotonically decreasing function of 1 /B that saturatesto λ → / /B → ∞ . The prediction of McLean andSaffman is unique in 1 /B , so the role of the aspect ratiois precluded from their theory.The relevance of three-dimensional effects was first sug-gested by Park and Homsy[14], who pointed out that athin film of fluid in the channel-thickness direction wouldcontribute to the pressure drop at the leading interface.Using perturbation methods for slightly curved leadinginterfaces(small ǫ ), they found that for low Ca the pres-sure drop varies as Ca / , a result that matched the earlyprediction of Bretherton[15] for capillary tubes.Sarkar and Jasnow[16] used the modified pressure dropto solve the steady state finger. Their results agreed bet-ter with experiments but were restricted to low valuesof Ca . It was shown by Tabeling and Libchaber[11]that corrections to the pressure drop can be used toreduce three-dimensional experimental data to the two-dimensional results of McLean and Saffman. A modifiedpressure drop can be accounted for as an effective sur-face tension. Using the correction of Park and Homsy,Tabeling and Libchaber were able to reduce their datato McLean and Saffman results for moderately low val-ues of 1 /B , where fitting parameters were used to esti-mate the correction terms. In an experimental study[12],Tabeling, Zocchi and Libchaber observed that, contraryto the McLean-Saffman prediction, the finger width cango below the one-half limit for sufficiently high 1 /B andsufficiently large ǫ .Reinelt extended the expansion of the pressure dropup to O (1) in Ca and included the effect of larger as-pect ratios[17]. Computation of the steady state fingeryielded solutions that better agreed with experiments for O (1) values of Ca . For small ǫ Reinelt observed a betteragreement between numerics and experiments. However,for relatively large ǫ this agreement is lost.Higher Ca values have only been explored in the caseof flat leading interfaces by Halpern and Gaver[18]. Theirnumerical results are consistent with results found byReinelt and Saffman[19] for Ca O (1) and ǫ = 0, andshow that the pressure drop is insensitive to the capil-lary number for Ca > et al [20, 21] used a set of coupledevolution equations for the velocity potential and orderparameter that describes accurately the early stages ofdestabilization of the leading interface, and approachesMcLean and Saffman results as the viscosity of the dis-placing fluid is made negligible. The strict one-sided sit-uation, were one of the fluids is inviscid, was studied byHern´andez-Machado et al [22]. They used a single evo-lution equation for the concentration that includes dy-namic effects in the form of chemical potential gradients and described correctly the steady state finger.In a preceding paper[23], we have shown that a detailedthree-dimensional description of fluid-flow in a channelcan be done by means of a mesoscopic model whichwe implement numerically via a Lattice-Boltzmann al-gorithm. The model considers a fluid-fluid interface incontact with solid boundaries. In contrast to classic ap-proaches, it allows for slip at the contact line by meansof a diffusive mechanism inherent to the mesoscopic na-ture of the interface. This circumvents the complica-tions of contact line dynamics in the classic formulation,associated to the viscous dissipation singularity[24]. InRef.[23], we focused on the case of a flat leading inter-face. We showed that depending on the velocity of thecontact lines, which we control by modifying the diffusionstrength, the interface can either advance as a meniscusor develop as a finger. In the latter case, a thin film ofdisplaced fluid is left adhered to the walls of the channel.In this paper we extend our Lattice-Boltzmann simu-lations to the case of a non-flat leading interface, wherea viscous finger is expected to appear. Our aim is toprovide a detailed description of the mechanisms thataffect the steady finger and that cause deviations fromtwo-dimensional results. To do so, we study fingers thatform in the meniscus and thin film regimes separately.We cover values of Ca up to O (10) and explore variousaspect ratios.The paper is organized in the following manner. InSec. II we present the governing equations of the systemwhich we solve numerically by means of a Lattice Boltz-mann algorithm, presented in the preceding paper[23].Results are presented in Sec. III. In Sec. III A we de-scribe the simulation strategy and parameter steeringprocedure. As a validation test, in Sec. III B we com-pute the dispersion relation of the interface in the two-dimensional limit and compare it to the analytic predic-tion of the Saffman-Taylor problem. Sec. III C is devotedto the study of stationary viscous fingers; in Sec III C 1 wefocus on fingers pertaining to the meniscus regime in thechannel thickness, which we found to be effectively two-dimensional, while in Sec. III C 2 fingers in the thin filmregime are studied. We find that fingers in the thin-filmregime are three-dimensional and cannot be describedby the two-dimensional theory in general. A discussionof our results where we compare with previous experi-ments is presented in Sec.IV. In Sec. V we present theconclusions of this work. II. GOVERNING EQUATIONS
We consider the motion of two viscous fluids, whosedynamics are governed by the Navier-Stokes equations, ρ (cid:16) ∂ t ~v + ~v ( ~ ∇ · ~v ) (cid:17) = − ~ ∇ P − φ~ ∇ µ + η ∇ ~v + ρ~g. (1)Here ~v is the fluid velocity, P is the pressure, ρ is thedensity, η is the fluid viscosity and ~g is the acceleration W fluid 2 ρg yxz b fluid 1 L FIG. 1: Schematic representation of the system. Dashed linesindicate projections of the fluid-fluid interface in the xy and xz planes. The leading interface corresponds to the xy pro-jection. of gravity. The extra term, φ~ ∇ µ , is mesoscopic and ac-counts for interfacial forces. Here, φ ( ~r, t ) is an order pa-rameter and µ ( φ ) is the chemical potential. φ has theproperty of being uniform in the volume of each phaseand non-uniform in an interfacial region of typical size ξ .In the present case, volume values are chosen as φ = ± φ = 0. The size of the interfaceis set to ξ = 0 . φ obey a convection-diffusion equa-tion, ∂ t φ + ~v · ~ ∇ φ = M ∇ µ, (2)where M is a mobility coefficient. In equilibrium, thepressure and chemical potential minimize a free energyfunctional, from which explicit expressions P ( ρ ) and µ ( φ )can be derived. For further details the reader is referredto Ref.[23].We work on a linear channel, composed by two solidplates of width W and length L parallel to the xy plane,separated by a distance b , as depicted in Fig. 1. Thereexist two principal directions in the system: a lateral di-rection, parallel to the xy plane, and a transverse oneparallel to the xz plane. We will denote these by sub-scripts k and ⊥ respectively.The impenetrability and stick boundary conditionsat the walls are enforced by setting ~v ( x, y, z = 0) = ~v ( x, y, z = b ) = 0 and φ~v ( x, y, z = 0) = φ~v ( x, y, z = b ) = 0. At both ends of the channel in the x directionthe flow is homogeneous. Hence, ∂ x ρ~v ( x = 0 , y, z ) = ∂ x ρ~v ( x = L, y, z ) = 0, and ∂ x φ~v ( x = 0 , y, z ) = ∂ x φ~v ( x = L, y, z ) = 0. Periodic boundary conditions are imposedin the y direction.As for the fluid-fluid boundary, the Gibbs-Thomsonrelation is recovered by integrating Eq.(1) across the in-terfacial region, ∆ P = σ (cid:18) R k + 1 R ⊥ (cid:19) , (3)where σ is the surface tension and R α is the radius ofcurvature of the interface in the direction α .We now briefly review the classic treatment of theproblem. First, v ⊥ is assumed to be much smaller than v k , which in turn is expected to be parabolic in z . As a result, Eq.(1) is recast in the form of an average veloc-ity field which holds in the volume of each fluid, calledDarcy’s Law, h v k i = − b η (cid:16) ~ ∇ P − ρ~g (cid:17) k , (4)where triangular brackets denote an average over thechannel thickness. Under these conditions, R ⊥ is ex-pected to be much larger than R k . Hence, in the two-dimensional theory the Gibbs-Thomson relation is sim-plified to ∆ P = σ/R k .Corrections to this expression arise whenever 1 /R ⊥ is not negligible. For such cases, Libchaber andTabeling[11] have proposed that thin film effects can beaccounted for by defining an effective surface tension σ ∗ = σ (cid:18) R k R ⊥ (cid:19) , so the two-dimensional form of the Gibbs-Thomson con-dition is recovered. For this purpose, they used the esti-mation of Park and Homsy[14] of the pressure drop for Ca → ǫ → P = σ (cid:18) π R k + 3 . b/ Ca / (cid:19) . (5)As a result, their experimental results collapsed to theMcLean-Saffman curve when using the correspondingdefinition of the control parameter, 1 /B ∗ = ( σ/σ ∗ )1 /B .We solve numerically Eqs. (1) and (2) by means of aLattice-Boltzmann algorithm. For further details of themethod, the reader is referred to the preceding paper[23]. III. RESULTSA. Simulation Parameters and Setup
The traditional description of the viscous fingeringproblem corresponds to situations in which the relevantforces at play are viscous stresses and capillarity. Forthe particular case of fingering in a Hele-Shaw cell theseforces are expressed in terms of a modified capillarynumber[10] 1 /B = 12( W/b ) (∆ ηU +∆ ρgb / /σ, where∆ η and ∆ ρ are the differences in viscosity and densitybetween the fluids.To ensure that capillarity and viscous forces dominatethe dynamics of the fluids, inertia must be small com-pared to both of these forces. We enforce this situationby neglecting the convective term in Eq. (1). As for com-pressibility, we consider low Mach number flows, whichwe achieve by keeping U ≪ c s . For our scheme, it sufficesto set U ≤ . /B . Due to computation resource limi-tations, ǫ is restricted to O (10) at most for the majorityof runs. To achieve high values of 1 /B , say O (10 ), Ca LBLBHele-Shaw k ′ ω ′ FIG. 2: Dispersion relation for the linear stability of theinterface. Simulation parameters (in simulation units) are σ = 0 . η = 0 . b = 11 . ρg = 3 . × − and ( (cid:4) ) ∆ ρg = 6 . × − must then be O (10). Our strategy is to keep the interfacevelocity and the viscosity in ranges of U = O (10 − ) and η = O (10 − ). Hence, Ca can be tuned by means of thesurface tension.The channel is implemented as follows. We set arectangular simulation box of dimensions N x × N y × N z.
Due to the flow symmetry, we simulate only onefourth of the real channel by setting boundary condi-tions as follows: ∂ y ρv y ( x, y = 0 , z ) = ∂ y ρv y ( x, y = W/ , z ) = 0, ∂ y φv y ( x, y = 0 , z ) = ∂ y φv y ( x, y = W/ , z )0, ∂ z ρv z ( x, y, z = 0) = ∂ z φv z ( x, y, z = 0) = 0 . B. Linear Stability in the two-dimensional limit
We first verify the linear stability of the interface, i.e. , the behavior of an initially flat interface that hasbeen subjected to a small perturbation. We study fluidsof equal viscosities, so the instability is gravitationallydriven. This is done by fixing the body force term inEq. (1) as ρg = 8 η/b U exp ( φ + 1) /
2, where U exp is themaximum expected velocity for a Poiseuille flow. Themodified capillary number reduces to 1 /B = W ∆ ρg/σ ,In this case, the linear stability analysis of the interfaceevolution of the averaged equations yields the dispersionrelation ω ( k ) = b η k (∆ ρg − σk ) , (6)where ω is the exponential growth rate of a sinu-soidal perturbation to the flat interface solution. Theperturbation is characterized by its wavelength, Λ =2 π/k . By considering dimensionless frequencies ω ′ = ω/ ( U/ W ) B / and wavenumbers, k ′ = W B / k , thedispersion relation becomes universal, i.e. , ω ′ ( k ′ ) = k ′ (1 − k ′ ) . We prepare a base flow corresponding to a flat inter-face in the xy plane that propagates at constant velocity.The interface in the xz plane is nearly flat throughout the simulation, so the system is effectively two-dimensional.Once the base flow is fully developed, the interface isshifted according to a single mode perturbation of wave-length Λ = W and an initial small amplitude. We followthe evolution of the amplitude, A ( t ), which is measuredas A ( t ) = x tip ( t ) − ¯ x ( t ) , where x tip and ¯ x ( t ) are the inter-face tip and mean interface positions respectively. Thegrowth rate, ω , is extracted as a linear fit of log( A ( t )) vs t . Fig. 2 shows a comparison between the universal dis-persion relation and our results. To quantify the degreeof accuracy of these results, we fit our data to the gen-eral form ak ′ ( b − ck ′ ). We find a most unstable modeat k ′ max ≃ .
56 and a first unstable mode at k ′ ≃ . C. Viscous Fingers
In a preceding study[23], we have shown that it is pos-sible to control the generation of a thin film in the channelby adjusting the diffusivity of the order parameter. Al-though for usual experimental conditions this is not a rel-evant parameter (it might be relevant in nano-channels),it gives the possibility of elucidating the role of a thinfilm in viscous fingers. Diffusivity is accounted for by aP´eclet number,
P e = U b/D , where D is the diffusion co-efficient. By combining the effects of P e and Ca , one caneither suppress or induce the formation of a thin film. Inparticular, a small value of the product CaP e results inthe suppression of thin films, while the contrary is ob-tained for high
CaP e . Results from the preceding workgive a penetration threshold of
CaP e ≃ − .The strategy is to study first fingers for which CaP e ≤ − and then extend this simulations to CaP e ≫ − .
1. Meniscus Regime
We first study fingers for which no film of displacedfluid develops in the xz plane of the channel. We carryout simulations with modified capillary numbers in the TABLE I: Control parameters and finger width for runs in ofthe meniscus regime.Run ǫ Ca CaP e /B λ (a) 0.17 0.11 0.08 99 0.709(b) 0.17 0.22 0.16 198 0.675(c) 0.17 0.45 0.19 290 0.640(d) 0.06 0.11 0.04 522 0.558(e) 0.06 0.19 0.02 1045 0.525(f) 0.06 0.23 0.03 1672 0.523(g) 0.06 0.48 0.11 2090 0.529(h) 0.06 0.68 0.21 3136 0.518(i) 0.04 0.74 0.26 4175 0.521(j) 0.05 0.76 0.27 6012 0.519 FIG. 3: Interface snapshots at two different times for ǫ = 0 . /B = 99 and CaP e = 0 .
08. The thickline parallel to the xy plane corresponds to the leading interface, while the thick line parallel to the xz plane corresponds tothe interface projection in the center of the channel. Thin lines correspond to the contact lines. The first snapshot correspondsto t = 0 . b/U , while the second snapshot, at t = 17 . b/U , corresponds to the steady state finger.FIG. 4: Collapsed interface profiles in the xy plane for themeniscus regime. Parameter values corresponding to eachsymbol can be consulted in Table I. The error (small bar atthe right bottom) corresponds to one lattice spacing. Thelarger bar indicates the size of the diffuse interface, approxi-mately 3 ξ . range 100 ≤ /B ≤ ǫ = 0 .
17 to ǫ = 0 .
04. The aspectratio is decreased by decreasing the channel thickness.We summarize the simulation parameters in Table I.For each run we observe the usual phenomenology forthe leading interface. During the early stages of interfaceevolution, the amplitude of the perturbation grows untila finger emerges and widens. This stage is followed by arelaxation of the interface shape, until a Saffman-Taylorfinger develops. The finger propagates with a steady ve-locity U , leaving behind a growing region where the fingerhas flat sides. In this region a constant finger width, λW ,can be defined. As for the channel thickness, we observethat the initially flat interface rapidly relaxes to a menis-cus, which also has a steady shape. In Fig. (3) we showa three-dimensional plot of the interface for run (a) inTable I at two different times. In the plot we show boththe contact lines and the leading interface; both contact lines follow the leading interface.To check for consistency in the steady state solutionwe use the semiempirical interface profile obtained byPitts[25], which reproduces experimental results accu-rately for a wide range of finger widths. The equationfor the interface shape reads,cos( πy ′ / λ ) = exp( πx ′ / λ ) . (7)where x ′ and y ′ measure the distance from the finger tipin units of half the channel width. The natural scalingsin this equation are πx ′ / λ and πy ′ / λ . Consequently,all profiles should collapse into a single curve if thesescalings are used. Fig. 4 shows plot of interface profilescorresponding to runs of Table I. As expected, all inter-face profiles fall in the same universal curve within error.In addition, our collapse is in fair agreement with Eq.(7).The selection rule in the viscous fingering problemis expressed as the functional dependence of the fingerwidth with the modified control parameter. We compareour results with the λ vs. /B prediction of McLeanand Saffman. We find that runs with ǫ = 0 .
17 showwider fingers than predicted, while runs with smaller ǫ agree better with the two-dimensional result. Even inthe absence of a thin film, the xz interface projectionhas a certain curvature. This can be accounted for bydefining an effective surface tension in terms of the radiiof curvature of the interface, which we are able to mea-sure directly. The effective surface tension then reads σ ∗ = σ (1 + R k /R ⊥ ). The correction factor in this ex-pression is given by the quantity in parentheses, whichincreases for strongly curved meniscus. The rescaled con-trol parameter then reads 1 /B ∗ = (1 /B ) / (1 + R k /R ⊥ ) . Of course this correction should be more evident in thelow 1 /B region, where λ varies rapidly with the modi-fied control parameter. In Fig. 5 we show a plot of λ vs. /B ∗ . Points fall on the McLean-Saffman curve for thewide range of 1 /B ∗ considered, regardless of the aspectratio.
2D Prediction ǫ = 0 . ǫ = 0 . ǫ = 0 . ǫ = 0 . /B ∗ λ FIG. 5: Finger width as a function of the rescaled controlparameter 1 /B ∗ in the meniscus regime.FIG. 6: Interface projections in the xy and xz planes for runswith different CaP e values. Plots correspond to the samesimulation time. (a)
CaP e = 0 .
85, (b)
CaP e = 4 .
2. Thin Film Regime
We now extend our simulations to fingers in the filmregime. Penetration in the xz plane occurs for high CaP e , so we choose to sample 1 /B at fixed D . Conse- TABLE II: Control parameters and finger width for thin filmregime runs.Run ǫ Ca CaP e /B λ (a) 0.25 2.80 12.32 835 0.592(b) 0.25 3.36 17.74 1002 0.589(c) 0.25 6.61 68.61 2004 0.569(d) 0.25 15.9 400.41 4003 0.558(e) 0.35 8.96 1515 1403 0.549(f) 0.49 34.7 4330 5247 0.527(g) 0.64 50.9 3973 5247 0.517(h) 0.78 68.5 7192 5247 0.508(i) 1.00 91.9 8019 5247 0.493(j) 1.00 131 156598 5430 0.494 quently, CaP e increases with increasing 1 /B . To resolvethe thin film correctly we must take into account thefinite size of the interface. As explained in Ref.[23], thethin film is insensitive to the channel thickness already for b = 23. We therefore choose sufficiently thick channels.We explore a wide range of aspect ratios, 0 . ≤ ǫ ≤ . ≤ /B ≤ CaP e ∼ O (1) region, close tothe penetration threshold. In Fig. 6 we show interfaceprojections in the xy and xz planes located at z = b/ y = W/ CaP e values;(a)
CaP e = 0 .
85 and (b)
CaP e = 4 .
44. In Fig. 6(a) the in-terface in the xz plane presents a penetrating structure,but a well developed film is absent. The finger in the xy plane is not well developed either, and it presents ananomalous tip. Conversely, in Fig. 6(b) both interfaceprojections describe well developed fingers. It is thenclear that deviations from the Saffman-Taylor finger inthe xy plane are correlated to the interface structure inthe xz plane. An interesting feature of the run corre-sponding to Fig. 6(a) is that the the xz interface structureis persistent. This means that the length of the finger inthe xz plane is constant in time, a consequence of the slipvelocity of the contact line. The diffusion strength is notlarge enough to maintain a meniscus, which in the onehand makes the slip velocity smaller than the channel ve-locity. Nevertheless, as the interface relaxes to a thin filmshape, curvature deviations from equilibrium increase theslip velocity, making the contact line advance to restorethe meniscus shape.We next explore the range CaP e ≥ O (10) for whichsimulation parameters and observed finger widths aresummarized in Table II. In Fig. 7 we present snapshots ofthe three dimensional interface at two different times forrun (b) in Table II. The first snapshot corresponds to theearly stage of the finger formation. Looking at the inter-face projections in the xy plane, we see that the contactline(light line) is close to the leading interface(dark line)and no film is present in the xz plane. In the next snap-shot the contact line has moved away from the tip, thusgiving rise to the growth of a wetting film. The shapeof the finger is in agreement with the typical morphol-ogy found in experiments. To illustrate this, in Fig. 8 wecompare the shape of the finger to Eq. (7). Within error,our profiles are consistent with Pitts shape.Fig. 9 shows the measured finger width as a functionof 1 /B . The lowest aspect ratio we have considered cor-responds to ǫ = 0 . /B values considered the finger width fallsabove the McLean-Saffman prediction. We increase theaspect ratio to ǫ = 0 . ǫ islarger confirm this tendency in a systematic way. Tests(f)-(j) in the same table correspond to a fixed value of1 /B with increasing ǫ . We find that for sufficiently large ǫ the finger width goes below the one-half theoretical limitof McLean and Saffman. FIG. 7: Interface snapshots at two different times for ǫ = 0 . /B = 1002 and CaP e = 17 .
74. Thicklines correspond to the xy and xz interface projections in the center of the channel. Thin lines correspond to the contact lines.Times are t = 0 . b/U and t = 28 . b/U .FIG. 8: Rescaled interface profiles for the thin film regime.Symbols correspond to data presented in Table II. The barsin the bottom at the right indicate the error bar and diffuseinterface size as in Fig. 4. IV. DISCUSSION
Our results show that the finger width decreases withincreasing aspect ratio. To maintain 1 /B fixed while
2D Prediction ǫ = 1 . ǫ = 1 . ǫ = 0 . ǫ = 0 . ǫ = 0 . ǫ = 0 . ǫ = 0 . /Bλ FIG. 9: Finger width as a function of 1 /B .
2D Prediction ǫ = 1 . ǫ = 1 . ǫ = 0 . ǫ = 0 . ǫ = 0 . ǫ = 0 . ǫ = 0 . /B ∗ λ FIG. 10: Finger width as a function of the rescaled controlparameter for the thin film regime. varying the aspect ratio of the channel, one has to vary Ca accordingly. As a consequence, the film thickness andthe capillary pressure are altered. If we increase the as-pect ratio(as in the high-1 /B region in Fig. 9), then Ca must decrease to keep 1 /B fixed. As a consequence, thefilm thickness and the capillary pressure decrease as ǫ increases, which is consistent with a narrower finger.This behavior has been observed, for instance, in ex-periments by Tabeling, Zocchi and Libchaber[12], andaddressed in numerical calculations by Reinelt[17] wherethe effect of the thin film was introduced perturbativelyin the two-dimensional model. Experiments suggest that,for high 1 /B , increasing the cell aspect ratio has a thin-ning effect on the finger, which is what we observe inour simulations. Results of Reinelt suggest the oppositetendency.The aforementioned experiments were done at small ǫ and Ca , and at high viscosity contrast, defined as c = ( η − η ) / ( η + η ). As we have shown in Ref.[23],the thin film thickens as c →
1. Under these conditions,experiments show that the finger width goes below theone-half limit even for cells with ǫ = 0 . Ca and high c values used in experiments.In our simulations the thin film is about t/b ≃ . t/b ≃ . ǫ . To achieve the experimen-tal regime thinner film should be considered. We haveconsidered a cell aspect ratio of ǫ ≃ .
05 and c = 0 . Ca is still too large, the film is then thickenough to keep us in the low 1 /B ∗ regime, where thefinger width is still larger than one half of the channelwidth. Due to computational limitations we do not ex-plore smaller ǫ .The fact that for a given 1 /B there exist different fin-ger widths for different aspect ratios raises the doubt of1 /B as being the only control parameter present in thesystem. To this end, we compute the rescaled surfacetension σ ∗ = σ (cid:0) R k /R ⊥ (cid:1) , where the radii of curva-ture are measured at the finger tip. We then rescale thecontrol parameter according to 1 /B ∗ = ( σ/σ ∗ )1 /B . InFig. 10 we show a plot of the finger width as a functionof the rescaled control parameter. At low 1 /B ∗ , resultsagree with McLean-Saffman results. We conclude that inthis region the finger is effectively two dimensional andthat three-dimensional effects can be accounted for evenat Ca ∼ /B ∗ . Hence, we concludethat deviations from two-dimensionality are caused bythe thin film.An important feature in the λ vs. /B ∗ plot is thatfinger width appears to be determined by 1 /B ∗ uniquely.This suggests that 1 /B ∗ is the only control parameter ofthe problem.We have explored a region of values of the aspect ra-tio between the Saffman-Taylor( ǫ →
0) and Rayleigh-Taylor( ǫ = 1) limits of the fingering instability. In bothlimits, the relevant control parameter appears to be aneffective modified capillary number. In addition, the in-terface shape is remarkably universal, as suggested byFigs.4 and 8. V. CONCLUSIONS
We have carried out a detailed study of the viscousfingering problem in three-dimensional channels for fluidsof different densities and viscosities. We have studied the single finger solution for systems in which either a thinfilm develops across the channel thickness or a meniscusis stabilized.For systems in which no thin film is present, McLean-Saffman two-dimensional results describe the dependencyof the finger width as a function of a rescaled modifiedcapillary number, 1 /B ∗ , which has a correction that de-pends curvature of the interface direction of the channelthickness. This holds for arbitrary high values of 1 /B ∗ ,evidencing that a complete displacement across the chan-nel thickness renders the problem two-dimensional.We have extended our studies to situations where athin film develops across the channel. We find differentvalues of the finger width when changing the channel as-pect ratios at fixed modified capillary number, an obser-vation that is consistent with previous experiments[12].This non-uniqueness seems to disappear as the controlparameter is corrected by curvature effects associated tothe thin film, i.e. , when the finger width is compared to1 /B ∗ .For low 1 /B ∗ , the finger width collapses to theMcLean-Saffman curve. However, at high 1 /B ∗ the fin-ger width deviates from this curve, and goes bellow thelimit of one half in units of the channel width.Our work is done at high values of the capillary num-ber. Consequently, the effective capillary pressure inour simulations is large enough to keep the finger abovethe one-half limit of the two-dimensional theory for highvalues of the channel aspect ratio. Experiments inRefs.[11, 12] were done in cells with ǫ ≃ .
03 and at Ca ≃ − , a regime that falls beyond the scope of thiswork for computational reasons. Nonetheless, for low1 /B ∗ , we recover the same results, indicating that thesame mechanisms hold, even if the actual aspect ratioand capillary number are not the same in experimentsand simulations.To our knowledge, experiments of fingering in highaspect ratio channels have not been conducted so far.Our results could be confirmed, for instance, in micro-channels, where the aspect ratio is typically large andin which the meniscus to thin film transition could beobserved. This is then an open question for experimen-talists to confirm. 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. [1] P. Tabeling,
Introduction to Microfluidics , Oxford Uni-versity Press, Oxford, 2005. [2] 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.[3] Y. Couder,
Perspectives in Fluid Dynamics , CambridgeUniversity Press, Cambridge, 2000.[4] P. Pelc´e, editor,
Dynamics of Curved Fronts , AcademicPress, San Diego, 1988.[5] D. Link, S. Anna, D. Weitz, and H. Stone, Geometricallymediated breakup of drops in microfluidic devices, Phys.Rev. E , 054503 (2004).[6] K. Hosokawa, T. Fujii, and I. Endo, Handling of pico-liter liquid samples in a poly(dimethylsiloxilane)-basedmicrofluidic device, Anal. Chem. , 4781 (1999).[7] H. Song, J. Tice, and F. Ismagilov, A microfluidic systemfor controlling reaction networks in time, Angew. Chem.Int. Ed. , 767 (2003).[8] H. Kinoshita, S. Kaneda, T. Fujii, and M. Oshima,Three-dimensional measurement and visualization of in-ternal flow of a moving droplet using confocal micro-PIV,Lab Chip , 338 (2007).[9] 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).[10] G. Tryggvason and H. Aref, Numerical experiments onHele-Shaw flow with a sharp interface, J. Fluid. Mech. , 1 (1983).[11] P. Tabeling and A. Libchaber, Film draining and theSaffman-Taylor problem, Phys. Rev. E , 794 (1986).[12] P. Tabeling, G. Zocchi, and A. Libchaber, An experi-mental study of the Saffman-Taylor instability, J. FluidMech. , 67 (1987).[13] J. McLean and P. Saffman, The effect of surface tensionon the shape of fingers in a Hele Shaw cell, J. Fluid Mech. , 455 (1981).[14] C. Park and G. Homsy, 2-Phase displacement in Hele-Shaw cells: Theory, J. Fluid Mech. , 291 (1984). [15] F. Bretherton, The motion of long bubbles in tubes, J.Fluid Mech. , 166 (1961).[16] S. Sarkar and D. Jasnow, Quantitative test of solvabilitytheory for the Saffman-Taylor problem, Phys. Rev. A , 4900 (1987).[17] D. Reinelt, The effect of thin film variations and trans-verse curvature on the shape of fingers in a Hele-Shawcell, Phys. Fluids , 2617 (1987).[18] D. Halpern and D. Gaver, Boundary element analysis ofthe time-dependent motion of a semi-infinite bubble in achannel, J. Comp. Phys. , 366 (1994).[19] D. A. Reinelt and P. G. Saffman, The penetration of afinger into a viscous fluid in a channel and tube, SIAMJ. Sci. Stat. Comp. , 542 (1985).[20] R. Folch, J. Casademunt, and A. Hern´andez-Machado,Phase-field model for Hele-Shaw flows with arbitrary vis-cosity contrast. I.Theoretical approach, Phys. Rev. E ,1724 (1999).[21] R. Folch, J. Casademunt, and A. Hern´andez-Machado,Phase-field model for Hele-Shaw flows with arbitrary vis-cosity contrast. II.Numerical study, Phys. Rev. E ,1734 (1999).[22] A. Hern´andez-Machado, A. Lacasta, E. Mayoral, andE. Corvera Poir´e, Phase-field model of Hele-Shaw flowsin the high-viscosity contrast regime, Phys. Rev. E ,046310 (2003).[23] R. Ledesma-Aguilar, A. Hern´andez-Machado, and I. Pag-onabarraga, Three dimensional aspects of fluid flows inchannels: I. Meniscus and thin film regimes, Submittedto Phys. Fluids.[24] P. de Gennes, Wetting: statics and dynamics, Rev. Mod.Phys. , 827 (1985).[25] E. Pitts, Penetration of a fluid into a Hele-Shaw cell:the Saffman-Taylor experiment, J. Fluid Mech.97