aa

*Corresponding author. Tel.: #44-0161-200-4353; fax: #44-0161- 200-4399. E-mail address: j.j.cilliers@umist.ac.uk (J. J. Cilliers). 1This work represents the last research contribution of Prof. Ted Woodburn, who passed away in October 1999 Chemical Engineering Science 55 (2000) 4021}4028 Prediction of the water distribution in a #owing foam S. J. Neethling, J. J. Cilliers*, E. T. Woodburn1 Froth and Foam Research Group, Department of Chemical Engineering, UMIST, P.O. Box 88, Manchester, M60 1QD, UK Received 28 July 1999; received in revised form 14 February 2000; accepted 22 February 2000 Abstract A two-dimensional mathematical model is presented which describes the time-averaged steady-state water distribution in a coalescing, #owing foam. The model uses previous work that predicts foam #ow velocity and bubble coalescence. Drainage is described using the transient model of Verbist et al., extended into two-dimensions. The e!ect of coalescence is accommodated by introducing the concept of Plateau border length for each bubble. The formulation and solution technique of the highly non-linear partial di!erential equation is discussed in some detail. Simulations are shown that compare the water distribution for non-coalescing and coalescing foams. The model will have application for determining the liquid content of #owing foams, the entrainment of solids in #otation systems and for designing foam process equipment. ( 2000 Elsevier Science Ltd. All rights reserved. Keywords: Foam; Simulation; Liquid; Motion; Drainage; Model 1. Introduction Foams commonly occur in the process industries, both as a problem (for example in distillation columns) and for performing separations (for example froth #otation). Pre- diction of the liquid content of #owing foams is of par- ticular importance, as it determines to a large extent the foam stability and the entrainment of undesirable mater- ial into the froth. This study refers to a two-phase disper- sion of air in water, which although of particular relevance to froth #otation, is referred to as a foam to distinguish it from the three-phase solid-bearing froths. The physical situation being considered, is the steady- state drainage of water through a #owing foam, in which the independent variables are the horizontal and vertical positions x and y. Bubbles continuously enter the foam through the pulp}foam interface and #ow upward to either burst on the surface of the foam, or #ow unbroken across a weir. This simulation assumes that the foam structure is symmetrical across a vertical section, reduc- ing the simulation to two dimensions (x, y). The x coordi- nate is lateral across the foam increasing from the vertical wall remote from the weir (x"0) and the y coordinate increases upwards through the foam from the pulp}foam interface (y"0). Work has previously been done on the transient drain- age of liquid out of a foam (Verbist, Weaire & Kraynik, 1996), which considered a static foam into which pulses of liquid were introduced from the top. Their work used the formulations reported by Leonard and Lemlich (1965). The stability of foams is discussed in texts of Bikerman (1973) and Derjaguin (1989). Recently, Goldshstein, Goldfard & Shreiber (1996) have discussed liquid drain- age through foams, and Monnereau and Vignes-Adler (1998) recently developed an experimental technique for the detection of bubble coalescence in foams. Murphy, Zimmerman & Woodburn (1996) calculated bubble velocities in a foam, v(x, y), from predicted streamlines, using the Laplace equation for irrotational #ow. Neethling and Cilliers (1998) extended this tech- nique to predict and visualise coalescence in #owing foams. 2. Foam drainage model * general basis The physical situation being considered, is the steady- state drainage of water through a #owing foam, in which the independent variables are the horizontal and vertical 0009-2509/00/$ - see front matter ( 2000 Elsevier Science Ltd. All rights reserved. PII: S 0 0 0 9 - 2 5 0 9 ( 0 0 ) 0 0 0 5 4 - 3 positions x and y. The bubbles deform as they pass upwards through the foam and they are characterised by their volume as a function of position. In the foam, the faces of two adjacent bubbles are separated by thin water lamellae. At the junction of three lamellae a single Plateau border is formed which is a channel along which liquid can move. The walls of this channel are de"ned by a small radius of curvature. In a coalescing foam, the Plateau borders will vary in length, diameter and orientation across the foam. At the junction of four bubbles, the Plateau borders meet at a vertex, from which the draining water #ow is redis- tributed. The local pressure of water p w in the Plateau borders is determined by the air pressure p a in the bubble, the interfacial tension c at the water}air interface and the local radius of curvature r #) of the Plateau border boundaries. Krotov (1981) was amongst the "rst people to recognise the importance of the pressure gradients within the Plateau borders on drainage. The pressure di!erence across an interface with a radius of curvature, r, in only one direction is given by the well-known Eq. (1), p a "p w #c r . (1) At the faces of the lamellae, the radius of curvature tends to R and the pressure in the lamella water equals the bubble pressure, while the pressure in the Plateau bor- ders is lower because of the smaller radius of curvature r #) of their walls. A pressure gradient develops across the water in the lamella "lms, causing drainage into the Plateau borders. This drainage continues until the forces between the two lamellae interface equal the force exerted from the Plateau border. It is assumed that this equilib- rium is rapidly attained. In the systems being considered, the amount of liquid present in the lamellae is much smaller than that in the Plateau borders. The water content of the froth is thus considered to be determined solely by the volume of the Plateau borders. 3. Simulation model The simulation is in two directions x and y, which requires symmetry in the z direction. The water content in a di!erential control volume is determined by the product of the average cross-sectional area, A(x, y), and the total channel length per unit volume j(x, y, z). The volumetric fraction of water, e(x, y, z), follows: e(x, y)"A(x, y)]j(x, y). (2) The fraction of the bounding areas of the control volume intersected by its Plateau borders is also assumed to be e(x, y). u(x, y) is the velocity vector of the liquid, while v(x, y) is the velocity vector of the gas bubbles. For clarity of presentation the variations of A(x, y), u(x, y), v(x, y) and j(x, y) are implicitly taken to be dependent on position x and y and are simply written as A, u, v and j. A steady-state mass balance of water through the con- trol volume can be computed in terms of the transla- tional #ow vector u of water #ow in the borders. In the limit this leads to the key equation L[ju x A] Lx dx#L[juyA] Ly dy"0, £ ) (jAu6 )"0. (3) To determine j, the length of Plateau border per unit volume of foam, the local bubble volume is expressed in terms of R, the radius of the circumscribing sphere of a dodecahedron. The ratio of the edge length of the dodecahedron to its volume leads to j j" 5J3 p]t]R2 , (4) where t, the ‘Golden Ratioa, "(1#J5)/2+1.618. The cross-sectional area A i of channel i is given in terms of the radius of curvature r #),i of its curved walls. Following Verbist et al. (1996) A i "AJ3! p 2B]r2#),i . (5) This is written as A i "C2r2 #),i where C"SJ3! p 2 . The three forces which determine the #ow of water in an individual channel are: (i) the capillary e!ect, due to local variations in the curvature of the Plateau border walls; (ii) the gravitational force causing downward #ow; (iii) viscous drag on the water by the wall of the Plateau channel which is moving upwards at the local velocity of the bubbles, v, taken to be equal to the local super"cial velocity of air. 4. Quantitative formulation of stress gradients (i) The capillary stress gradient, F #!1 , is generated by variations in the water pressure, p w , which is determined by the radius of curvature, r #) , of the walls of the Plateau channel and the bubble air pressure p a . The value of p a is assumed constant throughout the foam. F #!1 "£ ) p w , (6) where p w "p a ! c r #) . 4022 S. J. Neethling et al. / Chemical Engineering Science 55 (2000) 4021}4028 From Eq. (5) and (6) F #!1,x "!Cc 2 A~1.5 LA Lx , F #!1,y "!Cc 2 A~1.5 LA Ly . (7) (ii) The gravitational stress gradient F(y) G is considered to act only in the y direction. F G,x "0, F G,y "!o w g. (8) (iii) The gradient in the wall shear stress, F D , is caused by the viscous drag force and is proportional to the di!er- ence between the mean velocity of the draining water in the Plateau borders, u, and the motion of the wall, v, as well as the bulk liquid viscosity l. The velocity, v, of the surface of the Plateau border is taken as the super"cial velocity of the air through the foam. The proportionality constant, which for streamline #ow in a circular duct is 8p, is taken by Verbist et al. (1996) to have a value of 50. The Verbist analysis used an additional factor of 3 to allow for the increased distance of travel in the vertical y coordinate direction. Although his treatment is only for one spatial dimension, in the light of the uncertainties associated with the distribution of the angles of inclina- tion of the Plateau borders, we have also increased the proportionality constant by the factor 3. F D,x "!50]3]k](u6 x!v6 x) A , F D,y "!50]3]k](u6 y!v6 y ) A . (9) Under steady conditions these stress gradients are in balance, leading to F #!1,x #F D,x "0, F #!1,y #F G,y #F D,y "0. (10) Substituting from Eqs. (7)}(9) into Eq. (10) and solving for the local average water velocities, u6 x "!k 2 A1@2 LA Lx #v6 x , u6 y "!k 2 A1@2 LA Ly !k 1 A#v6 y , (11) where the constant parameters k 1 and k 2 are k 1 " owg 150k (m~1s~1) and k 2 " Cc 300k (ms~1). (12) Substituting the local average water velocities as deter- mined by Eq. (12) into the local water #ow conservation equation of Eq. (3) gives A in terms of the motion of the bubbles in the foam, v, and, j, the local Plateau border length per unit volume LjA[!k 2 A~1@2LA/Lx#v6 x ] Lx , #LjA[!k2A~1@2LA/Ly!k1A#v6 y] Ly , "0. (13) Expanding the partial derivative of the continuity equa- tion with respect to x and y in Eq. (13) gives jC!k2C 1 2 A~1@2A LA LxB 2#A1@2 L2A Lx2D#v6 x LA Lx #A Lv6 x Lx D #C!k2A1@2 LA Lx #v6 x AD Lj Lx , jC!k2C 1 2 A~1@2A LA Ly B 2#A1@2 L2A Ly2 D !2k 1 A LA Ly #v6 y LA Ly #A Lv6 y Ly D #C!k2A1@2 LA Ly !k 1 A2#v6 y A] Lj LyD"0. (14) Collecting terms in A and its derivatives, keeping in mind that the functional dependence on x and y is implied, Eq. (14) reduces to the governing equation L2A Lx2 ](!jk 2 A1@2) #L 2A Ly2 ](!jk 2 A1@2) #LA Lx ]A!j k 2 A~1@2 2 LA Lx #jv x !Lj Lx k 2 A1@2B #LA Ly ]A!j] k 2 A~1@2 2 LA Ly #jv y !2]jk 1 A!Lj Ly ]k 2 A1@2B #AAj] Mv x Lx #Lj Lx v x #j]Lvy Ly !Lj Ly ]k 1 A#Lj Ly ]v yB "0 (15) It should be remembered that in Eq. (15) the air velocity components v x (x, y) and v y (x, y) are obtained from the previous numerical simulation of Neethling and Cilliers (1998) as are the radii R(x, y) of the circumscribing spheres of the polyhedral bubbles. From this j(x, y) and its partial derivatives are also available. Solving Eq. (15) is a boundary value problem in the single variable A, the cross-sectional area of the Plateau border, which has to be solved numerically. S. J. Neethling et al. / Chemical Engineering Science 55 (2000) 4021}4028 4023 5. Boundary values The boundary conditions have to be expressed such that A, or an equation containing A and its di!erential, are the only dependent variables. 5.1. Boundary condition 1: pulp}froth interface At the pulp}froth interface, the bubbles entering the foam are nearly spherical. The value of A at this bound- ary, A 0 , is determined from the void fraction of close- packed, uniform spheres. For closest packing jA 0 "0.26. 5.2. Boundary condition 2: cell walls These may be either vertical or tilted at an angle h to the vertical. There is no #ow of water or air through these surfaces. Substituting from Eq. (11), the water and air #ows at an impervious boundary inclined at h to the vertical are u6 x cos h"A!k2A~1@2 LA Lx #v6 xBcos h !u6 y sin h"!A!k2A1@2 LA Ly !k 1 A#v6 yBsin h, (16) which when combined give 0"A!k2A1@2 LA Lx Bcos h !A!k2A~1@2 LA Ly !k 1 ABsin h, since both v x cos h!v y sin h and u x cos h!u y sin h are equal to zero. 5.3. Boundary condition 3: upper froth surface The top bursting surface is not impermeable to air, but it is impermeable to the #ow of water. Taking u y to be zero at the upper foam surface then from Eq. (11) v y,501 "k 2 A1@2 501A LA Ly 501 B#k1A501 . (17) 5.4. Boundary condition 4: concentrate overyow weir The assumption made for this simulation is that at the weir, the liquid is carried over solely due to the motion of the foam, which has only a horizontal component velo- city. Hence u 9,8%*3 "v 9,8%*3 which implies: LA Lx "0. (18) 6. Numerical integration of the governing equation (Eq. (15)) Eq. (15) can be written in the form of a general linear second-order elliptic equation a(x, y) L2A Lx2 #b(x, y) LA Lx #c(x, y) L2A Ly2 #d(x, y) LA Ly #e(x, y)A"0, (19) where the coe$cients are a(x, y)"c(x, y)"(!jxk 2 A1@2), b(x, y)"A!jx k 2 A~1@2 2 LA Lx #jxv x !Lj Lx xk 2 A1@2B, d(x, y)"A!jx k 2 A~1@2 2 LA Ly #jxv y !jx2k 1 A!Lj Ly xk 2 A1@2B e(x,y)"Ajx Lv x L #Lj Lx xv x #jxLvy Ly !Lj Ly xk 1 A#Lj Ly xv yB. (20) The simulation requires an estimate of the gradient of A and the value of A between two points. Even for the highly concave sections of the curve, the slope between two points gives a good approximation of the gradient midway between them. However, the average of the values of A at both points, does not provide a good estimate of the value of A midway between the points. To improve the estimate of A, a parabolic approxima- tion using the two neighbouring points and a further point in the area of interest, was used. The "nite di!erence form of Eq. (21) is aA j‘1,l #bA j~1,l #sA j,l‘1 #/A j,l~1 #uA j,l "0, where the coe$cients are a"Aaj,l# b j,l * 2 B, b"Aaj,l! b j,l * 2 B, s"Acj,l# d j,l * 2 B, /"Acj,l! d j,l * 2 B, u"!(2a j,l #2c j,l !e j,l *2). (21) 7. Preliminary calculation of bubble size distributions Neethling and Cilliers (1998) have previously reported a visualisation simulation of the foam #ow allowing for coalescence. From this simulation, the super"cial gas 4024 S. J. Neethling et al. / Chemical Engineering Science 55 (2000) 4021}4028 Fig. 1. A single frame from a visualisation simulation of a #owing, coalescing froth, including dimensions velocity v(x,y) and bubble size distributions R(x, y) are available, for di!erent values of the bubble bursting para- meter a, and weir angles h. In this visualisation simulation, coalescence is said to occur when a criterion C b ’‚ %2 /‚(x, y). ‚ %2 is the side length of an undeformed bubble, and ‚(x, y) is the actual observed side length. Although this is an arbitrary heuris- tic, it has physical signi"cance, since the lamella between two adjacent bubbles will thin when it is extended and hence the probability of failure increases, leading to co- alescence. C b can vary between 0 and 1. C b "0 representing no coalescence, as an in"nite extension of the lamella would be required before coalescence occurs, while when C b "1 virtually no extension will cause lamella failure and hence the foam is considered too unstable to exist. Many researchers have noted the importance of "lm thickness on coalescence (e.g. Narsimhan & Ruckenstein, 1986a, b, and Bhakta & Ruckenstein, 1996). Narsimhan and Ruckenstein (1986b), for instance, assume that only thermal disturbances occur in the "lms, and theoretically predict the critical thickness for "lm failure. This assump- tion is not appropriate for the rapidly #owing foams being considered here. If it is assumed that the pressure exerted by the curva- ture of the Plateau border is at equilibrium with the forces exerted within the "lm (van de Waals, electric double layer, etc.), then the cross-sectional area of the Plateau borders and the thickness of the lamellae are directly related. Simulations of Plateau border dimensions within a #owing aqueous foam show that the cross-sectional areas vary only slightly over most of the upper portion of the foam. This indicates that the lamellae in this region of the foam and at the upper surface are of approximately equal thickness. Since the bubbles at the top surface are bursting, this indicates that most of the lamellae within the foam are very near to the critical thickness required for coalescence. The added disturbance to the lamellae due to elongation can therefore reasonably be used as a coalescence criterion. The cross-sectional area of the Plateau borders at the bursting surface, and thus over most of the foam, can be estimated from the #ux of air through the top surface. Appropriate simpli"cation of Eqs. (11) and (12), by ignor- ing the capillary pressure gradient, shows that A 501 +v y,501 (k]150)/o w g. This implies that the #ux of air released through the top surface of the foam due to bursting is directly related to the critical thickness of the "lms at rupture. Although the visualisation reports bubble sections as irregular polyhedra (Fig. 1 shows a single frame from a typical visualisation), the bubble volumes are reported in terms of the cross-sectional areas A B (x, y) of the equiv- alent spheres. The bubble size data is reported by "tting the previously obtained sizes, using an equation of the following form: A B (x, y)"(a B #b B x#c B x2#d B x3) (e B #f B y#g B y2#h B y3). (22) Once the bubble sizes have been established the total length of the Plateau borders per unit volume j t (x, y) is calculated using Eq. (4). 8. Computational procedure The highly non-linear character of Eq. (15) requires a double iteration procedure for its numerical integra- tion. The initial estimate of A j,l for the whole "eld were taken to be A j,l "150k o w g xv y,501 04j4500, 14l4360, A j,0 "0.26 j 04j4500, l"0. (23) These initial estimates were then used to obtain an initial estimate of the coe$cients of Eq. (21). These are explicit functions of A j,l and the known values of v x , v y , j, k 1 and k 2 , as reported in Eq. (21). The coe$cient determination in terms of a current set of A values can be thought of as the outer iterative loop. Once a set of coe$cients has been calculated, they are used in an inner loop to determine an improved set of A values, using sequential over relaxation (SOR) proced- ures (Press, et al., 1990). The inner loop iteration proced- ures were aimed at minimising residual values m j,l at each point. These are de"ned in terms of Eq. (21) as m j,l "aA j‘1.l #bA j~1,l #sA j,l‘1 #/A j,l~1 #uA j,l . (24) S. J. Neethling et al. / Chemical Engineering Science 55 (2000) 4021}4028 4025 The inner loop procedure iterates on A j,l however hold- ing an unchanged coe$cient set over the entire "eld. Within this the improved set of A values is determined by the SOR procedure A/%8 j,l "A0-$ j,l !u]mj,l u . (25) The convergence criterion was taken as m j,l /u]A/%8 j,l . The convergence factor u"1.9875 was determined using recommendations of the Press et al. (1990). The inner loop was considered to close when the convergence criterion was less than 10~4 over the whole "eld. Once the inner loop is satis"ed then the new A j,l values are used in the outer loop to compute a new set of coe$cient