Langmuir 1991, 7, 1827-1832
1827
A Network Simulation for Drainage of Static Foam Columns Ashok R. Bhakta and Kartic C. Khilar' Department of Chemical Engineering, Indian Institute of Technology, Powai, Bombay 400 076, India Received December 10, 1990. In Final Form: March 20, 1991
A simple cylindrical network consistingof regular hexagonsis used to simulate the flow during drainage in a static foam column. The dimensions of the network are calculated by consideringthe foam to consist of regular pentagonal dodecahedral bubbles. Simulations offer several new insights concerning the microscopicorigin of the nonlinearity and inflection in the drainage plot. It is shown that the contraction of Plateau border channels in the presence of bubbles of varying sizes causes the observed nonlinearity in the drainage plot. Inflection in the drainage plot can be observed only when the liquid volume fraction increases with column height. Other calculations concerning the effects of bubble size and column height are consistent with experimental observations.
I. Introduction The study of foam drainage and its subsequent collapse is of both scientific and industrial importance. In the design of many foam processes such as foam bed reactions, foam bed separations, enhanced oil recovery, etc., the knowledge of foam drainage and collapse is essential. The major processes that take place during foam collapse are (a)capillary-inducedlamella drainageand gravity-induced Plateau border (PB) drainage, (b) interbubble diffusion of gas, and (c) coalescence and rupture of bubbles. The simultaneous occurrence of all these processes makes the analysis of foam drainage and collapse really complicated. The most common method of experimentally studying foam collapse is to carry out drainage and decay measurements by using a static foam column. In this measurement procedure, the cumulativevolume fraction ( VD/ Vm)( VD,cumulativeof volume of liquid drained; V m , initial volume of liquid present) of initial liquid content that has drained to the bottom of the column is measured at different intervals of time (tD). The drainage curve is then obtained by plotting ( V D /V,) vs tD. The slope of this plot is a measure of the rate of drainage. A parameter, called drainage half-life, 7112 (time required for half of the initial liquid to drain out), is usually used to compare the persistence and stability of different foams.' Numerous drainage plots have appeared in literature for both aqueous and nonaqueous systems containing synthetic, natural, and biologicalsurfactants.l+ Here, we show two typical drainage plots in Figure 1to define the primary objective of this work. It may be noted from Figure 1,that while the two drainage plots are similar at longer times, there is a fundamental difference at shorter times. In plot 1,the drainage rate is initially constant and then decreases,while in curve 2, the drainage rate initially increases to a maximum and then decreases showing a definite inflection. Such nonlinearity in drainage plots
Figure 1. Typical drainage curves.
and the fact that some drainage plots show an inflection while others do not, have not yet been systematically analyzed. Early studies aimed at describingthe drainage plot were entirely empirical in nature and empiricalexpressionswere developed to fit the drainage curves.2-1 Later, models were developed based on gravity-induced flow in film^.^*^ Although these models were far from being rigorous, they neverthelesshave laid down the foundation for the analysis of gravity-induced flow through channels. In the later years, some systematic studies of the flow in Plateau borders were carried 0 ~ t . 7 - g However, among these only Haas and Johnson8considered un-steady-state drainage in stationary foams. Though valuable, their model was valid only for cases in which the liquid fraction decreased with height. Besides, their model could not predict the inflection and did not account for the surface mobility. On the other hand, some investigators have developed models in which all the Plateau channels and films are assumed to drain directly into the liquid pool.1o These models strongly underpredict the half-life since the delay due to flow through the network of PB channels is not accounted for.
* Author to whom all the correspondence may be directed.
(1)Bikerman, J. J. Foam; Spri er-Verlag: New York, 1973. (2)Amerine, N.A.; Martin, L. P.&Mattei, W. Ind. Eng. Chem. 1942, 34,162. (3)Clark, G. L.;Rae, S. Ind. Eng. Chem. 1940,32,1594. (4) Miles, G. D.; Shedlovrky, L.; Rose,J. J. Phys. Chem. 1946,49,93. (6) Jacobi, W.M.; Woodcock, K. E.; Groove, C. S. Ind. Eng. Chem. 1966,48,11. (6) Brady, A.; Row,S. J. Am. Chem. SOC.1944,66,1348.
(7)Leonard, R. A.; Lamlich, R. AIChE J. 1966,11,18. (8)Haas, P.A.; Johneon, H. F. Ind. Eng. Chem. Fundam. 1967,6,225. (9)Desai, D.;Kumar, R. Chem. Eng. Sci. 1988,88,1529. (10)Sarma, D.S.H. R.; Pandit, J.; Khilar, Kartic C. J.Colloid Interface Sci. 1988,124 (l),339. (0
1991 American Chemical Societv
1828 Langmuir, Vol. 7, No. 8,1991
n
A II
-
Bhakta and Khilar
*
2TR,
row 1
row 2
/-
'.',
/e--
Figure 3. Two-dimensional network (opened up cylindrical
la)
network).
(C
1
Figure 2. (a) Schematic of an idealized cylindrical network. (b) Schematic of a typical tetrahedral link. (c) Schematic of an idealized two-dimensional link. In this work, using a simple regular hexagonal network of Plateau border (PB) channels of equilateral triangular cross section, simulations for the Hagen-Poiseuille type flow were carried out to study the drainage from a static foam column. It is shown that the nonlinearity in the drainage plot is due to the movement of a wave of Plateau border contraction from the top of the column to the bottom in the presence of bubbles of varying sizes. In cases of monosize foam bubbles, most of the drainage occurs before the wave reaches the bottom to bring in nonlinearity in the drainage plot. The existence of the inflection can be attributed to the variation of liquid volume fraction along the height of the column. Simulations carried out to study the effects of bubble size and column height yield results in qualitative agreement with experiments. 11. Formulation of the Network Model A static foam column of cylindrical cross section is considered in this work. Foam bubbles can be considered as regular pentagonal dodecahedra. Three pentagonal films of adjacent bubbles meet at 120° to form Plateau borders and four such borders meet at a tetrahedral angle of 109O 28/16'!. Figure 2b shows the basic tetrahedral link in the actual network. The liquid in the foam column is distributed between the films (lamellae)and the Plateau borders. Lamellae drain into the Plateau borders and these borders drain into the borders situated below. In order to model this complex flow path, we adopt a twodimensional equivalent of tetrahedral link as shown in Figure 2c. The foam column can be idealizedto consistof a number of concentric cylindrical networks (see Figure 2a) containing such two-dimensional links. The characteristics of the network can be understood more clearlyby opening out the cylindrical network into a rectangular sheet. The following are the major considerations used in this network analysis: (1) It is a network of regular hexagons similar to those used by Princenll and Khan and Armstrong12 to study foam rheology. (11) Princen, H. M. J. Colloid Interface Sci. 1988, 91, 160. (12) Khan,5.A.; Annatrong, R. C. J . Non-Newtonian Fluid Mech. 1986,22, 1.
(2) The links are made of Plateau border channels of length 1pB and cross-sectionalarea ape. While ape changes with time, Ipg is assumed to remain constant. The crosssectional area of a PB channel is assumed to be uniform along its entire length. (3) The liquid volume fraction is evenly distributed throughout the network, unless otherwise stated. (4) The dynamic liquid pressure is assumed to be uniform throughout the entire network. The flow occurs due to static pressure gradient. 1I.A. Geometrical Dimension of Network. The length (Lpg) of a PB channel is taken to be equal to that of the side of regular pentagonal fiim. If "a" is the equivalent radius of the regular pentagonal dodecahedron, then it can be shown that IpB
= 0.817a
(1)
The verticalheight (h)occupied by a single link row (shown in Figure 3) is given as h = 1pB + IpB cos 60 = 1.225a
(2)
Therefore, the number of link rows to be used for a column of heights H is given as NL = H / h
(3)
The width (W)occupied by a single link unit (shown in Figure 3) is given as
w = 21,
cos 30 = 31'21pB
The number of links units per link row in a cylindrical network of radius R, is given as NR= 2rR,/ W
(4)
The total number of PB channels in a network is thus given as N = 3NLNR
(5)
1I.B. Plateau Border Cross-sectional Area (8pB). The area (ape) of the Plateau border cross section needs to be related to a macroscopic parameter, gas volume fraction for simulation purposes. This can be done in the following manner. A pentagonal dodecahedron has 30 edges and each edge is a part of a PB channel which is shared by three bubbles. So in effect, each bubble is associated with 10 PB channels. The gas volume fraction T can be written as T'
(4/3)ra3 ( 4 1 3 ) ~+ ~ '(1oapBlpB + 4 d 6 )
(6)
where S is half the thickness of a lamella. The terms inside the parentheses in the denominator represent the volume of liquid in the PB channels and in the lamellae.
Drainage of Static Foam Columns
Langmuir, Vol. 7, No.8,1991 1829 pips
A
Flow from
P B chonncls
above
= Flow to channels below Figure 5. Liquid flow in a single channel. qod
b (7) Figure 4. (a) Actut. shape of Plateau ,order. (b-d) 7, i o u s idealized shapes of a PB channel cross section.
Rearranging, we obtain a p B =7.17 q ( * -37 e ) u -34 7 r 6 ]
=
(7)
1I.C. Dimensionless Cumulative Drainage Volume
( VD/V-).The initial volume of liquid present can be found
from the centroid of the equilateral triangle. The correction for the mobility of PB walls is accounted in the factor 8, which is a function of the ratio ( p r / k ) ,where A is the surface viscosity. The changes in the PB cross-sectional area (ape) are considered next. As shown in Figure 5, a PB channel receives liquid from the films attached to it as well as from the PB channels linked to it from above. At the same time, the channel drains into the channels below it. An unsteady-state volume balance for the channel shown in Figure 5 can be written as
as where where apBi and 6i represent initial values of Plateau border cross-sectional area and half the lamella thickness, respectively. Rb is the radius of the circle encircling the regular pentagonal lamella. The second term is obtained by accounting for the fact that each channel has three films draining into it and one-fifth of the liquid in each of these films drains into a channel. In order to compute VD,it must be clearly noted that only the bottom moat PB channels (see Figure 3) contribute to the cumulative drainage volume VD. Therefore, the dimensions of the cross sections and the flow from these channels primarily determine the nature of the drainage plot. Thus
Where qobi is the flow rate from the ith link in the bottom most link row. IID. Flow and Cross-Sectional Area Changes in a Single PB Channel. A Plateau border is formed at the intersection of three films approachingeach other at 120° as shown in Figure 4a. In order to calculate the flow through a PB channel, various idealized cross sections have been ad0pted.~$~J3 Some of these are shown in parts b-d in Figure 4. We have adopted the equilateral triangular shape and the corresponding expression for flow rate.13 In this expression, effect of mobility of the wall has been taken into account. It may be stated that any one of the other two shapes would have been equally good for our purpose since the expressions for flow rate differ only in the numerical constant. The average velocity o through a PB channel of equilateral triangular cross section is given as (10) Where p is the density of liquid, p is the viscosity of the liquid, and r is the perpendicular distance to any side (13) De&,
D.;K w ,R. Chem. Eng. Sci. 1985,97,1361.
and
The flow rate of liquid entering the channel (Qi) can be written as Qi
= QiF + QiPB
The rate of flow of liquid from the lamellae qg can be obtained as
The factor, 3/5 shows up as three films drain into the periphery of a pentagon of which only one-fifth constitute a channel. &b2 is the area of the lamella that is thinning at the rate of -2(d6/dt). The rate of lamella thinning can be obtained from Reynold’s equation14
where 6 can be found from integration as 1
(14)
It may be noted that this expression is valid only if the bubble radius does not change with time. Here, Rp is the radius of the Plateau border and y is the interfacialtension. Reynold’s equation has been derived assuming that the gas-liquid interface ia tangentially immobile, a condition that may be valid at high concentration of surfactant. At low or moderate surfactant concentration, a correction factor accounting for the partial mobility of the interface (14) Reynold, O., Philoa. 7’”.R. Soc. London, A ISM, 177, 17.
Bhakta and Khilar
1830 Langmuir, Vol. 7, No. 8, 1991 A
M2
0.0
I
I
I
2 tD
IS)
-
1.0”
I
I
3
4
I 5
Figure 7. Typical computed drainage plots for monosized networks.
(3 ‘0
-
Figure 6. Typical drainage plots for a single channel.
can be used.16 In this simulation, this correction,however, has not been made. The network model is now completely defined with necessary hydrodynamicand mass balance equations.The simulation was carried out by numerically solving the relevant equations (11-14) for each PB channel using the second-order Runge-Kutta technique, which gives good results with a time step of 0.01 s. The fourth-order technique was not used since it would require significantly more computation due to the fact that qi in eq 11is not an explicit function of 2” and has to be computed at each stage by solving the same equations for the channel immediately above the one considered. The network was generated by using the data structures facility available in Pascal. The programming was simplified by the symmetry present in the network. As a result of this the drainage characteristics of each link in a given link row are identical. Once a bubble size ( a ) and a gas volume fraction (7)are chosen, all the network dimensions can be found out. In addition the following values of various parameters were used in this simulation
R, = 10” m 7
= 0.9 p,
p
= lo-’ Paos
6, = 10-~ m
= 2.5 X
Pawm
111. Simulation Results and Discussions
1II.A. Single PB Channel. Prior to considering a network, it would be worthwhile to consider the dynamics of the drainage of a single channel. In fact, some insight into the nature of the drainage can be obtained from a single channel simulation. Consideringa single channel, we suppose that the only liquid entering the channel is that from the films, i.e. qi = qiF. This assumption is true for a channel right at the top of the column. The simulation results are shown in Figure 6 for two cases: (i) qplt-o > q o p o and (ii) qiF/t-0 C qoltlo. As one observes from Figure 6, case i predicts (15) Iranov, I. B.;Dimitrov, D. S. Colloid Polym. Sci. 1974,252,982.
a drainage curve with an inflection, while case ii predicts a curve with a monotonically decreasing slope. For the initial condition of case i, eq 11dictates that apB and hence the drainage rate should initially increase. However, since q g is decreasing, at a certain stage qo becomes equal to qg. From then onward, since q g never exceeds qo (i.e. dr/dt C 0), qo goes on decreasing, thus giving a drainage curve with an inflection corresponding to the time at which q p first equals qo. In case ii, since qF never exceeds qo, UPB decreases continually. As a result, the drainage rate keeps on decreasing to yield a monotonicplot. The values of qg/t-O and qo t-0 depend on the initial relative distribution of liquid Letween the lamella and PB. Thus the nature of the drainage plot depends to a great extent on this distribution. 1II.B. Network of Channels. Simulations were carried out by using the cylindrical network described in the previous section for two cases,viz., networks containing monosized bubbles and those containing bubbles of different sizes (polysize). The results of the monosize simulation are presented first. III.B.1. Monosized Network. Figure 7 presents the simulation results for a monosize network for two different bubble sizes, viz., 0.5 and 1mm. One observes from these plots that they are linear for a major portion of the drainage. The nonlinearity (if any) appears only when most of the liquid has drained out ( VD/V , > 0.8). This implies that the cross section (apB)and hence the drainage rates of the channels situated in the bottom-most link row remain constant for most of the time. For a monosized network, all channels initially have identical dimensions. As a result the rate of flow in (qpB) and the rate of flow out (qo)are almost equalfor all channels except those that are situated in the top-most row. It may be noted here, that in these calculations, q p is negligible as compared to Q p B . For these channels, qpB is zero and hence the channels shrink in size (ape decreases) continuously resulting in a continuous decrease in qo. Such a decrease in qo in the topmost row causes the channels in the second row from the top of contract. Consequently, qo from the second channel decreases. In this manner, the channels present in link rows from top to bottom contract in succession. One may think of this process as a motion of a contraction wave from top row to the bottom row. The rate of drainage alters only when the crosssectional area of the channels in the bottom most link row changes and this happens when the contraction wave reaches the bottom row. As a result, the nonlinearity in drainage plot appears after a rather long period of time. One further observes from Figure 7 that the drainage plot for a larger bubble size becomes nonlinear at an earlier time. This can be explained from the motion of the
Langmuir, Vol. 7, No. 8,1991 1831
Drainage of Static Foam Columns Table I. Variation of Half-Life with Height bubble size a, m 5 X lo-' 5 X lo-' 5Xlo-'
column height
HIm 1.22X 10-2 3.66X 10-2 6.1OX10-2
initial gas vol fraction T 0.9 0.9 0.9
column radius
1.01
drainage half-life
R,, m
LR
qp, 8
3X 3 X 1C2 3X1W2
20 60 100
0.85 2.5 4
Table 11. Variation of Half-Life with Bubble Size (e) bubble size a, m IXlo-' 5 X lo-' 1 X lC3
column height
HI m 3.66X10-2 3.66 X 10-2 3.66 X
initial gas vol
fraction T 0.9 0.9 0.9
column radius RE,m 3X10-* 3 X lo-* 3X
LR 300 60 30
drainage half-life 71/29 s 100 2.5 0.45
Table 111. Details of Bubble Size Distribution H = 2.5 X 10-2 m 7 = 0.9 ag2 = 2.5 X lo-' name P2B
0.1, m
lo-' P ~ C io-' lo-' P3A lo-' P3B
NI 100500 100500 100600 100500
azl m 10" 2x10-3 10-3 10-3
N2 03, m 402 43 25 5 X lo-' 138 5 x lo-'
71/2,8 0.5 1.25 2016 3 201 0.75 N3
contraction wave in the following manner. For a given height of the column, the number of link rows are fewer if the bubbles are of larger size. Therefore, the contraction wave, having to travel a lesser number of rows, reaches the bottom row in a shorter period of time. In addition, the differencebetween in flow and out flow at the top is greater for a large channel causing a larger rate of contraction. Simulations were carried out to compute the drainage half-life(7112)for columns of different height and different bubble size. Table I presents the results obtained when height of the column was varied from 1.22 X 10-2to 6.10 X 1W2 m, keeping all other parameters the same. It can be noted that large heights are not taken to avoid a rather large amount of computations. One observes from this table that 7112 increases with column height. This follows from the fact that while the drainage rate from the bottom channels remains nearly same for all the three cases considered, the initial liquid content (V,)increases with column height and therefore VD/V, decreases causing a decrease in 71p. The drainage half-lifes, 7112, at different bubble sizes are presented in Table 11. It can be observed from Table I1 that 7112 decreases sharply with the bubble size. With the increase in bubble size, the channel size increases proportionately. The out flow,qo,however, is proportional to the fourth power of the channel size, and therefore a sharp increase in drainage rate is observed. As a result, 7112 decreases sharply. These observations are in qualitative agreement with experimental measurements. III.B.2. Polysized Network. A bubble size distribution has been incorporated into the model by approximating the foam column to consist of a number of concentric cylindrical networks each of which comprises bubblesof a different size. For a given height, the number of links in a network and their channel dimensions will be different for different bubble sizes. Simulations were carried out for two simple size distributions, one consisting of two sizes and the other consisting of three sizes with their Sauter mean size being equal. The mean size can be calculated as follows:
The number of bubbles (Ni) can be calculatedas ~N&R/
+D (SI
-
Figure 8. Computed plots for bisized and trisized networks.
10, where the total number channels is 3 N f l ~and each bubble accounts for 10 channels. The drainage plot is obtained by superimposing flow of liquid following different paths, i.e., through networks having PB channels of different size. The drainage plots are shown in Figure 8, while the corresponding size distributions are presented in Table 111. One observes from Figure 8 that unlike the monosized case,the drainage plots are nonlinear. As indicated earlier, for larger bubbles, the nonlinearityappears much earlier. Thus, when a group of networks is considered,the network consistingof larger bubbles will drain at a decreasingrate at a very early stage of drainage, thus rendering nonlinearity in the drainage plot. Therefore, when the flows are superimposed, the nonlinearity is seen quite early. The larger size of the bubbles present earlier is the nonlinearity observed. This can be seen from a comparison between plots P2B and P2C in Figure 8. The plot P2C shows greater and earlier nonlinearitybecause the largest bubbles it contains (2 mm) are larger than those in P2B (1mm). A similar explanation holds for the trisized networks P3A and P3B. The nature of the distribution also influences the half-life of the system. P3A has a half-life almost double that of P3B. This is because the latter contains a larger number of the bigger bubbles. III.B.3. Inflection in the Drainage Plots. As mentioned previously, it has been experimentallyobserved that some systems give drainage plots with an inflection while others do not. Drainage plots with inflection were obtained by Arbuzov and Grebenshchikov,lBwho prepared foam by sucking saponin solution through a porous disk at the top of the column. As Jacobi et a1.5 have rightly mentioned, the liquid fraction i n such systems will be higher at the top of the column as compared to that at the bottom. On the basis of the above premise and an arbitrary decreasing variation of gas volume fraction with column height as shown in Table IV, simulationswere carried out. The drainageplot obtained from such a simulationis shown (16) Arbuzov, K.N.; Grebenschikov, B. N., J . Phye. Chem. (Moscow)
1937, 10, 32.
B h k t a and Khilar
1832 Langmuir, Vol. 7, No. 8, 1991
Table IV. Variation of Gar Fraction with Height T height, m m T
height, m m
0.95 0.93 0.91
Oto5 6to10 10 to 15
15 to 20 20 to 25
0.89 0.85
1.0
downward if they contract, and linear if their dimensions remain unchanged. 2. A monosize distribution of bubbles with uniform liquid fraction gives a predominantly linear drainage plot. 3. The observed nonlinearity in most drainage plots is due to the presence of bubbles of different sizes, while the observed inflection in some drainage plots is probably due to a variation (increase) of liquid volume fraction with column height. 4. The qualitative variation in drainage half-life with system variables such as bubble size, column height, etc., are consistent with measurements.
G1ossary
0.0
L I
2
3
1D (51
-
L
5
Figure 9. Computed drainage plot for a column with varying liquid fraction with height. in Figure 9. As can be observed from Figure 9, there exists a definite inflection in the drainage plot. The microscopicorigin of inflection can be explained as follows. The network consists of channels of larger sizes at the top and, therefore, liquid drains at a higher rate into the bottom-most channels that control the drainage rate. These channels,which receive more liquid than they can drain, expand, and thus the drainage rate initially increases with time. Eventually, the contraction wave reaches the bottom-most row and the channels shrink in size. Consequently,the drainage rate decreaseswhen this happens and the plot acquires an inflection feature.
IV. Conclusions The following conclusions emerge from this study: 1. The factor influencing the nature of the drainage plot is the size of cross section of the channelsin the bottommost row which are in contact with the liquid pool. The curve is concave upward if these channels expand, concave
bubble radius area of PB channel cross section gravitational acceleration height occupied by one link half the lamella thickness length of a PB channel number of link rows number of links per link row number of bubbles gas pressure total flow rate of liquid into a channel flow rate of liquid leaving a channel flow rate of liquid enteringa channelfrom the films flow rate of liquid entering a channel from PB’s above it flow rate from the bottom-most links perpendicular distance from the centroid to a side of the triangle radius of PB radius of a f i i variable distance from center of film radius of cylindrical network time average velocity in a PB channel total gas volume width occupied by link volume drained at any time volume drained after complete collapse Greek Letters a inverse of dimensionless surface viscosity P density P velocity coefficient
Acknowledgment. A.B.thanks Mr. M.Pradhan for many useful discussions.
