Surfactant Adsorption and Marangoni Flow in Liquid Jets. 3. Modeling

Ind. Eng. Chem. Res. 2006, 45, 2235-2248

2235

MATERIALS AND INTERFACES Surfactant Adsorption and Marangoni Flow in Liquid Jets. 3. Modeling in the Presence of Micellar Surfactant Michael Weiss† and Richard C. Darton* Department of Engineering Science, UniVersity of Oxford, Parks Road, Oxford OX1 3PJ, United Kingdom

Turgut Battal and Colin D. Bain Department of Chemistry, UniVersity of Oxford, Chemistry Research Laboratory, Mansfield Road, Oxford OX1 3TA, United Kingdom

We present a hybrid computational fluid dynamics (CFD) model of surfactant adsorption and Marangoni flow at the rapidly expanding surface of a gravity-driven, laminar liquid jet that takes account of the presence of micelles. Micellar diffusion and finite rates of micellar formation and degradation are included in the model. This hybrid CFD model, which is developed within the CFD code FIDAP, is an extension of one recently published hybrid CFD model for adsorption and Marangoni flow in the absence of micelles. A boundary-layer treatment of the jet flow near the nozzle exitswhich assumes that the stress at the nozzle exit is not exceeded at any point on the jet surfacesshows that the adsorption of surfactant causes the surface velocity, and also the surface concentration, to increase linearly from the point where the jet detaches from the nozzle. The two major results of this theory are used in the hybrid CFD model; these are that (i) the rate of surface expansion of the jet remains finite at the point of detachment and (ii) there is a finite surface concentration at the nozzle exit. The analytical solution for transport of surfactant monomers to the jet surface in the presence of micelles assumes infinitely fast micellar breakdown. In the numerical model, finite micellar breakdown rates can be accommodated. Only one micellar species is considered. The hybrid CFD model is used to establish the order of magnitude of the micellar disintegration rate constant and to study the effect of variations in the mean aggregation number of the micellar species. Computed results are presented for a micellar solution of aqueous C16TAB. These results are compared with experiments. Good agreement between the model and the experiments can be obtained when the rate constant for micellar breakdown takes a value O(102) s-1, for an aggregation number of 90. 1. Introduction Surfactant monomers have a tendency to form surfactant aggregates, called micelles, above a certain critical bulk concentration. The micelles, which are generally considered surface-inactive, coexist with surfactant monomers in the bulk liquid. If the equilibrium between micelles and monomers is disturbedsfor example, through the expansion (or compression) of the free surface, which causes monomers to diffuse toward (or away from) the surfacesthe micelles release (or take up) monomers to reestablish the state of equilibrium. In addition, micelles diffuse in the bulk liquid at a rate that differs from that of the monomers. Thus, depending on the rate of micellar dissociation (or formation) and on the rate of micellar diffusion, micelles can affect the adsorption of surfactant monomers at the free surface. In two recent publications, we have presented a theoretical and experimental study of the dynamics of surfactant adsorption at the continually expanding surface of a laminar liquid jet.1,2 While the experimental study describes experiments under both * To whom correspondence should be addressed. Tel.: +44 1865 273002. Fax: +44 1865 283310. E-mail: [email protected]. † Present address: Anjou Recherche Veolia Water, Chemin de la Digue, B.P. 76, 78603 Maisons-Laffitte Cedex, France.

micellar and nonmicellar conditions, the modeling work concerns only monomer adsorption. In the present paper, we extend our previous modeling work toward higher bulk concentrations above the critical micelle concentration (cmc) and study the effect of micelles on the adsorption of surfactant monomers in the jet. As before in the nonmicellar case, the modeling work contains an analytical part and a computational fluid dynamics (CFD) part, and both parts are combined to form a hybrid CFD model. This hybrid CFD model is developed within the commercial CFD code FIDAP, one of FLUENT’s general-purpose CFD solvers that uses the finite-element method. We compare computed results with measured surface velocity and surface concentration data for the surfactant C16TAB. The velocity data were measured using laser Doppler velocimetry (LDV), and ellipsometry was employed to obtain surface concentration data.2 The analytical part was originally developed to describe the adsorption of surfactant from micellar solutions in the center of an overflowing cylinder.3 Before describing this analytical model and its adaptation to the jet, we recapitulate the most important features and results of our previous analytical modeling work of the nonmicellar case in what remains of this introduction section, since the model of the micellar case rests largely on the nonmicellar case.

10.1021/ie050931p CCC: $33.50 © 2006 American Chemical Society Published on Web 03/01/2006

2236

Ind. Eng. Chem. Res., Vol. 45, No. 7, 2006

Figure 2. Schematic representation of the boundary-layer development in the jet. Our boundary-layer treatment assumes fully developed parabolic flow outside the boundary layer. Inside the boundary layer, the velocity is constant in the absence of surfactant (dotted line) and increases linearly from the surface in the presence of surfactant (solid line).

Figure 1. Schematic drawing of the axisymmetric jet flow showing the lower portion of the nozzle and the free section of the jet. R0 is the nozzle radius, g is the acceleration due to gravity, and uj0 is the mean nozzle exit velocity.

The liquid jet issues from a circular capillary. The inner diameter of the capillary pipe is on the order of 1 mm, and its length was chosen to be 100 times the diameter, sufficiently long to ensure Poiseuille flow at the nozzle exit. The capillary pipe is attached to a liquid reservoir in which the liquid height is maintained at a constant level, so that gravity provides the driving force for the liquid flow. The liquid leaves the nozzle with a mean velocity that is on the order of 1 m s-1. The velocity profile at the nozzle exit is that of fully developed laminar flow with the usual no-slip condition at the wall. This no-slip condition is immediately relaxed when the fluid leaves the tube. Figure 1 shows a schematic drawing of the axisymmetric jet flow. In the absence of surface shear, when no surfactant is present, a region of zero-vorticity appears at the surface of the jet. As the jet travels downward, this zero-vorticity region diffuses toward the center, eventually equalizing the velocity across the jet (see Figure 2). In our boundary-layer treatment of the jet flow, we considered its initial growth, when its thickness is much less than the radius of the nozzle or the jet. Near to the nozzle,

the effect of gravity may be neglected. Fluid in the zero-vorticity layer is accelerated by the adjacent layer of faster moving liquid. The velocity profile of the liquid outside the zero-vorticity layer remains parabolic, and the velocity within the growing layer of zero-vorticity is a function of the axial coordinate only. We set the condition that the velocity be continuous at the edge of the layer. Applying continuity and performing a momentum balance on the boundary layer led to

(Rez* )

u/s ) 4

1/3

(1)

where u/s ) us/uj0 is the dimensionless axial surface velocity and z* ) z/R0 is the dimensionless axial jet coordinate.1 The mean nozzle exit velocity, uj0, and the nozzle radius, R0, serve as characteristic parameters, and Re ≡ 2R0Fuj0/µ is the Reynolds number. The cube-root dependence of us on z in liquid jets has been found by other authors, most noticeably by Scriven and Pigford,4 Goren,5 and Tillett.6 We have discussed their work in part 2 of this series.1 In the derivation of eq 1, we neglect gravity, which in practice serves to accelerate the jet flow. The cube-root dependence of the surface velocity on the axial distance, eq 1, is thus valid provided that gravity does not affect the vorticity-free layer, that is, us(dus/dz) . g, where g is the acceleration due to gravity. We also require the bulk of the jet flow to be unaffected by gravity, so that uj02 . 2gz. Both criteria are fulfilled for short axial distances from the nozzle, which is important, as we use the theoretical solution of the jet flow near the nozzle to

Ind. Eng. Chem. Res., Vol. 45, No. 7, 2006 2237

construct our hybrid CFD model in this paper. In the CFD model, gravity is fully accounted for.1 Our boundary layer treatment assumes further that the flow at the nozzle exit is that of fully developed laminar flow. Tillett’s boundary-layer treatment of emerging flows, which does not confine modifications of the radial velocity profile to the boundary layer, has demonstrated that, in the high Reynolds number limit, this simplification is fully justified and, for moderate downstream distances from the nozzle, the interaction between the boundary layer and the core of the jet flow can be neglected.6 It has been shown that, for creeping jet flows (Re ≈ 0), this assumption is inadequate.7 We note that our CFD model of the jet flow does not assume parabolic flow at the nozzle exit. The inlet of the computational flow domain, where fully developed laminar flow is applied, is situated upstream of the nozzle exit. Our CFD simulations of the water jet show that, at high Reynolds numbers (Re > 1000), the flow at the nozzle exit is indeed fully developed laminar. The comparison with our experimental data has further confirmed the existence of a parabolic flow profile at the nozzle exit.1,2 When a surfactant is present, its nonuniform surface concentration will give rise to a surface tension gradientsthe Marangoni stress. This Marangoni stress causes a velocity gradient in the layer adjacent to the surface, which is thus a region of reduced rather than zero-vorticity (see Figure 2). The velocity profile in the core of the jet, however, remains parabolic. From the boundary-layer treatment, we obtained for the surface velocity

Ma (Rez* ) (1 - 4z* ) (1 + 41 dz*d (Ca1 ))

u/s ) 4

1/3

1/3

1/3

(2)

where Ma ≡ |Ca-1 - Ca0-1| is the Marangoni number, Ca ≡ µuj0/σ is the Capillary number, and Ca0 ≡ µuj0/σ0 is the Capillary number at the nozzle exit (z* ) 0). Equation 2 includes terms for the mean Marangoni stress, Ma/z*, and its local value, d/dz* (Ca-1). This stress cannot be determined from the hydrodynamics, as it depends on the surface concentration, which in turn depends on the diffusion of surfactant to the surface. Diffusion to the surface occurs when expansion of the surface, du/s /dz*, reduces the subsurface concentration of surfactant below its equilibrium (bulk) level, and it takes place within a thin diffusion layer that moves with a velocity u/s ) u/s (z*) throughout its depth. The convectiondiffusion equation, in which we neglect axial diffusion and convection, is given in eq 14. The boundary conditions are that w ) ws at y* ) 0 and w f wb as y* f ∞, where y* ) y/R0 is the dimensionless coordinate perpendicular to the surface, w ) cM/F is the mass fraction of surfactant, and the subscripts s and b refer to the values at the subsurface and in the bulk, respectively. The solution is

{x

w - ws y* ) erf wb - ws 2

ReSc

du/s dz*

}

(3)

where Sc ≡ µ/FD is the Schmidt number. The dimensionless diffusion mass flux, q/diff ) qdiffR0/(FD), to the surface is then

q/diff

x

∂w | ) ) (wb - ws) ∂y* y*)0

/

ReSc dus π dz*

(4)

and the convective mass flux, q/conv ) qconvR0/(FD), in the surface, which balances diffusion to the surface, is

Muj0Γsat d / (u Γ*) FD dz* s

q/conv )

(5)

where Γ* ) Γ/Γsat is the concentration at the surface. We divided the jet flow into two zonessthe near-nozzle region, which we termed the region of detachment, and the region downstream of the region of detachment. In the downstream region, where the Marangoni stresses are small compared with the maximum value that is found at the nozzle exit, we assumed that the surface velocity has the surfactant-free value, eq 1. Assuming that wb . ws and balancing diffusion and convection, eqs 4 and 5, lead to

Γ* ) 0.244

(

)

FD Re2/3Sc1/2wbz*1/3 Muj0Γsat

(6)

We note that eq 6 is mistyped in part 2 of this series of papers, where eq 48 should be identical to eq 6 in the present paper. For low concentrations, the variation in surface tension was obtained from Henry’s law, and when applied to ionic surfactants with monovalant ions, such as C16TAB, this gives

RmTFD d 1 )Re2/3Sc1/2wbz*-2/3 2 dz* Ca 6.14µMuj

( )

(7)

0

From eq 7, it follows that d/dz* (1/Ca) f ∞ (and also dΓ*/ dz*) as z* f 0. However, at z* ) 0, dΓ*/dz* must remain finite and sufficiently small that the maximum Marangoni stress, which equals the stress value at the wall inside the nozzle, is not exceeded. In the presence of surfactant, our boundary-layer treatment of the hydrodynamics shows that the cube-root dependence of eq 1 must break down as z* f 0. On the basis of this stress limitation, we have shown that in the detachment region, very near the nozzle, the surface velocity must scale linearly with the axial distance. Using Henry’s law and applying the stress limitation at the nozzle exit, we have shown that, to leading order, the surface concentration exhibits the same functional dependence on the axial distance in the region of detachment. Matching the surface velocity and the rate of surface expansion from the detachment region solution with that from the downstream solution gives the linear equation for the surface velocity in the region of detachment that contains the length of the detachment region, which we defined as λ* ) λ/R0, namely

u/s )

64 (3Re )

z* λ*2/3

1/3

(8)

The matching of the surface velocity in the region of detachment and in the downstream region is illustrated in Figure 3. Equation 8 is given by curve b, and curve a is the surfactant-free jet, eq 1. To carry out the matching, we have considered the reduction in surface velocity outside the region of detachment that is caused by the Marangoni stress equivalent to a small axial displacement ∆ of the nozzle in a downward sense in eq 1, as shown by curve c.1 The solutions for the surface concentration in the downstream region and in the region of detachment were matched using a procedure similar to that used for the matching of the surface velocity.1 The result of this second matching procedure was a quadratic equation for the surface concentration at the point of detachment from the nozzle, Γ/0, and the length of the region of detachment, λ*.

2238


Figure 3. Matching of the surface velocity, u/s ) u/s (z*), in the region of detachment, curve b, with the surface velocity in the downstream region in the presence of surfactant, curve c. Curve a is the surface velocity of the surfactant-free jet. The term λ* is the length of the detachment region, and ∆ is a matching parameter. The data are for Re ) 1950 and submicellar C16TAB at a bulk concentration of cb ) 0.3 mol m-3.

A mass balance around the surface at the point of detachment, in which we balanced diffusion and convection, eqs 4 and 5, and in which we assumed that wb . ws,0 ) ws(z* ) 0), and using the detachment region solution for u/s , eq 8, enabled us to find a second relationship for Γ/0 and λ*. Combining the two relationships for the length of the region of detachment, λ*, and for the surface concentration at z* ) 0, Γ/0, gives

(

FD λ* ) 2.543 × 10-3 E Muj0Γsat

)

3/2

ReSc3/4wb3/2

(9)

where E ≡ 2RmTΓsat/(µuj0) is the Elasticity number, in which Rm ) 8.314 J mol-1 K-1 is the universal gas constant, T is the temperature, and Γsat is the saturation surface concentration, and

(

Γ/0 ) 0.0462E1/2

FD Muj0Γsat

)

3/2

ReSc3/4wb3/2

(10)

2. Analytical Solution for Infinitely Fast Micellar Break Down Micelles are surfactant aggregates that can be made up of 50-150 or more surfactant molecules, and they may adopt various forms, including spheres and cylinders and disklike and lamellar shapes.8 The formation of micelles takes place if the bulk concentration exceeds a critical value, the critical micelle concentration, ccmc. The formation of micelles is usually understood as a series of reaction steps that lead to the micelle species, as described by Aniansson et al.,9 involving formation and dissociation rate constants at each step. Two surfactant monomers react to form a dimer, which finds another monomer to form a trimer, and so on, until the first stable micelle species is formed. Stable micelles are also formed through coalescence of smaller agglomerates. A range of stable micelle species is generally assumed. The process of micellization is governed by two time scaless one that describes the release (or takeup) of single monomers by the micelles and one that relates to the complete dissociation (or formation) of the micellessusually referred to as the fast and the slow processes, respectively. The time scale of the slow process may be ∼100 times longer than that of the fast one.10 However, the relaxation times may show large variations, depending upon the chemical nature of the surfactant polar head and the electrolyte concentration and on the molecular weight

Figure 4. Schematic drawing of the concentration profiles of micelles and monomers adjacent to the surface for the case of infinitely fast micellar break down.

of the surfactant. For example, the characteristic time of the slow process for amphiphilic block copolymers, with a molecular weight that is 2 orders of magnitude larger than that of C16TAB, can exceed 10 h and that for a fast process can reach several minutes. On the other hand, for a surfactant with a short hydrocarbon chain (heptylammonium chloride), the relaxation time of the fast process can be less than 0.1 µs.11 The fast process affects the mean aggregation number, whereas the slow process adjusts the concentration of micelles.12 The simplest model of micellar aggregation, which considers only the relaxation time of the slow process, assumes that N monomers, S, form one micelle, Sm. Muller noted that the formation of a micelle in one reaction step is highly improbable, as it requires the simultaneous collision of a large number of monomers.13 This simple aggregation model was employed by Miller,14 by Danov et al.,15 and very recently by Breward and Howell.16 It forms the underlying basis of both the analytical model, in which micellar breakdown is assumed to be infinitely fast, and the numerical model, which considers finite micellar breakdown rates, presented in this paper. 2.1. Instantaneous Monomer-Micelle Equilibrium. We assume that all the surfactant present in excess of the cmc is in the form of micelles, and we further assume that only one micellar species is present that has the mean aggregation number N. The concentration of monomers in the bulk is taken to be that at the cmc. Additionally, we assume that the equilibrium between monomers and micelles is instantaneous, so that the monomer concentration can only fall below the value at the cmc when there are no micelles present. There is thus a layer, of dimensionless thickness h* ) h/R0, in which monomers diffusing to the surface are found but no micelles. The concentration profiles of monomers, w(y*), and micelles, wm(y*), in this diffusion boundary layer h* adjacent to the surface are schematically shown in Figure 4. We note that the concentrations are given as mass fractions, w ) cM/F and wm ) cmNM/F, where M is the molecular weight of the surfactant species, F is the density of water, and c and cm are the molar concentrations of monomers and micelles, respectively. Both the monomer and micelle profiles have a finite slope at the edge of the monomer boundary layer, y* ) h*, since, at this plane, the micelles are converted into monomers. The


dimensionless diffusion mass fluxes of monomers and micelles, q/diff and q/m,diff, respectively, at y* ) h*, are related by

q/m,diff|y*)h* )

D / q | Dm diff y*)h*

(11)

in which

q/diff|y*)h* )

∂w | ∂y* y*)h*

(12)

The error function in eq 17 takes on values between unity and zero, and a comparison of eq 17 with eq 4 shows that the presence of micelles enhances the flux of monomers to the surface above what it would be at the cmc. We show in section 2.5 that the flux is also greater than what it would be if all the surfactant were present as monomers, above the cmc. 2.3. Micelle Concentration Profile and Micelle Diffusion Mass Flux. The convection-diffusion equation for the micelles, again neglecting axial convection and axial diffusion, reads

and

∂2wm ∂wm q/m,diff|y*)h* ) | ∂y* y*)h*

∂y

(13)

Since D/Dm > 1, where D and Dm are the diffusivities of monomers and micelles, respectively, we have that (∂wm/∂y*)y*)h* > (∂w/∂y*)y*)h*, which is indicated in Figure 4. Expressions for the concentration gradients (∂w/∂y*)y*)h* and (∂wm/∂y*)y*)h* are obtained from the solution of the respective convection-diffusion equation. 2.2. Monomer Concentration Profile and Monomer Diffusion Mass Flux. The diffusion boundary layers are much thinner than the hydrodynamic boundary layer, so that we can assume that the velocity parallel to the surface is constant throughout the region in which micelle and monomer concentrations are changing and the velocity perpendicular to the surface is given by -y*f ′(z*), where f ′(z*) ) du/s /dz* is the rate of surface expansion.1 The convection-diffusion equation for the monomers, neglecting axial convection and axial diffusion, is given by

*2

+

∂wm Re Scmy*f ′(z*) )0 2 ∂y*

where Scm ≡ µ/FDm is the micelle Schmidt number. The boundary conditions are that wm ) 0 at y* ) h* and wm ) wm,b as y* f ∞ (see Figure 4). We obtain the solution

{ x {x

erfc y* wm )1wm,b erfc h*

(14)

The boundary conditions are that w ) ws(z*) at y* ) 0 and w ) wcmc at y* ) h*, where ws(z*) is the subsurface concentration and wcmc is the concentration in the bulk at the cmc (see Figure 4). With these boundary conditions, we obtain the solution

w - ws(z*) wcmc - ws(z*)

)

{x {x

} }

erf y*

Re Scf ′(z*) 4

erf h*

Re Scf ′(z*) 4

(15)

q/m,diff|y*)h* )

∂wm | ) ∂y* y*)h*

x

wm,b

∂w ) | ∂y* y*)h*

(wcmc - ws)

{

Re exp -h*2 Scf ′(z*) 4

xReπ Scf ′(z*) erf h*

{x

Re Scf ′(z*) 4

} (16)

}

Similarly, we obtain the solution for the diffusion mass flux of monomers at the jet surface, y* ) 0,

q/diff|y*)0

x { xRe4 Scf ′(z*)} (17)

∂w ) ) (wcmc - ws) | ∂y* y*)0 erf h*

Re Scf ′(z*) π

Re Sc f ′(z*) 4 m

(19)

Re Sc f ′(z*) π m

{

exp -h*2

Re Sc f ′(z*) 4 m

{ x

}

Re Sc f ′(z*) 4 m

erfc h*

}

(20)

The concentration of micelles in the bulk, wm,b, can be found from a mass balance of the whole surfactant inventory,

wm,b )

NM M c ) (c - ccmc) ) wT - wcmc F m,b F T

(21)

where the subscript T indicates the total concentration of surfactant molecules. 2.4. Matching of the Diffusion Fluxes of Monomers and Micelles. At the dividing plane, y* ) h*, the diffusion mass fluxes of monomers and micelles are related by eq 11. With eqs 16 and 20 for q/diff|y*)h* and q/m,diff|y*)h*, respectively, we obtain

Using eq 15, we can calculate the diffusion mass flux of monomers toward the surface at the plane y* ) h*,

q/diff|y*)h* )

} }

Re Sc f ′(z*) 4 m

The mass flux of micelles to the surface at y* ) h* is given by

2

∂ w Re ∂w + Scy*f ′(z*) )0 *2 2 ∂y* ∂y

(18)

W)

{x }

Scm Scm exp{-Ψ2} erfc Ψ Sc Sc erf{Ψ} Sc m exp -Ψ2 Sc

x

{

}

(22)

where

W)

wT - wcmc wcmc - ws(z*)

(23)

and

Ψ ) h*

xRe4 Scf ′(z*)

(24)

We have thus found a relationship that describes the dependence of the monomer concentration at the jet surface, in the form of a normalized concentration W, on the ratio of the diffusivities

2240


5, and the pseudodiffusive flux given by eq 25. The mass balance then reads

Muj0Γsat d / (u Γ*)|z*)0 ) (wT - ws,0) FD dz* s

Figure 5. Scpseudo, as a function of the normalized concentration, W, according to eqs 22 and 26 for a total concentration of cT ) 1.25 mol m-3 (wT ) 455.625 × 10-6) C16TAB, which has a critical micellar concentration of ccmc ) 0.92 mol m-3 (wcmc ) 335.34 × 10-6), and varying ratios of Scm/Sc ) 1.0 (a), 2.0 (b), 3.0 (c), 5.0 (d), and 7.5 (e).

of monomers and micelles, Scm/Sc, and a parameter Ψ that contains the position of the dividing plane, h*, and the surface expansion rate, f ′(z*). 2.5. Pseudodiffusivity. Another way of regarding the enhancing effect of the presence of micelles on the diffusion mass flux of monomers to the surface is through the definition of a pseudodiffusivity. We compare the flux of monomers to the surface, eq 17, with a pseudoflux calculated for the hypothetical case that all surfactant is present as monomer. The monomers would then diffuse to the surface with a dimensionless pseudodiffusivity, Scpseudo-1, and the concentration difference that drives the mass transfer would be wT - ws. From eq 4, we then have the flux

q/diff|y*)0 ) (wT - ws)

xReπ Sc

pseudof

′(z*)

(25)

where, this time, we have accounted for the subsurface concentration at the point of detachment in the expression for the pseudodiffusive flux, ws(z* ) 0) ) ws,0. The surface velocity is zero at the point of detachment, u/s (z* ) 0) ) 0. The surface expansion rate at that point is a constant value, since the surface velocity increases linearly within the region of detachment. From eq 8, we have

du/s 64 1/3 -2/3 |z*)0 ) λ* dz* 3Re

( )

x

Scpseudo

)

x

Dpseudo ) (1 + W) erf{Ψ} D

(

Γ/0 ) 0.0462E1/2

FD Muj0Γsat

)

3/2

ReScpseudo3/4(wT - ws,0)3/2 (29)

where Scpseudo is defined in eq 26 and wT - ws,0 is the driving concentration force of the diffusion of surfactant to the surface at z* ) 0, wT being the total concentration, as defined in eq 21. Equation 29 may be contrasted with the equation for Γ/0 in the absence of micelles, eq 10. Using Langmuir’s equilibrium isotherm, we can replace ws,0 with

ws,0 ) (26)

The ratio Scpseudo/Sc is plotted in Figure 5 as a function of the normalized concentration W, where we have used eq 22 to eliminate the unknown parameter Ψ from eq 26. We see from Figure 5 that Scpseudo/Sc e1, so that the value of the pseudodiffusivity, Dpseudo, is greater than the monomer diffusivity, D. The presence of micelles enhances the diffusive flux of surfactant. 2.6. Estimation of Γ/0 in the Presence of Micelles. As before in the absence of micelles, we determine Γ/0 ) Γ*(z* ) 0), the surface concentration at the point of detachment from the nozzle, from the matching of the solutions for Γ* in the downstream region and in the region of detachment and from a mass balance around the surface at the point of detachment. To a first approximation, we assume that Henry’s law, which states that the surface tension depends linearly on the surface concentration, suffices to describe the link between the two variables, and we assume that wb - ws ∼ wb in finding the downstream solution for Γ*.1 Although this assumption is likely to be invalid at a monomer bulk concentration that equals the cmc, it enables us to obtain an estimate for Γ/0 in the presence of micelles. The matching procedure gives a relationship between Γ/0 and λ*. In the mass balance around the surface at the point of detachment, we balance the convective flux in the surface, eq

(28)

Solving the mass balance around the surface at the point of detachment, eq 27, gives a second relationship for Γ/0 and λ* that is used to eliminate the length of the detachment region from the first relationship. The result of this procedure is

We compare eq 25 with eq 17 and obtain

Sc

x

du/s Re Sc | π pseudo dz* z*)0 (27)

k*Γ/0

(30)

(1 - Γ/0)

in which k* ) kM/F is the dimensionless Langmuir constant. At the point of detachment, from eq 26, Scpseudo is given as

Scpseudo )

Sc [(1 + W0) erf{Ψ0}]2

(31)

where

W0 )

wT - wcmc wcmc - ws,0

(32)

and

Ψ0 ) h/0

xRe4 Scf ′(z*)|

z*)0

(33)

Equation 29 becomes implicit in Γ/0 if we substitute for ws,0, W0, and Scpseudo, and it must be solved iteratively. It is plotted in Figure 6 for different total concentrations of C16TAB, given by curves b-e. Equation 22 is also plotted in Figure 6, given by curve a. The intersection of the curves for eq 29 with the curve for eq 22 provides us with the value pairs for Γ/0 (through W0) and Ψ0 at the various total concentrations wT. These values for Γ/0 are used in the hybrid CFD model in the


the nozzle exit gave an estimate for this contact angle. A normal stress condition on the jet surface, the Young-Laplace equation, accounts for the varying pressure along the jet axis due to the curvature of the jet surface. The tangential stress condition at the surface represents a balance between viscous traction and the Marangoni stress due to a gradient in surface tension. This stress condition couples the fluid mechanics and the adsorption of surfactant. The surrounding air phase is neglected.

Figure 6. Dimensionless parameter Ψ as a function of the normalized concentration W at Re ) 1950 according to eq 29. The curves are for cT ) 1.25 (b), 1.5 (c), 2.5 (d), and 4.0 mol m-3 (e) C16TAB. Curve a is given by eq 22 for Scm/Sc ) 7.5.

Figure 7. Surface concentration at the point of detachment, Γ/0, for C16TAB at Re ) 1950 as a function of the bulk concentration, cb, in the region below the cmc (4) and as a function of the total concentration, cT, in the region above the cmc (b).

presence of micelles, to fix the surface concentration at the point of detachment. The values thus calculated are plotted in Figure 7. 3. Hybrid CFD Model 3.1. Fluid Mechanics. The numerical treatment of the hydrodynamics of the jet flow in the presence of micelles is completely analogous to the case of monomer adsorption below the cmc. We have detailed the numerical fluid mechanics model in part 2 of this series1 and give a brief summary here only. We have used FIDAPsone of FLUENT’s general-purpose CFD solvers that employs a finite-element method to discretize the governing set of partial differential equationssfor our numerical investigation of the jet flow. The fluid mechanics of the jet flow are governed by the Navier-Stokes equations and by the equation of continuity, which we solve in their dimensionless forms. Gravity points in the positive axial direction, and Poiseuille flow is applied at the inlet of the flow domain. We have a no-slip condition at the nozzle wall, and symmetry conditions on the axial center line of the jet flow. A kinematic condition, which sets the normal velocity on the free surface to zero, is used to locate the position of the free surface, and a contact angle is applied at the downstream position. The contact angle, which is a free surface boundary condition required by the CFD model, defines the angle between the tangent to the surface at the downstream position and the vertical outlet plane. Applying Bernoulli’s law and the equation of continuity and accounting for the parabolic velocity profile at

Both the normal and the tangential stress condition account for the locally changing surface tension along the jet surface, which changes as a function of the surface concentration. Although the CFD code FIDAP is capable of predicting Marangoni flows that arise from applied surface tension gradients, it requires the user to specify additional relationships in order to predict the Marangoni flow that results from the adsorption of surfactant along the jet surface. We have implemented a surface equation of state (eq 42) that links surface tension to surface concentration in the CFD code, and the Langmuir isotherm is employed to couple the surface concentration and the bulk value. Fortran user-defined subroutines were used to implement the surface equilibrium thermodynamics.1 3.2. Computational Domain and Meshing. We use the same computational domain and mesh that was already employed for the computation of the jet flow in the submicellar case.1 The axisymmetric computational domain considers a short part of the long capillary nozzle, one nozzle radius in length, and we chose a length of 100 × R0 for the free jet section. The grid is generated using FLUENT’s mesh-generating program GAMBIT, and a regular mesh is used. The mesh around the point of detachment is dense so as to capture the steep velocity gradients in the axial flow direction with sufficient accuracy. Since the diffusion boundary layer is about 2 orders of magnitude thinner than the hydrodynamic boundary layer, dense meshing underneath the surface is also required. The distance from the row of surface nodes to the first row of nodes is thus only about 40 nm (for R0 ) 0.79 mm). The total number of elements of the mesh is 12 012. In contrast to the computations without micelles, the fully coupled Newton-Raphson iteration scheme, which updates the free surface position at each iteration step, required significantly more computational time (up to 10 h on a Windows 2000 PC with 512 Mb RAM and a 1.0 GHz processor) and memory, when the concentration field of the micelles in the jet was computed along with the distribution of monomers in the jet. The free surface position was adjusted by moving the nodes on the surface and underneath the surface along the so-called spines to satisfy the kinematic condition on the surface. 3.3. Bulk Diffusion and Convection of Monomers and Micelles. 3.3.1. Governing Conservation Equations. In addition to the convection-diffusion equation for the surfactant monomers, a second conservation equation for the surfactant micelles must be considered. Both equations contain source terms that serve to describe the disintegration of micelles and the resulting supply of additional monomers. We thus have

(

) (

)

Re 1 ∂ ∂w ∂w ∂w ∂2w Sc V* + u* ) r* + + s* 2 ∂r* ∂z* r* ∂r* ∂r* ∂z*2 (34) for the monomers and

(

)

2242


(

) ( ( )

∂wm ∂wm Re Sc V* + u* ) 2 ∂r* ∂z*

)

∂wm ∂2wm Dm 1 ∂ r* + + s/m (35) 2 D r* ∂r* ∂r* ∂z*

for the micelles. The dimensionless radial and axial velocities, V* and u*, respectively, are scaled using the mean nozzle exit velocity, uj0, and the radial and axial coordinates, r* and z*, respectively, are made dimensionless using the nozzle radius, R0.The dimensionless source terms are defined as s* ) sR02/ FD and s/m ) smR02/FD. 3.3.2. Micelle-Monomer Kinetics. The source terms in eqs 34 and 35 are formulated using kinetic expressions. To formulate these expressions, we assume that N monomers, S, form one micelle, Sm, according to kf

N ‚ STSm

(36)

kd

where kf and kd are the rate constants of the formation and disintegration of the micelles, respectively. From the principle of mass action applied to the reaction mechanism, eq 36, we have

dc ) N(kdcm - kfcN) dt

(37)

for the monomers and

dcm ) kfcN - kdcm dt

(38)

for the micelles. The source terms then read

s* )

MR02 dc NMR02 ) (kdcm - kfcN) FD dt FD

(39)

and

s/m

NMR02 dcm NMR02 N ) (kfc - kdcm) ) FD dt FD

(40)

where t is the time. The molar concentrations of monomers and micelles, c and cm, are converted into mass fractions, w and wm, and eqs 39 and 40 are then implemented into the CFD code through a Fortran user subroutine. 3.3.3. Micelle-Monomer Equilibrium in the Bulk. The equilibrium between monomers and micelles in the bulk is described by the equilibrium constant, K, which we employ to couple the rate constants of the formation and disintegration of the micelles, kf and kd. From the reaction mechanism, eq 36, we obtain

K)

cm,b cT - ccmc kf ) ) kd c N Nc N cmc

(41)

cmc

where K has the units of mol(1-N) m3(N-1). Using eq 41, we are able to relate the rate constants and micelle characteristics. For example, for C16TAB at cT ) 1.5 mol m-3 and N ) 90, we find a value of K ) 11.7026 mol-89 m267 for the equilibrium constant. Choosing kd ) 1.0 s-1, we thus have kf ) 11.7026 m267 s-1 mol-89. Our micelle model uses explicit expressions for the formation and dissociation of micelles, since we compute both the

monomer and micelle bulk concentration fields. Weinheimer et al. have not used explicit rate constants in the course of their studies on monomer and micelle diffusion.17 Instead, they have combined the mass transfer equations for monomers and micelles to obtain an equation which relates the total flux of monomers and micelles to the gradient of the total surfactant concentration. This equation contains an apparent diffusion coefficient, in which both the diffusion coefficient of monomers and that of micelles and the equilibrium constant appear. Using this flux equation, the authors were able to determine the apparent diffusion coefficient as a function of the concentration by fitting their flux equation to measurements. The intercept of the fitted flux equation gives the micellar diffusion coefficient, and the slope of the (linear) curve is related to the diffusion coefficient of the monomers. Chan et al. have used a similar approach.18 The authors have studied the solubilization kinetics of fatty acids in micellar detergent solutions. A five-step mechanism was considered that describes the initial solubilization processsdiffusion of micelles toward the fatty acid surface (1), adsorption of micelles onto the surface (2), reaction to form a mixed micelle containing both detergent and solubilized fatty acid (3), mixed micelle desorption (4), and mixed micelle diffusion away from the surface (5). Assuming that steps 4 and 5 are rate controlling, an equation was derived that gives a linear relationship between the reciprocal of the solubilization rate and the reciprocal of the solubility. This equation contains combined rate constants, which in turn contain the individual reaction rate and equilibrium constants of the various steps of the solubilization mechanism. Solubilization experiments confirmed that steps 4 and 5 are rate controlling, and the authors used these experiments to find the values of the combined rate constants. 3.3.4. Boundary Conditions for the Monomers. At the inlet of the jet flow domain, sufficiently far upstream from the nozzle exit plane, bulk concentration is assumed everywhere in the fluid, which, in the presence of micelles, is equal to the value at the cmc. The nozzle wall is made impenetrable for surfactant; that is, a zero-flux condition is applied at the nozzle wall. A flux condition at the free surface of the jet allows the surface to take up surfactant from the bulk liquid. This surface flux boundary condition balances diffusion to the surface by convection in the surface, so that the rate of surface expansion determines the rate of diffusion to the surface.1 We assume diffusion-controlled adsorption and use the Langmuir equilibrium adsorption isotherm to link the surface concentration to the bulk value at the subsurface. At the downstream position, the surface experiences continuous expansion due to the action of gravity. The establishment of equilibrium conditions at the downstream position of the surface is thus not permitted, and consequently, no constraints to the monomer concentration at that position of the flow domain are applied. 3.3.5. Boundary Conditions for the Micelles. The equilibrium between monomers and micelles is undisturbed at the inlet of the flow domain, sufficiently far upstream from the nozzle exit plane, where the concentration of micelles takes on the bulk value. We can determine this value from the mass balance of the whole surfactant inventory, eq 21. As before in the case of the surfactant monomers, the nozzle wall is made impenetrable for micelles also. Micelles are considered surface-inactive, which is expressed through a zero-flux condition that is directed along the outward unit normal of the surface. Recent work by Colgate and Bain19 has suggested that there may indeed be a flux of micelles to the subsurface and from there adsorption of monomer at the surface. In the case of charged micelles, a barrier to


adsorption was provided by the accumulation of surface charge. We have not taken such effects into account. Due to the nonequilibrium conditions at the downstream position of the jet, no constraints to the micelle concentration at that position are applied. 3.4. Equilibrium Thermodynamics. The surface tension, σ, is linked to the surface concentration, Γ, by means of the Frumkin surface equation of state. The coupling between the concentration at the surface, Γ, and the concentration at the subsurface, ws, is established through the introduction of Langmuir’s equilibrium isotherm. The dimensionless form of the surface equation of state then reads

(

)

ws 1 1 ) + E ln 1 Ca Caw k* + ws

(42)

3.5. Limiting Marangoni Stress and Finite Surface Concentration at the Point of Detachment. We have previously reported that direct coupling of surfactant adsorption and fluid mechanics in the CFD modelsthe link between the shear stress condition at the surface and the surface equation of statesled to failure of the numerical iteration scheme.1 The Marangoni stress, d/dz* (1/Ca), very near to z* ) 0 was found to exceed the shear stress at the nozzle wall, which results in the appearance of a closed vortex flow in the nozzle exit plane. With the occurrence of such a toroidal flow, convergence for the convection-diffusion equation could no longer be attained. The failure of the numerical scheme results from its inability to handle the step change in the subsurface concentration that should occur at z* ) 0. Axial diffusion and numerical upwinding cause this step change in subsurface concentration to be smeared out over a few nodes, resulting in large concentration gradients immediately downstream of the point of detachment. When these concentration gradients are converted to surface tension gradients, large accelerations result, and the numerical computation of the convection-diffusion equation diverges. This problem was tackled by combining the theoretical solution for the subsurface concentration at the nozzle exit with the far-field solution that is obtained from the numerical computation, and the same procedure is applied in the presence of micellar surfactant. On the basis of the supposition that reverse flow does not occur, we apply a limiting Marangoni stress at the point of detachment. We also apply the finite surface concentration at the point of detachment, Γ/0, which this time has been calculated using the extended theory presented in this paper to account for the presence of micelles. We then derive a numerical fit to the subsurface concentration, eq 46, which is used in the calculation of the surface tension profile. The limiting Marangoni stress is implemented into the numerical computation through the subsurface concentration profile, ws(z*), at the point of detachment. It may be shown that the solution for ws(z*), which is valid as z* f 0, is given by

4 ws ) ws,0 + (k* + ws,0)z* E

solution for ws(z*) with a two-parameter power-law equation,

ws ) C(z*)m

which was matched asymptotically with eq 43. The resulting relation, which gave a smooth solution for ws(z*) that holds everywhere along the jet surface, contains a dimensionless length scale, l* ) l/R0, that is related to the length of the detachment region, λ* ) λ/R0, and that marks the point of transition from the (linear) upstream solution to the downstream power-law solution (eqs 108 and 109 in part 2 of this series1). The simple power-law relation ceases to accurately represent the numerical downstream solution for ws(z*) at higher bulk concentrations, most noticeably above the cmc. In this case, the monomer subsurface concentration approaches the bulk value at the downstream position (that is, ws f wb as z* becomes large), and this behavior must be accounted for. We have therefore developed the following method, which provides a better fit of ws(z*) for higher concentrations. We define a second length scale, L* ) L/R0, which indicates the transition from the second to the third region of the downstream jet flow, in which ws ∼ wb. This length scale is formally defined as

L* )

() wb C

1/m

(45)

where C and m are the fitting parameters of the simple powerlaw fit. The matched solution for ws(z*) is then extended to include the second length scale, that is,

{

4 ws ) ws,0 + (k* + ws,0)z* E

} {

}{

l*2 l*2 + z*2

(1-m)/2

}

m L* L* + z* (46)

and is used to describe the monomer subsurface concentration profile if micelles are present. Equation 46 is implemented into the CFD code together with the surface equation of state through a Fortran user subroutine, so that the surface equation of state is linked to eq 46. In this way, a smooth function Ca-1(z*) is obtained for the surface tension profile that accounts for the limiting Marangoni stress at the point of detachment. For z* , L*, the second length scale term in eq 46 takes a value of unity, and the solutions for ws are given by eqs 43 and 44 in the region of detachment and in the downstream region, respectively. For z* . L*, which implies that z* . l*, we have that ws f wb. If z* ∼ L* and large enough that z* . l*, we find that ws is given by

ws ) wb(z*)m

{( ) wb C

1/m

}

-m

+ z*

(47)

or

(43)

where ws,0 ) ws(z* ) 0) is the subsurface concentration at the point of detachment, which is calculated from the value for Γ/0 using the Langmuir isotherm. In the absence of micellar surfactant and at low to moderate monomer bulk concentrations, the surface is far from reaching a state of equilibrium with the bulk, even at large distances from the nozzle exit, so that ws(z*) , wb everywhere at the surface. In this case, we have approximated the numerical downstream

(44)

ws-1/m )

C-1/m + wb-1/m z*

(48)

Defining

X ) C-1/m

(49)

Y ) wb-1/m

(50)

and

2244


results in

ws-1/m )

X +Y z*

(51)

where the dimensionless parameters 1/m, X, and Y are obtained from a least-squares fit of eq 51 to the numerical downstream solution for ws(z*) in the region where z* . l*. 4. Results and Discussion We present computed results that we compare with experimental data for the cationic surfactant C16TAB. The hydrodynamic flow parameters, uj0 and R0, are kept constant, as well as the total concentration of surfactant in the bulk, cT. We use the hybrid CFD model to study the effects of variations in the dissociation rate constant and mean aggregation number of the micelle species, kd and N, respectively, on the adsorption behavior of C16TAB in the jet. We note that values for the parameters of the micelle kinetic relations, in particular for the dissociation constant, kd, are rare, and we seek to establish the order of magnitude of this parameter. Furthermore, our model assumes that only one micelle species is present, having a mean aggregation number N. It is desirable to study the sensitivity of the model to variations in N. 4.1. Parameters of the Hybrid CFD Model. The physical and flow parameters that were employed in the computations are given in Table 1. Bain et al.,20 using the Stokes-Einstein relation, calculated the coefficient of diffusivity of the micelle species to give Dm ) 1.0 × 10-10 m2 s-1 at T ) 298.15 K. Their calculation assumed spherical micelles, with a radius of rm ) 25 Å. We use this value in our calculation. Lianos and Zana measured values ranging from 88 to 92 for the mean aggregation number of the micelle species.21 In this work, we assume a value of N ) 90, which represents the lower limit for the mean aggregation number used in our calculations, the upper limit being N ) 120. The resulting values for the equilibrium constant, K, at the given total bulk concentration of cT ) 1.5 mol m-3 C16TAB are summarized in Table 2, together with the values for the Table 1. Values of Physical and Flow Parameters Used in the Numerical Computations with C16TAB above the cmc R0 uj0 cT ccmc θ F µ σw M k Γsat D Dm Rm T g

0.79 × 10-3 m 1.24 m s-1 1.5 mol m-3 0.92 mol m-3 89.97° 103 kg m-3 10-3 kg m-1 s-1 72.8 × 10-3 N m-1 364.5 × 10-3 kg mol-1 0.23 mol m-3 3.9 × 10-6 mol m-2 7.5 × 10-10 m2 s-1 1.0 × 10-10 m2 s-1 8.314 J mol-1 K-1 293.15 K 9.81 m s-2

nozzle radius mean nozzle exit velocity total molar conc of C16TAB critical micelle conc of C16TAB downstream contact angle at z* ) 100 density of water viscosity of water surface tension of water molecular weight of C16TAB Langmuir constant of C16TAB saturation surface conc of C16TAB diffusivity of C16TAB (monomers) diffusivity of C16TAB (micelles) universal gas constant temperature gravitational acceleration

dissociation and formation rate constants of the micelle species, kd and kf, respectively. K was calculated using eq 41. We have chosen a value for kd, which, through the given values for K and N, determines the value for kf. In this way, parameter studies for kd and N could be carried out. Several values for Γ/0 at different total bulk concentrations for C16TAB at Re ) 1950 are plotted in Figure 7. These values are shown together with the estimates for Γ/0 at submicellar bulk concentrations (see part 2 of this series for details1) and illustrate that micelles enhance the mass transfer of surfaceactive surfactant monomers to the surface. The value for Γ/0 that corresponds to cT ) 1.5 mol m-3 C16TAB at Re ) 1950, from Figure 7, is Γ/0 ) 0.2212. This value is transformed to the corresponding subsurface concentration value, ws,0, by means of the Langmuir isotherm. We find the value ws,0 ) 23.8128 × 10-6, which is used in the hybrid CFD model. We have outlined before that at high concentrations in the bulk, the surface approaches a state of equilibrium at the downstream position. The simple power-law relation, eq 44, which approximated the numerical downstream solution ws(z*) well when only small to moderate amounts of surfactant were dissolved, does not yield a satisfactory representation of the computed downstream data for ws(z*) at high bulk concentrations. A second length scale, L*, which is defined in eq 45, was introduced to take the equilibration of the surface into account, when ws f wb as z* becomes large. Equation 51, which contains three fitting parameters X, Y, and 1/m, gave the approximation of the numerical downstream data for ws(z*). With the values for these three fitting parameters known, the matched solution for ws(z*) that spans the entire length of the subsurface is given by eq 46. The surface equation of state was linked to eq 46, and the computation was carried out. An iterative procedure was applied, during which several computations, usually fewer than ten, were necessary to attain stable values for the three fitting parameters X, Y, and 1/m. Their values are summarized in Table 3 for varying dissociation rate constants, kd, and mean aggregation numbers, N. The values for the two length scales l* and L*, which were computed using the values for X, Y, and 1/m, are also displayed in Table 3. 4.2. Surface Concentration. 4.2.1. Effect of the Dissociation Rate Constant. The expanding jet surface causes monomers to diffuse to the surface. This diffusion process in turn disturbs the equilibrium between the monomers and the micelles. As a consequence of this disruption, micelles release monomers in an attempt to reestablish the equilibrium between monomers and micelles. The whole process is thus governed by two time scalessthe time scale of surface expansion, given by the reciprocal of the surface expansion rate, (dus/dz)-1, and the time scale over which micelles disintegrate, given by the reciprocal of the dissociation rate constant, 1/kd. The surface expansion rate varies along the jet surface and has its maximum value within the region of detachment. Using our theory, eq 28, we can estimate this maximum value to give

Table 2. Values of the Mean Aggregation Numbers, N, the Equilibrium Constants, K, and the Dissociation and Formation Rate Constants, kd and kf, Respectivelya

a

Γ*(z*) Figures (curves)

N

K (m3(N-1) mol1-N)

kd (s-1)

kf (m3(N-1) mol1-N s-1)

Figure 8(a) Figure 8(b) Figures 8(c) and 9(a) Figure 8(d) Figure 9(b)

90 90 90 90 120

11.7026 11.7026 11.7026 11.7026 107.0798

1.0 10.0 100.0 1000.0 100.0

11.7026 117.026 1170.26 11702.6 10707.98

The total bulk concentration is cT ) 1.5 mol m-3 C16TAB (ccmc ) 0.92 mol m-3).

Ind. Eng. Chem. Res., Vol. 45, No. 7, 2006 2245 Table 3. Values of the Parameters of the Least-Squares Fit to the Numerical Downstream Solution ws(z*), Equation 51, and the Resulting Values for the Length Scales l* and L*a kd (s-1)

N

X

Y

1/m

l*

L*

1.0 10.0 100.0 100.0 1000.0

90 90 90 120 90

3.4516 × 109 6.7004 × 108 9.0086 × 105 9.0096 × 105 3.0616 × 107

2.0488 × 107 1.9998 × 106 2.7207 × 104 2.7098 × 104 7.8812 × 105

2.066 1.906 1.285 1.285 1.700

0.7498 0.6835 0.4276 0.4548 2.3741

168.47 335.05 33.11 33.25 38.85

a The total bulk concentration is c ) 1.5 mol m-3 C TAB (c T 16 cmc ) 0.92 mol m-3).

Figure 9. Surface concentration curves, Γ*(z*), for N ) 90 (a) and 120 (b), in comparison with experimental data (4). The total concentration is cT ) 1.5 mol m-3 C16TAB at Re ) 1950. The dissociation rate constant is kd ) 100.0 s-1 in both cases. Curve c is the relative deviation factor, F, defined in eq 52.

Figure 8. Surface concentration curves, Γ*(z*), for varying dissociation rate constants, kd, in comparison with experimental data (4). The total concentration is cT ) 1.5 mol m-3 C16TAB at Re ) 1950. The lines are for N ) 90 and kd ) 1.0 (a), 10.0 (b), 100.0 (c), and 1000.0 s-1 (d).

dus/dz ∼ 103 s-1, from which it follows that (dus/dz)-1 ∼ 10-3 s. With increasing distance from the nozzle exit, the surface expansion rate decreases rapidly and the time scale of surface expansion becomes larger. At a constant mean aggregation number of N ) 90, we have varied the value of the dissociation rate constant, kd, to investigate its effect on the adsorption behavior in the jet. The graphs in Figure 8, where four Γ*(z*) curves are shown in comparison with one set of experimental data, were computed for dissociation rate constants of 1 e kd e 1000 s-1, that is, 1 g 1/kd g 10-3 s. Thus, for the disintegration of micelles to enhance the transport of monomers to the surface in the detachment region (z* < 1), the dissociation rate constant, kd, must be on the order of 103 s-1. This case is shown by curve d in Figure 8, where both the rate of surface expansion and the breakdown of micelles contribute to the transport of surfactant monomers to the surface. At kd ) 100 s-1, curve c, the breakdown of micelles is too slow to contribute to the transport of monomers to the surface in the region of detachment and the adsorption process is controlled by surface expansion. Decreasing the value of kd to 10 and 1 s-1, curves b and a, respectively, reduces the enhancing effect of micelle disintegration near the nozzle further and the adsorption of surfactant is dominated by surface expansion up to a distance of z* ) 10 from the nozzle exit. At longer distances from the nozzle, z* > 10, where the rate of surface expansion is greatly reduced, the disintegration of micelles enhances surfactant adsorption at small values for kd, as shown by curve b. The transport of surfactant monomers to the surface at z* > 10 is dominated by the disintegration of micelles in the case of fast micellar break down, curves c and d. In these cases, the numerical data at the downstream position (z* > 50) seem to plateau at a surface concentration value of Γ* ≈ 0.75, and the surface appears to be saturated with surfactant.

We note here that in our analytical treatment of monomer adsorption in the presence of micelles, we assumed the case of infinitely fast micellar breakdown kinetics; that is, kd f ∞ and 1/kd f 0. A value for kd that is on the order of 100 s-1 gives the best agreement between computation and experiment (see also curve a in Figure 9, which is identical to curve c in Figure 8). However, the experimental data near the nozzle disagree with all four computed curves, which may be the result of an overestimation of the value for Γ/0 by the theory. Using the stopped-flow method and conductometric registration, Czerniawski22 measured the rate of micellization of C16TAB in water and gave a value of kd ) 47.6 s-1 at a temperature of 25 °C and for total concentrations of 1-10 mol m-3. This value agrees reasonably well with the value found from our CFD computation. Chan et al.,23 on the other hand, obtained a value of kd ) 0.4 s-1 at 30 °C for a total concentration of 1.5 mol m-3 C16TAB in water. This value is much lower than would be consistent with our data, shown in Figure 8. The asymptotic values for the surface concentration in the experiment and in the computation at large z* disagree. Above the cmc and in the case of C16TAB, we would expect an equilibrium value of Γ/eq ≈ 0.8 at larger distances from the nozzle. This value seems to have been observed in the experiments, but the computations exhibit a lower value. The Langmuir isotherm, which we use to describe the local equilibrium between subsurface and surface, is only approximate, since it neglects both electrostatic effects and interactions between molecular chains. The agreement between experiments and the numerical model in the downstream region of the jet flow, where the surface becomes saturated and these additional effects are important, would be improved with a more sophisticated isotherm,3 but at the expense of significant additional complexity in the numerical computations. 4.2.2. Effect of Mean Aggregation Number. The effect of the mean aggregation number on the adsorption of C16TAB in the jet was investigated at a fixed micelle dissociation rate constant of kd ) 100 s-1. The computations were carried out with two different mean aggregation numbers of N ) 90 and 120 (see Table 2 for the values of the equilibrium constants and Table 3 for the values of the downstream solution fitting parameters and length scales). The results of this parameter study for Γ*(z*) are shown in Figure 9, curves a and b, where curve a is identical to curve c in Figure 8.

2246


Figure 10. Monomer bulk concentration profiles, w(r*), at three different axial downstream positions, z* ) 50 (a), 25 (b), and 10 (c), for cT ) 1.5 mol m-3 C16TAB and Re ) 1950. The mean aggregation number is N ) 90, and the dissociation rate constant is kd ) 100.0 s-1.

Figure 12. Bulk concentration profiles, w(r*) and wm(r*), of monomers (4) and micelles (O) in the boundary layer at z* ) 10. The calculations are for cT ) 1.5 mol m-3 C16TAB and Re ) 1950. The mean aggregation number is N ) 90, and the dissociation rate constant is kd ) 100.0 s-1.

Figure 11. Micelle bulk concentration profiles, wm(r*), at three different axial downstream positions, z* ) 50 (a), 25 (b), and 10 (c), for cT ) 1.5 mol m-3 C16TAB and Re ) 1950. The mean aggregation number is N ) 90, and the dissociation rate constant is kd ) 100.0 s-1.

Figure 13. Bulk concentration profiles, w(r*) and wm(r*), of monomers (4) and micelles (O) in the boundary layer at z* ) 50. The calculations are for cT ) 1.5 mol m-3 C16TAB and Re ) 1950. The mean aggregation number is N ) 90, and the dissociation rate constant is kd ) 100.0 s-1.

It is evident from Figure 9 that the effect of the mean aggregation number on the adsorption behavior is negligibly small, as both curves a and b are virtually indistinguishable. Also plotted in Figure 9 is a relative deviation factor F(z*), curve c, given by

our model, the micelle species is treated as surface-inactive through the zero-flux boundary condition at the surface. The micelle concentration profiles in Figures 12 and 13 reflect this boundary condition at the surface, where the gradients ∂wm/∂r* exhibit a slope of zero. Comparing the two radial micelle concentration profiles in Figures 12 and 13, we see that at z* ) 10, the micelle concentration at the surface is higher than further downstream at z* ) 50. In this region of the jet flow, the transport of monomers to the surface is dominated by the disintegration of micelles due to the low rate of surface expansion (see curve c in Figure 8). The micelle concentration is thus diminished at z* ) 50. As a consequence of the release of monomers by the micelles, the monomer concentration is increased at z* ) 50. Thus, if the micelles disintegrate at a finite rate kd, then a boundary layer region exists over which both micelles and monomers are present together, as the monomer concentration falls below the cmc. The numerical solutions found by Breward and Howell exhibit the same behavior.16 This case may be contrasted with the assumption of infinitely fast micelle break down (kd f ∞), which we assumed in our theoretical treatment. As the rate of disintegration decreases, the concentration of micelles within the diffusion boundary layer increases. As a consequence of this increase in the micelle concentration at low kd, the monomer concentration decreases, since less monomers are supplied from the disintegrating micelles. Curves a and b in Figure 8 show this situation at low kd. When the rate of

F)

|Γ*(N ) 120) - Γ*(N ) 90)| × 100% Γ*(N ) 90)

(52)

which is a measure of the relative deviation of the two solutions along the jet surface. In eq 52, Γ*(N ) 90) is given by curve a and Γ*(N ) 120) is given by curve b. Although the relative deviation increases with increasing jet length, it is very small, reaching a value of only little more than 0.05% at the outlet of the computational domain. 4.3. Radial Bulk Concentration Profiles. Calculated radial monomer and micelle bulk concentration profiles, w(r*) and wm(r*), are shown in Figures 10 and 11, respectively, for three different axial downstream positions and for N ) 90 and kd ) 100 s-1. The profiles shown correspond to the Γ*(z*) curve c in Figure 8. The data shown in Figures 12 and 13, where the bulk concentration profiles w(r*) and wm(r*) within the boundary layer at z* ) 10 and 50 are displayed, correspond to curves a and c in Figures 10 and 11. The profiles of the micelle concentration show that the boundary layer for this speciessthe distance from the surface to the point at which the bulk value, wm,b, is reachedshas about the same thickness as that of the surfactant monomers. Within


Figure 14. Surface tension parameter, Ca-1(z*), for varying dissociation rate constants, kd. The total concentration is cT ) 1.5 mol m-3 C16TAB at Re ) 1950. The lines are for N ) 90 and kd ) 1 (a), 10 (b), 100 (c), and 1000 s-1 (d).

Figure 15. Surface velocity curves, u/s (z*), for varying dissociation rate constants, kd, in comparison with experimental data (4). The lines are for cT ) 1.5 mol m-3 C16TAB at Re ) 1950, N ) 90, and kd ) 1 (a), 10 (b), 100 (c), and 1000 s-1 (d).

disintegration is very low, so that kd f 0, the bulk concentration of micelles exists up to the surface. The surfactant monomer profiles are smooth and wellbehaved and do not show any oscillations in the solution. The profiles of the micelle species, on the other hand, exhibit some oscillatory behavior, which is particularly pronounced at the edge of the boundary layer within the bulk liquid. This behavior may be related to the discretization of the computational flow domain. Although the mesh allows for sufficient accuracy when solving the monomer mass transfer, this may not be the case for solving the micelle mass transfer. Computations on a mesh that has a greater density toward the edge of the boundary layer would be necessary to establish the impact of the discretization of the flow domain on the solution for wm(r*,z*). The data displayed in Figures 12 and 13 show that as the boundary layer grows, with increasing distance from the nozzle exit, fewer nodes underneath the surface are necessary to accurately compute the gradients at the surface. Redistributing the nodes within the boundary layer would allow for greater density at the inner edge of the boundary layer without increasing the total number of elements within the computational flow domain, which ensures that the computational efforts remain the same. Also, the upwinding scheme employed to suppress the occurrence of such oscillations, which adds a small amount of artificial diffusion in the streamwise direction only, may not sufficiently work for the micelle species transport problem. 4.4. Surface Tension. The curves of the computed surface tension parameter, Ca-1(z*), which are plotted in Figure 14, correspond to the surface concentration curves in Figure 8. We see that, for larger values of kd, the surface tension curves reach a plateau at long axial distances from the nozzle. This plateau defines the state of equilibrium at the surface for this particular surfactant. 4.5. Surface Velocity. The computed surface velocity curves, u/s (z*), that correspond to the computed surface concentration data in Figure 8 are plotted in Figure 15, curves a-d, along with one set of experimentally obtained surface velocity data. The four curves show the effect of the micelle dissociation rate constant, kd, on the velocity retardation at the surface. In general, we can say that the point of transition from the linear to the far-field region of the jet moves toward larger values of z* with an increase in kd. In the cases of curves a and b, where kd ) 1 s-1 and 10 s-1, respectively, this point of transition lies at about z* ) 1. At the highest value for kd, curve d, this point has moved downstream, to about z* ) 5. Curve c shows an intermediate behavior.

Experimental velocity data are available to within a distance of about one nozzle radius from the nozzle exit. At this downstream distance, the velocity reduction as a consequence of the surface shear is no longer significant and the surface velocity is approaching the behavior of the water jet. We were thus not able to generate experimental data that show the reduction in surface velocity at a total bulk concentration of 1.5 mol m-3 C16TAB. Larger quantities of surfactant had to be dissolved in order to identify the Marangoni flow in the jet experimentally.2 5. Conclusions In part 3 of this series of publications on the dynamics of surfactant adsorption in liquid jets, we have presented a hybrid CFD model that takes the presence of micelles into account. This model has been developed within the CFD code FIDAP and is an extension of the hybrid CFD model of adsorption from submicellar surfactant solutions that was presented in part 2 of this series. In the analytical model, which has been used to complement the CFD model in the flow region very near the nozzle exit, we have assumed infinitely fast micellar breakdown, so that the equilibrium between micelles and monomers is attained instantaneously. Using this analytical model, we have demonstrated that the breakdown of micelles enhances the transport of monomers to the surface. In the CFD model, on the other hand, kinetic expressions have been employed and the disintegration rate of the micellar species is finite. Good agreement between computed and measured data was found for kd ∼ O(102 s-1) and N ) 90. The effect of variations in N on the computed result was shown to be negligibly small. Acknowledgment We thank the EPSRC (GR/M73194) for financial support and Dr. Dimitrina Valkovska for helpful discussions. Nomenclature Dimensionless Groups Ca ) Capillary number, µuj0/σ E ) Elasticity number, 2RmTΓsat/(µuj0) Ma ) Marangoni number, |Ca-1 - Ca0-1| Re ) Reynolds number, 2R0Fuj0/µ Sc ) Schmidt number, µ/(FD) Latin Symbols C ) dimensionless power-law fit parameter c ) molar concentration (mol m-3)

2248


D ) coefficient of diffusivity (m2 s-1) F ) relative deviation function (%) f(z) ) dependence of us on z (m s-1) f ′(z) ) rate of surface expansion, dus(z)/dz (s-1) g ) gravitational acceleration (9.81 m s-2) h ) monomer concentration boundary layer thickness (m) K ) equilibrium constant (m3(N-1) mol1-N) k ) Langmuir constant (mol m-3) kd ) micelle disintegration rate constant (s-1) kf ) micelle formation rate constant (m3(N-1) mol1-N s-1) L ) transition length scale (m) l ) transition length scale (m) M ) molecular weight (kg mol-1) m ) dimensionless power-law fit parameter N ) mean aggregation number q ) mass flux (kg m-2 s-1) Rm ) universal gas constant (8.314 J mol-1 K-1) R ) jet/nozzle radius (m) r ) radial coordinate (m) rm ) micelle radius (Å) S ) surfactant molecule s ) general source/sink term (kg m-3 s-1) T ) temperature (K) uj ) mean jet velocity (m s-1) u ) axial jet velocity component (m s-1) V ) radial velocity component (m s-1) W ) normalized dimensionless concentration w ) surfactant mass fraction, cM/F X ) dimensionless downstream solution fit parameter Y ) dimensionless downstream solution fit parameter y ) transverse coordinate in boundary layer (m) z ) axial coordinate (m) Greek Symbols θ ) downstream contact angle (deg) λ ) length of detachment region (m) µ ) dynamic viscosity (kg m-1 s-1) F ) density (kg m-3) σ ) surface tension (N m-1) ∆ ) axial nozzle displacement (m) Γ ) surface concentration (mol m-2) Subscripts 0 ) point of detachment b ) bulk cmc ) critical micelle concentration conv ) convection d ) disintegration diff ) diffusion eq ) equilibrium f ) formation m ) micelles pseudo ) pseudodiffusivity s ) surface sat ) saturation T ) total concentration w ) water Superscript * ) dimensionless variables

Literature Cited (1) Weiss, M.; Darton, R. C.; Battal, T.; Bain, C. D. Surfacatant Adsorption and Marangoni Flow in Liquid Jets. 2. Modeling. Ind. Eng. Chem. Res. 2004, 43, 5203-5220. (2) Battal, T.; Bain, C. D.; Weiss, M.; Darton, R. C. Surfactant Adsorption and Marangoni Flow in Liquid Jets: I. Experiments. J. Colloid Interface Sci. 2003, 263, 250-260. (3) Valkovska, D. S.; Shearman, G. C.; Bain, C. D.; Darton, R. C.; Eastoe, J. Adsorption of Ionic Surfactants at an Expanding Air-Water Interface. Langmuir 2004, 20, 4436-4445. (4) Scriven, L. E.; Pigford, R. L. Fluid Dynamics and Diffusion Calculations for Laminar Liquid Jets. AIChE J. 1959, 5, 397-402. (5) Goren, S. L. Development of the Boundary Layer at a Free Surface from a Uniform Shear Flow. J. Fluid Mech. 1966, 25, 87-95. (6) Tillett, J. P. K. On the Laminar Flow in a Free Jet of Liquid at High Reynolds Numbers. J. Fluid Mech. 1968, 32, 273-292. (7) Dutta, A.; Ryan, M. E. Dynamics of a Creeping Newtonian Jet with Gravity and Surface Tension: A Finite Difference Technique for Solving Steady Free-Surface Flows Using Orthogonal Curvilinear Coordinates. AIChE J. 1982, 28, 220-232. (8) Everett, D. H. Basic Principles of Colloid Science; Royal Society of Chemistry: London, 1988. (9) Aniansson, E. A. G.; Wall, S. N.; Almgren, M.; Hoffmann, H.; Kielmann, I.; Ulbricht, W.; Zana, R.; Lang, J.; Tondre, C. Theory of the Kinetics of Micellar Equilibria and Quantitative Interpretation of Chemical Relaxation Studies of Micellar Solutions of Ionic Surfactants. J. Phys. Chem. 1976, 80, 905-922. (10) Danov, K. D.; Kralchevsky, P. A.; Ivanov, I. B. Equilibrium and Dynamics of Surfactant Adsorption Monolayers and Thin Films. In Handbook of Detergents. Part A: Properties; Broze, G., Ed.; Marcel Dekker Inc.: New York, 1999. (11) Noskov, B. A. Kinetics of Adsorption from Micellar Solutions. AdV. Colloid Interface Sci. 2002, 95, 237-293. (12) Hoffmann, H.; Nagel, R.; Platz, G.; Ulbricht, W. Zur Kinetik der Mizellbildung von Alkylpyridiniumhalogeniden. Colloid Polym. Sci. 1976, 254, 812-834. (13) Muller, N. Kinetics of Micellization. In Solution Chemistry of Surfactants; Mittal, K. L., Ed.; Plenum Press: New York and London, 1979; Vol. 1. (14) Miller, R. Adsorption Kinetics of Surfactants from Micellar Solutions. Colloid Polym. Sci. 1981, 259, 1124-1128. (15) Danov, K. D.; Vlahovska, P. M.; Horozov, T.; Dushkin, C. D.; Kralchevsky, P. A.; Mehreteab, A.; Broze, G. Adsorption from Micellar Surfactant Solutions: Nonlinear Theory and Experiment. J. Colloid Interface Sci. 1996, 183, 223-235. (16) Breward, C. J. W.; Howell, P. D. Straining Flow of a Micellar Surfactant Solution. Eur. J. Appl. Math. 2004, 15, 511-531. (17) Weinheimer, R. M.; Evans, D. F.; Cussler, E. L. Diffusion in Surfactant Solutions. J. Colloid Interface Sci. 1981, 80, 357-368. (18) Chan, A. F.; Evans, D. F.; Cussler, E. L. Explaining Solubilization Kinetics. AIChE J. 1976, 22, 1006-1012. (19) Colgate, D. M.; Bain, C. D. Adsorption Kinetics in Micellar Solutions of Nonionic Surfactants. Phys. ReV. Lett. 2005, 95, 198302. (20) Bain, C. D.; Manning-Benson, S.; Darton, R. C. Rates of Mass Transfer and Adsorption of CTAB at an Expanding Air-Water Interface. J. Colloid Interface Sci. 2000, 229, 247-256. (21) Lianos, P.; Zana, R. Fluorescence Probe Studies of the Effect of Concentration on the State of Aggregation of Surfactants in Aqueous Solutions. J. Colloid Interface Sci. 1981, 84, 100-107. (22) Czerniawski, M. The Double Layer Structure of Colloidal Electrolytes. VI. Rate of Micellization Process of Colloidal Electrolytes in Aqueous Solutions. Chem. Abstr. 1966, 64, 13430-13431. (23) Chan, S.-K.; Herrmann, U.; Ostner, W.; Kahlweit, M. On the Kinetics of the Formation of Micelles in Aqueous Solutions III. Ber. BunsenGes. 1978, 82, 380-384.

ReceiVed for reView August 11, 2005 ReVised manuscript receiVed January 6, 2006 Accepted January 24, 2006 IE050931P

Surfactant Adsorption and Marangoni Flow in Liquid Jets. 3. Modeling

Recommend Documents