The formation of drops through viscous instability

zyxwvutsrqponm zyxwvutsrq zyxwvutsrq zyx 351 J. Fluid Mech. (1995), vol. 289, p p . 351-378 Copyright @ 1995 Cambridge University Press The formation of drops through viscous instability By S I L V A N A S. S. C A R D O S O AND A N D R E W W. W O O D S Institute of Theoretical Geophysics, Department of Applied Mathematics and Theoretical Physics, University of Cambridge, Silver Street, Cambridge CB3 9EW, UK (Received 8 June 1994 and in revised form 5 December 1994) The stability of an immiscible layer of fluid bounded by two other fluids of different viscosities and migrating through a porous medium is analysed, both theoretically and experimentally. Linear stability analyses for both one-dimensional and radial flows are presented, with particular emphasis upon the behaviour when one of the interfaces is highly stable and the other is unstable. For one-dimensional motion, it is found that owing to the unstable interface, the intermediate layer of fluid eventually breaks up into drops. However, in the case of radial flow, both surface tension and the continuous thinning of the intermediate layer as it moves outward may stabilize the system. We investigate both of these stabilization mechanisms and quantify their effects in the relevant parameter space. When the outer interface is strongly unstable, there is a window of instability for an intermediate range of radial positions of the annulus. In this region, as the basic state evolves to larger radii, the linear stability theory predicts a cascade to higher wavenumbers. If the growth of the instability is sufficient that nonlinear effects become important, the annulus will break up into a number of drops corresponding to the dominant linear mode at the time of rupture. In the laboratory, a Hele-Shaw cell was used to study these processes. New experiments show a cascade to higher-order modes and confirm quantitatively the prediction of drop formation. We also show experimentally that the radially spreading system is stabilized by surface tension at small radii and by the continual thinning of the annulus at large radii. zyxw zyxwvuts 1. Introduction When a fluid is horizontally displaced by another in a porous medium, the interface between them may be either stable or unstable, depending on the relative viscosities of the fluids and on their miscibility. The basic mechanism of this instability was first described by Hill (1952) and later by Saffman & Taylor (1958). Consider the rectilinear displacement of a fluid of viscosity p2 by another fluid of viscosity p1 in a homogeneous porous medium (figure 1). We assume the motion of the fluids through the porous medium is governed by Darcy’s law K u = --Vp = -MV P P (1.1) where u denotes velocity, K is the permeability of the medium and p is pressure. M = ~ / isp the mobility of the fluid. Suppose the interface between the two fluids is deformed slightly such that there is a perturbation of the interface of thickness 6x 352 zy zyxwvutsr zyxwvuts zyx S . S . S . Cardoso and A. W Woods P2 t X FIGURE 1. Rectilinear displacement of a fluid by another fluid of different viscosity. zyxwvutsrq FIGURE 2. Dispersion relation for the rectilinear displacement of two immiscible fluids. (figure lb). Then, from (l.l), the anomalous pressure across the displaced fluid is 6 p = (1/M2 - l/M])UGx, where U is the steady state velocity. If the net pressure difference is positive, then any small perturbation to the interface will grow, leading to an instability. Hence, the interface will be unstable when a less viscous fluid displaces a more viscous fluid (MI > M2). ‘Fingers’ of the displacing fluid will then develop and penetrate into the more viscous fluid ahead. We have assumed above that the interface between the two fluids is sharp, i.e. that the fluids are immiscible. A detailed stability analysis should therefore take into account the effect of surface tension. For disturbances in the form of normal modes proportional to exp(at+iky), the growth rate of the instability is found to be (Chouke, Meurs & Poel 1959) Here T denotes the surface tension coefficient between the two fluids and k is the wavenumber of the instability. The dispersion relation (1.2) shows that the effect of surface tension (second term) is to stabilize short-wavelength perturbations. The competition between the viscous instability and the surface tension leads to a most unstable mode (figure 2). In many situations, there is a continuous gradient of viscosity in the intermediate region between the two fluids, resulting from a certain degree of miscibility. Although in this case there are no interfaces, since molecular diffusion and dispersion will lead to a smooth variation of the viscosity, instability can still occur (Paterson 1985; zyxwvuts zyx zyxwvut The formation of drops through viscous instability 353 Manickam & Homsy 1993). The mathematical analysis of these miscible displacement problems is often complicated because of the time-dependent nature of the base state. Usually various approximations are made which lead to simpler problems. The phenomenon of viscous fingering in porous media is important in many natural and industrial situations. An example of the former arises in the geological context. Sequential batches of magmas with different chemical compositions, and hence different viscosities, often rise from depth along fissures in the Earth's crust, leading to magmatic intrusions. The extent of mixing occurring during such ascent has a key role in controlling the composition of the magma, which is an issue of central importance in petrology (McBirney 1984). In the industrial context, viscous fingering occurs, for example, in oil recovery processes, fixed bed regeneration in chemical processes and underground storage of gas (Homsy 1987). In many of these applications the instability is not desired and there have been a number of attempts to suppress it. These usually involve the reduction of the unfavourable viscosity gradient by the addition of a third fluid, of variable viscosity, at the boundary between the host and displacing fluids. Mungan (1971) showed experimentally that a spatially varying viscosity in this intermediate layer can stabilize the interface and these observations have been confirmed by a theoretical analysis of Gore11 & Homsy (1983). However, the stability of the immiscible displacement of an intermediate layer of fluid, bounded by two other fluids of different viscosities, has not been analysed in detail. The purpose of the present work is to investigate the stability of such flow, both theoretically and experimentally. A linear theoretical model is developed and quantitatively compared with new laboratory experiments. This simplified model allows a deeper understanding of some of the processes present in the more complex variable-viscosity problems. We analyse the interaction between the two interfaces and the effect of the thickness of the intermediate layer, in both rectilinear and circular geometries. The latter case is particularly interesting since surface tension can stabilize a single interface at small radii (Paterson 1981). However, here we identify a new mechanism of stabilization of a circularly spreading annulus which operates at large radii. Finally, we examine situations in which the intermediate layer is unstable and may break up to form drops. In 0 2, the rectilinear flow of an immiscible layer of fluid is analysed. In 6 3, we review and extend previous work on the linear stability of the radial displacement of a single interface. Since the base state is evolving with time, the most unstable mode also changes as the interface moves outward. New experimental results are presented which confirm this theoretical approach. Combining the results of &j2 and 3, we then consider the linear instability of a radially spreading annulus of fluid. A theoretical model is developed and a quantitative comparison with experimental observations is presented. In 0 5, we consider the stabilization of the flow by suitable variation of the displacement velocity. Finally, in 0 6, we draw some conclusions. zyx 2. Rectilinear displacement We begin by analysing the stability of the rectilinear displacement of a layer of fluid in order to gain fundamental insight on the behaviour of a three-fluid system. Later in the paper, we shall build upon these results when we analyse the more complex problem of a spreading annulus. Consider the displacement in a porous medium of an intermediate layer of fluid, bounded by two other fluids, as sketched in figure 3. The three fluids are assumed to have different viscosities and to be immiscible. The motion of the fluids is described by Darcy's law (cf. (Ll)), which may be written in 354 zy zyxwvuts zyxwv zyxwv zyxwvutsr zyxwv S. S. S. Cardoso and A. FK Woods 1(3 1(2 't X D FIGURE 3. Rectilinear displacement of an intermediate layer of fluid. terms of a velocity potential 4 = Mp as u= -v+ (2.1) and mass conservation v*u=o. (2.2) The velocity potential therefore satisfies Laplace's equation v24 = 0. (2.3) Across each interface the velocity is continuous and the pressure changes by an amount corresponding to the effect of surface tension. For a planar interface, these equations admit the steady state solution 4; = -Ux + cj, j = 1,2,3, (2.4) where the subscripts 1, 2 and 3 refer to the trailing, intermediate and leading fluids respectively (figure 3). The constants c j may be determined by specifying the magnitude of the pressure at some point in the flow. In order to examine the stability of the interfaces, we consider a wavelike perturbation a = Aexp(iky a t ) at the interface between fluids 1 and 2, and a similar perturbation b = Bexp(iky a t ) at the interface between fluids 2 and 3, where (r is the growth rate of the instability and k its wavenumber. The solution of (2.3) has the form 4, = 4; 4;, with + + + zyxwvut (2.5a) (2.5b) (2%) where a, p, y and E are constants to be determined. The continuity of velocity and pressure across the interface between fluids 1 and 2 (at x = U t ) are expressed, respectively, by - u21 -- at ul1 and ($) 4l (2.6) a+--Tu,,. 4; (2.7) M1 X M2 Note that the last term in (2.7) arises from the surface tension at the interface. We adopt here the simple formulation of Chouke et al. (1959) for this effect. Later in the paper, we shall discuss the inclusion of a more precise model for the surface tension. Similar velocity and pressure conditions may be written for the interface between fluids 2 and 3 (at x = Ut D). The resulting system of equations may be X + zyxwv zyx zyxwv zyxwvuts zyxw 355 The formation of drops through viscous instability solved for the growth rate of the instability, a, and for the ratio of the amplitudes of the perturbations at the two interfaces, A / B . The solution is zyxwv a0 ($)2 + U l ( 5 ) + u2 = 0 where ao, a1 and a2 are given by a0 = (1 + M1 M2 %) M2 sinh(kD) + a1 = ((2 - 2)+ (2+ 2),,(kD)’) TM2 a2 (- = (1 - 2)+ %(kD)’) ((1 sinh(kD), - 2)+ ,(kD)’) TM2 UD sinh(kD) and alUk A -- sinh(kD) + sinh(kD) + cosh(kD) (2.9) We have assumed here that the surface tension coefficient, T , is the same at both interfaces. Depending upon the relative magnitudes of the mobilities of the fluids, each interface may be either locally stable or locally unstable. For example, if M1 > M2 then the interface between fluids 1 and 2 will be locally unstable. The global stability is determined by the relative viscosity of the bounding fluids. Thus, in the absence of the intermediate layer, if M I < M3, the displacement would be stable. Figure 4 shows a plot of the dispersion relation for a case in which both interfaces are locally unstable, M1 > M2 and M2 > M3, (and hence, the flow is also globally unstable). The configuration of the intermediate layer, obtained from the ratio of amplitudes A / B , has been sketched beside each curve. It may be seen that there are two different modes, namely a sinuous mode in which the instabilities at the two interfaces grow in phase (solid curve), and a varicose mode in which the instabilities grow in antiphase (dashed curve). Both modes are unstable to perturbations of long wavelengths, with the most unstable mode being that for which the interfaces are in phase; short wavelengths are stabilized by surface tension. In figure 5, the dispersion relation for a case in which only the leading interface is unstable is represented; the flow is globally stable since M1 < M3. The varicose mode is now stable for all wavelengths. The sinuous mode is unstable for long wavelengths, but stabilized by surface tension for short wavelengths. Similar results are obtained for a globally unstable flow with only one unstable interface: only one of the modes is unstable. The limit of a thin intermediate layer We now consider the limit of a thin intermediate layer, as this will enable a simplification of equations (2.8) and (2.9) and a fuller understanding of the dynamics of the different modes. 356 zyxwvut zyxw zyxwvutsrq zyxwvut zyxwvuts S . S . S . Cardoso and A. u/: Woods 4 2 (TD - 0 U -2 zyxwvuts zy \ -4 5 10 zyxwvu 15 20 kD FIGURE 4.Dispersion relation when the two interfaces are locally unstable, Ml/M2 MZ/M3 = 2.5. =2 and When kD << 1, we obtain (2.10a) (2.10b) and ($) (;)2 = 1- 2 T M1M3 (-UD2 Mi M3 (1-$)2+ _ M 2- 1 .~M3 =1-- M2 M1 (2.1la) - ($4)2 (2-2)(1-2)(2 - 1) $(kD)2) . (2.11b) zyxw zyxwvuts zyxwv zy zyxwvutsr The formation of drops through viscous instability zyxwvutsrq '1 41 2 -4 357 1 zyxwvutsrqponmlkjihg ! I I 0 5 10 i 15 kD FIGURE 5. Dispersion relation when only the leading interface is unstable, M2/M3 = 2 ; the flow is globally stable, Mi/M3 = 0.5. FIGURE 6. Breakup of the intermediate layer into drops. The eventual separation of the drops will occur owing to the effect of surface tension, which tends to smooth the pointing edges at the rear at the moment of rupture. Equation ( 2 . 1 0 ~describes ) the stability of a mode with a behaviour similar to the Saffman-Taylor mode (compare with (1.2)). It is determined only by the properties of the bounding fluids, 1 and 3, and when unstable, is the dominant mode (figure 4). The other solution, given by (2.10b), has a smaller growth rate and represents an internal mode. The distortion of the interfaces corresponding to the Saffman-Taylor mode is always sinuous (in phase, A / B > 0). In contrast, the internal mode leads to a varicose deformation (antiphase) of the intermediate layer when both interfaces are locally unstable ( M I > M2 > M3), whereas a sinuous deformation is observed when one of the interfaces is highly stable (figure 5). The results in figure 5 suggest that even when the flow is globally stable, so that the Saffman-Taylor mode (GI)is stable, the internal mode remains unstable because of the local instability at one of the interfaces. This result suggests why the intermediate layer may eventually break up into drops. A sketch of the expected evolution of this layer is shown in figure 6. The linear displacement problem, described above, is of fundamental interest in that it illustrates in a simple way the interaction between the local and global stabilities and 358 zyxwvuts zyxw zyxwv S. S. S. Cardoso and A . W Woods FIGURE 7. Radial displacement of a fluid by another of different viscosity. their contribution to the resulting flow pattern. However, in some practical situations the fluid may spread radially from a relatively small localized area (Paterson 1981). In the following sections of this paper, we generalize the analysis to consider such radial displacements. 3. The dynamics of the radial displacement of a single interface 3.1. Theoretical description Before extending the analysis of 5 2, by considering the radial spreading of an annulus of fluid, we first examine the behaviour of a single interface in radial source flow. This will identify some of the differences between the rectilinear and radial displacement problems, which arise because the base state evolves with time in the latter case. We will also test a simple method of determining how this time dependence influences the wavelength of the instability by comparison with experimental observations. The linear stability of a circular interface spreading in a Hele-Shaw cell was modelled by Wilson (1975), and thereafter by Paterson (1981). Figure 7 depicts the immiscible displacement of fluid 2 by fluid 1. The flow is described by Laplace’s equation (cf. (2.3)), which in planar polar coordinates reads zyxwvuts zyxw The radial velocity component is u = -& For a point source with volume flow rate per unit depth Q, the velocity potential of the steady flow satisfying (3.1) is Q 1nr + c j , j = 1,2. 271 This solution satisfies the continuity of velocity at each interface. The constants cj may be determined by prescribing both the pressure discontinuity which results from the surface tension at the interface and the magnitude of the pressure at some point in the flow. Assuming a quasi-steady flow, in which the velocity at a certain radius is locally constant, and examining the stability of the interface to wavelike perturbations of the form A,(t) exp(inO), Paterson (1981) found that the growth rate of these perturbations is given by +? = -- ’ where A , denotes the amplitude of the sinusoidal mode with discrete azimuthal zyxwvuts zyxwv zy The formation of drops through viscous instability 25 359 20 15 10 5 10 5 0 zyxw zyxw zyxwvu zy w 15 20 25 FIGURE8. Growth rate of each mode as a function of the radial position of the interface. wavenumber n and R is the radial position of the interface. If the mobility of the inner fluid is much larger than that of the outer fluid, M I >> M2, the interface is very unstable and (3.3) may be simplified to the form A,, n-1 R2 An n(n +R . (3.4) The relation above shows that for sufficiently small radius, all modes are stabilized by surface tension. However, as the interface progresses outward, the different modes become unstable at different radial positions. We may define the critical radius at which instability first occurs, &, as the radius where mode two begins to grow (since mode one corresponds only to a displacement of the circle of inner fluid relative to the point source, with no change of shape). From (3.4), & is given by 12nM2T &= (3.5) Q and (3.4) may be rewritten as In figure 8, we show the variation of the growth rate of each mode as a function of the non-dimensional radial position 9 = R/&. The minimum radius for instability, for each mode n, is R = n(n 1)&/6. As the radial position of the interface increases, the growth rate of each mode increases to a maximum value + + at R = n(n 1)&/4. The growth rate of this mode subsequently decays, although remaining positive. Therefore, as R increases, the most unstable mode n gradually increases. Analysis of second-order terms shows that the linear theory remains valid 360 zyxwv zyxwvutsr zyxwvut zyxwvu zyxw S . S. S . Cardoso and A. u! Woods lo4 10-2 1 1 I 1 10 9 , zy zyxwv 8 1 100 1000 FIGURE 9. Amplitude of each mode as a function of the radial position of the interface, for an instantaneous disturbance of amplitude A0 at the critical radius (model A). while the amplitude of the perturbation remains small compared to its wavelength, nA,/R << 1. As a result, if the initial disturbance is sufficiently small, there may be a cascade through several modes before the amplitude of the disturbance has grown sufficiently for nonlinear effects to dominate. In this case, the development of the linear instability may be described by integrating the growth rate of each mode as a function of time (or interface radius). If we assume that owing to random noise in the system, the nth mode is perturbed with an amplitude A0 at the radial position 90,= &,/&, then since Rt = Q/2nR, we have for each mode from (3.6), The evolution of the system thus depends critically upon the initial disturbance. We may envisage two different model situations. Firstly, one may consider a situation in which all modes are perturbed instantaneously with the same amplitude A0 at a particular radius l70 - model A. Figure 9 shows the predicted magnitude for each of the modes as a function of 9 when & is taken to coincide with the critical radius & (equation (3.5)), that is 90"= 1 for all n. The higher modes are initially stable and hence A, initially decays for n 2 3; these modes then become unstable at larger R and A , begins to grow. As the radius increases, each mode in turn becomes dominant for a range of radii, and hence there is a cascade to higher modes. The precise details of this cascade depend upon the magnitude of &. The cascade predicted by linear theory only continues while nonlinear effects are negligible. As an alternative initial condition, one may consider a constant low level of noise in the system throughout the experiment - model B. Each mode n is assumed to be perturbed with an amplitude A0 until the radius where it first becomes unstable, gon= n(n 1)/6; subsquently, A , begins to grow. The magnitude of each mode as a function of the radius, as predicted for this initial condition, is graphed in figure 10. + z The formation of drops through viscous instability 104 103 zyx z 361 zyxwvutsrqp 102 10' 100 1 10 9 zyx zyx zyxwv 100 FIGURE 10. Amplitude of each mode as a function of the radial position of the the interface, for a constant disturbance of amplitude A0 (model B). Comparison of figures 9 and 10 shows that model B leads to a faster cascade to higher modes as 9 increases. Note that, as A , increases, nonlinear effects become increasingly important and the linear theory developed here ceases to apply. However, we may compare our linear theory, in particular the prediction of a cascade to higher modes, with experimental observations. We first describe the experimental method. 3.2. Experimental method Two-dimensional flow in a porous medium may be modelled in the laboratory using a Hele-Shaw cell. This consists of two closely spaced parallel plates. The analogue mobility of the fluid in this system is related to this gap thickness, b, by M = b2/12p (Saffman & Taylor 1958). The Hele-Shaw cell used here was made of two 13 mm thick perspex disks, 420 mm in diameter. The disks were screwed together and separated by washers to maintain a uniform gap spacing. Most experiments were conducted with a plate spacing of 0.75 mm. The fluids were injected through a nozzle 10 mm in diameter, at the centre of the top disk. The rate of addition of the fluid was controlled with a high-precision needle valve. The pressure difference across this valve was set to approximately 2 atm by applying air pressure ; this was sufficiently large that the increasing pressure difference across the Hele-Shaw cell during an experiment was negligble. Therefore, the flow rate was virtually constant. The experiment involved the injection of air into the cell previously filled with glycerine. The flow rate of air was measured with a soap film meter (Isenberg 1978). The range of flow rates (per unit depth) covered was 0.73 to 3.2 cm2 s-l. The kinematic viscosity of the liquid was measured with a U-tube, direct flow, capillary viscometer (see table 1). 362 zyxwvuts zyxwvu zyxwv zyxwv S. S. S. Cardoso and A. W Woods Fluid glycerine diluted glycerine cooking oil natrosol solution silicone oil DClOOO silicone oil DC12500 Viscosity, p (g cm-' s-l) 7.81 1.27 0.65 93.8 11.0 128 Surface tension, T (g sp2) 63 (with air)' 16 (with silicone oil DC1000) 32 (with air)' zyxw 16 (with diluted glycerine) TABLE1. Properties of the fluids. * Tennent (1971) 3.3. Comparison with experiments Paterson (1981) observed the rapid growth of eight fingers, as air was blown at a relatively large flow rate, Q = 9.3 cm2 s-l, into glycerine in a circular Hele-Shaw cell. That author suggested that the most unstable mode at onset of the instability grew very rapidly and therefore determined the number of fingers. Here, we present new experimental results where much smaller flow rates were used. As a result, the growth rate of the instability is sufficiently small so that the linear theory remains valid over a considerable range of radial positions. This enables the observation of a cascade to higher modes. In figure 11, a sequence of photographs of an experiment with a very low rate of injection of air is shown. The glycerine has been dyed for visualization. The mobilities of air and glycerine were, respectively, M1=2.6 cm3 s g-' and M2 = 6.0 x cm3 s g-l. These values satisfy our earlier theoretical assumption that M1 >> M2. Here Q = 0.73 cm2 s-', so that R, = 0.20 cm (see table 1 and equation (3.5)). The range of 98 in this experiment is approximately 1 through to 50. The transition from mode four to mode five is clearly visible. It is difficult to say if mode three dominates at an earlier stage. These observations are in excellent accord with the theoretical predictions of model B, shown in figure 10. It seems that the perturbation cannot be seen for A,/& smaller than approximately 10. Model A (figure 9) clearly underpredicts the observed wavenumbers. In figure 12, a sequence of photographs of an experiment with a larger flow rate of air is presented. The critical radius is now R, = 0.045 cm (Q = 3.2 cm2 s-l). The range of in the photographs is approximately 1 to 133. Mode six appears to dominate from the moment the perturbation becomes visible until the end of the experiment. This is again in good agreement with the predictions of model B (figure 10). However, the rate of growth of the fingers is now quite large, and so at the later stages (photographs ( d ) and ( e ) )the linear theory is no longer valid. The results above support model B, that is, the initial disturbance was of constant magnitude resulting, for example, from irregularities on the perspex surface ; this is in contrast to model A, of an instantaneous disturbance. In all our experiments, the capillary number Ca = pRJT was, approximately, in the range 10-3-10-1. These relatively small values of Ca allow us to use the simple formulation in (2.7) and (4.4) for the pressure jump resulting from the surface tension at the interface between the two fluids (Saffman & Taylor 1958; Chouke et al. 1959). Indeed, as shown by Maxworthy (1989), more detailed theories which include both the curvature of the interface in the thin gap of a Hele-Shaw cell (Park & Homsy 1984) and the effect of the viscous film left behind on the plate as the interface moves - wetting effect (Schwartz 1986; Reinelt 1988) lead to differences in the wavenumber zyxwv zyxwv zyxwvu zyx The formation of drops through viscous instability 363 zyxwvuts zyxw zyx zyxwvuts zyxwv zyxwv FIGURE 11. Displacement of glycerine by air at a flow rate Q = 0.73 cm2 s-l: ( a ) 9 = 8.5, ( b ) 9 = 15, (c) W = 21, ( d ) W = 23 and (e) W = 44. The scale is 1 cm. and the growth rate of the most unstable mode of the order of the precision of the current experimental techniques. Also, for the range of C a in this study, the correction of the radial displacement velocity for the plate wetting effect is small, less than 10% (Saffman 1986; Maxworthy 1989), and we have therefore neglected it. The good agreement between our theoretical predictions and experimental observations supports this simple approach. 4. The dynamics of the radial displacement of an annulus of fluid 4.1. Stability analysis We now combine the methods and results of § 5 2 and 3 to consider the stability of a spreading annulus of fluid, as shown in figure 13. As described previously in 0 3, the flow is governed by Laplace’s equation (cf. (3.1)). The steady state velocity potential is $?=--I Q j = 1,2,3 (4.1) 2n nr + cj, zyxwvutsr zyxw 364 S . S. S . Cardoso and A. u/: Woods zyxw zyx zyxwvut zyxw zyxwvut FIGURE 12. Displacement of glycerine by air at a flow rate Q = 3.2 cm2 s-l ( a ) W = 38, ( b ) W = 56, (c) 9 = 78, ( d ) W = 89 and ( e ) 9 = 133. where subscripts 1, 2 and 3 refer to the inner, annulus and outer fluids, respectively. The constants c j may be determined by prescribing both the pressure difference due to surface tension at each interface and the magnitude of the pressure at some point in the flow. Suppose that at some time the positions of the inner and outer interfaces are, respectively, r = R1 and r = R2 (figure 13) and that both interfaces are perturbed from the steady state circular configuration. We assume that the interface at r = R1 undergoes a wavelike perturbation a = A,(t)exp(inO) and the interface at r = R2 a perturbation b = B,(t)exp(inO). The amplitudes A , and B, are functions of time, t, and n = 1,2,3.. . is the quantized azimuthal wavenumber. The solution of (3.1) corresponding to this linear perturbation has the form 6 j = 47 6; with zyxwv + (4.2a) zyxwv zyxwvut zyx The formation of drops through viscous instability 365 zyxwvuts FIGURE 13. Radial displacement of an annulus of fluid. (4.2b) zyxwvuts zyx (4.2~) The functions of time in the equations above may be determined by the boundary conditions at each interface. The continuity of radial velocity at the interface between fluids 1 and 2 (at r = R1) requires that + vyra v: = v;ra + ui = a,. (4.3) The pressure condition at this interface is expressed by Similar continuity conditions may be written for the interface between fluids 2 and 3 (at r = R2). The resulting system of equations may be solved for the growth rate of the instabilities of the two interfaces, A,, and Brit. As found previously for the rectilinear displacement problem (see 4 2), there is an internal mode associated with the local stability of each of the two interfaces, and a global mode determined by the properties of only fluids 1 and 3. The general solution for this situation is complex. Generally, if the global mode is unstable, then fingers of the less viscous, trailing fluid will penetrate into the fluid ahead, as in the single interface problem. If the global mode is stable, but the internal mode is unstable, then the annulus may break up into a series of drops. It is this novel latter situation which we investigate in more detail. We consider the limit in which the viscosity of the displacing fluid (fluid 1) is very large, so that its mobility is very small, M I + 0, and hence A,, A,, -+ 0. Physically, this means that the inner interface is nearly rigid and so it remains almost circular as it moves outward with time. The growth rate of the instability of mode n at the 366 zyxwvut zyxwv zyxwvutsr S. S. S. Cardoso and A. W Woods outer interface is then given approximately by zyxwvut If we assume that M2 >> M3 and also that the annulus of fluid 2 is thin compared to the wavelength of the instability (n6 << l), we have M? B,. where 6 = (R2- R 1 ) / R 1is a non-dimensional thickness of the annulus. The first term on the right-hand side of the expression above represents the opposing effects of the viscosity difference between fluids 2 and 3 (destabilizing effect) and of the surface tension at this interface (stabilizing effect). The numerator of this term is similar to that found for the single interface in radial source flow (cf. (3.4)). The last term in (4.6) represents a new stabilizing effect which results from the thinning of the intermediate layer as the radial position increases. There is no analogous term in the single interface case, equation (3.4). A decrease in the magnitude of 6 tends to stabilize any particular perturbation to the interface between fluids 2 and 3. We deduce that at small radii, a perturbation of given wavenumber may be stabilized by surface tension, whereas at very large radius, the same mode may be stabilized by the rapid thinning of the annulus. As a result, for each mode there may only be a range of radii over which it is unstable. As the system evolves with time, the thickness of the intermediate layer, 6, decreases. However, the total volume of annulus fluid is fixed. Let us therefore define the nondimensional variable zy zyxwvu zyxwv where I/ = 2nR56 is the volume per unit depth of the annulus fluid and R, = 1 2 n M 3 T / Q ; A is thus a volume ratio, that of annulus fluid to the volume enclosed by the critical radius. Equation (4.6) may be written in terms of the non-dimensional variables A and 92 = R2/& as + - A g 3 / n ( n - 1 ) g - in(n2 - 1) 2~R,2 9 2 3 (1 A92'2/n) Q B,, - B, + (4.8) We have seen that there is a range of radii for which each mode is unstable, B,,/B, > 0 (n = 1,2,3...), for prescribed input conditions. Since n can only take integer values, the envelope of the unstable domain is given as a sequence of the arcs of the neutral curves B,, = 0, n = 1,2,3 ... (4.9) Bn with the end points of each arc defined by the condition Bn+l,/Bn+l= B,,/B,. In the limit of large n, this may be approximated by the solution of ;(2) =O coupled with equation (4.9). (4.10) zy zyxwv zyx zyxwvutsrq zyxwvuts 367 The formation of drops through viscous instability 0.6 - M 3 / M 2= 0. I 0.01 0.4 - Stable n=2 : A 0.2 I I n I .., '. '.'.. ,, .... ".. . .... . zyxwvuts zyxwvutsrqpo zyxwvutsrqponmlkjihgfe I I 8 0I ' '. -< I ............ ._ - _...._,"*............. # ...... I v.. . I --__ ' . . I 7 - _ _ ---.1';.:.; ....... ............. I........ L I "' -- - - - -..s.A..:.:.2..2 .-.. ...................i'................................................................. I 0 ----__ ------- _ _ _ _ _ _ ................. I 2 6 4 8 ~ ..c.*.I.I."I% 7 10 L% FIGURE 14. Stability diagram for the radial displacement of an annulus of fluid. These equations reduce, respectively, to 933 - n(n - 1) 1 (4.11) zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQPONMLKJIHGFEDCB 9 -n2 (n2 - 1) = 0, A 6A ~ + 2n (4.12) (2n - 1)g- - (2n2 - I) = 0. 6 The solution of this system of equations in A-93 space is shown as a solid curve in figure 14. The dashed lines define the stability domain for each mode (cf. (4.9)). It may be seen that for A greater than approximately 0.3, the system will always be stable, independent of 93. This corresponds to a relatively small volume of annulus fluid or small rate of injection, for a particular system of fluids. Physically, it means that stabilization at small radii by surface tension overlaps with the stabilization at higher radii caused by the thinning of the annulus . For lower values of A , there is a window of instability for a certain range of intermediate radial positions. For such values of A, as the radius increases, the most unstable mode may increase in wavenumber. While the perturbations remain small so that the linear theory applies, there may be a cascade to higher modes. Ultimately, at sufficiently large 9,the system becomes totally stable. Furthermore, for some values of A , a certain mode may grow and actually start decaying before the subsquent higher mode begins to grow and dominate. In the asymptotic limit of very large radius 9, total stabilization is predicted to occur for 3 (4.13) 9 2 9 --. - 24 The highest unstable mode has wavenumber given by (4.14) 368 zyxwvuts zyxwv zyxw zyxwvu zyxwv S. S. S. Cardoso and A. W Woods and hence this is the maximum number of drops which may form. The dotted curves in figure 14 represent a lower bound on the region of validity of the thin annulus approximation (cf. (4.5) and (4.6)). As seen, this is not a very restrictive simplification, particularly at small mobility ratio M3 /M2. Integrated linear growth rates As before (see 0 3.1), for very small initial perturbations, the most unstable mode may evolve owing to the change in radius of the base state. The growth rate of the perturbation should thus take into account the change in radial position of the interface with time. Including this displacement, the amplitude of the perturbation has the time-integrated form BO (4.15) zyxw zyxw Figure 15 shows the predicted evolution of the system for two different values of A . We have assumed here that noise in the system perturbs all modes with an initial disturbance of fixed amplitude Bo when they become unstable (equivalent to model B of 0 3). In figure 15(a), the magnitude of A is relatively large, and hence the system is unstable only over a very narrow range of 9.The transition from mode two to mode three occurs at small 9, but then the system rapidly becomes totally stable at larger B. The amplitude of the perturbation remains very small throughout. For a smaller value of A , the window of instability is wider, figure 15(b). A cascade of modes from two to three, then four and eventually five is predicted. The largest amplitude of the instability occurs for mode three, after which there is a gradual decay until complete stabilization at larger B. The amplitude of the perturbation is larger than in the previous case, but still very small in comparison with those observed in the single interface problem (see figure 10). In the limit of very small A, corresponding to a large volume of annulus fluid or a large injection rate, the behaviour of the system will approximate that of the single interface, analysed in 0 3. The quantitative criterion for this asymptotic limit may be derived by comparing (3.8) and (4.15) and is found to be A B 2 / n << 1. A further interesting result concerns the breakup of the intermediate layer and the consequent formation of drops. If the initial disturbance grows sufficiently for nonlinear interaction to become significant, then the perturbation will lead to the formation of a number of droplets. However, particularly for small n, there is a relatively large range of radii over which each linear mode is dominant, and the linear theory gives a good estimate of the number of drops: it corresponds to the dominant linear mode when nonlinear effects become important. These drops then gradually move and separate from each other as they are pushed outward by the radial flow. In 0 4.3, we compare these theoretical predictions on (i) total stabilization and (ii) drop formation with some experimental observations. The formation of drops through viscous instability 3 O ! 1 zyx 369 zyxwvuts zyxwvutsrqponmlk zyxwvutsrqpon zy zyxwvuts zyxwvu I I 10 100 w FIGURE 15. Amplitude of each mode as a function of the radial position of the annulus for ( a ) A = 0.2 and ( b ) A = 0.05. 4.2. Experiments A set of experiments involving three fluids was performed. The Hele-Shaw cell (see 0 3.2) was initially filled with a viscous fluid. A less viscous fluid was then slowly injected, so that it formed a small circular drop. Finally, the third highly viscous fluid was injected at a prescribed flow rate. Two different systems were used. The first consisted of very viscous natrosol solution as inner fluid, air as intermediate and cooking oil as outer fluid. The second system was silicone oil DC12500 as inner fluid, diluted glycerine as intermediate and silicone oil DClOOO as outer fluid (the 370 zyxwvutsr zyxwvuts zyxwv zyx S . S . S . Cardoso and A. W Woods numbers refer to nominal viscosity in centistokes). The fluids were dyed with different colourings so that the developing interfacial patterns could be observed. The relevant physical properties of the fluids used are given in table 1. The surface tension coefficient between the diluted glycerine solution and silicone oil was measured by the sessile drop method (Isenberg 1978). These different systems and the magnitude of the flow rates were selected so that different dynamical regimes could be observed. 4.3. Comparison with experiments 4.3.1. Suppression of the instability by annulus thinning Figure 16 shows a sequence of photographs of the displacement of an annulus of diluted glycerine by very viscous silicone oil, DC12500. The outer fluid is a less viscous silicone oil, DC1000, with mobility M3 = 4.28 x cm3 s g-’ (see table 1); this is approximately nine times smaller than that of the diluted glycerine, and an order of magnitude larger than the displacing silicone oil, so that our earlier theoretical simplifications are valid. The outer interface has therefore an unstable viscosity jump, while the inner interface is highly stable. The injection rate is Q = 1.38 cm2 s-l, the critical radius is R, = 0.019 cm and A = 3.2 x It may be seen that during the outward flow, there is a continuous thinning of the annulus of glycerine. The outer interface is apparently stable and no growing perturbation is visible. Let us compare this situation with the behaviour of a single diluted glycerine-silicone oil DClOOO interface. This is shown in figure 17; the flow rate is the same as in the experiment in figure 16. The difference in the dynamics is remarkable. The displacement is clearly unstable, and twelve fingers develop, growing to a length of approximately 2 cm. In figure 18, we present the theoretical predictions for the growth of the instabilities for the experiments in figures 16 (dashed lines) and 17 (solid lines). A permanent initial disturbance of amplitude Bo is assumed. The range of 9 shown corresponds to that where an instability is visible in the single interface experiment in figure 17. The observation of mode twelve throughout the experiment in 17 is well predicted by our calculations. These suggest that mode thirteen might become visible during the last frames, but this was not observed, possibly owing to nonlinear effects; these effects may also cause the amplitude to grow to 2 cm more rapidly. The dashed curves predict that the outer interface of the spreading annulus (figure 16) is not totally stabilized by the presence of the nearby inner stable interface. There is still a growing instability; however, the theory suggests that it grows much more slowly than in the single interface experiment. Indeed, there is almost an order of magnitude difference in the amplitudes of the instabilities (during the stage in which the instability is visible in figure 17); and this large difference accounts for the fact that a growing instability was not clearly seen in figure 16. We should mention that it is very difficult to achieve total stabilization experimentally. This would require a very large value of A, which we were unable to get with common immiscible fluids and feasible flow rates. Nevertheless, the results above show an O(1) effect of the thinning of the annulus on the growth of the instability, which results in the suppression of the instability on the scale of our experiments. zyx 4.3.2. Drop formation A further interesting experimental observation concerns the formation of drops. The eventual breakup of the intermediate layer of glycerine in the silicone oil DC12500, diluted glycerine, silicone oil DClOOO system would require a larger Hele-Shaw cell zyxwv zyxwvut zyx The formation of drops through viscous instability 371 zyxwvu zyxwv zy FIGURE 16. Displacement of silicone oil DC12500 (inner fluid), diluted glycerine (annulus fluid) and silicone oil DClOOO (outer fluid) at a flow rate Q = 1.38 cm2 s-’: ( a ) &? = 116, ( b ) 9 = 168, ( c ) 92 = 232, ( d ) W = 289 and ( e ) &? = 405. than the one used here. We have therefore used a system of fluids for which A is much smaller to show this effect. Figure 19 shows the displacement of an annulus of air by a very viscous natrosol solution. The outer fluid is cooking oil. It may be seen that a wavelike pattern develops at the unstable air-oil interface. The instability grows until its amplitude reaches the thickness of the air annulus; at this moment, the troughs of the outer fluid merge with the inner viscous fluid, and small bubbles of air form and separate. The breakup of the annulus of air results in the formation of fourteen bubbles which are then displaced by the outward radial flow. and R, = 0.033 cm ( Q = 27 cm2 s-l). Since 11%’ << 1, In this experiment, A the inner stable interface has only a secondary effect on the dynamics of the outer interface. The evolution of the air-oil interface should therefore be well described by the integrated growth rate for a single interface, equation (3.8) which is the limiting form of equation (4.15). The corresponding theoretical predictions are shown in figure 20. A constant noise level of amplitude A0 is assumed to trigger each - zyx 372 zyxwvu zyxwv S. S. S. Cardoso and A. W Woods zyxwvutsr zyxwvutsrq zyxwvut zyxwvu zyxwvut zyxw FIGURE17. Displacement of silicone oil DClOOO by diluted glycerine at a flow rate Q = 1.38 cm2 SKI: (a)9 2 = 84, ( b ) 92 = 163, ( c ) 92 = 221, ( d ) 92 = 279 and (e) B ! = 400. mode. The experiment began at W i= 61 (Ri= 2 cm) and hence the integration was performed for 90, > Bi.Our calculations suggest that mode 14 or 15 should dominate at the observed moment of rupture of the annulus. This is in excellent agreement with the laboratory observations above. The initial annulus of air in this experiment seems to be sufficiently thin for rupture to occur still within the linear growth regime. 5. Stabilization by controlling the rate of injection 5.1. Single interface We have seen that surface tension between the two fluids has a stabilizing effect on the viscous instability. In the case of a radial displacement, if the velocity of the interface is sufficiently small, the effect of surface tension is dominant and there is total stabilization. The interface may therefore be stable, even at large radial positions, by controlling the rate of injection of fluid. The upper bound on the flow rate such that the flow is stable may be calculated by requiring that the radial position of the The formation of drops through viscous instability 105 zyx z 373 zyxwvutsrqp zyxwvu zy I 200 250 zyxwvu 300 I I 350 400 9 FIGURE18. Amplitude of the dominant mode as a function of the radial position, for a constant disturbance of amplitude Bo for 2 > 1. Dashed line - theoretical prediction of (4.15) for the experiment in figure 16; solid line - theoretical prediction of (3.8) for the experiment in figure 17. interface is equal to the critical radius for instability (see 0 3). Using equation (33, we have 12~M2T Q= . For small radii, the flow rate can be relatively large whilst the interface is stable. This is because of the large curvature of the interface and hence a large stabilization effect by surface tension. However, as the radial position of the interface increases, the flow rate must decrease for the flow to remain stable. The variation of this critical flow rate with time, from the moment injection was initiated, may be obtained by integrating (5.1), where Q = ~ z R Rand , R = 0 at t = 0. We have Q = ( l 2 7 ~ M 2 T )(~&/ ~ t)-ll3. The injection rate should therefore decrease as t-'l3 for the flow to be stable. 5.2. Spreading annulus We have seen that in the case of a spreading annulus, there are two stabilizing mechanisms: surface tension operates at small radii and the continual thinning of the annulus at large radii. The requirement that at each radial position all modes are stable is expressed by equations (4.9) and (4.10). The numerical solution to these equations, represented in A-$2 space in figure 14, determines the maximum flow rate for which the flow is still stable. It is convenient to define the new non-dimensional variables zy 374 zyxwvu zyxwv S. S. S. Cardoso and A. W Woods zyxwvutsrq zyxwvuts zyxwvu zyxwvu FIGURE 19. Displacement of very viscous natrosol (inner fluid), air (annulus fluid) and cooking oil (outer fluid) at a flow rate Q = 27 cm2 s-l: ( a ) 92 = 61, ( b ) 92 = 115, ( c ) 3 = 161, ( d ) 92 = 173 and ( e ) 92 = 211. Figure 21 shows the variation of the maximum injection rate Q' with radial position R'. As in the case of the displacement of a single interface, at small radii, the large stabilizing effect of surface tension allows the injection rate to be relatively large. Then, as the radius increases, the flow rate must decrease. However, at still larger radii, the flow rate may increase again. This is because the stabilization caused by the thinning of the annulus becomes increasingly important. The allowed variation of the flow rate with time, from the moment injection was initiated, is shown in figure 22, where zyxwv zyxw (--) t* = M3 271 M2 v 3/2 (1271M3T)t. ( 5.5) The formation of drops through viscous instability zyx zyx 375 15 zyxwv zyx zyxwvuts zyxwvuts 61 111 161 21 1 9 FIGURE 20. Amplitude of the dominant mode as a function of the radial position, for a constant disturbance of amplitude A0 for W > 61. Theoretical prediction of (3.8) for the experiment in figure 19. I 0.1 I 1 10 R' zyxw FIGURE 21. Maximum injection rate for a stable flow as a function of the radial position. The limit of small radial positions At sufficiently small radial positions, the requirement that all modes are stable is equivalent to that of mode two being stable, since all higher modes are stable independently of the magnitude of the flow rate. It is then possible to determine 376 zyxwvutsr zyxwvuts zyx S . S. S. Cardoso and A. LK Woods 10-3 zyxwvutsrq zyxwvut 10-1 10-2 100 10' t* FIGURE 22. Maximum injection rate for a stable flow as a function of time. analytically the critical flow rate by solving (4.9) with n = 2. For R' << 1, we have Q. = -1 R' ' The time dependency is given by -113 Q'= ($t*) . (5.7) As expected, these solutions are similar to the corresponding ones for a single interface, (5.1) and (5.2). Solutions (5.6) and (5.7) are represented as dotted curves in figures 21 and 22, respectively. It may be seen that this asymptotic solution is a good approximation to the full solution for R' < 0.2, t' < 0.01. zyxwv The limit of large radial positions Our asymptotic analysis in 5 4.1 showed that for sufficiently large radial position, the stability domain is defined by equation (4.13). In terms of our new starred variables, we have in this limit Q * =-R' 2 (5.8) 3 The integration in time leads to (5.9) Solutions (5.8) and (5.9) are represented as dashed curves in figures 21 and 22, respectively. The asymptotic solution is a good approximation for R' > 10, t' > 50. The flow rate may therefore increase linearly with time for the flow to be stable. zyxwvuts zyxwvu zyx z The formation of drops through viscous instability 377 6. Conclusions We have investigated the stability of the immiscible displacement of an intermediate layer of fluid, bounded by two other fluids of different viscosities. The linear stability analyses for a thin intermediate layer shows that there are two different modes: a global mode which is determined by the viscosities of only the bounding fluids and an internal mode associated with the local stability of the two interfaces. When the displacement is globally unstable, viscous fingers of the rear, less viscous fluid eventually penetrate across the intermediate layer and into the more viscous leading fluid; this process is shown to be similar to the classical Saffman-Taylor instability described in many previous studies on viscous fingering (Homsy 1987). In this paper, we have focussed on the the dynamics of the new internal mode. A globally stable flow but with a locally unstable leading interface was investigated. It was found that in a rectilinear geometry, this displacement is always unstable and, as a result, the intermediate layer will eventually break up. However, in a radial source flow, there are two stabilizing mechanisms. For small radii, surface tension stabilizes the flow (as in the displacement of a single interface, Paterson 1981) and for large radii, the continual thinning of the annulus suppresses the instability of each mode, leading eventually to stability. In this case, there is a window of instability at intermediate radii for a constant rate of injection of fluid. As the radius increases, each of these modes is stabilized in turn; hence, a cascade to higher-order modes is observed, as the annulus moves outward. If the growth of the instability during this stage is sufficiently large, nonlinear effects become important and the annulus of fluid will rupture to form a number of drops; this number is given by the dominant linear mode at that time. In contrast, if the instability remains sufficiently small, then it may be totally stabilized once the radius becomes sufficiently large, 9 2 3/(2A), where A is inversely proportional to the volume of annulus fluid (see equation (4.7)). We also show that by suitably varying the magnitude of the rate of injection of fluid, the flow may remain stable for all radial positions. A series of experiments was carried out, using a circular Hele-Shaw cell. The different qualitative behaviours of a three-fluid system described above were confirmed. Quantitative comparison between the experimental results and the theoretical predictions for the cascade of modes and the number of drops formed shows remarkable agreement. zyx This work was begun while the authors were at the 1993 GFD summer program at the Woods Hole Oceanographic Institution. We would like to thank Professor J. Keller for helpful suggestions on the instability calculations, Dr J. R. Lister for useful discussion, and Dr J. Whitehead and Mr R. Frazel for use of the GFD laboratory at WHOI. The financial support of the EEC Science Program, WHOI and NERC are gratefully acknowledged. zyxwvutsrq REFERENCES CHOUKE, R. L., MEURS,P. VAN & POEL,C. VAN DER 1959 The instability of slow immiscible, viscous liquid-liquid displacements in permeable media. Trans. AZME 216, 188-194. GORELL,S. & HOMSY,G. M. 1983 A theory of optimal policy of oil recovery by secondary displacement processes. SZAM J. Appl. Maths 43, 79-98. HILL,S. 1952 Channelling in packed columns. Chem. Engng Sci. 1, 247-253. HOMSY,G. M. 1987 Viscous fingering in porous media. Ann. Rev. Fluid Mech. 19, 271-311. ISENBERG, C . 1978 The Science of Soap Films and Soap Bubbles. Clevedon: Tieto Ltd. zyxwvu 378 zyxwvutsrq zyxwvutsr zyxwv S . S . S . Cardoso and A. W Woods MANICKAM, 0. & HOMSY,G. M. 1993 Stability of miscible displacements in porous media with nonmonotonic viscosity profiles. Phys. Fluids A 5, 1356-1367. MAXWORTHY, T. 1989 Experimental study of interface instability in a Hele-Shaw cell. Phys. Rev. A 39, 5863-5866. MCBIRNEY, A. R. 1984 Zgneous Petrology. San Francisco: Freeman, Cooper and Co. MUNGAN,N. 1971 Improved waterflooding through mobility control. Can. J. Chem. Engng 49,32-37. PARK,C. W. & HOMSY,G. M. 1984 Two-phase displacement in Hele-Shaw cells: theory. J. Fluid Mech. 139, 291-308. PATERSON, L. 1981 Radial fingering in a Hele-Shaw cell. J. Fluid Mech. 11, 513-529. L. 1985 Fingering with miscible fluids in a Hele-Shaw cell. Phys. Fluids 28, 26-30. PATERSON, REINELT,D. A. 1988 The effect of thin film variations and transverse curvature on the shape of fingers in a Hele-shaw cell. Phys. Fluids 30, 2617-2623. SAFFMAN, P. G. 1986 Viscous fingering in Hele-Shaw cells. J. Fluid Mech. 173, 73-94. SAFFMAN, P. G. & TAYLOR,G. I. 1958 The penetration of a fluid into a porous medium or Hele-Shaw cell containing a more viscous liquid. Proc. R. SOC. Lond. A 245, 312-329. L. 1986 Stability of Hele-Shaw cell flows. Phys. Fluids 29, 30863088. SCHWARTZ, TENNENT, R. M. 1971 Science Data Book. Oliver and Boyd. WILSON,S. D. R. 1975 A note on the measurement of dynamic contact angles. J. Colloid Znterface Sci. 51, 532-534.